# The constants used in all functions are defined.

using JuMP, Ipopt, Clp
# properties of air
Cp_air = 1.005              # heat capacity of air kJ/kg.K
rho_air = 1.204             # density of air kg/m^3
# costs of various processes [refer to Table 1]
C_heat = 0.0170;             # unit hourly cost of heating
C_mechcool = 0.0333;         # unit hourly cost of mechanical cooling
C_mist = 0.00475;            # unit hourly cost of misting
C_dehumid = 0.0752;          # unit hourly cost of dehumidifying
C_CO2 = .339;                # unit hourly cost of injecting CO2
C_AXH = 0.0137;              # unit hourly cost of 
# facility parameters
U = 2                       # conductance of farm wall and roof [W/m^2/K]
SA = 2560                   # surface area of farm walls and roof m^2
V_farm = 9600;              # volume of farm m^3
q_lights = 992;              # [kW] heat being generated by lights
q_transpiration = 120.7;     # [g/s] loss of water into air
q_respiration = 2.73;        # [g/s] used by plants  

function temp70LP(RH_outside, TF)
    C_max = 18.03               # water holding capacity if air at 70 F [g/m^3] 
    
    logspaced = logspace(0, 2, 100);
    # various air exchange rates [m^3/s]
    V_dot = zeros(Array{Float64}(100)) 
    for i = 1:100
        V_dot[i] = logspaced[i]*2.675
    end 

    p = zeros(Array{Float64}(100))
    heatvar = zeros(Array{Float64}(100))
    mechcoolvar = zeros(Array{Float64}(100))
    evapcoolvar = zeros(Array{Float64}(100))
    dehumidvar = zeros(Array{Float64}(100))
    humidvar = zeros(Array{Float64}(100))
    CO2var = zeros(Array{Float64}(100))    
    RHvar = zeros(Array{Float64}(100))   

    for i = 1:100

            m = Model(solver = ClpSolver()) #solver = GurobiSolver(OutputFlag=

            # setpoint parameters
            T = (70 + 459.67)*5/9;                  # setpoint temperature, 70 [F] in [K]
            @variable(m, 50 <= RH <= 80)  # solve for the optimal relative humidity in the chamber
            @variable(m, CO2 >= 0.7368 ) # CO2 maintained at a minimum of 0.7368 [g/m^3], which is equal to 400 ppm
        
            # outdoor weather parameters
            T_outside = (TF + 459.67)*5/9;           # temperature, 80 [F] in [K] 299.817
            CO2_outside = 0.7368;           # carbon dioxide concentration in [g CO2/m^3 air]
            TD = (T_outside-273.15) - ((100 - RH_outside)/5) + 273.15;
            RH_newtemp = 5*((TD - 273.15) - (T-273.15)) + 100;
            
            # HVAC equipment variables
            @variable(m,q_heat >= 0) # quantity of heating
            @variable(m,q_mechcool >= 0) # quantity of cooling
            @variable(m,q_CO2 >= 0) # quantity of heating
            @variable(m,q_mist >= 0) # quantity of humidification
            @variable(m,q_dehumid >= 0) # quantity of dehumidification
            
            # energy and mass balances
            @constraint(m, rho_air*Cp_air*T*V_dot[i] == rho_air*Cp_air*T_outside*V_dot[i] + q_lights + q_heat - q_mechcool  - 2.257*(q_mist + q_transpiration) + 19.87*q_CO2 + U*SA*(T_outside - T)/1000) # SS energy balance 
            @constraint(m, C_max*RH/100*V_dot[i] == C_max*RH_newtemp/100*V_dot[i] + q_transpiration + q_mist - q_dehumid) # Moisture balance     
            @constraint(m, CO2*V_dot[i] == CO2_outside*V_dot[i] - q_respiration + q_CO2)                    # SS CO2 balance

            @objective(m, Min, q_heat*C_heat + q_mechcool*C_mechcool + q_CO2*C_CO2 + V_dot[i]*C_AXH + q_mist*C_mist + q_dehumid*C_dehumid)
            solve(m)

            p[i] = getvalue(getobjective(m))   
            heatvar[i] = getvalue(q_heat)
            mechcoolvar[i] = getvalue(q_mechcool)
            dehumidvar[i] = getvalue(q_dehumid)
            humidvar[i] = getvalue(q_mist)
            CO2var[i] = getvalue(q_CO2)    
            RHvar[i] = getvalue(RH)  

    end    

    # determine the minimum-cost ventilation rate
    ind = findmin(p)
    optcost = ind[1]
    optheat= heatvar[ind[2]]
    optmechcool= mechcoolvar[ind[2]]
    optmist= humidvar[ind[2]]
    optdehumid = dehumidvar[ind[2]]
    optCO2= CO2var[ind[2]] 
    optRH= RHvar[ind[2]]
    optVent = V_dot[ind[2]]  
    return(p, V_dot, heatvar, mechcoolvar, dehumidvar, humidvar, CO2var, RHvar, optcost, optheat, optmechcool, optmist, optdehumid, optCO2, optRH, optVent)
end;

function tempRangeNLP(RH_outside, TF, CO2_L, MAXTF)
    
    m = Model(solver = IpoptSolver(print_level=0))
    T_outside = (TF + 459.67)*5/9;# temperature, 80 [F] in [K] 299.817
    CO2_outside = 0.7368;         # carbon dioxide concentration in [g CO2/m^3 air]

    # setpoint variables
    @variable(m, (70 + 459.67)*5/9 <= T <= (MAXTF + 459.67)*5/9)
    @variable(m, 50 <= RH <= 80)  # solve for the optimal relative humidity in the chamber
    @variable(m, CO2 >= CO2_L )  # CO2 maintained at a minimum of 0.79 [g/m^3], which is equal to 400 ppm

    # HVAC equipment variables
    @variable(m, 2.675 <= V_dot <= 2.65*100) # quantity of heating
    @variable(m,q_heat >= 0) # quantity of heating
    @variable(m,q_mechcool >= 0) # quantity of cooling
    @variable(m,q_CO2 >= 0) # quantity of heating
    @variable(m,q_mist >= 0) # quantity of humidification
    @variable(m,q_RHdown >= 0) # quantity of dehumidification

    # Water holding capacity variables and constraints
    @variable(m, TD)
    @variable(m, RH_newtemp)
    @variable(m, C_max)
    @constraint(m, TD == (T_outside-273.15) - ((100 - RH_outside)/5))
    @constraint(m, RH_newtemp == 5*((TD) - (T-273.15)) + 100)
    @constraint(m,C_max == 0.7278*(T*9/5 - 459.67) - 32.189)

    # mass and energy balances
    @constraint(m, rho_air*Cp_air*T*V_dot == rho_air*Cp_air*T_outside*V_dot + q_lights + q_heat - q_mechcool  - 2.257*(q_mist + q_transpiration) + 19.87*q_CO2 + U*SA*(T_outside - T)/1000) # Energy balance 
    @NLconstraint(m, C_max*RH/100*V_dot == C_max*RH_newtemp/100*V_dot + q_transpiration + q_mist - q_RHdown) # Moisture balance     
    @constraint(m, CO2*V_dot == CO2_outside*V_dot - q_respiration + q_CO2)                    # CO2 balance

    @NLobjective(m, Min, q_heat*C_heat + q_mechcool*C_mechcool + q_CO2*C_CO2 + V_dot*C_AXH + q_mist*C_mist + q_RHdown*C_dehumid)
    solve(m)
    
    optcost = getvalue(q_heat*C_heat + q_mechcool*C_mechcool + q_CO2*C_CO2 + V_dot*C_AXH + q_mist*C_mist + q_RHdown*C_dehumid)
    optheat = getvalue(q_heat)
    optmechcool = getvalue(q_mechcool)
    optdehumid = getvalue(q_RHdown)
    optmist = getvalue(q_mist)
    optCO2 = getvalue(q_CO2)    
    optRH = getvalue(RH)  
    optT = getvalue(T)
    optVent = getvalue(V_dot)
    return(optcost,optheat, optmechcool, optdehumid, optmist, optCO2, optRH, optT, optVent) 
end;

raw = readcsv("726410TYA.csv");
T_TMY = raw[:,1];     # Temperature in [C]
RH_TMY = raw[:,3];    # Relative humidity in [%]

function temp70LPTMY()
# variables to hold optimum solutions at each hour of year
optcost_TMY = zeros(Array{Float64}(8760))
optheat_TMY = zeros(Array{Float64}(8760))
optmechcool_TMY = zeros(Array{Float64}(8760))
optmist_TMY = zeros(Array{Float64}(8760))
optdehumid_TMY = zeros(Array{Float64}(8760))
optCO2_TMY = zeros(Array{Float64}(8760))
optAXH_TMY = zeros(Array{Float64}(8760))
hourOfYear = zeros(Array{Float64}(8760))
optRH_TMY = zeros(Array{Float64}(8760))

for j = 2:8761

    # outdoor weather parameters
    T_outside = T_TMY[j] + 273.15;               # temperature in [K]
    CO2_outside = 0.7368;                        # carbon dioxide concentration in [g/m^3]
    RH_outside = RH_TMY[j];                      # relative humidity [%]
    TF = T_outside*9/5 - 459.67                  # temperature in [F]
    
    # call temp70LP
    (p, V_dot, heatvar, mechcoolvar, dehumidvar, humidvar, CO2var, RHvar, optcost, optheat, optmechcool, optmist, optdehumid, optCO2, optRH, optVent) = temp70LP(RH_outside, TF);
    
    # store optimum solutions
    optcost_TMY[j-1] = optcost
    hourOfYear[j-1] = j-1    
    optheat_TMY[j-1] = optheat
    optmechcool_TMY[j-1] = optmechcool
    optmist_TMY[j-1] = optmist
    optdehumid_TMY[j-1] = optdehumid
    optCO2_TMY[j-1]= optCO2
    optRH_TMY[j-1]= optRH
    optAXH_TMY[j-1] = optVent
    
end   
return(optcost_TMY, hourOfYear, optheat_TMY, optmechcool_TMY, optmist_TMY, optdehumid_TMY, optCO2_TMY, optRH_TMY, optAXH_TMY)    
end;

function tempRangeNLPTMY(CO2_L, MAX_TF)
# variables to hold optimum solutions at each hour of year
optcost_TMY = zeros(Array{Float64}(8760))
optheat_TMY = zeros(Array{Float64}(8760))
optmechcool_TMY = zeros(Array{Float64}(8760))
optmist_TMY = zeros(Array{Float64}(8760))
optdehumid_TMY = zeros(Array{Float64}(8760))
optCO2_TMY = zeros(Array{Float64}(8760))
optAXH_TMY = zeros(Array{Float64}(8760))
hourOfYear = zeros(Array{Float64}(8760))
optRH_TMY = zeros(Array{Float64}(8760))
optT_TMY = zeros(Array{Float64}(8760))
optVent_TMY = zeros(Array{Float64}(8760))
    
for j = 2:8761

    # outdoor weather parameters    
    T_outside = (T_TMY[j] + 273.15)*1;           # temperature in [F]
    CO2_outside = 0.7368;                        # carbon dioxide concentration in [ppm]
    RH_outside = RH_TMY[j];                      # relative humidity, (vapor pressure/saturation vapor pressure)*100 [%]

    TF = T_outside*9/5 - 459.67   
    
    # call tempRangeNLP
    (optcost,optheat, optmechcool, optdehumid, optmist, optCO2, optRH, optT, optVent)= tempRangeNLP(RH_outside, TF, CO2_L, MAX_TF);

    # store optimum solutions
    optcost_TMY[j-1] = optcost
    hourOfYear[j-1] = j-1    
    optheat_TMY[j-1] = optheat
    optmechcool_TMY[j-1] = optmechcool
    optmist_TMY[j-1] = optmist
    optdehumid_TMY[j-1] = optdehumid
    optCO2_TMY[j-1]= optCO2 
    optRH_TMY[j-1]= optRH 
    optAXH_TMY[j-1] = optVent 
    optT_TMY[j-1] = optT
    optVent_TMY[j-1] = optVent
    
end
return(optcost_TMY, hourOfYear, optheat_TMY, optmechcool_TMY, optmist_TMY, optdehumid_TMY, optCO2_TMY, optRH_TMY, optAXH_TMY, optT_TMY, optVent_TMY)       
end;

using PyPlot
# temp70LP at 50% RH, 80 F
(p5080, Vdot5080, heatvar5080, mechcoolvar5080, dehumidvar5080, humidvar5080, CO2var5080, RHvar5080, optcost5080, optheat5080, optmechcool5080, optmist5080, optdehumid5080, optCO25080, optRH5080, optVent5080) = temp70LP(50, 80);
plot(Vdot5080, p5080, label = "80 [F]")
# temp70LP at 50% RH, 60 F
(p5060, Vdot5060, heatvar5060, mechcoolvar5060, dehumidvar5060, humidvar5060, CO2var5060, RHvar5060, optcost5060, optheat5060, optmechcool5060, optmist5060, optdehumid5060, optCO25060, optRH5060, optVent5060) = temp70LP(50, 60);
plot(Vdot5060, p5060, label = "60 [F]")
# temp70LP at 50% RH, 0 F
(p500, Vdot500, heatvar500, mechcoolvar500, dehumidvar500, humidvar500, CO2var500, RHvar500, optcost500, optheat500, optmechcool500, optmist500, optdehumid500, optCO2500, optRH500, optVent500) = temp70LP(50, 0);
plot(Vdot500, p500, label = "0 [F]")
xlabel("Ventilation rate [m^3/s]")
ylabel("Hourly HVAC cost [dollars]")
ylim(0,100)
legend(loc="lower right");

using PyPlot
# caclulate costs at each hour of year using LP and NLP methods 
(optcost_TMY70, hourOfYear70, optheat_TMY70, optmechcool_TMY, optmist_TMY, optdehumid_TMY, optCO2_TMY, optRH_TMY, optAXH_TMY) =temp70LPTMY();
(optcost_TMYR, hourOfYearR, optheat_TMYR, optmechcool_TMYR, optmist_TMYR, optdehumid_TMYR, optCO2_TMYR, optRH_TMYR, optAXH_TMYR, optT_TMYR, optVentR) = tempRangeNLPTMY(0.7368, 70); 

# plot Typical Meteorlogical year data for Madison WI
plot(hourOfYear70, RH_TMY[2:8761], label = "RH_outside [%]")
plot(hourOfYear70, (T_TMY[2:8761]*(9/5)+32), label = "T_outside [F]")
xlabel("Hour of year")
ylabel("T/RH")
legend(loc = "lower right");

plot(hourOfYear70, optcost_TMY70, label = "LP")
plot(hourOfYearR, optcost_TMYR, label = "NLP")
xlabel("Hour of year")
ylabel("Hourly HVAC cost")
legend(loc = "right");

# call tempRangeNLPTMY to calculate annual costs vs. at various max indoor air temps and CO2 concentrations
(optcost_TMYR80, hourOfYearR80, optheat_TMYR80, optmechcool_TMYR80, optmist_TMYR80, optdehumid_TMYR80, optCO2_TMYR80, optRH_TMYR80, optAXH_TMYR80, optT_TMYR80, optVentR80) = tempRangeNLPTMY(0.7368, 80); 
(optcost_TMYR80080, hourOfYearR80080, optheat_TMYR80080, optmechcool_TMYR80080, optmist_TMYR80080, optdehumid_TMYR80080, optCO2_TMYR80080, optRH_TMYR80080, optAXH_TMYR80080, optT_TMYR80080, optVentR80080) = tempRangeNLPTMY(2*0.7368, 80); 
(optcost_TMYR120080, hourOfYearR120080, optheat_TMYR120080, optmechcool_TMYR120080, optmist_TMYR120080, optdehumid_TMYR120080, optCO2_TMYR120080, optRH_TMYR120080, optAXH_TMYR120080, optT_TMYR120080, optVentR120080) = tempRangeNLPTMY(3*0.7368, 80)

(optcost_TMYR75, hourOfYearR75, optheat_TMYR75, optmechcool_TMYR75, optmist_TMYR75, optdehumid_TMYR75, optCO2_TMYR75, optRH_TMYR75, optAXH_TMYR75, optT_TMYR75, optVentR75) = tempRangeNLPTMY(0.7368, 75); 
(optcost_TMYR80075, hourOfYearR80075, optheat_TMYR80075, optmechcool_TMYR80075, optmist_TMYR80075, optdehumid_TMYR80075, optCO2_TMYR80075, optRH_TMYR80075, optAXH_TMYR80075, optT_TMYR80075, optVentR80075) = tempRangeNLPTMY(2*0.7368, 75); 
(optcost_TMYR120075, hourOfYearR120075, optheat_TMYR120075, optmechcool_TMYR120075, optmist_TMYR120075, optdehumid_TMYR120075, optCO2_TMYR120075, optRH_TMYR120075, optAXH_TMYR120075, optT_TMYR120075, optVentR120075) = tempRangeNLPTMY(3*0.7368, 75);

(optcost_TMYR70, hourOfYearR70, optheat_TMYR70, optmechcool_TMYR70, optmist_TMYR70, optdehumid_TMYR70, optCO2_TMYR70, optRH_TMYR70, optAXH_TMYR70, optT_TMYR70, optVentR70) = tempRangeNLPTMY(0.7368, 70); 
(optcost_TMYR80070, hourOfYearR80070, optheat_TMYR80070, optmechcool_TMYR80070, optmist_TMYR80070, optdehumid_TMYR80070, optCO2_TMYR80070, optRH_TMYR80070, optAXH_TMYR80070, optT_TMYR80070, optVentR80070) = tempRangeNLPTMY(2*0.7368, 70); 
(optcost_TMYR120070, hourOfYearR120070, optheat_TMYR120070, optmechcool_TMYR120070, optmist_TMYR120070, optdehumid_TMYR120070, optCO2_TMYR120070, optRH_TMYR120070, optAXH_TMYR120070, optT_TMYR120070, optVentR120070) = tempRangeNLPTMY(3*0.7368, 70);

plot(T_TMY[2:8761]*(9/5) + 32, optcost_TMYR80, "r.", label = " 80 [F]")
plot(T_TMY[2:8761]*(9/5) + 32, optcost_TMYR75, "b.", label = " 75 [F]")
plot(T_TMY[2:8761]*(9/5) + 32, optcost_TMYR70, "k.", label = " 70 [F]")
xlabel("Air temperature [F]")
ylabel("Hourly HVAC cost [dollars]");
legend(loc="upper left");

plot(T_TMY[2:8761]*(9/5) + 32, optcost_TMYR80, "r.", label = "[400 ppm CO2]")
plot(T_TMY[2:8761]*(9/5) + 32, optcost_TMYR80080, "b.", label = "[800 ppm CO2]")
plot(T_TMY[2:8761]*(9/5) + 32, optcost_TMYR120080, "k.", label = "[1200 ppm CO2]")
xlabel("Air temperature [F]")
ylabel("Hourly HVAC cost [dollars]");
legend(loc="upper left");

# Code used to generate data for Table 3.

# 80 [F]
println(optcost5080)
println(optVent5080)
println(optheat5080)
println(optmechcool5080)
println(optmist5080)
println(optdehumid5080)
println(optCO25080)

# 60 [F]
println(optcost5060)
println(optVent5060)
println(optheat5060)
println(optmechcool5060)
println(optmist5060)
println(optdehumid5060)
println(optCO25060)

# 0 [F]
println(optcost500)
println(optVent500)
println(optheat500)
println(optmechcool500)
println(optmist500)
println(optdehumid500)
println(optCO2500)

# compare the results from the LP and NLP functions (Table 4)

# compare results at T_outside = 60 F, RH_outside = 50 %
(p, V_dot, heatvar, mechcoolvar, dehumidvar, humidvar, CO2var, RHvar, optcost, optheat, optmechcool, optmist, optdehumid, optCO2, optRH, optVent)= temp70LP(50, 60)
println("At T_outside = 60 F, RH_outside = 50 %, optimal hourly cost as found with temp70LP: ")
println(optcost)
(optcost,optheat, optmechcool, optdehumid, optmist, optCO2, optRH, optT, optVent) = tempRangeNLP(50, 60, 0.7368, 70);
println("At T_outside = 60 F, RH_outside = 50 %, Optimal hourly cost as found with tempRangeNLP: ")
println(optcost)    
println()   

# compare results at T_outside = 80 F, RH_outside = 50 %
(p, V_dot, heatvar, mechcoolvar, dehumidvar, humidvar, CO2var, RHvar, optcost, optheat, optmechcool, optmist, optdehumid, optCO2, optRH, optVent)= temp70LP(50, 80)
println("At T_outside = 80 F, RH_outside = 50 %, optimal hourly cost as found with temp70LP: ")
println(optcost)
(optcost,optheat, optmechcool, optdehumid, optmist, optCO2, optRH, optT, optVent) = tempRangeNLP(50, 80, 0.7368, 70);
println("At T_outside = 80 F, RH_outside = 50 %, Optimal hourly cost as found with tempRangeNLP: ")
println(optcost)  
println()   

# compare results at T_outside =0 F, RH_outside = 50 %
(p, V_dot, heatvar, mechcoolvar, dehumidvar, humidvar, CO2var, RHvar, optcost, optheat, optmechcool, optmist, optdehumid, optCO2, optRH, optVent)= temp70LP(50, 0)
println("At T_outside = 0 F, RH_outside = 50 %, optimal hourly cost as found with temp70LP: ")
println(optcost)
(optcost,optheat, optmechcool, optdehumid, optmist, optCO2, optRH, optT, optVent) = tempRangeNLP(50, 0, 0.7368, 70);
println("At T_outside = 0 F, RH_outside = 50 %, Optimal hourly cost as found with tempRangeNLP: ")
println(optcost)   


# print data for Table 5
annualcost40080 = sum(optcost_TMYR80)
annualcost80080 = sum(optcost_TMYR80080)
annualcost120080 = sum(optcost_TMYR120080)
println(annualcost40080)
println(annualcost80080)
println(annualcost120080)
println(sum(optcost_TMYR75))
println(sum(optcost_TMYR80075))
println(sum(optcost_TMYR120075))
println(sum(optcost_TMYR70))
println(sum(optcost_TMYR80070))
println(sum(optcost_TMYR120070))

# print data for Table 6.
println("Breakdown of annual costs when max indoor temp is 80 and min CO2 is 400 ppm")
println(sum(optheat_TMYR80)*C_heat/sum(optcost_TMYR80))
println(sum(optmechcool_TMYR80)*C_mechcool/sum(optcost_TMYR80))
println(sum(optmist_TMYR80)*C_mist/sum(optcost_TMYR80))
println(sum(optdehumid_TMYR80)*C_dehumid/sum(optcost_TMYR80))
println(sum(optCO2_TMYR80)*C_CO2/sum(optcost_TMYR80))
println(sum(optVentR80)*C_AXH/sum(optcost_TMYR80))
println()
println("Breakdown of annual costs when max indoor temp is 70 and min CO2 is 400 ppm")
println(sum(optheat_TMYR70)*C_heat/sum(optcost_TMYR70))
println(sum(optmechcool_TMYR70)*C_mechcool/sum(optcost_TMYR70))
println(sum(optmist_TMYR70)*C_mist/sum(optcost_TMYR70))
println(sum(optdehumid_TMYR70)*C_dehumid/sum(optcost_TMYR70))
println(sum(optCO2_TMYR70)*C_CO2/sum(optcost_TMYR70))
println(sum(optVentR70)*C_AXH/sum(optcost_TMYR70))
println()
println("Breakdown of annual costs when max indoor temp is 75 and min CO2 is 400 ppm")
println(sum(optheat_TMYR75)*C_heat/sum(optcost_TMYR75))
println(sum(optmechcool_TMYR75)*C_mechcool/sum(optcost_TMYR75))
println(sum(optmist_TMYR75)*C_mist/sum(optcost_TMYR75))
println(sum(optdehumid_TMYR75)*C_dehumid/sum(optcost_TMYR75))
println(sum(optCO2_TMYR75)*C_CO2/sum(optcost_TMYR75))
println(sum(optVentR75)*C_AXH/sum(optcost_TMYR75))

# sum of hourly costs shown in figure 7
println("The estimated annual cost of using the LP is: ")
println(sum(optcost_TMY70))
println("The estimated annual cost of using the NLP is: ")
println(sum(optcost_TMYR))