From f8f56c886f8321432a8f83ff8811296eab83323b Mon Sep 17 00:00:00 2001 From: "David P. Chassin" Date: Wed, 21 Feb 2024 17:25:37 -0800 Subject: [PATCH] Enabled OPF solver --- converters/py2glm.py | 21 ++- module/powerflow/substation.cpp | 13 +- module/pypower/autotest/case14.glm | 129 ++++++++++++------ .../autotest/test_case14_substation.glm | 3 +- module/pypower/gencost.cpp | 92 +++++-------- module/pypower/gencost.h | 23 +--- module/pypower/pypower.cpp | 19 ++- module/pypower/pypower_solver.py | 38 +++--- 8 files changed, 196 insertions(+), 142 deletions(-) diff --git a/converters/py2glm.py b/converters/py2glm.py index 5db48a084..bdc92fe9c 100644 --- a/converters/py2glm.py +++ b/converters/py2glm.py @@ -71,13 +71,30 @@ def convert(ifile,ofile,py_type): bus = "bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin", gen = "bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf", branch = "fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax", - # gencost = "TODO" ).items(): + glm.write(f"{NL}//{NL}// {name}{NL}//{NL}") for line in data[name]: - glm.write(f"""object {name} + glm.write(f"""object pypower.{name} {{ {NL.join([f" {x} {line[n]};" for n,x in enumerate(spec.split())])} }} +""") + if 'gencost' in data: + glm.write("\n//\n// gencost\n//\n") + for line in data['gencost']: + model = line[0] + startup = line[1] + shutdown = line[2] + count = line[3] + costs = line[4:] + assert(len(costs)==count) + glm.write(f"""object pypower.gencost +{{ + model {int(model)}; + startup {startup}; + shutdown {shutdown}; + costs "{','.join([str(x) for x in costs])}"; +}} """) if __name__ == '__main__': diff --git a/module/powerflow/substation.cpp b/module/powerflow/substation.cpp index ca2b93761..8e44c3709 100644 --- a/module/powerflow/substation.cpp +++ b/module/powerflow/substation.cpp @@ -245,13 +245,22 @@ int substation::init(OBJECT *parent) } else if ( gl_object_isa(parent,"load","pypower") ) { - // TODO: link to pypower bus object fetch_complex(&pPositiveSequenceVoltage,"V",parent); fetch_complex(&pConstantPowerLoad,"P",parent); fetch_complex(&pConstantCurrentLoad,"I",parent); fetch_complex(&pConstantImpedanceLoad,"Z",parent); fetch_double(&pTransNominalVoltage,"Vn",parent); - // throw "substation does not support pypower bus linkage yet"; + if (fabs(*pTransNominalVoltage-nominal_voltage)>0.001) + { + gl_error("pypower load bus nominal voltage (Vn %.1f V) and substation (nominal_voltage %.1f V) do not match to within 0.001 V",*pTransNominalVoltage,nominal_voltage); + return 0; + } + if (bustype != SWING) + { + warning("substation attached to pypower load and not a SWING bus - forcing bustype SWING"); + bustype = SWING; + } + has_parent = 1; } else //Parent isn't a pw_load, so we just become a normal node - let it handle things { diff --git a/module/pypower/autotest/case14.glm b/module/pypower/autotest/case14.glm index 0b0dba84c..7f7225daa 100644 --- a/module/pypower/autotest/case14.glm +++ b/module/pypower/autotest/case14.glm @@ -4,7 +4,11 @@ module pypower version 2; baseMVA 100.0; } -object bus + +// +// bus +// +object pypower.bus { bus_i 1.0; type 3.0; @@ -20,7 +24,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 2.0; type 2.0; @@ -36,7 +40,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 3.0; type 2.0; @@ -52,7 +56,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 4.0; type 1.0; @@ -68,7 +72,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 5.0; type 1.0; @@ -84,7 +88,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 6.0; type 2.0; @@ -100,7 +104,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 7.0; type 1.0; @@ -116,7 +120,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 8.0; type 2.0; @@ -132,7 +136,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 9.0; type 1.0; @@ -148,7 +152,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 10.0; type 1.0; @@ -164,7 +168,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 11.0; type 1.0; @@ -180,7 +184,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 12.0; type 1.0; @@ -196,7 +200,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 13.0; type 1.0; @@ -212,7 +216,7 @@ object bus Vmax 1.06; Vmin 0.94; } -object bus +object pypower.bus { bus_i 14.0; type 1.0; @@ -228,7 +232,11 @@ object bus Vmax 1.06; Vmin 0.94; } -object gen + +// +// gen +// +object pypower.gen { bus 1.0; Pg 232.4; @@ -252,7 +260,7 @@ object gen ramp_q 0.0; apf 0.0; } -object gen +object pypower.gen { bus 2.0; Pg 40.0; @@ -276,7 +284,7 @@ object gen ramp_q 0.0; apf 0.0; } -object gen +object pypower.gen { bus 3.0; Pg 0.0; @@ -300,7 +308,7 @@ object gen ramp_q 0.0; apf 0.0; } -object gen +object pypower.gen { bus 6.0; Pg 0.0; @@ -324,7 +332,7 @@ object gen ramp_q 0.0; apf 0.0; } -object gen +object pypower.gen { bus 8.0; Pg 0.0; @@ -348,7 +356,11 @@ object gen ramp_q 0.0; apf 0.0; } -object branch + +// +// branch +// +object pypower.branch { fbus 1.0; tbus 2.0; @@ -364,7 +376,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 1.0; tbus 5.0; @@ -380,7 +392,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 2.0; tbus 3.0; @@ -396,7 +408,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 2.0; tbus 4.0; @@ -412,7 +424,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 2.0; tbus 5.0; @@ -428,7 +440,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 3.0; tbus 4.0; @@ -444,7 +456,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 4.0; tbus 5.0; @@ -460,7 +472,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 4.0; tbus 7.0; @@ -476,7 +488,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 4.0; tbus 9.0; @@ -492,7 +504,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 5.0; tbus 6.0; @@ -508,7 +520,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 6.0; tbus 11.0; @@ -524,7 +536,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 6.0; tbus 12.0; @@ -540,7 +552,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 6.0; tbus 13.0; @@ -556,7 +568,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 7.0; tbus 8.0; @@ -572,7 +584,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 7.0; tbus 9.0; @@ -588,7 +600,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 9.0; tbus 10.0; @@ -604,7 +616,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 9.0; tbus 14.0; @@ -620,7 +632,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 10.0; tbus 11.0; @@ -636,7 +648,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 12.0; tbus 13.0; @@ -652,7 +664,7 @@ object branch angmin -360.0; angmax 360.0; } -object branch +object pypower.branch { fbus 13.0; tbus 14.0; @@ -668,3 +680,42 @@ object branch angmin -360.0; angmax 360.0; } + +// +// gencost +// +object pypower.gencost +{ + model 2; + startup 0.0; + shutdown 0.0; + costs "0.0430293,20.0,0.0"; +} +object pypower.gencost +{ + model 2; + startup 0.0; + shutdown 0.0; + costs "0.25,20.0,0.0"; +} +object pypower.gencost +{ + model 2; + startup 0.0; + shutdown 0.0; + costs "0.01,40.0,0.0"; +} +object pypower.gencost +{ + model 2; + startup 0.0; + shutdown 0.0; + costs "0.01,40.0,0.0"; +} +object pypower.gencost +{ + model 2; + startup 0.0; + shutdown 0.0; + costs "0.01,40.0,0.0"; +} diff --git a/module/pypower/autotest/test_case14_substation.glm b/module/pypower/autotest/test_case14_substation.glm index 8535e3421..e0a5f840d 100644 --- a/module/pypower/autotest/test_case14_substation.glm +++ b/module/pypower/autotest/test_case14_substation.glm @@ -18,10 +18,9 @@ object pypower.load object substation { phases ABC; + bustype SWING; nominal_voltage 12.5 kV; base_power ${pypower::baseMVA}; power_convergence_value 0.01; }; } - - diff --git a/module/pypower/gencost.cpp b/module/pypower/gencost.cpp index c18c4d353..84afeabcf 100644 --- a/module/pypower/gencost.cpp +++ b/module/pypower/gencost.cpp @@ -22,70 +22,24 @@ gencost::gencost(MODULE *module) defaults = this; if (gl_publish_variable(oclass, - PT_int32, "gencost_i", get_gencost_i_offset(), - PT_DESCRIPTION, "gencost number (1 to 29997)", + PT_enumeration, "model", get_model_offset(), + PT_DESCRIPTION, "cost model (1=piecewise linear, 2=polynomial)", + PT_KEYWORD, "UNKNOWN", (enumeration)CM_UNKNOWN, + PT_KEYWORD, "PIECEWISE", (enumeration)CM_PIECEWISE, + PT_KEYWORD, "POLYNOMIAL", (enumeration)CM_POLYNOMIAL, - PT_enumeration, "type", get_type_offset(), - PT_DESCRIPTION, "gencost type (1 = PQ, 2 = PV, 3 = ref, 4 = isolated)", - PT_KEYWORD, "UNKNOWN", (enumeration)0, - PT_KEYWORD, "PQ", (enumeration)1, - PT_KEYWORD, "PV", (enumeration)2, - PT_KEYWORD, "REF", (enumeration)3, - PT_KEYWORD, "NONE", (enumeration)4, + PT_double, "startup[$]", get_startup_offset(), + PT_DESCRIPTION, "startup cost ($)", - PT_double, "Pd[MW]", get_Pd_offset(), - PT_DESCRIPTION, "real power demand (MW)", + PT_double, "shutdown[$]", get_shutdown_offset(), + PT_DESCRIPTION, "shutdown cost($)", - PT_double, "Qd[MVAr]", get_Qd_offset(), - PT_DESCRIPTION, "reactive power demand (MVAr)", + PT_char1024, "costs", get_costs_offset(), + PT_DESCRIPTION, "cost model (comma-separate values)", - PT_double, "Gs[MW]", get_Gs_offset(), - PT_DESCRIPTION, "shunt conductance (MW at V = 1.0 p.u.)", - - PT_double, "Bs[MVAr]", get_Bs_offset(), - PT_DESCRIPTION, "shunt susceptance (MVAr at V = 1.0 p.u.)", - - PT_int32, "area", get_area_offset(), - PT_DESCRIPTION, "area number, 1-100", - - PT_double, "baseKV[kV]", get_baseKV_offset(), - PT_DESCRIPTION, "voltage magnitude (p.u.)", - - PT_double, "Vm[pu*V]", get_Vm_offset(), - PT_DESCRIPTION, "voltage angle (degrees)", - - PT_double, "Va[deg]", get_Va_offset(), - PT_DESCRIPTION, "base voltage (kV)", - - PT_int32, "zone", get_zone_offset(), - PT_DESCRIPTION, "loss zone (1-999)", - - PT_double, "Vmax", get_Vmax_offset(), - PT_DESCRIPTION, "maximum voltage magnitude (p.u.)", - - PT_double, "Vmin", get_Vmin_offset(), - PT_DESCRIPTION, "minimum voltage magnitude (p.u.)", - - PT_double, "lam_P", get_lam_P_offset(), - PT_DESCRIPTION, "Lagrange multiplier on real power mismatch (u/MW)", - PT_ACCESS, PA_REFERENCE, - - PT_double, "lam_Q", get_lam_Q_offset(), - PT_DESCRIPTION, "Lagrange multiplier on reactive power mismatch (u/MVAr)", - PT_ACCESS, PA_REFERENCE, - - PT_double, "mu_Vmax", get_mu_Vmax_offset(), - PT_DESCRIPTION, "Kuhn-Tucker multiplier on upper voltage limit (u/p.u.)", - PT_ACCESS, PA_REFERENCE, - - PT_double, "mu_Vmin", get_mu_Vmin_offset(), - PT_DESCRIPTION, "Kuhn-Tucker multiplier on lower voltage limit (u/p.u.)", - PT_ACCESS, PA_REFERENCE, - - NULL)<1){ - char msg[256]; - snprintf(msg,sizeof(msg)-1, "unable to publish properties in %s",__FILE__); - throw msg; + NULL)<1) + { + throw "unable to publish properties in pypower gencost"; } } } @@ -108,5 +62,23 @@ int gencost::create(void) int gencost::init(OBJECT *parent) { + if ( model == CM_UNKNOWN ) + { + error("cost model must be PIECEWISE or POLYNOMIAL"); + return 0; + } + + if ( startup < 0 ) + { + error("startup cost must be non-negative"); + return 0; + } + + if ( shutdown < 0 ) + { + error("shutdown cost must be non-negative"); + return 0; + } + return 1; } diff --git a/module/pypower/gencost.h b/module/pypower/gencost.h index febe127a3..fabdcb930 100644 --- a/module/pypower/gencost.h +++ b/module/pypower/gencost.h @@ -8,26 +8,15 @@ class gencost : public gld_object { +private: + typedef enum {CM_UNKNOWN=0,CM_PIECEWISE=1,CM_POLYNOMIAL} COSTMODEL; public: // published properties - GL_ATOMIC(int32,gencost_i); - GL_ATOMIC(enumeration,type); - GL_ATOMIC(double,Pd); - GL_ATOMIC(double,Qd); - GL_ATOMIC(double,Gs); - GL_ATOMIC(double,Bs); - GL_ATOMIC(int32,area); - GL_ATOMIC(double,baseKV); - GL_ATOMIC(double,Vm); - GL_ATOMIC(double,Va); - GL_ATOMIC(int32,zone); - GL_ATOMIC(double,Vmax); - GL_ATOMIC(double,Vmin); - GL_ATOMIC(double,lam_P); - GL_ATOMIC(double,lam_Q); - GL_ATOMIC(double,mu_Vmax); - GL_ATOMIC(double,mu_Vmin); + GL_ATOMIC(enumeration,model); + GL_ATOMIC(double,startup); + GL_ATOMIC(double,shutdown); + GL_ATOMIC(char1024,costs); public: // event handlers diff --git a/module/pypower/pypower.cpp b/module/pypower/pypower.cpp index 49c772af7..90e53ae09 100644 --- a/module/pypower/pypower.cpp +++ b/module/pypower/pypower.cpp @@ -133,10 +133,12 @@ EXPORT bool on_init(void) if ( enable_opf ) { - // TODO: required for OPF solution - throw "OPF not supported yet"; - // gencostdata = PyList_New(7); - // PyDict_SetItemString(data,"gencost",gencostdata); + gencostdata = PyList_New(ngencost); + PyDict_SetItemString(data,"gencost",gencostdata); + for ( size_t n = 0; n < ngencost ; n++ ) + { + PyList_SetItem(gencostdata,n,PyList_New(4)); + } } return true; @@ -223,6 +225,15 @@ EXPORT TIMESTAMP on_sync(TIMESTAMP t0) PyList_SetItem(pyobj,24,PyFloat_FromDouble(obj->get_mu_Qmin())); } } + for ( size_t n = 0 ; n < ngencost ; n++ ) + { + gencost *obj = gencostlist[n]; + PyObject *pyobj = PyList_GetItem(gencostdata,n); + PyList_SetItem(pyobj,0,PyLong_FromLong(obj->get_model())); + PyList_SetItem(pyobj,1,PyFloat_FromDouble(obj->get_startup())); + PyList_SetItem(pyobj,2,PyFloat_FromDouble(obj->get_shutdown())); + PyList_SetItem(pyobj,3,PyUnicode_FromString(obj->get_costs())); + } // run solver PyObject *result = PyObject_CallOneArg(solver,data); diff --git a/module/pypower/pypower_solver.py b/module/pypower/pypower_solver.py index 280e6ddab..c392b0b8c 100644 --- a/module/pypower/pypower_solver.py +++ b/module/pypower/pypower_solver.py @@ -2,10 +2,9 @@ # Copyright (C) 2024 Regents of the Leland Stanford Junior University import os, sys -from pypower.api import case14, ppoption, runpf +from pypower.api import case14, ppoption, runpf, runopf from numpy import array -with_opf = False save_case = False debug = False verbose = False @@ -49,23 +48,26 @@ def solver(pf_case): # copy from model for name in ['bus','gen','branch']: - casedata[name] = array(pf_case[name]) + if name in pf_case: + casedata[name] = array(pf_case[name]) + if 'gencost' in pf_case: + costdata = [] + for cost in pf_case['gencost']: + costs = [float(x) for x in cost[3].split(',')] + costdata.append([int(cost[0]),cost[1],cost[2],len(costs)]) + costdata[-1].extend(costs) + casedata['gencost'] = array(costdata) + if save_case: with open("pypower_casedata.py","w") as fh: fh.write(str(casedata)) - # TODO: call solver - # print(f"solver(pf_case={pf_case})",file=sys.stderr,flush=True) - # stdout = sys.stdout - # stderr = sys.stderr - # devnull = open("/dev/null","w") - # sys.stdout = devnull - # sys.stderr = devnull - results,success = runpf(casedata,options) - # sys.stdout = stdout - # sys.stderr = stderr - # devnull.close() + if 'gencost' in casedata: + results = runopf(casedata,options) + success = results['success'] + else: + results,success = runpf(casedata,options) if save_case: with open("pypower_results.py","w") as fh: @@ -84,9 +86,13 @@ def solver(pf_case): # print(" --> FAILED:",results,file=sys.stderr,flush=True) return False - except Exception as err: + except Exception: - print(" --> EXCEPTION:",err,file=sys.stderr,flush=True) + e_type,e_value,e_trace = sys.exc_info() + print(" --> EXCEPTION:",e_type,e_value,file=sys.stderr,flush=True) + if debug: + import traceback + traceback.print_tb(e_trace,file=sys.stderr) return False