Optimal Power Flow¶
See the module pypsa.opf
.
NonLinear Optimal Power Flow¶
Optimisation with the full nonlinear power flow equations is not yet supported.
Linear Optimal Power Flow¶
Optimisation with the linearised power flow equations for (mixed) AC and DC networks is fully supported.
All constraints and variables are listed below.
Overview¶
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, 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 InputOutput 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 (Nonestring) – 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 timeseries 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
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 setup 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.
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 LOPF, so the network equations may no longer be valid. Note also that all the expansion is continuous.
Optimising dispatch only: a market model¶
Capacity optimisation can be turned off so that only the dispatch is optimised, like a shortrun 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 pointtopoint HVDC link).
Optimising total annual system costs¶
To minimise longrun 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
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 N1\}\) label the buses
\(t \in T = \{0,\dots T1\}\) label the snapshots
\(l \in L = \{0,\dots L1\}\) label the branches
\(s \in S = \{0,\dots S1\}\) 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}\) startup cost if generator with unit commitment is started at time \(t\)
\(sdc_{n,s,t}\) shutdown 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\) CO2equivalenttonneperMWh 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
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:
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
timedependent restrictions on the dispatch (per unit of nominal
power) due to e.g. wind availability or power plant derating.
For generators with timevarying 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.
Generator unit commitment constraints¶
These are defined in pypsa.opf.define_generator_variables_constraints(network,snapshots)
.
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:
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
i.e. if the generator has just started up at time \(t\) then \(u_{n,s,t1} = 0\), \(u_{n,s,t} = 1\) and \(u_{n,s,t}  u_{n,s,t1} = 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 uptime before snapshots
starts is less than the minimum uptime, 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 uptime 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 uptime 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 uptime 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
You can also defined generator.down_time_before
for periods before network.snapshots
, analagous to the up time.
For nonzero 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
so that it is only nonzero if \(u_{n,s,t}  u_{n,s,t1} = 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
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
for \(t \in \{1,\dots T1\}\).
For generators with unit commitment you can also specify ramp limits at startup \(rusu_{n,s}\) and shutdown \(rdsd_{n,s}\)
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 timedependent variable, each storage unit has three:
The storage dispatch \(h_{n,s,t}\) (when it depletes the state of charge):
The storage uptake \(f_{n,s,t}\) (when it increases the state of charge):
and the state of charge itself:
where \(r_{n,s}\) is the number of hours at nominal power that fill the state of charge.
The variables are related by
\(\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=T1}\).
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 timedependent variables:
The store dispatch \(h_{n,s,t}\):
and the energy:
The variables are related by
\(\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=T1}\).
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\)
(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
Note that 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:
Iterate the LOPF again with the updated impedances (see e.g. http://www.sciencedirect.com/science/article/pii/S0360544214000322#).
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
andcycles
are two new formulations based on a graphtheoretic 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.
Controllable branch flows: links¶
See pypsa.opf.define_controllable_branch_flows(network,snapshots)
and pypsa.opf.define_branch_extension_variables(network,snapshots)
.
For links, whose power flow is controllable, there is simply an optimisation variable for each component which satisfies
If the link flow is positive \(f_{l,t} > 0\) then it withdraws
\(f_{l,t}\) from bus0
and feeds in \(\eta_l f_{l,t}\) to
bus1
, where \(\eta_l\) is the link efficiency.
If additional output buses busi
for \(i=2,3,\dots\) are
defined (i.e. bus2
, bus3
, etc) and their associated
efficiencies efficiencyi
, i.e. \(\eta_{i,l}\), then at
busi
the feedin is \(\eta_{i,l} f_{l,t}\). See also
Link with multiple outputs or inputs.
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\).
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).
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 carrierspecific 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 noncyclic
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
) CO2equivalenttonneperMWh and the
generator with carrier \(s\) at node \(n\) has efficiency
\(\eta_{n,s}\) then the CO2 constraint is
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¶
PyPSA uses the Python optimisation language pyomo to construct the OPF problem. You can easily
extend the optimisation problem constructed by PyPSA using the usual
pyomo syntax. To do this, 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.
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.
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 nonzero, 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}