Distribution_Network_Optimization/opf_losses_v2.py

### 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
# <editor-fold desc="Reading the data from Excel file">
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')
# </editor-fold>
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()
# ##
#