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

Add working voltage regulation at the generator busses

parent abec30c4
......@@ -32,7 +32,7 @@ test1: app
#./grid2ode ../data/test3.txt
#./grid2ode ../data/test2.txt
./m2c circuit.mo
#${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp euler.cpp event.cpp
#${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp euler.cpp event.cpp guess.cpp
${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp cvode.cpp event.cpp guess.cpp -lsundials_cvode
app: ${OBJS}
......
......@@ -77,7 +77,7 @@ int main(int argc, char** argv)
dq = new double[numStateVars];
q = new double[numStateVars];
// Calculate initial values
init_state(q);
init_state(q,alg);
for (int i = 0; i < numStateVars; i++)
{
NV_Ith_S(y,i) = q[i];
......
......@@ -12,7 +12,7 @@ using namespace std;
// Current simulation time
static double tNow = 0.0;
// Maximum quanta of time or state
static const double step = 1E-6;
static const double step = 1E-4;
// Simulation ending time
static double tEnd = 0.0;
// Reporting interval
......@@ -47,7 +47,7 @@ int main(int argc, char** argv)
dq = new double[numStateVars];
q = new double[numStateVars];
// Calculate initial values
init_state(q);
init_state(q,alg);
event(dq,q,tNow);
func(dq,q,alg,0.0);
event(dq,q,tNow);
......
#include <cmath>
#include <iostream>
#include <cassert>
#include <map>
#include <cfloat>
#include <cstdlib>
using namespace std;
static const double pi = 3.141592653589793;
......@@ -17,7 +17,7 @@ static void Jinv(double JJ[2][2], double V, double theta, double wt, double va,
JJ[1][1] = dfa_dV;
JJ[0][1] = -dfa_dtheta;
JJ[1][0] = -dfb_dV;
double det = JJ[0][0]*JJ[1][1]-JJ[0][1]*JJ[1][1];
double det = JJ[0][0]*JJ[1][1]-JJ[0][1]*JJ[1][0];
JJ[0][0] /= det;
JJ[0][1] /= det;
JJ[1][0] /= det;
......@@ -30,13 +30,26 @@ static void f(double* ff, double V, double theta, double wt, double va, double v
ff[1] = V*cos(wt+theta+2.0*pi/3.0)-vb;
}
double guessV(double V, double theta, double w, double t, double va, double vb)
static double uniform(double a, double b)
{
double u = double(rand())/double(RAND_MAX);
return u*a+(1.0-u)*b;
}
void guessV(double& V, double& theta, double w, double t, double va, double vb)
{
double V0 = V;
double t0 = theta;
int repeat = 0;
double JJ[2][2];
double ff[2];
double err = DBL_MAX;
double step = 1E-4;
do
repeat:
double step = 1E-3;
f(ff,V,theta,w*t,va,vb);
double err = fabs(ff[0])+fabs(ff[1]);
double Vprev = V;
double tprev = theta;
while (err > 1E-4 && step > 1E-8 && V > 0.0)
{
Jinv(JJ,V,theta,w*t,va,vb);
f(ff,V,theta,w*t,va,vb);
......@@ -44,12 +57,29 @@ double guessV(double V, double theta, double w, double t, double va, double vb)
assert(isfinite(ff[1]) && !isnan(ff[1]));
V -= step*(JJ[0][0]*ff[0]+JJ[1][0]*ff[1]);
theta -= step*(JJ[0][1]*ff[0]+JJ[1][1]*ff[1]);
f(ff,V,theta,w*t,va,vb);
double newErr = fabs(ff[0])+fabs(ff[1]);
if (newErr > err) step *= 0.5;
if (newErr > err)
{
V = Vprev;
theta = tprev;
step *= 0.1;
}
else
{
Vprev = V;
tprev = theta;
}
err = newErr;
}
while (err > 1E-4 && step > 1E-8);
if (err > 1E-4 || V < 0.0)
{
V = V0*uniform(0.5,1.5);
theta = t0*uniform(0.5,1.5);
repeat++;
goto repeat;
}
if (repeat > 0) cerr << "guess repeated " << repeat << " times" << endl;
assert(err <= 1E-4);
return fabs(V);
}
......@@ -10,7 +10,7 @@
enum
{
DELAY,
GUESS,
DER,
INIT,
NONE
......@@ -23,6 +23,7 @@ struct init_val_t
};
struct init_val_t to_init;
std::set<std::string> keywords;
std::set<std::string> alg_vars;
std::set<std::string> state_vars;
std::list<init_val_t> state_vars_init_vals;
......@@ -58,7 +59,7 @@ void lhs_to_rhs()
void identifier(const char* ID)
{
// Build list of variables
if (lhs && whichPass == 0)
if (whichPass == 0)
{
if (lhs)
{
......@@ -70,13 +71,16 @@ void identifier(const char* ID)
state_vars.insert(ID);
else if (whichFunc == INIT)
to_init.name = ID;
else alg_vars.insert(ID);
else if (strcmp(ID,"guessV") == 0)
whichFunc = GUESS;
else if (whichFunc != GUESS)
alg_vars.insert(ID);
}
}
// Generate code
else if (whichPass == 1)
{
if (strcmp(ID,"initialize") == 0 || strcmp(ID,"fix") == 0)
if (strcmp(ID,"initialize") == 0)
skip = true;
if (skip) return;
if (newLine)
......@@ -84,14 +88,8 @@ void identifier(const char* ID)
func << " ";
newLine = false;
}
// First argument of delay uses variable identifier
if (whichFunc == DELAY)
{
whichFunc = NONE;
func << ID;
}
// Otherwise use value for a variable
else if (state_vars.find(ID) != state_vars.end())
if (state_vars.find(ID) != state_vars.end())
{
if (lhs) func << "d";
func << "q[" << ID << "]";
......@@ -99,14 +97,23 @@ void identifier(const char* ID)
else if (alg_vars.find(ID) != alg_vars.end())
{
func << "alg[" << ID << "]";
}
}
else if (strcmp(ID,"guessV") == 0)
{
lhs = false;
func << ID;
whichFunc = GUESS;
}
// Otherwise its a function call. der only exists on lhs
// to identify state variables
else if (strcmp(ID,"der") != 0)
else if (!lhs && keywords.find(ID) != keywords.end())
{
func << ID;
if (strcmp(ID,"delay") == 0)
whichFunc = DELAY;
}
else if (strcmp(ID,"der") != 0)
{
alg_vars.insert(ID);
func << "alg[" << ID << "]";
}
}
}
......@@ -183,6 +190,11 @@ void make_gnuplot(const char* filename, const std::list<plot_var>& vars)
int main(int argc, char* argv[])
{
keywords.insert("initialize");
keywords.insert("time");
keywords.insert("cos");
keywords.insert("pi");
std::list<plot_var> power;
std::list<plot_var> speed;
std::list<plot_var> volts;
......@@ -232,7 +244,7 @@ int main(int argc, char* argv[])
volts.push_back(plot_var(alg_var,var_idx));
else if (alg_var.find("pe") != alg_var.npos)
power.push_back(plot_var(alg_var,var_idx));
else if (alg_var.find("Amp") != alg_var.npos)
else if (alg_var.find("vt") != alg_var.npos)
volt_amp_est.push_back(plot_var(alg_var,var_idx));
vars_out << "static const int " << alg_var << " = " << (k++) << "; // " << (var_idx++) << std::endl;
}
......@@ -243,10 +255,15 @@ int main(int argc, char* argv[])
func_out << func.str() << std::endl;
func_out << "}" << std::endl;
func_out << "void init_state(double* const q)";
func_out << "void init_state(double* const q, double* const alg)";
func_out << std::endl << "{" << std::endl;
for (auto val: state_vars_init_vals)
func_out << " q[" << val.name << "] = " << std::setprecision(20) << val.value << ";" << std::endl;
{
if (state_vars.find(val.name) != state_vars.end())
func_out << " q[" << val.name << "] = " << std::setprecision(20) << val.value << ";" << std::endl;
else
func_out << " alg[" << val.name << "] = " << std::setprecision(20) << val.value << ";" << std::endl;
}
func_out << "}" << std::endl << std::endl;
func_out.close();
vars_out << std::endl << "#endif" << std::endl;
......
......@@ -14,30 +14,29 @@ static const double pi = 3.14159265358979;
static double sys_freq_hz = 60.0;
// Enumeration of system phases
static const char* const phases[3] = { "a", "b", "c"};
// Set of voltage amplitude measurement variables
static set<string> voltageAmplitudeEstimators;
// Default initial values for voltages and currents
// Initial value data for a state variable
struct state_var_init_t
// Initial value data for variables
struct var_init_t
{
// Variable name
string name;
// Initial value
double init_value;
state_var_init_t(string name, double init_value):
var_init_t(string name, double init_value):
name(name),init_value(init_value){}
state_var_init_t(){}
var_init_t(){}
// These get stored in a set that is sorted by name
// to avoid repetitions
bool operator<(const state_var_init_t& other) const
bool operator<(const var_init_t& other) const
{
return name < other.name;
}
};
// Initial variable values for all state variables
set<state_var_init_t> stateVars;
set<var_init_t> stateVars;
set<var_init_t> algVars;
// Number to string
string number(double val)
......@@ -140,11 +139,12 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
string C = "c_"+busName;
string Pref = "pref_"+busName;
string Vguess = "vt_"+busName;
string ThetaGuess = "theta_t_"+busName;
string powerDraw;
// Power demand at the terminal is sum of demand for each phase
for (auto phase: phases)
{
stateVars.insert(state_var_init_t(generatorCurrent(busName,phase),data.genr.i[i++]));
stateVars.insert(var_init_t(generatorCurrent(busName,phase),data.genr.i[i++]));
result += generatorSourceVoltage(busName,phase);
result += "der("+generatorCurrent(busName,phase)+")=("+
generatorVoltage(busName,phase)+"-"+voltage(busName,phase)+"-"+
......@@ -155,17 +155,17 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
powerDraw += "+"+generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
}
result += Pe+"="+powerDraw+";\n";
stateVars.insert(state_var_init_t(generatorAngle(busName),data.genr.Etheta));
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));
stateVars.insert(state_var_init_t(Pref,data.genr.Pm0));
stateVars.insert(state_var_init_t(Ef,data.genr.Eamp));
/* result += Vguess+"="+"guessV("+Ef+","+generatorAngle(busName)+","+
number(120.0*pi)+",time,"+voltage(busName,phases[0])+","+voltage(busName,phases[1])+");\n"; */
//result += "der("+Ef+")="+number(data.genr.G)+"*("+number(data.vamp)+"-"+Vguess+");\n";
result += "der("+Ef+")=0.0;\n";
stateVars.insert(var_init_t(generatorAngle(busName),data.genr.Etheta));
stateVars.insert(var_init_t(w,data.genr.w0));
stateVars.insert(var_init_t(Pm,data.genr.Pm0));
stateVars.insert(var_init_t(C,data.genr.Pm0));
stateVars.insert(var_init_t(Pref,data.genr.Pm0));
stateVars.insert(var_init_t(Ef,data.genr.Eamp));
algVars.insert(var_init_t(Vguess,data.vamp));
algVars.insert(var_init_t(ThetaGuess,data.theta));
result += "guessV("+Vguess+","+ThetaGuess+","+
number(120.0*pi)+",time,"+voltage(busName,phases[0])+","+voltage(busName,phases[1])+");\n";
result += "der("+Ef+")="+number(data.genr.G)+"*("+number(data.vamp)+"-"+Vguess+");\n";
result += "der("+generatorAngle(busName)+")="+w+";\n";
result += "der("+w+")=("+Pm+"-"+Pe+"-"+number(data.genr.D)+"*"+w+")/"+Mstr+";\n";
result += "der("+C+")="+number(data.genr.Tspd1)+"*("+Pref+"-"+number(data.genr.R)+
......@@ -249,18 +249,18 @@ string createBus(
isum += "+"+generatorCurrent(busName,phase);
isum += ")";
// Voltage is a state variable
stateVars.insert(state_var_init_t(v,bus.v[phase_num]));
stateVars.insert(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)
{
stateVars.insert(state_var_init_t(iload,bus.iload[phase_num]));
stateVars.insert(var_init_t(iload,bus.iload[phase_num]));
result += string("der(")+iload+")=("+v+"-"+iload+"*"+number(bus.R)+")/"+number(bus.L)+";\n";
}
// Add load voltage as a state variable
else if (bus.C > 0.0)
{
stateVars.insert(state_var_init_t(vload,bus.vload[phase_num]));
stateVars.insert(var_init_t(vload,bus.vload[phase_num]));
result += "der("+vload+")=(("+v+"-"+vload+")/"+number(bus.R)+")/"+number(bus.C)+";\n";
}
// Next phase
......@@ -280,7 +280,7 @@ string createLine(ElectricalData::line_t& line)
string va = voltage(bus_name(line.from),phase);
string vb = voltage(bus_name(line.to),phase);
string i = current(line_name(line),phase);
stateVars.insert(state_var_init_t(i,line.i[phase_num]));
stateVars.insert(var_init_t(i,line.i[phase_num]));
result += string("der(")+i+string(")=(")+va+"-"+vb;
if (line.R > 0.0)
result += string("-")+number(line.R)+string("*")+i;
......@@ -307,6 +307,10 @@ void make_model(ElectricalData& data)
{
fout << "initialize " << stateVar.name << "=" << number(stateVar.init_value) << ";" << endl;
}
for (auto algVar: algVars)
{
fout << "initialize " << algVar.name << "=" << number(algVar.init_value) << ";" << endl;
}
fout.close();
}
......
......@@ -4,7 +4,7 @@
extern const int numStateVars;
extern const int numAlgVars;
void init_state(double* const q);
void init_state(double* const q, double* const alg);
void func(double* const dq, const double* const q, double* const alg, const double t);
void event(double* const dq, double* const q, const double t);
/**
......@@ -14,6 +14,6 @@ void event(double* const dq, double* const q, const double t);
* is the present frequency of the generator. The voltages va
* and vb are the present voltages at the bus to estimate.
*/
double guessV(double V, double theta, double w, double t, double va, double vb);
void guessV(double& V, double& theta, double w, double t, double va, double vb);
#endif
......@@ -23,7 +23,8 @@ int main(int argc, char** argv)
double vb = V0*cos(w*t+theta+2.0*pi/3.0);
double V = uniform(0.5,1.5);
theta = uniform(-pi,pi);
cout << guessV(V,theta,w,t,va,vb) << " =? " << V0 << endl;
guessV(V,theta,w,t,va,vb);
cout << V << " =? " << V0 << endl;
return 0;
}
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