Stable model with proper voltages. Solver is very very slow. Much simplified...

Stable model with proper voltages. Solver is very very slow. Much simplified generator speed control.

src/ElectricalData.cpp

@@ -15,14 +15,12 @@ ElectricalData::genr_t::genr_t():
 	Tspd2(1.0),
 	M(5.0),
 	R(1.0/0.05),
-	Te(100.0),
+	Te(10.0),
 	Vref(1.0),
-	D(1.0),
 	w0(0.0),
 	Ef0(0.0),
 	T0(0.0),
 	Pm0(0.0),
-	C0(0.0),
 	genr_mw(0.0),
 	genr_mvar(0.0),
 	agc(R)
@@ -39,29 +37,25 @@ void ElectricalData::configureGenerator(genr_t& genr_data, double mvar_base, dou
 	if (genr_data.genr_mw * mvar_base > 600.0) // Big nuke plant
 	{
 		genr_data.M = 3.5*genr_data.genr_mw;
-		genr_data.D = genr_data.M/3.5;
 	}
 	else if (genr_data.genr_mw * mvar_base > 300.0) // Big fossil plant
 	{
 		genr_data.M = 4.0*genr_data.genr_mw;
-		genr_data.D = genr_data.M/4.0;
 	}
 
 	else if (genr_data.genr_mw * mvar_base > 150.0) // Smaller steam plant
 	{
 		genr_data.M = 6.0*genr_data.genr_mw;
-		genr_data.D = genr_data.M/6.0;
 	}
 	// If the plant is too small then it does not have automatic controls
 	else 
 	{
 		genr_data.M = 2.0;
-		genr_data.D = 1.0;
 	}
 	// Try to get close to an appropriate set of initial conditions
 	genr_data.Ef0 = 1.0;
 	genr_data.w0 = 1.0;
-	genr_data.C0 = genr_data.Pm0 = genr_data.genr_mw;
+	genr_data.Pm0 = sqrt(genr_data.genr_mw*genr_data.genr_mw+genr_data.genr_mvar*genr_data.genr_mvar);
 }
 
 Complex ElectricalData::getAdmittance(int node)

src/ElectricalData.h

@@ -57,8 +57,6 @@ class ElectricalData
 			double Te;
 			/// Voltage magnitude reference
 			double Vref;
-			/// Resistance to off-nominal operating speeds
-			double D;
 			/// Initial angular velocity (p.u.)
 			double w0;
 			/// Initial exiciter field voltage (p.u.)
@@ -67,8 +65,6 @@ class ElectricalData
 			double T0;
 			/// Initial mechanical power output (p.u.)
 			double Pm0;
-			/// Initial governor signal
-			double C0;
 			/// Power for plant sizing information
 			double genr_mw, genr_mvar;
 			/// AGC control gain

src/Makefile

@@ -4,7 +4,7 @@
 CXX = g++
 CFLAGS = -Wall -Ofast -std=c++11
 LIBS = 
-INCLUDE = 
+INCLUDE = -I/usr/include/eigen3
 
 .SUFFIXES: .cpp
 .cpp.o:
@@ -21,10 +21,10 @@ OBJS = \
 test1: app 
 	flex m2c.lex
 	g++ -Wall -o m2c lex.yy.c
-	#./grid2ode ../data/ieee14cdf.txt
-	./grid2ode ../data/test1.txt
+	./grid2ode ../data/ieee14cdf.txt
+	#./grid2ode ../data/test1.txt
 	./m2c circuit.mo
-	${CXX} ${CFLAGS} func.cpp sim.cpp 
+	${CXX} ${CFLAGS} ${INCLUDE} func.cpp sim.cpp 
 
 app: ${OBJS}
 	${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}

src/m2c.lex

@@ -12,7 +12,6 @@ enum
 	DELAY,
 	DER,
 	INIT,
-	FIX,
 	NONE
 };
 
@@ -20,21 +19,18 @@ struct init_val_t
 {
 	std::string name;
 	double value;
-	bool fixed;
 };
 
 struct init_val_t to_init;
 std::set<std::string> alg_vars;
 std::set<std::string> state_vars;
 std::list<init_val_t> state_vars_init_vals;
-std::ostringstream func, opt;
+std::ostringstream func;
 bool lhs = true;
 bool newLine = true;
 bool skip = false;
-bool residualExpr = false;
 int whichFunc = NONE;
 int whichPass = 0;
-int numFixed = 0;
 
 void end_eqn()
 {
@@ -44,7 +40,6 @@ void end_eqn()
 	if (!skip && whichPass == 1)
 	{
 		func << ";\n";
-		if (residualExpr) opt << ";\n";
 	}
 	skip = false;
 }
@@ -56,7 +51,6 @@ void lhs_to_rhs()
 	if (whichPass == 1)
 	{
 		func << "=";
-		if (residualExpr) opt << "=";
 	}
 }
 
@@ -71,14 +65,9 @@ void identifier(const char* ID)
 				whichFunc = DER;
 			else if (strcmp(ID,"initialize") == 0)
 				whichFunc = INIT;
-			else if (strcmp(ID,"fix") == 0)
-			{
-				whichFunc = FIX;
-				numFixed++;
-			}
 			else if (whichFunc == DER)
 				state_vars.insert(ID);
-			else if (whichFunc == INIT || whichFunc == FIX)
+			else if (whichFunc == INIT)
 				to_init.name = ID;
 			else alg_vars.insert(ID);
 		}
@@ -92,9 +81,7 @@ void identifier(const char* ID)
 		if (newLine)
 		{
 			func << "    ";
-			if (residualExpr) opt << "        ";
 			newLine = false;
-			residualExpr = false;
 		}
 		// First argument of delay uses variable identifier
 		if (whichFunc == DELAY)
@@ -107,29 +94,10 @@ void identifier(const char* ID)
 		{
 			if (lhs) func << "d";
 			func << "q[" << ID << "]";
-			if (lhs)
-			{
-				for (auto var: state_vars_init_vals)
-				{
-					if (var.name == ID && var.fixed)
-					{
-						opt << "residue[j++]";
-						residualExpr = true;
-						break;
-					}
-				}
-			}
-			else if (residualExpr) opt << ID ;
 		}
 		else if (alg_vars.find(ID) != alg_vars.end())
 		{
 			func << "alg[" << ID << "]";
-			if (lhs)
-			{
-				residualExpr = true;
-				opt << "T " << ID;
-			}
-			else if (residualExpr) opt << ID;
 		}	
 		// Otherwise its a function call. der only exists on lhs
 		// to identify state variables
@@ -138,7 +106,6 @@ void identifier(const char* ID)
 			func << ID;
 			if (strcmp(ID,"delay") == 0)
 				whichFunc = DELAY;
-			if (residualExpr) opt << ID;
 		}
 	}
 }
@@ -147,10 +114,7 @@ void mathop(const char* op)
 {
 	if (skip) return;
 	if (whichPass == 1 && !lhs)
-	{
 		func << op;
-		if (residualExpr) opt << op;
-	}
 }
 
 void number(const char* num)
@@ -158,17 +122,15 @@ void number(const char* num)
 	if (skip) return;
 	if (!lhs)
 	{
-		if (whichPass == 0 && (whichFunc == INIT || whichFunc == FIX))
+		if (whichPass == 0 && whichFunc == INIT)
 		{
 			to_init.value = atof(num);
-			to_init.fixed = whichFunc == FIX;
 			state_vars_init_vals.push_back(to_init);
 			whichFunc = NONE;
 		}
 		else if (whichPass == 1)
 		{
 			func << num;
-			if (residualExpr) opt << "T(" << num << ")";
 		}
 	}
 }
@@ -189,12 +151,12 @@ OP [-+*/(),]
 {OP} {mathop(yytext);}
 {ID} {identifier(yytext);}
 {NUMBER} {number(yytext);}
+[.|\n] {}
 
 %%
 
 int main(int argc, char* argv[])
 {
-	opt << "        ";
 	for (whichPass = 0; whichPass < 2; whichPass++)
 	{
 		FILE* fp = fopen(argv[1],"r");
@@ -202,6 +164,7 @@ int main(int argc, char* argv[])
 		yylex();
 		fclose(fp);
 	}
+	int var_idx = 2;
 	int k = 0;
 	std::ofstream func_out("func.cpp");
 	func_out << "#include \"sim.h\"" << std::endl;
@@ -210,26 +173,13 @@ 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)
-		func_out << "static const int " << state_var << " = " << (k++) << ";" << std::endl;
+		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)
-		func_out << "static const int " << alg_var << " = " << (k++) << ";" << std::endl;
+		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 << "const int numFixed = " << numFixed << ";" << std::endl << std::endl;
-
-	func_out << "void init_state(double* const q, const double* const x)";
-	func_out << std::endl << "{" << std::endl << "    int k = 0;" << std::endl;
-	for (auto val: state_vars_init_vals)
-	{
-		if (!val.fixed)
-			func_out << "    q[" << val.name << "] = x[k++];" << std::endl; 
-		else
-			func_out << "    q[" << val.name << "] = " << val.value << ";" << std::endl;
-	}
-	func_out << "}" << std::endl << std::endl;
-
 	func_out << "void func(double* const dq, const double* const q, double* const alg, const double time)";
 	func_out << std::endl << "{" << std::endl;
 	func_out << func.str() << std::endl;
@@ -242,25 +192,5 @@ int main(int argc, char* argv[])
 	func_out << "}" << std::endl << std::endl;
 	func_out.close();
 
-	std::ofstream opt_out("opt.h");
-	opt_out << "struct CostFunctor" << std::endl << "{" << std::endl;
-	opt_out << "    template<typename T>" << std::endl;
-	opt_out << "    bool operator()(const T* const x, T* residue) const" << std::endl;
-	opt_out << "    {" << std::endl;
-	opt_out << "        int k = 0, j = 0;" << std::endl;
-	opt_out << "        static const T time = T(0.0);" << std::endl;
-	opt_out << "        static const T pi =  T(" << 3.141592653589793 << ");" << std::endl;
-	for (auto val: state_vars_init_vals)
-	{
-		if (!val.fixed)
-			opt_out << "        T "  << val.name << " = T(x[k++]);" << std::endl;
-		else
-			opt_out << "        T "  << val.name << " = T(" << val.value << ");" << std::endl;
-	}
-	opt_out << opt.str() << std::endl;
-	opt_out << "        return true;" << std::endl;
-	opt_out << "    }" << std::endl;
-	opt_out << "};" << std::endl;
-	opt_out.close();
 	return EXIT_SUCCESS;
 }

src/main.cpp

@@ -8,33 +8,44 @@
 #include <set>
 using namespace std;
 
-static const double pi = 3.14159265358979;
-int lineSegments = 1;
+const double pi = 3.14159265358979;
+// 60 Hertz systems default. Set at command line
 double Freq = 2.0*pi*60.0; 
-vector<string> phases;
+// Number of segments into which a transmission line will be sliced
+// This is set by a command line argument
+int lineSegments = 1;
+// Enumeration of system phases
+const char* const phases[3] = { "a", "b", "c"};
+// Set of voltage amplitude measurement variables
 set<string> voltageAmplitudeEstimators;
 
-static const double TfRate = 1.0;
-static const double v0 = 0.0;
-static const double i0 = 0.0;
+// Default rate of change for a transformer tap
+const double TfRate = 1E-3;
+// Default initial values for voltages and currents
+const double v0 = 0.0;
+const double i0 = 0.0;
 
+// Initial value data for a state variable
 struct state_var_init_t
 {
+	// Variable name
 	string name;
+	// Initial value
 	double init_value;
-	bool fixed;
 	state_var_init_t(string name, double init_value):
-		name(name),init_value(init_value),fixed(false){}
-	state_var_init_t(string name, double init_value, bool fixed):
-		name(name),init_value(init_value),fixed(fixed){}
+		name(name),init_value(init_value){}
 	state_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
 	{
 		return name < other.name;
 	}
 };
+// Initial variable values for all state variables
 set<state_var_init_t> stateVars;
 
+// Number to string
 string number(double val)
 {
 	ostringstream sin;
@@ -42,6 +53,7 @@ string number(double val)
 	return sin.str();
 }
 
+// Bus ID number to string
 string bus_name(int i)
 {
 	ostringstream sin;
@@ -49,6 +61,7 @@ string bus_name(int i)
 	return sin.str();
 }
 
+// Line ID to string
 string line_name(const ElectricalData::line_t& data)
 {
 	ostringstream sin;
@@ -56,59 +69,61 @@ string line_name(const ElectricalData::line_t& data)
 	return sin.str();
 }
 
+// Bus voltage to string
 string voltage(string busName, string phase)
 {
 	return string("v")+phase+string("_")+busName;
 }
 
+// Generator rotational speed to a string
 string omega(string genName)
 {
 	return string("w_")+genName;
 }
 
-string control_var(string genName)
-{
-	return string("c_")+genName;
-}
-
+// Generator mechanical power to a string
 string mech_power(string genName)
 {
 	return string("pm_")+genName;
 }
 
+// Current to a string
 string current(string flowName, string phase)
 {
 	return string("i")+phase+string("_")+flowName;
 }
 
+// Transformer modulus at a bus and phase to a string
 string transformerModulus(string bus_name, string phase, ElectricalData::bus_data_t& bus)
 {
 	return string("tf_")+bus_name+"_"+phase+"_"+number(bus.ID);
 }
 
+// Line transformer modulus to a string
 string transformerModulus(string bus_name, string phase, ElectricalData::line_t& line)
 {
 	return string("tf_")+bus_name+"_"+phase+"_"+number(line.to)+"_"+number(line.from);
 }
 
+// Generator phase angle to a string
 string generatorAngle(string busName)
 {
 	return string("theta_")+busName;
 }
 
-string generatorSourceVoltage(string busName, string phase)
+// Internal generator excitation voltage
+string generatorVoltage(string busName, string phase)
 {
-	string freq = number(Freq);
-	string w = omega(busName);
-	string theta = generatorAngle(busName);
-	if (phase == phases[0])
-		return voltage(busName,phase)+"=cos("+freq+"*"+w+"*time+"+theta+")";
-	else if (phase == phases[1])
-		return voltage(busName,phase)+"=cos("+freq+"*"+w+"*time-2.0*pi/3.0+"+theta+")";
-	else
-		return voltage(busName,phase)+"=cos("+freq+"*"+w+"*time-4.0*pi/3.0+"+theta+")";
+	return string("vg_")+busName+"_"+phase;
+}
+
+// Generator output current
+string generatorCurrent(string busName, string phase)
+{
+	return string("ig_")+busName+"_"+phase;
 }
 
+// Generate functions to estimate the amplitude of a sinuisoid
 string sineAmplitudeMeasurement(string busName, string phase, string genName, string& estimate)
 {
 	/**
@@ -135,54 +150,101 @@ string sineAmplitudeMeasurement(string busName, string phase, string genName, st
 	return y1_gen+"\n"+y2_gen+"\n"+y3_gen+"\n"+Aest+"\n";
 }
 
-string createGenerator(string busName, string powerDraw, ElectricalData::genr_t data)
+// Generator excitation voltage
+string generatorSourceVoltage(string busName, string phase)
 {
-	string result("");
+	string E = "ef_"+busName+"_"+phase;
+	string freq = number(Freq);
+	string w = omega(busName);
+	string theta = generatorAngle(busName);
+	if (phase == phases[0])
+		return generatorVoltage(busName,phase)+"="+E+"*cos("+freq+"*"+w+"*time+"+theta+")";
+	else if (phase == phases[1])
+		return generatorVoltage(busName,phase)+"="+E+"*cos("+freq+"*"+w+"*time+2.0*pi/3.0+"+theta+")";
+	else
+		return generatorVoltage(busName,phase)+"="+E+"*cos("+freq+"*"+w+"*time+4.0*pi/3.0+"+theta+")";
+}
+
+// Create dynamic model of a generator. The cap argument is the capacitance of the
+// node to which the generator is attached. The internal inductance of the generator
+// is set to make sure the resulting LRC circuit is critically damped.
+string createGenerator(string busName, ElectricalData::genr_t data, double cap)
+{
+	string result;
 	string varName;
 	string Vref = number(data.Vref);
 	string Hstr = number(data.M*Freq);
-	string Tp1str = number(data.Tspd1*Freq);
-	string Tp2str = number(data.Tspd2*Freq);
-	string Rstr = number(data.R*Freq);
+	string Tp1str = number(data.Tspd1);
+	string Tp2str = number(data.Tspd2);
+	string Rstr = number(data.R);
 	string w = omega(busName);
-	string C = control_var(busName);
 	string Pm = mech_power(busName);
+	string Pe = "pe_"+busName;
 	string Pref = "pref_"+busName;
+	string powerDraw;
+	// Small internal impedance to damp oscillations
+	static const double resistance = 1E-3;
+	// Set the inductance to its critically damped value
+	const double inductance = (resistance*resistance*cap/4.0);
 	for (auto phase: phases)
 	{
+		string E = "ef_"+busName+"_"+phase;
+		string Amp;
+		stateVars.insert(state_var_init_t(E,data.Ef0));
+		result += sineAmplitudeMeasurement(busName,phase,busName,Amp);	
+		result += "der("+E+")=("+number(data.Vref)+"-"+Amp+")*"+number(data.Te)+";\n";
+	}
+	// Power demand at the terminal is sum of demand for each phase
+	for (auto phase: phases)
+	{	
+		if (powerDraw.empty())
+			powerDraw = generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
+		else 
+			powerDraw += "+"+generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
 		result += generatorSourceVoltage(busName,phase)+";\n";
+		string i = generatorCurrent(busName,phase);
+		stateVars.insert(state_var_init_t(i,0.0));
+		result += "der("+i+")=("+generatorVoltage(busName,phase)+"-"+voltage(busName,phase)+"-"+number(resistance)+"*"+i+")/"+number(inductance)+"\n;";
 	}
-	stateVars.insert(state_var_init_t(Pref,data.Pm0,true));
-	stateVars.insert(state_var_init_t(generatorAngle(busName),data.T0*pi/180.0,true));
-	stateVars.insert(state_var_init_t(w,data.w0,true));
+	result += Pe+"="+powerDraw+";\n";
+	stateVars.insert(state_var_init_t(generatorAngle(busName),data.T0*pi/180.0));
+	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 += string("der(")+w+string(")=(")+Pm+
-		string("-("+powerDraw+"))/")+Hstr+string(";\n");
-	stateVars.insert(state_var_init_t(Pm,data.Pm0,true));
-	result += string("der(")+Pm+string(")=")+Tp2str+
-		string("*(")+C+string("-")+Pm+string(");\n"); 
-	stateVars.insert(state_var_init_t(C,data.C0,true));
-	result += string("der(")+C+string(")=")+Tp1str+
-		string("*(")+Pref+"-"+Rstr+string("*")+
-		string("(")+w+string("-1.0)-")+C+string(");\n");
-	result += string("der(")+Pref+")="+number(data.agc)+"*(1.0-"+w+");\n";
+		string("-"+Pe+")/")+Hstr+string(";\n");
+	result += string("der(")+Pm+string(")="+Tp2str+"*("+
+		Rstr+"*(1.0-"+w+")+"+Pe+"-"+Pm+");\n"); 
+	// Useful to make the generators sit still when we are testing the network model
+/*	result += "der("+generatorAngle(busName)+")=0.0;\n";
+	result += string("der(")+w+string(")=0.0;\n");
+	result += string("der(")+Pm+string(")=0.0;\n"); 
+	result += string("der(")+C+string(")=0.0;\n");
+	result += string("der(")+Pref+")=0.0;\n"; */ 
 	return result;
 }
 
+// Create an electrical bus, possible with a load
 string createBus(
 		int nodeID,
 		ElectricalData& data)
 {
 	ElectricalData::bus_data_t bus = data.getBusInfo(nodeID);
 	string busName = bus_name(bus.ID);
-	static const double Cmin = 1E-2;
-	double C = 0.0, R = 0.0;//real(data.getAdmittance(nodeID));
+	// Small, minimum value for the capacitance relative to ground
+	static const double Cmin = 1E-5;
+	Complex Z = 1.0/data.getAdmittance(nodeID);
+	double C = 0.0, R = real(Z), Rimg = imag(Z);
 	string result = "";
+	bool hasGenr = (bus.genr.genr_mw != 0.0 || bus.genr.genr_mvar != 0.0);
+	// Half of the line charging goes to each bus attached to the line
 	for (auto line: data.getLines())
 	{
 		if (line.to == bus.ID || line.from == bus.ID)
 			C += fabs(imag(line.B))/2.0;
 	}
+	// Check for valid resistance
 	if (isnan(R) || !isfinite(R) || R < 0.0)
 	{
 		cout << "Bad load resistance " << R 
@@ -190,6 +252,15 @@ string createBus(
 		cout << "Replacing with zero" << endl;
 		R = 0.0;
 	}
+	// Check for valid reactance
+	if (isnan(Rimg) || !isfinite(Rimg))
+	{
+		Rimg = 0.0;
+	}
+	// Load capacitor is in parallel with the shunt/line charging capacitor
+	if (R == 0.0 && Rimg < 0.0)
+		C -= Rimg;
+	// Check for valid capacitance
 	if (isnan(C) || !isfinite(C) || C < Cmin)
 	{
 		cout << "Bad bus capacitance " << C
@@ -199,61 +270,20 @@ string createBus(
 	}
 	string Cstr = number(C/Freq);
 	string Rstr = number(R);
-	// If the voltage at this bus is regulated by a generator
-	if (bus.genr.genr_mw != 0.0 || bus.genr.genr_mvar != 0.0)
-	{
-		// Power is sum of power on phases
-		string NetPower;
-		for (auto phase: phases)
-		{
-			string PowerDraw;
-			for (auto line: data.getLines())
-			{
-				if (line.to == bus.ID)
-				{
-					if (PowerDraw.size() == 0)
-						PowerDraw = current(line_name(line),phase);
-					else
-						PowerDraw = "+"+current(line_name(line),phase);
-				}
-			}
-			for (auto line: data.getLines())
-			{
-				if (line.from == bus.ID)
-				{
-					PowerDraw = "-"+current(line_name(line),phase);
-				}
-			}
-			if (R > 0.0)
-			{
-				PowerDraw += "-"+voltage(busName,phase)+"/"+number(R);
-			}
-			PowerDraw = "("+voltage(busName,phase)+"*("+PowerDraw+"))";
-			if (NetPower.empty())
-				NetPower = PowerDraw;
-			else
-				NetPower += "+"+PowerDraw;
-		}
-		result += createGenerator(busName,NetPower,bus.genr);
-	}
-	// Otherwise this is just a regular load node
-	else for (auto phase: phases)
+	if (hasGenr)
+		result += createGenerator(busName,bus.genr,C/Freq);
+	// Generator the node equations
+	for (auto phase: phases)
 	{
-		// If the voltage can sag, include a transformer to regulate it
-/*		if (R > 0.0)
-		{