Optimal Power Flow with AMPL and Python - conventional Power Flow#
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#
\[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(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 pf.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 PQ within N; # set of load buses
set PV within N; # set of voltage-controlled 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
# generator data
param PG {GEN}; # real power generation
param QG {GEN}; # reactive power generation
# branch indexed 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]
);
# variables
var V {i in N} := V0[i]; # voltage magnitude
var delta {i in N} := delta0[i]; # voltage angle
# 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])
);
# constraints
s.t. p_flow {i in (PQ union PV)}:
P[i] == (if (i in GEN) then PG[i] else 0) - PL[i];
s.t. q_flow {i in PQ}:
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 union PV)}:
V[i] == V0[i];
Overwriting pf.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"])
gen = [1, 3, 4]
ref = [1]
pq = [3, 4]
pv = [2, 5]
#print(df_bus)
#print(df_branch)
# data preprocessing
ampl_bus = pd.DataFrame()
cols = ["PL", "QL", "g_s", "b_s"]
for col in cols:
ampl_bus.loc[:, col] = df_bus.loc[:, col]
ampl_bus["V0"] = 1.0
ampl_bus["delta0"] = 0.0
ampl_branch = pd.DataFrame()
ampl_branch = df_branch.copy()
ampl_gen = pd.DataFrame()
ampl_gen.index = gen
ampl_gen["PG"] = 0.0
ampl_gen["QG"] = 0.0
# convert degrees to radians
ampl_bus["delta0"] = ampl_bus["delta0"].apply(radians)
ampl_branch["phi"] = ampl_branch["phi"].apply(radians)
print(ampl_bus)
print(ampl_branch)
print(ampl_gen)
PL QL g_s b_s V0 delta0
bus
1 0.000 0.000 0.00 0.0 1.0 0.0
2 0.000 0.000 0.00 0.3 1.0 0.0
3 0.000 0.000 0.05 0.0 1.0 0.0
4 0.900 0.400 0.00 0.0 1.0 0.0
5 0.239 0.129 0.00 0.0 1.0 0.0
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 QG
1 0.0 0.0
3 0.0 0.0
4 0.0 0.0
def pf_run(bus, branch, gen, ref, pq, pv):
# initialyze AMPL and read the model
ampl = AMPL()
ampl.read("pf.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.set["PQ"] = pq
ampl.set["PV"] = pv
ampl.solve(solver=SOLVER)
solve_result = ampl.get_value("solve_result")
if solve_result != "solved":
print("WARNING: solver exited with %s status." %(solve_result,))
return ampl.get_data("V", "delta").to_pandas(), solve_result
df_res, solver_status = pf_run(ampl_bus, ampl_branch, ampl_gen, ref, pq, pv)
# convert radians back to degrees
df_res["delta"] = df_res["delta"].apply(degrees)
# print results
print("solver status:", solver_status)
print(df_res)
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...: 30
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 18
Total number of variables............................: 6
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 6
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 0.0000000e+00 7.00e-01 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 4.87e-02 0.00e+00 -1.7 1.71e-01 - 1.00e+00 1.00e+00h 1
2 0.0000000e+00 2.17e-04 0.00e+00 -2.5 4.08e-03 - 1.00e+00 1.00e+00h 1
3 0.0000000e+00 4.11e-09 0.00e+00 -5.7 1.76e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 3
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 4.1074473979318246e-09 4.1074473979318246e-09
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 4.1074473979318246e-09 4.1074473979318246e-09
Number of objective function evaluations = 4
Number of objective gradient evaluations = 4
Number of equality constraint evaluations = 4
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 4
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 3
Total CPU secs in IPOPT (w/o function evaluations) = 0.002
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;
solver status: solved
V delta
1 1.000000 0.000000
2 1.000000 -8.657929
3 0.981536 -5.893046
4 0.983056 -9.440548
5 1.000000 -9.950946
Conclusion#
Bibliography#
Stephen Frank & Steffen Rebennack (2016) An introduction to optimal power flow: Theory, formulation, and examples, IIE Transactions, 48:12, 1172-1197.