# 1. Introduction¶

Indoor agriculture is a rapidly growing industry driven by demand for local produce, cannabis legalization, and incentives to reduce water consumption. However, growing food indoors with artificial light is energy intensive. The heat generated by the lights also needs to be exhausted [1]. In addition, humidity needs to be controlled to prevent pathogen growth, and allow transpiration and nutrient uptake. Plant growth can be improved by increasing the CO2 concentration in the growing space [2].

Recent suvey results report that the two greatest challenges facing growers are (1) managing and (2) predicting/stabilizing operating costs [3]. While current guidelines give good rules-of-thumb for growing plants indoors, they do not ensure cost-optimal HVAC control. The primary reason is that there are multiple ways to achieve the same environmental conditions. For example, ‘cooling’ can be achieved by simply ventilating if the outside air is cooler than the setpoint, by evaporative cooling if the wet-bulb temperature of the outside air is cooler than the setpoint, or by using refrigeration in any scenario. Ventilation is the cheapest, followed by evaporative cooling, and mechanical cooling (i.e. equipment using vapor-compression refrigeration) is the most expensive. Humidity control is a similar problem. Ventilation is the cheapest way to exhaust moisture to the atmosphere, but if a facility is heating, cooling, or CO2-enriching the indoor air, or the outside air is excessively humid, then ventilation can be counterproductive. Therefore, some level of mechanical dehumidification makes sense in many cases.

An optimization program is developed herein to solve this problem. The program determines the cost-optimal HVAC system operation for maintaining the indoor growing space within acceptable ranges for temperature, relative humidity, and CO2 concentration given the outdoor weather conditions. That is, of the possible combinations of HVAC equipment operation that could bring the indoor farm's conditions within their desired ranges, this program determines the cheapest combination.

This report will explain the mathematical model developed, including the facility characteristics and assumptions, the environmental requirements of the crop being grown (temperature, relative humidity, and CO2 concentration), the variables and constants used to represent each piece of HVAC equipment, and laws of conservation of mass and energy that constrain the solution. The objective function is formulated to minimize the the operating costs of the HVAC equipment. The solution to the optimization problem is implemented both as a linear program (LP) and as a non-linear program (NLP). The LP requires fixing some of the variables, but is useful to give physical intuition for the problem and to demonstrate that the NLP is correctly solving the optimization problem. NLPs are susceptible to finding local solutions. The programs developed are then used, along with Typical Meteorlogical Year weather data [4] for Madison, Wisconsin to calculate annual operating costs for the facility considered at several levels of CO2-enrichment and maximum indoor temperature.

[1] Storey, N. (2016). Operating Expenses of Indoor Growing | Let's Talk Indoor Farming! https://www.youtube.com/watch?v=Qwy8jrZ1eNE

[2] Sabeh, N. (2014). Tomato Greenhouse Roadmap. Chap. 4: Climate Control Systems.

[3] Agrilyst. State of Indoor Farming. White paper. http://stateofindoorfarming.agrilyst.com/#cta

[4] National Renewable Energy Laboratory. National Solar Radiation Data Base. TMY 3. http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/tmy3/.

# 2. Mathematical model¶

## 2a. Facility characteristics and key assumptions¶

The indoor farm considered has the characteristics presented in Table 1. The justification for this set of assumptions is contained in [5]:

Table 1. Facility characteristics, including generation of heat and moisture and consumption of CO2.

Characteristic Constant/value
Dimensions (l x w x h) 40 [m] x 40 [m] x 6 [m]
Wall/roof conductivity $U$ = 2 [$W/m^{2} \cdot ^{\circ}K$]
Heat generated by lighting $\dot{q}_{lights}$ = 992 $[kW]$
Moisture generated evaporation and evapotranspiration $\dot{q}_{transpiration}$ = 120.7 [$g_{water}/s$]
CO2 consumed by the plants $\dot{q}_{photosynthesis}$ = 2.73 [$g_{CO2}/s$]

The following assumptions are also important to the model:

• All processes occur in steady state in one hour intervals.
• The air is well-mixed (i.e. uniform). This means the properties of air inside the farm and air leaving the farm are equal.

## 2b. Grow space requirements¶

The conditions required inside the farm, along with the variable names used to represent them, are shown in figure 1.

Figure 1. Grow space requirements. ($T$) is air temperature, ($RH$) is relative humidity, and ($CO2$) is CO2 concentration. In the LP solution, T is set to 70 F.

## 2c. Representation of HVAC equipment¶

The equipment the indoor farm has, along with the variables representing their operation, are presented in Table 2.

Table 2. Indoor farm HVAC equipment, variable representation, and cost.

Equipment Variable Cost constant Unit hourly cost
Ventilation $\dot{V}$ [$m^{3}$/s] $C_{vent}$ 0.0137 $[\$/((m^{3}/s)h)]$Gas heater$\dot{q}_{heat}[kW]C_{heat}$0.0171$[\$/kWh]$
Water misting $\dot{q}_{mist}$ [$g_{water}/s$] $C_{mist}$ 0.00475 $[\$/(g_{water}/s)h)]$Mechanical cooling$\dot{q}_{mechcool}$[$kW$]$C_{mechcool}$0.0333$[\$/kWh]$
Dehumidification $\dot{q}_{dehumid}$ [$g_{water}/s$] $C_{dehumid}$ 0.0752 $[\$/(g_{water}/s)h)]$CO2 generation$\dot{q}_{CO2}$[$g/s$]$C_{CO2}$0.339$[\$/(g_{CO2}/s)h)]$

The justifications for cost assumptions are given in [5]. The minimum ventilation rate is set to 2.65 [m^3/s], which corresponds to approximately 1 air exchange per hour, and the maximum ventilation rate is set to 265 [m^3/s], which is approximately 100 air exchanges per hour. All other equipment can operate at any quantity greater than 0.

## 2d. Conservation of mass and energy¶

To satisfy the laws of conservation of mass and energy, the solution is contrainted by an energy balance on the (sensible) heat in the grow facility air, a moisture balance on the water content of the grow facility air, and a CO2 balance on the concentration of CO2 in the grow facility air. These are equality constraints the optimization problem must satisfy.

### Energy balance:¶

The energy balance is shown in figure 2 and equation 1. All terms have units [kW].

Figure 2. Energy balance of growing space.

$$(Eq. ~ 1) ~ ~ ~ ~ ~ T_{outside}\cdot \dot{V} \cdot \rho_{air} \cdot C_{p,air} + \dot{q}_{lights} + \dot{q}_{heat} - \dot{q}_{mechcool} - 2.26 \cdot (\dot{q}_{mist} + \dot{q}_{transpiration}) + 19.9 \cdot \dot{q}_{CO2} + U \cdot SA \cdot (T_{outside} - T)/1000 = T \cdot \dot{V} \cdot \rho_{air} \cdot C_{p,air}$$

where,

$T_{outside}$ is the outdoor air temperature $[^{\circ}K]$, $T$ is the air temperature inside the grow facility $[^{\circ}K]$, $\dot{V}$ is the ventilation rate $[m^{3}/s]$, $\rho_{air}$ is the density fo air $[kg/m^{3}]$, $C_{p,air}$ is the heat capacity of air $[kJ/(kg\cdot ^{\circ}K)]$, $\dot{q}_{lights}$ is heat generation from the grow lights $[kW]$, $\dot{q}_{heat}$ is heat generation from gas heating $[kW]$, $\dot{q}_{mechcool}$ is heat removed using mechanical cooling $[kW]$, $-2.257 \cdot \dot{q}_{mist}$ is the sensible heat removed from the air due to misting $[kW]$, $-2.257 \cdot \dot{q}_{transpiration}$ is the sensible heat removed from the air due to transpiration $[kW]$, $19.87 \cdot \dot{q}_{CO2}$ is the heat generated when CO2 is "injected" through burning natural gas $[kW]$, $U$ is the conductance of the farm roof and walls $[W/(m^{2}\cdot ^{\circ}K)]$, and SA is the surface area of the roof and walls $[m^{2}]$.

$\dot{q}_{mist}$ and $\dot{q}_{transpiration}$ have the units $[g_{water}/s]$. Multiplying by factor of 2.257 $[kJ/g_{water}]$, which is the latent heat of vaporization of water, converts their units to $[kW]$. $\dot{q}_{CO2}$ has the units $[g_{CO2}/s]$ and the heat this process gives off $[kW]$ is accounted for by the factor 19.87, which has units $[kJ/g_{CO2}]$.

### Moisture balance:¶

The moisture balance is shown in figure 3 and equation 2. All terms have the units $[g_{water}/s]$.

Figure 3. Moisture balance of growing space.

$$(Eq. ~ 2) ~ ~ ~ ~ ~ C_{max} \cdot (RH_{insidetemp}/100) \cdot \dot{V} + \dot{q}_{transpiration} + \dot {q}_{mist} + \dot {q}_{dehumid} = C_{max} \cdot (RH/100) \cdot \dot{V}$$

where,

$C_{max}$ is the water holding capacity of air at the indoor temperature [$g_{water}/m^{3}$], $RH$ is the relative humidity inside the farm $[\%]$, $\dot{V}$ is the ventilation rate $[m^{3}/s]$, $RH_{insidetemp}$ is the relative humdity of the outside air at the indoor air temperature $[\%]$ (can be greater than 100%), $\dot{q}_{transpiration}$ is the rate of evapotranspiration from the plants $[g_{water}/s]$, and $\dot{q}_{mist}$ is the rate of misting into the indoor farm air $[g_{water}/s]$.

$RH_{insidetemp}$ is calculated by using equations from [6], which approximate the relationship between air temperature, relative humidity, and dewpoint temperature. First, the dewpoint temperature of the outside air is calculated.

$$(Eq. ~ 3) ~ ~ ~ ~ ~ T_{dewpoint} = (T - 273.15) - (100 - RH_{outside})/5 + 273.15$$

Then, the dewpoint temperature is used, along with the inside temperature, $T$, to calculate $RH_{insidetemp}$.

$$(Eq. ~ 4) ~ ~ ~ ~ ~ RH_{insidetemp} = 100 - 5 \cdot ((T - 273.15)- (T_{dewpoint} - 273.15))$$

The water holding capacity of air, $C_{max} [g_{water}/m^{3}]$, can be linearly approximated within the indoor temperature range explored in this report.

$$(Eq. ~ 5) ~ ~ ~ ~ ~ C_{max} = (0.7278 \cdot (T \cdot 9/5 - 459.67) - 32.189)$$

### CO2 balance:¶

The CO2 balance is shown in figure 4 and equation 5.

Figure 4. CO2 balance of growing space.

$$(Eq. ~ 5) ~ ~ ~ ~ ~ CO2 \cdot \dot{V} = CO2_{outside} \cdot \dot{V} - \dot{q}_{photosynthesis} + \dot {q}_{CO2}$$

where,

$CO2$ is the concentration of CO2 in the inside air [$g_{CO2}/m^{3}$], $CO2_{outside}$ is the concentration of CO2 is the air outside [$g_{CO2}/m^{3}$], $\dot{q}_{photosynthesis}$ is the rate of CO2 use by the plants, and $\dot {q}_{CO2}$ is the rate of CO2 generation from a CO2 generator [$g_{CO2}/s$].

## 2e. Standard form of optimization problem¶

Minimize:

$\dot{q}_{heat} \cdot C_{heat} + \dot{q}_{mechcool} \cdot C_{mechcool} + \dot{q}_{CO2} \cdot C_{CO2} + \dot{V} \cdot C_{vent} + \dot{q}_{mist} \cdot C_{mist} + \dot{q}_{dehumid} \cdot C_{dehumid}$

Subject to:

Conservation of mass and energy:

$~ 1) ~ T \cdot \dot{V} \cdot \rho_{air} \cdot C_{p,air} = T_{outside}\cdot \dot{V} \cdot \rho_{air} \cdot C_{p,air} + \dot{q}_{heat} + \dot{q}_{lights} - \dot{q}_{mechcool} - 2.26 \cdot (\dot{q}_{mist} + \dot{q}_{transpiration}) + 19.9 \cdot \dot{q}_{CO2} + U \cdot SA \cdot (T_{outside} - T)/1000$

$~ 2) ~ C_{max} \cdot (RH/100) \cdot \dot{V} = C_{max} \cdot (RH_{outside}/100) \cdot \dot{V} + \dot{q}_{transpiration} + \dot {q}_{mist}$

$~ 3) ~ C_{CO2} \cdot \dot{V} = C_{CO2,outside} \cdot \dot{V} - \dot{q}_{respiration} + \dot {q}_{CO2}$

Grow space requirements:

$~ 4) ~ 50 \leq RH \leq 80$

$~ 5) ~ 294.26 \leq T \leq 299.82 \ast$

$~ 6) ~ CO2 \geq 0.7368 \ast$

Changing water holding capacity of air vs. temperature:

$~ 7) ~ T_{dewpoint} = (T_{outside}-273.15) - ((100 - RH_{outside})/5) + 273.15$

$~ 8) ~ RH_{newtemp} = 5 \cdot ((T_{dewpoint} - 273.15) - (T-273.15)) + 100$

$~ 9) ~ C_{max} = (0.7278 \cdot (T \cdot 9/5 - 459.67) - 32.189)$

Limitations of HVAC equipment operation:

$~ 10) ~ \dot{q}_{heat}, \dot{q}_{mechcool}, \dot{q}_{CO2}, \dot{q}_{mist}, \dot{q}_{dehumid} \geq 0$

$~ 11) ~ \dot{V} \geq 2.65$

The constraints 1,2, and 3 are non-linear, but can be made linear by fixing $T$ (which fixes $C_{max}$) and $\dot{V}$ and solving the optimization problem over a range of $\dot{V}$. The $\dot{V}$ with the lowest-cost solution is the approximate solution to the overall optimization problem. The solution is limited to the discrete values of $\dot{V}$ used. Constraints 5, 7, 8, and 9 become free of variables and become simply equations. This strategy of formulating the LP is implemented in the function named $temp70LP$. The solution to the NLP with all variables in the above formulation remaining as variables is implemented in the function named $tempRangeNLP$.

$\ast$ Minimum CO2 and temperature range are varied as specified in the results. The LP sets temperature temperature equal to 294.26 $[^{\circ}K]$ (70 $[^{\circ}F]$).

[6] Lawrence, M. G. (2005). The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air: A Simple Conversion and Applications. AMERICAN METEOROLOGICAL SOCIETY. http://journals.ametsoc.org/doi/pdf/10.1175/BAMS-86-2-225

# 3. Solution¶

This section of the report shows the implementation of the solution to the optimization problem using 4 functions. The first function $temp70LP$ is an LP and was developed to determine the optimal operational parameters at a fixed temperature of 70 $[^{\circ}F]$ by varying $\dot{V}$ at discrete intervals. The function $tempRangeNLP$ is an NLP and was developed to determine the optimal operational parameters leaving all variables in the constraints as variables. The function $temp70LPTMY$ uses the weather data and the function $temp70LP$ to calculate the optimal operating parameters at each hour of the year. The function $tempRangeNLPTMY$ uses the weather data and the function $tempRangeNLP$ to calculate the optimal operating parameters at each hour of a typical year in Madison. $tempRangeNLPTMY$ takes the arguments CO2_L, which is the minimum CO2 level $[g/m^{3}]$ in the indoor farm and MAX_TF $[^{\circ}F]$, which is the maximum temperature in the farm, to allow the calculation of annual operating costs with these parameters varied.

In [1]:
# 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


## 3a. LP at a single hour "temp70LP"¶

In [3]:
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;

WARNING: Method definition temp70LP(Any, Any) in module Main at In[2]:2 overwritten at In[3]:2.


## 3b. NLP at a single hour "tempRangeNLP"¶

In [4]:
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;


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


## 3c. LP for each hour in year "temp70LPTMY"¶

In [6]:
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;


## 3d. NLP for each hour in year "tempRangeNLPTMY"¶

In [7]:
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;


# 4. Results and discussion¶

## 4a. LP results at specified outdoor conditions¶

The $temp70LP$ function was run at temperatures of 80, 60 and 0 $[^{\circ}F]$ at $RH$ = 50 $[\%]$. Ventilation rate was varied logarithmically from 2.675 $[m^{3}/s]$ to 267.5 $[m^{3}/s]$ (or approximately from 1 to 100 air exchanges per hour). Logarithmic spacing of possible ventilation rates was chosen because preliminary results were more sensitive at low ventilation rates, which occur at relatively hot or cold outdoor temperatures.

The cost as a function of ventilation rate at each temperature is plotted is figure 5, and the optimal operating conditions at each temperature are presented in Table 3. At 80 $[^{\circ}F]$, the optimal operation occurs at the lowest ventilation rate allowed, 2.675 $[m^{3}/s]$. This is because at an outside temperature higher than the inside temperature, the ventilation air needs to be cooled. At 60 $[^{\circ}F]$, the minimum HVAC cost occurs at 55 $[m^{3}/s]$, and at 0 $[^{\circ}F]$, the minimum HVAC cost occurs at 6.47 $[m^{3}/s]$. These ventilation rates exhaust the heat generated from the lighting, while minimizing the other inputs. At lower ventilation rates, mechanical cooling and dehumidification are needed. At higher ventilation rates heating and misting are needed to maintain temeperature and $RH$.

In [8]:
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");


Figure 5. Hourly cost vs. ventilation rate at various outside temperatures, $RH_{outside}$ = 50 [%].

Table 3. Optimal costs and operational parameters at 80, 60, and 0 $[^{\circ}F]$ and RH = 50 $[\%]$.

80 $[^{\circ}F]$ / RH = 50 $[\%]$ 60 $[^{\circ}F]$ / RH = 50 $[\%]$ 0 $[^{\circ}F]$ / RH = 50 $[\%]$
NLP \$1.57623 \$ 2.46952 \$37.26929 ## 4c. Annual cost and operation in Madison, Wisconsin¶ Annual Typical Meteorlogical Year weather data for Madison is plotted in figure 6. In [9]: 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");  ****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit http://projects.coin-or.org/Ipopt ******************************************************************************  Figure 6. Typical Meterorlogical Year weather in Madison, Wisconsin. The costs at each hour of the typical meteorlogical year were estimated using both the LP and NLP formulations with the farm's indoor air temperature set to 70$[^{\circ}F]$and CO2 level to a minimum of 0.739$[g/m^{3}]$(400 ppm). The results are plotted in figure 7. For these indoor conditions, the sum of the annual costs is estimated at \$72,368.79 using the LP and \$72,096.62 using the NLP. In [10]: plot(hourOfYear70, optcost_TMY70, label = "LP") plot(hourOfYearR, optcost_TMYR, label = "NLP") xlabel("Hour of year") ylabel("Hourly HVAC cost") legend(loc = "right");  Figure 7. HVAC costs at each hour of the year for indoor temperature of 70$[^{\circ}F]$and CO2 level of 400$[ppm]$. The function tempRangeNLPTMY was used to vary the maximum indoor temperature and minimum CO2 levels to determine the changes in cost from allowing a higher temperature in the growing facility and from requiring a higher level of CO2. Table 5 presents the annual HVAC costs at various temperatures and CO2 concentrations. Allowing the temperature to increase to 80$[^{\circ}F]$instead of 70$[^{\circ}F]$decreases the HVAC costs while increasing the minimum level of CO2 increases the costs. In [11]: # 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);  Table 5. Annual HVAC costs at various maximum temperatures and CO2 levels. CO2/Max Temp 70$[^{\circ}F]$75$[^{\circ}F]$80$[^{\circ}F]$400 [ppm] \$ 72,096.62 \$43,132.32 \$ 25,018.07
800 [ppm] \$166,523.51 \$ 136,426.39 \$106,227.07 1200 [ppm] \$ 231,764.79 \$198,037.52 \$ 166,855.58

Table 6. Percent breakdown of annual HVAC costs at various maximum temperatures and CO2 = 400 ppm.

Cost type/Max Temp 70 $[^{\circ}F]$ 75 $[^{\circ}F]$ 80 $[^{\circ}F]$
Ventilation 7.0 $\%$ 11.6 $\%$ 17.6 $\%$
Heat 0.0 $\%$ 0.0 $\%$ 0.0 $\%$
Misting 5.9 $\%$ 11.2 $\%$ 19.9 $\%$
Mechanical cooling 54.4 $\%$ 42.2 $\%$ 22 .8 $\%$
Dehumidification 21.5 $\%$ 16.2 $\%$ 7.2 $\%$
CO2 generation 11.2 $\%$ 18.8 $\%$ 32.4 $\%$

To gain further insight into the cost behavior of the HVAC system in Madison, the hourly HVAC costs are plotted vs. outside air temperature for the cases in the first row and third column of Table 5. Figure 8 presents the hourly HVAC costs vs. outside air temperature for several maximum temperatures (T = 70 $[^{\circ}F]$, 75 $[^{\circ}F]$, and 80 $[^{\circ}F]$), while holding CO2 constant at 400 [ppm]. Figure 9 presents the results for several levels of CO2 enrichment (CO2 = 400 [ppm], 800 [ppm], and 1200 [ppm]), while holding maximum temperature constant at 80 $[^{\circ}F]$. As can be seen from figure 8, cost increases rapidly as the outdoor temperature approaches the maximum indoor temperature. This is because ventilation can no longer be used for cooling and refrigeration must be used instead. As can be seen in figure 9, a higher concentration of CO2 decreases the outdoor air temperature at which cost rises sharply. This is due to the heat generated when CO2 is burned.

In [12]:
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");


Figure 8. Hourly HVAC costs vs. outside air temperature for hours in Madison Typical Meteorlogical Year for several maximum indoor air temperatures.

In [13]:
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");


Figure 9. Hourly HVAC costs vs. outside air temperature for hours in Madison Typical Meteorlogical Year for several levels of CO2 enrichment.

## 5. Conclusion¶

This report develops a mathematical model and solution for determining the optimal operational parameters of an indoor farm. This program could be used in the following ways:

1. Estimating HVAC costs to aid in determining indoor farm profitability.
2. Determining which pieces and what sizes of HVAC equipment are needed for optimal operation.
3. Controlling an indoor farm's HVAC system in real-time.

The following insights were drawn from the simulations run in this report:

1. Mechanical cooling/dehumidification is expensive. Ventilation and misting (i.e. evaporative cooling) should be used instead of mechanical cooling when conditions allow.
2. Increasing the allowable temperature in the farm from 70 [F] to 80 [F] very significantly decreases the HVAC operating costs by reducing the usage of mechanical cooling/dehumidification.
3. In Madison, Wisconsin, the farm considered gets plenty of heat from lights and CO2 injection and would not need a dedicated "heating" system.

One future improvement to this optimization program is the incorporation of a model for evapotranspiration, instead of assuming a constant rate of moisture generation. This may affect how much dehumidification is needed. Another improvement to the model is to optimize the lighting and CO2-injection schedules. This program assumes the farm's grow lights are always on and that CO2-enrichment is always occuring. However, grow lights are not typically on 24/7 and CO2-enrichment can be done intermittently at times when the ventilation rate is low. A schedule could be optimized based on the expected weather conditions and price of electricity to optimize the lighting and CO2-injection strategies.

## 4. Appendix¶

In [14]:
# 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)

37.27254720083333
2.675
0.0
820.251886111111
0.0
119.62821666666662
2.73
2.469549167974722
55.011704248326424
0.0
0.0
166.40396205740018
0.0
2.729999999999997
1.5825595928836813
6.473843108324972
0.0
0.0
119.66272469465879
0.0
2.73

In [15]:
# 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)

At T_outside = 60 F, RH_outside = 50 %, optimal hourly cost as found with temp70LP:
2.469549167974722
At T_outside = 60 F, RH_outside = 50 %, Optimal hourly cost as found with tempRangeNLP:
2.469522859336779

At T_outside = 80 F, RH_outside = 50 %, optimal hourly cost as found with temp70LP:
37.27254720083333
At T_outside = 80 F, RH_outside = 50 %, Optimal hourly cost as found with tempRangeNLP:
37.26929729899658

At T_outside = 0 F, RH_outside = 50 %, optimal hourly cost as found with temp70LP:
1.5825595928836813
At T_outside = 0 F, RH_outside = 50 %, Optimal hourly cost as found with tempRangeNLP:
1.5762340306769056

In [18]:
# 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))

25018.07056485681
106227.07626708364
166855.5813255608
43132.320206197095
136426.3860972927
198037.52962807103
72096.61739405835
166523.51957196055
231764.79016878596

In [19]:
# 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))

Breakdown of annual costs when max indoor temp is 80 and min CO2 is 400 ppm
6.716536238769925e-10
0.2284256903180522
0.1994431384891465
0.07215203115362598
0.32405041420653774
0.17592872516098415

Breakdown of annual costs when max indoor temp is 70 and min CO2 is 400 ppm
2.1598848962720688e-10
0.543630180858043
0.059077280916909836
0.21471310880698882
0.11244793804395783
0.07013149115811217

Breakdown of annual costs when max indoor temp is 75 and min CO2 is 400 ppm
3.75242851953314e-10
0.42168102020178405
0.11215791867665517
0.1623501201497735
0.1879591902881049
0.1158517503084396

In [20]:
# 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))

The estimated annual cost of using the LP is:
72368.79197067142
The estimated annual cost of using the NLP is:
72096.61739405835