Optimal Power Flow¶
Alternating current optimal power flow (ACOPF) is a fundamental nonlinear optimization problem in power systems. It is used to determine the optimal power flow of a power system to minimize the generation cost while satisfying the power flow equations and system constraints.
In this example, we will show how to use PyOptInterface to solve an single-period optimal power flow problem using the structured nonlinear programming formulation.
Problem Formulation¶
\[\begin{split}
\min \quad & \sum_{i \in G} a_i P_{i}^2 + b_i P_{i} + c_i \\
\textrm{s.t.} \quad & \theta_r = 0 \quad & \forall r \in R \\
\quad & P_{min,i} \leq P_{i} \leq P_{max,i} \quad & \forall i \in G \\
\quad & Q_{min,i} \leq Q_{i} \leq Q_{max,i} \quad & \forall i \in G \\
\quad & V_{min,b} \leq V_{b} \leq V_{max,b} \quad & \forall b \in BUS \\
\quad & \sum_{i \in G_b} P_{i} - \sum_{i \in D_b} L^P_{i} - G_{sh,b} V_b^2 = \sum_{(i, j) \in BRANCH} P_{ij} \quad & \forall b \in BUS \\
\quad & \sum_{i \in G_b} Q_{i} - \sum_{i \in D_b} L^Q_{i} + B_{sh,b} V_b^2 = \sum_{(i, j) \in BRANCH} Q_{ij} \quad & \forall b \in BUS \\
\quad & P_{ij} = G_{ij} V_i^2 - V_i V_j (G_{ij}\cos(\theta_i-\theta_j)+B_{ij}\sin(\theta_i-\theta_j)) \quad & \forall (i, j) \in BRANCH \\
\quad & Q_{ij} = -B^C_{ij} V_i^2 - B_{ij} V_i^2 - V_i V_j (G_{ij}\sin(\theta_i-\theta_j)-B_{ij}\cos(\theta_i-\theta_j)) \quad & \forall (i, j) \in BRANCH \\
\quad & P_{ji} = G_{ij} V_j^2 - V_i V_j (G_{ij}\cos(\theta_j-\theta_i)+B_{ij}\sin(\theta_j-\theta_i)) \quad & \forall (i, j) \in BRANCH \\
\quad & Q_{ji} = -B^C_{ij} V_j^2 - B_{ij} V_j^2 - V_i V_j (G_{ij}\sin(\theta_j-\theta_i)-B_{ij}\cos(\theta_j-\theta_i)) \quad & \forall (i, j) \in BRANCH \\
\quad & -S_{max,ij} \leq P_{ij}^2 + Q_{ij}^2 \leq S_{max,ij} \quad & \forall (i, j) \in BRANCH \\
\quad & -\Delta\theta_{min,ij} \leq \theta_i - \theta_j \leq -\Delta\theta_{max,ij} \quad & \forall (i, j) \in BRANCH
\end{split}\]
The decision variables are the active power output of the generators \(P_{i}\), reactive power output of the generators \(Q_{i}\), voltage magnitude of the buses \(V_{b}\), phase angle of the buses \(\theta_{b}\) for \(b \in BUS\) and the branch power flows \(P_{ij}\) and \(Q_{ij}\) for \((i, j) \in BRANCH\).
The objective function is the total cost of the generators, which is the sum of the quadratic cost, the linear cost, and the constant cost. The first constraint ensures that the phase angle of the reference bus is zero. The second and third constraints are the active and reactive power output limits of the generators. The fourth constraint is the voltage magnitude limits of the buses. The fifth and sixth constraints are the power balance equations of the buses. The seventh and eighth constraints are the power flow equations of the branches. The ninth and tenth constraints are the power flow equations of the branches in the opposite direction. The eleventh and twelfth constraints are the apparent power limits of the branches. The last constraint is the phase angle difference limits of the branches.
Implementation¶
We notice that the branch power flow equations with respect to \(P_{ij}\), \(Q_{ij}\), \(P_{ji}\) and \(Q_{ji}\) are nonlinear equations and have the same structure for each branch. We will use a single function to model the nonlinear power flow equations of one branch.
import math
import pyoptinterface as poi
from pyoptinterface import nlfunc, ipopt
model = ipopt.Model()
def branch_flow(vars, params):
G, B, Bc = params.G, params.B, params.Bc
Vi, Vj, theta_i, theta_j, Pij, Qij, Pji, Qji = (
vars.Vi,
vars.Vj,
vars.theta_i,
vars.theta_j,
vars.Pij,
vars.Qij,
vars.Pji,
vars.Qji,
)
sin_ij = nlfunc.sin(theta_i - theta_j)
cos_ij = nlfunc.cos(theta_i - theta_j)
Pij_eq = G * Vi**2 - Vi * Vj * (G * cos_ij + B * sin_ij) - Pij
Qij_eq = -(B + Bc) * Vi**2 - Vi * Vj * (G * sin_ij - B * cos_ij) - Qij
Pji_eq = G * Vj**2 - Vi * Vj * (G * cos_ij - B * sin_ij) - Pji
Qji_eq = -(B + Bc) * Vj**2 - Vi * Vj * (-G * sin_ij - B * cos_ij) - Qji
return [Pij_eq, Qij_eq, Pji_eq, Qji_eq]
bf = model.register_function(branch_flow)
Here the nonlinear function takes two arguments: vars and params. Although the power flow constraints take the same form for each branch, the parameters G B and Bc, namely the \(\pi\) circuit parameters of the branch, are different for each branch.
Next, we will use PJM 5-bus system as an example to demonstrate the implementation of the optimal power flow problem. The PJM 5-bus system is a small power system with 5 buses and 6 branches. The system data is shown below.
branches = [
# (from, to, R, X, B, angmin, angmax, Smax)
(0, 1, 0.00281, 0.0281, 0.00712, -30.0, 30.0, 4.00),
(0, 3, 0.00304, 0.0304, 0.00658, -30.0, 30.0, 4.26),
(0, 4, 0.00064, 0.0064, 0.03126, -30.0, 30.0, 4.26),
(1, 2, 0.00108, 0.0108, 0.01852, -30.0, 30.0, 4.26),
(2, 3, 0.00297, 0.0297, 0.00674, -30.0, 30.0, 4.26),
(3, 4, 0.00297, 0.0297, 0.00674, -30.0, 30.0, 2.40),
]
buses = [
# (Pd, Qd, Gs, Bs, Vmin, Vmax)
(0.0, 0.0000, 0.0, 0.0, 0.9, 1.1),
(3.0, 0.9861, 0.0, 0.0, 0.9, 1.1),
(3.0, 0.9861, 0.0, 0.0, 0.9, 1.1),
(4.0, 1.3147, 0.0, 0.0, 0.9, 1.1),
(0.0, 0.0000, 0.0, 0.0, 0.9, 1.1),
]
generators = [
# (bus, Pmin, Pmax, Qmin, Qmax, a, b, c)
(0, 0.0, 0.4, -0.300, 0.300, 0.0, 1400, 0.0),
(0, 0.0, 1.7, -1.275, 1.275, 0.0, 1500, 0.0),
(2, 0.0, 5.2, -3.900, 3.900, 0.0, 3000, 0.0),
(3, 0.0, 2.0, -1.500, 1.500, 0.0, 4000, 0.0),
(4, 0.0, 6.0, -4.500, 4.500, 0.0, 1000, 0.0),
]
slack_bus = 3
Then we declare the variables:
N_branch = len(branches)
N_bus = len(buses)
N_gen = len(generators)
Pbr_from = model.add_variables(range(N_branch))
Qbr_from = model.add_variables(range(N_branch))
Pbr_to = model.add_variables(range(N_branch))
Qbr_to = model.add_variables(range(N_branch))
V = model.add_variables(range(N_bus), name="V")
theta = model.add_variables(range(N_bus), name="theta")
for i in range(N_bus):
Vmin, Vmax = buses[i][4], buses[i][5]
model.set_variable_bounds(V[i], Vmin, Vmax)
model.set_variable_bounds(theta[slack_bus], 0.0, 0.0)
P = model.add_variables(range(N_gen), name="P")
Q = model.add_variables(range(N_gen), name="Q")
for i in range(N_gen):
model.set_variable_bounds(P[i], generators[i][1], generators[i][2])
model.set_variable_bounds(Q[i], generators[i][3], generators[i][4])
Next, we add the constraints:
# nonlinear constraints
for k in range(N_branch):
branch = branches[k]
R, X, Bc2 = branch[2], branch[3], branch[4]
G = R / (R**2 + X**2)
B = -X / (R**2 + X**2)
Bc = Bc2 / 2
i = branch[0]
j = branch[1]
Vi = V[i]
Vj = V[j]
theta_i = theta[i]
theta_j = theta[j]
Pij = Pbr_from[k]
Qij = Qbr_from[k]
Pji = Pbr_to[k]
Qji = Qbr_to[k]
model.add_nl_constraint(
bf,
vars=nlfunc.Vars(
Vi=Vi,
Vj=Vj,
theta_i=theta_i,
theta_j=theta_j,
Pij=Pij,
Qij=Qij,
Pji=Pji,
Qji=Qji,
),
params=nlfunc.Params(G=G, B=B, Bc=Bc),
eq=0.0,
)
# power balance constraints
P_balance_eq = [poi.ExprBuilder() for i in range(N_bus)]
Q_balance_eq = [poi.ExprBuilder() for i in range(N_bus)]
for b in range(N_bus):
Pd, Qd = buses[b][0], buses[b][1]
Gs, Bs = buses[b][2], buses[b][3]
Vb = V[b]
P_balance_eq[b] -= poi.quicksum(
Pbr_from[k] for k in range(N_branch) if branches[k][0] == b
)
P_balance_eq[b] -= poi.quicksum(
Pbr_to[k] for k in range(N_branch) if branches[k][1] == b
)
P_balance_eq[b] += poi.quicksum(P[i] for i in range(N_gen) if generators[i][0] == b)
P_balance_eq[b] -= Pd
P_balance_eq[b] -= Gs * Vb * Vb
Q_balance_eq[b] -= poi.quicksum(
Qbr_from[k] for k in range(N_branch) if branches[k][0] == b
)
Q_balance_eq[b] -= poi.quicksum(
Qbr_to[k] for k in range(N_branch) if branches[k][1] == b
)
Q_balance_eq[b] += poi.quicksum(Q[i] for i in range(N_gen) if generators[i][0] == b)
Q_balance_eq[b] -= Qd
Q_balance_eq[b] += Bs * Vb * Vb
model.add_quadratic_constraint(P_balance_eq[b], poi.Eq, 0.0)
model.add_quadratic_constraint(Q_balance_eq[b], poi.Eq, 0.0)
for k in range(N_branch):
branch = branches[k]
i = branch[0]
j = branch[1]
theta_i = theta[i]
theta_j = theta[j]
angmin = branch[5] / 180 * math.pi
angmax = branch[6] / 180 * math.pi
model.add_linear_constraint(theta_i - theta_j, poi.In, angmin, angmax)
Smax = branch[7]
Pij = Pbr_from[k]
Qij = Qbr_from[k]
Pji = Pbr_to[k]
Qji = Qbr_to[k]
model.add_quadratic_constraint(Pij * Pij + Qij * Qij, poi.Leq, Smax * Smax)
model.add_quadratic_constraint(Pji * Pji + Qji * Qji, poi.Leq, Smax * Smax)
Finally, we set the objective function:
cost = poi.ExprBuilder()
for i in range(N_gen):
a, b, c = generators[i][5], generators[i][6], generators[i][7]
cost += a * P[i] * P[i] + b * P[i] + c
model.set_objective(cost)
After optimization, we can retrieve the optimal solution:
model.optimize()
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 152
Number of nonzeros in inequality constraint Jacobian.: 33
Number of nonzeros in Lagrangian Hessian.............: 60
Total number of variables............................: 43
variables with only lower bounds: 0
variables with lower and upper bounds: 15
variables with only upper bounds: 0
Total number of equality constraints.................: 34
Total number of inequality constraints...............: 18
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 6
inequality constraints with only upper bounds: 12
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0059989e+02 3.99e+00 2.88e+01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 8.2642475e+03 2.47e+00 1.01e+02 -1.0 2.79e+00 - 4.08e-03 3.81e-01h 1
2 7.3953828e+03 2.42e+00 9.90e+01 -1.0 1.89e+01 - 6.33e-02 1.87e-02f 2
3 6.5016499e+03 2.36e+00 9.57e+01 -1.0 1.58e+01 - 2.57e-01 2.73e-02f 1
4 6.3495650e+03 2.27e+00 9.10e+01 -1.0 1.48e+01 - 4.57e-01 3.80e-02f 3
5 6.2696718e+03 2.07e+00 8.07e+01 -1.0 1.38e+01 - 7.59e-01 8.54e-02f 2
6 7.1293922e+03 1.63e+00 6.07e+01 -1.0 1.07e+01 - 8.60e-01 2.13e-01h 2
7 9.2308978e+03 9.26e-01 3.44e+01 -1.0 8.92e+00 - 4.35e-01 4.33e-01h 1
8 1.1477927e+04 6.73e-01 2.82e+01 -1.0 2.96e+00 - 1.27e-01 2.73e-01h 1
9 1.5885484e+04 1.86e-01 1.12e+01 -1.0 2.20e+00 - 5.73e-01 7.24e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.7582902e+04 7.44e-03 1.39e+00 -1.0 1.54e+00 - 8.05e-01 1.00e+00h 1
11 1.7581523e+04 5.37e-04 6.38e-02 -1.0 8.19e-01 - 1.00e+00 1.00e+00h 1
12 1.7556846e+04 4.26e-03 1.53e-01 -2.5 3.55e-01 - 8.73e-01 9.76e-01f 1
13 1.7552853e+04 3.50e-03 1.39e-01 -2.5 8.11e-02 - 1.00e+00 7.89e-01h 1
14 1.7552881e+04 1.84e-05 1.96e-03 -2.5 1.30e-02 - 1.00e+00 1.00e+00h 1
15 1.7551955e+04 1.39e-05 6.19e-03 -3.8 1.78e-02 - 1.00e+00 9.72e-01f 1
16 1.7551939e+04 1.54e-08 1.13e-06 -3.8 6.68e-04 - 1.00e+00 1.00e+00f 1
17 1.7551891e+04 1.33e-08 1.02e-06 -5.7 9.92e-04 - 1.00e+00 1.00e+00f 1
18 1.7551891e+04 1.15e-12 1.20e-10 -8.6 1.29e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 18
(scaled) (unscaled)
Objective...............: 4.3879727096482048e+02 1.7551890838592819e+04
Dual infeasibility......: 1.1970113789061543e-10 4.7880455156246171e-09
Constraint violation....: 1.1473044736476368e-12 1.1473044736476368e-12
Complementarity.........: 2.5405436222116973e-09 1.0162174488846789e-07
Overall NLP error.......: 2.5405436222116973e-09 1.0162174488846789e-07
Number of objective function evaluations = 33
Number of objective gradient evaluations = 19
Number of equality constraint evaluations = 33
Number of inequality constraint evaluations = 33
Number of equality constraint Jacobian evaluations = 19
Number of inequality constraint Jacobian evaluations = 19
Number of Lagrangian Hessian evaluations = 18
Total CPU secs in IPOPT (w/o function evaluations) = 0.009
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
print(model.get_model_attribute(poi.ModelAttribute.TerminationStatus))
P_value = P.map(model.get_value)
print("Optimal active power output of the generators:")
for i in range(N_gen):
print(f"Generator {i}: {P_value[i]}")
TerminationStatusCode.LOCALLY_SOLVED
Optimal active power output of the generators:
Generator 0: 0.4
Generator 1: 1.7
Generator 2: 3.2449849431098143
Generator 3: 0.0
Generator 4: 4.7069359970171085