eleNet20220927.py

import gurobipy as gp
from gurobipy import GRB
import numpy as np
import sys
from matplotlib import pyplot as plt
Vbase,Sbase,zbase = 12660,1e7,16.02756
lft = np.array([[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15],[15,16],[16,17],[17,18],[2,19],[19,20],[20,21],[21,22],[3,23],[23,24],[24,25],[6,26],[26,27],[27,28],[28,29],[29,30],[30,31],[31,32],[32,33]])
r32 = 1/zbase*np.array([0.0922,0.4930,0.3661,0.3811,0.8190,0.1872,0.7115,1.0299,1.0440,0.1967,0.3744,1.4680,0.5416,0.5909,0.7462,1.2889,0.7320,0.1640,1.5042,0.4095,0.7089,0.4512,0.8980,0.8959,0.2031,0.2842,1.0589,0.8043,0.5074,0.9745,0.3105,0.3411])
x32 = 1/zbase*np.array([0.0470,0.2512,0.1864,0.1941,0.7070,0.6188,0.2351,0.7400,0.7400,0.0651,0.1298,1.1549,0.7129,0.5260,0.5449,1.7210,0.5739,0.1565,1.3555,0.4784,0.9373,0.3084,0.7091,0.7071,0.1034,0.1447,0.9338,0.7006,0.2585,0.9629,0.3619,0.5302])
P33 = 1e6*np.array([0,100,90,120,60,60,200,200,60,60,45,60,60,120,60,60,60,90,90,90,90,90,90,420,420,60,60,60,120,200,150,210,60])/Sbase/150 # 33 nodes
Q33 = 1e6*np.array([0,60,40,80,30,20,100,100,20,20,30,35,35,80,10,20,20,40,40,40,40,40,50,200,200,25,25,20,70,100,70,100,40])/Sbase/150

Pb = np.zeros_like(x32) # branch active power flow
Qb = np.zeros_like(x32) # branch reactive power flow
Pb[5:17] = P[17]
Qb[5:17] = Q[17]
Pb[24:32] = P[32]
Qb[24:32] = Q[32]
Pb[2:5] = P[17] + P[32]
Qb[2:5] = Q[17] + Q[32]
Pb[21:25] = P[24]
Qb[21:25] = Q[24]
Pb[1] = P[17] + P[32] + P[24]
Qb[1] = Q[17] + Q[32] + Q[24]
Pb[17:21] = P[21]
Qb[17:21] = Q[21]
Pb[0] = P[17] + P[32] + P[24] + P[21]
Qb[0] = Q[17] + Q[32] + Q[24] + Q[21]

m = gp.Model('vols')
V = m.addVars(33)
m.addConstr(V[0]==1.05)
for b in range(32):
    m.addConstr( V[lft[b,1]-1] == V[lft[b,0]-1] - (Pb[b]*r32[b] + Qb[b]*x32[b]) )
m.setObjective(0)
m.optimize()
if m.status != GRB.OPTIMAL:
    print('opt Fail >>>>>>>>>>>>>')
    sys.exit(3)
# for i in range(33):
#     print('%8d|%12g' % (i,V[i].X))
# for i in [17,21,24,32]:
#     print(V[i].X)


print(Pb[0],Qb[0])