Optimal Power Flow

See the module pypsa.opf and pypsa.linopf. Optimisation with the linearised power flow equations for (mixed) AC and DC networks is fully supported. Note that optimisation with the full non-linear power flow equations is not yet supported.

All constraints and variables are listed below.

Overview

  • The linear OPF module can optimise the dispatch of generation and storage and the capacities of generation, storage and transmission infrastructure.

  • It is assumed that the load is inelastic and must be met in every snapshot (this will be relaxed in future versions).

  • The optimisation currently uses continuous variables for most functionality; unit commitment with binary variables is also implemented for generators.

  • The objective function is the total system cost for the snapshots optimised.

  • Each snapshot can be given a weighting \(w_t\) to represent e.g. multiple hours.

  • This set-up can also be used for stochastic optimisation, if you interpret the weighting as a probability.

  • Each transmission asset has a capital cost.

  • Each generation and storage asset has a capital cost and a marginal cost.

Execute:

network.lopf(snapshots, solver_name="glpk", solver_io=None,
             extra_functionality=None, solver_options={},
             keep_files=False, formulation="angles",
             extra_postprocessing=None)``

where snapshots is an iterable of snapshots, solver_name is a string, e.g. “gurobi” or “glpk”, solver_io is a string, extra_functionality is a function of network and snapshots that is called before the solver (see below), extra_postprocessing is a function of network, snapshots and duals that is called after solving (see below), solver_options is a dictionary of flags to pass to the solver, keep_files means that the .lp file is saved and formulation is a string in ["angles","cycles","kirchhoff","ptdf"] (see Passive branch flow formulations for more details).

Network.lopf(snapshots=None, pyomo=True, solver_name='glpk', solver_options={}, solver_logfile=None, formulation='kirchhoff', keep_files=False, extra_functionality=None, **kwargs)

Linear optimal power flow for a group of snapshots.

Parameters
  • snapshots (list or index slice) – A list of snapshots to optimise, must be a subset of network.snapshots, defaults to network.snapshots

  • pyomo (bool, default True) – Whether to use pyomo for building and solving the model, setting this to False saves a lot of memory and time.

  • solver_name (string) – Must be a solver name that pyomo recognises and that is installed, e.g. “glpk”, “gurobi”

  • solver_options (dictionary) – A dictionary with additional options that get passed to the solver. (e.g. {‘threads’:2} tells gurobi to use only 2 cpus)

  • solver_logfile (None|string) – If not None, sets the logfile option of the solver.

  • keep_files (bool, default False) – Keep the files that pyomo constructs from OPF problem construction, e.g. .lp file - useful for debugging

  • formulation (string) – Formulation of the linear power flow equations to use; must be one of [“angles”,”cycles”,”kirchhoff”,”ptdf”]

  • extra_functionality (callable function) – This function must take two arguments extra_functionality(network,snapshots) and is called after the model building is complete, but before it is sent to the solver. It allows the user to add/change constraints and add/change the objective function.

Other Parameters
  • ptdf_tolerance (float) – Only taking effect when pyomo is True. Value below which PTDF entries are ignored

  • free_memory (set, default {‘pyomo’}) – Only taking effect when pyomo is True. Any subset of {‘pypsa’, ‘pyomo’}. Allows to stash pypsa time-series data away while the solver runs (as a pickle to disk) and/or free pyomo data after the solution has been extracted.

  • solver_io (string, default None) – Only taking effect when pyomo is True. Solver Input-Output option, e.g. “python” to use “gurobipy” for solver_name=”gurobi”

  • skip_pre (bool, default False) – Only taking effect when pyomo is True. Skip the preliminary steps of computing topology, calculating dependent values and finding bus controls.

  • extra_postprocessing (callable function) – Only taking effect when pyomo is True. This function must take three arguments extra_postprocessing(network,snapshots,duals) and is called after the model has solved and the results are extracted. It allows the user to extract further information about the solution, such as additional shadow prices.

  • skip_objective (bool, default False) – Only taking effect when pyomo is False. Skip writing the default objective function. If False, a custom objective has to be defined via extra_functionality.

  • warmstart (bool or string, default False) – Only taking effect when pyomo is False. Use this to warmstart the optimization. Pass a string which gives the path to the basis file. If set to True, a path to a basis file must be given in network.basis_fn.

  • store_basis (bool, default True) – Only taking effect when pyomo is False. Whether to store the basis of the optimization results. If True, the path to the basis file is saved in network.basis_fn. Note that a basis can only be stored if simplex, dual-simplex, or barrier with crossover is used for solving.

  • keep_references (bool, default False) – Only taking effect when pyomo is False. Keep the references of variable and constraint names withing the network. These can be looked up in n.vars and n.cons after solving.

  • keep_shadowprices (bool or list of component names) – Only taking effect when pyomo is False. Keep shadow prices for all constraints, if set to True. If a list is passed the shadow prices will only be parsed for those constraint names. Defaults to [‘Bus’, ‘Line’, ‘GlobalConstraint’]. After solving, the shadow prices can be retrieved using pypsa.linopt.get_dual() with corresponding name

  • solver_dir (str, default None) – Only taking effect when pyomo is False. Path to directory where necessary files are written, default None leads to the default temporary directory used by tempfile.mkstemp().

Returns

  • status (str) – Status of optimization. Either “ok” if solution is optimal, or “warning” if not.

  • termination_condition (str) – More information on how the solver terminated. One of “optimal”, “suboptimal” (in which case a solution is still provided), “infeasible”, “infeasible or unbounded”, or “other”.

Important

Since version v0.16.0, PyPSA enables optimisation without the use of pyomo by setting pyomo=False. This make the lopf function much more efficient in terms of memory usage and time. For this purpose two new module were introduced, pypsa.linopf and pypsa.linopt wich mainly reflect the functionality of pypsa.opf and pypsa.opt but without using pyomo. Note that when setting pyomo to False, the extra_functionality has to be adapted to the appropriate syntax (see guidelines below). Some unit commitment functionality is not yet implemented without pyomo.

Warning

If the transmission capacity is changed in passive networks, then the impedance will also change (i.e. if parallel lines are installed). This is NOT reflected in the ordinary LOPF, however pypsa.linopf.ilopf covers this through an iterative process as done in here.

Optimising dispatch only - a market model

Capacity optimisation can be turned off so that only the dispatch is optimised, like a short-run electricity market model.

For simplified transmission representation using Net Transfer Capacities (NTCs), there is a Link component which does controllable power flow like a transport model (and can also represent a point-to-point HVDC link).

Optimising total annual system costs

To minimise long-run annual system costs for meeting an inelastic electrical load, capital costs for transmission and generation should be set to the annualised investment costs in e.g. EUR/MW/a, marginal costs for dispatch to e.g. EUR/MWh and the weightings (now with units hours per annum, h/a) are chosen such that

\[\sum_t w_t = 8760\]

In this case the objective function gives total system cost in EUR/a to meet the total load.

Stochastic optimisation

For the very simplest stochastic optimisation you can use the weightings w_t as probabilities for the snapshots, which can represent different load/weather conditions. More sophisticated functionality is planned.

Variables and notation summary

\(n \in N = \{0,\dots |N|-1\}\)

label the buses

\(t \in T = \{0,\dots |T|-1\}\)

label the snapshots

\(l \in L = \{0,\dots |L|-1\}\)

label the branches

\(s \in S = \{0,\dots |S|-1\}\)

label the different generator/storage types at each bus

\(w_t\)

weighting of time \(t\) in the objective function

\(g_{n,s,t}\)

dispatch of generator \(s\) at bus \(n\) at time \(t\)

\(\bar{g}_{n,s}\)

nominal power of generator \(s\) at bus \(n\)

\(\bar{g}_{n,s,t}\)

availability of generator \(s\) at bus \(n\) at time \(t\) per unit of nominal power

\(u_{n,s,t}\)

binary status variable for generator with unit commitment

\(suc_{n,s,t}\)

start-up cost if generator with unit commitment is started at time \(t\)

\(sdc_{n,s,t}\)

shut-down cost if generator with unit commitment is shut down at time \(t\)

\(c_{n,s}\)

capital cost of extending generator nominal power by one MW

\(o_{n,s}\)

marginal cost of dispatch generator for one MWh

\(f_{l,t}\)

flow of power in branch \(l\) at time \(t\)

\(F_{l}\)

capacity of branch \(l\)

\(\eta_{n,s}\)

efficiency of generator \(s\) at bus \(n\)

\(\eta_{l}\)

efficiency of controllable link \(l\)

\(e_s\)

CO2-equivalent-tonne-per-MWh of the fuel carrier \(s\)

Further definitions are given below.

Objective function

See pypsa.opf.define_linear_objective(network,snapshots).

The objective function is composed of capital costs \(c\) for each component and operation costs \(o\) for generators

\begin{gather*} \sum_{n,s} c_{n,s} \bar{g}_{n,s} + \sum_{n,s} c_{n,s} \bar{h}_{n,s} + \sum_{l} c_{l} F_l \\ + \sum_{t} w_t \left[\sum_{n,s} o_{n,s,t} g_{n,s,t} + \sum_{n,s} o_{n,s,t} h_{n,s,t} \right] + \sum_{t} \left[suc_{n,s,t} + sdc_{n,s,t} \right] \end{gather*}

Additional variables which do not appear in the objective function are the storage uptake variable, the state of charge and the voltage angle for each bus.

Generator constraints

These are defined in pypsa.opf.define_generator_variables_constraints(network,snapshots).

Generator nominal power and generator dispatch for each snapshot may be optimised.

Each generator has a dispatch variable \(g_{n,s,t}\) where \(n\) labels the bus, \(s\) labels the particular generator at the bus (e.g. it can represent wind/gas/coal generators at the same bus in an aggregated network) and \(t\) labels the time.

It obeys the constraints:

\[\tilde{g}_{n,s,t}*\bar{g}_{n,s} \leq g_{n,s,t} \leq \bar{g}_{n,s,t}*\bar{g}_{n,s}\]

where \(\bar{g}_{n,s}\) is the nominal power (generator.p_nom) and \(\tilde{g}_{n,s,t}\) and \(\bar{g}_{n,s,t}\) are time-dependent restrictions on the dispatch (per unit of nominal power) due to e.g. wind availability or power plant de-rating.

For generators with time-varying p_max_pu in network.generators_t the per unit availability \(\bar{g}_{n,s,t}\) is a time series.

For generators with static p_max_pu in network.generators the per unit availability is a constant.

If the generator’s nominal power \(\bar{g}_{n,s}\) is also the subject of optimisation (generator.p_nom_extendable == True) then limits generator.p_nom_min and generator.p_nom_max on the installable nominal power may also be introduced, e.g.

\[\tilde{g}_{n,s} \leq \bar{g}_{n,s} \leq \hat{g}_{n,s}\]

Generator unit commitment constraints

These are defined in pypsa.opf.define_generator_variables_constraints(network,snapshots).

Important

Unit commitment constraints will only be build fully if pyomo is set to True. If pyomo is set to False a simplified version of the unit commitment is calculated by ignoring the parameters min_up_time, min_down_time, start_up_cost, shut_down_cost, up_time_before and down_time_before.

The implementation is a complete implementation of the unit commitment constraints defined in Chapter 4.3 of Convex Optimization of Power Systems by Joshua Adam Taylor (CUP, 2015).

Unit commitment can be turned on for any generator by setting committable to be True. This introduces a times series of new binary status variables \(u_{n,s,t} \in \{0,1\}\), saved in network.generators_t.status, which indicates whether the generator is running (1) or not (0) in period \(t\). The restrictions on generator output now become:

\[u_{n,s,t}*\tilde{g}_{n,s,t}*\bar{g}_{n,s} \leq g_{n,s,t} \leq u_{n,s,t}*\bar{g}_{n,s,t}*\bar{g}_{n,s} \hspace{.5cm} \forall\, n,s,t\]

so that if \(u_{n,s,t} = 0\) then also \(g_{n,s,t} = 0\).

Note that a generator cannot be both extendable (generator.p_nom_extendable == True) and committable (generator.committable == True) because of the coupling of the variables \(u_{n,s,t}\) and \(\bar{g}_{n,s}\) here.

If the minimum up time \(T_{\textrm{min_up}}\) (generator.min_up_time) is set then we have for generic times

\[\sum_{t'=t}^{t+T_\textrm{min_up}} u_{n,s,t'}\geq T_\textrm{min_up} (u_{n,s,t} - u_{n,s,t-1}) \hspace{.5cm} \forall\, n,s,t\]

i.e. if the generator has just started up at time \(t\) then \(u_{n,s,t-1} = 0\), \(u_{n,s,t} = 1\) and \(u_{n,s,t} - u_{n,s,t-1} = 1\), so that it has to run for at least \(T_{\textrm{min_up}}\) periods.

The generator may have been up for some periods before the snapshots simulation period. If the up-time before snapshots starts is less than the minimum up-time, then the generator is forced to be up for the difference at the start of snapshots. If the start of snapshots is the start of network.snapshots, then the up-time before the simulation is read from the input variable generator.up_time_before. If snapshots falls in the middle of network.snapshots, then PyPSA assumes the statuses for hours before snapshots have been set by previous simulations, and reads back the previous up-time by examining the previous statuses. If the start of snapshots is very close to the start of network.snapshots, it will also take account of generator.up_time_before as well as the statuses in between.

At the end of snapshots the minimum up-time in the constraint is only enforced for the remaining snapshots, if the number of remaining snapshots is less than \(T_{\textrm{min_up}}\).

Similarly if the minimum down time \(T_{\textrm{min_down}}\) (generator.min_up_time) is set then we have

\[\sum_{t'=t}^{t+T_\textrm{min_down}} (1-u_{n,s,t'})\geq T_\textrm{min_down} (u_{n,s,t-1} - u_{n,s,t}) \hspace{.5cm} \forall\, n,s,t\]

You can also defined generator.down_time_before for periods before network.snapshots, analagous to the up time.

For non-zero start up costs \(suc_{n,s}\) a new variable \(suc_{n,s,t} \geq 0\) is introduced for each time period \(t\) and added to the objective function. The variable satisfies

\[suc_{n,s,t} \geq suc_{n,s} (u_{n,s,t} - u_{n,s,t-1}) \hspace{.5cm} \forall\, n,s,t\]

so that it is only non-zero if \(u_{n,s,t} - u_{n,s,t-1} = 1\), i.e. the generator has just started, in which case the inequality is saturated \(suc_{n,s,t} = suc_{n,s}\). Similarly for the shut down costs \(sdc_{n,s,t} \geq 0\) we have

\[sdc_{n,s,t} \geq sdc_{n,s} (u_{n,s,t-1} - u_{n,s,t}) \hspace{.5cm} \forall\, n,s,t\]

Generator ramping constraints

These are defined in pypsa.opf.define_generator_variables_constraints(network,snapshots).

The implementation follows Chapter 4.3 of Convex Optimization of Power Systems by Joshua Adam Taylor (CUP, 2015).

Ramp rate limits can be defined for increasing power output \(ru_{n,s}\) and decreasing power output \(rd_{n,s}\). By default these are null and ignored. They should be given per unit of the generator nominal power. The generator dispatch then obeys

\[-rd_{n,s} * \bar{g}_{n,s} \leq (g_{n,s,t} - g_{n,s,t-1}) \leq ru_{n,s} * \bar{g}_{n,s}\]

for \(t \in \{1,\dots |T|-1\}\).

For generators with unit commitment you can also specify ramp limits at start-up \(rusu_{n,s}\) and shut-down \(rdsd_{n,s}\)

\begin{gather*} \left[ -rd_{n,s}*u_{n,s,t} -rdsd_{n,s}(u_{n,s,t-1} - u_{n,s,t})\right] \bar{g}_{n,s} \\ \leq (g_{n,s,t} - g_{n,s,t-1}) \leq \\ \left[ru_{n,s}*u_{n,s,t-1} + rusu_{n,s} (u_{n,s,t} - u_{n,s,t-1})\right]\bar{g}_{n,s} \end{gather*}

Storage Unit constraints

These are defined in pypsa.opf.define_storage_variables_constraints(network,snapshots).

Storage nominal power and dispatch for each snapshot may be optimised.

With a storage unit the maximum state of charge may not be independently optimised from the maximum power output (they’re linked by the maximum hours variable) and the maximum power output is linked to the maximum power input. To optimise these capacities independently, build a storage unit out of the more fundamental Store and Link components.

The storage nominal power is given by \(\bar{h}_{n,s}\).

In contrast to the generator, which has one time-dependent variable, each storage unit has three:

The storage dispatch \(h_{n,s,t}\) (when it depletes the state of charge):

\[0 \leq h_{n,s,t} \leq \bar{h}_{n,s}\]

The storage uptake \(f_{n,s,t}\) (when it increases the state of charge):

\[0 \leq f_{n,s,t} \leq \bar{h}_{n,s}\]

and the state of charge itself:

\[0\leq soc_{n,s,t} \leq r_{n,s} \bar{h}_{n,s}\]

where \(r_{n,s}\) is the number of hours at nominal power that fill the state of charge.

The variables are related by

\[soc_{n,s,t} = \eta_{\textrm{stand};n,s}^{w_t} soc_{n,s,t-1} + \eta_{\textrm{store};n,s} w_t f_{n,s,t} - \eta^{-1}_{\textrm{dispatch};n,s} w_t h_{n,s,t} + w_t\textrm{inflow}_{n,s,t} - w_t\textrm{spillage}_{n,s,t}\]

\(\eta_{\textrm{stand};n,s}\) is the standing losses dues to e.g. thermal losses for thermal storage. \(\eta_{\textrm{store};n,s}\) and \(\eta_{\textrm{dispatch};n,s}\) are the efficiency losses for power going into and out of the storage unit.

There are two options for specifying the initial state of charge \(soc_{n,s,t=-1}\): you can set storage_unit.cyclic_state_of_charge = False (the default) and the value of storage_unit.state_of_charge_initial in MWh; or you can set storage_unit.cyclic_state_of_charge = True and then the optimisation assumes \(soc_{n,s,t=-1} = soc_{n,s,t=|T|-1}\).

If in the time series storage_unit_t.state_of_charge_set there are values which are not NaNs, then it will be assumed that these are fixed state of charges desired for that time \(t\) and these will be added as extra constraints. (A possible usage case would be a storage unit where the state of charge must empty every day.)

Store constraints

These are defined in pypsa.opf.define_store_variables_constraints(network,snapshots).

Store nominal energy and dispatch for each snapshot may be optimised.

The store nominal energy is given by \(\bar{e}_{n,s}\).

The store has two time-dependent variables:

The store dispatch \(h_{n,s,t}\):

\[-\infty \leq h_{n,s,t} \leq +\infty\]

and the energy:

\[\tilde{e}_{n,s} \leq e_{n,s,t} \leq \bar{e}_{n,s}\]

The variables are related by

\[e_{n,s,t} = \eta_{\textrm{stand};n,s}^{w_t} e_{n,s,t-1} - w_t h_{n,s,t}\]

\(\eta_{\textrm{stand};n,s}\) is the standing losses dues to e.g. thermal losses for thermal storage.

There are two options for specifying the initial energy \(e_{n,s,t=-1}\): you can set store.e_cyclic = False (the default) and the value of store.e_initial in MWh; or you can set store.e_cyclic = True and then the optimisation assumes \(e_{n,s,t=-1} = e_{n,s,t=|T|-1}\).

Passive branch flows: lines and transformers

See pypsa.opf.define_passive_branch_flows(network,snapshots) and pypsa.opf.define_passive_branch_constraints(network,snapshots) and pypsa.opf.define_branch_extension_variables(network,snapshots).

For lines and transformers, whose power flows according to impedances, the power flow \(f_{l,t}\) in AC networks is given by the difference in voltage angles \(\theta_{n,t}\) at bus0 and \(\theta_{m,t}\) at bus1 divided by the series reactance \(x_l\)

\[f_{l,t} = \frac{\theta_{n,t} - \theta_{m,t}}{x_l}\]

(For DC networks, replace the voltage angles by the difference in voltage magnitude \(\delta V_{n,t}\) and the series reactance by the series resistance \(r_l\).)

This flow is the limited by the capacity :math:F_l of the line

\[|f_{l,t}| \leq F_l\]

Note

If \(F_l\) is also subject to optimisation (branch.s_nom_extendable -- True), then the impedance \(x\) of the line is NOT automatically changed with the capacity (to represent e.g. parallel lines being added).

There are two choices here:

  1. Iterate the LOPF again with the updated impedances, see e.g. http://www.sciencedirect.com/science/article/pii/S0360544214000322#, like done by pypsa.linopf.ilopf

  2. João Gorenstein Dedecca has also implemented a MILP version of the transmission expansion, see https://github.com/jdedecca/MILP_PyPSA, which properly takes account of the impedance with a disjunctive relaxation. This will be pulled into the main PyPSA code base soon.

Passive branch flow formulations

PyPSA implements four formulations of the linear power flow equations that are mathematically equivalent, but may have different solving times. These different formulations are described and benchmarked in the arXiv preprint paper Linear Optimal Power Flow Using Cycle Flows.

You can choose the formulation by passing network.lopf the argument formulation, which must be in ["angles","cycles","kirchhoff","ptdf"].

  • angles is the standard formulations based on voltage angles described above, used for the linear power flow and found in textbooks.

  • ptdf uses the Power Transfer Distribution Factor (PTDF) formulation, found for example in http://www.sciencedirect.com/science/article/pii/S0360544214000322#.

  • kirchhoff and cycles are two new formulations based on a graph-theoretic decomposition of the network flows into a spanning tree and closed cycles.

Based on the benchmarking in Linear Optimal Power Flow Using Cycle Flows for standard networks, kirchhoff almost always solves fastest, averaging 3 times faster than the angles formulation and up to 20 times faster in specific cases. The speedup is higher for larger networks with dispatchable generators at most nodes.

Nodal power balances

See pypsa.opf.define_nodal_balances(network,snapshots).

This is the most important equation, which guarantees that the power balances at each bus \(n\) for each time \(t\).

\[\sum_{s} g_{n,s,t} + \sum_{s} h_{n,s,t} - \sum_{s} f_{n,s,t} - \sum_{l} K_{nl} f_{l,t} = \sum_{s} d_{n,s,t} \hspace{.4cm} \leftrightarrow \hspace{.4cm} w_t\lambda_{n,t}\]

Where \(d_{n,s,t}\) is the exogenous load at each node (load.p_set) and the incidence matrix \(K_{nl}\) for the graph takes values in \(\{-1,0,1\}\) depending on whether the branch \(l\) ends or starts at the bus. \(\lambda_{n,t}\) is the shadow price of the constraint, i.e. the locational marginal price, stored in network.buses_t.marginal_price.

The bus’s role is to enforce energy conservation for all elements feeding in and out of it (i.e. like Kirchhoff’s Current Law).

_images/buses.png

Global constraints

See pypsa.opf.define_global_constraints(network,snapshots).

Global constraints apply to more than one component.

Currently only “primary energy” constraints are defined. They depend on the power plant efficiency and carrier-specific attributes such as specific CO2 emissions.

Suppose there is a global constraint defined for CO2 emissions with sense <= and constant \textrm{CAP}_{CO2}. Emissions can come from generators whose energy carriers have CO2 emissions and from stores and storage units whose storage medium releases or absorbs CO2 when it is converted. Only stores and storage units with non-cyclic state of charge that is different at the start and end of the simulation can contribute.

If the specific emissions of energy carrier \(s\) is \(e_s\) (carrier.co2_emissions) CO2-equivalent-tonne-per-MWh and the generator with carrier \(s\) at node \(n\) has efficiency \(\eta_{n,s}\) then the CO2 constraint is

\[\sum_{n,s,t} \frac{1}{\eta_{n,s}} w_t\cdot g_{n,s,t}\cdot e_{n,s} + \sum_{n,s}\left(e_{n,s,t=-1} - e_{n,s,t=|T|-1}\right) \cdot e_{n,s} \leq \textrm{CAP}_{CO2} \hspace{.4cm} \leftrightarrow \hspace{.4cm} \mu\]

The first sum is over generators; the second sum is over stores and storage units. \(\mu\) is the shadow price of the constraint, i.e. the CO2 price in this case. \(\mu\) is an output of the optimisation stored in network.global_constraints.mu.

Custom constraints and other functionality

Since PyPSA v0.16.0, the lopf function is provided by two different modules. The ordinary implementation based on the pypsa.opf module uses pyomo to set up the linear optimisation problem and passing it to the solver. The implementation without pyomo, based on the module pypsa.linopf, uses a straight-forward approach to write out the .lp file directly and explicitly running it from a solver’s interface. Therefore the application of custom constraints depends on whether pyomo is activated or not.

In general for a custom constraint, pass the function network.lopf a function extra_functionality as an argument. This function must take two arguments extra_functionality(network,snapshots) and is called after the model building is complete, but before it is sent to the solver. It allows the user to add, change or remove constraints and alter the objective function.

1. pyomo is set to True

You can easily extend the optimisation problem constructed by PyPSA using the usual pyomo syntax.

The CHP example and the example that replaces generators and storage units with fundamental links and stores both pass an extra_functionality argument to the LOPF to add functionality.

The function extra_postprocessing is called after the model has solved and the results are extracted. This function must take three arguments extra_postprocessing(network,snapshots,duals). It allows the user to extract further information about the solution, such as additional shadow prices for constraints.

2. pyomo is set to False

In general when pyomo is disabled, all variable and constraint references are stored in the network object itself. Thus every variable and constraint is attached to a component, e.g. the dispatch variable of network.generators.p is attached to the component ‘Generator’ and can be easily accessed by

>>> get_var(n, 'Generator', 'p')

An additional constraint can easily be implemented by using the functions

  • pypsa.linopt.get_var for getting the variables which should be included in the constraint

  • pypsa.linopt.linexpr for creating linear expressions for the left hand side (lhs) of the constraint. Note that only the lhs includes all terms with variables, the rhs is a constant.

  • pypsa.linopt.define_constraints for defining a network constraint.

The are functions defined as such:

linopt.get_var(c, attr, pop=False)

Retrieves variable references for a given static or time-depending attribute of a given component. The function looks into n.variables to detect whether the variable is a time-dependent or static.

Parameters
  • n (pypsa.Network) –

  • c (str) – component name to which the constraint belongs

  • attr (str) – attribute name of the constraints

Example

>>> get_var(n, 'Generator', 'p')
linopt.linexpr(*, as_pandas=True, return_axes=False)

Elementwise concatenation of tuples in the form (coefficient, variables). Coefficient and variables can be arrays, series or frames. Per default returns a pandas.Series or pandas.DataFrame of strings. If return_axes is set to True the return value is split into values and axes, where values are the numpy.array and axes a tuple containing index and column if present.

Parameters
  • tuples (tuple of tuples) –

    Each tuple must of the form (coeff, var), where

    • coeff is a numerical value, or a numerical array, series, frame

    • var is a str or a array, series, frame of variable strings

  • as_pandas (bool, default True) – Whether to return to resulting array as a series, if 1-dimensional, or a frame, if 2-dimensional. Supersedes return_axes argument.

  • return_axes (Boolean, default False) – Whether to return index and column (if existent)

Example

Initialize coefficients and variables

>>> coeff1 = 1
>>> var1 = pd.Series(['a1', 'a2', 'a3'])
>>> coeff2 = pd.Series([-0.5, -0.3, -1])
>>> var2 = pd.Series(['b1', 'b2', 'b3'])

Create the linear expression strings

>>> linexpr((coeff1, var1), (coeff2, var2))
0    +1.0 a1 -0.5 b1
1    +1.0 a2 -0.3 b2
2    +1.0 a3 -1.0 b3
dtype: object

For a further step the resulting frame can be used as the lhs of pypsa.linopt.define_constraints()

For retrieving only the values:

>>> linexpr((coeff1, var1), (coeff2, var2), as_pandas=False)
array(['+1.0 a1 -0.5 b1', '+1.0 a2 -0.3 b2', '+1.0 a3 -1.0 b3'], dtype=object)
linopt.define_constraints(lhs, sense, rhs, name, attr='', axes=None, spec='')

Defines constraint(s) for pypsa-network with given left hand side (lhs), sense and right hand side (rhs). The constraints are stored in the network object under n.cons with key of the constraint name. If multiple constraints are defined at ones, only using np.arrays, then the axes argument can be used for defining the axes for the constraints (this is espececially recommended for time-dependent constraints). If one of lhs, sense and rhs is a pd.Series/pd.DataFrame the axes argument is not necessary.

Parameters
  • n (pypsa.Network) –

  • lhs (pd.Series/pd.DataFrame/np.array/str/float) – left hand side of the constraint(s), created with pypsa.linot.linexpr().

  • sense (pd.Series/pd.DataFrame/np.array/str/float) – sense(s) of the constraint(s)

  • rhs (pd.Series/pd.DataFrame/np.array/str/float) – right hand side of the constraint(s), must only contain pure constants, no variables

  • name (str) –

    general name of the constraint (or component which the constraint is referring to). The constraint will then be stored under:

    • n.cons[name].pnl if the constraint is two-dimensional

    • n.cons[name].df if the constraint is one-dimensional

  • attr (str default '') – Specifying name of the constraint, defines under which name the constraint(s) are stored in n.cons[name].pnl if two-dimensional or in n.cons[name].df if one-dimensional

  • axes (pd.Index or tuple of pd.Index objects, default None) – Specifies the axes if all of lhs, sense and rhs are np.arrays or single strings or floats.

Example

Let’s say we want to constraint all gas generators to a maximum of 100 MWh during the first 10 snapshots. We then firstly get all operational variables for this subset and constraint there sum to less equal 100.

>>> from pypsa.linopt import get_var, linexpr, defin_constraints
>>> gas_i = n.generators.query('carrier == "Natural Gas"').index
>>> gas_vars = get_var(n, 'Generator', 'p').loc[n.snapshots[:10], gas_i]
>>> lhs = linexpr((1, gas_vars)).sum().sum()
>>> define_(n, lhs, '<=', 100, 'Generator', 'gas_power_limit')

Now the constraint references can be accessed by pypsa.linopt.get_con() using

>>> cons = get_var(n, 'Generator', 'gas_power_limit')

Under the hook they are stored in n.cons.Generator.pnl.gas_power_limit. For retrieving their shadow prices add the general name of the constraint to the keep_shadowprices argument.

Note that this is usefull for the extra_functionality argument.

The function extra_postprocessing is not necessary when pyomo is deactivated. For retrieving additional shadow prices, just pass the name of the constraint, to which the constraint is attached, to the keep_shadowprices parameter of the lopf function.

Inputs

For the linear optimal power flow, the following data for each component are used. For almost all values, defaults are assumed if not explicitly set. For the defaults and units, see Components.

  • network.{snapshot_weightings}

  • bus.{v_nom, carrier}

  • load.{p_set}

  • generator.{p_nom, p_nom_extendable, p_nom_min, p_nom_max, p_min_pu, p_max_pu, marginal_cost, capital_cost, efficiency, carrier}

  • storage_unit.{p_nom, p_nom_extendable, p_nom_min, p_nom_max, p_min_pu, p_max_pu, marginal_cost, capital_cost, efficiency*, standing_loss, inflow, state_of_charge_set, max_hours, state_of_charge_initial, cyclic_state_of_charge}

  • store.{e_nom, e_nom_extendable, e_nom_min, e_nom_max, e_min_pu, e_max_pu, e_cyclic, e_initial, capital_cost, marginal_cost, standing_loss}

  • line.{x, s_nom, s_nom_extendable, s_nom_min, s_nom_max, capital_cost}

  • transformer.{x, s_nom, s_nom_extendable, s_nom_min, s_nom_max, capital_cost}

  • link.{p_min_pu, p_max_pu, p_nom, p_nom_extendable, p_nom_min, p_nom_max, capital_cost}

  • carrier.{carrier_attribute}

  • global_constraint.{type, carrier_attribute, sense, constant}

Note

Note that for lines and transformers you MUST make sure that \(x\) is non-zero, otherwise the bus admittance matrix will be singular.

Outputs

  • bus.{v_mag_pu, v_ang, p, marginal_price}

  • load.{p}

  • generator.{p, p_nom_opt}

  • storage_unit.{p, p_nom_opt, state_of_charge, spill}

  • store.{p, e_nom_opt, e}

  • line.{p0, p1, s_nom_opt, mu_lower, mu_upper}

  • transformer.{p0, p1, s_nom_opt, mu_lower, mu_upper}

  • link.{p0, p1, p_nom_opt, mu_lower, mu_upper}

  • global_constraint.{mu}

Utility Functions (without pyomo)

Build optimisation problems from PyPSA networks without Pyomo. Originally retrieved from nomopyomo ( -> ‘no more Pyomo’).

pypsa.linopf.assign_solution(n, sns, variables_sol, constraints_dual, keep_references=False, keep_shadowprices=None)

Helper function. Assigns the solution of a succesful optimization to the network.

pypsa.linopf.define_dispatch_for_extendable_and_committable_variables(n, sns, c, attr)

Initializes variables for power dispatch for a given component and a given attribute.

Parameters
  • n (pypsa.Network) –

  • c (str) – name of the network component

  • attr (str) – name of the attribute, e.g. ‘p’

pypsa.linopf.define_dispatch_for_extendable_constraints(n, sns, c, attr)

Sets power dispatch constraints for extendable devices for a given component and a given attribute.

Parameters
  • n (pypsa.Network) –

  • c (str) – name of the network component

  • attr (str) – name of the attribute, e.g. ‘p’

pypsa.linopf.define_dispatch_for_non_extendable_variables(n, sns, c, attr)

Initializes variables for power dispatch for a given component and a given attribute.

Parameters
  • n (pypsa.Network) –

  • c (str) – name of the network component

  • attr (str) – name of the attribute, e.g. ‘p’

pypsa.linopf.define_fixed_variable_constraints(n, sns, c, attr, pnl=True)

Sets constraints for fixing variables of a given component and attribute to the corresponding values in n.df(c)[attr + ‘_set’] if pnl is True, or n.pnl(c)[attr + ‘_set’]

Parameters
  • n (pypsa.Network) –

  • c (str) – name of the network component

  • attr (str) – name of the attribute, e.g. ‘p’

  • pnl (bool, default True) – Whether variable which should be fixed is time-dependent

pypsa.linopf.define_global_constraints(n, sns)

Defines global constraints for the optimization. Possible types are

  1. primary_energy

    Use this to constraint the byproducts of primary energy sources as CO2

  2. transmission_volume_expansion_limit

    Use this to set a limit for line volume expansion. Possible carriers are ‘AC’ and ‘DC’

  3. transmission_expansion_cost_limit

    Use this to set a limit for line expansion costs. Possible carriers are ‘AC’ and ‘DC’

pypsa.linopf.define_kirchhoff_constraints(n, sns)

Defines Kirchhoff voltage constraints

pypsa.linopf.define_nodal_balance_constraints(n, sns)

Defines nodal balance constraint.

pypsa.linopf.define_nominal_for_extendable_variables(n, c, attr)

Initializes variables for nominal capacities for a given component and a given attribute.

Parameters
  • n (pypsa.Network) –

  • c (str) – network component of which the nominal capacity should be defined

  • attr (str) – name of the variable, e.g. ‘p_nom’

pypsa.linopf.define_objective(n, sns)

Defines and writes out the objective function

pypsa.linopf.define_ramp_limit_constraints(n, sns)

Defines ramp limits for generators wiht valid ramplimit

pypsa.linopf.define_storage_unit_constraints(n, sns)

Defines state of charge (soc) constraints for storage units. In principal the constraints states:

previous_soc + p_store - p_dispatch + inflow - spill == soc

pypsa.linopf.define_store_constraints(n, sns)

Defines energy balance constraints for stores. In principal this states:

previous_e - p == e

pypsa.linopf.ilopf(n, snapshots=None, msq_threshold=0.05, min_iterations=1, max_iterations=100, track_iterations=False, **kwargs)

Iterative linear optimization updating the line parameters for passive AC and DC lines. This is helpful when line expansion is enabled. After each sucessful solving, line impedances and line resistance are recalculated based on the optimization result. If warmstart is possible, it uses the result from the previous iteration to fasten the optimization.

Parameters
  • snapshots (list or index slice) – A list of snapshots to optimise, must be a subset of network.snapshots, defaults to network.snapshots

  • msq_threshold (float, default 0.05) – Maximal mean square difference between optimized line capacity of the current and the previous iteration. As soon as this threshold is undercut, and the number of iterations is bigger than ‘min_iterations’ the iterative optimization stops

  • min_iterations (integer, default 1) – Minimal number of iteration to run regardless whether the msq_threshold is already undercut

  • max_iterations (integer, default 100) – Maximal numbder of iterations to run regardless whether msq_threshold is already undercut

  • track_iterations (bool, default False) – If True, the intermediate branch capacities and values of the objective function are recorded for each iteration. The values of iteration 0 represent the initial state.

  • **kwargs – Keyword arguments of the lopf function which runs at each iteration

pypsa.linopf.network_lopf(n, snapshots=None, solver_name='cbc', solver_logfile=None, extra_functionality=None, skip_objective=False, extra_postprocessing=None, formulation='kirchhoff', keep_references=False, keep_files=False, keep_shadowprices=['Bus', 'Line', 'Transformer', 'Link', 'GlobalConstraint'], solver_options=None, warmstart=False, store_basis=False, solver_dir=None)

Linear optimal power flow for a group of snapshots.

Parameters
  • snapshots (list or index slice) – A list of snapshots to optimise, must be a subset of network.snapshots, defaults to network.snapshots

  • solver_name (string) – Must be a solver name that pyomo recognises and that is installed, e.g. “glpk”, “gurobi”

  • pyomo (bool, default True) – Whether to use pyomo for building and solving the model, setting this to False saves a lot of memory and time.

  • solver_logfile (None|string) – If not None, sets the logfile option of the solver.

  • solver_options (dictionary) – A dictionary with additional options that get passed to the solver. (e.g. {‘threads’:2} tells gurobi to use only 2 cpus)

  • solver_dir (str, default None) – Path to directory where necessary files are written, default None leads to the default temporary directory used by tempfile.mkstemp().

  • keep_files (bool, default False) – Keep the files that pyomo constructs from OPF problem construction, e.g. .lp file - useful for debugging

  • formulation (string) – Formulation of the linear power flow equations to use; must be one of [“angles”,”cycles”,”kirchhoff”,”ptdf”]

  • extra_functionality (callable function) – This function must take two arguments extra_functionality(network,snapshots) and is called after the model building is complete, but before it is sent to the solver. It allows the user to add/change constraints and add/change the objective function.

  • skip_objective (bool, default False) – Skip writing the default objective function. If False, a custom objective has to be defined via extra_functionality.

  • extra_postprocessing (callable function) – This function must take three arguments extra_postprocessing(network,snapshots,duals) and is called after the model has solved and the results are extracted. It allows the user to extract further information about the solution, such as additional shadow prices.

  • warmstart (bool or string, default False) – Use this to warmstart the optimization. Pass a string which gives the path to the basis file. If set to True, a path to a basis file must be given in network.basis_fn.

  • store_basis (bool, default False) – Whether to store the basis of the optimization results. If True, the path to the basis file is saved in network.basis_fn. Note that a basis can only be stored if simplex, dual-simplex, or barrier with crossover is used for solving.

  • keep_references (bool, default False) – Keep the references of variable and constraint names withing the network. These can be looked up in n.vars and n.cons after solving.

  • keep_shadowprices (bool or list of component names) – Keep shadow prices for all constraints, if set to True. If a list is passed the shadow prices will only be parsed for those constraint names. Defaults to [‘Bus’, ‘Line’, ‘GlobalConstraint’]. After solving, the shadow prices can be retrieved using pypsa.linopt.get_dual() with corresponding name

pypsa.linopf.prepare_lopf(n, snapshots=None, keep_files=False, skip_objective=False, extra_functionality=None, solver_dir=None)

Sets up the linear problem and writes it out to a lp file.

Returns

  • Tuple (fdp, problem_fn) indicating the file descriptor and the file name of

  • the lp file

Utility Functions (with pyomo)

Optimal Power Flow functions.

pypsa.opf.define_nodal_balances(network, snapshots)

Construct the nodal balance for all elements except the passive branches.

Store the nodal balance expression in network._p_balance.

pypsa.opf.define_passive_branch_flows_with_kirchhoff(network, snapshots, skip_vars=False)

define passive branch flows with the kirchoff method

pypsa.opf.define_sub_network_cycle_constraints(subnetwork, snapshots, passive_branch_p, attribute)

Constructs cycle_constraints for a particular subnetwork

pypsa.opf.network_lopf(network, snapshots=None, solver_name='glpk', solver_io=None, skip_pre=False, extra_functionality=None, solver_logfile=None, solver_options={}, keep_files=False, formulation='angles', ptdf_tolerance=0.0, free_memory={}, extra_postprocessing=None)

Linear optimal power flow for a group of snapshots.

Parameters
  • snapshots (list or index slice) – A list of snapshots to optimise, must be a subset of network.snapshots, defaults to network.snapshots

  • solver_name (string) – Must be a solver name that pyomo recognises and that is installed, e.g. “glpk”, “gurobi”

  • solver_io (string, default None) – Solver Input-Output option, e.g. “python” to use “gurobipy” for solver_name=”gurobi”

  • skip_pre (bool, default False) – Skip the preliminary steps of computing topology, calculating dependent values and finding bus controls.

  • extra_functionality (callable function) – This function must take two arguments extra_functionality(network,snapshots) and is called after the model building is complete, but before it is sent to the solver. It allows the user to add/change constraints and add/change the objective function.

  • solver_logfile (None|string) – If not None, sets the logfile option of the solver.

  • solver_options (dictionary) – A dictionary with additional options that get passed to the solver. (e.g. {‘threads’:2} tells gurobi to use only 2 cpus)

  • keep_files (bool, default False) – Keep the files that pyomo constructs from OPF problem construction, e.g. .lp file - useful for debugging

  • formulation (string) – Formulation of the linear power flow equations to use; must be one of [“angles”,”cycles”,”kirchhoff”,”ptdf”]

  • ptdf_tolerance (float) – Value below which PTDF entries are ignored

  • free_memory (set, default {'pyomo'}) – Any subset of {‘pypsa’, ‘pyomo’}. Allows to stash pypsa time-series data away while the solver runs (as a pickle to disk) and/or free pyomo data after the solution has been extracted.

  • extra_postprocessing (callable function) – This function must take three arguments extra_postprocessing(network,snapshots,duals) and is called after the model has solved and the results are extracted. It allows the user to extract further information about the solution, such as additional shadow prices.

Returns

Return type

None

pypsa.opf.network_lopf_build_model(network, snapshots=None, skip_pre=False, formulation='angles', ptdf_tolerance=0.0)

Build pyomo model for linear optimal power flow for a group of snapshots.

Parameters
  • snapshots (list or index slice) – A list of snapshots to optimise, must be a subset of network.snapshots, defaults to network.snapshots

  • skip_pre (bool, default False) – Skip the preliminary steps of computing topology, calculating dependent values and finding bus controls.

  • formulation (string) – Formulation of the linear power flow equations to use; must be one of [“angles”,”cycles”,”kirchhoff”,”ptdf”]

  • ptdf_tolerance (float) – Value below which PTDF entries are ignored

Returns

Return type

network.model

pypsa.opf.network_lopf_prepare_solver(network, solver_name='glpk', solver_io=None)

Prepare solver for linear optimal power flow.

Parameters
  • solver_name (string) – Must be a solver name that pyomo recognises and that is installed, e.g. “glpk”, “gurobi”

  • solver_io (string, default None) – Solver Input-Output option, e.g. “python” to use “gurobipy” for solver_name=”gurobi”

Returns

Return type

None

pypsa.opf.network_lopf_solve(network, snapshots=None, formulation='angles', solver_options={}, solver_logfile=None, keep_files=False, free_memory={'pyomo'}, extra_postprocessing=None)

Solve linear optimal power flow for a group of snapshots and extract results.

Parameters
  • snapshots (list or index slice) – A list of snapshots to optimise, must be a subset of network.snapshots, defaults to network.snapshots

  • formulation (string) – Formulation of the linear power flow equations to use; must be one of [“angles”,”cycles”,”kirchhoff”,”ptdf”]; must match formulation used for building the model.

  • solver_options (dictionary) – A dictionary with additional options that get passed to the solver. (e.g. {‘threads’:2} tells gurobi to use only 2 cpus)

  • solver_logfile (None|string) – If not None, sets the logfile option of the solver.

  • keep_files (bool, default False) – Keep the files that pyomo constructs from OPF problem construction, e.g. .lp file - useful for debugging

  • free_memory (set, default {'pyomo'}) – Any subset of {‘pypsa’, ‘pyomo’}. Allows to stash pypsa time-series data away while the solver runs (as a pickle to disk) and/or free pyomo data after the solution has been extracted.

  • extra_postprocessing (callable function) – This function must take three arguments extra_postprocessing(network,snapshots,duals) and is called after the model has solved and the results are extracted. It allows the user to extract further information about the solution, such as additional shadow prices.

Returns

Return type

None

pypsa.opf.network_opf(network, snapshots=None)

Optimal power flow for snapshots.