import os
import sys
import math
import cvxopt as cvx
import picos as pic
import pandas as pd
import matplotlib.pyplot as plt
from gurobipy import *
from statsmodels.tsa.arima_model import ARIMA
import numpy as np
from scipy import *
#import DeferableLoad
OPTmodel = Model('OPTIMIZER')
#general parameters
Tamb =22
N = 1440 # maximum iteration
i = range(1, N)
COP= 3.4 # Coeffient of performance
'''
Prediction need to be added here
'''
# Datacenter room defintion
R = 10 #length of room
B = 7
H = 9 #Height of room
L = 10
dT = 60
A = 2*((L*B)+(B*H)+(H*L))
Thick = 0.33 # thickness of wall
k = 0.7 # thermal conductivity of wall
mAir = 1.2 * (L * B * H)
C = 718
landa = k * A / Thick
a0 = 0.05 / dT
a1 = 1
ki = math.exp(-(landa * 60) / (mAir * C)) # value that constant and its related to property of room
kc = (1 - ki) * a0
ko = (1 - ki) * a1
kp = (1 - ki) * (COP / landa)
Tmin= 18
Tmax= 27
Tamb= 22
PcoolingRated = 100
Pbess_rated = 30.462
Pbess_ratedN = -30.462
Ebess_min = 0
Ebess_max = 300
with open ('Pcooling.csv','r') as f:
Pcooling = []
for line in f:
Pcooling.append(line)
f.close()
with open ('ITpower.csv','r') as f1:
ITload = []
for line1 in f1:
ITload.append(line1)
f1.close()
with open ('DR.csv','r') as f2:
DR =[]
for line2 in f2:
DR.append(line2)
f2.close()
print ITload
print Pcooling
print DR
for i in range(1,200):
for it in range(1, 1440):
Tm = np.empty(1440)
Tm.fill(18)
TmA = np.empty(1440)
TmA.fill(27)
Phvac_flex = {}
Phvac_up = {}
Phvac_down_= {}
Phvac_up_ = {}
Pbess_out_ = {}
Pbess_in_ = {}
Phvac_down = {}
Pbess_flex_ = {}
Pbess_flex = {}
Phvac_flex_ = {}
Pbess_in = {}
Pdc = {}
Pdc_base = {}
Pflex_i = {}
Tdc_i = {}
Pbess_out ={}
Ebess_i = {}
Phvac_flex[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name="PHVAC_flex"+str(i))
Phvac_up[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PHVAC_up" + str(i))
Phvac_up_[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PHVAC_up_" + str(i))
Phvac_down_[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PHVAC_down_" + str(i))
Pbess_out_[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PBESS_out_" + str(i))
Pbess_in_[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PBESS_in_" + str(i))
Phvac_down[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PHVAC_down" + str(i))
Pbess_flex_[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PBESS_flex_" + str(i))
Pbess_flex[i] = OPTmodel.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PBESS_flex" + str(i))
Phvac_flex_[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PHVAC_flex_" + str(i))
Pbess_in[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PBESS_in" + str(i))
Pdc[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PDC" + str(i))
Pdc_base[i] = OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PDC_base" + str(i))
Pflex_i[i]= OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="Pflex_i" + str(i))
Tdc_i[i]= OPTmodel.addVar(ub=GRB.INFINITY,vtype = GRB.CONTINUOUS, name = "Tdc_i" + str(i))
Pbess_out[i] = OPTmodel.addVar(lb=-GRB.INFINITY,ub=GRB.INFINITY,vtype=GRB.CONTINUOUS, name="PBESS_out" + str(i))
Ebess_i[i]= OPTmodel.addVar(ub=GRB.INFINITY,vtype=GRB.CONTINUOUS,name="Ebess_i" + str(i))
Pflex_i[1] = 0
Pflex_i[1] = 0
Tdc_i[0] = 18
Phvac_flex[1] = 0
# Phvac_flex_[1] = 0
Phvac_down[1] = 0
Phvac_up[1] = 0
Phvac_down_[1] = 0
Phvac_up_[1] = 0
# Phvac_down_pos[1] = 0
# Phvac_up_pos(1) = 0;
Pbess_flex[1] = 0
# Pbess_flex_[1] = 0
Pbess_out[1] = 0
Pbess_in[1] = 0
# Pbess_out_[1] = 0
Pbess_in_[1] = 0
# Pbess_out_pos[1] = -250
# Pbess_in_pos(1) = 250;
Ebess_i[1] = 150
OPTmodel.update()
'''
if float(DR[i]) > 0:
Phvac_down_[i] = 0
Phvac_up_[i] = float(DR[i])
Pbess_out_[i] = 0
Pbess_in_[i] = float(DR[i])
#Pbess_flex_[i] = Pbess_in_[i] + Pbess_out_[i]
#Phvac_flex_[i] = Phvac_down_[i] + Phvac_up_[i]
OPTmodel.update()
elif float(DR[i]) < 0:
Phvac_down_[i] = float(DR[i])
Phvac_up_[i] = 0
#Phvac_flex_[i] = Phvac_down_[i] + Phvac_up_[i]
Pbess_out_[i] = float(DR[i])
Pbess_in_[i] = 0
#Pbess_flex_[i] = Pbess_in_[i] + Pbess_out_[i]
OPTmodel.update()
else:
Phvac_down_[i] = 0
Phvac_up_[i] = 0
Phvac_flex_[i] = Phvac_down_[i] + Phvac_up_[i]
Pbess_out_[i] = 0
Pbess_in_[i] = 0
Pbess_flex_[i] = Pbess_in_[i] + Pbess_out_[i]
OPTmodel.update()
'''
#print Phvac_up.values()
#print Phvac_flex_[i]
print OPTmodel
OPTmodel.update()
ConHVAC1 = OPTmodel.addConstr(Phvac_flex[i] == Phvac_up[i] + Phvac_down[i], name='ConHVAC1')
ConHVAC2 = OPTmodel.addConstr(0 <= Phvac_flex[i] , name='ConHVAC2')
ConHVAC3 = OPTmodel.addConstr(Phvac_flex[i] <= PcoolingRated, name='ConHVAC3')
PH = pd.read_csv('Pcooling.csv')
PHVAC = PH.values
newList2 = map(lambda x: x / 1000, PHVAC)
p=[]
p=PcoolingRated-newList2[i]
#CONHVAC4 = OPTmodel.addConstr(Phvac_up[i]==np.minimum((Phvac_up_[i]),(float(newList2[i]))))
#Phvac_u(1:MaxIter) == min(Phvac_u_(1:MaxIter), (repelem(Phvac_max, MaxIter) - (Pcooling(1:MaxIter)'/1000)))
ConTemp1 = OPTmodel.addConstr(Tm[it] <= Tdc_i[i] <= TmA[it], name='ConTemp1')
ConBESS1 = OPTmodel.addConstr(Pbess_ratedN <= Pbess_flex[i] <= Pbess_rated, name='ConBESS1')
ConBESS2 = OPTmodel.addConstr(Pbess_flex[i] == Pbess_in[i] + Pbess_out[i], name='ConBESS2')
ConBESS3 = OPTmodel.addConstr(0 <= Pbess_in[i] <= min(Pbess_rated, Pbess_in_[i]), name='ConBESS3')
ConBESS4 = OPTmodel.addConstr(np.maximum(Pbess_ratedN,Pbess_out_[i]) <= Pbess_out[i]<=0 , name='ConBESS4') # need to modifty
ConEBESS1 = OPTmodel.addConstr(Ebess_min <= Ebess_i[i], name='ConEBESS1')
ConEBESS2 = OPTmodel.addConstr(Ebess_i[i] <= Ebess_max, name='ConEBESS2')
D = pd.read_csv('DR.csv').values
DRN = map(lambda x: x / 1000, D)
PDRN=map(lambda x: x / 4.8, DRN)
if float((PDRN[i])) > 0:
CON1 = OPTmodel.addConstr(Pbess_flex_[i] == Pbess_in_[i] + Pbess_out_[i],'CON1')
CON2 = OPTmodel.addConstr(Phvac_flex_[i] == Phvac_up_[i] + Phvac_down_[i],'CON2')
CON3=OPTmodel.addConstr(Phvac_down_[i] == 0, name='CON3')
CON4=OPTmodel.addConstr(Phvac_up_[i] == float((PDRN[i])),name='CON4')
CON5=OPTmodel.addConstr(Pbess_out_[i] == 0,name='CON5')
CON6=OPTmodel.addConstr(Pbess_in_[i] == float((PDRN[i])),name='CON6')
elif float(np.transpose(PDRN[i])) < 0:
CON7=OPTmodel.addConstr(Phvac_down_[i] == float(np.transpose(PDRN[i])),name='CON7')
CON8=OPTmodel.addConstr(Phvac_up_[i] == 0,name='CON8')
# Phvac_flex_[i] = Phvac_down_[i] + Phvac_up_[i]
CON9=OPTmodel.addConstr(Pbess_out_[i] == float((PDRN[i])),name='CON9')
CON10=OPTmodel.addConstr(Pbess_in_[i] == 0,name='CON10')
else:
CON11=OPTmodel.addConstr(Phvac_down_[i] == 0,name='CON11')
CON12=OPTmodel.addConstr(Phvac_up_[i] == 0,name='CON12')
CON13=OPTmodel.addConstr(Phvac_flex_[i] == Phvac_down_[i] + Phvac_up_[i],name='CON13')
CON14=OPTmodel.addConstr(Pbess_out_[i] == 0)
CON15=OPTmodel.addConstr(Pbess_in_[i] == 0,name='CON15')
CON16=OPTmodel.addConstr(Pbess_flex_[i] == Pbess_in_[i] + Pbess_out_[i],name='CON16')
OPTmodel.update()
ConPDC = OPTmodel.addConstr(Pdc[i] == Pflex_i[i] + float(ITload[i]), name='ConPDC')
# OPTmodel.addConstr(Tdc_i[i]==(ki*Tdc_i[i-1]+(ko*Tamb)))
#for x in Ebess_i:
#ConEBESS2 = OPTmodel.addConstr(Ebess_i[i] ==((Pbess_in[i] / 0.75) + (Pbess_out[i] * 0.75)))
cooling = np.array(pd.read_csv('Pcooling.csv'))
DRR = pd.read_csv('DR.csv')
DR = DRR.values
IT = pd.read_csv('ITpower.csv')
ITload = IT.values
newList = map(lambda x: x / 1000, ITload)
PH = pd.read_csv('Pcooling.csv')
PHVAC = PH.values
newList2 = map(lambda x: x / 1000, PHVAC)
#for y in Tdc_i:
T=pd.read_csv('TT.csv').values
OPTmodel.addConstr(Tdc_i[i]==((ki*float(T[i]))+(ko*Tamb)+(kc*float(newList[i]))-((kp*(float(newList2[i])))+(Phvac_flex[i]*3.14))))
print Tdc_i.values()
OPTmodel.addConstr(Pbess_out_[i]<=Phvac_flex[i] + Pbess_flex[i]<=Pbess_in_[i])
# Tdc_i[1:len(i)]==(Ki*Tdc_i[1:1438])+(Kc*array2[1:1438])+(Ko*Tamb))
ConBESS5 = OPTmodel.addConstr(Pbess_flex[i] == Pbess_in[i] + Pbess_out[i], name='ConBESS5')
#OPTmodel.addConstr(defIT[i]==DeferableLoad.j2 + DeferableLoad.j3)
# OPTmodel.addConstr(Pdc_base[i]==predictions[i])
ConFLEX = OPTmodel.addConstr(Pflex_i[i] == Pbess_flex[i] + Phvac_flex[i], name='ConFLEX')
PcoolingPredicted = pd.read_csv('PcoolingPredictionResult.csv')
PcoolingPredictedValue = PcoolingPredicted.values
ITPredicted = pd.read_csv('ITpredictionResult.csv')
ITPredictedValue = ITPredicted.values
ConPDCbase = OPTmodel.addConstr(Pdc_base[i] == np.transpose(ITPredictedValue[i]) + np.transpose(PcoolingPredictedValue[i]))
OPTmodel.update()
# OPTmodel.addConstr(Pdc_base[i]==prediction[i])
OPTmodel.setObjective((np.transpose(Pdc_base[i])-float(DR[i]) - (Pdc[i]) ), GRB.MINIMIZE)
OPTmodel.update()
OPTmodel.optimize()
print Pdc_base[i].X
#print Ebess_i[i].X
#print Phvac_flex[i].X
print Tdc_i[i]
print Pdc[i]
print Phvac_flex[i]
print Pbess_flex[i]
print Pbess_out[i]
print Pbess_in[i]
print Ebess_i[i]
print Pbess_flex_[i]
print Phvac_down[i]
print Phvac_up[i]
'''
def get_results(self):
"""
This function gets the results of the current optimization model
Returns
-------
"""
HVACresult = np.zeros(1,N)
BatteryResult = np.zeros(1,N)
SOC = np.zeros(1,N)
#r_Q_dot = np.zeros((self.gp.N_H, self.N_S))
#r_P = np.zeros((self.gp.N_H, self.N_S))
#r_P_self = np.zeros((self.gp.N_H, self.N_S))
#r_P_ex = np.zeros((self.gp.N_H, self.N_S))
#r_Q_dot_gas = np.zeros((self.gp.N_H, self.N_S))
#Load = np.zeros((self.gp.N_H, self.N_S))
try:
for t in range(1,N):
HVACresult[t]= Phvac_flex[t].X
BatteryResult[t]=Pbess_flex[t].X
SOC[t] = Ebess_i[t].X / Ebess_max
except:
pass
return { 'SOC' : SOC , 'BatteryResult': BatteryResult }
print OPTmodel.getVars()
# get results
Temp = {}
Battery = {}
Ebess_result = {}
ITloadd = {}
for t in range(1,N):
Temp[t] = OPTmodel.getVarByName("Tdc_i" )
Battery[t] = OPTmodel.getVarByName("PBESS_flex" )
Ebess_result[t] = OPTmodel.getVarByName("Ebess_i" )
#r_P_e[t] = model.getVarByName("P_export_%s_0" % t).X
fig, axes = plt.subplots(4, 1)
# plot elctricity
ax5 = axes[2]
ax6 = ax5.twinx()
ax5.plot( [Temp[t] for t in range(1,N)], 'g-')
ax6.plot([Ebess_result[t] for t in range(1,N)], 'b-')
ax5.set_xlabel('Time index')
ax5.set_ylabel('Power Import [W]', color='g')
ax6.set_ylabel('Power CHP [W]', color='b')
ax7 = axes[3]
ax7.plot([Battery[t] for t in range(1,N)], 'g-')
ax7.set_ylabel('Power Export [W]', color='g')
'''
print Pflex_i.values()
# print OPTmodel.getVars()
print OPTmodel.feasibility()
print OPTmodel.getObjective()
print Pdc_base.values()
'''
b = map(float, Phvac_flex)
plt.plot(b)
plt.show()
'''
#c = map(float, Pbess_flex_)
#plt.plot(c)
#plt.show()
print OPTmodel
print Tdc_i.values()
# get results
print OPTmodel.getVars()
# print OPTmodel.getAttr('EBESS_i')
status = OPTmodel.status
print status
# print Con10,Con12
print Phvac_flex.values()
print Pbess_flex.values()
print Ebess_i.values()
print OPTmodel.objval
print Tdc_i
print Pbess_in
print Pbess_out.values()
# print Pbess_flex
# print Phvac_flex
# print Ebess_i
print Pflex_i.values()
print Pbess_flex_.values()
#print OPTmodel.getVars()
print OPTmodel.feasibility()
print OPTmodel.getObjective()
print Ebess_i.values()
if OPTmodel.status == GRB.Status.INF_OR_UNBD:
# Turn presolve off to determine whether model is infeasible
# or unbounded
OPTmodel.setParam(GRB.Param.Presolve, 0)
OPTmodel.optimize()
OPTmodel.write("mymodel.lp")
if OPTmodel.status == GRB.Status.OPTIMAL:
print('Optimal objective: %g' % OPTmodel.objVal)
OPTmodel.write('model.sol')
exit(0)
elif OPTmodel.status != GRB.Status.INFEASIBLE:
print('Optimization was stopped with status %d' % OPTmodel.status)
exit(0)
# Model is infeasible - compute an Irreducible Inconsistent Subsystem (IIS)
print('')
print('Model is infeasible')
OPTmodel.computeIIS()
OPTmodel.write("model.ilp")
print("IIS written to file 'model.ilp'")
I want to plot the computed values from gurobi but when I want to get the X attribute of gurobi variable it says that AttributeError: it has no attribute 'X' and the when I cast the value from float to int it just showed me the empty plot but at the lp file I could see the result of each iteration I am anxiously waiting for your response cherrs
