Fixed error in calculting generator voltages. Added start on a voltage control for the generator.

src/ElectricalData.cpp

@@ -28,27 +28,27 @@ void ElectricalData::configureGenerator(genr_t& genr_data)
 	/**
 	 * Inertia estimates are calculated from data in Table 5.1 of
 	 * "Computer-Aided Power System Analysis" by Ramasamy Natarajan,
-	 * published by Marcel Dekker, Inc., c. 2002.
+	 * published by Marcel Dekker, Inc., c. 2002. Define:
 	 */
 	if (genr_data.genr_mw * SBase > 600.0) // Big nuke plant
 	{
-		genr_data.H = 3.5;
+		genr_data.M = 3.5*600.0/SBase;
 		genr_data.Lg = 0.26/(2.0*pi*freq);
 	}
 	else if (genr_data.genr_mw * SBase > 300.0) // Big fossil plant
 	{
-		genr_data.H = 4.0;
+		genr_data.M = 4.0*300.0/SBase;
 		genr_data.Lg = 0.25/(2.0*pi*freq);
 	}
 
-	else if (genr_data.genr_mw * SBase > 150.0) // Smaller steam plant
+	else if (genr_data.genr_mw * SBase > 120.0) // Smaller steam plant
 	{
-		genr_data.H = 6.0;
+		genr_data.M = 3.3*120.0/SBase;
 		genr_data.Lg = 0.26/(2.0*pi*freq);
 	}
 	else 
 	{
-		genr_data.H = 2.0;
+		genr_data.M = 6.0*50.0/SBase;
 		genr_data.Lg = 0.2/(2.0*pi*freq);
 	}
 	// Try to get close to an appropriate set of initial conditions
@@ -56,8 +56,8 @@ void ElectricalData::configureGenerator(genr_t& genr_data)
 	genr_data.Tspd1 = 10.0;
 	genr_data.Tspd2 = 1.0;
 	genr_data.R = 1.0/0.05;
-	genr_data.Agc = genr_data.R/10.0;
-	genr_data.K = 100.0;
+	genr_data.Agc = genr_data.R*0.1;
+	genr_data.D = 1.0;
 	genr_data.G = 400.0;
 	genr_data.Rg = 0.0;	
 }
@@ -249,6 +249,7 @@ void ElectricalData::initialize_circuit()
 			}
 		}
 		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;

src/ElectricalData.h

@@ -56,16 +56,17 @@ class ElectricalData
 		struct genr_t
 		{
 			/// Inertia constant
-			double H;
+			double M;
 			/// These are the speed control parameters
-			double Tspd1, Tspd2, R, Agc; // Speed gain
+			double Tspd1, Tspd2, R, Agc, D; // Speed gain
 			/// Reference output current amplitude and phase
 			double Eamp, Etheta, Lg, Rg;
 			double i[3];
 			/// Voltage control
-			double K, G;
+			double G;
 			/// Power for plant sizing information
 			double genr_mw, genr_mvar;
+			Complex S;
 			/// Initial speed
 			double w0;
 			/// Initial power

src/Makefile

@@ -19,6 +19,7 @@ OBJS = \
 	PSSE_Raw_Data.o \
 	AdmittanceNetwork.o \
 	ComplexMatrix_LAPACK.o \
+	guess.o \
 	main.o 
 
 test1: app 
@@ -32,7 +33,7 @@ test1: app
 	#./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 cvode.cpp event.cpp -lsundials_cvode
+	${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp cvode.cpp event.cpp guess.cpp -lsundials_cvode
 
 app: ${OBJS}
 	${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}

src/PSSE_Raw_Data.cpp

@@ -85,7 +85,7 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 			bus_data.genr.genr_mw = 0.0;
 			bus_data.genr.genr_mvar = 0.0;
 			// Add this bus only if it has a non-zero voltage
-			if (bus_data.v != 0)
+			if (bus_data.vamp != 0)
 			{
 				bus_id_to_idx[bus_data.ID] = nodes.size();
 				nodes.push_back(bus_data);
@@ -190,11 +190,9 @@ skip_line:
 			GETTOK(NULL);
 			GETTOK(NULL);
 			// Transformer data, ignore it
-			float turns = 0.0; // atof(GETTOK(NULL));
+			float turns = atof(GETTOK(NULL));
 			if (turns != 0.0)
-			{
 				line_data.tside = line_data.from;
-			}
 			else GETTOK(NULL); // Skip the phase shift field if there is no transformer
 			// Get shunts
 			Complex BI(atof(GETTOK(NULL)),atof(GETTOK(NULL)));

src/cvode.cpp

@@ -17,8 +17,8 @@
 using namespace std;
 
 // Numerical tolerace. Picked these arbitrarily.
-static const double tolerance = 1E-10;
-static const double maxStep = 1E-8;
+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;

src/event.cpp

@@ -11,5 +11,13 @@ void event(double* const dq, double* const q, double tNow)
 	dq[ia_3_2_4] = 0.0;
 	dq[ib_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;
+	q[ib_102_66_67] = 0.0;
+	q[ic_102_66_67] = 0.0;
+	dq[ia_102_66_67] = 0.0;
+	dq[ib_102_66_67] = 0.0;
+	dq[ic_102_66_67] = 0.0; */
 }
 

src/guess.cpp

+#include <cmath>
+#include <iostream>
+#include <cassert>
+#include <map>
+#include <cfloat>
+using namespace std;
+
+static const double pi = 3.141592653589793;
+
+static void Jinv(double JJ[2][2], double V, double theta, double wt, double va, double vb)
+{
+	double dfa_dV = cos(wt+theta);
+	double dfa_dtheta = -V*sin(wt+theta);
+	double dfb_dV = cos(wt+theta+2.0*pi/3.0);
+	double dfb_dtheta = -V*sin(wt+theta+2.0*pi/3.0);
+	JJ[0][0] = dfb_dtheta;
+	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];
+	JJ[0][0] /= det;
+	JJ[0][1] /= det;
+	JJ[1][0] /= det;
+	JJ[1][1] /= det;
+}
+
+static void f(double* ff, double V, double theta, double wt, double va, double vb)
+{
+	ff[0] = V*cos(wt+theta)-va;
+	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)
+{
+	double JJ[2][2];
+	double ff[2];
+	double err = DBL_MAX;
+	double step = 1E-4;
+	do
+	{
+		Jinv(JJ,V,theta,w*t,va,vb);
+		f(ff,V,theta,w*t,va,vb);
+		assert(isfinite(ff[0]) && !isnan(ff[0]));
+		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]);
+		double newErr = fabs(ff[0])+fabs(ff[1]);
+		if (newErr > err) step *= 0.5;
+		err = newErr;
+	}
+	while (err > 1E-4 && step > 1E-8);
+	assert(err <= 1E-4);
+	return fabs(V);
+}
+

src/main.cpp

@@ -111,20 +111,18 @@ string generatorExcitor(string busName)
 }
 
 // Generator source voltage
-string generatorSourceVoltage(string busName, string phase, ElectricalData::genr_t genr)
+string generatorSourceVoltage(string busName, string phase)
 {
 	string ef = generatorExcitor(busName);
 	string result;
 	string freq = number(2.0*pi*sys_freq_hz);
-	string w = omega(busName);
-	string g = number((2.0*pi*sys_freq_hz)/genr.Eamp);
-	string theta = generatorAngle(busName)+"+"+number(genr.Etheta);;
+	string theta = generatorAngle(busName);
 	if (phase == phases[0])
-		result += generatorVoltage(busName,phase)+"=("+number(genr.Eamp)+"+"+ef+")*cos(("+freq+"+"+w+")*time+"+theta+");\n";
+		result += generatorVoltage(busName,phase)+"="+ef+"*cos("+freq+"*time+"+theta+");\n";
 	else if (phase == phases[1])
-		result += generatorVoltage(busName,phase)+"=("+number(genr.Eamp)+"+"+ef+")*cos(("+freq+"+"+w+")*time+2.0*pi/3.0+"+theta+");\n";
+		result += generatorVoltage(busName,phase)+"="+ef+"*cos("+freq+"*time+2.0*pi/3.0+"+theta+");\n";
 	else
-		result += generatorVoltage(busName,phase)+"=("+number(genr.Eamp)+"+"+ef+")*cos(("+freq+"+"+w+")*time+4.0*pi/3.0+"+theta+");\n";
+		result += generatorVoltage(busName,phase)+"="+ef+"*cos("+freq+"*time+4.0*pi/3.0+"+theta+");\n";
 	return result;
 }
 
@@ -135,20 +133,19 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
 	string result;
 	string varName;
 	string Ef = generatorExcitor(busName);
-	string Hstr = number(data.genr.H);
-	string Pt = "pt_"+busName;
-	string wt = "wt_"+busName;
+	string Mstr = number(data.genr.M);
 	string w = omega(busName);
 	string Pm = mech_power(busName);
 	string Pe = "pe_"+busName;
 	string C = "c_"+busName;
 	string Pref = "pref_"+busName;
+	string Vguess = "vt_"+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++]));
-		result += generatorSourceVoltage(busName,phase,data.genr);
+		result += generatorSourceVoltage(busName,phase);
 		result += "der("+generatorCurrent(busName,phase)+")=("+
 			generatorVoltage(busName,phase)+"-"+voltage(busName,phase)+"-"+
 			generatorCurrent(busName,phase)+"*"+number(data.genr.Rg)+")/"+number(data.genr.Lg)+";\n";
@@ -158,28 +155,28 @@ 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),0.0));
+	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(wt,data.genr.w0));
 	stateVars.insert(state_var_init_t(Pm,data.genr.Pm0));
-	stateVars.insert(state_var_init_t(Pt,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,0.0));
+	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";
 	result += "der("+generatorAngle(busName)+")="+w+";\n";
-	result += "der("+w+")=("+Pt+"-"+Pe+")/"+Hstr+";\n";
-	result += "der("+wt+")=("+Pm+"-"+Pt+")/"+Hstr+";\n";
-	result += "der("+Pt+")=("+wt+"-"+w+"-1E-1*"+Pt+")/1E-5;\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)+
-		"*"+wt+"-"+C+");\n";
+		"*"+w+"-"+C+");\n";
 	result += string("der(")+Pm+")="+
 		number(data.genr.Tspd2)+"*("+C+"-"+Pm+");\n";
-	result += "der("+Pref+")=-"+number(data.genr.Agc)+"*"+wt+";\n";
+	result += "der("+Pref+")=-"+number(data.genr.Agc)+"*"+w+";\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+")/"+Hstr+";\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+")="+

src/sim.h

@@ -4,9 +4,16 @@
 extern const int numStateVars;
 extern const int numAlgVars;
 
-double delay(int var, double h, double t);
 void init_state(double* const q);
 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);
+/**
+ * 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
+ * be the same at the attached generator to be close. The frequency w
+ * 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);
 
 #endif

src/solve.cpp

+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+#include "sim.h"
+using namespace std;
+
+static const double pi = 3.14;
+
+static double uniform(double a, double b)
+{
+	double u = double(rand())/double(RAND_MAX);
+	return a+u*(b-a);
+}
+
+int main(int argc, char** argv)
+{
+	srand(atoi(argv[1]));
+	double w = 120.0*pi;
+	double t = uniform(0.0,1000.0);
+	double theta = uniform(-pi,pi);
+	double V0 = uniform(0.5,1.5);
+	double va = V0*cos(w*t+theta);
+	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;
+	return 0;
+}
+