SymEnergy

Symbolic modelling of energy systems, based on SymPy.

SymEnergy provides a framework for the structured solution of cost-optimization of energy systems using the standard Lagrange multiplier approach. The result consists in the close-form analytical solutions to the optimization problem. For example, the energy production from a power plant is expressed as a function of the vector of symbolic parameters. These solutions are evaluated for certain parameters in order to identify relevant constraint combinations.

Installation

pip install symenergy

Repository

https://github.com/mcsoini/symenergy

Publications

Content

Minimal example

Investigate the dispatch of two power plants during a single time slot in dependence on decreasing residual load (increasing variable renewable energy production).

Initialize a SymEnergy model instance:

from symenergy.core import model
m = model.Model()

Add a time slot 't0' with 3000 load and 1 VRE production. VRE production is varied later; any finite value of this parameter is acceptable at this point to initialize the parameter.

m.add_slot('t0', load=3000, vre=1)

Add two power plants: Cheap baseload power 'n' with limited capacity, expensive peaker plants 'g' with unconstrained power production:

m.add_plant('n', vc0=1, capacity=1000)
m.add_plant('g', vc0=2)

All Lagrange problems corresponding to this model are generated and solved.

m.generate_solve()

All results of this process are stored in the model attribute df_comb, indexed by the constraint combinations (see below). The table below shows a single line of this table. It corresponds to the case where

m.df_comb.iloc[[0]].T

===================  ==================================================================================================================================================
act_lb_n_pos_p_t0    True
act_lb_n_p_cap_C_t0  False
act_lb_g_pos_p_t0    False
lagrange             g_p_t0*vc0_g_none*w_none + lb_n_pos_p_t0*n_p_t0 + n_p_t0*vc0_n_none*w_none + pi_supply_t0*w_none*(-g_p_t0 + l_t0 - n_p_t0 - vre_scale_none*vre_t0)
variabs_multips      [g_p_t0, lb_n_pos_p_t0, n_p_t0, pi_supply_t0]
result               (-(-l_t0*w_none + vre_scale_none*vre_t0*w_none)/w_none, -(-vc0_g_none*w_none**2 + vc0_n_none*w_none**2)/w_none, 0, vc0_g_none)
tc                   -vc0_g_none*(-l_t0*w_none + vre_scale_none*vre_t0*w_none)
===================  ==================================================================================================================================================

The lagrange value in the table above corresponds to the Lagrange function constructed for this combination of active and inactive constraints. It is equal to the total cost (g_p_t0*vc0_g_none*w_none + n_p_t0*vc0_n_none*w_none), the equality constraint (energy balance/supply constraint -g_p_t0 + l_t0 - n_p_t0 - vre_scale_none*vre_t0 with multiplier (shadow price) pi_supply_t0), plus the only binding constraint with the corresponding multiplier lb_n_pos_p_t0*n_p_t0.

The variabs_multips are a list of all symbols/variables for which the problem was solved. In this case, these are the peaker plant power output g_p_t0, the shadow price of the baseload power posivity constraint (lb_n_pos_p_t0), the baseload power production n_p_t0, and the shadow price of the supply constraint (the “electricity price”) pi_supply_t0. The result column contains a list of the same variables and multipliers. Not surprisingly, the solution (not fully simplified) of the peaker power production is equal to the residual load l_t0 - vre_scale_none*vre_t0; the electricity price is equal to the variable cost of peaker power production vc0_g_none. The column tc is the total cost, generated by substituting the variable solutions into the total cost expression.

The SymEnergy evaluator calculates numerical results from the closed-form solutions for certain parameter values. This allows the identify the relevant constraint combinations. In the example below the model is evaluated for increasing VRE production: The value of the parameter vre of the time slot 't0' is varied in 31 steps between 0 and 3000 (residual load 0).

from symenergy.evaluator import evaluator
import numpy as np

x_vals_vre = {m.slots['t0'].vre: np.linspace(0, 3000, 31)}
ev = evaluator.Evaluator(m, x_vals=x_vals_vre)
ev.get_evaluated_lambdas_parallel()
ev.expand_to_x_vals_parallel()

The final results are stored in the attribute ev.df_exp. The dataframe ev.df_bal is a convenient way of accessing the optimal solution’s energy balance:

(ev.df_bal.query('slot not in ["global"]')
   .pivot_table(index='vre_t0', columns='func',
                   values='lambd')[['n_p_t0', 'g_p_t0', 'vre_t0', 'l_t0']].plot.area())
minimal balance

The optimal solution thus corresponds to the case where the cheap baseload power plant produces power at maximum output, while the peaker plants cover the remaining residual load. Once the VRE production is large enough to fully replace the peaker plants, the production from baseload plants is reduced.

Additional analyses can be performed using the ev.df_exp table, which contains all evaluated variables. For example, we can identify for which ranges of the varied parameter vre_t0 the individual constraint combinations are active:

(ev.df_exp.query('is_optimum and func == "tc"')
        .groupby(['idx'] + m.constrs_cols_neq)
        .vre_t0.apply(lambda x: (x.min(), x.max()))
        .reset_index())
idx act_lb_n_pos_p_t0 act_lb_n_p_cap_C_t0 act_lb_g_pos_p_t0 vre_t0
1 True False False (3000.0, 3000.0)
3 False True False (0.0, 2000.0)
4 False False True (2100.0, 2900.0)

See also

Section Theoretical discussion of the minimal example
Theoretical description of the minimal example’s solution process

Theoretical description of the solution process

The figure below provides an overview of the solution process:

flowchart

Components and their constraints and costs

Time slots

Time slots \(t\) are characterized by parameters representing the demand level \(\Phi_{\mathrm{Load},t}\) and the VRE production \(\Phi_{\mathrm{VRE},t}\). Optionally, a power variable \(p_{\mathrm{Curt},t}\) allows for surplus production. These variables must be greater than zero, i.e. they are subject to positivity constraintsfootnote{Equality constraint expressions are denoted by \(\Pi\), inequality constraint expressions are represented by the symbol \(\Lambda\).} \(\Lambda_{\mathrm{pos}, p_{\mathrm{Curt},t}}\geqslant 0\) with \(\Lambda_{\mathrm{pos}, p_{\mathrm{Curt},t}} = p_{\mathrm{Curt},t}\).

Power plants

Power plants pp produce electric power (variable \(p_{\mathrm{pp},t}\)) at a linear specific variable cost \(c_{\mathrm{v},t} = \gamma_\mathrm{v,0,pp} + \gamma_\mathrm{v,1,pp}\cdot p_{\mathrm{pp},t}\) with parameters \(\gamma_\mathrm{v,0,pp}\) and \(\gamma_\mathrm{v,1,pp}\) (compare section ref{sect_method_dispatch}). Optional parameters are fixed O&M costs \(\gamma_\mathrm{f,pp}\) (EUR/MW) and installed capacities \(\Gamma_\mathrm{pp}\) (MW). The optional retired capacity variable \(C_\mathrm{pp,ret}\) (MW) is subject to optimization. The total cost associated with power plant operation is thus

\[\begin{split}c_\mathrm{Total} = \sum_{\mathrm{pp}, t}\big[ & p_{\mathrm{pp},t}\left( \gamma_\mathrm{v,0,pp} + 0.5\cdot \gamma_\mathrm{v,1,pp} \cdot p_{\mathrm{pp},t} \right) \\ & + \gamma_\mathrm{f,pp} \left(\Gamma_\mathrm{pp} - C_\mathrm{pp,ret}\right) \big]. \label{equ_total_cost}\end{split}\]

The feasible solution space of the power plant variables is reduced by the capacity constraints \(\Lambda_{\mathrm{cap},p_{\mathrm{pp},t}} = \left(\Gamma_\mathrm{pp} - C_\mathrm{pp,ret}\right) - p_{\mathrm{pp},t}\geqslant 0\) and \(\Lambda_{\mathrm{cap},C_\mathrm{pp,ret}} = \Gamma_\mathrm{pp} - C_\mathrm{pp,ret}\geqslant 0\), as well as the positivity constraints \(\Lambda_{\mathrm{pos},p_{\mathrm{pp},t}} = p_{\mathrm{pp},t} \geqslant 0\) and \(\Lambda_{\mathrm{pos},C_\mathrm{pp,ret}} = C_\mathrm{pp,ret} \geqslant 0\).

Storage

Storage assets s charge and discharge with power \(p_{\mathrm{chg/dch,s},t}\). The round-trip efficiency is \(\eta\). Optionally, the subsets of time slots during which the storage is charging and/or discharging can be selected. This allows for model simplifications. The internal structure of the storage representation depends on the amount of simplification:

  • Strong simplification: Assume a two time slot system with \(\mathcal{H}={t_0, t_1}\). If the charging and discharging slots are defined in such a way that the two corresponding slot subsets are disjoint \(\mathcal{H}_\mathrm{chg}\left\{t_0\right\}\), \(\mathcal{H}_\mathrm{dch}\left\{t_1\right\}\), the maximum stored energy \(e_\mathrm{s}\) is a time-independent variable. In this case, energy conservation is enforced through the charging-discharging equality constraint \(\Pi_{\mathrm{cd,s}} = \sum_t p_{\mathrm{chg,s},t} \eta_\mathrm{s} - \sum_\mathrm{t} p_{\mathrm{dch,s},t} = 0\) (round-trip efficiency \(\eta_\mathrm{s}\)); the charged energy follows from \(\Pi_\mathrm{e,s} = \sum_t p_{\mathrm{chg,s},t}\eta_\mathrm{s}^{1/2} - e_\mathrm{s} = 0\).
  • Weak simplification: For the same set of time slots, if charging-discharging sets overlap (e.g. \(\mathcal{H}_\mathrm{chg}=\left\{t_0\right\}\), \(\mathcal{H}_\mathrm{dch}=\left\{t_1\right\}\)), or no specifications are made, the stored energy \(e_{t, \mathrm{s}}\) is time-dependent. In this case, the stored energy is defined through the equality-constraint \(e_t = e_{t-1} - \eta^{-1/2} p_{\mathrm{dch,t} + \eta^{1/2} p_{\mathrm{chg},t}\).

All variables are subject to positivity constraints and optional capacity constraints (power capacity \(\Gamma_\mathrm{s}\) and energy capacity \(\Sigma_\mathrm{s}\)), analogous to the power plant variables.

Model constraints

For each \(t\), the demand equality constraint is defined as

\[\Pi_{\mathrm{d},t} = \Phi_{\mathrm{Load},t} + p_{\mathrm{Curt},t} + \sum_\mathrm{s}p_{\mathrm{chg,s},t} - \Phi_{\mathrm{VRE},t} - \sum_\mathrm{s}p_{\mathrm{dch,s},t} - \sum_\mathrm{pp}p_{\mathrm{pp}, t} = 0.\]

This describes the energy balance, i.e. all power consumption (load, curtailment, charging) must equal all power production (power plant production, storage discharging).

See also

Add model components
API documentation on adding components to a model instance

Theoretical discussion of the minimal example

The theoretical background of the basic example in section Minimal example is further explained. This example allows to demonstrate the optimization using SymEnergy. The considered system is simple enough for its optimal solution to be intuitively accessible. Energy is supplied during a single time slot (load parameter \(\Phi_\mathrm{Load}\), VRE production parameter \(\Phi_\mathrm{VRE}\)) by cheap baseload power and expensive peaker plants:

  • Base plant production \(p_\mathrm{Base}\) (specific variable cost \(\gamma_\mathrm{Base}\)) must be positive (constraint \(\Lambda_{\mathrm{pos},p_{\mathrm{Base}}}\geqslant 0\), multiplier \(\lambda_{\mathrm{pos},p_{\mathrm{Base}}}\)); it is constrained by the base capacity \(\Gamma_\mathrm{Base}\) (\(\Lambda_{\mathrm{cap},p_{\mathrm{Base}}}\geqslant 0\) with multiplier \(\lambda_{\mathrm{cap},p_\mathrm{Base}}\)).
  • Peaker plant production \(p_\mathrm{Peak}\) (specific variable cost \(\gamma_\mathrm{Peak} > \gamma_\mathrm{Base}\)) must be positive (\(\Lambda_{\mathrm{pos},p_{\mathrm{Peak}}}\geqslant 0\), multiplier \(\lambda_{\mathrm{pos},p_{\mathrm{Peak}}}\)) and is otherwise unconstrained.

The total cost is \(p_\mathrm{Base}\gamma_\mathrm{Base} + p_\mathrm{Peak}\gamma_\mathrm{Peak}\), the demand constraint is \(\Pi_{\mathrm{d}}= \Phi_\mathrm{Load}-p_\mathrm{Base}-p_\mathrm{Peak} - \Phi_\mathrm{VRE} = 0\) with multiplier \(\pi_\mathrm{d}\).

example_toy

Without further analysis, the optimal solution can be expected to correspond to the case where \(p_\mathrm{Base}\) covers a maximum of the residual load \(\Phi_\mathrm{Load} - \Phi_\mathrm{VRE}\), with \(p_\mathrm{Peak}\) supplying the rest. This is shown in panel (a) of the figure above as a function of the VRE share. SymEnergy starts from the 3 binding or non-binding inequality constraints. Consequently, the total number of potential general solutions (constraint combinations) is \(2^3=8\):

  • Two of these constraint combinations—\((\Lambda_{\mathrm{pos},p_{\mathrm{Base}}},\Lambda_{\mathrm{cap},p_{\mathrm{Base}}},\Lambda_{\mathrm{pos},p_{\mathrm{Peak}}})\in \left\{\left(0,0,\mathbb{R}_{\neq0}\right),\left(0,0,0\right)\right\}\)—correspond to the case where baseload power production is simultaneously zero and equal its installed power capacity. They can be excluded a priori, without loss of generality.
  • Three of the remaining constraints would be valid only for specific parameter value combinations. For example, all constraint combinations with zero peak production and binding baseload capacity constraint (\(\Lambda_{\mathrm{pos},p_{\mathrm{Peak}}}=\Lambda_{\mathrm{cap},p_{\mathrm{Base}}}=0\)), would require the baseload capacity parameter to be exactly equal the residual load parameter: \(\Gamma_\mathrm{Base} \equiv \Phi_\mathrm{Load} - \Phi_\mathrm{VRE}\). Therefore, no general solution exists.
  • The remaining three combinations
\[\begin{split} (\Lambda_{\mathrm{pos},p_{\mathrm{Base}}}, \Lambda_{\mathrm{cap},p_{\mathrm{Base}}}, \Lambda_{\mathrm{pos},p_{\mathrm{Peak}}}) \in \{ &(\mathbb{R}_{\neq0},\mathbb{R}_{\neq0},0),\\ &(\mathbb{R}_{\neq0},0,\mathbb{R}_{\neq0}),\\ &(0,\mathbb{R}_{\neq0},\mathbb{R}_{\neq0})\}\end{split}\]

are potentially optimal solutions, depending on the parameter values.

For illustration, the solution for the constraint combination \((\mathbb{R}_{\neq0},0,\mathbb{R}_{\neq0})\) is calculated. It implies a binding base capacity constraint (\(\Lambda_{\mathrm{cap},p_\mathrm{Base}} = \Gamma_\mathrm{Base} - p_\mathrm{Base} = 0\)) and both plants producing positive power output (\(\Lambda_{\mathrm{pos},p_\mathrm{Base}} = p_\mathrm{Base} > 0\), \(\Lambda_{\mathrm{pos},p_\mathrm{Peak}} = p_\mathrm{Peak} > 0\)). The resulting Lagrange function is

\[\begin{split}\mathcal{L}_{(\mathbb{R}_{\neq0},0,\mathbb{R}_{\neq0})} & = p_\mathrm{Base}\gamma_\mathrm{Base} + p_\mathrm{Peak}\gamma_\mathrm{Peak} \\ &+ \lambda_{\mathrm{cap},p_\mathrm{Base}}\Lambda_{\mathrm{cap},p_\mathrm{Base}}+ \pi_\mathrm{d}\Pi_\mathrm{d}\end{split}\]

The KKT conditions \(\nabla_{\mathbf{v},\mathbf{\pi},\mathbf{\lambda}}\mathcal{L}_{(\mathbb{R}_{\neq0},0,\mathbb{R}_{\neq0})}(\mathbf{p},\mathbf{v},\mathbf{\pi},\mathbf{\lambda})=0\) with respect to all variables and multipliers yield a system of linear equations:

\[\begin{split}\left[\begin{matrix} 0 & 0 & 1 & 1\\ 0 & 0 & 1 & 0\\ 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0 \end{matrix}\right] \cdot \left[\begin{matrix} p_\mathrm{Base} \\ p_\mathrm{Peak} \\ \pi_\mathrm{d} \\ \lambda_{\mathrm{cap},p_\mathrm{Base}} \\ \end{matrix}\right] = \left[\begin{matrix} \gamma_\mathrm{Base}\\ \gamma_\mathrm{Peak}\\ \Gamma_\mathrm{Base}\\ \Phi_\mathrm{Load} - \Phi_\mathrm{VRE} \\ \end{matrix}\right].\end{split}\]

Its closed-form analytical solution \(p_\mathrm{Base} = \Gamma_\mathrm{Base}\), \(p_\mathrm{Peak} = \Phi_\mathrm{Load} - \Phi_\mathrm{VRE} - \Gamma_\mathrm{Base}\), \(\pi_\mathrm{d} =\gamma_\mathrm{Peak}\), \(\lambda_{\mathrm{cap},p_\mathrm{Base}} = \gamma_\mathrm{Base} - \gamma_\mathrm{Peak}\) corresponds to the case where the baseload plants are operating at the capacity limit, and the peaker plant is used to cover the remaining residual load. With the respective definitions of the other Lagrange functions \(\mathcal{L_\mathrm{CC}}\), the general solutions for the remaining constraint combinations above are found analogously.

In the figure above, b-d, the solutions of all constraint combinations are evaluated for 51 discrete VRE production values from 0 to \(\Phi_\mathrm{Load}\).

  • Optimal: By definition, the optimal solution for each set of parameter values is determined by the feasible solution with the lowest total cost (panel b). For low VRE shares, this implies maximum baseload operation—\((\mathbb{R}_{\neq0},0, \mathbb{R}_{\neq0})\); meanwhile, peak plants cover the remaining power demand. For higher VRE shares, the residual demand is smaller than the base capacity, such that the output from the latter gradually decreases; in this case, the peak output is zero and the base plant is not capacity constrained—\((\mathbb{R}_{\neq0},\mathbb{R}_{\neq0},0)\).
  • Non-optimal: Another feasible solution exists for all VRE shares—\((0, \mathbb{R}_{\neq0},\mathbb{R}_{\neq0})\): It corresponds to the base plant delivering zero output and all demand being covered by the peaker plants (panel c). The higher variable cost of the latter makes this option prohibitively expensive in all cases.
  • Infeasible: Finally, the constraint combinations comprising the optimal solutions can also be infeasible, depending on the VRE share: At high VRE shares, this is the case if the base plants produce at full power output and peaker plants produce negative output to compensate for the overproduction. At low VRE shares, base plant capacity violations are equally cost-effective yet infeasible. For any given parameter set, these infeasible solutions must be identified by evaluating all of the model’s constraints.

See also

Minimal example
Description of the implementation of this example in SymEnergy
Components and their constraints and costs
Discussion of model component properties
https://doi.org/10.1016/j.eneco.2019.104495
This example is included in the publication Does bulk electricity storage assist wind and solar in replacing dispatchable power production?, Energy Economics, Volume 85, January 2020, 104495

Indices and tables