From 3e96b4222919e94d61aa658d9e01d8910712703b Mon Sep 17 00:00:00 2001 From: Sadam Hussain Date: Sat, 18 May 2024 11:53:01 -0400 Subject: [PATCH] Upload files to "/" --- opf_losses_v2.py | 405 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 405 insertions(+) create mode 100644 opf_losses_v2.py diff --git a/opf_losses_v2.py b/opf_losses_v2.py new file mode 100644 index 0000000..9d55747 --- /dev/null +++ b/opf_losses_v2.py @@ -0,0 +1,405 @@ +### 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 + +np.seterr(divide='ignore', invalid='ignore') +############################################################First Step############################################# +# impedance matrix +# importing the data from the excel for the bus +# +data = pd.ExcelFile('opf_data_time.xlsx') +databus = data.parse('bus') +datademand = data.parse('demand') +datacoeff = data.parse('coefficent') +datagen = data.parse('gen') +datatime = data.parse('time') +# +dataimpedence = data.parse('impedance') +Imp_matrix = np.zeros((5, 5)).astype(complex) +for y in range(len(dataimpedence)): + Imp_matrix[dataimpedence.loc[y, "From"] - 1, dataimpedence.loc[y, "To"] - 1] = dataimpedence.loc[y, "R"] + 1j * \ + dataimpedence.loc[y, "X"] + Imp_matrix[dataimpedence.loc[y, "To"] - 1, dataimpedence.loc[y, "From"] - 1] = Imp_matrix[ + dataimpedence.loc[y, "From"] - 1, dataimpedence.loc[y, "To"] - 1] + +print(Imp_matrix) + +# shunt admitance/ line charging +shunt_matrix = np.zeros((5, 5)).astype(complex) +for y in range(len(dataimpedence)): + shunt_matrix[dataimpedence.loc[y, "From"] - 1, dataimpedence.loc[y, "To"] - 1] = 1j * dataimpedence.loc[y, "B"] + shunt_matrix[dataimpedence.loc[y, "To"] - 1, dataimpedence.loc[y, "From"] - 1] = shunt_matrix[ + dataimpedence.loc[y, "From"] - 1, dataimpedence.loc[y, "To"] - 1] +print(shunt_matrix) + +# converting the impedance to admittance + +yb = -1 / Imp_matrix +yb[np.isnan(yb)] = 0 + +# positive of yb +pyb = -(yb) + +# sum the row of the shunt admittance to preparid for the diaginal element +shunt_matrix_sum = np.sum(shunt_matrix, axis=1) + +# ist i make the diagonal matrix to zero and then add row wise +pyb[np.diag_indices_from(pyb)] = 0 +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 + +yii = np.add(shunt_matrix_sum, y_offdiagonal_sum) # get the diaginal matrix +# print(yii) + + +np.fill_diagonal(yb, yii) +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) +print('Y bus angle (Rad) : ', y_rad) +y_ang = np.degrees(y_rad) + +# convert the ybus matrix into real (Conductance) and imaginary ( substance) matrix +ybus_real = ybus.real +print('Conductance: ', ybus_real) + +ybus_imag = ybus.imag + +print('Susseptance: ', ybus_imag) + +#########################################################Second Step################################################### + + +# optimal power flow + +model = py.ConcreteModel() +# Set +# model.b = py.Set(initialize=[0, 1, 2, 3, 4], doc='bus') +model.b = py.Set(initialize=databus.bus, doc='bus') +model.index_i = py.Set(initialize=databus.bus, doc='index i') +model.index_j = py.Set(initialize=databus.bus, doc='index j') +model.t = py.Set(initialize=datatime.time, doc='time') + +# Parameter + +# demand form the homes orginal +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.append(datademand.active_power.values, len(model.t.data())))), + 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.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) +model.v[3, :].fix(1.03) +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)), + doc='active power upper bound') + + +def pgenbound(model, i, t): + return model.plb[i], model.pub[i] + + +# model.pgen = py.Var(model.b, model.t, bounds=pgenbound) +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 + +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') + + +def qgenbound(model, i, t): + return model.qlb[i], model.qub[i] + + +# model.qgen = py.Var(model.b, model.t, bounds=qgenbound) +model.qgen = py.Var(model.b, model.t) +model.qgen[4, :].fix(0) +model.qgen[5, :].fix(0) +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') + + +# the plus sign is because of the P-Qj = ..... +def reactive_power(model, i, t): + return model.qgen[i, t] - model.qdmand[i, t]+ sum( + [y_mag[i - 1, j - 1] * (model.v[i, t]) * (model.v[j, t]) * py.sin( + 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 + model.index_i for n in model.index_j if m != n) for l in model.t) + + +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') +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): + print('optimal solution') # Do something when the solution in optimal and feasible +elif result.solver.termination_condition == py.TerminationCondition.infeasible: + print('solution is not feasible') # Do something when model in infeasible +else: + # Something else is wrong + 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] +pl_qdem = [[] for j in model.b] +pl_v = [[] for j in model.b] +pl_qgen = [[] for j in model.b] +pl_theta = [[] for j in model.b] + +for j in model.b: + for k, v in model.pgen.items(): + if k[0] == j: + pl_pgen[j - 1].append(py.value(v)) + +for j in model.b: + for k, v in model.qgen.items(): + if k[0] == j: + pl_qgen[j - 1].append(py.value(v)) + +for j in model.b: + for k, v in model.v.items(): + if k[0] == j: + pl_v[j - 1].append(py.value(v)) + +for j in model.b: + for k, v in model.theta.items(): + if k[0] == j: + pl_theta[j - 1].append(py.value(v)) + +for j in model.b: + for k, v in model.pdmand.items(): + if k[0] == j: + pl_pdem[j - 1].append(py.value(v)) + +for j in model.b: + for k, v in model.qdmand.items(): + if k[0] == j: + pl_qdem[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') +# # fig = plt.figure() +# ax = fig.add_subplot(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() + + + +# all graphs +fig2, ax = plt.subplots(5, len(model.b), constrained_layout=True, figsize=(10, 10)) +for k in range(len(model.b)): + ax[0, k].plot(t, pl_v[k], color="blue") + # ax[1, k].plot(t, pl_theta[k], color="blue") + ax[1, k].bar(t, pl_pdem[k], color="green") + ax[2, k].bar(t, pl_qdem[k], color="black") + ax[3, k].bar(t, pl_pgen[k], color="red") + ax[4, k].bar(t, pl_qgen[k], color="brown") + + ax[0, 0].set_ylabel('Voltage (PU)') + # ax[1, 0].set_ylabel('Angle') + ax[1, 0].set_ylabel('P_dem (PU)') + ax[2, 0].set_ylabel('Q_dem (PU)') + ax[3, 0].set_ylabel('P_gen (PU)') + ax[4, 0].set_ylabel('Q_gen (PU)') + + ax[0, k].set_title(f"Bus {k + 1}") + fig2.suptitle('Grid Data', fontsize=10) +pl.show() +# ## +#