I want to convert this python code to matlab code . How can i do it?

1 次查看（过去 30 天）
DHARTI PATEL
DHARTI PATEL 2021-3-15
评论： DHARTI PATEL 2021-5-4
from __future__ import division, print_function
from gurobipy import Model as GurobiModel, GRB, quicksum
from numpy import array, sqrt, real, imag, pi
from math import asin
from scipy.sparse import dok_matrix, hstack
from collections import defaultdict, deque
from loadcase import load_case
def build_U_matrices(G, B):
S2 = sqrt(2)
n = G.shape[0]
Ureal = dok_matrix((n-1, n))
Ureac = dok_matrix((n-1, n))
S2 = 2**.5
for i in range(1, n):
Ureal[i-1, i] = S2*G[i, :].sum()
Ureac[i-1, i] = -S2*B[i, :].sum()
return Ureal, Ureac
def build_R_matrices(G, B, branch_map):
""" rows are buses 2 to n; cols are branches """
n = G.shape[0]
Rreal = dok_matrix((n-1, n-1))
Rreac = dok_matrix((n-1, n-1))
for fbus in range(1, n):
for tbus in B[fbus, :].nonzero()[1]:
branch = branch_map[(fbus, tbus)]
Rreal[fbus-1, branch] = -G[fbus, tbus]
Rreac[fbus-1, branch] = B[fbus, tbus]
return Rreal, Rreac
def build_I_matrices(G, B, branch_map):
""" rows are buses 2 to n; cols are branches """
n = B.shape[0]
Ireal = dok_matrix((n-1, n-1))
Ireac = dok_matrix((n-1, n-1))
for fbus in range(1, n):
for tbus in B[fbus, :].nonzero()[1]:
branch = branch_map[(fbus, tbus)]
s = -1 if fbus < tbus else 1
Ireal[fbus-1, branch] = s*B[fbus, tbus]
Ireac[fbus-1, branch] = s*G[fbus, tbus]
return Ireal, Ireac
def build_constraint_matrix(G, B, branch_map):
Ureal, Ureac = build_U_matrices(G, B)
Rreal, Rreac = build_R_matrices(G, B, branch_map)
Ireal, Ireac = build_I_matrices(G, B, branch_map)
Areal = hstack([Ureal, Rreal, Ireal])
Areac = hstack([Ureac, Rreac, Ireac])
return Areal, Areac
def build_gurobi_model(case):
G, B = case.G, case.B
P = real(case.demands)
Q = imag(case.demands)
branches = case.branch_list
n = len(case.demands)
vhat = case.vhat
s2 = 2**.5
gens = {bus: gen.v for bus, gen in case.gens.items()}
del gens[0]
m = GurobiModel("jabr")
u = [m.addVar(name='u_%d'%i) for i in range(n)]
R = {(i, j): m.addVar(name='R_%d_%d' % (i, j)) for i, j in branches}
I = {(i, j): m.addVar(lb=-GRB.INFINITY, name='I_%d_%d' % (i, j)) for i, j in branches}
for i, j in branches:
R[j, i] = R[i, j]
I[j, i] = I[i, j]
m.update()
m.addConstr(u[0] == vhat*vhat/s2, 'u0')
for gen, v in gens.iteritems():
m.addConstr(u[gen] == v*v/s2, 'u%d' % gen)
for i, j in branches:
m.addQConstr(2*u[i]*u[j] >= R[i,j]*R[i,j] + I[i,j]*I[i,j], 'cone_%d_%d' % (i, j))
k = lambda i: (j for j in B[i, :].nonzero()[1])
s = lambda i, j: 1 if i < j else -1
for i in range(1, n):
m.addConstr(-s2*u[i]*G[i, :].sum() + quicksum(G[i,j]*R[i,j] + B[i,j]*s(i,j)*I[i,j] for j in k(i)) == P[i],
'real_flow_%d_%d' % (i, j))
if i in gens:
continue
m.addConstr(s2*u[i]*B[i, :].sum() + quicksum(-B[i,j]*R[i,j] + G[i,j]*s(i,j)*I[i,j] for j in k(i)) == Q[i],
'reac_flow_%d_%d' % (i, j))
m.setObjective(quicksum(R[i,j] for i, j in branches), sense=GRB.MAXIMIZE)
m.params.outputFlag = 0
#m.params.barQCPConvTol = 5e-10
m.optimize()
if m.status != 2:
raise ValueError("gurobi failed to converge: %s (check log)" % m.status)
u_opt = [x.getAttr('x') for x in u]
R_opt = {(i, j): x.getAttr('x') for (i, j), x in R.items()}
I_opt = {(i, j): x.getAttr('x') for (i, j), x in I.items()}
return u_opt, R_opt, I_opt
def recover_original_variables(u, I):
""" given Jabr variables u and I, return bus voltages and angles """
V = sqrt(sqrt(2) * array(u))
theta_branch = {(i, j): asin(I[(i, j)]/(V[i]*V[j])) for (i, j) in I}
theta_bus = recover_bus_angles(theta_branch)
return V, theta_bus
def recover_bus_angles(theta_branch):
""" assumes 0 is the root. traverses the tree and converts branch voltage
angles to bus voltage angles """
# TODO clarify what you're doing here
A = defaultdict(list) # adjacency matrix
theta_bus_dict = {0: 0}
for i, j in theta_branch.keys():
A[i] += [j]
q = deque([0])
while q:
parent = q.popleft()
for child in A[parent]:
if child not in theta_bus_dict:
theta_ij = theta_branch[(parent, child)]
if parent > child:
theta_ij *= -1
theta_bus_dict[child] = theta_bus_dict[parent] - theta_ij
q.append(child)
theta_bus = [theta_bus_dict[i] for i in range(max(theta_bus_dict.keys())+1)]
return theta_bus
def solve(casefile):
""" given a matpower casefile, solves the power flow using the Jabr method
and returns a dictionary mapping bus number to
(voltage magnitude, voltage angle) tuples. angles are in radians
"""
case = load_case(casefile)
u, R, I = build_gurobi_model(case)
V, theta = recover_original_variables(u, I)
i2e = case.i2e
answer = {i2e[i]: (v, t) for (i, (v, t)) in enumerate(zip(V, theta))}
return answer
if __name__ == '__main__':
answer = solve('cases/case5_renumber_tree.m')
for bus in sorted(answer.keys()):
v, t = answer[bus]
print('%3d %7.3f %7.3f' % (bus, v, 180/pi*t))

请先登录，再回答此问题。

回答（1 个）

Monika Jaskolka
Monika Jaskolka 2021-5-4

类别

在 Help Center 和 File Exchange 中查找有关 Data Type Identification 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by