""" This is the code for optimizing NajadAbad Pump Station by Pedram Jazayeri
    with 1 houre periods 
    in a day
    without electricity tarriff"""
""" v3.1.1 = delete Heights of twins and Suction revervoirs
             """

###############################################################################################
import matplotlib.pyplot as plt
from docplex.mp.model import Model
import pandas as pd
#import math
import numpy as np


###############################################################################################
#VARIABELS 
Res_10000Area=(33.57-2*0.35)*(67.27-3*0.35)    
Res_20000Area=(54.65-2*0.35)*(78.55-3*0.35) 

Path2='D:/BMN/1401/04/vilashahr.xlsx'
Names2=['Name1', 'Name2', 'Value', 'Unit', 'Set', 'Date']
InFlowRead=pd.read_excel(Path2, sheet_name='InFlow Discharge', header=None, names=Names2)
InFlow=list(InFlowRead['Value'][0:288])
NumberOfDaysForAvailableData=int(len(InFlow)/288) 

InFlow_1Hr=list(range(24*NumberOfDaysForAvailableData))
InFlow_2Hr=list(range(12*NumberOfDaysForAvailableData))
InFlow_3Hr=list(range(8*NumberOfDaysForAvailableData))
InFlow_4Hr=list(range(6*NumberOfDaysForAvailableData))
InFlow_6Hr=list(range(4*NumberOfDaysForAvailableData))


TimeFrame=InFlow_4Hr          #SELECT DESIRED TIME FRAMES(1/2/3/4/5/6 HOURS)
num=int(len(TimeFrame))
print('')
print(f"You selected time period='{int(24/len(TimeFrame))} hours'. ")
print('')
period=int(24/num)

for i in range(num):
    c=0
    for j in range(12*period*i, 12*period*(i+1)):
        c=c+InFlow[j]
    TimeFrame[i]=c/(12*period)
      
Pump=list(range(num+1))
Pump[0]=np.nan

Net= list(range(num+1))
Net[0]=np.nan

H10= list(range(num+1))
H10[0]=3.89713526

H20= list(range(num+1))
H20[0]=3.84223247

Multipliers=[0.88, 0.63, 0.45, 0.36, 0.35, 0.44, 0.59, 0.79,\
             0.98, 1.13, 1.25, 1.30, 1.32, 1.32, 1.30, 1.30,\
             1.30, 1.30, 1.30, 1.28, 1.22, 1.15, 1.07, 0.99]
Sum_Multipliers=sum(Multipliers)

Daily_ConsumePredicted=34750  #cubic meters per day


###############################################################################################
#CREATE OPTIMIZATION MODEL
m = Model(name='Pump Station Optimization')
NumberOfTurnedOnPumps=list(range(num))
H10000=list(range(num))
H20000=list(range(num))
NetUsage=list(range(num))
NetConstraint=list(range(num))
H10000Constraint=list(range(num))
H20000Constraint=list(range(num))

for i in range(num):
    NumberOfTurnedOnPumps[i]=m.integer_var(name='NumberOfTurnedOnPumps'+f'{i}', lb=0, ub=3)
    H10000[i]=m.continuous_var(name='H10000'+f'{i}', lb=0.5 , ub=4.5)
    H20000[i]=m.continuous_var(name='H20000'+f'{i}', lb=0.5 , ub=4.5)
    NetUsage[i]=m.continuous_var(name='NetUsage'+f'{i}')


###############################################################################################
#CONSTRAINTS
H10000Constraint[0]=m.add_constraint((H10000[0]*Res_10000Area)==(H10[i]*Res_10000Area)+\
                                   NumberOfTurnedOnPumps[0]*100*0.3-0.35*NetUsage[0])
H20000Constraint[0]=m.add_constraint((H20000[0]*Res_20000Area)==(H20[i]*Res_20000Area)+\
                                   0.9*300*0.3-NumberOfTurnedOnPumps[0]*100*0.3-0.5*NetUsage[0])
   
if num==24:
    NetConstraint[0]=m.add_constraint( NetUsage[0]==( (Daily_ConsumePredicted/Sum_Multipliers) *  (Multipliers[0])))
    convert_index=3.6
elif num==12:
    NetConstraint[0]=m.add_constraint( NetUsage[0]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[0]+Multipliers[1])/2)*2))
    convert_index=3.6*2
elif num==8:
    NetConstraint[0]=m.add_constraint( NetUsage[0]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[0]+Multipliers[1]+Multipliers[2])/3)*3))
    convert_index=3.6*3
elif num==6:
    NetConstraint[0]=m.add_constraint( NetUsage[0]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[0]+Multipliers[1]+Multipliers[2]+Multipliers[3])/4)*4))
    convert_index=3.6*4
elif num==4:
    NetConstraint[0]=m.add_constraint( NetUsage[0]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[0]+Multipliers[1]+Multipliers[2]+Multipliers[3]+Multipliers[4]+Multipliers[5])/6)*6))
    convert_index=3.6*6


for i in range(num-1):
    H10000Constraint[i+1]=m.add_constraint((H10000[i+1]*Res_10000Area)==(H10000[i]*Res_10000Area)+\
                                         NumberOfTurnedOnPumps[i+1]*100*convert_index-0.35*NetUsage[i+1])
    H20000Constraint[i+1]=m.add_constraint((H20000[i+1]*Res_20000Area)==(H20000[i]*Res_20000Area)+\
                                             0.9*300*convert_index-NumberOfTurnedOnPumps[i+1]*100*convert_index-0.5*NetUsage[i+1])    
        
    if   num==24:
        NetConstraint[i+1]=m.add_constraint( NetUsage[i+1]==( (Daily_ConsumePredicted/Sum_Multipliers) *  (Multipliers[i+1])))       

    elif num==12:
        NetConstraint[i+1]=m.add_constraint( NetUsage[i+1]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[2*(i+1)]+Multipliers[2*(i+1)+1]) /2) *2) )  
    
    elif num==8:
        NetConstraint[i+1]=m.add_constraint( NetUsage[i+1]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[3*(i+1)]+Multipliers[3*(i+1)+1]+Multipliers[3*(i+1)+2])/3)*3))  
    
    elif num==6:
        NetConstraint[i+1]=m.add_constraint( NetUsage[i+1]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[4*(i+1)]+Multipliers[4*(i+1)+1]+Multipliers[4*(i+1)+2]+Multipliers[4*(i+1)+3])/4)*4))
    
    elif num==4:
        NetConstraint[i+1]=m.add_constraint( NetUsage[i+1]==( (Daily_ConsumePredicted/Sum_Multipliers) * ((Multipliers[6*(i+1)]+Multipliers[6*(i+1)+1]+Multipliers[6*(i+1)+2]+Multipliers[6*(i+1)+3]+Multipliers[6*(i+1)+4]+Multipliers[6*(i+1)+5])/6)*6))    


###############################################################################################
#OBJECTIVE FUNCTION
if num==24:     
    m.minimize( NumberOfTurnedOnPumps[0] +NumberOfTurnedOnPumps[1] +NumberOfTurnedOnPumps[2] +NumberOfTurnedOnPumps[3]\
               +NumberOfTurnedOnPumps[4] +NumberOfTurnedOnPumps[5] +NumberOfTurnedOnPumps[6] +NumberOfTurnedOnPumps[7]\
               +NumberOfTurnedOnPumps[8] +NumberOfTurnedOnPumps[9] +NumberOfTurnedOnPumps[10]+NumberOfTurnedOnPumps[11]\
               +NumberOfTurnedOnPumps[12]+NumberOfTurnedOnPumps[13]+NumberOfTurnedOnPumps[14]+NumberOfTurnedOnPumps[15]\
               +NumberOfTurnedOnPumps[16]+NumberOfTurnedOnPumps[17]+NumberOfTurnedOnPumps[18]+NumberOfTurnedOnPumps[19]\
               +NumberOfTurnedOnPumps[20]+NumberOfTurnedOnPumps[21]+NumberOfTurnedOnPumps[22]+NumberOfTurnedOnPumps[23] )

elif num==12:    
    m.minimize( NumberOfTurnedOnPumps[0] +NumberOfTurnedOnPumps[1] +NumberOfTurnedOnPumps[2] +NumberOfTurnedOnPumps[3]\
               +NumberOfTurnedOnPumps[4] +NumberOfTurnedOnPumps[5] +NumberOfTurnedOnPumps[6] +NumberOfTurnedOnPumps[7]\
               +NumberOfTurnedOnPumps[8] +NumberOfTurnedOnPumps[9] +NumberOfTurnedOnPumps[10]+NumberOfTurnedOnPumps[11] )

elif num==8:    
    m.minimize( NumberOfTurnedOnPumps[0] +NumberOfTurnedOnPumps[1] +NumberOfTurnedOnPumps[2] +NumberOfTurnedOnPumps[3]\
               +NumberOfTurnedOnPumps[4] +NumberOfTurnedOnPumps[5] +NumberOfTurnedOnPumps[6] +NumberOfTurnedOnPumps[7] )
        
elif num==6:    
    m.minimize( NumberOfTurnedOnPumps[0] +NumberOfTurnedOnPumps[1] +NumberOfTurnedOnPumps[2] +NumberOfTurnedOnPumps[3]\
               +NumberOfTurnedOnPumps[4] +NumberOfTurnedOnPumps[5] )

elif num==4:    
    m.minimize( NumberOfTurnedOnPumps[0] +NumberOfTurnedOnPumps[1] +NumberOfTurnedOnPumps[2] +NumberOfTurnedOnPumps[3] )


###############################################################################################
#UPDATE PARAMETERS
sol=m.solve()
for i in range(num):
    H10[i+1]= sol.get_value(H10000[i])
    H20[i+1]= sol.get_value(H20000[i])
    Pump[i+1]=sol.get_value(NumberOfTurnedOnPumps[i])
    Net[i+1]= sol.get_value(NetUsage[i])
sol.display()


###############################################################################################  
#PLOT    
plt.plot(H20, label='Res20000')
plt.legend()
plt.show()

plt.plot(H10, label='Res10000')
plt.legend()
plt.show()

plt.plot(Multipliers, label='Hourly Multipliers', color='red')
plt.legend()
plt.show()

plt.plot(Net, label='Network Usage', color='red')
plt.legend()
plt.show()

plt.bar(list(range(1,num+1)), Pump[1:num+1], label='Number of Turnd On Pumps', color='c')
plt.legend()
plt.show()


###############################################################################################