Background

This section gives a brief overview of the concept behind COMANDO, the problem formulations that can be addressed with it, as well as its usage.

Concept

The aim of COMANDO is to enable object-oriented modeling of energy systems and their constituting components for the purpose of creating optimization problems related to the system’s design and operation. Unlike other energy system modeling frameworks, COMANDO does not impose restrictions to maintain the resulting optimization problems in a particular class such as linear programming (LP) or mixed-integer linear programming (MILP). Instead, component and system models may contain a wide range of nonlinear expressions, including nonconvex and nonsmooth functions. Furthermore ordinary differential equations describing system dynamics can be included directly, without manual discretization. In this way COMANDO allows users to create component and system models in a natural fashion, without having to obscure them with optimization-specific implementation details.

By specifying the connections between components, different component models can be combined into a system model. Based on a finished system model, different optimization problems can be created by specifying objective functions, the operational scenarios and their time-structure to be considered, as well as associated parameter data.

To bring the resulting problem formulations into a form that solvers understand, certain reformulations may be necessary. COMANDO provides routines for substitution and reformulation of expressions, as well as for automatic approximation via linearization and time-discretization. Once an acceptable problem formulation is created, it can be translated to different algebraic modeling languages (AMLs) or be passed directly to an appropriate solver for solution.

The details of the different steps of this process are explained in the following.

Usage

The usage of COMANDO can be split into three phases

  1. Modeling phase

    • Component model creation

    • System model creation

  2. Problem formulation phase

    • problem generation

      • objective selection

      • time-structure selection

      • scenario-structure selection

      • providing data

    • problem reformulation

      • time discretization

      • linearization

      • …

  3. Problem solution phase

    • via a solver interface (e.g., GUROBI, BARON)

    • via an algebraic modeling language (AML) interface (e.g., GAMS, PYOMO)

    • via custom algorithm

      • initialization

      • decomposition

      • …

Modeling phase

In the modeling phase the system behavior is specified in terms of the component behavior and system structure.

Components

In COMANDO a component is the basic building block of an energy system. It represents a model of a generic real-world component, specified via a collection of mathematical expressions that are given in a symbolic form.

Take for example the following simple description of some generic component C:

\begin{align} C_{\mathrm{out}} &= C_{\mathrm{in}} C_{\eta} \\ C_{\eta} &= (C_0 + C_1 \, C_{\mathrm{pl}} + C_2 \, C_{\mathrm{pl}}^2) \\ C_{\mathrm{pl}} &= C_{\mathrm{out}} / C_{\mathrm{out,max}} \end{align}

It contains different symbols that represent different design and operational quantities of C and equations that correlate these symbols with each other. Here, we may assume that \(C_{\mathrm{out,max}}\) is a design quantity, while the symbols \(C_{\mathrm{out}}\), \(C_{\mathrm{in}}\), \(C_{\eta}\), and \(C_{\mathrm{pl}}\) are operational quantities and \(C_0\), \(C_1\), and \(C_2\) are some numerical coefficients.

We can create a COMANDO model of this component by deciding which role the different symbols play for some design or operation process. In COMANDO a distinction is made between quantities that are given as input data and those that are to be determined by the use of the model, i.e., by formulating some optimization problem incorporating the model and solving that problem.

Symbols that act as placeholder for input data can be modeled as a comando.Parameter, while those that are to be determined can be modeled as a comando.Variable or a comando.VariableVector, depending on whether they represent design- or an operational-related quantity.

Note

Objects of type Variable / VariableVector always represent scalar / vector values while those of type Parameter may represent both, depending on the data that is assigned to them. In the context of design and operation, design quantities are always modeled by symbols or composed expressions that represent scalar values, while operational quantities are always modeled by symbols or composed expressions that represent vector values, with each entry representing the value of the corresponding quantity for a certain operational scenario and timestep.

We may for example model the coefficients \(C_0\), \(C_1\), and \(C_2\) as objects C_0, C_1, and C_2 of type Parameter, \(C_{\mathrm{out,max}}\) as a Variable C_out_max and \(C_{\mathrm{out}}\), \(C_{\mathrm{in}}\) as objects C_out, C_in of type VariableVector.

For the definition of component COMANDO provides the Component class. This class defines methods to create and add the corresponding symbols to Component instances:

# Define a new type of component
class MyComponent(comando.Component):
    """A generic component."""

    def __init__(self, label, *coeffs):
        """Initialize the component.

        Arguments
        ---------
        label : str
            The name of the component
        coeffs : Iterable[float]
            The polynomial coefficients
        """
        super().__init__(label)  # Initialize the parent class (Component)

        # Create parameters
        # >>> help(comando.Parameter)
        # NOTE: Parameters may hold time-constant or time-variables
        #       data, i.e., they may be `design' or `operational'
        #       parameters.
        #       Which type they are is not specified at the modeling
        #       stage, but only in the subsequent problem formulation
        #       stage.
        C_coeffs = [
            self.make_parameter(
              name=i,
              value=val,
              # Internal argument not intended for direct use
              # parent=None
            )
            for i, val in enumerate(coeffs)
        ]

        # Create design and operational variables
        # Define variables, also see
        # >>> help(comando.Variable)
        # for design variables
        # and
        # >>> help(comando.VariableVector)
        # for operational variables
        C_out_max = self.make_design_variable('out_max')
        C_out = self.make_operational_variable('out')
        C_in = self.make_operational_variable('in')

        # Expressions formed with the defined symbols

        # Expressions that are to be remembered can be stored
        C_pl = self.add_expression('part_load', C_out / C_out_max)
        # Intermediate expressions that don't need to be remembered
        # can also be created as temporary objects
        C_eta = sum(C_i * C_pl ** i for i, C_i in enumerate(C_coeffs))

        # Add constraints (left-hand-side, right-hand-side, identifier)
        self.add_eq_constraint(C_out, C_in * C_eta, 'conversion')
        self.add_le_constraint(C_out, C_out_max, 'output_limit')
        # ... or `self.add_ge_constraint(...)'
        # The resulting constraints can be accessed via the dict
        # self.constraints

        # Add connectors
        self.add_input('INPUT', C_in)
        self.add_output('OUTPUT', C_out)
        # also see
        # >>> help(comando.Component.add_connector)
        # >>> help(comando.Component.add_connectors)

# Create an instance of our new class
c = comando.Component('C')
c_pl = c.get_expression('part_load')  # Access previously defined expression

Note how instead of introducing additional VariableVector instances for \(C_{\eta}\) and \(C_{\mathrm{pl}}\), we directly used their definitions and create expressions for these quantities using overloaded Python functions! The expression for C_pl is considered interesting and is therefore stored in C under the name ‘part_load’. At a later point, this expression can be accessed via the Component.get_expression() method.

We may have similarly defined C_out as an expression, but instead we decided to make it an operational variable and introduce the correlation between input and output as an equality constraint. Similarly, we introduced an inequality constraint specifying that the component’s output must be smaller than it’s maximum value.

It is also possible to declare operational variables to be ‘states’, i.e., quantities whose time-derivative is given by some algebraic expression. This allows the consideration of dynamic effects. A typical example for this are storage units:

class Storage(Component):
    """General storage unit.

    Arguments
    ---------
    - label : str
        Unique sting that serves as an identifier of the storage unit.
    - max_cap: (None) upper bound for nominal size
    - charge_eff: (1) the charging efficiency of the storage unit
    - discharge_eff: (1) the discharging efficiency of the storage unit
    """

    def __init__(
        self,
        label,
        max_cap=None,
        charge_eff=1,
        discharge_eff=1,
    ):
        super().__init__(label)
        # design variables
        cap = self.make_design_variable(
            "cap", bounds=(0, max_cap), init_val=max_cap
        )
        # define operational variables
        soc = self.make_operational_variable("soc", bounds=(0, max_cap))
        self.add_le_constraint(soc, cap, name="cap_limit")

        # input at timestep t
        inp = self.make_operational_variable("inp", bounds=(0, max_cap))
        # output at timestep t
        out = self.make_operational_variable("out", bounds=(0, max_cap))
        # define storage constraints and declare state variables
        state_change = inp * charge_eff - out / discharge_eff

        # declare the operational variable a state
        # also see the alternative method
        # >>> help(comando.Component.make_state)
        self.declare_state(
            var=soc,
            rate_of_change=state_change,
            init_state=0,
            # Default values for derivative variable
            # der_bounds=(None, None),
            # der_init_val=None
        )

        # make sure, we don't charge and discharge at the same time
        self.add_eq_constraint(inp * out, 0, name="charge_or_discharge")

        self.add_expression("input", inp)
        self.add_expression("output", out)

        # set connectors
        self.add_input("INPUT", inp)
        self.add_output('OUTPUT', out)

In both examples above we defined ‘INPUT’ and ‘OUTPUT’ connectors that can be used to connect the component to others within a comando.System model. Components may assign individual algebraic expressions to connectors, allowing them to be interfaced with other components. Connectors may be specified as inputs, outputs or bidirectional connectors. The former two restrict the value of the corresponding expression to be nonnegative and nonpositive, respectively, while the latter does not impose additional restrictions.

Systems

Systems are modeled as collections of interconnected components, i.e., a system model consists of a set of components and the specification of how their connectors are connected. For this purpose COMANDO provides the comando.System class. As the System class inherits from the Component class, systems can define additional expressions, constraints and connectors and be incorporated as subsystems within other systems. As with the Component class, the System class may also be specialized via inheritance. For simple cases, however, we may also use it directly as in the example below:

# Instantiate some previously defined components
source = ...
component = MyComponent("C", 1.2, 3, 4)
storage = Storage('S')
consumer = ...

from comando import system
ES = System(
    'Energy_System',
    components=[source, component, storage, consumer],
    connections=dict(
      Source_Commodity_Bus=[source.OUTPUT, component.INPUT]
      Main_Commodity_Bus=[
          component.OUTPUT,
          storage.INPUT,
          storage.OUTPUT,
          consumer.INPUT
        ]
    )
)

Problem formulation phase

Given a system model, COMANDO can currently be used to create a comando.Problem object, representing a mathematical optimization problem (OP) of the form:

\[ % DEFINITIONS \renewcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \renewcommand{\min}[1][]{ {\underset{#1}{\text{min}\,}} } \renewcommand{\st}{\mathrm{s.\,t.}} \renewcommand{\dv}{\bm{x}} \renewcommand{\ov}{\bm{y}_{s}(\cdot)} % operational variables \renewcommand{\ovfunc}{\bm{y}_{s}(t)} % operational variables \renewcommand{\derov}{\dot{\bm{y}}^\mathrm{d}_s(t)} % derivative of ov \renewcommand{\roc}{\bm{f}} % specified rate of change for ov \renewcommand{\iv}{\bm{y}^\mathrm{d}_{s}(t=0)} % initial value \renewcommand{\design}{I} \renewcommand{\operation}{{II}} \renewcommand{\of}{F} % objective function \renewcommand{\dobj}{\of_\design} % design objective \renewcommand{\oobj}{\of_{\operation,s}} % operational objective \renewcommand{\moobj}{\dot{\of}_\operation} % momentary operational objective \renewcommand{\myargs}{\big(\dv, \ovfunc, \bm{p}_s(t)\big)} % \begin{alignedat}{4} & \min[\dv] \;\; && \dobj(\dv) + \sum_{s \in \mathcal{S}} w_s \, \oobj^*(\dv) \\[-1mm] & \underset{\hphantom{\dv, \bm{y}}}{\st} && \bm{g}_\design(\dv) \leq \bm{0} \\[-2mm] & && \bm{h}_\design(\dv) = \bm{0} \\ & && \left.\kern-0.75ex \begin{array}{lll} \oobj^*(\dv) = & \min[\bm{y}_s(\cdot)] & \oobj(\dv, \ov) = \!\! \int_{\mathcal{T}_s} \!\! \moobj\myargs \, \text{d}t \\ & \st & \bm{y}^\text{d}_{s}(\cdot) \in \bm{y}_{s}(\cdot) \\ & & \iv = \bm{y}^\text{d}_{s, 0} \\ & & \left.\kern-0.3ex \begin{array}{ll} \derov = \roc\myargs \!\!\!&\\ \bm{g}_\operation\myargs \leq \bm{0} &\\ \bm{h}_\operation\myargs = \bm{0} &\\ \bm{y}_{s}(t) = [\bm{y}^\text{d}_{s}(t), ...] \\ \bm{y}_{s}(t) \in \mathcal{Y}_{s}(t) \subset \mathbb{R}^{n_y} \!\times \mathbb{Z}^{m_y} \!\!\!\!\!\!\!\!\! \\ \end{array} \right\} \; \forall t \in \mathcal{T}_s \\ & & \mathcal{T}_s = \left[0, T_s\right] \\[-2mm] \end{array} \right\} \; \forall s \in \mathcal{S} \\[1mm] %= \{s_1, s_2, \cdots, s_{\card{\mathcal{S}}}\}\\ & && \dv \in \mathcal{X} \subset \mathbb{R}^{n_x} \!\times \mathbb{Z}^{m_x} \\[-2mm] & && \mathcal{S} = \{s_1, s_2, \cdots, s_{|\mathcal{S}|}\} \end{alignedat} \]

Where \(\dv\) and \(\ov\) are the vectors of design- and operation-variables, respectively.

Problem generation

\(F_I\) and \(F_{II}\) are user-specified scalar and indexed expressions corresponding to one-time and momentary costs, \(\mathcal{T}_s\) are time horizons for scenarios \(s\) in from the set \(\mathcal{S}\), with corresponding probabilities \(w_s\).

The constraints for the OP are automatically generated from the system model, i.e., all scalar relational expressions are taken as constraints and all indexed relational expressions are taken as constraints, parametrized by t, and s.

The dependence of the objective and constraint functions on time and scenario can be expressed in terms of the parameter values, which are user input. This data can currently be given only in discrete form, i.e., for a discrete time and scenario. The time-steps at which data is available

If any states were defined in the model the corresponding differential equations are currently discretized by default. An exception is the use of the pyomo_dae interface; here the time-continuous representation is passed and discretization via collocation can be performed. Data for the parameter values and initial guesses for the variable values can be provided based on the user-chosen sets T and S.

Note

The time- and scenario structure is only specified during Problem creation and is therefor not available during the modeling phase! While this may be unintuitive at first, it allows for a clean separation between component and system behavior on the one hand, and problem-specific implementation details on the other. An important consequence of this is, that certain aspects of component behavior such as, e.g., ramping constraints must either be defined in a more general way (e.g., via limiting the derivative of the ramped quantity), or must be deferred and carried out after Problem creation (e.g. via appropriate callbacks).

Given the System instance above, we can define a problem as follows:

# Access the design variable out_max from the component
# and use it as the design-part of the objective
d_obj = component["C_out_max"]

# Minimize required size of component for one week of operation
ts = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]
scenarios = {
    "low_demand": 0.2,
    "medium_demand": 0.7,
    "high_demand": 0.1
}
P = ES.create_problem(
    design_objective=d_obj,
    operational_objective=0,
    timesteps=(ts, 169),
    scenarios=scenarios,
    name='Example_Problem'
)

# Parameter values can be set as follows:
P["<name of some parameter>"] = values
# where values can be:
# - a scalar (for design parameters)
# - an iterable (e.g., numpy array, pandas Series) of length:
#   - len(ts) (for time-dependent, but scenario-independent data)
#   - len(scenarios) (for scenario-dependent, but time-independent data)
#   - len(ts) * len(scenarios) (for time- and scenario-dependent data)

For a more complete example, have a look at examples/example_system.ipynb

Todo

Examples for creating problem-specific constraints such as:

  • Constraints that restrict cumulative quantities (time-integrals of quantities)

  • Constraints that only consider a certain time range

  • Constraints explicitly referencing quantities corresponding to multiple time-points and/or scenarios

Problem reformulation

It is possible to use manual or automated reformulations of the original problem formulation.

Todo

Example reformulations, linearization module and default discretization scheme in the interfaces.

Problem solution

Given a Problem, the user can chose to pass it directly to a solver capable of handling the corresponding problem type, transform the COMANDO Problem to a representation in an AML (currently Pyomo or GAMS), or work directly with the COMANDO representation in a custom algorithm to preprocess or solve the Problem. Currently available interfaces as well as their installation requirements are described in Interfaces to solvers and algebraic modeling languages (AMLs).

If a solution is obtained, it is loaded back into the COMANDO Problem after the solver terminates.