# Optimal Power Flow¶

The operation of power systems relies on a number of optimization tasks known as Optimal Power Flow (OPF) problems. The objective of a standard OPF problem is to minimize operational costs such that the underlying grid constraints on generation, demand, and voltage limits are satisfied.

This Mod considers the Alternating Current (AC) and Direct Current (DC) OPF formulations. The ACOPF problem, in its natural form, requires the introduction of complex numbers to formulate voltage constraints. This Mod uses a cartesian-coordinates formulation of ACOPF that reformulates complex-valued terms via nonconvex quadratic relationships. The DCOPF problem approximates the ACOPF problem by making additional assumptions to produce linear constraints. While the additional assumptions result in a potential loss of solution accuracy, they make the DCOPF problem much easier to solve. This is especially useful if the solution accuracy can be neglected in favour of solution time and problem size. Please refer to the Full Problem Specification for full details of the formulations.

Here we assume basic familiarity with concepts such as *voltage* (potential
energy), *current* (charge flow), and *power* (instantaneous energy generation
or consumption). The engineering community also uses the terms *bus* to refer to
nodes in a network and branch to refer to arcs in a network (a connection
between two buses, typically a line or a transformer). Please refer to the
Recommended Literature section for more details and comprehensive
descriptions of power systems and the underlying problems.

## MATPOWER Case Input Format¶

This Mod has multiple API functions, each of which accepts a dictionary as input
describing an OPF case. This case dictionary follows the MATPOWER Case Format
conventions
and holds all essential information about the underlying network: buses, branch
connections, and generators. Several pre-defined MATPOWER cases can be loaded
directly using the `gurobi_optimods.datasets`

module. The following code loads
an example case with 9 buses and displays the key fields:

`baseMVA`

holds the base units of voltage for this case.

```
>>> from gurobi_optimods import datasets
>>>
>>> case = datasets.load_opf_example("case9")
>>> case['baseMVA']
100.0
```

Each entry in the list `case['bus']`

represents a bus in the network.

```
>>> from pprint import pprint
>>> from gurobi_optimods import datasets
>>>
>>> case = datasets.load_opf_example("case9")
>>> pprint(case['bus'][0]) # The first bus
{'Bs': 0.0,
'Gs': 0.0,
'Pd': 0.0,
'Qd': 0.0,
'Va': 0.0,
'Vm': 1.0,
'Vmax': 1.1,
'Vmin': 0.9,
'area': 1.0,
'baseKV': 345.0,
'bus_i': 1,
'type': 3,
'zone': 1.0}
```

Each entry in the list `case['branch']`

represents a branch in the
network. Note that the `fbus`

and `tbus`

fields must refer to a bus
by its `bus_i`

value.

```
>>> from pprint import pprint
>>> from gurobi_optimods import datasets
>>>
>>> case = datasets.load_opf_example("case9")
>>> pprint(case['branch'][0]) # The first branch
{'angle': 0.0,
'angmax': 360.0,
'angmin': -360.0,
'b': 0.0,
'fbus': 1,
'r': 0.0,
'rateA': 250.0,
'rateB': 250.0,
'rateC': 250.0,
'ratio': 0.0,
'status': 1.0,
'tbus': 4,
'x': 0.0576}
```

Each entry in the list `case['gen']`

represents a generator in the
network. Note that the `bus`

field must refer to a bus by its
`bus_i`

value.

```
>>> from pprint import pprint
>>> from gurobi_optimods import datasets
>>>
>>> case = datasets.load_opf_example("case9")
>>> pprint(case['gen'][0]) # The first generator
{'Pc1': 0,
'Pc2': 0,
'Pg': 0,
'Pmax': 250,
'Pmin': 10,
'Qc1max': 0,
'Qc1min': 0,
'Qc2max': 0,
'Qc2min': 0,
'Qg': 0,
'Qmax': 300,
'Qmin': -300,
'Vg': 1,
'apf': 0,
'bus': 1,
'mBase': 100,
'ramp_10': 0,
'ramp_30': 0,
'ramp_agc': 0,
'ramp_q': 0,
'status': 1}
```

Each entry in the list `case['gencost']`

corresponds to the generator
at the same position in `case['gen']`

, and defines the cost function
of that generator.

```
>>> from pprint import pprint
>>> from gurobi_optimods import datasets
>>>
>>> case = datasets.load_opf_example("case9")
>>> pprint(case['gencost'][0]) # Cost function for the first generator
{'costtype': 2.0,
'costvector': [0.11, 5.0, 150.0],
'n': 3.0,
'shutdown': 0.0,
'startup': 1500.0}
```

Warning

The Mod only supports generator costs with `costtype = 2`

and ```
n
<= 3
```

(i.e. linear and quadratic cost functions) in the `gencost`

structure, as these are the most commonly used settings in practice.
If a different costtype or larger value of `n`

value is provided,
an error will be issued.

Cases can also be read directly from MATPOWER format using
`read_case_matpower()`

. This function reads a standard
MATLAB `.mat`

data file holding the case data into the native Python format
accepted by the Mod. For example:

```
case = opf.read_case_matpower("my_case.mat")
```

## Solving an OPF Problem¶

After reading in or otherwise generating a MATPOWER case in the format described
above, we can solve an OPF problem defined by the given network data. For this
task, we use the `solve_opf()`

function. We can define
the type of the OPF problem that we want to solve by defining the `opftype`

argument when calling the function. Currently, the available options are `AC`

,
`ACrelax`

, and `DC`

.

The

`AC`

setting solves an ACOPF problem defined by the given network data. The ACOPF problem is formulated as a nonconvex bilinear model as described in the ACOPF section of the Full Problem Specification. This setting yields the most accurate model of the physical power system. However, it is also the most difficult problem to solve and thus usually leads to the longer runtimes.The

`ACrelax`

setting solves a Second Order Cone (SOC) relaxation of the nonconvex bilinear ACOPF problem formulation defined by the given network data. The relaxation is constructed by dropping nonconvex bilinear terms but simultaneously keeping the convex JABR inequalities, see JABR Relaxation for more details. This setting often yields a good approximation of the physical power system and is of moderate difficulty.The

`DC`

setting solves a DCOPF problem defined by the given network data. The DCOPF problem is a linear approximation of the ACOPF problem, see DCOPF section of the Full Problem Specification for more details. This setting only yields a crude approximation of the physical power system, but is usually an easy problem that can be solved very quickly even for large networks.

The `solve_opf`

function solves an `AC`

problem unless the `opftype`

argument specifies otherwise.

```
from gurobi_optimods import opf
from gurobi_optimods import datasets
case = datasets.load_opf_example("case9")
result = opf.solve_opf(case, opftype="AC")
```

```
...
Optimize a model with 73 rows, 107 columns and 208 nonzeros
...
Optimal solution found...
...
Objective value = 5296...
...
```

ACOPF and Branch Switching models are most often very hard to solve to
optimality. For this reason, it is recommended to specify a solver time limit
using the `time_limit`

parameter. If the problem has not been solved to
optimality within the time limit, the best known solution will be returned.

```
result = opf.solve_opf(case, opftype="AC", time_limit=60)
```

## Solution Format¶

The Mod returns the result as a dictionary, following the same MATPOWER Case Format conventions as the case dictionary. However, in the result dictionary some object entries are modified compared to the input case dictionary. These modified fields hold the solution values of the optimization. There are also additional fields to store the solution information, as specified below.

The `bus`

entries `result['bus'][i]['Vm']`

and
`result['bus'][i]['Va']`

store the voltage magnitude (Vm) and voltage
angle (Va) values in the OPF solution for bus i.

If a DCOPF problem was solved, additional fields
`result['bus'][i]['mu']`

hold the shadow prices for balance
constraints at bus i.

```
>>> result['bus'][0]
{... 'Vm': 1.09..., 'Va': 0, ...}
```

The additional `branch`

entries populated in the solution data are:

`Pf`

: real power injection at the “from” end of the branch,`Pt`

: real power injection at the “to” end of the branch,`Qf`

: reactive power injection at the “from” end of the branch, and`Qt`

: reactive power injection at the “to” end of the branch.

The `switching`

field defines whether a branch is turned on (1) or off
(0) in the given result.

```
>>> result['branch'][1]
{... 'Pf': 35.2..., 'Pt': -35.0..., 'Qf': -3.8..., 'Qt': -13.8..., ...}
```

Then `gen`

entries `result['gen'][i]['Pg']`

and
`result['gen'][i]['Qg']`

hold the real and reactive power injection
values, respectively, for the \(i^{th}\) generator.

```
>>> result['gen'][2]
{... 'Pg': 94.1..., 'Qg': -22.6..., ...}
```

## Plotting Feasible Solutions¶

In addition to solving an OPF problem, this Mod also provides plotting functions to display graphical representation of the network and the OPF result. There are already very involved graphical tools to represent OPF solutions provided by other packages such as:

thus the graphical representation provided by this Mod is kept intentionally
simple. In order to use this functionality, it is necessary to install the
`plotly`

package as follows:

```
pip install plotly
```

In order to plot a previously obtained result, you must provide \((x, y)\) coordinates for all buses in the network. Coordinates are provided as a dictionary mapping bus IDs to coordinates. The OptiMods datasets module provides an example set of coordinates for plotting the 9 bus test case:

```
>>> coordinates = datasets.load_opf_extra("case9-coordinates")
>>> from pprint import pprint
>>> pprint(coordinates)
{1: (44.492, -73.208),
2: (41.271, -73.953),
3: (41.574, -73.966),
4: (40.814, -72.94),
5: (43.495, -76.451),
6: (42.779, -78.427),
7: (44.713, -73.456),
8: (43.066, -76.214),
9: (43.048, -78.854)}
```

Given a solution and a coordinate mapping, the plotting functions return plotly figures which can be displayed in a web browser. In the following code, we solve the DCOPF problem for a real-world dataset for the city of New York and produce a solution plot:

```
case = datasets.load_opf_example("caseNY")
solution = opf.solve_opf(case, opftype='DC')
coordinates = datasets.load_opf_extra("caseNY-coordinates")
fig = opf.solution_plot(case, coordinates, solution)
fig.show() # open plot in a browser window
```

The above image shows the grid solution generated from the given network data, plotted using the provided coordinates. The colored circles depict generators and the amount of power they generate:

Black bus: Power generation \(\leq 75\) and load \(< 50\)

Blue bus: Power generation \(\leq 75\) and load \(\geq 50\)

Purple bus: Power generation \(> 75\)

Orange bus: Power generation \(> 150\)

Red bus: Power generation \(> 500\)

## Branch Switching¶

An important extension of the OPF problem is Branch Switching, where branches may be turned off. Note that already turning off a single branch changes the whole power flow through the network. Thus in practice, it is rare that branches are turned off at all, even if this option is enabled. If any are turned off, then it is usually only a small fraction of the overall power grid. For the mathematical formulation, please refer to the Branch Switching subsection of the Full Problem Specification.

To enable branch switching in a given OPF problem, set the `branch_switching`

argument to `True`

when calling `gurobi_optimods.opf.solve_opf()`

. The Mod
additionally offers the possibility to control the number of branches that must
remain switched on via the `min_active_branches`

argument. In practice, it is
expected that only a very small fraction of branches are turned off. Thus, the
default value of the `min_active_branches`

argument is 0.9 (90%). In the
following example, we solve a modified version of the 9 bus network to see
whether branch switching allows a better solution.

```
case = datasets.load_opf_example("case9-switching")
result = opf.solve_opf(
case, opftype="AC", branch_switching=True, min_active_branches=0.1
)
```

Plotting the resulting solution shows that one branch has been turned off in the
optimal solution. Please note, that the used examplary network has been
artificially adjusted to achieve this result and this is **not** the usual
behavior in a realistic power grid of such small size.

## Violations for Pre-defined Voltage Values¶

In practice we may have voltage magnitudes and voltage angles for each bus at
hand and would like to know whether these values are actually feasible within a
given network. To tackle this problem, we can use the
`gurobi_optimods.opf.compute_violations()`

function. This function takes a
set of bus voltages in addition to the case data and returns the computed
violations in this voltage solution.

Bus voltages are provided as a dictionary mapping the Bus ID to a pair \((V_m, V_a)\) where \(V_m\) is the voltage magnitude and \(V_a\) is the voltage angle. An example is provided for the 9 bus case:

```
>>> voltages = datasets.load_opf_extra("case9-voltages")
>>> pprint(voltages)
{1: (1.089026, 0.0),
2: (1.099999, 20.552543),
3: (1.090717, 16.594399),
4: (1.084884, -2.408447),
5: (1.096711, 2.43001),
6: (1.099999, 11.859651),
7: (1.072964, 9.257936),
8: (1.066651, 11.200108),
9: (1.08914, 2.847507)}
```

Using this voltage data, we can check for possible model violations by calling
the `gurobi_optimods.opf.compute_violations()`

function. The function
returns a dictionary which follows the MATPOWER Case Format with
additional fields storing the violations for particular buses and branches.

The following fields in the violations dictionary are added to store violations data:

`violation['bus'][i]['Vmviol']`

Voltage magnitude violation at bus i`violation['bus'][i]['Pviol']`

real power injection violation at bus i`violation['bus'][i]['Qviol']`

reactive power injection violation at bus i`violation['branch'][i]['limitviol']`

branch limit violation at branch i

```
volts_dict = datasets.load_opf_extra("case9-voltages")
case = datasets.load_opf_example("case9")
violations = opf.compute_violations(case, volts_dict)
```

```
>>> print(violations['branch'][6]['limitviol'])
66.33435016796234
>>> print(violations['bus'][3]['Pviol'])
-318.8997836192236
```

In this case, the limit at branch 6 and the real power injection at bus 3 are violated by the given input voltages.

## Plotting Violations¶

Similar to generating a graphical representation of a feasible solution, it is
also possible to generate a figure representing the violations within a given
power grid. We can use the `gurobi_optimods.opf.violation_plot()`

for this.
The violations result is provided as input, along with the case and coordinate
data, to produce this plot.

```
voltages = datasets.load_opf_extra("case9-voltages")
case = datasets.load_opf_example("case9")
coordinates = datasets.load_opf_extra("case9-coordinates")
violations = opf.compute_violations(case, voltages)
fig = opf.violation_plot(case, coordinates, violations)
```

The above image shows the power grid generated from the given network data together with coordinate and violation data. The red circles depict buses where the voltage magnitude or real or reactive power injections are violated. Red marked branches depict branches with violated limits.

## Recommended Literature¶

Power systems and the optimal power flow problem are well studied. For a more comprehensive descrition, we recommend the following literature.

G. Andersson.

*Modelling and Analysis of Electric Power Systems*. Power Systems Laboratory, ETH Zürich, 2004.A.R. Bergen and V. Vittal.

*Power Systems Analysis*. Prentice-Hall, 1999.D. Bienstock.

*Electrical Transmission Systems Cascades and Vulnerability, an Operations Research viewpoint*. SIAM, 2015. ISBN 978-1-61197-415-7.D.K. Molzahn and I.A. Hiskens.

*A survey of relaxations and approximations of the power flow equations*. Foundations and Trends in Electric Energy Systems, 4:1–221, 2019.J.D. Glover, M.S. Sarma, and T.J. Overbye.

*Power System Analysis and Design*. CENGAGE Learning, 2012.

## Formulation Details¶

Full details of the formulations used in the OPF Mod are given in the following sections: