This notebook demonstrates the use of semidefinite relaxation techniques for AC optimal power flow (ACOPF) problems. A description of the problem formulation and its implementation can be found in the following papers:
M. S. Andersen, A. Hansson, L. Vandenberghe, "Reduced-Complexity Semidefinite Relaxations of Optimal Power Flow Problems", IEEE Transactions on Power Systems, 29 (4), pp. 1855–1863, 2014.
A. Eltved, J. Dahl, M. S. Andersen, "On the Robustness and Scalability of Semidefinite Relaxation for Optimal Power Flow Problems", arXiv (math.OC), 2018.
To run this notebook, you will need CVXOPT, CHOMPACK, MOSEK, and Matplotlib.
MATPOWER provides a number of test cases that are basically MATLAB m-files, most of which (those that do not perform any computations) can be loaded directly from a file or from a url. For example, the following code
from opfsdr import opf
prob = opf('https://raw.githubusercontent.com/MATPOWER/matpower/master/data/case300.m')downloads the case300 test case from the MATPOWER repository on Github and generates a semidefinite relaxation problem. Alternatively, using the Github API, a list of test cases and their urls can be retrieved from the MATPOWER repository as follows:
In [1]:
import json, re
import requests
testcases = {}
clist = []
# Retrieve list of MATPOWER test cases
response = requests.get('https://api.github.com/repos/MATPOWER/matpower/contents/data')
clist += json.loads(response.text)
# Retrieve list of pglib-opf test cases
response = requests.get('https://api.github.com/repos/power-grid-lib/pglib-opf/contents/')
clist += json.loads(response.text)
response = requests.get('https://api.github.com/repos/power-grid-lib/pglib-opf/contents/api')
clist += json.loads(response.text)
response = requests.get('https://api.github.com/repos/power-grid-lib/pglib-opf/contents/sad')
clist += json.loads(response.text)
# Build list of test cases
for c in clist:
if not c['name'].endswith('.m'): continue
casename = c['name'].split('.')[0]
testcases[casename] = c['download_url']
In [2]:
from opfsdr import opf
#case = 'pglib_opf_case1354_pegase__sad'
case = 'case1354pegase'
print("Test case: %s" % case)
options = {'conversion': True, # apply chordal conversion? (default: False)
'tfill': 16, # use clique merging heuristic based on fill (default: 0)
'tsize': 0, # use clique merging heuristic based on clique size (default: 0)
'branch_rmin': 1e-5, # minimum transmission line resistance (default: -inf)
'line_constraints': True, # include apparent power flow constraints in problem formulation
'pad_constraints': True, # include phase angle difference constraints
'truncate_gen_bounds': 1e3, # reduce large generator bounds (default: None)
'verbose': 1, # print info, progress, etc. (default: 0)
'scale': False # apply scaling heuristic to cone LP? (default: False)
}
%time prob = opf(testcases[case], **options)
print(prob)
The semidefinite relaxation can be solved with MOSEK as follows:
In [3]:
# Set MOSEK tolerances
from opfsdr import msk
msk.options[msk.mosek.dparam.intpnt_co_tol_pfeas] = 1e-9 # default: 1e-8
msk.options[msk.mosek.dparam.intpnt_co_tol_dfeas] = 1e-9 # default: 1e-8
msk.options[msk.mosek.dparam.intpnt_co_tol_rel_gap] = 1e-8 # default: 1e-7
In [4]:
# Solve semidefinite relaxation
%time sol = prob.solve(solver="mosek")
In [5]:
print("Solver status: %s" % sol['status'])
The objective value associated with the semidefinite relaxation does not include constant cost, and hence it is not necessarily equal to the cost of generation (the cost of generation is available as sol['cost']):
In [6]:
print("Generation cost: %.2f USD/hour (fixed cost: %.2f USD/hour)"%(sol['cost'],prob.const_cost))
In [7]:
%pylab inline
from pylab import plot,xlabel,ylabel,arange,linspace,hist,xscale,rcParams
rcParams['figure.figsize'] = (10, 6)
if len(sol['eigratio']) == 1:
print("Eigenvalue ratio: %.2e" % (sol['eigratio'][0]))
else:
import numpy as np
hist(sol['eigratio'],bins=[10.**a for a in linspace(0,9,100)])
xscale('log')
xlabel('Clique eigenvalue ratio')
ylabel('Frequency')
print("Cliques with eigenvalue ratio less than 1e5: %i of %i"\
% (len([1 for evr in sol['eigratio'] if evr < 1e5]),len(sol['eigratio'])))
In [8]:
pmin = np.array([gen['Pmin'] for gen in prob.generators])*prob.baseMVA
pmax = np.array([gen['Pmax'] for gen in prob.generators])*prob.baseMVA
hist((np.array(sol['Sg'].real().T).squeeze()-pmin)/(pmax-pmin),20)
xlabel(r'$\frac{P_{\mathrm{g}}-P_{\mathrm{g}}^{\min}}{P_{\mathrm{g}}^{\max} - P_{\mathrm{g}}^{\min}}$')
ylabel('Number of generators');
In [9]:
pmin = np.array([gen['Qmin'] for gen in prob.generators])*prob.baseMVA
pmax = np.array([gen['Qmax'] for gen in prob.generators])*prob.baseMVA
hist((np.array(sol['Sg'].imag().T).squeeze()-pmin)/(pmax-pmin),20)
xlabel(r'$\frac{Q_{\mathrm{g}}-Q_{\mathrm{g}}^{\min}}{Q_{\mathrm{g}}^{\max} - Q_{\mathrm{g}}^{\min}}$')
ylabel('Number of generators');
In [10]:
minvm = np.array([b['minVm'] for b in prob.busses])
maxvm = np.array([b['maxVm'] for b in prob.busses])
hist((np.array(sol['Vm']).squeeze()-minvm)/(maxvm-minvm),20)
xlabel(r'$\frac{|V| - V_{\min}}{V_{\max}-V_{\min}}$')
ylabel('Number of busses');
In [11]:
if prob.branches_with_flow_constraints():
Sapp = np.array(abs(sol['St'])).squeeze()
Smax = np.array([br['rateA'] for _,br in prob.branches_with_flow_constraints()])
hist(Sapp/Smax,20)
pyplot.xlabel(r'$\frac{|S_{i,j}|}{S_{i,j}^{\max}}$')
pyplot.ylabel('Number of transmission lines')
else:
print("No flow constraints.")