Commit de88c7e5 authored by Jim Nutaro's avatar Jim Nutaro
Browse files

More physically reasonable model uses a current source. The extra capacitors...

More physically reasonable model uses a current source. The extra capacitors seem to help with stability. The new, long range issue may be slow inter generator oscillations that are not damped.
parent 60796adf
......@@ -54,8 +54,6 @@ void ElectricalData::configureGenerator(genr_t& genr_data)
genr_data.R = 1.0/0.05;
genr_data.Agc = genr_data.R/100.0;
genr_data.D = 1.0;
genr_data.ga = 10000.0;
genr_data.gb = 50000.0;
}
void ElectricalData::initialize_circuit()
......@@ -218,7 +216,7 @@ void ElectricalData::initialize_circuit()
{
if (nodes[entry.second].genr.genr_mw == 0.0 && nodes[entry.second].genr.genr_mvar == 0.0)
continue;
Complex Igenr = nodes[entry.second].V*nodes[entry.second].y;
Complex Igenr = nodes[entry.second].V*(nodes[entry.second].B+nodes[entry.second].y);
for (auto line: lines)
{
if (line.from == nodes[entry.second].ID)
......@@ -240,6 +238,8 @@ void ElectricalData::initialize_circuit()
double vg = abs(Vnode[entry.second])*cos(arg(Vnode[entry.second])+phase_angle[i]);
nodes[entry.second].genr.Pm0 += ig*vg;
}
nodes[entry.second].genr.Iamp = abs(Igenr);
nodes[entry.second].genr.Itheta = arg(Igenr);
configureGenerator(nodes[entry.second].genr);
}
delete Y;
......
......@@ -59,7 +59,8 @@ class ElectricalData
double H;
/// These are the speed control parameters
double Tspd1, Tspd2, R, Agc; // Speed gain
double ga, gb; // Hueristic exictation control parameters
/// Reference output current amplitude and phase
double Iamp, Itheta;
double D; // Friction for off nominal speed
/// Power for plant sizing information
double genr_mw, genr_mvar;
......
......@@ -28,9 +28,9 @@ test1: app
#./grid2ode ../data/ieee118cdf.txt
./grid2ode ../data/ieee14cdf.txt
#./grid2ode ../data/test1.txt
#./grid2ode ../data/test3.txt
#./grid2ode ../data/test2.txt
./m2c circuit.mo
#${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp cvode.cpp -lsundials_cvode
${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp sim.cpp poly.cpp event.cpp
app: ${OBJS}
......
......@@ -3,11 +3,13 @@
void event(double* const dq, double* const q, double tNow)
{
/* q[ia_2_2_3] = 0.0;
// Drop the the line from bus 2 to bus 3 in the IEEE 14 bus model
// Setting the line current to zero is the same as removing it.
q[ia_2_2_3] = 0.0;
q[ib_2_2_3] = 0.0;
q[ic_2_2_3] = 0.0;
dq[ia_2_2_3] = 0.0;
dq[ib_2_2_3] = 0.0;
dq[ic_2_2_3] = 0.0; */
dq[ic_2_2_3] = 0.0;
}
......@@ -16,8 +16,6 @@ static double sys_freq_hz = 60.0;
static const char* const phases[3] = { "a", "b", "c"};
// Set of voltage amplitude measurement variables
static set<string> voltageAmplitudeEstimators;
// Set of gyrator parameters
static set<string> gyrators;
// Default rate of change for a transformer tap
static const double TfRate = 1E-1;
......@@ -109,12 +107,6 @@ string generatorAngle(string busName)
return string("theta_")+busName;
}
// Internal generator excitation voltage
string generatorVoltage(string busName, string phase)
{
return voltage(busName,phase);
}
// Generator output current
string generatorCurrent(string busName, string phase)
{
......@@ -129,7 +121,7 @@ string sineAmplitudeMeasurement(string busName, string phase, string genName, st
* Five-Point Amplitude Estimation ofSinusoidal Signals: With Applicationto LVDT Signal Conditioning,
* IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 59, NO. 3, MARCH 2010
*
* This is almost certainly good enough if we 60 hertz as the signal frequency. We should never be
* This is almost certainly good enough if we use 60 hertz as the signal frequency. We should never be
* more than 1.5% away from that because the generators would be tripping offline otherwise.
* For busses attached to the generator, we'll use the generator frequency.
*/
......@@ -148,25 +140,24 @@ string sineAmplitudeMeasurement(string busName, string phase, string genName, st
return y1_gen+"\n"+y2_gen+"\n"+y3_gen+"\n"+Aest+"\n";
}
// Generator excitation voltage
string generatorSourceVoltage(string busName, string phase, double Vref, ElectricalData::bus_data_t data)
string generatorExcitor(string busName, string phase)
{
return "E_"+busName+"_"+phase;
}
// Generator injected current
string generatorSourceCurrent(string busName, string phase)
{
string result;
string freq = number(2.0*pi*sys_freq_hz);
string w = omega(busName);
string theta = generatorAngle(busName);
string gyrator = "g_"+busName;
if (gyrators.find(gyrator) == gyrators.end())
{
result = gyrator+"="+number(data.genr.ga)+"*exp(-"+number(data.genr.gb)+"*fabs("+w+"))*"+w+";\n";
gyrators.insert(gyrator);
}
if (phase == phases[0])
result += generatorVoltage(busName,phase)+"="+number(Vref)+"*(1.0+"+gyrator+")*cos(("+freq+"+"+w+")*time+"+theta+")";
result += generatorCurrent(busName,phase)+"="+generatorExcitor(busName,phase)+"*cos(("+freq+"+"+w+")*time+"+theta+");\n";
else if (phase == phases[1])
result += generatorVoltage(busName,phase)+"="+number(Vref)+"*(1.0+"+gyrator+")*cos(("+freq+"+"+w+")*time+2.0*pi/3.0+"+theta+")";
result += generatorCurrent(busName,phase)+"="+generatorExcitor(busName,phase)+"*cos(("+freq+"+"+w+")*time+2.0*pi/3.0+"+theta+");\n";
else
result += generatorVoltage(busName,phase)+"="+number(Vref)+"*(1.0+"+gyrator+")*cos(("+freq+"+"+w+")*time+4.0*pi/3.0+"+theta+")";
result += generatorCurrent(busName,phase)+"="+generatorExcitor(busName,phase)+"*cos(("+freq+"+"+w+")*time+4.0*pi/3.0+"+theta+");\n";
return result;
}
......@@ -184,14 +175,22 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
string powerDraw;
// Power demand at the terminal is sum of demand for each phase
for (auto phase: phases)
{
{
string estimate;
string amp = sineAmplitudeMeasurement(busName,phase,busName,estimate);
if (!amp.empty())
result += amp;
string Ie = generatorExcitor(busName,phase);
stateVars.insert(state_var_init_t(Ie,data.genr.Iamp));
result += "der("+Ie+")=10.0*("+number(data.vamp)+"-"+estimate+");\n";
result += generatorSourceCurrent(busName,phase);
if (powerDraw.empty())
powerDraw = generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
powerDraw = voltage(busName,phase)+"*"+generatorCurrent(busName,phase);
else
powerDraw += "+"+generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
powerDraw += "+"+voltage(busName,phase)+"*"+generatorCurrent(busName,phase);
}
result += Pe+"="+powerDraw+";\n";
stateVars.insert(state_var_init_t(generatorAngle(busName),data.theta));
stateVars.insert(state_var_init_t(generatorAngle(busName),data.genr.Itheta));
stateVars.insert(state_var_init_t(w,data.genr.w0));
stateVars.insert(state_var_init_t(Pm,data.genr.Pm0));
stateVars.insert(state_var_init_t(C,data.genr.Pm0));
......@@ -207,7 +206,8 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
result += "der("+Pref+")=-"+number(data.genr.Agc)+"*"+w+";\n";
/* result += "der("+generatorAngle(busName)+")=0;\n";
result += "der("+w+")=0;\n";
result += "der("+Pm+")=0;\n"; */
result += "der("+Pm+")=0;\n";
result += "der("+C+")=0;\n";*/
return result;
}
......@@ -222,13 +222,12 @@ string createBus(
// Small, minimum value for the capacitance relative to ground
string result = "";
bool hasGenr = (bus.genr.genr_mw != 0.0 || bus.genr.genr_mvar != 0.0);
if (hasGenr)
result += createGenerator(busName,bus);
// Generator the node equations
for (auto phase: phases)
{
if (hasGenr)
result += generatorSourceVoltage(busName,phase,bus.vamp,bus)+";\n";
string isum = "(";
string igenr = generatorCurrent(busName,phase);
// Load current to use if this bus has an inductive load
string iload = "il_"+busName+"_"+phase;
// Load voltage it the load is capacitive
......@@ -270,18 +269,12 @@ string createBus(
{
isum += "-"+v+"/"+number(bus.R);
}
isum += ")";
// Generator current, bus voltage is regulated by the generator
if (hasGenr)
{
result += generatorCurrent(busName,phase)+"=-"+isum+";\n";
}
isum += "+"+generatorCurrent(busName,phase);
isum += ")";
// Voltage is a state variable
else
{
stateVars.insert(state_var_init_t(v,bus.v[phase_num]));
result += "der("+v+")="+isum+"/"+number(bus.Cshunt)+string(";\n");
}
stateVars.insert(state_var_init_t(v,bus.v[phase_num]));
result += "der("+v+")="+isum+"/"+number(bus.Cshunt)+string(";\n");
// Add load current through the inductor as a state variable
if (bus.L > 0.0)
{
......@@ -297,8 +290,6 @@ string createBus(
// Next phase
phase_num++;
}
if (hasGenr)
result += createGenerator(busName,bus);
return result;
}
......
......@@ -75,7 +75,7 @@ double delay(int var, double h, double t)
int main(int argc, char** argv)
{
int reports = 1, delays = 1;
double hmax = 0, h, hReport, hEnd, hDelay;
double h, hReport, hEnd, hDelay;
if (argc != 7)
{
cout << "Must supply quantum, cint, and tend" << endl;
......@@ -129,7 +129,6 @@ int main(int argc, char** argv)
q[i] += p[i].extrapolate(h);
}
// Advance the clock
hmax = ::max(hmax,h);
tNow += h;
// Store new delay values
if (h == hDelay)
......@@ -151,8 +150,6 @@ int main(int argc, char** argv)
if (h == hReport || h == hEnd)
{
reports++;
cerr << hmax << endl;
hmax = 0.0;
print();
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment