diff --git a/opf_losses_v2.py b/opf_losses_v2.py index 9d55747..1ba2b35 100644 --- a/opf_losses_v2.py +++ b/opf_losses_v2.py @@ -1,21 +1,15 @@ -### this code is for OPF with voltage graph import itertools - -import numpy as np import pyomo.environ as py -import matplotlib as plt -import matplotlib.pyplot as plt import pandas as pd import pylab as pl -# to ignore the warning for dividing by 0 - -import matplotlib.pyplot as plt import numpy as np from matplotlib import pyplot as plt -from matplotlib import cm -from pyomo.core import Objective + + +# to ignore the warning for dividing by 0 np.seterr(divide='ignore', invalid='ignore') + ############################################################First Step############################################# # impedance matrix # importing the data from the excel for the bus @@ -46,7 +40,6 @@ for y in range(len(dataimpedence)): print(shunt_matrix) # converting the impedance to admittance - yb = -1 / Imp_matrix yb[np.isnan(yb)] = 0 @@ -62,8 +55,7 @@ y_offdiagonal = pyb y_offdiagonal_sum = np.sum(y_offdiagonal, axis=1) # print(y_offdiagonal_sum) -# now i will add the y_offdiagonal_sum with shunt_matrix_sum for diagonal matrix - +# now adding the y_offdiagonal_sum with shunt_matrix_sum for diagonal matrix yii = np.add(shunt_matrix_sum, y_offdiagonal_sum) # get the diaginal matrix # print(yii) @@ -73,7 +65,6 @@ ybus = yb print('Y bus: ', ybus) # calculating the magnitude and angle of ybus matrix - y_mag = abs(ybus) print('Y bus magnitude : ', y_mag) y_rad = np.angle(ybus) @@ -92,7 +83,6 @@ print('Susseptance: ', ybus_imag) # optimal power flow - model = py.ConcreteModel() # Set # model.b = py.Set(initialize=[0, 1, 2, 3, 4], doc='bus') @@ -112,20 +102,7 @@ model.qdmand = py.Param(model.b, model.t, initialize=dict( zip(list(itertools.product(range(1, len(databus.bus.values) + 1), range(1, len(model.t.data()) + 1))), np.append(datademand.reactive_power.values, len(model.t.data())))), doc='Reactive power demand') -# fixed demand -# model.pdmand = py.Param(model.b, model.t, initialize=dict( -# zip(list(itertools.product(range(1,len(databus.bus.values)+1), range(1,len(model.t.data())+1))), -# np.repeat(datademand.active_power1.values, len(model.t.data())))),mutable=True, -# doc='Active power demand') -# model.qdmand = py.Param(model.b, model.t, initialize=dict( -# zip(list(itertools.product(range(1,len(databus.bus.values)+1), range(1,len(model.t.data())+1))), -# np.repeat(datademand.reactive_power1.values, len(model.t.data())))),mutable=True, -# doc='Reactive power demand') - - # Variable - -# todo do we need to put some bound on the voltage and angle need to fixed it model.v = py.Var(model.index_i, model.t, initialize=1, bounds=(0.9, 1.07)) model.v[1, :].fix(1.06) model.v[2, :].fix(1.045) @@ -134,8 +111,6 @@ model.theta = py.Var(model.index_i, model.t, initialize=0, bounds=(-np.pi, np.pi model.theta[1, :].fix(0) # bounds for generator active power - - model.plb = py.Param(model.b, initialize=dict(zip(databus.bus, datagen.pl)), doc='active power lower bound') model.pub = py.Param(model.b, initialize=dict(zip(databus.bus, datagen.pu)), @@ -150,9 +125,8 @@ def pgenbound(model, i, t): model.pgen = py.Var(model.b, model.t,bounds=(0, None)) model.pgen[4, :].fix(0) model.pgen[5, :].fix(0) -# bounds for generator reactive power -# Todo should i need to put j with each element +# bounds for generator reactive power model.qlb = py.Param(model.b, initialize=dict(zip(databus.bus, datagen.ql)), doc='reactive power lower bound') model.qub = py.Param(model.b, initialize=dict(zip(databus.bus, datagen.qu)), doc='reactive power upper bound') @@ -169,17 +143,12 @@ print('model data is ok') # Constraints -# Todo check the error for indexing def active_power(model, i, t): return model.pgen[i, t] - model.pdmand[i, t]- sum( [y_mag[i - 1, j - 1] * (model.v[i, t]) * (model.v[j, t]) * py.cos( y_rad[i - 1, j - 1] + model.theta[j, t] - model.theta[i, t]) for j in model.index_j]) ==0 - -# i is not equal to j - - model.Const_Active_power = py.Constraint(model.index_i, model.t, rule=active_power, doc='active power injection') @@ -190,24 +159,17 @@ def reactive_power(model, i, t): y_rad[i - 1, j - 1] + model.theta[j, t] - model.theta[i, t]) for j in model.index_j])==0 - model.Const_Reactive_power = py.Constraint(model.index_i, model.t, rule=reactive_power, doc='Reactive power injection') print('model constraint are ok') -# Objective Function * .25 - - for k in model.index_i: for l in model.index_j: ybus_real[k - 1, l - 1] = (-(ybus_real[k - 1, l - 1])) def objective_rule(model): - # for m in model.index_i: - # for n in model.index_j: - # if m != n: return sum( sum(0.5 * (ybus_real[m - 1, n - 1]) * ((model.v[m, l]) ** 2 + (model.v[n, l]) ** 2 - 2 * ( (model.v[m, l]) * (model.v[n, l])) * py.cos(model.theta[n, l] - model.theta[m, l])) for m in @@ -216,7 +178,6 @@ def objective_rule(model): model.objective = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function') -print(' model objective is ok ') # Solver opt = py.SolverFactory('ipopt') @@ -224,7 +185,6 @@ result = opt.solve(model) model.write('model_opf.nl', io_options={'symbolic_solver_labels': True}) print('objective function = ', py.value(model.objective) * 100) -# model.display() if (result.solver.status == py.SolverStatus.ok) and ( result.solver.termination_condition == py.TerminationCondition.optimal): @@ -236,92 +196,6 @@ else: print('Solver Status: ', result.solver.status) # ploting: - -# pl_pgen = [] -# pl_pdem = [] -# pl_qdem = [] -# -# pl_qgen = [] -# pl_v = [[] for j in model.b] -# pl_theta = [[] for j in model.b] -# -# for ll in range( -# len([model.pgen, model.qgen, model.v, model.theta, model.pdmand, model.qdmand])): -# if ll == 0: -# for i in model.pgen: -# pl_pgen.append(py.value(model.pgen[i])) -# if ll == 1: -# for i in model.qgen: -# pl_qgen.append(py.value(model.qgen[i])) -# -# # if ll == 2: -# # for i in model.v: -# # pl_v.append(py.value(model.v[i])) -# # if ll == 3: -# # for i in model.theta: -# # pl_theta.append(py.value(model.theta[i])) -# if ll == 3: -# for i in model.pdmand: -# pl_pdem.append(py.value(model.pdmand[i])) -# if ll == 3: -# for i in model.qdmand: -# pl_qdem.append(py.value(model.qdmand[i])) -# for j in model.b: -# for k, v in model.v.items(): -# if k[1] == j: -# pl_v[j -1].append(py.value(v)) -# for j in model.b: -# for k, v in model.theta.items(): -# if k[1] == j: -# pl_theta[j-1 ].append(py.value(v)) -# # width = 0.2 -# b = np.arange(len(model.b)) -# t = np.arange(len(model.t)) -# x = b # number of buss -# y = t # time -# z = pl_v # votlage -# X, Y = np.meshgrid(y, x) -# Z = np.reshape(z, X.shape) # Z.shape must be equal to X.shape = Y.shape -# fig = plt.figure() -# ax = fig.gca(projection='3d') -# ax.plot_surface(X, Y, Z, cmap=cm.coolwarm) -# ax.set_xlabel('time') -# ax.set_ylabel('bus') -# ax.set_zlabel('voltage') - -# plt.show() -# fig = plt.figure() -# ax = fig.add_subplot(111, projection='3d') -# ax.plot_surface(X, Y, Z) -# ax.set_xlabel('time') -# ax.set_ylabel('bus') -# ax.set_zlabel('voltage') -# plt.show() - -# ax1.set_xlabel('x axis') -# ax1.set_ylabel('y axis') -# ax1.set_zlabel('z axis') -# fig2, ax = plt.subplots() -# v = ax.bar(x1, pl_v, label='voltage-bus_i') -# ax.set_ylabel('Voltages') -# ax.set_title('Bus') -# ax.set_xticks(x) -# ax.set_xticklabels(x) -# ax.bar_label(v) -# pl.ylabel('Voltage (PU)') -# pl.title('Voltages of the buses') -# fig2.tight_layout() -# pl.show() - -# fig2, ax = plt.subplots(2, len(b), constrained_layout=True, figsize=(10, 10)) -# for k in range(len(b)): -# ax[0, k].bar(t, pl_v[k], color="blue") -# ax[1, k].bar(t, pl_theta[k], color="blue") -# ax[0, 0].set_ylabel('Voltage') -# ax[0, k].set_title(f"Bus {k + 1}") -# fig2.suptitle('Grid Data', fontsize=10) -# pl.show() - pl_pdem = [[] for j in model.b] pl_pgen = [[] for j in model.b] time = [i for i in model.t] @@ -366,7 +240,9 @@ t = np.arange(len(model.t)) x = b # number of buss y = t # time z = pl_v # votlage -# + + +#3d voltage profile # X, Y = np.meshgrid(y, x) # Z = np.reshape(z, X.shape) # Z.shape must be equal to X.shape = Y.shape # fig = plt.figure()