# VSG Group 3: CEA HVAC System

Team members

* Adam Waterman (problem presenter) [e-mail]
* Ann Smith (facilitator) [e-mail]
* Andrew Lacey [e-mail]
* Nam Nguyen [e-mail]
* Chris Breward [e-mail]

# Chat
I've dumped a lot of stuff in Overleaf verbatim, feel free to edit .
Nam Great comment! please do let me have outputs etc for the report. Directly into Overleaf or here. Thanks.(AS)

**Nam**: Hi Ann, I forgot what comment was ... There is not a lot of outputs at the moment because the modelling/control optimisation are still in progress. Nevertheless, here is a GitHub link to the code: https://github.com/namhhnguyen/Maths4Food/blob/main/PlantAirControl.ipynb
Cheers Nam.

I've just dumped all the Hack MD work in Overleaf as a basis for the final report (AS)

# Food Security – Virtual Study Group
Tuesday 16th- Thursday 18th January 2024

# Introduction
This challenge focusses on a vertical farm or CEA hvac system. The aim is to maintain target
temperature, relative humidity and vapor pressure deficit within the chamber while using the
minimum amount of energy. All of the target factors affect each other so changing one changes the others. Also, methods of control affect multiple target factors in different ways. This leads to a complex multivariate system with multiple target parameters and multiple interlinked control
methods

# The hvac system
The hvac (heating ventilation and air condition) system is a CEA (controlled Environment Agriculture) controlled growing room.

The system is heated and cooled by opening and closing vents and heat exchanger grates. To reduce the temperature the lower cooling grate is opened to increase air flow from outside the system. Humidity can be controlled either by increasing air flow in the chamber, pumping in fresh air to displace the higher humidity air. Also, humidity can be reduced by reducing the temperature. However, if the temperature is reduced below the dew point condensation will be introduced to the system which will itself affect the humidity. An excess of water can be relesed through a drain in the base of the unit. The dew point is the atmospheric temperature (varying according to pressure and humidity) below which water droplets begin to condense and dew can form.

Control is currently governed by simple proportional integral differential (PID) loops. A well established method to control process variables via a feedback and updating mechanism. Once steady state is achieved, targets are generally maintained with good acuracy. Current control protocols hold conditions accurate to 0.2C. However, the practice is not efficient as there can be large shifts in conditions and responding dynamically to temperature and humidity corrections can lead to high energy costs.

In short, the hvac system is a complex system of few but highly correlated variables, under the current control sytem high levels of power effort are required to maintain equilibrium. Thus control via PID loops, although effective, is not an environmentally sound solution.

# 

Figure: example of a simple vertical farm or CEA hvac system. It should be assumed that there are a significant amount of plants, lighting and irrigation equipment within the chamber. This causes the temperature to rise and increases the humidity in the air.
## Assumptions
Assume the chamber is very well insulated. So heat is not gained or lost through walls.
The volume and pressure of the room is constant.

Equipment such as lighting in the room uses energy which ends up as heat.

The crops and the irrigation system generate humidity.

## Aim
To maintain target temperature, $T$, relative humidity, $R$, and vapour pressure deficit, $VPD$, within the chamber while using the minimum amount of energy.
All of the target factors affect each other so changing one changes the others. Also all of the different methods of control usually affect multiple target factors in different ways.
This leads to a complex multivariate system with multiple target parameters and multiple interlinked control methods.

## Factors at play
Temperature of the chamber and fresh / outside in air (degree c)
Relative humidity (RH) (expressed as a percent) also measures water vapor, but RELATIVE to the temperature of the air. In other words, it is a measure of the actual amount of water vapor in the air compared to the total amount of vapor that can exist in the air at its current temperature.
Absolute humidity (expressed as grams of water vapor per cubic meter volume of air) is a measure of the actual amount of water vapor (moisture) in the air, regardless of the air's temperature. The higher the amount of water vapor, the higher the absolute humidity.
Dew point is the temperature the air needs to be cooled to (at constant pressure) in order to achieve a relative humidity (RH) of 100%. At which point condensation occurs.
Vapor pressure deficit (VPD) is calculated as the difference between the amount of moisture that's actually in the air and the amount of moisture that air could hold at saturation.
Energy. (joules) the energy used to increase or decrease the air temperature using the heating and cooling grates

## Controls
**Temperature** can be control by:
Increasing the cooling grate will cause the air temperature to drop.
Increasing the heating grate will cause the air temperature to rise
Adding fresh air will either increase or decrease the air temperature depending on the relative air temperature outside.
Relative humidity can be control by:
Raising the temperature will cause the RH to drop.
Lowing the temperature will cause RH to increase.
Reducing the absolute humidity will cause the RH to drop.
Increase the absolute humidity will cause the RH to increase.
Absolute humidity can be control by:
Lowing the temperature of the cooling grate to below the dew point of the air leading to condensation. (the heating grate can then be used to raise the air temperature back up to its original temperature)
Adding fresh air fresh air will either increase or decrease the absolute humidity depending on if the fresh air has a high or lower absolute humidity. The fresh air temperature can then be brought up or down as required.
A humidifier can be added to increase humidity but in reality is almost never needed because all of the plants and irrigation systems create a lot of humidity. So only a decrease in humidity is required.

## Solution

1. Identify the underlying dynamical system as a mathematical ODE model.
2. Parameterise the model (data? reinforcement learning? wishful thinking??)
3. Code and run scenarios.

# Mathematical Model
**Growing chamber** conditions:
**(Absolute) Humidity:** ($g/m^3$)
Amount of water vapour in the air. Grams per meter cubed.
**Temperature** (°C)
**Relative Humidity** (RH)(%):
measure of proportion of actual amount of water vapour in air: total that can exist at a given temperature.
**Vapour Pressure Deficit (VPD)** difference? between the moisture in the air and the saturation amount.

$VPD=VP_{sat} -VP_{air}$

Example chart for growing Basil, ideal region highlighted by the white box.


------
Chris and Andrew notes

Let $T$ be the temperature in the air and $t$ be time. We'll write down a model for what's going on without any controls other than the lights going on and off:

Cons of temp in air:
$$\rho c_p V \frac{\mbox{d}T}{\mbox{d}t}=\lambda(t)=\lambda_0(1+\lambda_1f(t)),$$
where $V$ is the volume of the room (not measured in cubic feet!), $\lambda_0$ is the rate of energy production when the lights are off and $\lambda_1$ is the rate of energy production when the lights are on, and
$f(t)=\Theta(t-t_1)\Theta(t_2-t)$ is the lights schedule, where $\Theta(\cdot)$ is the Heaviside function, and $t_1$ and $t_2$ are the times that the lights are turned on/off, respectively (we only consider $24$ hrs; could repeat as needed). (Comment from AAL: Looking at humidity and latent heat, and comparing with specific heat of air x temperature change of about 20$^\circ$C, it appears heat contents are a vey similar size; we shall need to allow for latent heat of evaporation. We should also perhaps be allowing for heat content for soil etc. If this soil (etc.) has a low heat capacity, we may neglect it; if the heat transfer is ``fast'', we can simply include the soil's heat capacity on the left-hand side. Otherwise, an extra heat equation will be needed.)

We denote by $H$ the absolute water content, measured in mol m$^{-3}$ (Comment from AAL: Might be better to use kg m$^{-3}$?), and by $R$ the relative humidity given by $$R=\frac{H}{S(T)},$$ where $S(T)=C e^{\nu T}$ is the saturation level ($C$ and $\nu$ are constants), which is the most water vapour that the air can hold at a given temperature. Conservation of water in the air gives:
$$V \frac{\mbox{d}H}{\mbox{d}t}=\kappa_s(T,R)+\kappa_P(T,R,f)$$
where $$\kappa_s(T,R)=A_0(1+A_1T)(1-R)$$ is the supply rate of water due to evaporation from all non-plant surfaces (eg soil), which goes up as the temperature goes up and down as the relative humidity goes up (we picked linear for convenience); $A_0$ and $A_1$ are constants, and $$\kappa_p(T,R,f)=A_2(1+A_3T)(1-R)(1+A_4 f)$$

NB could eliminate one of $R$ and $H$ from the equations.

Now include the introduction of fresh air, with external temperature $T_e$ and absolute humidity $H_e$, with volume flow rate $Q_v$ (mass flow rate $\rho Q_v$, taking the air to have approximately constant density), assumed much less than the recirculation rate, say $Q_r$ (which ensures that the air doesn't suddenly become supersaturated), and equal venting rate. The net rate of change of energy per unit time is then (input less output) $c_p \rho Q_v (T_e - T)$ and corresponding rate of changes of vapour content is $Q_v ( H_e - H )$.

Combining with the earlier equations for a closed chamber (still ignoring latent heat of evaporation and heat capacity of soil):
$$
\rho c_p V \frac{\mbox{d} T}{\mbox{d}t}=\lambda(t) +c_p \rho Q_v (T_e - T),
$$
$$
V \frac{\mbox{d}H}{\mbox{d}t}=\kappa_s(T,R)+\kappa_P(T,R,f) + Q_v ( H_e - H ).
$$
The mixed chamber and fresh air,
with temperature and humidity approximately those of the chamber, flows from the inlet to the grates, which we must now model.

To get an idea as to how to model the two heat-exchanging grates, we first consider the hot grate as a counter-flow heat exchanger (although it is not one!). Taking the air to flow with mass flux $q_1$, to have specific heat $c_1$, and local temperature $T_1$, going from $x=0$ to $x=l$, the water to have corresponding quantities $q_2$, $c_2$, and $T_2$, it going from $x=l$ to $x=0$, and a heat-exchange coefficient $h$ (J K$^{-1}$ m$^{-1}$), the temperatures satisfy
$$
c_1 q_1 \frac{\partial T_1}{\partial x} = h (T_2 - T_1) ,
$$
$$
-c_2 q_2 \frac{\partial T_2}{\partial x} = h (T_2 - T_1) .
$$
For convenience we partially non-dimensionalise, writing $x = ly$, $\alpha_j = (hl)/(q_j c_j)$ for $j = 1, 2$:
$$
{T_1}' = \alpha_1 (T_2 - T_1) ; {T_2}' = \alpha_2 (T_2 - T_1) ;
$$
with $'=$ d/d$y$, and given inlet temperatures $T_1(0) = T_{1i}$ and $T_2(1) = T_{2i}$. Solving these gives (apart from the special case of $\alpha_1 = \alpha_2$) an outlet air temperature of
$$
T_{1O}= T_1(1) = \frac{E (1 - (\alpha_2/\alpha_1))T_{1i} + (1 - E)T_{2i}}{1 - (\alpha_2/\alpha_1)E}
$$
with $E= \exp(\alpha_2 - \alpha_1)$. (Comment from AAL: Anyone care to check??!!)will do, just dumping in Overleaf-feel free to cut as you want or just scroll through (AS)

From these we can see that (as should be expected):
for $\alpha_1$ small (air whistling through), $T_{1O} \sim T_{1i}$;
for $\alpha_1$ large (air hardly moving), $T_{1O} \sim T_{2i}$;
for $\alpha_2$ small (lots of heating water, so its temperature remains at $T_{2i}$), $T_{1O} \sim T_{2i} (1 - \exp(-\alpha_1)) + T_{1i}\exp(-\alpha_1)$;
for $\alpha_2$ large (hardly any heating water), $T_{1O} \sim T_{1i}$.

All these suggest that for the heating grate, the warmed air temperature, $T_O$, should be related to the incoming air temperature, $T_i$, by a law of the form
$$
T_O = h_h T_h + (1 - h_h) T_i ,
$$
with $T_h$ the temperature of the heating water (as it comes in), and $h_h$ a coefficient depending upon the flow rates: $h_h \to 0$ for water flow rate going to zero, $h_h \to 1$ for water flow rate going to infinity. For the present, we might take $h_h$ to vary with the hot-water flow rate as
$$
h_h = 1 - \exp(-Q_h/Q_h^*)
$$
with $Q_h^*$ some constant value.
N.B. If the hot water is recirculated, given that we want to minimise power use, we should also give the out-put hot-water temperature here (got by conservation of energy).

The cold grate might be expected to work in the same sort of way, except for occrrence of condensation. Looking at a counter-current representation, the same equations as for the hot grate might apply where the air is not saturated with water vapour, and latent heat has to be allowed for where the air is saturated:
$$
c_1 q_1 \frac{\partial T_1}{\partial x} = h (T_2 - T_1) \mbox{ and } H=const. \mbox{ where } H<R(T_1);
$$
$$
q_1 \left( c_1 + L S'(T_1) \right) \frac{\partial T_1}{\partial x}= h (T_2 - T_1) \mbox{ where } H=R(T_1);
$$
$$
-c_2 q_2 \frac{\partial T_2}{\partial x} = h (T_2 - T_1);
$$
$H$ and $T_1$ are given at $x=0$, and $T_2$ at $x=l$.
In general this free-boundary problem doesn't have an explicit solution, so, for the present, we proceed as at the final stage of the hot-grate modelling and assume a form of a law for its operation. We then take
the cooled air temperature, $T_O$, to be related to the incoming air temperature, $T_i$, by a law of the form
$$
T_O = h_c T_h + (1 - h_c) T_i + (L/(\rho c_p))(H_i - H_O),
$$
$$
H_O = S(T_O),
$$
so
$$
T_O + (L/(\rho c_p)) S(T_O)= h_c T_h + (1 - h_c) T_i + (L/(\rho c_p))
H_i,
$$
with $T_c$ the temperature of the cooling water (as it comes in), and $h_c$ a coefficient depending upon the flow rates: $h_c \to 0$ for water flow rate going to zero, $h_c \to 1$ for water flow rate going to infinity. (We could take $h_c=h_h$.) As long as $h_c$ has the expected these expected properties, along with it being increasing with flow rate, we again will have $T_O$ and $H_O$ increasing with $H_i$ and $T_i$ and decreasing with with water flow rate -- we're assuming that $T_c < T_i$ and that condensation will occur at the grate. For the present, we might again take $h_c$ to vary with the cold-water flow rate as
$$
h_h = 1 - \exp(-Q_h/Q_h^*)
$$
with $Q_h^*$ some constant value, despite any effects of water condensation.
N.B. If the cold water is recirculated, given that we want to minimise power use, we should also give the out-put hot-water temperature here (got by conservation of energy).

Putting the effects of the grates into the model:

The cold grate sees input temperature and humidity (approximately) $T$ and $H$, so its outputs will be $T_i$ and $H_i = S(T_i)$ given by

$$
T_i + (L/(\rho c_p)) S(T_i)= h_c T_h + (1 - h_c) T + (L/(\rho c_p))
H.
$$
These then come into the hot grate which outputs temperature and humidity $T_O$ and $H_O = H_i$, with
$$
T_O = h_h T_h + (1 - h_h) T_i .
$$
Air with this temperature and humidity (re)enters the chamber, with mass flow rate $Q_r \rho$, while the same flow rate takes air of temperature $T$ and $H$ out. Accounting for the further flow differences we then get
$$
\rho c_p V \frac{\mbox{d} T}{\mbox{d}t}=\lambda(t) +c_p \rho Q_v (T_e - T) +c_p \rho Q_r (T_O - T),
$$
$$
V \frac{\mbox{d}H}{\mbox{d}t}=\kappa_s(T,R)+\kappa_P(T,R,f) + Q_v ( H_e - H ) + Q_r ( H_O - H ).
$$
We still have to put numbers in -- and look at power!!

------
### CODING (Nam)
Coding up a struture for a proof of concept Hvac model, e.g., set up observations, control actions, and a reward function that we wish to optimise for. Next stage will be to work with Andrew, Ann, and Chris to bring actual maths in.

# Signal Outputs

Output signals captured during operation between 13:28 and 22:59 hours, March 2023.
Signals:
Date and time stamp (minutes), H_Target (target humidity value), H_Hvac Avg (average humidity value per minute), H_Effort (effort measured in Joules to maintain target humidity value), T_Effort (effort in Joules to maintain target termperature value), T_Hvac Avg (average temperature value per minute, degress Celcius) , T_target (tartget temperature value).

Plants like other living organisms require periods of rest, figures abd...

Effort to maintain target temperatures and humidity values is most pronounced during hours of daylight. Overnight, especially pre dawn (fig ?? 03.40 to 05.00 hours) there is

Figures:

Figure

Figure

Figure




### DEWPOINT AND VAPOUR PRESSURE DEFICIT EQUATIONS
From Tetens Formula, the relation between temperature and the partial pressure of water vapour:

$e(millibars) = 6.1078 exp( (17.269*T) / (237.3+T) )$
where,
$e$ is saturated vapour pressure in millibars
$T$ is temperature in degrees C and the equation for relative humidity:

$Rh=(ea/es)*100$
where,
$ea$ is the actual vapour pressure or vapor pressure at dewpoint temperature
$es$ is the saturation vapour pressure or vapor pressure at air temperature it can be shown that:

**Dew point temperature = **
$Tdew = (237.3*X) / (17.269‐X)$

where $X=ln(Rh/100)+((17.269*Tair)/(237.3+Tair))$
$Tair$ is air temperature
$Tdew$ is dew point temperature

Vapor Pressure Deficit = es‐ea = es‐(Rh*es/100) at any instant.

********************************

# DEFINITIONS:

SATURATION VAPOUR PRESSURE: water evaporating and condensing through the air‐water interface yields zero net transport.

Ideal gas law to sat vap pr:
e = pRT
e=vapor pressure in millibars
p=vapor density in mass/unit volume (g/cubic meter)
T=absolute temperature (degrees Kelvin)
R=vapoUr gas constant

If we use the universal gas constant for R, the formula becomes
$e = 1.61pRT$

which becomes
$p = 0.622e/RT$

This is known as VAPOUR DENSITY or ABSOLUTE HUMIDITY of the atmosphere.

RELATIVE HUMIDITY is the ratio of the vapor density (or pressure) to the saturation vapor density (or pressure) at the same temperature.

‐‐‐‐‐‐‐‐‐‐‐‐‐‐

# APPROXIMATIONS:

For saturated vapor pressure in temperature range 25‐55 deg Fahrenheit:
$e == 6.11 + 0.339(D‐32)$

$e$ = sat vap pr in millibars
$D$ = dew point temperature in Fahrenheit

For saturated vapor pressure in temperature range 25‐55 deg Fahrenheit:
$e == 0.18 + 0.01(D‐32)$

$e$ = sat vap pr in inches of mercury
$D$ = dew point temperature in Fahrenheit

Relative humidity can be approximated from air and dew point by

$rh = 100 ((112 ‐ 0.1T + D) / 112 + 0.9T)**8$

$T$ = temperature
$D$ = dew point
(Range of valid temperatures is unknown).

Dew point can be approximated in temp range ‐40 to 50 deg C by

$T ‐ D = (14.55+0.114T)*x + ((2.5+0.007T)*x)**3 + ((15.9+0.117T)*x)**14$

$T$ = temperature
$D$ = dewpt
x = complement of relative humidity in decimal form, ie. x = 1.0 ‐ (rh/100)
-----