Optimal Power Flow with AMPL and Python - Bus Injection Model (BIM) with controllable-phase shifting transformers and tap-changing transformers#
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, \varphi, T) = P_i^G - P_i^L, \forall i \in N \\
& Q_i(V, \delta, \varphi, T) = 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 \\
& \varphi_i^{min} \leq \varphi_i \leq \varphi_i^{max}, \forall i \in L \\
& T_i^{min} \leq T_i \leq T_i^{max}, \forall i \in L \\
\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 opfadv.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 T0 {L}; # initial voltage ratio
param phi0 {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
# lower and upper bounds
param T_min {L};
param T_max {L};
param phi_min {L};
param phi_max {L};
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
# variables
# 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
var T {(i,k,m) in L} >= T_min[i,k,m], <= T_max[i,k,m] := T0[i,k,m]; # voltage ratio
var phi {(i,k,m) in L} >= phi_min[i,k,m], <= phi_max[i,k,m] := phi0[i,k,m]; # phase angle
# auxiliary variables
# bus admittance matrix real part
var 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
var 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]
);
# 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 opfadv.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")
h = pd.DataFrame(
[
[4, 3, 4, -30, 30]
],
columns=["row", "from", "to", "phi_min", "phi_max"]
).set_index(["row", "from", "to"])
k = pd.DataFrame(
[
[5, 3, 5, 0.95, 1.05]
],
columns=["row", "from", "to", "T_min", "T_max"]
).set_index(["row", "from", "to"])
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_branch = ampl_branch.rename({"T":"T0", "phi":"phi0"}, axis=1)
ampl_branch["T_min"] = ampl_branch["T0"]
ampl_branch["T_max"] = ampl_branch["T0"]
ampl_branch["phi_min"] = ampl_branch["phi0"]
ampl_branch["phi_max"] = ampl_branch["phi0"]
ampl_branch.loc[(k.index), ["T_min", "T_max"]] = k
ampl_branch.loc[(h.index), ["phi_min", "phi_max"]] = h
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["phi0"] = ampl_branch["phi0"].apply(radians)
ampl_branch["phi_min"] = ampl_branch["phi_min"].apply(radians)
ampl_branch["phi_max"] = ampl_branch["phi_max"].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 T0 phi0 T_min T_max phi_min \
row from to
1 1 2 0.000 0.300 0.0 0.00 1.00 0.00000 1.00 1.00 0.000000
2 1 3 0.023 0.145 0.0 0.04 1.00 0.00000 1.00 1.00 0.000000
3 2 4 0.006 0.032 0.0 0.01 1.00 0.00000 1.00 1.00 0.000000
4 3 4 0.020 0.260 0.0 0.00 1.00 -0.05236 1.00 1.00 -0.523599
5 3 5 0.000 0.320 0.0 0.00 0.98 0.00000 0.95 1.05 0.000000
6 4 5 0.000 0.500 0.0 0.00 1.00 0.00000 1.00 1.00 0.000000
phi_max
row from to
1 1 2 0.000000
2 1 3 0.000000
3 2 4 0.000000
4 3 4 0.523599
5 3 5 0.000000
6 4 5 0.000000
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 opfadv_run(bus, branch, gen, ref):
# initialyze AMPL and read the model
ampl = AMPL()
ampl.read("opfadv.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
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()
bus_res = ampl.get_data("V", "delta").to_pandas()
branch_res = ampl.get_data("T", "phi").to_pandas()
gen_res = ampl.get_data("PG", "QG").to_pandas()
return obj, bus_res, branch_res, gen_res
obj, ma, pt, pg = opfadv_run(ampl_bus, ampl_branch, ampl_gen, ref)
# convert radians back to degrees
ma["delta"] = ma["delta"].apply(degrees)
pt["phi"] = pt["phi"].apply(degrees)
print("generation cost:", obj)
print(ma)
print(pt)
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...: 70
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 55
Total number of variables............................: 16
variables with only lower bounds: 0
variables with lower and upper bounds: 14
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.0608254e-01 3.89e-02 5.32e-01 -1.0 8.36e-01 - 9.78e-01 1.00e+00f 1
2 4.0912162e-01 7.90e-04 1.85e-02 -1.7 5.41e-02 - 1.00e+00 1.00e+00h 1
3 4.0618651e-01 1.70e-04 1.84e-02 -2.5 3.18e-02 - 9.83e-01 1.00e+00h 1
4 4.0306632e-01 2.17e-03 5.31e-03 -3.8 6.81e-02 - 9.38e-01 1.00e+00h 1
5 4.0188992e-01 4.78e-03 2.00e-04 -3.8 9.15e-02 - 1.00e+00 1.00e+00h 1
6 4.0194463e-01 2.52e-04 1.14e-05 -3.8 3.07e-02 - 1.00e+00 1.00e+00h 1
7 4.0172007e-01 6.48e-04 1.43e-03 -5.7 7.34e-02 - 7.79e-01 1.00e+00h 1
8 4.0166464e-01 3.97e-04 3.95e-05 -5.7 6.21e-02 - 9.97e-01 1.00e+00h 1
9 4.0166242e-01 8.53e-05 2.08e-06 -5.7 2.97e-02 - 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.0166299e-01 1.11e-05 6.07e-07 -5.7 5.88e-03 - 1.00e+00 1.00e+00h 1
11 4.0165991e-01 1.52e-05 4.75e-05 -8.6 7.30e-03 - 9.03e-01 1.00e+00h 1
12 4.0165958e-01 4.30e-06 3.94e-08 -8.6 9.72e-04 - 1.00e+00 1.00e+00h 1
13 4.0165956e-01 6.09e-09 4.53e-11 -8.6 3.62e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 4.0165956437302702e-01 4.0165956437302702e-01
Dual infeasibility......: 4.5338876966059241e-11 4.5338876966059241e-11
Constraint violation....: 6.0877359264743802e-09 6.0877359264743802e-09
Complementarity.........: 2.5163793223846341e-09 2.5163793223846341e-09
Overall NLP error.......: 6.0877359264743802e-09 6.0877359264743802e-09
Number of objective function evaluations = 14
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 14
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 13
Total CPU secs in IPOPT (w/o function evaluations) = 0.003
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
Ipopt 3.12.13: Optimal Solution Found
suffix ipopt_zU_out OUT;
suffix ipopt_zL_out OUT;
generation cost: 0.401659564373027
V delta
1 1.000000 0.000000
2 0.980824 -12.584346
3 0.956703 -1.671734
4 0.967651 -13.859993
5 0.958942 -9.133640
T phi
index0 index1 index2
1 1 2 1.000000 0.000000
2 1 3 1.000000 0.000000
3 2 4 1.000000 0.000000
4 3 4 1.000000 12.375281
5 3 5 0.950005 0.000000
6 4 5 1.000000 0.000000
PG QG
1 0.946724 0.386688
3 0.191521 -0.126627
4 0.053072 0.199994
Conclusion#
Bibliography#
Stephen Frank & Steffen Rebennack (2016) An introduction to optimal power flow: Theory, formulation, and examples, IIE Transactions, 48:12, 1172-1197.