Updating code of HEMS connected to transformer_model.py

This commit is contained in:
Sadam93 2024-05-31 11:40:34 -04:00
parent 69d363b01d
commit 1e8509714d
2 changed files with 3053 additions and 3050 deletions

File diff suppressed because it is too large Load Diff

View File

@ -1,228 +1,231 @@
m3 = py.ConcreteModel()
from HEMS import *
# tranformer parameters
m3.Transformer_Limit = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Transformer_limit)), m3 = py.ConcreteModel()
doc='Maximum power limit on transformer')
m3.tran_life = py.Param(mutable=True, doc='Transformer life expectancy') # tranformer parameters
m3.HST = py.Var(m1.t, bounds=(0, None), doc="windign hottest spot temp") m3.Transformer_Limit = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Transformer_limit)),
m3.LOL = py.Var(m1.t, bounds=(0, None), doc="lost of life") doc='Maximum power limit on transformer')
# m3.FAA = py.Var(m1.t, bounds=(0, None), doc="HST") m3.tran_life = py.Param(mutable=True, doc='Transformer life expectancy')
m3.DTOT = py.Var(m1.t, bounds=(0, None), doc="top oil rise temp") m3.HST = py.Var(m1.t, bounds=(0, None), doc="windign hottest spot temp")
m3.DHST = py.Var(m1.t, bounds=(0, None), doc="winding HST rise over the top") m3.LOL = py.Var(m1.t, bounds=(0, None), doc="lost of life")
m3.DTOTU = py.Var(m1.t, bounds=(0, None), doc="ultimate top-oil rise over the ambient temperature") # m3.FAA = py.Var(m1.t, bounds=(0, None), doc="HST")
m3.DHSTU = py.Var(m1.t, bounds=(0, None), doc="ultimate top-oil rise over the ambient") m3.DTOT = py.Var(m1.t, bounds=(0, None), doc="top oil rise temp")
m3.Power_P_T = py.Var(m1.t, bounds=(0, None), doc="power from home to Transformer") m3.DHST = py.Var(m1.t, bounds=(0, None), doc="winding HST rise over the top")
m3.Power_T_P = py.Var(m1.t, bounds=(0, None), doc="power from transformer") m3.DTOTU = py.Var(m1.t, bounds=(0, None), doc="ultimate top-oil rise over the ambient temperature")
m3.S_P_T = py.Var(m1.h, m1.t, within=py.Binary) m3.DHSTU = py.Var(m1.t, bounds=(0, None), doc="ultimate top-oil rise over the ambient")
m3.FAA = py.Var(m1.t, bounds=(0, None)) # The output of piecewise. Equivalent to model.y=piecewise(model.x) m3.Power_P_T = py.Var(m1.t, bounds=(0, None), doc="power from home to Transformer")
m3.n = py.Param(initialize=0.9, m3.Power_T_P = py.Var(m1.t, bounds=(0, None), doc="power from transformer")
doc='top oil rise constant for transformer ') # 1, 0.9=ONOF tranformer, 0.8=ONON Transformer, it was 1 m3.S_P_T = py.Var(m1.h, m1.t, within=py.Binary)
m3.m = py.Param(initialize=1.6, m3.FAA = py.Var(m1.t, bounds=(0, None)) # The output of piecewise. Equivalent to model.y=piecewise(model.x)
doc='winding constant for transformer') # 1, 0.8 , 0.8 so for 2m will 2, 1.6 , 1.6 it was 2 m3.n = py.Param(initialize=0.9,
# <editor-fold desc="transformer parameter"> doc='top oil rise constant for transformer ') # 1, 0.9=ONOF tranformer, 0.8=ONON Transformer, it was 1
m3.g1 = py.Param(initialize=0.1535, m3.m = py.Param(initialize=1.6,
doc='time duration ') # i t was 0.4866, but 0.1535 if we consider delta t to 15 min if we consider 0.25 detla t then 0.0027 doc='winding constant for transformer') # 1, 0.8 , 0.8 so for 2m will 2, 1.6 , 1.6 it was 2
# m3.g1 = py.Param(initialize=0.0027,doc='time duration ') # i t was 0.4866, if we consider 0.25 detla t then 0.0027 # <editor-fold desc="transformer parameter">
m3.g2 = py.Param(initialize=0.8826, m3.g1 = py.Param(initialize=0.1535,
doc='time duration ') # it was 0.9998 but 0.8826 if we consider delta t to 15 min if we consider 0.25 detla t then 0.038 doc='time duration ') # i t was 0.4866, but 0.1535 if we consider delta t to 15 min if we consider 0.25 detla t then 0.0027
# m3.g2 = py.Param(initialize=0.038, doc='time duration ') # it was 0.9998 but if we consider 0.25 detla t then 0.038 # m3.g1 = py.Param(initialize=0.0027,doc='time duration ') # i t was 0.4866, if we consider 0.25 detla t then 0.0027
m3.TR = py.Param(initialize=6, doc='Resistance ') m3.g2 = py.Param(initialize=0.8826,
m3.DTOR = py.Param(initialize=56, doc='e top-oil rise over ambient at the rated load ') # it was 56 doc='time duration ') # it was 0.9998 but 0.8826 if we consider delta t to 15 min if we consider 0.25 detla t then 0.038
m3.DHR = py.Param(initialize=80, doc='hottest-spot rise over top-oil at the rated load ') # 24 and 80 it was 24 # m3.g2 = py.Param(initialize=0.038, doc='time duration ') # it was 0.9998 but if we consider 0.25 detla t then 0.038
m3.HSTL = py.Param(initialize=110, doc='time duration ') m3.TR = py.Param(initialize=6, doc='Resistance ')
m3.t_TO = py.Param(initialize=180, doc=' in min ') # it was 90min in one paper and other paper 180 min m3.DTOR = py.Param(initialize=56, doc='e top-oil rise over ambient at the rated load ') # it was 56
m3.t_w = py.Param(initialize=48, doc=' in min ') # it was 7min in one paper and other paper 48 min m3.DHR = py.Param(initialize=80, doc='hottest-spot rise over top-oil at the rated load ') # 24 and 80 it was 24
m3.Ambient_Transformer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_tr.celsius1)), mutable=True, m3.HSTL = py.Param(initialize=110, doc='time duration ')
doc='Buying Price') m3.t_TO = py.Param(initialize=180, doc=' in min ') # it was 90min in one paper and other paper 180 min
m3.P_Buy_extra = py.Param(m1.h, m1.t, initialize=0) m3.t_w = py.Param(initialize=48, doc=' in min ') # it was 7min in one paper and other paper 48 min
m3.P_sell_extra = py.Param(m1.h, m1.t, initialize=0) m3.Ambient_Transformer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_tr.celsius1)), mutable=True,
delta_t = 15 doc='Buying Price')
m3.P_Buy_extra = py.Param(m1.h, m1.t, initialize=0)
m3.P_sell_extra = py.Param(m1.h, m1.t, initialize=0)
# <editor-fold desc="Transformer thermal Modeling"> delta_t = 15
def HST_A(m3, i):
return m3.HST[i] == m3.Ambient_Transformer[i] + m3.DTOT[i] + m3.DHST[i]
# <editor-fold desc="Transformer thermal Modeling">
def HST_A(m3, i):
m3.Const_Tr_a = py.Constraint(m1.t, rule=HST_A, doc='Hotest spot temperature') return m3.HST[i] == m3.Ambient_Transformer[i] + m3.DTOT[i] + m3.DHST[i]
def HST_B(m3, i): m3.Const_Tr_a = py.Constraint(m1.t, rule=HST_A, doc='Hotest spot temperature')
# return m3.DTOT[i] == m3.DTOTU[i] * m3.g1
return m3.DTOT[i] == m3.DTOTU[i] * (1 - py.exp(-delta_t / m3.t_TO.value))
def HST_B(m3, i):
# return m3.DTOT[i] == m3.DTOTU[i] * m3.g1
m3.Const_Tr_b = py.Constraint([m1.t_first.value], rule=HST_B, doc='Hotest spot temperature') return m3.DTOT[i] == m3.DTOTU[i] * (1 - py.exp(-delta_t / m3.t_TO.value))
def HST_C(m3, i): m3.Const_Tr_b = py.Constraint([m1.t_first.value], rule=HST_B, doc='Hotest spot temperature')
return m3.DTOT[i] == (m3.DTOTU[i] - m3.DTOT[i - 1]) * (1 - py.exp(-delta_t / m3.t_TO.value)) + m3.DTOT[i - 1]
# return m3.DTOT[i] == (m3.DTOTU[i] - m3.DTOT[i - 1]) * m3.g1 + m3.DTOT[i - 1]
def HST_C(m3, i):
return m3.DTOT[i] == (m3.DTOTU[i] - m3.DTOT[i - 1]) * (1 - py.exp(-delta_t / m3.t_TO.value)) + m3.DTOT[i - 1]
m3.Const_Tr_c = py.Constraint(m1.t_second, rule=HST_C, doc='Hotest spot temperature') # return m3.DTOT[i] == (m3.DTOTU[i] - m3.DTOT[i - 1]) * m3.g1 + m3.DTOT[i - 1]
def HST_D(m3, i): m3.Const_Tr_c = py.Constraint(m1.t_second, rule=HST_C, doc='Hotest spot temperature')
return m3.DHST[i] == m3.DHSTU[i] * (1 - py.exp(-delta_t / m3.t_w.value))
# return m3.DHST[i] == m3.DHSTU[i] * m3.g2.value
def HST_D(m3, i):
return m3.DHST[i] == m3.DHSTU[i] * (1 - py.exp(-delta_t / m3.t_w.value))
m3.Const_Tr_d = py.Constraint([m1.t_first.value], rule=HST_D, doc='Hotest spot temperature') # return m3.DHST[i] == m3.DHSTU[i] * m3.g2.value
def HST_E(m3, i): m3.Const_Tr_d = py.Constraint([m1.t_first.value], rule=HST_D, doc='Hotest spot temperature')
# return m3.DHST[i] == (m3.DHSTU[i] - m3.DHST[i - 1]) * m3.g2.value + m3.DHST[i - 1]
return m3.DHST[i] == (m3.DHSTU[i] - m3.DHST[i - 1]) * (1 - py.exp(-delta_t / m3.t_w.value)) + m3.DHST[i - 1]
def HST_E(m3, i):
# return m3.DHST[i] == (m3.DHSTU[i] - m3.DHST[i - 1]) * m3.g2.value + m3.DHST[i - 1]
m3.Const_Tr_e = py.Constraint(m1.t_second, rule=HST_E, doc='Hotest spot temperature') return m3.DHST[i] == (m3.DHSTU[i] - m3.DHST[i - 1]) * (1 - py.exp(-delta_t / m3.t_w.value)) + m3.DHST[i - 1]
# m = 0.8, n =0.9 or .8 ,.8, or 1,1 m3.Const_Tr_e = py.Constraint(m1.t_second, rule=HST_E, doc='Hotest spot temperature')
def HST_F(m3, i):
# return m3.DTOTU[i] == (
# m3.DTOR.value * ((m3.Power_P_T[i].value ** 2) * m3.TR.value + 1)) / ( # m = 0.8, n =0.9 or .8 ,.8, or 1,1
# m3.TR.value + 1) def HST_F(m3, i):
return m3.DTOTU[i] == m3.DTOR.value * ( # return m3.DTOTU[i] == (
(((m1.Power_Total_power[i].value / m3.Transformer_Limit[i]) ** 2) * m3.TR.value + 1) / ( # m3.DTOR.value * ((m3.Power_P_T[i].value ** 2) * m3.TR.value + 1)) / (
m3.TR.value + 1)) ** m3.n.value # m3.TR.value + 1)
return m3.DTOTU[i] == m3.DTOR.value * (
(((m1.Power_Total_power[i].value / m3.Transformer_Limit[i]) ** 2) * m3.TR.value + 1) / (
m3.Const_Tr_f = py.Constraint(m1.t, rule=HST_F, doc='Hotest spot temperature') m3.TR.value + 1)) ** m3.n.value
def HST_G(m3, i): m3.Const_Tr_f = py.Constraint(m1.t, rule=HST_F, doc='Hotest spot temperature')
return m3.DHSTU[i] == m3.DHR.value * (
((m1.Power_Total_power[i].value) / m3.Transformer_Limit[i]) ** m3.m.value) # was 2 , 2*0.8 = 1.6
def HST_G(m3, i):
return m3.DHSTU[i] == m3.DHR.value * (
m3.Const_Tr_g = py.Constraint(m1.t, rule=HST_G, doc='Hotest spot temperature') ((m1.Power_Total_power[i].value) / m3.Transformer_Limit[i]) ** m3.m.value) # was 2 , 2*0.8 = 1.6
# def HSTL(m3, i): m3.Const_Tr_g = py.Constraint(m1.t, rule=HST_G, doc='Hotest spot temperature')
# return m3.HST[i] <= m3.HSTL.value
#
# # def HSTL(m3, i):
# m3.Const_Tr_i = py.Constraint(m1.t, rule=HSTL, doc='Hotest upeer limit') # return m3.HST[i] <= m3.HSTL.value
# #
#
def Lol1(m3, i): # m3.Const_Tr_i = py.Constraint(m1.t, rule=HSTL, doc='Hotest upeer limit')
return m3.FAA[i] == py.exp(15000 / 383 - (15000 / (m3.HST[i] + 273))) #
def Lol1(m3, i):
m3.Const_Tr_j = py.Constraint(m1.t, rule=Lol1, doc='Aging factor') return m3.FAA[i] == py.exp(15000 / 383 - (15000 / (m3.HST[i] + 273)))
def Lol(m3, i): m3.Const_Tr_j = py.Constraint(m1.t, rule=Lol1, doc='Aging factor')
return m3.LOL[i] == (sum(m3.FAA[y] for y in m1.t) * m1.time_d.value) / 180000 # it was m1.time_d.value
# return m3.LOL[i] == (m3.FAA[i] * m1.time_d.value) / 180000
def Lol(m3, i):
return m3.LOL[i] == (sum(m3.FAA[y] for y in m1.t) * m1.time_d.value) / 180000 # it was m1.time_d.value
m3.Const_Tr_k = py.Constraint(m1.t, rule=Lol, doc='LOST OF LIFE transformer') # return m3.LOL[i] == (m3.FAA[i] * m1.time_d.value) / 180000
# LDRv(t)=E= SUM(td, LDR(t,td)) ; m3.Const_Tr_k = py.Constraint(m1.t, rule=Lol, doc='LOST OF LIFE transformer')
# MP(t).. Ptot(t)=E= P(t)-LLC(t)-b*PNSS-SUM(td, LDR(t,td))+SUM(td, LDR(td,t)) ;
# C_HST(t).. HST(t)=l=HSTL;# GAMS code
# C_LC(t).. LLC(t)=l=Ptot(t)-LLC(t);# HST1(t).. HST(t)=e=AT(t)+DTOT(t)++(t); # LDRv(t)=E= SUM(td, LDR(t,td)) ;
# MP(t).. Ptot(t)=E= P(t)-LLC(t)-b*PNSS-SUM(td, LDR(t,td))+SUM(td, LDR(td,t)) ;
# C_HST(t).. HST(t)=l=HSTL;# GAMS code
# OF .. Cost =E= SUM(t, LLC(t))*CLC + b*CSW + SUM(t, LDRv(t))*CDR ; # C_LC(t).. LLC(t)=l=Ptot(t)-LLC(t);# HST1(t).. HST(t)=e=AT(t)+DTOT(t)++(t);
def objective_rule(m3): # OF .. Cost =E= SUM(t, LLC(t))*CLC + b*CSW + SUM(t, LDRv(t))*CDR ;
return m3.LOL[i]
# return m3.LOL[i]
def objective_rule(m3):
return m3.LOL[i]
m3.obj_HEMS = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function') # return m3.LOL[i]
# <editor-fold desc="Solver">
# m3.write('mtr12.lp', io_options={'symbolic_solver_labels': True})
opt = py.SolverFactory('ipopt') m3.obj_HEMS = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function')
# opt = py.SolverFactory('mindtpy').solve(m3, mip_solver='cplex', nlp_solver='ipopt', tee=True) # <editor-fold desc="Solver">
# opt = py.SolverFactory('couenne ') # m3.write('mtr12.lp', io_options={'symbolic_solver_labels': True})
# opt.options["mipgap"] = 1 # 0.155 for load and 0.8 for other load new one: 0.045 opt = py.SolverFactory('ipopt')
result = opt.solve(m3, tee=True) # opt = py.SolverFactory('mindtpy').solve(m3, mip_solver='cplex', nlp_solver='ipopt', tee=True)
# result = opt.solve(m2) # opt = py.SolverFactory('couenne ')
# opt.options["mipgap"] = 1 # 0.155 for load and 0.8 for other load new one: 0.045
# this is for calling the neos result = opt.solve(m3, tee=True)
# os.environ['NEOS_EMAIL'] = 'sadamengr15@gmail.com' # result = opt.solve(m2)
# opt = py.SolverManagerFactory('neos')
# result = opt.solve(m2, opt='cplex') # this is for calling the neos
# os.environ['NEOS_EMAIL'] = 'sadamengr15@gmail.com'
# m2.pprint() # opt = py.SolverManagerFactory('neos')
# print(result) # result = opt.solve(m2, opt='cplex')
# print( "results.solver.status =", result.solver.status) # m2.pprint()
# print("results.solver.termination_condition =", result.solver.termination_condition ) # print(result)
# print("SolverStatus.warning =", result.warning)
print(result) # print( "results.solver.status =", result.solver.status)
# print("results.solver.termination_condition =", result.solver.termination_condition )
# print("SolverStatus.warning =", result.warning)
print(result)
# # # <editor-fold desc="Ploting">
# #
# # # Plotting for automatic ploting
time = [i for i in m1.t] # # # <editor-fold desc="Ploting">
# #
HST = [] # cost = [i for i in m2.c.values()] <<< check this # # # Plotting for automatic ploting
DHST = [] # cost = [i for i in m2.c.values()] <<< check this time = [i for i in m1.t]
DTOT = [] # cost = [i for i in m2.c.values()] <<< check this
Transformer_Limit = [] # cost = [i for i in m2.c.values()] <<< check this HST = [] # cost = [i for i in m2.c.values()] <<< check this
Power_P_T = [] # cost = [i for i in m2.c.values()] <<< check this DHST = [] # cost = [i for i in m2.c.values()] <<< check this
Power_T_P = [] # cost = [i for i in m2.c.values()] <<< check this DTOT = [] # cost = [i for i in m2.c.values()] <<< check this
Load_Transformer = [] # cost = [i for i in m2.c.values()] <<< check this Transformer_Limit = [] # cost = [i for i in m2.c.values()] <<< check this
Ambient_Transformer = [] # cost = [i for i in m2.c.values()] <<< check this Power_P_T = [] # cost = [i for i in m2.c.values()] <<< check this
FAA = [] # cost = [i for i in m2.c.values()] <<< check this Power_T_P = [] # cost = [i for i in m2.c.values()] <<< check this
LOL = [] # cost = [i for i in m2.c.values()] <<< check this Load_Transformer = [] # cost = [i for i in m2.c.values()] <<< check this
FAA_1 = [1] * 96 Ambient_Transformer = [] # cost = [i for i in m2.c.values()] <<< check this
HST_ref = [110] * 96 FAA = [] # cost = [i for i in m2.c.values()] <<< check this
for i in m3.LOL: LOL = [] # cost = [i for i in m2.c.values()] <<< check this
LOL.append(py.value(m3.LOL[i])) FAA_1 = [1] * 96
for i in m3.FAA: HST_ref = [110] * 96
FAA.append(py.value(m3.FAA[i])) for i in m3.LOL:
for i in m3.Ambient_Transformer: LOL.append(py.value(m3.LOL[i]))
Ambient_Transformer.append(py.value(m3.Ambient_Transformer[i])) for i in m3.FAA:
for i in m3.HST: FAA.append(py.value(m3.FAA[i]))
HST.append(py.value(m3.HST[i])) for i in m3.Ambient_Transformer:
for i in m3.DTOT: Ambient_Transformer.append(py.value(m3.Ambient_Transformer[i]))
DTOT.append(py.value(m3.DTOT[i])) for i in m3.HST:
for i in m3.DHST: HST.append(py.value(m3.HST[i]))
DHST.append(py.value(m3.DHST[i])) for i in m3.DTOT:
# for i in m1.Power_P_T: DTOT.append(py.value(m3.DTOT[i]))
# Power_P_T.append(py.value(m1.Power_P_T[i])) for i in m3.DHST:
for i in m3.Transformer_Limit: DHST.append(py.value(m3.DHST[i]))
Transformer_Limit.append(py.value(m3.Transformer_Limit[i])) # for i in m1.Power_P_T:
for i in m1.Power_Total_power: # Power_P_T.append(py.value(m1.Power_P_T[i]))
Load_Transformer.append(py.value(m1.Power_Total_power[i])) for i in m3.Transformer_Limit:
fig_T1, ax = plt.subplots(3, constrained_layout=True, ) # , sharey='row' Transformer_Limit.append(py.value(m3.Transformer_Limit[i]))
gs = fig_T1.add_gridspec(hspace=0, wspace=0) for i in m1.Power_Total_power:
# ax[0].plot(time, Transformer_Limit, color="red") Load_Transformer.append(py.value(m1.Power_Total_power[i]))
# ax[0].bar(time, Power_P_T, color="blue") fig_T1, ax = plt.subplots(3, constrained_layout=True, ) # , sharey='row'
# ax[0].legend(['Tranformer_limit']) gs = fig_T1.add_gridspec(hspace=0, wspace=0)
ax[0].plot(time, Transformer_Limit, color="red",label="Transformer limit") # ax[0].plot(time, Transformer_Limit, color="red")
ax[0].bar(time, Load_Transformer, color="blue",label="Load on Transformer") # ax[0].bar(time, Power_P_T, color="blue")
ax[0].legend() # ax[0].legend(['Tranformer_limit'])
ax[1].plot(time, HST, label='HST', color="red") ax[0].plot(time, Transformer_Limit, color="red",label="Transformer limit")
# ax[1].plot(time, DHST, label='DHST', color="blue") ax[0].bar(time, Load_Transformer, color="blue",label="Load on Transformer")
# ax[1].plot(time, DTOT, label='DTOT', color="green") ax[0].legend()
ax[1].plot(time, HST_ref, label='HST_ref',linestyle=':', color="black") ax[1].plot(time, HST, label='HST', color="red")
ax[1].legend() # ax[1].plot(time, DHST, label='DHST', color="blue")
ax[0].set_ylabel('Load_Transformer(kW)') # ax[1].plot(time, DTOT, label='DTOT', color="green")
ax[2].plot(time, FAA, label='FAA', color="red") ax[1].plot(time, HST_ref, label='HST_ref',linestyle=':', color="black")
ax[2].plot(time, FAA_1, label='FAA=1',linestyle=':' ,color="black") ax[1].legend()
ax[2].legend() ax[0].set_ylabel('Load_Transformer(kW)')
# ax[3].plot(time, LOL, label='% LOL', color="blue") ax[2].plot(time, FAA, label='FAA', color="red")
# ax[4].bar(HST, FAA, label='LOL', color="blue") ax[2].plot(time, FAA_1, label='FAA=1',linestyle=':' ,color="black")
# ax[5].bar(HST, LOL, label='LOL', color="blue") ax[2].legend()
ax[0].set_ylabel('Load_Transformer (kW)') # ax[3].plot(time, LOL, label='% LOL', color="blue")
ax[1].set_ylabel('Temp') # ax[4].bar(HST, FAA, label='LOL', color="blue")
ax[2].set_ylabel('Acc aging factor') # ax[5].bar(HST, LOL, label='LOL', color="blue")
# ax[3].set_ylabel('LoL %') ax[0].set_ylabel('Load_Transformer (kW)')
# ax[4].set_ylabel('HST Vs FAA') ax[1].set_ylabel('Temp')
# ax[5].set_ylabel('HST Vs %LoL') ax[2].set_ylabel('Acc aging factor')
fig_T1.suptitle('Transformer connected to HEMS', fontsize=16) # ax[3].set_ylabel('LoL %')
# # # ax[4].set_ylabel('HST Vs FAA')
# pl.tight_layout() # ax[5].set_ylabel('HST Vs %LoL')
fig_T1.suptitle('Transformer connected to HEMS', fontsize=16)
# #
# pl.tight_layout()
plt.show() plt.show()