commit 30bd12534a50f7f2d001492b9684f8479c2c30e3 Author: Sadam Hussain Date: Wed May 1 11:38:07 2024 -0400 Upload files to "/" diff --git a/transformer_model.txt b/transformer_model.txt new file mode 100644 index 0000000..98ed955 --- /dev/null +++ b/transformer_model.txt @@ -0,0 +1,228 @@ +m3 = py.ConcreteModel() + +# tranformer parameters +m3.Transformer_Limit = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_data.Transformer_limit)), + doc='Maximum power limit on transformer') +m3.tran_life = py.Param(mutable=True, doc='Transformer life expectancy') +m3.HST = py.Var(m1.t, bounds=(0, None), doc="windign hottest spot temp") +m3.LOL = py.Var(m1.t, bounds=(0, None), doc="lost of life") +# m3.FAA = py.Var(m1.t, bounds=(0, None), doc="HST") +m3.DTOT = py.Var(m1.t, bounds=(0, None), doc="top oil rise temp") +m3.DHST = py.Var(m1.t, bounds=(0, None), doc="winding HST rise over the top") +m3.DTOTU = py.Var(m1.t, bounds=(0, None), doc="ultimate top-oil rise over the ambient temperature") +m3.DHSTU = py.Var(m1.t, bounds=(0, None), doc="ultimate top-oil rise over the ambient") +m3.Power_P_T = py.Var(m1.t, bounds=(0, None), doc="power from home to Transformer") +m3.Power_T_P = py.Var(m1.t, bounds=(0, None), doc="power from transformer") +m3.S_P_T = py.Var(m1.h, m1.t, within=py.Binary) +m3.FAA = py.Var(m1.t, bounds=(0, None)) # The output of piecewise. Equivalent to model.y=piecewise(model.x) +m3.n = py.Param(initialize=0.9, + doc='top oil rise constant for transformer ') # 1, 0.9=ONOF tranformer, 0.8=ONON Transformer, it was 1 +m3.m = py.Param(initialize=1.6, + 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.1535, + 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.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.g2 = py.Param(initialize=0.8826, + 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.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.TR = py.Param(initialize=6, doc='Resistance ') +m3.DTOR = py.Param(initialize=56, doc='e top-oil rise over ambient at the rated load ') # it was 56 +m3.DHR = py.Param(initialize=80, doc='hottest-spot rise over top-oil at the rated load ') # 24 and 80 it was 24 +m3.HSTL = py.Param(initialize=110, doc='time duration ') +m3.t_TO = py.Param(initialize=180, doc=' in min ') # it was 90min in one paper and other paper 180 min +m3.t_w = py.Param(initialize=48, doc=' in min ') # it was 7min in one paper and other paper 48 min +m3.Ambient_Transformer = py.Param(m1.t, initialize=dict(zip(sheet_data.time, sheet_tr.celsius1)), mutable=True, + 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) +delta_t = 15 + + +# +def HST_A(m3, i): + return m3.HST[i] == m3.Ambient_Transformer[i] + m3.DTOT[i] + m3.DHST[i] + + +m3.Const_Tr_a = py.Constraint(m1.t, rule=HST_A, doc='Hotest spot temperature') + + +def HST_B(m3, i): + # 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)) + + +m3.Const_Tr_b = py.Constraint([m1.t_first.value], rule=HST_B, doc='Hotest spot temperature') + + +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] + # return m3.DTOT[i] == (m3.DTOTU[i] - m3.DTOT[i - 1]) * m3.g1 + m3.DTOT[i - 1] + + +m3.Const_Tr_c = py.Constraint(m1.t_second, rule=HST_C, doc='Hotest spot temperature') + + +def HST_D(m3, i): + 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 + + +m3.Const_Tr_d = py.Constraint([m1.t_first.value], rule=HST_D, doc='Hotest spot temperature') + + +def HST_E(m3, i): + # 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] + + +m3.Const_Tr_e = py.Constraint(m1.t_second, rule=HST_E, doc='Hotest spot temperature') + + +# m = 0.8, n =0.9 or .8 ,.8, or 1,1 +def HST_F(m3, i): + # return m3.DTOTU[i] == ( + # m3.DTOR.value * ((m3.Power_P_T[i].value ** 2) * m3.TR.value + 1)) / ( + # 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.TR.value + 1)) ** m3.n.value + + +m3.Const_Tr_f = py.Constraint(m1.t, rule=HST_F, doc='Hotest spot temperature') + + +def HST_G(m3, i): + 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 + + +m3.Const_Tr_g = py.Constraint(m1.t, rule=HST_G, doc='Hotest spot temperature') + + +# def HSTL(m3, i): +# return m3.HST[i] <= m3.HSTL.value +# +# +# m3.Const_Tr_i = py.Constraint(m1.t, rule=HSTL, doc='Hotest upeer limit') +# + +def Lol1(m3, i): + return m3.FAA[i] == py.exp(15000 / 383 - (15000 / (m3.HST[i] + 273))) + + +m3.Const_Tr_j = py.Constraint(m1.t, rule=Lol1, doc='Aging factor') + + +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 + # return m3.LOL[i] == (m3.FAA[i] * m1.time_d.value) / 180000 + + +m3.Const_Tr_k = py.Constraint(m1.t, rule=Lol, doc='LOST OF LIFE transformer') + + +# 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 +# C_LC(t).. LLC(t)=l=Ptot(t)-LLC(t);# HST1(t).. HST(t)=e=AT(t)+DTOT(t)++(t); + + +# OF .. Cost =E= SUM(t, LLC(t))*CLC + b*CSW + SUM(t, LDRv(t))*CDR ; + + +def objective_rule(m3): + return m3.LOL[i] + # return m3.LOL[i] + + +m3.obj_HEMS = py.Objective(rule=objective_rule, sense=py.minimize, doc='Definition of objective function') +# +# m3.write('mtr12.lp', io_options={'symbolic_solver_labels': True}) +opt = py.SolverFactory('ipopt') +# opt = py.SolverFactory('mindtpy').solve(m3, mip_solver='cplex', nlp_solver='ipopt', tee=True) +# opt = py.SolverFactory('couenne ') +# opt.options["mipgap"] = 1 # 0.155 for load and 0.8 for other load new one: 0.045 +result = opt.solve(m3, tee=True) +# result = opt.solve(m2) + +# this is for calling the neos +# os.environ['NEOS_EMAIL'] = 'sadamengr15@gmail.com' +# opt = py.SolverManagerFactory('neos') +# result = opt.solve(m2, opt='cplex') + +# m2.pprint() +# 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) + + + +# # # +# # +# # # Plotting for automatic ploting +time = [i for i in m1.t] + +HST = [] # cost = [i for i in m2.c.values()] <<< check this +DHST = [] # cost = [i for i in m2.c.values()] <<< check this +DTOT = [] # cost = [i for i in m2.c.values()] <<< check this +Transformer_Limit = [] # cost = [i for i in m2.c.values()] <<< check this +Power_P_T = [] # cost = [i for i in m2.c.values()] <<< check this +Power_T_P = [] # cost = [i for i in m2.c.values()] <<< check this +Load_Transformer = [] # cost = [i for i in m2.c.values()] <<< check this +Ambient_Transformer = [] # cost = [i for i in m2.c.values()] <<< check this +FAA = [] # cost = [i for i in m2.c.values()] <<< check this +LOL = [] # cost = [i for i in m2.c.values()] <<< check this +FAA_1 = [1] * 96 +HST_ref = [110] * 96 +for i in m3.LOL: + LOL.append(py.value(m3.LOL[i])) +for i in m3.FAA: + FAA.append(py.value(m3.FAA[i])) +for i in m3.Ambient_Transformer: + Ambient_Transformer.append(py.value(m3.Ambient_Transformer[i])) +for i in m3.HST: + HST.append(py.value(m3.HST[i])) +for i in m3.DTOT: + DTOT.append(py.value(m3.DTOT[i])) +for i in m3.DHST: + DHST.append(py.value(m3.DHST[i])) +# for i in m1.Power_P_T: +# Power_P_T.append(py.value(m1.Power_P_T[i])) +for i in m3.Transformer_Limit: + Transformer_Limit.append(py.value(m3.Transformer_Limit[i])) +for i in m1.Power_Total_power: + Load_Transformer.append(py.value(m1.Power_Total_power[i])) +fig_T1, ax = plt.subplots(3, constrained_layout=True, ) # , sharey='row' +gs = fig_T1.add_gridspec(hspace=0, wspace=0) +# ax[0].plot(time, Transformer_Limit, color="red") +# ax[0].bar(time, Power_P_T, color="blue") +# ax[0].legend(['Tranformer_limit']) +ax[0].plot(time, Transformer_Limit, color="red",label="Transformer limit") +ax[0].bar(time, Load_Transformer, color="blue",label="Load on Transformer") +ax[0].legend() +ax[1].plot(time, HST, label='HST', color="red") +# ax[1].plot(time, DHST, label='DHST', color="blue") +# ax[1].plot(time, DTOT, label='DTOT', color="green") +ax[1].plot(time, HST_ref, label='HST_ref',linestyle=':', color="black") +ax[1].legend() +ax[0].set_ylabel('Load_Transformer(kW)') +ax[2].plot(time, FAA, label='FAA', color="red") +ax[2].plot(time, FAA_1, label='FAA=1',linestyle=':' ,color="black") +ax[2].legend() +# ax[3].plot(time, LOL, label='% LOL', color="blue") +# ax[4].bar(HST, FAA, label='LOL', color="blue") +# ax[5].bar(HST, LOL, label='LOL', color="blue") +ax[0].set_ylabel('Load_Transformer (kW)') +ax[1].set_ylabel('Temp') +ax[2].set_ylabel('Acc aging factor') +# ax[3].set_ylabel('LoL %') +# ax[4].set_ylabel('HST Vs FAA') +# ax[5].set_ylabel('HST Vs %LoL') +fig_T1.suptitle('Transformer connected to HEMS', fontsize=16) +# # +# pl.tight_layout() +plt.show() \ No newline at end of file