centralized.py

#!/usr/bin/env python

import cvxpy as cp
from cvxpy.settings import OPTIMAL, OPTIMAL_INACCURATE
from math import sqrt
from typing import *


Vbase = 231
Sbase = 1e6
Zbase = Vbase**2 / Sbase

buses: List[dict] = [
    {'id':1, 'pc':0,    'pg':0,    'qc':0,        'sg':0},
    {'id':2, 'pc':5e03, 'pg':3e03,  'qc':5e03,      'sg':10e03},
    {'id':3, 'pc':4e03,  'pg':3e03,  'qc':6.113e03, 'sg':10e03},
    {'id':4, 'pc':10e03, 'pg':3e03,  'qc':0e03,      'sg':10e03},
]

lines = [
    {'r':0.09323, 'x':0.05585},
    {'r':0.19495, 'x':0.0812},
    {'r':0.11738, 'x':0.02925},
]


V0=231.0
epsilon=0.05

##############################################################################

V0 = V0/Vbase

r = [line['r']/Zbase for line in lines]
x = [line['x']/Zbase for line in lines]

Q = cp.Variable((len(buses)))
P = cp.Variable((len(buses)))
U = cp.Variable((len(buses)))

pc = [bus['pc']/Sbase for bus in buses]
qc = [bus['qc']/Sbase for bus in buses]

qg = [cp.Parameter(value=0), cp.Variable(), cp.Variable(), cp.Variable()]
pg = [bus['pg']/Sbase for bus in buses]

sg = [bus['sg']/Sbase for bus in buses]
s_h = [sqrt(sg[j]**2 - pg[j]**2) for j in range(len(buses))]

constraints = []

for j in range(len(buses)-1):
    jp1 = j+1
    constraints.extend((
        P[jp1] == (P[j] - pc[jp1] + pg[jp1]),
        Q[jp1] == (Q[j] - qc[jp1] + qg[jp1]),
        U[jp1] == (U[j] - 2 * (r[j] * P[j] + x[j] * Q[j])),
    ))
constraints.extend((
    P[len(buses)-1] == 0,
    Q[len(buses)-1] == 0,
    U[0] == 0,
))

for j in range(1, len(buses)):
    constraints.extend((
        V0**2 * (epsilon**2 - 2*epsilon) <= U[j],
                                            U[j] <= V0**2 * (epsilon**2 + 2*epsilon),
        cp.abs(qg[j]) <= s_h[j],
    ))

prob = cp.Problem(
    cp.Minimize(
        cp.sum(
            cp.multiply(
                r,
                (P[:-1]**2 + Q[:-1]**2) / V0**2
            )
        )
    ),
    constraints
)

prob.solve(verbose=True, solver=cp.ECOS)

if not prob.status in (OPTIMAL, OPTIMAL_INACCURATE):
    raise Exception(prob.status)

assert 6600 < qg[1].value * Sbase < 6900
assert 7600 < qg[2].value * Sbase < 7800
assert 2200 < qg[3].value * Sbase < 2500