Add weather sensitivities (#201)

docs/Module/Pypower.md

@@ -57,10 +57,10 @@ GridLAB-D `pypower` module classes.
 ~~~
 class bus {
     int32 bus_i; // bus number (1 to 29997)
-    complex S[MVA]; // base load demand not counting child objects, copied from Pd,Qd by default (MVA)
+    complex S[MVA]; // base load demand not counting child objects, including weather sensitivities, copied to Pd,Qd by default (MVA)
     enumeration {PQREF=1, NONE=4, REF=3, PV=2, PQ=1, UNKNOWN=0} type; // bus type (1 = PQ, 2 = PV, 3 = ref, 4 = isolated)
-    double Pd[MW]; // real power demand (MW)
-    double Qd[MVAr]; // reactive power demand (MVAr)
+    double Pd[MW]; // (OUTPUT) real power demand (MW)
+    double Qd[MVAr]; // (OUTPUT) reactive power demand (MVAr)
     double Gs[MW]; // shunt conductance (MW at V = 1.0 p.u.)
     double Bs[MVAr]; // shunt susceptance (MVAr at V = 1.0 p.u.)
     int32 area; // area number, 1-100
@@ -74,9 +74,34 @@ class bus {
     double lam_Q; // Lagrange multiplier on reactive power mismatch (u/MVAr)
     double mu_Vmax; // Kuhn-Tucker multiplier on upper voltage limit (u/p.u.)
     double mu_Vmin; // Kuhn-Tucker multiplier on lower voltage limit (u/p.u.)
-}
+    char1024 weather_file; // Source object for weather data
+    char1024 weather_variables; // Weather variable column names (col1,col2,...)
+    double weather_resolution[s]; // Weather time downsampling resolution (s)
+    double Sn[W/m^2]; // Solar direct normal irradiance (W/m^2)
+    double Sh[W/m^2]; // Solar horizontal irradiance (W/m^2)
+    double Sg[W/m^2]; // Solar global irradiance (W/m^2)
+    double Wd[deg]; // Wind direction (deg)
+    double Ws[m/2]; // Wind speed (m/2)
+    double Td[degC]; // Dry-bulb air temperature (degC)
+    double Tw[degC]; // Wet-bulb air temperature (degC)
+    double RH[%]; // Relative humidity (%)
+    double PW[in]; // Precipitable_water (in)
+    double HI[degF]; // Heat index (degF)
+    char1024 weather_sensitivity; // Weather sensitivities {PROP: VAR[ REL VAL],SLOPE[; ...]}
+}}
 ~~~
 
+Weather sensitivities are cumulative. The following relations may be used to
+construct the piecewise linear sensitivity curve:
+
+* `>` the cutoff value applies the slope when the weather variable is greater than the cutoff value
+
+* `<` the cutoff value applies the slope when the weather variable is less than the cutoff value
+
+* `@` the cutoff value is not used for comparison.
+
+Note that the cutoff value is always used to identify the intercept point, i.e., $adjust = (source - cutoff) * slope$.
+
 ## Branch Objects
 
 ~~~

module/pypower/autotest/test_case14_sensitivity.glm

@@ -0,0 +1,40 @@
+#define CASE=14
+#ifexists "../case.glm"
+#define DIR=..
+#endif
+#include "${DIR:-.}/case.glm"
+
+module pypower
+{
+	solver_method NR;
+}
+
+modify pp_bus_1.weather_file ${DIR:-.}/test_case14_sensitivity.csv;
+modify pp_bus_1.weather_variables Sh,Sn,Sg,Wd,Ws,Td,Tw,RH,PW;
+modify pp_bus_1.weather_sensitivity "S.real:Td<10,-0.1;S.real:Td<50,-0.1;S.imag:Td>80,+0.1;S.real:Sh,-0.01";
+
+module tape
+{
+	csv_header_type NAME;
+}
+
+object recorder
+{
+	parent pp_bus_1;
+	property S,Sh,Sn,Sg,Wd,Ws,Td,Tw,RH,PW;
+	file "test_case14_sensitivity_pp_bus_1_record.csv";
+	interval 3600;
+}
+
+#set suppress_repeat_messages=FALSE
+
+clock
+{
+	timezone "PST+8PDT";
+	starttime "2020-01-01 00:00:00 PST";
+	stoptime "2021-01-01 00:00:00 PST";
+}
+
+#ifexists ../case.glm
+#on_exit 0 diff ../test_case14_sensitivity_pp_bus_1_record.csv test_case14_sensitivity_pp_bus_1_record.csv > gridlabd.diff
+#endif

module/pypower/autotest/test_case14_weather.glm

@@ -20,7 +20,7 @@ module tape
 object recorder
 {
 	parent pp_bus_1;
-	property Sh,Sn,Sg,Wd,Ws,Td,Tw,RH,PW;
+	property S,Sh,Sn,Sg,Wd,Ws,Td,Tw,RH,PW;
 	file "test_case14_weather_pp_bus_1_record.csv";
 	interval 3600;
 }

module/pypower/bus.cpp

@@ -32,7 +32,7 @@ bus::bus(MODULE *module)
 				PT_DESCRIPTION, "bus number (1 to 29997)",
 
 			PT_complex, "S[MVA]", get_Pd_offset(),
-				PT_DESCRIPTION, "base load demand not counting child objects, copied from Pd,Qd by default (MVA)",
+				PT_DESCRIPTION, "base load demand not counting child objects, including weather sensitivities, copied to Pd,Qd by default (MVA)",
 
 			PT_enumeration, "type", get_type_offset(),
 				PT_DESCRIPTION, "bus type (1 = PQ, 2 = PV, 3 = ref, 4 = isolated)",
@@ -136,6 +136,9 @@ bus::bus(MODULE *module)
 			PT_double, "HI[degF]", get_HI_offset(),
 				PT_DESCRIPTION, "Heat index (degF)",
 
+			PT_char1024, "weather_sensitivity", get_weather_sensitivity_offset(),
+				PT_DESCRIPTION, "Weather sensitivities {PROP: VAR[ REL VAL],SLOPE[; ...]}",
+
 			NULL)<1)
 		{
 				throw "unable to publish bus properties";
@@ -158,6 +161,7 @@ int bus::create(void)
 
 	// initialize weather data
 	current = first = last = NULL;
+	sensitivity_list = NULL;
 
 	return 1; // return 1 on success, 0 on failure
 }
@@ -184,13 +188,85 @@ int bus::init(OBJECT *parent)
 		return 0;
 	}
 
+	char buffer[strlen(weather_sensitivity)+1];
+	strcpy(buffer,weather_sensitivity);
+	char *next=NULL, *last=NULL;
+	// fprintf(stderr,"weather_sensitivity = '%s'\n",(const char*)get_weather_sensitivity());
+	while ( (next=strtok_r(next?NULL:buffer,";",&last)) != NULL )
+	{
+		char propname[65], varname[65], cutoff_test='@';
+		double cutoff_value=0, slope_value;
+		// PROP: VAR[ REL VAL],SLOPE
+		// fprintf(stderr,"next = '%s'\n",next);
+		if ( sscanf(next,"%64[^:]:%64[^<>]%c%lf,%lf",propname,varname,&cutoff_test,&cutoff_value,&slope_value) == 5 
+			|| sscanf(next,"%64[^:]:%64[^,],%lf",propname,varname,&slope_value) == 3 )
+		{
+			// fprintf(stderr,"sensitivity: %s += %s * %lf if %s %c %lf else 0\n",propname,varname,slope_value,varname,cutoff_test,cutoff_value);
+			gld_property source(my(),varname);
+			if ( ! source.is_valid() )
+			{
+				error("weather_sensitivity source '%s' is not valid",varname);
+				return 0;
+			}
+			else if ( source.get_type() != PT_double )
+			{
+				error("weather_sensitivity source '%s' is not a double",varname);
+				return 0;
+			}
+			SENSITIVITY *sensitivity = new SENSITIVITY;
+			if ( strcmp(propname,"Pd") == 0 || strcmp(propname,"S.real") == 0 )
+			{
+				sensitivity->value = &S.Re();
+			}
+			else if ( strcmp(propname,"Qd") == 0 || strcmp(propname,"S.imag") == 0 )
+			{
+				sensitivity->value = &S.Im();
+			}
+			else
+			{
+				error("property '%s' is not valid",propname);
+			}
+			sensitivity->def = strdup(next);
+			sensitivity->source = (double*)source.get_addr();
+			sensitivity->slope = slope_value;
+			sensitivity->cutoff_test = cutoff_test;
+			sensitivity->cutoff_value = cutoff_value;
+			sensitivity->last_adjustment = 0.0;
+			sensitivity->next = sensitivity_list;
+			sensitivity_list = sensitivity;
+		}
+		else
+		{
+			error("weather_sensitivity '%s' in not valid",next);
+			return 0;
+		}
+	}
+
 	return 1; // return 1 on success, 0 on failure, 2 on retry later
 }
 
 TIMESTAMP bus::precommit(TIMESTAMP t0)
 {
 	get_weather(t0);
 
+	// adjust values with sensitivities
+	for ( SENSITIVITY *sensitivity = sensitivity_list ; sensitivity != NULL ; sensitivity = sensitivity->next )
+	{
+		*(sensitivity->value) -= sensitivity->last_adjustment;
+		if ( ( sensitivity->cutoff_test == '<' && *sensitivity->source < sensitivity->cutoff_value )
+			|| (sensitivity->cutoff_test == '>' && *sensitivity->source > sensitivity->cutoff_value )
+			|| sensitivity->cutoff_test == '@'
+			)
+		{
+			sensitivity->last_adjustment = (*sensitivity->source - sensitivity->cutoff_value) * sensitivity->slope;
+			*(sensitivity->value) += sensitivity->last_adjustment;
+		}
+		else
+		{
+			sensitivity->last_adjustment = 0.0;
+		}
+	}
+
 	return current && current->next ? current->next->t : TS_NEVER;	
 }
 

module/pypower/bus.h

@@ -12,6 +12,19 @@ class bus : public gld_object
 
 	static char256 timestamp_format;
 
+public:
+	typedef struct s_sensitivity
+	{
+		char *def;
+		double *value; 
+		double *source;
+		double slope;
+		char cutoff_test;
+		double cutoff_value;
+		double last_adjustment;
+		struct s_sensitivity *next;
+	} SENSITIVITY;
+
 public:
 	// published properties
 	GL_ATOMIC(int32,bus_i);
@@ -48,6 +61,8 @@ class bus : public gld_object
 	GL_ATOMIC(double,HI);
 #define N_WEATHERDATA 10 // adjust if adding more weather data items
 
+	GL_ATOMIC(char1024,weather_sensitivity);
+
 private:
 
 	bool load_weather(void);
@@ -62,6 +77,8 @@ class bus : public gld_object
 	gld_property *weather_mapper[N_WEATHERDATA];
 	WEATHERDATA *first, *last, *current;
 
+	SENSITIVITY *sensitivity_list;
+
 public:
 
 	complex V;