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

Fixed error in calculation of the initial conditions. Changed the matrix...

Fixed error in calculation of the initial conditions. Changed the matrix calculates to clearly indicate the input is a current and not a voltage. Add a max step size argument to the cvode simulator.
parent 3def0974
......@@ -32,9 +32,9 @@ void AdmittanceNetwork::add_self(int i, Complex y)
impl->incr(i,i,y);
}
void AdmittanceNetwork::solve(const Complex* Vsrc, Complex* V)
void AdmittanceNetwork::solve(const Complex* Isrc, Complex* V)
{
impl->solve(Vsrc,V);
impl->solve(Isrc,V);
}
......@@ -26,7 +26,7 @@ class AdmittanceNetwork
// Add self-impedence at node i
void add_self(int i, Complex y);
// Solve for voltage given a current
void solve(const Complex* Vsrc, Complex* V);
void solve(const Complex* Isrc, Complex* V);
// Get matrix entry
Complex get(int row, int col) { return impl->get(row,col); }
// Destructor
......
......@@ -14,7 +14,7 @@ ComplexMatrix_LAPACK::ComplexMatrix_LAPACK(int N):
N(N),
Y(new double[2*N*N]),
LUY(new double[2*N*N]),
v_lapack(new double[2*N]),
iv_lapack(new double[2*N]),
piv_lapack(new int[N])
{
for (int i = 0; i < 2*N*N; i++)
......@@ -49,7 +49,7 @@ Complex ComplexMatrix_LAPACK::get(int i, int j)
return Complex(Y[real_part],Y[imag_part]);
}
void ComplexMatrix_LAPACK::solve(const Complex* Vsrc, Complex* V)
void ComplexMatrix_LAPACK::solve(const Complex* Isrc, Complex* V)
{
int i,info,SIZE=N,lda=N,ldb=N,nrhs=1;
char type = 'N';
......@@ -65,11 +65,11 @@ void ComplexMatrix_LAPACK::solve(const Complex* Vsrc, Complex* V)
// Put the current into LAPACK form
for (i = 0; i < N; i++)
{
v_lapack[2*i] = real(Vsrc[i]);
v_lapack[2*i+1] = imag(Vsrc[i]);
iv_lapack[2*i] = real(Isrc[i]);
iv_lapack[2*i+1] = imag(Isrc[i]);
}
// Solve for the voltage
zgetrs_(&type,&SIZE,&nrhs,LUY,&lda,piv_lapack,v_lapack,&ldb,&info);
zgetrs_(&type,&SIZE,&nrhs,LUY,&lda,piv_lapack,iv_lapack,&ldb,&info);
if (info != 0)
{
cerr << "zgetrs_ returned " << info << endl;
......@@ -79,7 +79,7 @@ void ComplexMatrix_LAPACK::solve(const Complex* Vsrc, Complex* V)
// Get the current into the complex return vector
for (i = 0; i < N; i++)
{
V[i] = Complex(v_lapack[2*i],v_lapack[2*i+1]);
V[i] = Complex(iv_lapack[2*i],iv_lapack[2*i+1]);
}
}
......@@ -87,7 +87,7 @@ ComplexMatrix_LAPACK::~ComplexMatrix_LAPACK()
{
delete [] Y;
delete [] LUY;
delete [] v_lapack;
delete [] iv_lapack;
delete [] piv_lapack;
}
......@@ -14,7 +14,7 @@ class ComplexMatrix_LAPACK:
public:
// Create an N bus network
ComplexMatrix_LAPACK(int N);
void solve(const Complex* Vsrc, Complex* V);
void solve(const Complex* Isrc, Complex* V);
// Get the value at entry i,j
Complex get(int row, int col);
void set(int row, int col, Complex y);
......@@ -26,7 +26,7 @@ class ComplexMatrix_LAPACK:
const int N;
double *Y, *LUY;
// Scratch space for the solvers
double *v_lapack;
double *iv_lapack;
int* piv_lapack;
};
......
......@@ -146,45 +146,36 @@ void ElectricalData::initialize_circuit()
}
// Find initial generator currents and voltages
AdmittanceNetwork* Y = new AdmittanceNetwork(new ComplexMatrix_LAPACK(nodes.size()));
Complex *Vsrc = new Complex[nodes.size()];
Complex *Isrc = new Complex[nodes.size()];
Complex *Vnode = new Complex[nodes.size()];
// Node to ground currents or source voltages for generators
// Injected currents at generators
for (unsigned i = 0; i < nodes.size(); i++)
{
if (nodes[i].genr.genr_mw != 0.0 || nodes[i].genr.genr_mvar != 0.0)
{
Vsrc[i] = polar(nodes[i].vamp,nodes[i].theta);
Y->add_self(i,1.0);
Complex S = Complex(nodes[i].genr.genr_mw,nodes[i].genr.genr_mvar);
Complex V = polar(nodes[i].vamp,nodes[i].theta);
Isrc[i] = conj(S/V);
}
else
{
Vsrc[i] = Complex(0,0);
Y->add_self(i,nodes[i].y+nodes[i].B);
}
Isrc[i] = Complex(0,0);
Y->add_self(i,nodes[i].y+nodes[i].B);
}
// Node to node currents
for (auto line: lines)
{
int to_idx = node_to_idx[line.to];
int from_idx = node_to_idx[line.from];
bool to_genr = nodes[to_idx].genr.genr_mw != 0.0 || nodes[to_idx].genr.genr_mvar != 0.0;
bool from_genr = nodes[from_idx].genr.genr_mw != 0.0 || nodes[from_idx].genr.genr_mvar != 0.0;
if (!to_genr)
{
if (line.turns == 0.0)
Y->add_line(to_idx,from_idx,line.y);
else
Y->add_line(to_idx,from_idx,line.y,line.turns,0.0,line.tside);
}
if (!from_genr)
{
if (line.turns == 0.0)
Y->add_line(from_idx,to_idx,line.y);
else
Y->add_line(from_idx,to_idx,line.y,line.turns,0.0,line.tside);
}
if (line.turns == 0.0)
Y->add_line(to_idx,from_idx,line.y);
else
Y->add_line(to_idx,from_idx,line.y,line.turns,0.0,line.tside);
if (line.turns == 0.0)
Y->add_line(from_idx,to_idx,line.y);
else
Y->add_line(from_idx,to_idx,line.y,line.turns,0.0,line.tside);
}
Y->solve(Vsrc,Vnode);
Y->solve(Isrc,Vnode);
// Node initial conditions
for (auto entry: node_to_idx)
{
......@@ -236,24 +227,12 @@ void ElectricalData::initialize_circuit()
if (nodes[entry.second].genr.genr_mw == 0.0 && nodes[entry.second].genr.genr_mvar == 0.0)
continue;
configureGenerator(nodes[entry.second].genr);
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)
{
Igenr += line.I;
}
else if (line.to == nodes[entry.second].ID)
{
Igenr -= line.I;
}
}
Complex Igenr = Isrc[entry.second];
Complex S = Vnode[entry.second]*conj(Igenr);
nodes[entry.second].genr.S = S;
nodes[entry.second].genr.genr_mw = real(S);
nodes[entry.second].genr.genr_mvar = imag(S);
nodes[entry.second].genr.Pm0 = 0.0;
nodes[entry.second].genr.Rg = 0.0;
// Critically damp the generator
Complex X(nodes[entry.second].genr.Rg,nodes[entry.second].genr.Lg*2.0*pi*freq);
Complex Vgenr = Igenr*X+nodes[entry.second].V;
......@@ -267,6 +246,6 @@ void ElectricalData::initialize_circuit()
nodes[entry.second].genr.Etheta = arg(Vgenr);
}
delete Y;
delete [] Vsrc;
delete [] Isrc;
delete [] Vnode;
}
......@@ -45,7 +45,7 @@ void print(double t)
int main(int argc, char** argv)
{
double t = 0, cint = 1E-4, tend = 10.0, tReport = cint;
double t = 0, hmax = 1E-5, cint = 1E-4, tend = 10.0, tReport = cint;
N_Vector y, abstol;
SUNMatrix A;
SUNLinearSolver LS;
......@@ -57,9 +57,12 @@ int main(int argc, char** argv)
{
cint = atof(argv[++i]);
tReport = cint;
hmax = ::min(cint,hmax);
}
else if (strcmp("tend",argv[i]) == 0)
tend = atof(argv[++i]);
else if (strcmp("hmax",argv[i]) == 0)
hmax = atof(argv[++i]);
}
if (cint < 0.0 || tend <= 0.0)
{
......@@ -87,7 +90,7 @@ int main(int argc, char** argv)
}
// Setup cvode
cvode_mem = CVodeCreate(CV_BDF,CV_NEWTON);
CVodeSetMaxStep(cvode_mem,cint);
CVodeSetMaxStep(cvode_mem,hmax);
CVodeInit(cvode_mem,f,0.0,y);
CVodeSVtolerances(cvode_mem,tolerance,abstol);
A = SUNDenseMatrix(numStateVars,numStateVars);
......
......@@ -5,13 +5,13 @@ void event(double* const dq, double* const q, double* const alg, double tNow)
{
// Lose a generator
q[ig_2_a] = 0;
/* 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;
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;
......
......@@ -176,20 +176,12 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
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";
result += "der("+w+")=("+Pm+"-"+Pe+")/"+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"; */
/* result += "der("+Ef+")=0.0;\n";
result += "der("+generatorAngle(busName)+")=0;\n";
result += "der("+w+")=0;\n";
result += "der("+Pm+")=0;\n";
result += "der("+C+")=0;\n";
result += "der("+Pref+")=0;\n"; */
result += "der("+Pref+")=0;\n";*/
return result;
}
......
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