Updated PSSE raw parser. Picked generator control parameterse that should be...

Updated PSSE raw parser. Picked generator control parameterse that should be stable, and appear to be so in practice.

src/ElectricalData.cpp

@@ -11,8 +11,12 @@ using namespace std;
  * Natarajan.
  */
 ElectricalData::genr_t::genr_t():
-	M(5.0),
-	R(1.0/0.05),
+	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),
@@ -20,8 +24,7 @@ ElectricalData::genr_t::genr_t():
 	T0(0.0),
 	Pm0(0.0),
 	genr_mw(0.0),
-	genr_mvar(0.0),
-	agc(R)
+	genr_mvar(0.0)
 {
 }
 
@@ -34,21 +37,21 @@ 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.H = 3.5*genr_data.genr_mw;
 	}
 	else if (genr_data.genr_mw * mvar_base > 300.0) // Big fossil plant
 	{
-		genr_data.M = 4.0*genr_data.genr_mw;
+		genr_data.H = 4.0*genr_data.genr_mw;
 	}
 
 	else if (genr_data.genr_mw * mvar_base > 150.0) // Smaller steam plant
 	{
-		genr_data.M = 6.0*genr_data.genr_mw;
+		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.M = 2.0;
+		genr_data.H = 2.0;
 	}
 	// Try to get close to an appropriate set of initial conditions
 	genr_data.Ef0 = 1.0;

src/ElectricalData.h

@@ -48,9 +48,11 @@ class ElectricalData
 			/// initial condition
 			genr_t();
 			/// Inertia constant
-			double M;
-			/// inverse of Droop
-			double R;
+			double H;
+			/// These are the speed control parameters
+			double gp; // Proportional gain
+			double gd; // Derivative gain
+			double gi; // Integral gain
 			/// Time constant for the exicitation contoller
 			double Te;
 			/// Voltage magnitude reference
@@ -65,8 +67,6 @@ class ElectricalData
 			double Pm0;
 			/// Power for plant sizing information
 			double genr_mw, genr_mvar;
-			/// AGC control gain
-		 	double agc;
 		};
 		struct bus_data_t
 		{

src/Makefile

@@ -16,15 +16,16 @@ INCLUDE =
 OBJS = \
 	ElectricalData.o \
 	IEEE_CDF_Data.o \
+	PSSE_Raw_Data.o \
 	main.o 
 
 test1: app 
 	flex m2c.lex
-	g++ -Wall -o m2c lex.yy.c
-	./grid2ode ../data/ieee14cdf.txt
+	#g++ -Wall -o m2c lex.yy.c
+	#./grid2ode ../data/ieee14cdf.txt
 	#./grid2ode ../data/test1.txt
-	./m2c circuit.mo
-	${CXX} ${CFLAGS} ${INCLUDE} func.cpp sim.cpp poly.cpp
+	#./m2c circuit.mo
+	#${CXX} ${CFLAGS} ${INCLUDE} func.cpp sim.cpp poly.cpp
 
 app: ${OBJS}
 	${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}

src/PSSE_Raw_Data.cpp

@@ -36,7 +36,7 @@ const int PSSE_Raw_Data::LINE_LEN = 5000;
 PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 	ElectricalData(Freq)
 {
-	map<int,int> busID_to_index;
+	std::map<int,unsigned> bus_id_to_idx;
 	char* line = new char[LINE_LEN];
 	ifstream fin(data_file);
 	if (fin.bad() && !fin)
@@ -47,11 +47,10 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 	fin.getline(line,LINE_LEN,'\r');
 	stringstream sin(line);
 	double Sbase; int IC; sin >> IC >> Sbase;
-	cout << "Sbase = " << Sbase << endl;
 	// Toss the next two lines to get at the bus data
 	fin.getline(line,LINE_LEN,'\r'); fin.getline(line,LINE_LEN,'\r');
 	// Read the bus information
-	int flag, index = 0;
+	int flag;
 	do
 	{
 		fin.getline(line,LINE_LEN,'\r');
@@ -59,7 +58,6 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 		bus_data_t bus_data;
 		if (flag != 0)
 		{
-			// Make sure indices are contiguous
 			bus_data.ID = flag;
 			// Skip the name
 			GETTOK(NULL);
@@ -68,20 +66,12 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 			// Skip the bus type code
 			GETTOK(NULL);
 			// Get the complex shunt admittance. PSSE provides these as power at 1 ang 0 volts
+			// Skip these. Sometimes they are inductive. Not sure what is going on with this.
+			// Maybe add them to the load?
 			double B_real = atof(GETTOK(NULL));
 			double B_imag = atof(GETTOK(NULL));
-			bus_data.B = Complex(B_real,B_imag);
-			#define ABS(_x) (_x < 0 ? -(_x) : _x)
-			#define ISSMALL(_x) (ABS(_x) < 1e-6)
-			if( ISSMALL(Sbase) )
-			{
-			    cerr << "Ignoring small Sbase " << Sbase << endl;
-			    bus_data.B = 0;
-			}
-			else
-			{
-			    bus_data.B /= Sbase;
-			}
+			//bus_data.B = Complex(B_real,B_imag);
+			//bus_data.B /= Sbase;
 			// Skip the area and zone numbers
 			GETTOK(NULL);
 			GETTOK(NULL);
@@ -91,17 +81,14 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 			// Zero the flow of power
 			bus_data.load_mw = 0.0;
 			bus_data.load_mvar = 0.0;
-			bus_data.genr_mw = 0.0;
-			bus_data.genr_mvar = 0.0;
+			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)
 			{
-				busID_to_index[bus_data.ID] = index;
-				bus_data.ID = index++;
+				bus_id_to_idx[bus_data.ID] = nodes.size();
 				nodes.push_back(bus_data);
 			}
-			cout << "BUS " << flag << " V = " << bus_data.v <<
-				" ang. " << bus_data.theta/DEGS_TO_RADS << endl;
 		}
 	}
 	while (flag != 0);
@@ -111,10 +98,10 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 		fin.getline(line,LINE_LEN,'\r');
 		flag = atoi(GETTOK(line)); // First token is the bus ID
 		// Make sure the bus exists; if not, skip it
-		if (flag != 0 && busID_to_index.find(flag) != busID_to_index.end())
+		if (flag != 0 && bus_id_to_idx.find(flag) != bus_id_to_idx.end())
 		{
 			// Get the bus data
-			bus_data_t& bus_data = nodes[busID_to_index[flag]];
+			bus_data_t& bus_data = nodes[bus_id_to_idx[flag]];
 			// Skip the bus name, status, area, and zone
 			GETTOK(NULL);
 			GETTOK(NULL);
@@ -130,20 +117,6 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 			// This is the constant admittance part
 			bus_data.load_mw += atof(GETTOK(NULL))/Sbase;
 			bus_data.load_mvar += atof(GETTOK(NULL))/Sbase;
-			// No negative loads, turn it into a generator
-			if (bus_data.load_mw < 0.0)
-			{
-				bus_data.genr_mw -= bus_data.load_mw;
-				bus_data.genr_mvar -= bus_data.load_mvar;
-				bus_data.load_mw = 0.0;
-				bus_data.load_mvar = 0.0;
-				// Update the data for this generator
-				assert(bus_data.genr_mw > 0.0);
-				addGenr(bus_data,Sbase);
-			}
-			cout << "BUS " << flag << " Load P,Q = " << bus_data.load_mw*Sbase << "," 
-				<< bus_data.load_mvar*Sbase << " , Genr P,Q = " <<
-				bus_data.genr_mw*Sbase << "," << bus_data.genr_mvar*Sbase << endl;
 		}
 	}
 	while (flag != 0);
@@ -152,7 +125,6 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 	bool hasshuntdata=false; if(getenv("PSSEHASSHUNTDATA")) hasshuntdata=true;
 	if(hasshuntdata)
 	{
-		cout << "*** Skipping shunt data ***" << endl;
         	do
 		{
 			fin.getline(line,LINE_LEN,'\r');
@@ -161,7 +133,6 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 	}
 	//------------------------------------------------------------------------
 	// Read generator data
-	index = 0; // For reporting the generator ID number
 	char savedline[10000]; int maxsavedline = 10000; //XXX-Hack
 	do
 	{
@@ -169,43 +140,18 @@ PSSE_Raw_Data::PSSE_Raw_Data(const char* data_file, double Freq):
 		strncpy( savedline, line, maxsavedline-1 ); savedline[maxsavedline-1] = 0; //XXX-Hack
 		flag = atoi(GETTOK(line)); // First token is the bus ID
 		// Make sure the bus exists
-		if (flag != 0 && busID_to_index.find(flag) != busID_to_index.end())
+		if (flag != 0 && bus_id_to_idx.find(flag) != bus_id_to_idx.end())
 		{
 			// Get the bus data
-			bus_data_t& bus_data = nodes[busID_to_index[flag]];
+			bus_data_t& bus_data = nodes[bus_id_to_idx[flag]];
 			// Skip the bus name
 			GETTOK(NULL);
 			// Real and reactive power
 			double mw = atof(GETTOK(NULL))/Sbase;
 			double mvar = atof(GETTOK(NULL))/Sbase;
-			// Generators must produce power
-			if (mw <= 0.0)
-			{
-				bus_data.load_mw -= mw;
-				bus_data.load_mvar -= mvar;
-			}
-			else 
-			{
-				bus_data.genr_mw += mw;
-				bus_data.genr_mvar += mvar;
-				// If this is a new bus, update the index for reporting
-				// if (genrs.find(bus_data.ID) == genrs.end()) index++;
-				// Update the data for this generator
-				if( !(bus_data.genr_mw > 0.0) )
-				{
-				    cerr << "Generator turned off? MW=" << bus_data.genr_mw << endl;
-				    cerr << "    -- Entry: \"" << savedline << "\"" << endl;
-				    if(0)assert(bus_data.genr_mw > 0.0);
-				}
-				else
-				{
-				    addGenr(bus_data,Sbase);
-				}
-			}
-			cout << "BUS " << flag << "(" << index << ") Genr P,Q = " << bus_data.genr_mw*Sbase << "," 
-				<< bus_data.genr_mvar*Sbase << endl;
-			cout << "BUS " << flag << " Load P,Q = " << bus_data.load_mw*Sbase << "," 
-				<< bus_data.load_mvar*Sbase << endl;
+			bus_data.genr.genr_mw += mw;
+			bus_data.genr.genr_mvar += mvar;
+			configureGenerator(bus_data.genr,Sbase);
 		}
 	}while (flag != 0);
 	int line_count = 0;
@@ -223,14 +169,14 @@ skip_line:
 			// end of the bus (see the RAW file format documentation)
 			int to_bus = abs(atoi(GETTOK(NULL)));
 			// Remove lines attending busses that are out of service
-			if (busID_to_index.find(flag) == busID_to_index.end() ||
-					busID_to_index.find(to_bus) == busID_to_index.end())
+			if (bus_id_to_idx.find(flag) == bus_id_to_idx.end() ||
+					bus_id_to_idx.find(to_bus) == bus_id_to_idx.end())
 				goto skip_line;
 			// Assign a number to this bus
 			line_data.number = line_count;
 			// Get the starting and ending bus IDs
-			line_data.from = busID_to_index[flag];
-			line_data.to = busID_to_index[to_bus];
+			line_data.from = flag;
+			line_data.to = to_bus;
 			// Skip the circuit identifier
 			GETTOK(NULL);
 			// Read the line impedence
@@ -262,13 +208,34 @@ skip_line:
 				// Add the line
 				lines.push_back(line_data);
 			}
-			cout << "LINE " << line_data.from  << "->" << line_data.to << " y = " << line_data.y << endl;
 			line_count++;
 		}
 	}while (flag != 0);
-	cout << "Got " << line_count << " transmission lines" << endl;
 	// Done
 	delete [] line;
+	// Delete nodes that don't have any lines attached to them
+	unsigned idx = 0;
+	while (idx < nodes.size())
+	{
+		bool found_it = false;
+		for (auto line: lines)
+		{
+			if (line.to == nodes[idx].ID || line.from == nodes[idx].ID)
+			{
+				found_it = true;
+				break;
+			}
+		}
+		if (found_it)
+			idx++;
+		else
+		{
+			nodes[idx] = nodes.back();
+			nodes.pop_back();
+		}
+	}
+	cout << "Got " << line_count << " transmission lines" << endl;
+	cout << "Got " << nodes.size() << " busses" << endl;
 }
 
 PSSE_Raw_Data::~PSSE_Raw_Data()

src/m2c.lex

@@ -184,6 +184,7 @@ 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> volt_amp_est;
 	std::list<plot_var> angle;
 
 	for (whichPass = 0; whichPass < 2; whichPass++)
@@ -223,6 +224,8 @@ int main(int argc, char* argv[])
 			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));
+		else if (alg_var.find("Amp") != alg_var.npos)
+			volt_amp_est.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;
@@ -243,6 +246,7 @@ int main(int argc, char* argv[])
 	make_gnuplot("power.gnu",power);
 	make_gnuplot("speed.gnu",speed);
 	make_gnuplot("angle.gnu",angle);
+	make_gnuplot("volt_amplitude.gnu",volt_amp_est);
 
 	return EXIT_SUCCESS;
 }

src/main.cpp

-//#include "PSSE_Raw_Data.h"
+#include "PSSE_Raw_Data.h"
 #include "IEEE_CDF_Data.h"
 #include <iostream>
 #include <fstream>
@@ -173,8 +173,7 @@ 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 Rstr = number(data.R);
+	string Hstr = number(data.H*pi*Freq);
 	string w = omega(busName);
 	string Pm = mech_power(busName);
 	string Pe = "pe_"+busName;
@@ -205,15 +204,17 @@ string createGenerator(string busName, ElectricalData::genr_t data, double cap)
 		result += "der("+i+")=("+generatorVoltage(busName,phase)+"-"+voltage(busName,phase)+"-"+number(resistance)+"*"+i+")/"+number(inductance)+";\n";
 	}
 	result += Pe+"="+powerDraw+";\n";
-	stateVars.insert(state_var_init_t(generatorAngle(busName),data.T0*pi/180.0));
+	stateVars.insert(state_var_init_t(generatorAngle(busName),data.T0));
 	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)+")="+number(Freq)+"*("+w+"-1.0);\n";
 	result += string("der(")+w+string(")=(")+Pm+
 		string("-"+Pe+")/")+Hstr+string(";\n");
-	result += string("der(")+Pm+string(")="+
-		Rstr+"*(1.0-"+w+"+"+Pe+"-"+Pm+");\n"); 
+	result += string("der(")+Pm+")="+
+		number(data.gp)+"*(1.0-"+w+")+"+
+		number(data.gd)+"*("+Pe+"-"+Pm+")-"+
+		number(data.gi)+"*"+generatorAngle(busName)+";\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");
@@ -445,9 +446,9 @@ int main(int argc, char** argv)
 	ElectricalData* data;
 	string fn = argv[argc-1];
 	string ext = fn.substr(fn.find_last_of(".")+1);
-/*	if (ext == "raw")
-		data = new PSSE_Raw_Data(fn.c_str(),Freq); */
-	if (ext == "txt")
+	if (ext == "raw")
+		data = new PSSE_Raw_Data(fn.c_str(),Freq);
+	else if (ext == "txt")
 		data = new IEEE_CDF_Data(fn.c_str(),Freq);
 	else
 	{

src/sim.cpp

@@ -118,14 +118,12 @@ 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++)
 		{
 			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] += p[i].extrapolate(h);