optimal_power_flow.md

file_format	mystnb
kernelspec	name python3

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

$$ \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 $$

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 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:

import math
import pyoptinterface as poi
from pyoptinterface import nl, ipopt

model = ipopt.Model()

N_branch = len(branches)
N_bus = len(buses)
N_gen = len(generators)

Pbr_from = model.add_m_variables(N_branch)
Qbr_from = model.add_m_variables(N_branch)
Pbr_to = model.add_m_variables(N_branch)
Qbr_to = model.add_m_variables(N_branch)

V = model.add_m_variables(N_bus, name="V")
theta = model.add_m_variables(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):
    with nl.graph():
        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]

        sin_ij = nl.sin(theta_i - theta_j)
        cos_ij = nl.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

        model.add_nl_constraint(
            Pij_eq == 0.0,
        )
        model.add_nl_constraint(
            Qij_eq == 0.0,
        )
        model.add_nl_constraint(
            Pji_eq == 0.0,
        )
        model.add_nl_constraint(
            Qji_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, (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()

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]}")