Optimal Feedback Control for HVAC Systems: An Integral Sliding Mode Control Approach Based on Barrier Function

Designing of each HVAC system’s control is to provide a suitable and desirable environment for human life and achieve comfortable and justifiable indoor air quality (IAQ). Main part of HVAC is Air-handling units AHU. AHU which provides a desired supply air with specific temperature and humidity. This units constituting of various equipment and mechanical parts. Therefore, the possibility of complete identification and derivation of equations and designing an effective control system by classical methods is impossible and unattainable and it is difficult to obtain an accurate mathematical dynamic model of AHU [1].

The dynamic mathematical model for AHUs is essential to present an effective control method which can operate in an acceptable way taking into account all practical constraints, it is possible to have an overview of the different modelling techniques to best improve the control strategy of the HVAC systems developed in studies [2-4]. Temperature and humidity are the most important controlled parameters (state variables) in air handling units. Humidity is a significant factor affecting both thermal comfort and indoor air quality. Energy efficiency is a major challenge in buildings for maintaining comfortable conditions. Reducing energy losses and improving the dynamic behavior of the system like comfort conditions for temperature and humidity are the primary objectives of HVAC control systems. In order to achieve these objectives, numerous control approaches have been used, such as the decoupling method, which aims to overcome the problem of humidity coupling by eliminating the temperature disturbance completely and not focusing on reducing the recovery time of humidity [5]. A neural fuzzy structure of decoupled parameter self-tuning fuzzy neural PID controller was proposed in the study [6]. Adaptive and learning-based approaches have also been introduced in the study [7]. Seong et al. [8] the authors performed studies of optimal HVAC control methods using Genetic Algorithms that were implemented to both a variable air volume (VAV) air-conditioning system and chilled water system for optimization of the control variables in each system. Robust control system is crucial due to the dynamic uncertainties and external disturbance. Sliding mode control (SMC) is known to be an appropriate strategy for its robustness against external disturbances and little sensitivity to parameter variations under matching conditions. The SMC needs a suitable control law that is a trajectory moves toward the sliding surface and reaches in a finite time and stays on. SMC is introduced [9, 10] which compared the performance of the system with proposed controller versus performance with PID controller to demonstrate the benefit of using SMC to enhance the desired tracking. Two strategies are implemented on nonlinear model nonlinear decoupled sliding mode control and Linear optimal robust H¥ technique, the result showed that the nonlinear decoupled sliding mode control has given better performance [11].

The fuzzy algorithm is used to investigate the stability of a sliding control system by adjusting the parameters in the approach rate, reducing the switching frequency, and weakening chatter [12]. In the study [13], adaptive sliding mode control to overcome the overvalued of controller gain in sliding mode control is proposed. In spite of the specifications of sliding mode controller it has many drawbacks like the chattering effectiveness, the reaching phase, and sensitiveness to matchless uncertainties [14, 15]. Integral sliding mode (ISMC) is one of the SMC strategies that have been proposed to solve these problems, which is looking to remove the reaching phase by enforcing sliding mode during the entire system response [16]. Integral sliding mode based on barrier function is proposed [17] was applied for a DC servo actuator system containing friction, which is enabled the controller to compensate the external disturbance and uncertainty from the first instant. The ISMC is also proposed [18] for unmatched disturbance. Finally, for the purpose of designing an optimal control the LQR was proposed to provide optimally controlled feedback gains to enable the high performance [19].

The structure of this article is as follows: Section 2 introduces the problem statement, while Section 3 describes the system and the dynamic model of the AHU system. Section 4 illustrates the sliding mode controller, and Section 5 presents the proposed control design. The simulation results are demonstrated in Section 6, and Section 7 concludes the paper.

Air handling unit (AHU) is a device used for the control of the air flow quality in HVAC systems. Figure 1 which presents a schematic view of an air-handling unit having one zone (indoor) in HVAC. HVAC system for single thermal load consist of a heat exchanger, a chiller, which provides chilled water to the heat exchanger; a circulating air fan; the thermal space; connecting ductwork; dampers; and mixing air components [9]. A mixed air of 25% of fresh air with 75% of recirculating air passes through the heat exchanger. The temperature and humidity ratio of the hot and humid air are reduced as it passes through the heat exchanger.  Finally, the desired supply air is supplied and delivered to the thermal space through the connecting ductwork.

Many assumptions must be considered [9]: 1) Ideal gas behavior; 2) Mixing air should be perfect; 3) Fixed pressure process; 4) Wall and thermal storage must be neglected; 5) thermal losses between components should be negligible; 6) Infiltration and exfiltration effects is neglected; and 7) negligible transient effects in the flow splitter and mixer.

1.png

Figure 1. AHU diagram in HVAC systems

The dynamic modeling of the AHU is given in studies [9, 11, 21], by the following differential equations.

$\dot{W}_{i a}=\frac{\dot{M}_l}{\rho_{a d} V_{i a}}-\frac{\dot{f}_{a r}}{V_{i a}}\left(W_{i a}-W_{s a}\right)$

$\dot{T}_{\mathrm{ia}}=\frac{1}{\rho_{a d} C_a V_{i a}}\left(\dot{H}_l-h_v \dot{M}_l\right)+\frac{h_v \dot{f}_{a r}}{C_a V_{i a}}\left(W_{i a}-W_{\mathrm{sa}}\right)-\frac{\dot{f}_{a r}}{V_{i a}}\left(T_{\mathrm{ia}}-T_{\mathrm{sa}}\right)$

$\dot{T}_{s a}=\frac{\dot{f}_{a r}}{V_{c u}}\left(T_{i a}-T_{s a}\right)+\frac{0.25 \dot{f}_{a r}}{V_{c u}}\left(T_{o a}-T_{i a}\right)- \frac{\dot{f}_{a r} h_s}{C_a V_{c u}}\left(0.25 W_{o a}+0.75 W_{i a}-W_{s a}\right)-\frac{\rho_{w d} C_w \delta T_{c u}}{\rho_{a d} C_a V_{c u}} \dot{f}_{w r}$              (1)

The parameters, variables, state and control are described in Tables 1 and 2 [22].

Table 1. Thermo-fluid parameters of AHU

Parameters
Via	Thermmeal space volume	Woa	Humidity ratio of outdoor fresh air
Vcu	Cooling unit volume	Wsa	Humidity proportion of Supply air
ρwd	Water mass density	Wia	Thermal zone humidity ratio
ρad	Air mass density	Toa	Outdoor fresh air temperature
hv	Enthalpy of vapors	Tsa	Conditioned air supply temperature
hs	Enthalpy of saturated water	Tia	Thermal zone air temperature
Cw	Specific heat of water	δTcu	Cooling unit temperature gradient
Ca	Specific heat of air	$\dot{M}_l$	Humidity source strength
$\dot{f}_{a r}$	Air flow rate	$\dot{H}_l$	Heat load
$\dot{f}_{w r}$	Water flow rate

Table 2. Thermo-fluid parameters values in AHU around an operating point

In state space form the nonlinear AHU model is given by;

$\begin{aligned} & \dot{x}_1=f_1+g_1 u_1 \\ & \dot{x}_2=f_2 \\ & \dot{x}_3=f_3-g_2 u_2\end{aligned}$             (2)

where,

$f_1=a_1, g_1=-a_2 x_1+a_3$,

$f_2=\left(b_2 x_1-a_2 x_2+a_2 x_3-b_3\right) u_1+b_1$

$f_3=\left(-0.75 c_2 x_1+0.75 c_1 x_2-c_1 x_3+c_3\right) u_1, \quad g_2=c_4$

$u_1=\dot{f}_{a r}, u_2=\dot{f}_{w r}, y_1=W_{i a}, y_2=T_{i a}$

$x_1=W_{i a}, x_2=T_{i a}, x_3=T_{s a}$

$a_1=\frac{1}{\rho_{a d} V i a} \dot{M}_l, a_2=\frac{1}{V_{i a}}, a_3=a_2 W_{s a}$

$b_1=\frac{1}{\rho_{a d} C_a V_{i a}}\left(\dot{H}_l-h_v \dot{M}_l\right), b_2=\frac{h_v}{C_a V_{i a}}$,

$b_3=b_2 * W_{s a}, c_1=\frac{1}{V_{c u}}, c_2=\frac{h_s}{C_a V_{c u}}$,

$c_3=0.25 c_1 T_{o a}-0.25 c_2 W_{o a}+c_2 W_{s a}, c_4=\frac{\rho_{w d} C_w \delta T}{\rho_{a d} C_a V_{c u}}$

Remark (1): According to the physical characteristics for the HVAC system which is represented in Eq. (2), the control inputs are positive quantities. So, the system is non-affine.

Firstly, the dynamic model is splitting into two parts, upper part and lower part. Eq. (7) is the upper part model which describe the humidity ratio dynamic where u₁ is the control input;

$\dot{x}_1=f_1+g_1 u_1$            (7)

while the lower part (indoor temperature) model is given by;

$\begin{gathered}\dot{x}_2=f_2 \\ \dot{x}_3=f_3-g_2 u_2\end{gathered}$              (8)

where, u₂ is the control input.

Before designing the controllers, the following remarks must be considered

Remark 2: According to the Eq. (7) and Eq. (8), u₁&u₂ must be positive (u₁>0 and u₂>0) to achieve the characteristic of the cooling mode (decreasing humidity ratio and indoor temperature).

Remark 3: From Eq. (2), the system will be open loop when u₁=0, which means:

$\begin{gathered}x_1 \rightarrow \infty \text { as } t \rightarrow \infty \text { where } x_1=W_{i a} \\ x_2 \rightarrow \infty \text { as } t \rightarrow \infty, \quad \text { where } x_2=T_{i a}\end{gathered}$

where, a₁, b₁>0.

5.1 Designing the humidity ratio control u₁

To design u₁, the error function for x₁ is defined as:

$e_1=x_1-x_{1 d}$              (9)

where, x_1d is the desired humidity ratio

In order to design the classical ISMC, the input output model with respect to e₁ from Eq. (7) is given by:

$\dot{e}_1=\dot{x}_1$             (10)

Eq. (10) can be written in terms of the nominal and perturbation form as follows:

$\dot{e}_1=\left[-a_{2o} x_1+a_{3o}\right] u_1+\delta_1$         (11)

where,

$\delta_1=\left[-\Delta a_2 x_1+\Delta a_3\right] u_1+a_1$

which denoted the perturbation of the upper part model.

Based on Eq. (11), the control law for u₁ can be taken as

$u_1=\frac{-1}{-a_{2o} x_1+a_{3o}}\left(u_{1 n}+u_{1 s}\right)$            (12)

where, (a_2o, a_3o) refer to the nominal parameters, u_1n is the linear control, while u_1s is the discontinuous control which designed to reject the perturbation term.

According to the ISMC design the sliding variable of ISMC is defined as [23, 27]:

$s_1=s_{1 o}(e)+z_1$              (13)

where, s₁ includes two parts: The first one s_1o(e) which can be designed as a linear combination of the system states, like in the conventional sliding mode. The second part z₁ is the integral term for the first controller and will be derived below.

So

$s_{1 o}=e_1$               (14)

and, the derivative of the integral sliding variable is given by:

$\dot{s}_1=\dot{e}_1+\dot{z}_1$           (15)

with z₁(0)=-e₁(0), which ensure that s₁(0)=0, $\forall t \geq 0$.

By substituting Eq. (11) and Eq. (12) into Eq. (15), we obtain

$\dot{s}_1=-\left(u_{1 n}+u_{1 s}\right)+\delta_1+\dot{z}_1$               (16)

Define the dynamic of the integral part as

$\dot{z}_1=u_{1 n}$            (17)

Accordingly, we obtain

$\dot{s}_1=\left(-u_{1 s}\right)+\delta_1$            (18)

By applying the equivalent control [33] to Eq. (18) yield:

$\begin{gathered}\dot{s}_1=0=\left[-u_{1 s}\right]_{e q}+\delta_1 \\ \text { So }\left[u_{1 s}\right]_{e q}=\delta_1\end{gathered}$

From Eq. (11) the resultant error dynamic in the equivalent mode is given by:

$\dot{e}_1=-u_{1 n}$              (19)

Therefore, u_1n can be selected as in Eq. (20) which makes the origin of the error dynamics in Eq. (19) globally asymptotically stable.

$u_{1 n}=\alpha_1 e_1$           (20)

where, α₁ is a positive constant, therefore,

$\dot{z}_1=\alpha_1 e_1$              (21)

The discontinuous control term in Eq. (18) is simply given by:

$u_{1 s}=k_1 \operatorname{sign}\left(s_1\right)$            (22)

Selecting the discontinuous control gain k₁ such that it satisfies the inequity.

$k_1>\left|\delta_1\right|$             (23)

The first control, eventually, is given by:

$u_1=\frac{-1}{-a_{2 o} x_1+a_{3 o}}\left(\alpha_1 e_1+k_1 \operatorname{sign}\left(s_1\right)\right)$             (24)

5.2 Designing the indoor temperature control u₂

Defined the error function of the temperature state x₂ as:

$e_2=x_2-x_{2 \mathrm{~d}}$              (25)

where, x_2dis the desired indoor temperature.

By considering e₂ as the output, the relative degree is two. So, the input-output model for the lower part of the system model is given by:

$\begin{gathered}\dot{e}_2=\dot{x}_2=f_2=e_3 \\ \dot{e}_3=F(x)-G(x) u_2\end{gathered}$                  (26)

where, $F(x)=\left(b_2 \dot{x}_1-a_2 f_2+a_2 f_3\right) u_1+\left(b_2 x_1-a_2 x_2+a_2 x_3-b_3\right) \dot{u}_1$ and $G(x)=a_2 g_2 u_1$. e₂ & e₃ will be needed later to constract the sliding variable for the lower system model.

Remark 4: Since f₂ is uncertain function, e₃ can’t be computed exactly. Hence a robust differentiator is proposed to obtain e₃ (the derivative of e₂), where the ACSMD is used here as presented in section (5.4).

The error dynamic model for the lower system model can also be written as:

$\begin{gathered}\dot{e}_2=e_3 \\ \dot{e}_3=-G_o(x) u_2+F_o(x)+\delta_2\end{gathered}$                (27)

where, G_o(x) is nominal term with G_o(x)>0 $\forall x$, and $\delta_2$ is the perturbation term of the lower sub system.

$\delta_2=\Delta F(x)-\Delta G(x) u_2$                 (28)

Let the control law be taken as:

$u_2=G_o(x)^{-1}\left(u_{2 n}+u_{2 s}\right)$            (29)

and select the sliding variable of ISMC be chosen as

$s_2=s_{2 o}(e)+z_2$              (30)

Then, for

$s_{2 o}=e_3$                (31)

with z₂(0)=-e₃(0).

We obtain

$s_2=e_3+z_2$              (32)

where, z₂ is the integral term for the second controller.

The derivative of s₂ is:

$\dot{s}_2=\dot{e}_3+\dot{z}_2$              (33)

By substituting Eq. (27) into Eq. (33), we obtain

$\dot{s}_2=-\left(u_{2 n}+u_{2 s}\right)+F_o(x)+\delta_2+\dot{z}_2$              (34)

Let the integral part derivative be defined as

$\dot{z}_2=u_{2 n}-F_o(x)$               (35)

Accordingly, Eq. (34) becomes

$\dot{s}_2=-u_{2 s}+\delta_2$               (36)

Then, applying the equivalent control [33]

$\dot{s}_1=0$

Leads to                                                                  (37)

$\left[u_{2 s}\right]_{e q}=\delta_2$

As a result, Eq. (27) can be written as:

$\begin{gathered}\dot{e}_2=e_3 \\ \dot{e}_3=-u_{2 n}+F_o(x)\end{gathered}$              (38)

where, u_2n can be selected a linear controller as in Eq. (39) which grantees the asymptotically stability of the origin (e₂, e₃)=(0,0).

$u_{2 n}=F_o(x)+\alpha_2 e_2+\alpha_3 e_3$           (39)

where, α₂, α₃ are positive constants. Consequently:

$\dot{z}_2=\alpha_2 e_2+\alpha_3 e_3$            (40)

Finally, the discontinuous control term is given by:

$u_{2 s}=k_2 \operatorname{sign}\left(s_2\right)$           (41)

Then

$u_2=G_o(x)^{-1}\left(F_o(x)+\alpha_2 e_2+\alpha_3 e_3+k_2 \operatorname{sign}\left(s_2\right)\right)$               (42)

Remark (5): According to remark (2) and Eqns. (24) and (42) that u₁&u₂ are a positive quantity, then u₁&u₂ must follow the following rule:

$u_i=\left\{\begin{array}{cl}+v e & \text { if } u_i>0 \\ 0 & \text { if } u_i \leq 0\end{array}\right.$ where $i=1,2$           (43)

Remark (6): To use the barrier function instead of discontinuous term in the control law, the following PBFs h_ps(x) will be used [17]:

$u_s(s)=h_{p s}(s) * \operatorname{sign}(s)=\frac{s}{\epsilon-|s|}$               (44)

where, u_s(s) is a differentiable function of s(e).

Therefore, using the PBFs for u_s as in Eq. (44), u₁ and u₂ become:

$\begin{gathered}u_1=\frac{-1}{-a_{2 o} x_1+a_{3 o}}\left(\alpha_1 e_1+\frac{s_1}{\epsilon_1-\left|s_1\right|}\right) \\ u_2=G_o(x)^{-1}\left(F_o(x)+\alpha_2 e_2+\alpha_3 e_3+\frac{s_2}{\epsilon_2-\left|s_2\right|}\right)\end{gathered}$               (45)

where, $\epsilon_1$ is a constant, while $\epsilon_2$ is a time variable function which given by:

$\epsilon_2=\epsilon_{21} e x p^{-\lambda t}+\epsilon_{22}\left(1-e x p^{-\lambda t}\right)$              (46)

where, ϵ₂₁>ϵ₂₂>0.

Remark (7): Since u₁ is a differentiable function according to Eq. (45) $\forall e_1>0$, the time derivative of u₁ can be estimated using ACSMD.

5.3 Optimal of a linear control law

To guarantee an optimal ISMC, the LQR will be used to design a linear control for the control law in Eq. (39) [19]. The state space of the lower subsystem in Eq. (38) is described as:

$\dot{e}=\left[\begin{array}{ll}0 & 1 \\ 0 & 0\end{array}\right] e-\left[\begin{array}{l}0 \\ 1\end{array}\right] u_{2 n}$               (47)

The pair (A, B) is a controllable pair. Consequently, the linear control law which is described in Eq. (39) is derived from minimizing the cost function that mentioned in Eq. (3).

5.4 Chattering free sliding mode differentiator

In Eq. (32), the first step to implement the second controller u₂ is the constraction of the sliding variable s₂. Since e₂ is uncertain function, so we need to estimate its derivative value using an observer as mentioned in Remark (4). The second state e₃ is the time derivative of e₂. The ACSMD is proposed here to obtain e₃ which is introduced in studies [34, 35]. ACSMD is given by

$\left.\begin{array}{c}\sigma=e_2+\phi \\ \eta=\frac{2 \alpha_4}{\pi} \tan ^{-1}(\mu \sigma) \\ \dot{\phi}=-\eta \\ \dot{v}=\frac{1}{\tau}(-v+\eta)\end{array}\right\}$              (48)

where, σ is the sliding mode differentiator variable, α₄ and μ are the differentiator parameters. The fourth equation in Eq. (48), $\dot{v}=\frac{1}{\tau}(-v+\eta)$, is a low pass filter (LPF) with time constant $\tau$. The output of the LPF, v, is the estimated value of e₃. According to the study [34], the bound on the steady state estimation error is given by:

$\left|v(t)-e_3(t)\right| \leq \frac{2}{\tau \mu} \tan \left(\frac{\pi}{2 \alpha_4}\left|e_3\right|_{\max }\right)$                 (49)

where, the initial conditions ν(t_o)=0, ϕ(t_o)=-e₂(t_o), and α₄>|e₃|.