Commit 3def0974 authored by Jim Nutaro's avatar Jim Nutaro
Browse files

Added switch to effect loss of a generator. Fixed output so that is was at a...

Added switch to effect loss of a generator. Fixed output so that is was at a regular interval for easy use in an FFT
parent bdc17553
......@@ -18,18 +18,16 @@ using namespace std;
// Numerical tolerace. Picked these arbitrarily.
static const double tolerance = 1E-6;
static const double maxStep = 1E-3;
// Model variable values
static double *q, *dq, *alg;
static double hmin = 0.0;
static double *q, *qp, *dq, *alg;
static int f(realtype t, N_Vector y, N_Vector ydot, void* user_data)
{
for (int i = 0; i < numStateVars; i++)
q[i] = NV_Ith_S(y,i);
event(dq,q,t);
event(dq,q,alg,t);
func(dq,q,alg,t);
event(dq,q,t);
event(dq,q,alg,t);
for (int i = 0; i < numStateVars; i++)
NV_Ith_S(ydot,i) = dq[i];
return 0;
......@@ -39,15 +37,15 @@ void print(double t)
{
cout << t;
for (int i = 0; i < numStateVars; i++)
cout << " " << q[i];
cout << " " << qp[i];
for (int i = 0; i < numAlgVars; i++)
cout << " " << alg[i];
cout << " " << hmin << endl;
cout << endl;
}
int main(int argc, char** argv)
{
double t = 0, cint = 0.0, tend = 0.0, tReport = 0.0;
double t = 0, cint = 1E-4, tend = 10.0, tReport = cint;
N_Vector y, abstol;
SUNMatrix A;
SUNLinearSolver LS;
......@@ -56,7 +54,10 @@ int main(int argc, char** argv)
for (int i = 1; i < argc; i++)
{
if (strcmp("cint",argv[i]) == 0)
{
cint = atof(argv[++i]);
tReport = cint;
}
else if (strcmp("tend",argv[i]) == 0)
tend = atof(argv[++i]);
}
......@@ -76,48 +77,52 @@ int main(int argc, char** argv)
alg = new double[numAlgVars];
dq = new double[numStateVars];
q = new double[numStateVars];
qp = new double[numStateVars];
// Calculate initial values
init_state(q,alg);
for (int i = 0; i < numStateVars; i++)
{
NV_Ith_S(y,i) = q[i];
NV_Ith_S(y,i) = qp[i] = q[i];
NV_Ith_S(abstol,i) = tolerance;
}
// Setup cvode
cvode_mem = CVodeCreate(CV_BDF,CV_NEWTON);
CVodeSetMaxStep(cvode_mem,maxStep);
CVodeSetMaxStep(cvode_mem,cint);
CVodeInit(cvode_mem,f,0.0,y);
CVodeSVtolerances(cvode_mem,tolerance,abstol);
A = SUNDenseMatrix(numStateVars,numStateVars);
LS = SUNDenseLinearSolver(y,A);
CVDlsSetLinearSolver(cvode_mem,LS,A);
// Report initial state
event(dq,q,t);
event(dq,q,alg,t);
func(dq,q,alg,t);
event(dq,q,t);
event(dq,q,alg,t);
print(t);
hmin = DBL_MAX;
// Run simulation
while (t < tend)
{
double tPrev = t;
CVode(cvode_mem,t+maxStep,y,&t,CV_ONE_STEP);
hmin = ::min(hmin,t-tPrev);
if (t >= tReport+cint)
assert(t < tReport);
CVode(cvode_mem,t+cint,y,&t,CV_ONE_STEP);
if (t >= tReport)
{
tReport += cint;
double h = (tReport-tPrev)/(t-tPrev);
assert(h >= tReport-tPrev);
for (int i = 0; i < numStateVars; i++)
q[i] = NV_Ith_S(y,i);
event(dq,q,t);
func(dq,q,alg,t);
event(dq,q,t);
print(t);
hmin = DBL_MAX;
qp[i] = (1.0-h)*qp[i]+h*NV_Ith_S(y,i);
event(dq,qp,alg,tReport);
func(dq,qp,alg,tReport);
event(dq,qp,alg,tReport);
print(tReport);
tReport += cint;
}
for (int i = 0; i < numStateVars; i++)
qp[i] = NV_Ith_S(y,i);
}
// Cleanup and exit
delete [] dq;
delete [] q;
delete [] qp;
delete [] alg;
N_VDestroy(y);
......
......@@ -48,9 +48,9 @@ int main(int argc, char** argv)
q = new double[numStateVars];
// Calculate initial values
init_state(q,alg);
event(dq,q,tNow);
event(dq,q,alg,tNow);
func(dq,q,alg,0.0);
event(dq,q,tNow);
event(dq,q,alg,tNow);
// Record the initial state
print();
hmin = DBL_MAX;
......@@ -77,11 +77,11 @@ int main(int argc, char** argv)
// Advance the clock
tNow += h;
// Apply events to state variable values
event(dq,q,tNow);
event(dq,q,alg,tNow);
// Calculate new variable values for tNow
func(dq,q,alg,tNow);
// Apply derivatives values
event(dq,q,tNow);
event(dq,q,alg,tNow);
// Report if its time or at the end of the simulation
if (h == hReport || h == hEnd)
{
......
#include "sim.h"
#include "vars.h"
void event(double* const dq, double* const q, double tNow)
void event(double* const dq, double* const q, double* const alg, double tNow)
{
// Lose a generator
q[ig_2_a] = 0;
q[ig_2_b] = 0;
q[ig_2_c] = 0;
dq[ig_2_a] = 0;
dq[ig_2_b] = 0;
dq[ig_2_c] = 0;
alg[switch_2] = 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_3_2_4] = 0.0;
/* q[ia_3_2_4] = 0.0;
q[ib_3_2_4] = 0.0;
q[ic_3_2_4] = 0.0;
dq[ia_3_2_4] = 0.0;
dq[ib_3_2_4] = 0.0;
dq[ic_3_2_4] = 0.0;
dq[ic_3_2_4] = 0.0; */
// Drop a line in the ieee 188 bus model
/* q[ia_102_66_67] = 0.0;
......
......@@ -36,8 +36,9 @@ static double uniform(double a, double b)
return u*a+(1.0-u)*b;
}
void 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, double on)
{
if (on <= 0.0) return;
double V0 = V;
double t0 = theta;
int repeat = 0;
......
......@@ -223,7 +223,7 @@ int main(int argc, char* argv[])
vars_out << "// State variables " << std::endl << std::endl;
for (auto state_var: state_vars)
{
if (state_var.find('w') != state_var.npos)
if (state_var.find("w_") != state_var.npos)
speed.push_back(plot_var(state_var,var_idx));
else if (state_var.find("va") != state_var.npos)
volts.push_back(plot_var(state_var,var_idx));
......@@ -238,7 +238,7 @@ int main(int argc, char* argv[])
k = 0;
for (auto alg_var: alg_vars)
{
if (alg_var.find('w') != alg_var.npos)
if (alg_var.find("w_") != alg_var.npos)
speed.push_back(plot_var(alg_var,var_idx));
else if (alg_var.find("va") != alg_var.npos)
volts.push_back(plot_var(alg_var,var_idx));
......
......@@ -140,6 +140,7 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
string Pref = "pref_"+busName;
string Vguess = "vt_"+busName;
string ThetaGuess = "theta_t_"+busName;
string GenrSwitch = "switch_"+busName;
string powerDraw;
// Power demand at the terminal is sum of demand for each phase
for (auto phase: phases)
......@@ -163,16 +164,17 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
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));
algVars.insert(var_init_t(GenrSwitch,1.0));
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)+
"*"+w+"-"+C+");\n";
result += string("der(")+Pm+")="+
number(data.genr.Tspd2)+"*("+C+"-"+Pm+");\n";
result += "der("+Pref+")=-"+number(data.genr.Agc)+"*"+w+";\n";
number(120.0*pi)+",time,"+voltage(busName,phases[0])+","+voltage(busName,phases[1])+","+GenrSwitch+");\n";
result += "der("+Ef+")=("+number(data.genr.G)+"*("+number(data.vamp)+"-"+Vguess+"))*"+GenrSwitch+";\n";
result += "der("+generatorAngle(busName)+")=("+w+")*"+GenrSwitch+";\n";
result += "der("+w+")=(("+Pm+"-"+Pe+"-"+number(data.genr.D)+"*"+w+")/"+Mstr+")*"+GenrSwitch+";\n";
result += "der("+C+")=("+number(data.genr.Tspd1)+"*("+Pref+"-"+number(data.genr.R)+
"*"+w+"-"+C+"))*"+GenrSwitch+";\n";
result += string("der(")+Pm+")=("+
number(data.genr.Tspd2)+"*("+C+"-"+Pm+"))*"+GenrSwitch+";\n";
result += "der("+Pref+")=(-"+number(data.genr.Agc)+"*"+w+")*"+GenrSwitch+";\n";
// Comment these out and see the next comment block to keep the generators still
/* result += "der("+Ef+")="+number(data.genr.G)+"*("+Pm+"-"+Pe+"-"+number(data.genr.K)+"*"+Ef+");\n";
result += "der("+generatorAngle(busName)+")="+w+";\n";
......
......@@ -6,7 +6,7 @@ extern const int numAlgVars;
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);
void event(double* const dq, double* const q, double* const alg, const double t);
/**
* This estimates the amplitude of V at a node attached to a generator.
* The arguments are the initial guess for theta and V, which should
......@@ -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.
*/
void 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, double on);
#endif
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