$title myopic optimization model $onText myopic optimization model with heat storage and electrical level as parameters $offText $set matout "'gams_matlab.gdx' cost.l, R_m, R_p, Q_m, Q_p, t " Set t 'time periods' i 'thermal energy' /'c','h'/; *alias (t,lo) ; $gdxin gams_matlab.gdx $loadR t $gdxin display t,i ; *$exit $call GDXXRW data_opt_lookup_inc.xlsx Index=ParamsAndSets!A1 $GDXIN data_opt_lookup_inc.gdx ************************************************ parameters P_rer(t) 'Current amount of RER energy' a_grid (t) 'Current grid energy price' s_grid(t) 'Current grid supply disruption' E_L(t) 'Current energy demand of electrical load' Q_L(t,i) 'Current energy demand of thermal load' $LOAD P_rer a_grid s_grid E_L Q_L $gdxin display P_rer, a_grid, s_grid, E_L, Q_L; *********************************************** parameters P_rer_m 'model P_rer' a_grid_m 'model a_grid' s_grid_m 'model s_grid' E_L_m 'model E_L_m' Q_L_m(i) 'model thermal load' a_rer 'rer energy production cost'/0.03/ h_E ' electrical energy storage cost during period t' /0.0034/ b_E 'penalty cost for electrical energy '/1000/ h_T 'heating energy storage cost during period t'/0.00012/ b_T(i) 'penalty cost for thermal energy '/c 1000, h 1000 / beta_c_E 'electrical storage charging efficiency'/0.95 / beta_d_E 'electrical storage discharging efficiency'/0.95/ gamma_c_E 'electrical storage charging rate'/5/ gamma_d_E 'electrical storage discharging rate'/5/ beta_c_T 'Thermal storage charging efficiency'/0.95/ beta_d_T 'Thermal storage discharging efficiency'/0.95/ gamma_c_T 'Thermal storage charging rate'/7/ gamma_d_T 'Thermal storage discharging rate'/7/ P_grid_max'Grid maximal energy supply'/150/ Rmax 'Electrical storage maximal energy level'/100 / Qmax 'Thermal storage maximal energy level'/125/ Erl_max 'maximal electrical shedding load'/50/ Qrl_max(i)'maximal thermal shedding load'/c 150,h 150/ k1(i) 'binary number for heating load balance'/c 0,h 1/ k2(i) 'binary number for cooling load balance'/c 1,h 0/ PchL(i) 'load heating and cooling system'/c 200,h 200/ PchST 'storage heating system'/50/ k_per 'P_rer production control constant' /10/ k_grid 'grid price control constant' /1.5/ cop_rhp 'coefficient of performance heat pump' /3/ cop_ab 'coefficient of performance absorption chiller' /0.7/; parameter R_p 'previous electrical storage energy level at t=0' Q_p 'previous electrical storage energy level at t=0' R_m 'next energy level of electrical storage' Q_m 'next energy level of thermal storage' $gdxin gams_matlab.gdx $loadR Q_m, Q_p, R_m, R_p $gdxin display Q_m, Q_p, R_m, R_p; *$exit ************************************************ Positive Variable x1 'Energy trade-off from the grid to the electrical load' x2 'Energy trade-off from the grid to the electrical storage' x3(i) 'Energy trade-off from the grid to the thermal load chiller i' x4 'Energy trade-off from the grid to the thermal storage chiller' x5 'energy trade-off from the RER to the electrical load' x6 'Energy trade-off from the RER to the electrical storage' x7(i) 'Energy trade-off from the RER to the thermal load i' x8 'Energy trade-off from the RER to the thermal storage chiller' x9 'Energy trade-off from the electrical storage to the electrical load' x10(i)'Energy trade-off from the thermal load chiller to the thermal load i' x11(i)'Energy trade-off from the thermal load chiller to the thermal storage' x12(i)'Energy trade-off from the thermal storage to the thermal load i' QchST 'Thermal storage chiller production'; variables cost 'cost function'; ************************************************ Equations eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15 eq16,eq17,eq18,eq19,eq20,eq21; *cost function eq1.. cost=E=k_grid*a_grid_m*(x1+x2+sum(i, x3(i))+x4) +a_rer*(x5+x6+sum(i, x7(i))+x8) +h_E*(R_p+beta_c_E*(x2+x6)-x9) +h_T *(Q_p+beta_c_T*(QchST+sum(i,k1(i)*x11(i)))- sum(i,k1(i)*x12(i)+k2(i)*x12(i))) +b_E*(E_L_m-x1-x5-x9*beta_d_E) +sum(i,b_T(i)*(Q_L_m(i)-(x10(i)+beta_d_T*k1(i)*x12(i)+beta_d_T*k2(i)*x12(i)*cop_ab))); *transition function of electrical storage eq2..R_m =e= R_p + beta_c_E*(x2 + x6) - x9; *transition function of heating storage eq3..Q_m =e= Q_p + beta_c_T*(QchST + sum(i,k1(i) * x11(i)))-sum(i,k1(i)*x12(i)+k2(i)*x12(i)); *electrical load demand balance eq4.. x1+x5+beta_d_E*x9=l= E_L_m; *thermal load demand balance eq5(i).. x10(i)+ beta_d_T*k1(i)*x12(i)+beta_d_T*k2(i)*x12(i)*cop_ab=l=Q_L_m(i); *cop_absorption chiller *RER production and decision balance eq6.. x5+x6+sum(i, x7(i))+x8=l=k_per*P_rer_m; *grid energy supply capacity and and decision balance eq7.. x1+x2+sum(i, x3(i))+x4=l=s_grid_m*P_grid_max; *Electrical storage charging rate limitation eq8.. x2+x6=l=gamma_c_E; *Electrical storage discharging rate limitation eq9.. x9=l=gamma_d_E; eq10.. x9=l=R_p; *Thermal storage charging rate limitation eq11.. sum(i,k1(i)*x11(i))+QchST=l=gamma_c_T; *Thermal storage discharging rate limitation eq12.. sum (i,x12(i))=l=gamma_d_T; eq13.. sum (i,x12(i))=l=Q_p; *Electrical storage charging limitation eq14.. beta_c_E*(x2+x6)=l=Rmax-R_p; *Thermal storage charging limitation eq15.. beta_c_T*(QchST+sum(i,k1(i)*x11(i)))=l=Qmax-Q_p; *heating storage chiller heating production eq16..QchST=e=cop_rhp*(x4+x8); *load chiller power consumption limitation eq17(i).. x3(i)+x7(i)=l=PchL(i); *heating storage chiller power consumption limitation eq18.. x4+x8=l=PchST; *electrical load shedding limitation eq19.. E_L_m-(x1+x5+x9*beta_d_E)=l= Erl_max ; *Thermal load shedding limitation eq20(i).. Q_L_m(i)-(x10(i)+beta_d_T*k1(i)*x12(i)+ beta_d_T*k2(i)* x12(i)*cop_ab)=l= Qrl_max(i); *thermal load chiller energy balance eq21(i)..x10(i)+ k1(i)*x11(i)=e=cop_rhp*(x3(i)+x7(i)); ************************************************************ model myopic_lookup1 /all/; *********************************************************** Options LP=CPlex; Parameter report1(t,*),report2(t,i,*), report3(t,i,*), report4 ; $exit ********************************************************** loop(t, *execute_load 'gams_matlab.gdx' Q_m, Q_p, R_m, R_p; *display Q_m, Q_p, R_m, R_p; P_rer_m = P_rer(t); a_grid_m = a_grid(t); s_grid_m = s_grid(t); E_L_m = E_L(t); Q_L_m(i)=Q_L(t,i) ; solve myopic_lookup1 using LP minimizing cost; report1(t,'x1') = x1.l; report1(t,'x2') = x2.l; report1(t,'x4') = x4.l; report1(t,'x5') = x5.l; report1(t,'x6') = x6.l; report1(t,'x8') = x8.l; report1(t,'x9') = x9.l; report1(t,'QchST') = QchST.l; report1(t,'grid') = x1.l+x2.l+x4.l+sum(i,x3.l(i)); report1(t,'RER') = x5.l+x6.l+x8.l+sum(i,x7.l(i)); report1(t,'E_load') = x1.l+beta_d_E*x9.l+x5.l; report1(t,'E_shp') = x8.l+x4.l; report1(t,'RER_p') = k_per*P_rer(t); report1(t,'RER_re') = k_per*P_rer(t)-(x5.l+x6.l+x8.l+sum(i,x7.l(i))); report1(t,'E_shedding') = E_L(t)-(x1.l+beta_d_E*x9.l+x5.l); report1(t,'grid_cost(t)') = k_grid*a_grid_m*(x1.l+x2.l+sum (i, x3.l(i))+x4.l); report1(t,'rer_cost(t)') = a_rer*(x5.l+x6.l+sum (i, x7.l(i))+x8.l); report1(t,'Es_cost(t)') = h_E*(R_p+beta_c_E*(x2.l+x6.l)-x9.l); report1(t,'Qs_cost(t)') = h_T *(Q_p+beta_c_T*(QchST.l+sum(i,k1(i)*x11.l(i)))- sum(i,k1(i)*x12.l(i)+k2(i)*x12.l(i))); report1(t,'Eshed_cost(t)') = b_E*(E_L_m-x1.l-x5.l-x9.l*beta_d_E); report1(t,'Qshed_cost(t)') = sum(i,b_T(i)*(Q_L_m(i)-(x10.l(i)+beta_d_T*k1(i)*x12.l(i)+ beta_d_T*k2(i)* x12.l(i)*cop_ab))); report2(t,'c','x3') = x3.l('c'); report2(t,'c','x7') = x7.l('c'); report2(t,'c','x10') = x10.l('c'); report2(t,'c','x11') = x11.l('c'); report2(t,'c','x12') = x12.l('c'); report2(t,'c','Q_load_c') = x10.l('c')+beta_d_T*x12.l('c')*cop_ab; report2(t,'c','Q_shedding_c')= Q_L_m('c')-(x10.l('c')+beta_d_T*x12.l('c')*cop_ab); report3(t,'h','x3') = x3.l('h'); report3(t,'h','x7') = x7.l('h'); report3(t,'h','x10') = x10.l('h'); report3(t,'h','x11') = x11.l('h'); report3(t,'h','x12') = x12.l('h'); report3(t,'h','Q_load_h') = x10.l('h')+beta_d_T*x12.l('h'); report3(t,'h','Q_shedding_h')= Q_L_m('h')-(x10.l('h')+beta_d_T*x12.l('h')); report1(t,'cost') = k_grid*a_grid_m*(x1.l+x2.l+sum (i, x3.l(i))+x4.l) +a_rer*(x5.l+x6.l+sum (i, x7.l(i))+x8.l) +h_E*(R_p+beta_c_E*(x2.l+x6.l)-x9.l) +h_T *(Q_p+beta_c_T*(QchST.l+sum(i,k1(i)*x11.l(i)))- sum(i,k1(i)*x12.l(i)+k2(i)*x12.l(i))) +b_E*(E_L_m-x1.l-x5.l-beta_d_E*x9.l) +sum(i,b_T(i)*(Q_L_m(i)-(x10.l(i)+beta_d_T*k1(i)*x12.l(i)+ beta_d_T*k2(i)*x12.l(i)*cop_ab))); ); report4('total_cost')=sum(t, report1(t,'cost')); report4('grid_cost')=sum(t, report1(t,'grid_cost(t)')); report4('rer_cost')=sum(t, report1(t,'rer_cost(t)')); report4('rer_cost')=sum(t, report1(t,'rer_cost(t)')); report4('Es_cost')=sum(t, report1(t,'Es_cost(t)')); report4('Qs_cost')=sum(t, report1(t,'Qs_cost(t)')); report4('Eshed_cost')=sum(t, report1(t,'Eshed_cost(t)')); report4('Qshed_cost')=sum(t, report1(t,'Qshed_cost(t)')); report4('E_shedding')=sum(t, report1(t,'E_shedding')); report4('E_shp')=sum(t, report1(t,'E_shp')); report4('Q_shedding_h')=sum(t, report3(t,'h','Q_shedding_h')); report4('Q_shedding_c')=sum(t, report2(t,'c','Q_shedding_c')); display report1, report2,report3, report4; execute_unload %matout%; $exit