Add weather loader for pypower module (#188)

Signed-off-by: David P. Chassin <dchassin@slac.stanford.edu>
Signed-off-by: Alyona Teyber <Ivanova.alyona5@gmail.com>
Co-authored-by: Alyona Teyber <Ivanova.alyona5@gmail.com>

module/pypower/Makefile.mk

@@ -25,5 +25,6 @@ module_pypower_pypower_la_SOURCES += module/pypower/powerplant.cpp module/pypowe
 module_pypower_pypower_la_SOURCES += module/pypower/relay.cpp module/pypower/relay.h
 module_pypower_pypower_la_SOURCES += module/pypower/scada.cpp module/pypower/scada.h
 module_pypower_pypower_la_SOURCES += module/pypower/transformer.cpp module/pypower/transformer.h
+module_pypower_pypower_la_SOURCES += module/pypower/weather.cpp module/pypower/weather.h
 
 dist_pkgdata_DATA += module/pypower/pypower_solver.py

module/pypower/autotest/case.glm

@@ -10,7 +10,11 @@ clock
     stoptime "2020-02-01 00:00:00 PST";
 }
 
+#ifexists "case${CASE}.glm"
+#include "case${CASE}.glm"
+#else
 #input "case${CASE}.py" -t pypower
+#endif
 
 #gridlabd -C "case${CASE}.glm" -D "starttime=${starttime}" -o "case${CASE}_ref.json" 
 

module/pypower/autotest/test_case14_load.glm

@@ -9,7 +9,6 @@
 module pypower
 {
     maximum_timestep 3600;
-    save_case TRUE;
 }
 
 object pypower.load

module/pypower/autotest/test_case14_powerline.glm

@@ -7,7 +7,6 @@
 module pypower
 {
 	maximum_timestep 3600;
-	save_case TRUE;
 }
 
 object powerline

module/pypower/autotest/test_case14_powerplant.glm

@@ -9,7 +9,11 @@
 module pypower
 {
     maximum_timestep 3600;
-    save_case TRUE;
+}
+
+module tape
+{
+    csv_header_type NAME;
 }
 
 module tape

module/pypower/autotest/test_case14_ts.glm

@@ -6,6 +6,5 @@
 
 module pypower
 {
-	solver_method GS;
 	maximum_timestep 3600;
 }

module/pypower/autotest/test_case14_weather.glm

@@ -0,0 +1,39 @@
+#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_weather.csv;
+modify pp_bus_1.weather_variables Sh,Sn,Sg,Wd,Ws,Td,Tw,RH,PW;
+
+module tape
+{
+	csv_header_type NAME;
+}
+
+object recorder
+{
+	parent pp_bus_1;
+	property Sh,Sn,Sg,Wd,Ws,Td,Tw,RH,PW;
+	file "test_case14_weather_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_weather_pp_bus_1_record.csv test_case14_weather_pp_bus_1_record.csv > gridlabd.diff
+#endif

module/pypower/bus.cpp

@@ -2,15 +2,18 @@
 // Copyright (C) 2024 Regents of the Leland Stanford Junior University
 
 #include "pypower.h"
+#include <time.h>
 
 EXPORT_CREATE(bus);
 EXPORT_INIT(bus);
 EXPORT_SYNC(bus);
+EXPORT_PRECOMMIT(bus);
 
 CLASS *bus::oclass = NULL;
 bus *bus::defaults = NULL;
 
 static int last_i = 0;
+char256 bus::timestamp_format = ""; // "%Y-%m-%d %H:%M:%S%z"; // RFC-822/ISO8601 standard timestamp
 
 bus::bus(MODULE *module)
 {
@@ -93,6 +96,46 @@ bus::bus(MODULE *module)
 				PT_DESCRIPTION, "Kuhn-Tucker multiplier on lower voltage limit (u/p.u.)",
 				PT_ACCESS, PA_REFERENCE,
 
+			PT_char1024, "weather_file", get_weather_file_offset(),
+				PT_DESCRIPTION, "Source object for weather data",
+
+			PT_char1024, "weather_variables", get_weather_variables_offset(),
+				PT_DESCRIPTION, "Weather variable column names (col1,col2,...)",
+
+			PT_double, "weather_resolution[s]", get_weather_resolution_offset(),
+				PT_DEFAULT, "3600 s",
+				PT_DESCRIPTION, "Weather time downsampling resolution (s)",
+
+			PT_double, "Sn[W/m^2]", get_Sn_offset(),
+				PT_DESCRIPTION, "Solar direct normal irradiance (W/m^2)",
+
+			PT_double, "Sh[W/m^2]", get_Sh_offset(),
+				PT_DESCRIPTION, "Solar horizontal irradiance (W/m^2)",
+
+			PT_double, "Sg[W/m^2]", get_Sg_offset(),
+				PT_DESCRIPTION, "Solar global irradiance (W/m^2)",
+
+			PT_double, "Wd[deg]", get_Wd_offset(),
+				PT_DESCRIPTION, "Wind direction (deg)",
+
+			PT_double, "Ws[m/2]", get_Ws_offset(),
+				PT_DESCRIPTION, "Wind speed (m/2)",
+
+			PT_double, "Td[degC]", get_Td_offset(),
+				PT_DESCRIPTION, "Dry-bulb air temperature (degC)",
+
+			PT_double, "Tw[degC]", get_Tw_offset(),
+				PT_DESCRIPTION, "Wet-bulb air temperature (degC)",
+
+			PT_double, "RH[%]", get_RH_offset(),
+				PT_DESCRIPTION, "Relative humidity (%)",
+
+			PT_double, "PW[in]", get_PW_offset(),
+				PT_DESCRIPTION, "Precipitable_water (in)",
+
+			PT_double, "HI[degF]", get_HI_offset(),
+				PT_DESCRIPTION, "Heat index (degF)",
+
 			NULL)<1)
 		{
 				throw "unable to publish bus properties";
@@ -113,6 +156,9 @@ int bus::create(void)
 		throw "maximum bus entities exceeded";
 	}
 
+	// initialize weather data
+	current = first = last = NULL;
+
 	return 1; // return 1 on success, 0 on failure
 }
 
@@ -131,14 +177,29 @@ int bus::init(OBJECT *parent)
 		S.Im() = Qd;
 	}
 
+	// load weather data, if any
+	if ( get_weather_file()[0] != '\0' && ! load_weather() )
+	{
+		error("unable to load weather_file '%s'",(const char*)get_weather_file());
+		return 0;
+	}
+
 	return 1; // return 1 on success, 0 on failure, 2 on retry later
 }
 
+TIMESTAMP bus::precommit(TIMESTAMP t0)
+{
+	get_weather(t0);
+
+	return current && current->next ? current->next->t : TS_NEVER;	
+}
+
 TIMESTAMP bus::presync(TIMESTAMP t0)
 {
 	// reset to base load
 	Pd = S.Re();
 	Qd = S.Im();
+
 	return TS_NEVER;
 }
 
@@ -153,3 +214,130 @@ TIMESTAMP bus::postsync(TIMESTAMP t0)
 	exception("invalid postsync call");
 	return TS_NEVER;
 }
+
+bool bus::load_weather()
+{
+	// setup weather data
+	char *varlist = strdup(get_weather_variables());
+	char *next=NULL, *last=NULL;
+	int n_var = 0;
+	while ( (next=strtok_r(next?NULL:varlist,",",&last)) != NULL )
+	{
+		gld_property *prop = new gld_property(my(),next);
+		if ( prop == NULL || ! prop->is_valid() )
+		{
+			error("weather variable '%s' is not valid",next);
+			return 0;
+		}
+		weather_mapper[n_var++] = prop;
+	}
+	weather_mapper[n_var] = NULL; // mark end of weather variables with a NULL pointer
+	free(varlist);
+
+	// process weather file
+	bool result = false;
+	FILE *fp = fopen(get_weather_file(),"rt");
+	if ( fp == NULL )
+	{
+		error("file '%s' open for read failed",(const char*)get_weather_file());
+		return false;
+	}
+	char buffer[1024];
+
+	for ( int lineno = 0 ; !feof(fp) && fgets(buffer,sizeof(buffer)-1,fp) != NULL ; lineno++ )
+	{
+		char *data = strchr(buffer,',');
+		if ( data == NULL || ! isdigit(buffer[0]) )
+		{
+			// no/bad data -- ignore line
+			continue;
+		}
+		*data++ = '\0';
+		gld_clock t(buffer);
+		if ( ! t.is_valid() )
+		{
+			error("%s@%d: timestamp '%s' is not valid",(const char*)get_weather_file(),lineno,buffer);
+			result = false;
+			break;
+		}
+		time_t tu = (TIMESTAMP)t;
+		if ( get_weather_resolution() > 0 && tu % (long long)get_weather_resolution() == 0 )
+		{
+			if ( ! add_weather(tu,data) )
+			{
+				error("%s@%d: weather data read failed",(const char*)get_weather_file(),lineno);
+				result = false;
+				break;
+			}
+		}
+		result = true;
+	}
+	fclose(fp);
+
+	// start with first record
+	current = first;
+
+	return result;
+}
+
+bool bus::add_weather(TIMESTAMP t, char *buffer)
+{
+	gld_clock ts(t);
+
+	// check time
+	if ( last && t <= last->t )
+	{
+		error("duplicate or out-of-sequence timestamp (t='%s' <= last->t='%s')",
+			gld_clock(t).get_string().get_buffer(),gld_clock(last->t).get_string().get_buffer());
+		return false;
+	}
+
+	// add new weather record
+	WEATHERDATA *data = new WEATHERDATA;
+	data->t = t;
+	data->next = NULL;
+	if ( last != NULL )
+	{
+		last->next = data;
+	}
+	last = data;
+	if ( first == NULL )
+	{
+		first = data;
+	}
+
+	// parse buffer
+	char *next=NULL, *last=NULL;
+	int n_var = 0;
+	while ( (next=strtok_r(next?NULL:buffer,",",&last)) != NULL )
+	{
+		data->value[n_var++] = atof(next);
+	}
+	while ( n_var < N_WEATHERDATA )
+	{
+		data->value[n_var++] = 0.0;
+	}
+	return true;
+}
+
+bool bus::get_weather(TIMESTAMP t)
+{
+	if ( last == NULL )
+	{
+		return true; // no weather
+	}
+	while ( current && current->t < t )
+	{
+		current = current->next;
+	}
+	if ( current == NULL )
+	{
+		warning("weather data ended unexpectedly");
+		return true;
+	}
+	for ( int n_var = 0 ; n_var < N_WEATHERDATA && weather_mapper[n_var] != NULL ; n_var++ )
+	{
+		weather_mapper[n_var]->setp(current->value[n_var]);
+	}
+	return true;
+}

module/pypower/bus.h

@@ -8,6 +8,9 @@
 
 class bus : public gld_object 
 {
+public:
+
+	static char256 timestamp_format;
 
 public:
 	// published properties
@@ -29,6 +32,35 @@ class bus : public gld_object
 	GL_ATOMIC(double,mu_Vmax);
 	GL_ATOMIC(double,mu_Vmin);
 	GL_ATOMIC(complex,S);
+	GL_ATOMIC(char1024,weather_file);
+	GL_ATOMIC(char1024,weather_variables);
+	GL_ATOMIC(double,weather_resolution);
+
+	GL_ATOMIC(double,Sh);
+	GL_ATOMIC(double,Sn);
+	GL_ATOMIC(double,Sg);
+	GL_ATOMIC(double,Wd);
+	GL_ATOMIC(double,Ws);
+	GL_ATOMIC(double,Td);
+	GL_ATOMIC(double,Tw);
+	GL_ATOMIC(double,RH);
+	GL_ATOMIC(double,PW);
+	GL_ATOMIC(double,HI);
+#define N_WEATHERDATA 10 // adjust if adding more weather data items
+
+private:
+
+	bool load_weather(void);
+	bool add_weather(TIMESTAMP t,char *buffer);
+	bool get_weather(TIMESTAMP t);
+
+	typedef struct s_weatherdata {
+		TIMESTAMP t;
+		double value[N_WEATHERDATA];
+		struct s_weatherdata *next;
+	} WEATHERDATA;
+	gld_property *weather_mapper[N_WEATHERDATA];
+	WEATHERDATA *first, *last, *current;
 
 public:
 
@@ -39,6 +71,7 @@ class bus : public gld_object
 	TIMESTAMP presync(TIMESTAMP t0);
 	TIMESTAMP sync(TIMESTAMP t0);
 	TIMESTAMP postsync(TIMESTAMP t0);
+	TIMESTAMP precommit(TIMESTAMP t0);
 
 public:
 	// internal properties

module/pypower/pypower.cpp

@@ -66,6 +66,7 @@ EXPORT CLASS *init(CALLBACKS *fntable, MODULE *module, int argc, char *argv[])
     new relay(module);
     new scada(module);
     new transformer(module);
+    new weather(module);
 
     gl_global_create("pypower::version",
         PT_int32, &pypower_version,