Optimal Power Flow and Locational Marginal Prices

Optimal Power Flow and Locational Marginal Prices#

A fundamental component of the power system is the transmission network that allows the flow of power from generating resources to loads. A simple linear optimal power flow model is stated below.

1. Import packages

  • We are using the gurobipy package to formulate a mathematical model and solve it.

# import required packages

import pandas as pd
import numpy as np
import math

import gurobipy as gp
from gurobipy import GRB

import matplotlib.pyplot as plt

2. Define data

Symbol

Description

Data

c

Generation cost

[5, 20,100]

\(\overline{P}_i\)

Capacaity limit

[102, 100, 100]

\(\overline{F}_l\)

Flow capacity limit

[100, 100, 30]

D

Load

[0, 0, 130]

L

Number of Lines

3

G

Number of Generators

3

G

Number of Buses

3

\(B_{l,j}\)

Power Transfer Distribution Factor Matrix

\begin{bmatrix} 1/3 & -1/3 & 0 \ 2/3 & 1/3 & 0 \ 1/3 & 2/3 & 0 \end{bmatrix}

 
# problem parameters
c = [5, 20, 100] # generation cost
pmax = [102, 100, 100] # maximum generation
fmax = [100, 100, 30] # maximum flow

# Case 1
D = [0, 0, 130] # nodal demand
# Case 2
# D = [0, 30, 130] # nodal demand

n_bus = 3 # number of buses/nodes
n_gen = 3 # number of generators
n_line = 3 # number of lines

# ptdf matrix 
B = np.array([[1/3, -1/3, 0], 
              [2/3,  1/3, 0],
              [1/3,  2/3, 0]])

3. Define mathematical model

Objective function:

  • Minimize the cost of generation given by cost function \(c_i(p_i)\).

Decision variables:

  • \(p_{i}^{\text{c}}\) production of generator i.

  • \(f_{l}\) flow on line l.

Constraints:

  • Energy balance: the total sum of generation needs to equal the load.

  • Capacity constraints: each generator is limited by its capacity

  • Power flow constraint: the flow in each line is calculated using the PTDF matrix and the net injections between generation and load in each node

  • Flow capacity constraint: each flow is limited by the line capacity.

(10)#\[\begin{align} \min \quad & \sum_{i\in[G]} c_i(p_i) \\ \text{s.t} \quad &(\lambda):\qquad \sum_{i\in[G]} p_i = \sum_{i\in[N]} D_i \\ &(\overline{\sigma}_i):\qquad p_i \le \overline{P}_i && \forall i\in[G] \\ &(\beta_{l}):\qquad f_l = \sum_{i\in[N]} B_{l,i}\Big(\sum_{j\in[G]_i} p_j - D_i\Big)&& \forall l\in[L] \\ &(\mu_l^+):\qquad f_l \le \overline{F}_l && \forall l\in[L] \\ &(\mu_l^-):\qquad -f_l \le \overline{F}_l && \forall l\in[L] \end{align}\]
# define model
m = gp.Model()
m.setParam("OutputFlag", 0)

# variables
p = m.addVars(n_gen, lb=0, ub=GRB.INFINITY, name="p") # generation at bus
f = m.addVars(n_line, lb=-GRB.INFINITY, ub=GRB.INFINITY, name="f") # flow through line 

# energy balance
m.addConstr(sum(p[i] for i in range(n_gen)) == sum(D[i] for i in range(n_bus)), name="enerbal")

# generator constraints
m.addConstrs(p[i] <= pmax[i] for i in range(n_gen))

# power flow 
m.addConstrs((f[l] == sum(B[l,i]*(p[i] - D[i]) for i in range(n_bus)) for l in range(n_line)), name="flowdef")
m.addConstrs((f[l] <= fmax[l] for l in range(n_line)), name="flowlim_pos")
m.addConstrs((-f[l] <= fmax[l] for l in range(n_line)), name="flowlim_neg")

# objective
m.setObjective(sum(p[i]*c[i] for i in range(n_gen)), GRB.MINIMIZE)

# run
m.optimize()

Set parameter Username

Academic license - for non-commercial use only - expires 2025-11-15

# objective (cost of production)
print("System cost")
print(f"{m.ObjVal:.0f}")

System cost
4450

# production 
p_res = [m.getVarByName(f"p[{i}]").X for i in range(n_gen)]
print("Production:")
print(p_res)

Production:
[90.0, 0.0, 40.0]

# flow
f_res = [m.getVarByName(f"f[{l}]").X for l in range(n_line)]
print("Flows:")
print(f_res)

Flows:
[30.0, 60.0, 30.0]

# congestion prices
mu_pos = [m.getConstrByName(f"flowlim_pos[{l}]").Pi for l in range(n_line)]
mu_neg = [m.getConstrByName(f"flowlim_neg[{l}]").Pi for l in range(n_line)]
print("mu+")
print(mu_pos)
print()
print("mu-")
print(mu_neg)

mu+
[0.0, 0.0, -285.0]

mu-
[0.0, 0.0, 0.0]

# lmp 
lam = m.getConstrByName("enerbal").Pi
lmp = [lam + sum(B[l,i]*(mu_pos[l] - mu_neg[l]) for l in range(n_line)) for i in range(n_bus)]
print("lambda")
print(lam)
print()
print("LMPs:")
print(lmp)

lambda
100.0

LMPs:
[5.0, -90.0, 100.0]

4. Define mathematical model using voltage angle formulation

Objective function:

  • Minimize the cost of generation given by cost function \(c_i(p_i)\).

Decision variables:

  • \(p_{i}^{\text{c}}\) production of generator i.

  • \(f_{ij}\) flow from node i to j.

  • \(\theta_{i}\) voltage angle in node i.

Constraints:

  • Energy balance: the total sum of generation subtracted by the load equals the sum of incoming flows into the node.

  • Capacity constraints: each generator is limited by its capacity.

  • Power flow constraint: the flow in each line is the susceptance multiplied my the difference in voltage angles from node i to j.

  • Flow capacity constraint: each flow is limited by the line capacity.

(11)#\[\begin{align} \min \quad & \sum_{i\in[G]} c_i(p_i) \\ \text{s.t} \quad &(\lambda_i):\qquad \sum_{j\in[G]_i}p_j - D_i = \sum_{j:ij\in\set{L}}f_{ij} && \forall i \in [N] \\ &(\overline{\sigma}_i):\qquad p_i \le \overline{P}_i && \forall i\in[G] \\ &(\beta_{ij}):\qquad f_{ij} = b_{ij}(\theta_i - \theta_j )&& \forall ij\in\set{L} \\ &\hspace{5em} \theta_{i=slack} = 0 \\ &(\mu_{ij}^+):\qquad f_{ij} \le \overline{F}_{ij} && \forall ij\in\set{L} \\ &(\mu_{ij}^-):\qquad -f_{ij} \le \overline{F}_{ij} && \forall ij\in\set{L} \end{align}\]

Define new data

Symbol

Description

Data

c

Generation cost

[5, 20,100]

\(\overline{P}_i\)

Capacity limit

[102, 100, 100]

\(\overline{F}_{i,j}\)

Flow capacity limit

\begin{bmatrix} 0 & 100 & 50 \ 50 & 0 & 100 \ 30 & 50 & 0 \end{bmatrix}

\(D_i\)

Load in node i

[0, 0, 130]

L

Number of Lines

3

G

Number of Generators

3

G

Number of Buses

3

b

susceptance

500

 
# define model 
# assuming generator 1 is in node 1, generator 2 is in node 2, generator 3 is in node 3
# assuming constant susceptance of 500
b = 500
slack = 1
fmax = np.array([[0, 100, 50], 
              [50,  0, 100],
              [30,  50, 0]])

m = gp.Model()
m.setParam("OutputFlag", 0)

# variables
p = m.addVars(n_gen, lb=0, ub=GRB.INFINITY, name="p") # generation at bus
f = m.addVars(n_bus, n_bus, lb=-GRB.INFINITY, ub=GRB.INFINITY, name="f") # flow through line 
theta = m.addVars(n_line, lb=-GRB.INFINITY, ub=GRB.INFINITY, name="f") # flow through line 

# energy balance
for i in range(n_bus):
    m.addConstr(p[i] - D[i] == sum(f[i,j] for j in range(n_bus)), name="enerbal")

# generator constraints
m.addConstrs(p[i] <= pmax[i] for i in range(n_gen))

# power flow 
for i in range(n_bus):
    for j in range(n_bus):
        m.addConstr((f[i,j] == b*(theta[i] - theta[j])), name="flowdef")

for i in range(n_bus): 
    for j in range(n_bus):
        m.addConstr((f[i,j] <= fmax[i,j]), name="flowlim_pos")
        m.addConstr((-f[i,j] <= fmax[i,j]), name="flowlim_neg")

# slack bus
m.addConstr(theta[1] == 0,name="slack_bus")


# objective
m.setObjective(sum(p[i]*c[i] for i in range(n_gen)), GRB.MINIMIZE)

# run
m.optimize()

# objective (cost of production)
print(f"{m.ObjVal:.0f}")

6450