Optimal Power Flow with AMPL and Python - Bus Injection Model (BIM)#
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#
\[\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(V, \delta) = P_i^G - P_i^L, \forall i \in N \\
& Q_i(V, \delta) = Q_i^G - Q_i^L, \forall i \in N \\
& P_i^{G,min} \leq P_i^{G} \leq P_i^{G,max}, \forall i \in G \\
& Q_i^{G,min} \leq Q_i^{G} \leq Q_i^{G,max}, \forall i \in G \\
& V_i^{min} \leq P_i \leq V_i^{max}, \forall i \in N \\
& \delta_i^{min} \leq \delta_i \leq \delta_i^{max}, \forall i \in N \\
\end{split}
\end{split}\]
\[P_i(V, \delta) = V_i \sum_{k=1}^{N}V_k(G_{ik}\cos(\delta_i-\delta_k) + B_{ik}\sin(\delta_i-\delta_k)), \forall i \in N\]
\[Q_i(V, \delta) = V_i \sum_{k=1}^{N}V_k(G_{ik}\sin(\delta_i-\delta_k) - B_{ik}\cos(\delta_i-\delta_k)), \forall i \in N\]
AMPL model#
%%writefile opf.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 V0 {N}; # initial voltage magnitude
param delta0 {N}; # initial voltage angle
param PL {N}; # real power load
param QL {N}; # reactive power load
param g_s {N}; # shunt conductance
param b_s {N}; # shunt susceptance
# lower and upper bounds
param V_min {N};
param V_max {N};
param delta_min {N};
param delta_max {N};
# generator data
param PG0 {GEN}; # initial real power generation
param QG0 {GEN}; # initial reactive 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};
param QG_min {GEN};
param QG_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 g_sh {L}; # shunt conductance
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 real part
param G {(i,k) in YN} =
if (i == k) then (
g_s[i] +
sum{(l,i,u) in L} (g[l,i,u] + g_sh[l,i,u]/2)/T[l,i,u]**2 +
sum{(l,u,i) in L} (g[l,u,i] + g_sh[l,u,i]/2)
)
else (
-sum{(l,i,k) in L} ((
g[l,i,k]*cos(phi[l,i,k])-b[l,i,k]*sin(phi[l,i,k])
)/T[l,i,k]) -
sum{(l,k,i) in L} ((
g[l,k,i]*cos(phi[l,k,i])+b[l,k,i]*sin(phi[l,k,i])
)/T[l,k,i])
);
# 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 V {i in N} >= V_min[i], <= V_max[i] := V0[i]; # voltage magnitude
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
var QG {i in GEN} >= QG_min[i], <= QG_max[i] := QG0[i]; # reactive power generation
# real power injection
var P {i in N} =
V[i] * sum{(i,k) in YN} V[k] * (
G[i,k] * cos(delta[i] - delta[k]) +
B[i,k] * sin(delta[i] - delta[k])
);
# reactive power injection
var Q {i in N} =
V[i] * sum{(i,k) in YN} V[k] * (
G[i,k] * sin(delta[i] - delta[k]) -
B[i,k] * cos(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. q_flow {i in N}:
Q[i] == (if (i in GEN) then QG[i] else 0) - QL[i];
s.t. fixed_angles {i in REF}:
delta[i] == delta0[i];
s.t. fixed_voltages {i in REF}:
V[i] == V0[i];
Overwriting opf.mod
Numerical example#
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", "QL", "g_s", "b_s", "V_min", "V_max"]
for col in cols:
ampl_bus.loc[:, col] = df_bus.loc[:, col]
ampl_bus["V0"] = 1.0
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_gen = df_gen.copy()
ampl_gen["PG0"] = 0.0
ampl_gen["QG0"] = 0.0
# convert degrees to radians
ampl_bus["delta0"] = ampl_bus["delta0"].apply(radians)
ampl_bus["delta_min"] = ampl_bus["delta_min"].apply(radians)
ampl_bus["delta_max"] = ampl_bus["delta_max"].apply(radians)
ampl_branch["phi"] = ampl_branch["phi"].apply(radians)
print(ampl_bus)
print(ampl_branch)
print(ampl_gen)
print(ref)
PL QL g_s b_s V_min V_max V0 delta0 delta_min delta_max
bus
1 0.000 0.000 0.00 0.0 1.00 1.00 1.0 0.0 -3.141593 3.141593
2 0.000 0.000 0.00 0.3 0.95 1.05 1.0 0.0 -3.141593 3.141593
3 0.000 0.000 0.05 0.0 0.95 1.05 1.0 0.0 -3.141593 3.141593
4 0.900 0.400 0.00 0.0 0.95 1.05 1.0 0.0 -3.141593 3.141593
5 0.239 0.129 0.00 0.0 0.95 1.05 1.0 0.0 -3.141593 3.141593
R X g_sh b_sh T phi
row from to
1 1 2 0.000 0.300 0.0 0.00 1.00 0.00000
2 1 3 0.023 0.145 0.0 0.04 1.00 0.00000
3 2 4 0.006 0.032 0.0 0.01 1.00 0.00000
4 3 4 0.020 0.260 0.0 0.00 1.00 -0.05236
5 3 5 0.000 0.320 0.0 0.00 0.98 0.00000
6 4 5 0.000 0.500 0.0 0.00 1.00 0.00000
PG_min PG_max QG_min QG_max const linear quad PG0 QG0
bus
1 -inf inf -inf inf 0.0 0.35 0.0 0.0 0.0
3 0.10 0.4 -0.2 0.3 0.0 0.20 0.4 0.0 0.0
4 0.05 0.4 -0.2 0.2 0.0 0.30 0.5 0.0 0.0
[1]
def opf_run(bus, branch, gen, ref):
# initialyze AMPL and read the model
ampl = AMPL()
ampl.read("opf.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=SOLVER)
solve_result = ampl.get_value("solve_result")
if solve_result != "solved":
print("WARNING: solver exited with %s status." %(solve_result,))
obj = ampl.obj["generation_cost"].value()
ma = ampl.get_data("V", "delta").to_pandas()
pg = ampl.get_data("PG", "QG").to_pandas()
return obj, ma, pg
obj, ma, pg = opf_run(ampl_bus, ampl_branch, ampl_gen, ref)
# convert radians back to degrees
ma["delta"] = ma["delta"].apply(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...: 62
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 38
Total number of variables............................: 14
variables with only lower bounds: 0
variables with lower and upper bounds: 12
variables with only upper bounds: 0
Total number of equality constraints.................: 10
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 6.47e-01 4.52e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 4.0623606e-01 3.61e-02 5.17e-01 -1.0 8.37e-01 - 9.80e-01 1.00e+00f 1
2 4.1018715e-01 1.30e-04 1.37e-02 -1.7 6.02e-02 - 1.00e+00 1.00e+00h 1
3 4.0744338e-01 5.79e-05 2.05e-02 -2.5 3.18e-02 - 9.82e-01 1.00e+00h 1
4 4.0511863e-01 1.50e-04 2.57e-03 -3.8 5.31e-02 - 9.83e-01 1.00e+00h 1
5 4.0456520e-01 1.92e-04 6.80e-05 -3.8 4.51e-02 - 1.00e+00 1.00e+00h 1
6 4.0424502e-01 7.31e-04 2.73e-03 -5.7 8.44e-02 - 7.20e-01 1.00e+00h 1
7 4.0416214e-01 3.04e-04 1.39e-04 -5.7 5.66e-02 - 9.80e-01 1.00e+00h 1
8 4.0414750e-01 3.11e-05 2.23e-06 -5.7 1.97e-02 - 1.00e+00 1.00e+00h 1
9 4.0414754e-01 7.05e-08 1.56e-08 -5.7 8.23e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 4.0414383e-01 1.15e-06 2.09e-06 -8.6 3.49e-03 - 9.92e-01 1.00e+00h 1
11 4.0414383e-01 1.52e-10 3.74e-11 -8.6 4.41e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 4.0414383068667636e-01 4.0414383068667636e-01
Dual infeasibility......: 3.7426345610312982e-11 3.7426345610312982e-11
Constraint violation....: 1.5229577321473897e-10 1.5229577321473897e-10
Complementarity.........: 2.6076022452054029e-09 2.6076022452054029e-09
Overall NLP error.......: 2.6076022452054029e-09 2.6076022452054029e-09
Number of objective function evaluations = 12
Number of objective gradient evaluations = 12
Number of equality constraint evaluations = 12
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 12
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 11
Total CPU secs in IPOPT (w/o function evaluations) = 0.004
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.40414383068667636
V delta
1 1.000000 0.000000
2 0.982697 -7.498240
3 0.964150 -4.223082
4 0.969636 -8.204157
5 0.950000 -8.640602
PG QG
1 0.946145 0.248706
3 0.194654 -0.072208
4 0.057509 0.199998
Conclusion#
Bibliography#
Stephen Frank & Steffen Rebennack (2016) An introduction to optimal power flow: Theory, formulation, and examples, IIE Transactions, 48:12, 1172-1197.