Small models seem to work properly and larger models are stable. PSSE Raw...

Small models seem to work properly and larger models are stable. PSSE Raw files won't necessarily parse correctly. There is a lot of noise, almost certainly numerical, in the solutions but it is worse for larger models. Not sure what to do about that yet.

data/test1.txt

  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       3.45  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       1.0101  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
-   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    
+   1  2   1  1  1 0 0.01     0.01       0.0       0     0     0     0    0  0.0    0.0    0.0    0.0     0.0    0.0    0.0    
 -999
 LOSS ZONES FOLLOWS                     1 ITEMS
   1 ORNL Test Case

data/test2.txt

+ 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.0  0.1101   1.0     1.0   10.0    10.0    0.0     0.0      0   
+   2  Generator   1 1   2  1.011 0.5100 0.0      0.0       3.0  0.1101   1.0     1.0   10.0    10.0    0.0     0.0      0   
+   3  Load        1 1   2  1.000 0.0000 4.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
+   1  3   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    
+   2  3   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
+-999
+INTERCHANGE DATA FOLLOWS               0 ITEMS
+-999
+TIE LINES FOLLOWS                      0 ITEMS
+-999
+END OF DATA

src/AdmittanceNetwork.cpp

+#include "AdmittanceNetwork.h"
+using namespace std;
+
+AdmittanceNetwork::AdmittanceNetwork(ComplexMatrix* impl):
+	impl(impl)
+{
+}
+
+AdmittanceNetwork::~AdmittanceNetwork()
+{
+	delete impl;
+}
+
+void AdmittanceNetwork::add_line(int i, int j, Complex y)
+{
+	impl->incr(i,i,y);
+	impl->incr(i,j,-y);
+}
+
+void AdmittanceNetwork::add_line(int i, int j, Complex y, double turns, double phase_shift, int side)
+{
+	Complex a = polar(turns,phase_shift);
+	if (i == side)
+		impl->incr(i,i,a*a*y);
+	else
+		impl->incr(i,i,y);
+	impl->incr(i,j,-a*y);
+}
+
+void AdmittanceNetwork::add_self(int i, Complex y)
+{
+	impl->incr(i,i,y);
+}
+
+void AdmittanceNetwork::solve(const Complex* Vsrc, Complex* V)
+{
+	impl->solve(Vsrc,V);
+}
+
+

src/AdmittanceNetwork.h

+#ifndef _AdmittanceNetwork_h_
+#define _AdmittanceNetwork_h_
+#include "ElectricalData.h"
+
+class ComplexMatrix
+{
+	public:
+		virtual void set(int row, int col, Complex val) = 0;
+		virtual Complex get(int row, int col) = 0;
+		virtual void incr(int row, int col, Complex val) = 0;
+		virtual void solve(const Complex* Vsrc, Complex* V) = 0;
+		virtual ~ComplexMatrix(){}
+};
+
+class AdmittanceNetwork
+{
+	public:
+		// The impl object is owned by the network
+		AdmittanceNetwork(ComplexMatrix* impl);
+		// Add a line with a transformer from i to j with the transformer on the specified side (i or j,
+		// with the specified phase shift in radians and number of turns. This only adds the equations
+		// for the i side. 
+		void add_line(int i, int j, Complex y, double turns, double phase_shift, int side);
+		// Add a line from i to j 
+		void add_line(int i, int j, Complex y);
+		// 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);
+		// Get matrix entry
+		Complex get(int row, int col) { return impl->get(row,col); }
+		// Destructor
+		~AdmittanceNetwork();
+	private:
+		ComplexMatrix* impl;
+};
+
+#endif
+

src/ComplexMatrix_LAPACK.cpp

+#include "ComplexMatrix_LAPACK.h"
+#include <iostream>
+#include <cassert>
+using namespace std;
+
+// LAPACK function for solving a complex Ax = b
+extern "C"
+{
+	void zgetrs_(char*,int*,int*,double*,int*,int*,double*,int*,int*);
+	void zgetrf_(int*,int*,double*,int*,int*,int*);
+};
+
+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]),
+	piv_lapack(new int[N])
+{
+	for (int i = 0; i < 2*N*N; i++)
+		Y[i] = 0.0;
+}
+
+static int index(int row, int col, int size)
+{
+	return 2*(col*size+row);
+}
+
+void ComplexMatrix_LAPACK::incr(int i, int j, Complex y)
+{
+	int real_part = index(i,j,N);
+	int imag_part = real_part+1;
+	Y[real_part] += real(y);
+	Y[imag_part] += imag(y);
+}
+
+void ComplexMatrix_LAPACK::set(int i, int j, Complex y)
+{
+	int real_part = index(i,j,N);
+	int imag_part = real_part+1;
+	Y[real_part] = real(y);
+	Y[imag_part] = imag(y);
+}
+
+Complex ComplexMatrix_LAPACK::get(int i, int j)
+{
+	int real_part = index(i,j,N);
+	int imag_part = real_part+1;
+	return Complex(Y[real_part],Y[imag_part]);
+}
+
+void ComplexMatrix_LAPACK::solve(const Complex* Vsrc, Complex* V)
+{
+	int i,info,SIZE=N,lda=N,ldb=N,nrhs=1;
+	char type = 'N';
+	// Do the LU decomposition
+	for (int i = 0; i < 2*N*N; i++) LUY[i] = Y[i];
+	zgetrf_(&SIZE,&SIZE,LUY,&lda,piv_lapack,&info);
+	if (info != 0)
+	{
+		cerr << "zgetrf_ returned " << info << endl;
+		cerr.flush();
+	}
+	assert(info == 0);
+	// 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]);
+	}
+	// Solve for the voltage
+	zgetrs_(&type,&SIZE,&nrhs,LUY,&lda,piv_lapack,v_lapack,&ldb,&info);
+	if (info != 0)
+	{
+		cerr << "zgetrs_ returned " << info << endl;
+		cerr.flush();
+	}
+	assert(info == 0);
+	// 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]);
+	}
+}
+
+ComplexMatrix_LAPACK::~ComplexMatrix_LAPACK()
+{
+	delete [] Y;
+	delete [] LUY;
+	delete [] v_lapack;
+	delete [] piv_lapack;
+}
+

src/ComplexMatrix_LAPACK.h

+#ifndef _lapack_h_
+#define _lapack_h_
+#include "AdmittanceNetwork.h"
+
+/**
+ * This is the admittance matrix for the transmission network.
+ * It can be used to solve for current given the voltage or 
+ * vice versa. This implementation uses LAPACK to solve the
+ * Ax=b problem.
+ */
+class ComplexMatrix_LAPACK:
+	public ComplexMatrix
+{
+	public:
+		// Create an N bus network
+		ComplexMatrix_LAPACK(int N);
+		void solve(const Complex* Vsrc, Complex* V);
+		// Get the value at entry i,j
+		Complex get(int row, int col);
+		void set(int row, int col, Complex y);
+		void incr(int row, int col, Complex y);
+		// Destructor
+		~ComplexMatrix_LAPACK();
+	private:
+		// The admittance network stored for LAPACK
+		const int N;
+		double *Y, *LUY;
+		// Scratch space for the solvers
+		double *v_lapack;
+		int* piv_lapack;
+};
+
+#endif

src/ElectricalData.cpp

 #include "ElectricalData.h"
+#include "AdmittanceNetwork.h"
+#include "ComplexMatrix_LAPACK.h"
 #include <iostream>
 #include <cassert>
 using namespace std;
 
-/**
- * This constructor sets the default parameters for the models
- * of the generators. The default data is designed to ensure
- * a stable system. The values are based on typical values
- * described in "Computer-Aided Power System Analysis" by
- * Natarajan.
- */
-ElectricalData::genr_t::genr_t():
-	H(5.0),
-	// This choice leads to a stable generator
-	// with a "good" damping
-	gp(200.0),
-	gd(1.0),
-	gi(0.01),
-	Te(10.0),
-	Vref(1.0),
-	w0(0.0),
-	Ef0(0.0),
-	T0(0.0),
-	Pm0(0.0),
-	genr_mw(0.0),
-	genr_mvar(0.0)
+static const double pi = 3.14159265358979;
+
+static const double phase_angle[3] =
+{
+	0.0,
+	2.0*pi/3.0,
+	4.0*pi/3.0
+};
+
+ElectricalData::ElectricalData(double freq):
+	freq(freq)
 {
 }
 
-void ElectricalData::configureGenerator(genr_t& genr_data, double mvar_base, double kv_base)
+ElectricalData::~ElectricalData()
+{
+}
+
+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.
 	 */
-	if (genr_data.genr_mw * mvar_base > 600.0) // Big nuke plant
+	if (genr_data.genr_mw * SBase > 600.0) // Big nuke plant
 	{
 		genr_data.H = 3.5*genr_data.genr_mw;
 	}
-	else if (genr_data.genr_mw * mvar_base > 300.0) // Big fossil plant
+	else if (genr_data.genr_mw * SBase > 300.0) // Big fossil plant
 	{
 		genr_data.H = 4.0*genr_data.genr_mw;
 	}
 
-	else if (genr_data.genr_mw * mvar_base > 150.0) // Smaller steam plant
+	else if (genr_data.genr_mw * SBase > 150.0) // Smaller steam plant
 	{
 		genr_data.H = 6.0*genr_data.genr_mw;
 	}
-	// If the plant is too small then it does not have automatic controls
 	else 
 	{
 		genr_data.H = 2.0;
 	}
 	// Try to get close to an appropriate set of initial conditions
-	genr_data.Ef0 = 1.0;
-	genr_data.w0 = 1.0;
-	genr_data.Pm0 = genr_data.genr_mw;
+	genr_data.w0 = 0.0;
+	genr_data.Tspd1 = 10.0;
+	genr_data.Tspd2 = 1.0;
+	genr_data.R = 1.0/0.05;
+	genr_data.Agc = genr_data.R/100.0;
+	genr_data.D = 1.0;
+	genr_data.ga = 10000.0;
+	genr_data.gb = 50000.0;
 }
 
-Complex ElectricalData::getAdmittance(int node)
+void ElectricalData::initialize_circuit()
 {
-	// Passive loads. Admittance is calculated from the
-	// assumed power draw.
-	bus_data_t b = nodes[node];
-	Complex S(b.load_mw,b.load_mvar);
-	Complex V = polar(b.v,b.theta);
-	Complex Ii = conj(S)/conj(V);
-	Complex Y = Ii/V + b.B;
-	return Y;
-}
-
-void ElectricalData::segmentLines(int segments)
-{
-	bus_data_t node;
-	node.v = 0.0;
-	node.theta = 0.0;
-	node.load_mw = 0.0;
-	node.load_mvar = 0.0;
-	node.B = Complex(0.0,0.0);
-	vector<line_t> initLines(lines);
-	lines.clear();
-	for (auto line: initLines)
+	map<int,int> node_to_idx;
+	// Move line charging to shunt capacitances at the busses and
+	// calculate the line electrical parameters
+	for (unsigned i = 0; i < lines.size(); i++)
+	{
+		assert(real(lines[i].B) == 0.0);
+		for (unsigned j = 0; j < nodes.size(); j++)
+		{
+			if (lines[i].from == nodes[j].ID)
+				nodes[j].B += lines[i].B/2.0;
+			else if (lines[i].to == nodes[j].ID)
+				nodes[j].B += lines[i].B/2.0;
+		}
+		// Calculate line parameters
+		lines[i].B = Complex(0.0,0.0);
+		lines[i].Z = 1.0/lines[i].y;
+		lines[i].R = real(lines[i].Z);
+		assert(imag(lines[i].Z) > 0.0);
+		lines[i].L = imag(lines[i].Z)/(2.0*pi*freq);
+	}
+	// Calculate load parameters
+	for (unsigned i = 0; i < nodes.size(); i++)
+	{
+		node_to_idx[nodes[i].ID] = i;
+		nodes[i].y = nodes[i].Z = Complex(0.0,0.0);
+		nodes[i].R = nodes[i].L = nodes[i].C = nodes[i].Cshunt = 0.0;
+		if (nodes[i].B != Complex(0.0,0.0))
+		{
+			assert(real(nodes[i].B) == 0.0);
+			if (imag(nodes[i].B) < 0.0)
+				nodes[i].y += nodes[i].B;
+			// Charging as capacitance
+			else
+				nodes[i].Cshunt = imag(nodes[i].B)/(2.0*pi*freq);
+		}
+		Complex V = polar(nodes[i].vamp,nodes[i].theta);
+		Complex S(nodes[i].load_mw,nodes[i].load_mvar);
+		if (S != Complex(0.0,0.0))
+		{
+			Complex I = conj(S/V);
+			nodes[i].Z = V/I;
+			nodes[i].y += 1.0/nodes[i].Z;
+			nodes[i].Z = 1.0/nodes[i].y;
+		}
+		nodes[i].R = real(nodes[i].Z);
+		if (imag(nodes[i].Z) > 0.0)
+			nodes[i].L = imag(nodes[i].Z)/(2.0*pi*freq);
+		else if (imag(nodes[i].Z) < 0.0)
+		{
+			nodes[i].C = -1.0/(imag(nodes[i].Z)*(2.0*pi*freq));
+			assert(nodes[i].C > 0.0);
+		}
+		// Lump a purely capacitive load into the shunt capacitor
+		if (nodes[i].R == 0.0 && nodes[i].L == 0.0)
+		{
+			nodes[i].Cshunt += nodes[i].C;
+			nodes[i].C = 0.0;
+			nodes[i].y = Complex(0.0,0.0);
+			nodes[i].Z = 1.0/nodes[i].y;
+		}
+		// Must have a small shunt capacitance
+		if (nodes[i].Cshunt == 0.0)
+		{
+			nodes[i].Cshunt = 1E-5;
+		}
+		nodes[i].B = Complex(0,2.0*pi*freq*nodes[i].Cshunt);
+	}
+	// Find initial generator currents and voltages
+	AdmittanceNetwork* Y = new AdmittanceNetwork(new ComplexMatrix_LAPACK(nodes.size()));
+	Complex *Vsrc = new Complex[nodes.size()];
+	Complex *Vnode = new Complex[nodes.size()];
+	// Node to ground currents or source voltages for generators
+	for (unsigned i = 0; i < nodes.size(); i++)
 	{
-		line.y = line.y*double(segments);
-		line.B = line.B/double(segments);
-		int end = line.to;
-		lines.push_back(line);
-		for (int j = 1; j < segments; j++)
+		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);
+		}
+		else
+		{
+			Vsrc[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);
+		}
+	}
+	Y->solve(Vsrc,Vnode);
+	// Node initial conditions
+	for (auto entry: node_to_idx)
+	{
+		nodes[entry.second].V = Vnode[entry.second];
+		nodes[entry.second].vamp = abs(Vnode[entry.second]);
+		nodes[entry.second].theta = arg(Vnode[entry.second]);
+		Complex S = Vnode[entry.second]*conj(Vnode[entry.second]*nodes[entry.second].y);
+		nodes[entry.second].load_mw = real(S);
+		nodes[entry.second].load_mvar = imag(S);
+		for (int i = 0; i < 3; i++)
+		{
+			nodes[entry.second].iload[i] = nodes[entry.second].vload[i] = 0.0;
+			nodes[entry.second].v[i] = nodes[entry.second].vamp*cos(phase_angle[i]+nodes[entry.second].theta);
+		}
+		nodes[entry.second].Iload = nodes[entry.second].Vload = Complex(0,0);
+		if (nodes[entry.second].C > 0.0)
+		{
+			assert(nodes[entry.second].L == 0.0);
+			Complex Iload = Vnode[entry.second]*nodes[entry.second].y;
+			nodes[entry.second].Vload = Iload*Complex(0.0,imag(nodes[entry.second].Z));
+			for (int i = 0; i < 3; i++)
+			{
+				nodes[entry.second].vload[i] =
+					abs(nodes[entry.second].Vload)*cos(arg(nodes[entry.second].Vload)+phase_angle[i]);
+			}
+		}
+		else if (nodes[entry.second].L > 0.0)
+		{
+			assert(nodes[entry.second].C == 0.0);
+			nodes[entry.second].Iload = Vnode[entry.second]*nodes[entry.second].y;
+			for (int i = 0; i < 3; i++)
+				nodes[entry.second].iload[i] =
+					abs(nodes[entry.second].Iload)*cos(arg(nodes[entry.second].Iload)+phase_angle[i]);
+		}
+	}
+	for (unsigned i = 0; i < lines.size(); i++)
+	{
+		lines[i].number = i;
+		Complex va = nodes[node_to_idx[lines[i].from]].V;
+		Complex vb = nodes[node_to_idx[lines[i].to]].V;
+		lines[i].I = (va-vb)*lines[i].y;
+		for (int j = 0; j < 3; j++)
+			lines[i].i[j] =
+				abs(lines[i].I)*cos(arg(lines[i].I)+phase_angle[j]);
+	}
+	// Calculate generator initial conditions
+	for (auto entry: node_to_idx)
+	{
+		if (nodes[entry.second].genr.genr_mw == 0.0 && nodes[entry.second].genr.genr_mvar == 0.0)
+			continue;
+		Complex Igenr = nodes[entry.second].V*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 S = Vnode[entry.second]*conj(Igenr);
+		nodes[entry.second].genr.genr_mw = real(S);
+		nodes[entry.second].genr.genr_mvar = imag(S);
+		nodes[entry.second].genr.Pm0 = 0.0;
+		for (int i = 0; i < 3; i++)
 		{
-			lines.back().to = node.ID = nodes.size();
-			ElectricalData::line_t seg;
-			seg.from = node.ID;
-			seg.y = line.y;
-			seg.B = line.B;
-			seg.tside = -1;
-			if (j == segments - 1)
-				seg.to = end;
-			lines.push_back(seg);
-			nodes.push_back(node);
+			double ig = abs(Igenr)*cos(arg(Igenr)+phase_angle[i]);
+			double vg = abs(Vnode[entry.second])*cos(arg(Vnode[entry.second])+phase_angle[i]);
+			nodes[entry.second].genr.Pm0 += ig*vg;
 		}
+		configureGenerator(nodes[entry.second].genr);
 	}
+	delete Y;
+	delete [] Vsrc;
+	delete [] Vnode;
+	
 }

src/ElectricalData.h

@@ -22,14 +22,13 @@ class ElectricalData
 		{
 			/// Default constructor
 			line_t():
-				number(0),from(0),to(0),y(),B(),tside(-1){}
+				from(0),to(0),y(0,0),B(0,0),tside(-1),turns(0),L(0),R(0),Z(1.0/y),I(0,0){}
+			int number;
 			/**
-			 * Create a line from not "from" to "to" with admittance y.
+			 * Create a line from "from" to "to" with admittance y.
 			 */
-			line_t(int number, int from, int to, Complex y):
-				number(number),from(from),to(to),y(y),B(),tside(-1){}
-			/// Unique line number
-			int number;
+			line_t(int from, int to, Complex y):
+				from(from),to(to),y(y),B(),tside(-1){}
 			/// Nodes connected by the line
 			int from, to;
 			/// Admittance of the line
@@ -38,69 +37,93 @@ class ElectricalData
 			Complex B;
 			/// Transformer side. This should be to or from or -1 if there is no transformer
 			int tside;
+			/// Transformer turns ratio and initial condition
+			double turns;
+			/// Inductance
+			double L;
+			/// Resistance
+			double R;
+			/// Impedance
+			Complex Z;
+			/// Initial current
+			Complex I;
+			/// Initial current by phase
+			double i[3];
 		};
 		/**