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

Generates working 14 and 1 bus models. The 14 bus model appears to suffer...

Generates working 14 and 1 bus models. The 14 bus model appears to suffer voltage collapse (?) due to sagging output current. I think this is caused by the rather ad hoc power output control law that the model uses.
parent 4ffebb38
16/11/09 THYME TEST CASE 1.0 2009 W ORNL Test Case
BUS DATA FOLLOWS 2 ITEMS
1 Generator 1 1 2 1.011 0.5100 0.0 0.0 1.0101 0.1101 1.0 1.0 10.0 10.0 0.0 0.0 0
1 Generator 1 1 2 1.011 0.5100 0.0 0.0 3.45 0.1101 1.0 1.0 10.0 10.0 0.0 0.0 0
2 Load 1 1 2 1.000 0.0000 1.0 0.1 0.0000 0.0000 1.0 1.0 10.0 10.0 0.0 0.0 0
-999
BRANCH DATA FOLLOWS 1 ITEMS
2 1 1 1 1 0 0.01 0.01 0.0 0 0 0 0 0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
1 2 1 1 1 0 0.01 0.01 0.0 0 0 0 0 2 1.0 0.0 0.0 0.0 0.0 0.0 0.0
-999
LOSS ZONES FOLLOWS 1 ITEMS
1 ORNL Test Case
......
......@@ -16,7 +16,7 @@ INCLUDE = -I/usr/include/eigen3
OBJS = \
ElectricalData.o \
IEEE_CDF_Data.o \
main.o
main.o
test1: app
flex m2c.lex
......@@ -24,11 +24,11 @@ test1: app
./grid2ode ../data/ieee14cdf.txt
#./grid2ode ../data/test1.txt
./m2c circuit.mo
${CXX} ${CFLAGS} ${INCLUDE} func.cpp sim.cpp
${CXX} ${CFLAGS} ${INCLUDE} func.cpp sim.cpp poly.cpp
app: ${OBJS}
${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}
clean:
rm -f *.o a.out grid2ode func.cpp soln.txt m2c lex.yy.c
rm -f *.o a.out grid2ode func.cpp soln.txt m2c lex.yy.c *.gnu
va_1=cos(376.991*w_1*time+theta_1);
vb_1=cos(376.991*w_1*time-2.0*pi/3.0+theta_1);
vc_1=cos(376.991*w_1*time-4.0*pi/3.0+theta_1);
der(theta_1)=(w_1-1.0);
der(w_1)=(pm_1-((va_1*(ia_0_2_1))+(vb_1*(ib_0_2_1))+(vc_1*(ic_0_2_1))))/753.982;
der(pm_1)=376.991*(c_1-pm_1);
der(c_1)=3769.91*(pref_1-7539.82*(w_1-1.0)-c_1);
der(pref_1)=20*(1.0-w_1);
der(va_2)=(-ia_0_2_1)/2.65258e-05;
der(vb_2)=(-ib_0_2_1)/2.65258e-05;
der(vc_2)=(-ic_0_2_1)/2.65258e-05;
der(ia_0_2_1)=(va_2-va_1-0.01*ia_0_2_1)/2.65258e-05;
der(ib_0_2_1)=(vb_2-vb_1-0.01*ib_0_2_1)/2.65258e-05;
der(ic_0_2_1)=(vc_2-vc_1-0.01*ic_0_2_1)/2.65258e-05;
fix c_1=1.0101;
initialize ia_0_2_1=0;
initialize ib_0_2_1=0;
initialize ic_0_2_1=0;
fix pm_1=1.0101;
fix pref_1=1.0101;
fix theta_1=0.000155276;
initialize va_2=0;
initialize vb_2=0;
initialize vc_2=0;
fix w_1=1;
......@@ -155,8 +155,37 @@ OP [-+*/(),]
%%
struct plot_var
{
plot_var(std::string name, int ID):name(name),ID(ID){}
std::string name;
int ID;
};
void make_gnuplot(const char* filename, const std::list<plot_var>& vars)
{
int k = 0;
std::ofstream func_out(filename);
func_out << "plot ";
for (auto var: vars)
{
if (k == 0) func_out << "'soln.txt' ";
else func_out << ",'' ";
func_out << "using 1:" << (var.ID) << " with lines";
k++;
}
func_out << std::endl;
func_out.close();
}
int main(int argc, char* argv[])
{
std::list<plot_var> power;
std::list<plot_var> speed;
std::list<plot_var> volts;
std::list<plot_var> angle;
for (whichPass = 0; whichPass < 2; whichPass++)
{
FILE* fp = fopen(argv[1],"r");
......@@ -173,11 +202,29 @@ int main(int argc, char* argv[])
func_out << "using namespace std;" << std::endl << std::endl;
func_out << "static const double pi = " << 3.141592653589793 << ";" << std::endl << std::endl;
for (auto state_var: state_vars)
{
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));
else if (state_var.find("pm") != state_var.npos)
power.push_back(plot_var(state_var,var_idx));
else if (state_var.find("theta") != state_var.npos)
angle.push_back(plot_var(state_var,var_idx));
func_out << "static const int " << state_var << " = " << (k++) << "; // " << (var_idx++) << std::endl;
}
func_out << std::endl << "const int numStateVars = " << state_vars.size() << ";" << std::endl << std::endl;
k = 0;
for (auto alg_var: alg_vars)
{
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));
else if (alg_var.find("pe") != alg_var.npos)
power.push_back(plot_var(alg_var,var_idx));
func_out << "static const int " << alg_var << " = " << (k++) << "; // " << (var_idx++) << std::endl;
}
func_out << std::endl << "const int numAlgVars = " << alg_vars.size() << ";" << std::endl << std::endl;
func_out << "void func(double* const dq, const double* const q, double* const alg, const double time)";
......@@ -192,5 +239,10 @@ int main(int argc, char* argv[])
func_out << "}" << std::endl << std::endl;
func_out.close();
make_gnuplot("volts.gnu",volts);
make_gnuplot("power.gnu",power);
make_gnuplot("speed.gnu",speed);
make_gnuplot("angle.gnu",angle);
return EXIT_SUCCESS;
}
......@@ -211,7 +211,7 @@ string createGenerator(string busName, ElectricalData::genr_t data, double cap)
stateVars.insert(state_var_init_t(w,data.w0));
stateVars.insert(state_var_init_t(Pm,data.Pm0));
// Comment these out and see the next comment block to keep the generators still
result += "der("+generatorAngle(busName)+")=("+w+"-1.0);\n";
result += "der("+generatorAngle(busName)+")="+number(Freq)+"*("+w+"-1.0);\n";
result += string("der(")+w+string(")=(")+Pm+
string("-"+Pe+")/")+Hstr+string(";\n");
result += string("der(")+Pm+string(")="+Tp2str+"*("+
......@@ -233,7 +233,7 @@ string createBus(
ElectricalData::bus_data_t bus = data.getBusInfo(nodeID);
string busName = bus_name(bus.ID);
// Small, minimum value for the capacitance relative to ground
static const double Cmin = 1E-5;
static const double Cmin = 1E-3;
Complex Z = 1.0/data.getAdmittance(nodeID);
double C = 0.0, R = real(Z), Rimg = imag(Z);
string result = "";
......
#include "sim.h"
#include "poly.h"
#include <cfloat>
#include <algorithm>
#include <vector>
......@@ -28,6 +29,7 @@ static const int dSize = 1+int(maxDelay/dint);
static double **dQueue;
// Model variable values
static double *q, *dq, *alg;
static poly* p;
void print()
{
......@@ -91,6 +93,7 @@ int main(int argc, char** argv)
alg = new double[numAlgVars];
dq = new double[numStateVars];
q = new double[numStateVars];
p = new poly[numStateVars];
dQueue = new double*[numStateVars];
for (int i = 0; i < numStateVars; i++)
{
......@@ -101,6 +104,8 @@ int main(int argc, char** argv)
// Calculate initial values
init_state(q);
func(dq,q,alg,0.0);
for (int i = 0; i < numStateVars; i++)
p[i].init(dq[i]);
// Record (the very dull) initial state
print();
// Simulate until requested end of interval
......@@ -113,15 +118,17 @@ int main(int argc, char** argv)
hDelay = dint*double(delays)-tNow;
hEnd = tEnd-tNow;
h = ::min(hDelay,::min(hReport,hEnd));
#pragma omp parallel for reduction(min:h)
for (int i = 0; i < numStateVars; i++)
{
if (fabs(dq[i]) > 0.0)
h = min(h,quant/fabs(dq[i]));
p[i].append(dq[i]);
h = min(h,p[i].time_until(quant));
}
// Update the state variables
#pragma omp parallel for
for (int i = 0; i < numStateVars; i++)
{
q[i] += dq[i]*h;
q[i] += p[i].extrapolate(h);
}
// Advance the clock
tNow += h;
......@@ -149,6 +156,7 @@ int main(int argc, char** argv)
delete [] dq;
delete [] q;
delete [] alg;
delete [] p;
for (int i = 0; i < numStateVars; i++)
delete [] dQueue[i];
delete [] dQueue;
......
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