Modeling
The following notebook provides a guide for how to model new power system entities (power flow, busses, components) in CommonPower.
[7]:
from typing import Callable, List
import pyomo.environ as pyo
from pyomo.core import Constraint, ConcreteModel, Expression
from commonpower.models.components import ESS
from commonpower.modeling.mip_builder import MIPExpressionBuilder
from commonpower.modeling import ModelElement, ElementTypes as et
from commonpower.modeling.robust_cost import CostScenario
from commonpower.modeling.robust_constraints import ConstraintScenario
Introduction
In CommonPower, every power system entity (power flow, busses, components) is described by a symbolic model. We use these models to compute optimal control actions or safety adjustments. The necessary optimization problems are automatically generated in the background and solved within the Python framework Pyomo. Pyomo defines models with variables, parameters, and constraints. We extend this classification to gain some nuance. In CommonPower, we use the following types of model elements (see below for examples):
INPUT - Input variable. Maps to the Pyomo type Var.
DATA - Exogenous input, which is read from a data provider. Maps to the Pyomo type Param.
STATE - State variable. Maps to the Pyomo type Var.
VAR - Generic variable. The difference to state variables is that VAR does not have to be initialized. Maps to the Pyomo type Var.
CONSTANT - Fixed Parameter. Parameters can either be constant across runs or be initialized in each run based on a specific logic. Maps to the Pyomo type Param.
COST - Cost variable. This is essentially a generic variable but explicitly defined to simplify downstream analysis. Maps to the Pyomo type Var.
CONSTRAINT - Constraint. Input coupling and dynamics functions are defined with this type. Maps to the Pyomo type Constraint.
SET - Set. Sets can be useful to specify the values a discrete variable can take. Maps to the Pyomo type Set.
We define dynamics and cost functions as constraints. These model elements are the basis of every model. Let’s have a look at the ESS model.
Note: The ‘M’ is the ‘big M constant’ used for mixed-integer constraints in optimization. Big-M constraints are typically used to propagate the implications of a binary, on-off decision to a continuous variable.
.info() method provides a fast overview over which model elements require configuration.[ ]:
ESS.info()
._get_model_elements() method of the ESS class, you will find that soc_init and cost are not defined.Node parent class defines an augmentation method which adds an init parameter for every defined state variable and one cost variable.._get_additional_constraints() method to keep ._get_model_elements() clean.MIPExpressionBuilder utility is used. It provides a convenient interface (based on the Big-M approach) to generate binary indicator variables.expbuilder.from_geq("p", 0, "p_ec").p_ec and the two indicator constraints you can see above. They make sure that p_ec == 1 iff p >= 0. We can then use p_ec in the dynamics function.Defining a custom ESS component
soc >= 0.8 * <soc upper bound>, and finally adjust the dynamic function accordingly.[9]:
def NewESS(ESS):
CLASS_INDEX = "new_ev"
def _get_additional_constraints(self) -> list[ModelElement]:
soc_rel = ModelElement(
"soc_rel",
et.VAR,
"relative state of charge",
bounds=(0, 1), # if we define it, the user will not be asked to provide bounds.
domain=pyo.NonNegativeReals,
)
def soc_rel_cf(scenario: ConstraintScenario):
def soc_rel_f(model, t):
scenario("soc_rel", model)[t] == scenario("soc", model, True)[t] / scenario("soc", model)[t].ub
return soc_rel_f
c_soc_rel = ModelElement(
"c_soc_rel",
et.ROBUST_CONSTRAINT,
"constraint to track the relative soc",
expr=soc_rel_cf,
)
mb = MIPExpressionBuilder(self)
mb.from_geq("p", 0, "p_ec")
mb.from_geq("soc_rel", 0.8, "soc80")
return [soc_rel, c_soc_rel] + mb.model_elements # the order here is important, as the elements are added to the pyomo model in this order
def _get_dynamic_fcn(self) -> list[ModelElement]:
def dynamics(scenario: ConstraintScenario):
def dynamic_fcn(model, t):
# p > 0 is charging
if t == self.horizon: # horizon+1 cannot have a constraint
return Constraint.Skip
else:
return (
scenario(self, "etas", model) * scenario(self, "soc", model)[t] * self.tau # self-discharge
) + (
#scenario(self, "etac", model)
(
scenario(self, "soc80") * 0.7 * scenario("etac", model) # if soc>=80% -> etac*0.7
+ (1 - scenario("soc80")) * scenario("etac", model) # rest similar to standard ESS component
)
* scenario("p", model)[t]
* scenario("p_ec", model)[t]
* self.tau
) + (
(1 / scenario("etad", model))
* scenario("p", model)[t]
* (1 - scenario("p_ec", model)[t])
* self.tau
) == scenario("soc", model, True)[t + 1]
return dynamic_fcn
dyn = ModelElement("dynamic_fcn", et.ROBUST_CONSTRAINT, "dynamic function", expr=dynamics)
return [dyn]
[10]:
def NewESS2(NewESS):
def cost_fcn(self, scenario: CostScenario, model: ConcreteModel, t: int = 0) -> Callable:
"""
The cost represents: cost = abs(p) * rho with an additional penalty factor of 1.2 if soc>=80%
"""
return (
(
scenario(self, "p", model)[t] * (scenario(self, "p_ec", model)[t])
- scenario(self, "p", model)[t] * (1 - scenario(self, "p_ec", model)[t])
)
* (
scenario(self, "soc80") * 1.2 * scenario(self, "rho", model)
+ (1 - scenario(self, "soc80")) * scenario(self, "rho", model)
)
* self.tau
)
Some notes
Node class also exposes an ._unmodeled_updates() method.CommonPower’s mixed-integer programming (MIP) approach is very flexible and allows modelling a large varieties of components. The use of indicator constraints generally produces quadratically constrained optimization problems, and often quadratically constrained quadratic programs (MIQCQPs). However, as most of the quadratic constraints are binary, our experience is that the resulting optimization problems can be solved reasonably fast. See the MIP tutorial for more information on this.