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=["highs"],  # 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];

Overwriting dcopf.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", "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(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  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=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("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(degrees)

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

HiGHS 1.7.1:HiGHS 1.7.1: optimal solution; objective 0.3843420339
0 simplex iterations
0 barrier iterations
generation cost: 0.3843420338524032
      delta
1  0.000000
2 -8.117763
3 -3.676181
4 -9.014099
5 -8.372766
         PG
1  0.903905
3  0.186916
4  0.051047

Conclusion#

Bibliography#

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