Optimal Power Flow with AMPL and Python - DC Power Flow#

opf3.ipynb Open In Colab Kaggle Gradient Open In SageMaker Studio Lab Hits

Description: Optimal Power Flow

Tags: AMPL, amplpy, Optimal Power Flow, Python

Notebook author: Nicolau Santos <nicolau@ampl.com>

# Install dependencies
%pip install -q amplpy pandas

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["coin"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Introduction#

Content will be available soon!

Problem description#

DC Power Flow assumptions

  • \(G_{ik}=0\)

  • \(\sin(\delta_i-\delta_k) \approx \delta_i-\delta_k\) and \(\cos(\delta_i-\delta_k) \approx 1\)

  • \(V_i \approx 1, \forall i \in N\)

  • Reactive Power Flow is neglected

\[ P_i(\delta) \approx \sum_{k=1}^{N}B_{ik}(\delta_i-\delta_k) \]
\[\begin{split} \begin{split} \min \enspace & \sum_{i \in G}(\text{const} + \text{linear}P_i^G + \text{quad}(P_i^G)^2) \\ \text{s.t.} \enspace & P_i(\delta) = P_i^G - P_i^L, \forall i \in N \\ & P_i^{G,min} \leq P_i^{G} \leq P_i^{G,max}, \forall i \in G \\ & \delta_i^{min} \leq \delta_i \leq \delta_i^{max}, \forall i \in N \\ \end{split} \end{split}\]

AMPL model#

%%writefile dcopf.mod
# data

set N;                              # set of buses in the network
param nL;                           # number of branches in the network
set L within 1..nL cross N cross N; # set of branches in the network
set GEN within N;                   # set of generator buses
set REF within N;                   # set of reference (slack) buses
set YN :=                           # index of the bus admittance matrix
    setof {i in N} (i,i) union
    setof {(i,k,l) in L} (k,l) union
    setof {(i,k,l) in L} (l,k);

# bus data

param delta0    {N}; # initial voltage angle
param PL        {N}; # real power load
param b_s       {N}; # shunt susceptance

# lower and upper bounds
param delta_min {N};
param delta_max {N};


# generator data

param PG0    {GEN}; # initial real power generation
param const  {GEN}; # constant cost of a given generator
param linear {GEN}; # linear cost of a given generator
param quad   {GEN}; # quadratic cost of a given generator

# lower and upper bounds
param PG_min {GEN};
param PG_max {GEN};


# branch data

param T    {L}; # initial voltage ratio
param phi  {L}; # initial phase angle
param R    {L}; # branch resistance
param X    {L}; # branch reactance
param b_sh {L}; # shunt susceptance

param g {(l,k,m) in L} := R[l,k,m]/(R[l,k,m]^2 + X[l,k,m]^2);  # series conductance
param b {(l,k,m) in L} := -X[l,k,m]/(R[l,k,m]^2 + X[l,k,m]^2); # series susceptance


# bus admittance matrix imaginary part
param B {(i,k) in YN} =
    if (i == k) then (
        b_s[i] +
        sum{(l,i,u) in L} (b[l,i,u] + b_sh[l,i,u]/2)/T[l,i,u]**2 +
        sum{(l,u,i) in L} (b[l,u,i] + b_sh[l,u,i]/2)
    )
    else (
        -sum{(l,i,k) in L} (
            g[l,i,k]*sin(phi[l,i,k])+b[l,i,k]*cos(phi[l,i,k])
        )/T[l,i,k] -
        sum{(l,k,i) in L} (
            -g[l,k,i]*sin(phi[l,k,i])+b[l,k,i]*cos(phi[l,k,i])
        )/T[l,k,i]
    );


# decision variables with lower bounds, upper bounds and initial guess

var delta {i in N} >= delta_min[i], <= delta_max[i] := delta0[i]; # voltage angle

var PG {i in GEN} >= PG_min[i], <= PG_max[i] := PG0[i]; # real power generation


# real power injection
var P {i in N} = sum{(i,k) in YN} B[i,k] * (delta[i] - delta[k]);


# objective

minimize generation_cost:
    sum{i in GEN} (const[i] + PG[i] * linear[i] + (PG[i] ** 2) * quad[i]);


# constraints

s.t. p_flow {i in N}:
    P[i] == (if (i in GEN) then PG[i] else 0) - PL[i];

s.t. fixed_angles {i in REF}:
    delta[i] == delta0[i];

Writing dcopf.mod

Numerical example#

import pandas as pd
import math

inf = float("inf")

df_bus = pd.DataFrame(
    [
        [1, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0],
        [2, 0.0, 0.0, 0.0, 0.3, 0.95, 1.05],
        [3, 0.0, 0.0, 0.05, 0.0, 0.95, 1.05],
        [4, 0.9, 0.4, 0.0, 0.0, 0.95, 1.05],
        [5, 0.239, 0.129, 0.0, 0.0, 0.95, 1.05],
    ],
    columns=["bus", "PL", "QL", "g_s", "b_s", "V_min", "V_max"],
).set_index("bus")

df_branch = pd.DataFrame(
    [
        [1, 1, 2, 0.0, 0.3, 0.0, 0.0, 1.0, 0.0],
        [2, 1, 3, 0.023, 0.145, 0.0, 0.04, 1.0, 0.0],
        [3, 2, 4, 0.006, 0.032, 0.0, 0.01, 1.0, 0.0],
        [4, 3, 4, 0.02, 0.26, 0.0, 0.0, 1.0, -3.0],
        [5, 3, 5, 0.0, 0.32, 0.0, 0.0, 0.98, 0.0],
        [6, 4, 5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0],
    ],
    columns=["row", "from", "to", "R", "X", "g_sh", "b_sh", "T", "phi"],
).set_index(["row", "from", "to"])

df_gen = pd.DataFrame(
    [
        [1, float(-inf), float(inf), float(-inf), float(inf), 0.0, 0.35, 0.0],
        [3, 0.10, 0.40, -0.20, 0.30, 0.0, 0.20, 0.40],
        [4, 0.05, 0.40, -0.20, 0.20, 0.0, 0.30, 0.50],
    ],
    columns=["bus", "PG_min", "PG_max", "QG_min", "QG_max", "const", "linear", "quad"],
).set_index("bus")


ref = [1]

# print(df_bus)
# print(df_branch)
# print(df_gen)

# data preprocessing

ampl_bus = pd.DataFrame()
cols = ["PL", "b_s"]
for col in cols:
    ampl_bus.loc[:, col] = df_bus.loc[:, col]
ampl_bus["delta0"] = 0.0
ampl_bus["delta_min"] = -180.0
ampl_bus["delta_max"] = 180.0

ampl_branch = pd.DataFrame()
ampl_branch = df_branch.copy()
ampl_branch = ampl_branch.drop("g_sh", axis=1)

ampl_gen = df_gen.copy()
ampl_gen["PG0"] = 0.0
ampl_gen = ampl_gen.drop(["QG_min", "QG_max"], axis=1)

# convert degrees to radians
ampl_bus["delta0"] = ampl_bus["delta0"].apply(math.radians)
ampl_bus["delta_min"] = ampl_bus["delta_min"].apply(math.radians)
ampl_bus["delta_max"] = ampl_bus["delta_max"].apply(math.radians)
ampl_branch["phi"] = ampl_branch["phi"].apply(math.radians)

print(ampl_bus)
print(ampl_branch)
print(ampl_gen)
print(ref)

        PL  b_s  delta0  delta_min  delta_max
bus                                          
1    0.000  0.0     0.0  -3.141593   3.141593
2    0.000  0.3     0.0  -3.141593   3.141593
3    0.000  0.0     0.0  -3.141593   3.141593
4    0.900  0.0     0.0  -3.141593   3.141593
5    0.239  0.0     0.0  -3.141593   3.141593
                 R      X  b_sh     T      phi
row from to                                   
1   1    2   0.000  0.300  0.00  1.00  0.00000
2   1    3   0.023  0.145  0.04  1.00  0.00000
3   2    4   0.006  0.032  0.01  1.00  0.00000
4   3    4   0.020  0.260  0.00  1.00 -0.05236
5   3    5   0.000  0.320  0.00  0.98  0.00000
6   4    5   0.000  0.500  0.00  1.00  0.00000
     PG_min  PG_max  const  linear  quad  PG0
bus                                          
1      -inf     inf    0.0    0.35   0.0  0.0
3      0.10     0.4    0.0    0.20   0.4  0.0
4      0.05     0.4    0.0    0.30   0.5  0.0
[1]

def dcopf_run(bus, branch, gen, ref):
    # initialyze AMPL and read the model
    ampl = AMPL()
    ampl.read("dcopf.mod")

    # load the data
    ampl.set_data(bus, "N")
    ampl.param["nL"] = branch.shape[0]
    ampl.set_data(branch, "L")
    ampl.set_data(gen, "GEN")
    ampl.set["REF"] = ref

    # uncomment to show expanded problem
    # ampl.eval("solexpand;")

    # uncomment to show admittance matrix
    # ampl.eval("display G,B;")

    ampl.solve(solver="ipopt")
    assert ampl.solve_result == "solved", ampl.solve_result

    if ampl.solve_result != "solved":
        print("WARNING: solver exited with %s status." % (ampl.solve_result,))

    obj = ampl.obj["generation_cost"].value()
    ma = ampl.get_data("delta").to_pandas()
    pg = ampl.get_data("PG").to_pandas()

    return obj, ma, pg


obj, ma, pg = dcopf_run(ampl_bus, ampl_branch, ampl_gen, ref)

# convert radians back to degrees
ma["delta"] = ma["delta"].apply(math.degrees)

# print results
print("generation cost:", obj)
print(ma)
print(pg)

Ipopt 3.12.13: 

******************************************************************************
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 is Ipopt version 3.12.13, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:       17
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

Total number of variables............................:        7
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        6
                     variables with only upper bounds:        0
Total number of equality constraints.................:        5
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.2324719e-02 8.47e-01 4.54e-02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.8944604e-01 6.66e-16 1.97e-03  -1.0 7.89e-01    -  9.78e-01 1.00e+00f  1
   2  3.8974111e-01 1.11e-16 7.79e-16  -1.7 9.79e-03    -  1.00e+00 1.00e+00f  1
   3  3.8688794e-01 4.44e-16 3.49e-03  -3.8 3.83e-02    -  8.90e-01 1.00e+00f  1
   4  3.8501142e-01 1.11e-16 2.68e-15  -3.8 4.48e-02    -  1.00e+00 1.00e+00f  1
   5  3.8450538e-01 4.44e-16 8.26e-04  -5.7 2.30e-02    -  9.55e-01 1.00e+00f  1
   6  3.8438280e-01 1.11e-15 1.26e-15  -5.7 9.33e-03    -  1.00e+00 1.00e+00f  1
   7  3.8435155e-01 1.44e-15 2.05e-15  -5.7 4.68e-03    -  1.00e+00 1.00e+00f  1
   8  3.8434431e-01 8.88e-16 2.21e-15  -5.7 2.23e-03    -  1.00e+00 1.00e+00f  1
   9  3.8434240e-01 3.33e-16 2.95e-06  -8.6 1.30e-03    -  9.98e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  3.8434207e-01 1.22e-15 3.01e-15  -8.6 5.92e-04    -  1.00e+00 1.00e+00f  1
  11  3.8434203e-01 6.66e-16 2.44e-15  -8.6 2.22e-04    -  1.00e+00 1.00e+00h  1
  12  3.8434203e-01 3.33e-16 9.28e-16  -8.6 4.31e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:   3.8434203386092713e-01    3.8434203386092713e-01
Dual infeasibility......:   9.2775847754871432e-16    9.2775847754871432e-16
Constraint violation....:   3.3306690738754696e-16    3.3306690738754696e-16
Complementarity.........:   4.3487950483208504e-09    4.3487950483208504e-09
Overall NLP error.......:   4.3487950483208504e-09    4.3487950483208504e-09


Number of objective function evaluations             = 13
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 13
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 12
Total CPU secs in IPOPT (w/o function evaluations)   =      0.013
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
 
Ipopt 3.12.13: Optimal Solution Found

suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
generation cost: 0.3843420338609271
      delta
1  0.000000
2 -8.117728
3 -3.676164
4 -9.014060
5 -8.372741
         PG
1  0.903901
3  0.186916
4  0.051051

Conclusion#

Bibliography#

Stephen Frank & Steffen Rebennack (2016) An introduction to optimal power flow: Theory, formulation, and examples, IIE Transactions, 48:12, 1172-1197.