Implementation of dual optimisation and solution recovery strategies. Download.
""" IIB Project: SDP relaxation dual optimisation and solution recovery method implementation. """
import sys
import pickle
import datetime
import time
import numpy as np
import scipy as sp
import cvxpy as cp
import matplotlib.pyplot as plt
sys.path.append('../') # add parent folder to Python path to allow import
from Datasets.Networks.test_case import network_structure, node_data # import network definition
from Datasets.Networks.haversine import haversine # custom implementation of Haversine function
# =============================================================================
# Grab power data and initialise required objects
# =============================================================================
# grab list of nodes
nodes = list(node_data.keys()) # this list defines the node indexing (i.e. the first is node 0, etc.)
N = len(nodes) # number of nodes in the system
# import node power data
fpath = '../Datasets/Analysis Datasets/test_case_2020-10-17.pkl'
with open(fpath, 'rb') as pfile: # read byte pickle file
power_data = pickle.load(pfile)
# create datetime objects of base samples for plotting
base_dt_objs = [datetime.datetime.strptime(d, '%Y-%m-%d %H:%M:%S') for d in power_data['timestamps']]
# specify initial fullness constants
alphas = 0.8*np.ones(N,dtype=float)
# define generation installed capacities and demand magnitude
wind_caps = {'lond':5, 'camb':10, 'bris':2.5, 'south':5, 'birm':0,
'liv':5, 'leed':0, 'newc':10, 'edin':10, 'inv':2.5} # wind installed capacities - GW
total_wind_cap = sum(wind_caps.values())
solar_caps = {'lond':2.5, 'camb':2.5, 'bris':3, 'south':5, 'birm':4,
'liv':1, 'leed':1, 'newc':1, 'edin':0, 'inv':0} # solar installed capacities - GW
total_solar_cap = sum(solar_caps.values())
assert [n for n in list(node_data.keys()) if node_data[n]['wind'] == [None]] == [n for n in list(wind_caps.keys()) if wind_caps[n] == 0]
peak_demand = 65 # peak national demand power - GW
# =============================================================================
# Produce node power deficit time series for analysis time window
# =============================================================================
start_time = '2018-03-01 09:00:00' # str time of start of analysis window
tstart_index = base_dt_objs.index(datetime.datetime.strptime(start_time, '%Y-%m-%d %H:%M:%S')) # index of window start
T = 6 # duration of window in time steps (of 1hr)
node_power_defs = np.zeros((N,T)) # initialise array to store node power deficit times series, (Pd - Pg)[t]
# time series are rows of the matrix, with the row index corresponding to the node number (row length = T)
# compute real power deficit series
for i in range(N):
n = nodes[i]
node_power_defs[i] += power_data[n]['demand'][tstart_index:tstart_index+T]*peak_demand
node_power_defs[i] -= power_data[n]['wind'][tstart_index:tstart_index+T]*wind_caps[n]
node_power_defs[i] -= power_data[n]['solar'][tstart_index:tstart_index+T]*solar_caps[n]
# if we are performing the analysis over a long duration, e.g. a year, at this stage we may
# wish to check the net energy deficit for the network over the time window
net_energy_deficit = node_power_defs.sum() # net energy deficit over time window - GWh
node_power_defs *= 1000 # convert powers to MW
# Note:
# The consistent numerical unit set to be used is (MVA,kV,kA,Ohms - hrs,MWh)
network_power_defs = np.sum(node_power_defs, axis=0) # net power deficit over network at each time instance - MWh
# =============================================================================
# Define parameters, matrices, and constraints required for the optimisation
# =============================================================================
Smax = 9e3 # line apparent power limit - MVA
Smax2 = Smax**2 # defined Smax-squared for convenience
# note, as some point we may wish to move to Smax being defined for each line
# for network node voltages, assume 400kV transmission level with \pm 10% operation tolerance
Vmin = 360 # minimum allowable node voltage - kV
Vmin2 = Vmin**2 # define Vmin-squared for convenience
Vmax = 440 # maximum allowable node voltage - kV
Vmax2 = Vmax**2 # define Vmax-squared for convenience
eta = 1 # round-trip efficiency of storage
dt = 1 # time increment - hrs
storage_flow_factor = np.sqrt(eta)/dt # define multiplier sqrt(eta)/dt
re_line_imp = 0.025 # real line impedance per unit length - Ohms/km
im_line_imp = 0.21 # imaginary line impedance per unit length - Ohms/km
pl_line_imp = re_line_imp + im_line_imp*1j # per unit length line impedance - Ohms/km
# construct the admittance matrix and set of lines
Y = np.zeros((N,N),dtype=complex) # initialise array for admittance matrix
lines = [] # initialise list to store connection node tuples
for k in range(N):
n = nodes[k] # base node
# note, in the model y_kk = 0, i.e. zero admittance-to-ground for the buses
for l in range(N):
m = nodes[l] # potentially connected node
if m in network_structure[n]['cons']: # m connected to n
lines.append((k,l)) # add connection to list of lines
# compute distance between the nodes - km
nm_dist = haversine(network_structure[n]['location'], network_structure[m]['location'])
z_nm = nm_dist*pl_line_imp # complex impedance of line n -> m - Ohms
y_nm = 1/z_nm # complex admittance of line n -> m - Siemens
Y[k,k] += y_nm # add term to Y[k,k] sum
Y[k,l] = -1*y_nm # set Y[k,l]
L = len(lines) # number of (directed) connections in the network
unique_lines = [p for p in lines if p[0]> 0]
# PSD constraints on A(lambda,R,t)
As = [] # initialise list to store A(lambda,R,t) matrices for later use
for t in range(T):
A = 0 # initialise A(lambda,R,t)
# construct A(lambda,R,t) matrix
for n in range(N): # for each node
A += lambda_P[n,t]*bYks[n] + (lambda_VU[n,t] - lambda_VL[n,t])*Mks[n]
for l,m in lines: # for each line
# note r^2 = R[0,1], r^3 = R[0,2]
A += 2*Rs[(l,m,t)][0,1]*bYlms[l,m] + 2*Rs[(l,m,t)][0,2]*bYbarlms[l,m]
constraints += [A >> 0]
As.append(A) # save A matrix for later observation
# KKT stationarity conditions/constraints for S_n
for n in range(N): # for each node
c = 1 - lambda_S[n] - alphas[n]*lambda_P[n,0]*storage_flow_factor
# note lambda_P is zero indexed in t, so this is indeed using
# lamdba_P(n,t=1) from the formulation expression
for t in range(T):
c += -1*lambda_eU[n,t]
constraints += [c == 0]
# KKT stationarity conditions/constraints for e_n[t]
for n in range(N): # for each node
for t in range(T-1): # for non-terminal time instances
c = storage_flow_factor*(lambda_P[n,t] - lambda_P[n,t+1]) + lambda_eU[n,t] - lambda_eL[n,t]
constraints += [c == 0]
# for the terminal time instance
c = storage_flow_factor*lambda_P[n,T-1] + lambda_eU[n,T-1] - lambda_eL[n,T-1]
constraints += [c == 0]
# epsilon-modification constraints - 2-norm chosen as vector norm to retain convexity in a simple manner
# epsilon = 1e-9 # modification parameter
# recepsilon = 1/epsilon
# for t in range(T): # for each time instance
# for n in range(N): # for each node
# lambda_comb = cp.vstack([lambda_P[n,t],lambda_VU[n,t],lambda_VL[n,t],lambda_eU[n,t],lambda_eL[n,t]])
# constraints += [cp.norm(lambda_comb,2) <= recepsilon]
# #constraints += [epsilon <= lambda_P[n,t], lambda_P[n,t] <= recepsilon]
# for l,m in lines: # for each line
# r_comb = cp.vec(Rs[(l,m,t)])
# constraints += [cp.norm(r_comb,2) <= recepsilon]
# Create dual problem objective function
# =============================================================================
obj = 0 # initialise objective
for t in range(T): # for each time instance
for n in range(N):
obj += lambda_P[n,t]*node_power_defs[n,t] + lambda_VL[n,t]*Vmin2 - lambda_VU[n,t]*Vmax2
for l,m in lines:
# note r^1 = R[0,0], r^4 = R[1,1], r^6 = R[2,2]
obj += -1*Rs[(l,m,t)][0,0]*Smax2 + -1*Rs[(l,m,t)][1,1] + -1*Rs[(l,m,t)][2,2]
# Define CVXPY problem
# =============================================================================
problem = cp.Problem(cp.Maximize(obj), constraints)
# Solve dual problem
# =============================================================================
print("\nOptimisation start.")
start = time.time()
# problem.solve(solver=cp.SCS,verbose=True,
# max_iters=int(5e4),eps=2e-4,
# use_indirect=True,alpha=1.0) # SCS solver
problem.solve(solver=cp.CVXOPT, verbose=True, kktsolver='robust',
refinement=5, reltol=1e-8, feastol=1e-7,
max_iters=250) # CVXOPT solver
# problem.solve(solver=cp.MOSEK, verbose=True,
# mosek_params={mosek.iparam.optimizer:mosek.optimizertype.conic,
# mosek.dparam.intpnt_co_tol_rel_gap:1e-12}) # MOSEK solver
end = time.time()
runtime = end-start
print("\n============================")
print("Op solve time: %ss %sms"%(int(runtime),int((runtime%1)*1000)))
print("============================\n")
# Verify conditions for exactness are satisfied
# =============================================================================
print("\nEvaluating exactness.")
failures = [] # initialise list to store faild time instances
rcond = 1e-4
# relative condition number - threshold for eigenvalue (really singular value)
# to be treated as zero - see https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.null_space.html
# hopefully all the eigenvalues are positive real, so they conincide with the singular values
for t in range(T): # for each time instance
eigvals = np.linalg.eigvals(As[t].value)
mev = np.max(np.abs(eigvals)) # find largest eigenvalue
numzeros = 0
for ev in eigvals:
if np.abs(ev) < mev*rcond:
numzeros += 1
if numzeros != 2: # criterion for exactness not satisfied
print("Fail; %s,%s"%(t,numzeros))
print(eigvals)
print("Constraint violation: %s"%constraints[L*T+t].violation())
failures.append(t) # store failed time index
def user_break(message):
answer = input(message)
if answer == 'y':
return # continue
elif answer == 'n':
sys.exit() # terminate program
else:
print("Invalid answer - choose [y]/n")
user_break(message) # ask again for a valid answer
if len(failures) != 0:
print("\n===================================")
print("Exactness condition failed at %s"%failures)
print("===================================")
user_break("Continue to solution recovery? [y]/n\n")
print("\n")
else:
print("\n==============================")
print("Exactness conditions satisfied")
print("==============================\n")
# =============================================================================
# Recover primal solution
# =============================================================================
Vopts = np.empty((T,N),dtype=np.complex) # initialise array to store optimised node voltage vectors
Vopts[:] = np.nan # set entries to nan to prevent confusion with unfilled time instances
Wopts = np.zeros((T,2*N,2*N)) # initialise array to store optimised W matrices
# Recover W matrices for time instances where exactness criterion is SATISFIED
# =============================================================================
# Method from Lavaei & Low 2012 / Gayme & Topcu 2012
print("Recovering voltage matrices for exact instances.")
for t in [h for h in range(T) if h not in failures]:
print(t)
ns_vecs = sp.linalg.null_space(As[t].value,rcond=rcond) # compute null space of A(lambda,R,t)
assert ns_vecs.shape == (2*N,2) # verify dimension of null space
# grab first null space basis vector and split into components
v1 = ns_vecs[:,0][:N] # first N entries
v2 = ns_vecs[:,0][N:] # last N entries
# identify the first node which has an exactly active constraint (i.e. the
# Lagrange multiplier is exactly zero), and set the boundary voltage
# variable accordingly
### for SCS which gives true zero values for lambda_VL/VU
# for k in range(N):
# if (lambda_VL.value[k,t] != 0) and (lambda_VU.value[k,t] == 0):
# Vbk2 = Vmin2
# break
# elif (lambda_VL.value[k,t] == 0) and (lambda_VU.value[k,t] != 0):
# Vbk2 = Vmax2
# break
### for CVXOPT which gives numerically zero values for lambda_VL/VU
VLU_tol = 1e-6
for k in range(N):
if (lambda_VL.value[k,t] > VLU_tol) and (lambda_VU.value[k,t] < VLU_tol):
Vbk2 = Vmin2
break
elif (lambda_VL.value[k,t] < VLU_tol) and (lambda_VU.value[k,t] > VLU_tol):
Vbk2 = Vmax2
break # note by breaking we retain the value of k at which this occurs
# calculate constants for constructing Vopt[t]
zeta1 = np.sqrt((Vbk2/(v1[k]**2 + v2[k]**2))/(1 + (v2[0]/v1[0])**2))
zeta2 = -1*(v2[0]/v1[0])*zeta1
# construct Vopt[t]
Vopt_t = (zeta1 + zeta2*1j)*(v1 + v2*1j)
# adjust sign of voltage vector to set phase angle of reference to 0 not 180
Vopt_t *= np.sign(Vopt_t[0])
Vopts[t] = Vopt_t # store constructed optimal voltage (row) vector
print(Vopt_t)
print(np.abs(Vopt_t))
# construct Wopt[t] and store
Xopt_t = np.hstack(( Vopt_t.real , Vopt_t.imag )) # extended voltage (row) vector
Xopt_t = np.reshape(Xopt_t, (2*N,1)) # convert to column vector
Wopt_t = Xopt_t@Xopt_t.T # produce Wopt[t] matrix
Wopts[t] = Wopt_t # store
print('\n')
# Recover W matrices for time instances which FAIL exactness criterion
# =============================================================================
print("Recovering voltage matrices for non-exact instances.")
redeemed_times = [] # initialise list to store time instances which satisfy exactness following recovery
for t in failures:
print(t)
# Set up the CVXPY problem for recovering Wopt
# =========================================================================
# create the optimisation variable W, specify symmetry
Wnamestr = "recW-%s"%t
W = cp.Variable(shape=(2*N,2*N),name=Wnamestr,symmetric=True)
rec_constraints = [] # initialise list of constraints for optimisation
# constraint W to be psd
rec_constraints += [W >> 0]
# apply KKT condition constraint - Tr(A(t)W[t]) = 0
rec_constraints += [cp.trace(As[t].value@W) == 0]
# apply node voltage constraints (primal)
for n in range(N):
rec_constraints += [Vmin2 <= cp.trace(Mks[n]@W)]
rec_constraints += [cp.trace(Mks[n]@W) <= Vmax2]
# apply line current constraints (primal)
for l,m in lines:
rec_constr_mat = cp.vstack((
cp.hstack(( -1*Smax2 , cp.trace(bYlms[l,m]@W) , cp.trace(bYbarlms[l,m]@W) )),
cp.hstack(( cp.trace(bYlms[l,m]@W), -1, 0 )),
cp.hstack(( cp.trace(bYbarlms[l,m]@W), 0 , -1 ))
))
rec_constraints += [rec_constr_mat << 0]
recovery_obj = cp.normNuc(W) # nuclear norm rank-regularisation objective
recovery_op = cp.Problem(cp.Minimize(recovery_obj), rec_constraints)
rec_max_iters = int(2e4) # solver max iterations for recovery
rec_tol = 1e-6 # solver tolerance for recovery
recovery_op.solve(solver=cp.SCS, verbose=False,
max_iters=rec_max_iters, eps=rec_tol,
use_indirect=True) # use solver SCS
# recovery_op.solve(solver=cp.CVXOPT, verbose=True, kktsolver='robust',
# refinement=5)
Weigvals = np.linalg.eigvals(W.value)
print(np.flip(sorted(Weigvals, key=lambda l:np.abs(l))))
print("PSD violation: %s"%rec_constraints[0].violation())
# note, the value of this violation must be considered relative to the eigenvalue size
Wrank = np.linalg.matrix_rank(W.value, tol=np.max(Weigvals)*rcond)
# note the same rcond as before is used for consistency
print("rank - %s"%Wrank)
if Wrank == 1:
print("Time %s redeemed!"%t)
redeemed_times.append(t)
# Store recovered W matrix
# =========================================================================
Wopts[t] = W.value
assert False not in [W.any() for W in Wopts] # ensure all W matrices calculated
### Check if power limit constraints are binding
# for t in range(T):
# for (l,m) in lines:
# vol = Smax2 - (np.trace(bYlms[l,m]@Wopts[t])**2 + np.trace(bYbarlms[l,m]@Wopts[t])**2)
# if vol < 1000:
# print(vol, (l,m), t)
# else:
# print(vol)
# Recover energy storage levels
# =============================================================================
print("\nRecovering storage energy levels")
"""
Work through the energy storage levels for each node. Start at
e_n[t] = alpha*S_n, and solve forwards in time in terms of S_n while
lambda^{eU/eL} = 0. When a non-zero value is reached, e_n[t] has a
determined value, so S_n can be calculated, allowing the previous values
to be back calculated. Solution forwards in time then proceeds from here.
"""
# initialise array to store store energy levels, e_n[t]
storage_energy_levels = np.zeros((N,T))
storage_capacities = np.zeros(N) # array to store storage capacities
### check Lagrange multipliers eL/eU for violations
# ztol = 1e-3
# for n in range(N):
# print(n)
# for t in range(T):
# fail = (lambda_eL.value[n,t] > ztol) and (lambda_eU.value[n,t] > ztol)
# if fail:
# print(t,lambda_eL.value[n,t],lambda_eU.value[n,t],"fail")
# else:
# print(t,lambda_eL.value[n,t],lambda_eU.value[n,t])
# print('\n')
assert not np.any(lambda_eL.value < 0) # check Lagrage multipliers non-negative
assert not np.any(lambda_eU.value < 0)
for n in range(N): # solve energy levels independently for each node
assert lambda_S.value[n] >= 0
lambda_S_tol = 1e-3 # tolerance for declaring lambda_S_n to be non-zero
if lambda_S.value[n] > lambda_S_tol: # if the storage Lagrange multiplier is non-zero
storage_capacities[n] = 0 # set storage capacity to zero
storage_energy_levels[n,:] = np.nan # declare energy levels undefined
else:
# set relatively small/effectively zero lambda_P values to zero
lambda_P_tol = 1e-5 # tolerance to zero declaration
lambdaPs = np.where(lambda_P.value[n] < lambda_P_tol, 0, lambda_P.value[n]) # screen lambda_Ps
# Compute energy levels of storage
# =====================================================================
nz_tol = 1e-3 # tolerance for lambda_eL/lambda_eU to be non-zero
### Method 1 - KKT, completementary slackness based method
### This is the method from Gayme & Topcu 2012
# Perform an initial pass to compute Sn
# =====================================================================
delta_e = 0 # initialise variable to store net energy accumulation
Sn_calculated = False # flag to show whether Sn have been calculated
for t in range(T):
if lambda_P.value[n,t] > nz_tol: # if lambda_P is non-zero
# increment energy accumulation
delta_e += -1*(1/storage_flow_factor)*(np.trace(bYks[n]@Wopts[t]) + node_power_defs[n,t])
print(t,delta_e)
if (lambda_eL.value[n,t] > nz_tol) and (lambda_eU.value[n,t] > nz_tol):
print("Warning! Invalid Lagrange multipliers encountered eL-eU before solution")
pass # ignore if both multipliers non-zero
# if lambda_eL > 0, then e_n[t] = 0
elif lambda_eL.value[n,t] > nz_tol:
S_n = -1*delta_e/alphas[n] # note delta_e should be negative
storage_capacities[n] = S_n # store determined storage capacity
print(S_n,n,t,1)
Sn_calculated = True # set flag
break # break to save computation
# if lambda_eU > 0, then e_n[t] = Sn
elif lambda_eU.value[n,t] > nz_tol:
if alphas[n] == 1:
pass # this level cannot be used to determine Sn
else:
S_n = delta_e/(1-alphas[n])
storage_capacities[n] = S_n # store determined storage capacity
Sn_calculated = True # set flag
print(S_n,n,t,2)
break
else:
pass
else: # if lambda_P is zero
if lambda_eU.value[n,t] > nz_tol:
S_n = delta_e/(1-alphas[n])
storage_capacities[n] = S_n # approximate storage capacity
Sn_calculated = True # set flag
print("Warning - Sn solution inaccurate, lambda_P = 0 encountered.")
break
else: # if lambda_P = 0 and e_n[t] not defined, we have a problem
raise ValueError("Indeterminate energy level")
### Is this is right way to deal with lambda_P = 0?
if not Sn_calculated:
raise ValueError("Storage capacity not calculated - no energy level constraints binding.")
### Method 2 - global consistency based method (more robust but less accurate)
# Perform an initial pass to compute Sn
# ====================================================================
# delta_es = np.zeros(T+1) # array to store net energy accumulation - with initial zero (t=0 instance)
# for t in range(T):
# if lambda_P.value[n,t] > nz_tol: # if lambda_P is non-zero
# delta_es[t+1] = delta_es[t] + -1*(1/storage_flow_factor)*(np.trace(bYks[n]@Wopts[t]) + node_power_defs[n,t])
# else: # if lambda_P is zero
# print("Warning - indeterminate energy level")
# delta_es[t+1] = delta_es[t] # approximate energy level
# #raise ValueError("Indeterminate energy level")
# if alphas[n] != 1:
# S_n = np.max([np.max(delta_es)/(1-alphas[n]),np.abs(np.min(delta_es))/alphas[n],np.ptp(delta_es)])
# # consider storage level change diagram
# else:
# if np.max(delta_es) <= 0:
# S_n = np.abs(np.min(delta_es)) # max energy extraction
# else:
# print("Warning - storage capacity violation")
# S_n = np.ptp(delta_es) # max energy extraction from excess peak
# storage_capacities[n] = S_n
# Perform a second pass to compute e_n[t]
# =====================================================================
for t in range(T):
if lambda_P.value[n,t] > nz_tol: # if lambda_P is non-zero
# compute storage levels using update rule
if t == 0:
storage_energy_levels[n,t] = alphas[n]*S_n + -1*(1/storage_flow_factor)*(np.trace(bYks[n]@Wopts[t]) + node_power_defs[n,t])
else:
storage_energy_levels[n,t] = storage_energy_levels[n,t-1] + -1*(1/storage_flow_factor)*(np.trace(bYks[n]@Wopts[t]) + node_power_defs[n,t])
else: # if lambda_P is zero
if lambda_eU.value[n,t] > nz_tol: # if storage capacity constraint active
storage_energy_levels[n,t] = S_n # set storage level to capacity
else: # if lambda_P = 0 and e_n[t] not defined, we have a problem
print("Warning - indeterminate energy level")
storage_energy_levels[n,t] = storage_energy_levels[n,t-1]
#raise ValueError("Indeterminate energy level")
print("\n\n===================")
print("Storage capacities:")
print(storage_capacities)
print("Total capacity: %sMWh - Obj fn value: %sMWh"%(int(np.sum(storage_capacities)),int(problem.value)))
print("%% error: %s%%"%(round((np.sum(storage_capacities)-problem.value)*100/problem.value,3)))
print("===================\n")
# Check energy level feasibility violations
for n in range(N):
if np.min(storage_energy_levels[n]) < 0:
print(n,"min neg",np.min(storage_energy_levels[n]))
if np.max(storage_energy_levels[n]) > storage_capacities[n]:
print(n,"max excess",np.max(storage_energy_levels[n]),storage_capacities[n])
# Plot energy storage levels
fig = plt.figure()
ax = plt.gca()
for n in range(N):
if storage_capacities[n] != 0:
plt.plot(range(T), storage_energy_levels[n]/1000,label=n+1)
plt.xlabel("Time instance, $t$")
plt.ylabel("Energy storage level (GWh)")
ax.grid(True, which='major', alpha=0.75, ls='--', zorder=0)
ax.set_axisbelow(True)
plt.legend(title="Node/Bus")
plt.xlim(0,T-1)
plt.tight_layout()
plt.savefig('solution.png', dpi=300)
plt.show()