Unveiling overlooked aspects of model predictive control for building air conditioning systems

Abstract

The ongoing research on model predictive control (MPC) for building air conditioning systems predominantly centers on improving the predictive capabilities of system models. In this paper, the impacts of three additional pivotal factors on MPC performance are assessed by examining a generic MPC design for a typical variable air volume (VAV) system that serves large commercial buildings. The three factors encompass the nuanced reformulation of optimization, the judicious relaxation of constraints, and the meticulous tuning of parameters. Detailed case studies with an integrated Modelica and EnergyPlus model of the US Department of Energy's Commercial Reference Building are conducted. The results confirm that the optimization formulation, along with relaxation methods, significantly affects MPC performance in terms of energy savings, zonal thermal comfort level, and computational demand. They also reveal that the impact of the MPC control parameters on the energy savings and thermal comfort may vary by season and can be non-monotonic.

KEYWORDS:

• Model predictive control
• optimization formulation
• constraint relaxation
• parameter tuning
• Modelica
• EnergyPlus

1. Introduction

Building operations encompass a multitude of objectives ranging from the enhancement of indoor air quality, provision of thermal comfort, and maximization of energy efficiency. This diversity gives rise to multifaceted challenges in achieving optimal building performance. Recently, model prediction control (MPC) (Allgöwer and Zheng 2012) has been emerging as a comprehensive approach that can effectively address the intricate complexity of building operations (Drgoňa et al. 2020). MPC orchestrates a sequence of control actions by optimizing an objective function across a future time horizon, followed by execution of the initial action. In comparison with the prevalent rule-based control that predominates in contemporary buildings, MPC has two major advantages. First, it reconciles conflicting control objectives through co-optimization, thus deriving control actions for balanced operational goals. Second, it foresees upcoming system conditions through prediction, thus refining control actions for optimal building performance.

Although pilot demonstrations (Drgoňa, Picard, and Helsen 2020) have revealed promising results, the practical implementation of MPC into real-world scenarios is contingent upon a multitude of factors. In the literature, it has been commonly believed that the predictive accuracy and computational efficiency of building system models hold paramount importance for the performance of MPC. Hence, considerable research efforts have been dedicated to reducing modeling error and computational demand of building system models. For example, random forests (Smarra et al. 2018) and latent heat information (Raman et al. 2020) are used to improve the predictive accuracy of the building system models. Linearization (Yang and Wan 2022), model order reduction (Chen and Li 2021), and rule extraction(Yang et al. 2023) are used to reduce the computation time for executing the building system models.

However, the predictive accuracy and computational efficiency of building system models are not the sole factors influencing MPC performance. It turns out that the intricate processes of formulating, solving, and configuring MPC also warrant examination. During MPC formulation, it becomes imperative to incorporate the domain knowledge of the studied system. This knowledge should be embedded into the definition of optimization objectives and constraints, which is referred to as optimization reformulation. Following the optimization formulation, constraint relaxation must be undertaken for effectively solving the optimization problem with high complexity. Constraint relaxation aims to reduce the complexity of the original problem and is often achieved by either removing certain constraints or converting hard constraints into soft ones. Lastly, to achieve better control performance, fine-tuning parameters can be attempted for different MPC configurations. Prediction horizon and control interval are the two most often adjusted parameters, which define how far into the future the control considers and how often the control actions are updated, respectively.

The use of optimization reformulation, constraint relaxation, and parameter fine-tuning for improving MPC performance has been recognized across various domains, including applications in industrial electronics (Bolognani et al. 2008), chemical reactors (Mhaskar 2006), power systems (Valverde and Van Cutsem 2013), spacecraft (Dang et al. 2020), unmanned aerial vehicles (Michel et al. 2017), and robotic manipulators (Elsisi et al. 2021). However, these crucial factors have not received adequate attention in the design of MPC for building systems as only a limited number of studies have considered MPC parameter tuning for building applications. For example, the cost and the computing time of an MPC for a residential house with varying prediction horizons was studied (Lefort et al. 2013). The effect of the prediction horizon and control interval on the energy performance of an MPC was examined (Gholamibozanjani et al. 2018). In previous studies (Huang, Lin, Chinde et al. 2021), we scrutinized the impacts of control interval on MPC performance for a commercial building in terms of cost savings and computational time. While quantitatively describing the sensitivity of the control performance to the selected MPC parameters, those studies ignore the impacts from the operational conditions of buildings, such as changing weather conditions, leading to partial capture of this sensitivity.

In this paper, we unveil the impacts of optimization reformulation, constraint relaxation, and parameter tuning on MPC performance for building systems. We consider a generic MPC for a typical variable air volume (VAV) system of large commercial buildings to reduce overall energy consumption and enhance zonal thermal comfort. We then define various customization options to explore the effects of those three factors, which specifically include

• two optimization reformulations for eliminating simultaneous cooling and heating and avoiding unrealistic fan operations, respectively;

• two options for relaxing the thermal comfort constraint;

• five values for the prediction horizon; and

• three values of the control interval.

Detailed case studies with an integrated Modelica (Fritzson and Engelson 1998) and EnergyPlus (Crawley et al. 2001) model of the US Department of Energy's (DOE's) Commercial Reference Building (Deru et al. 2011) (hereinafter referred to as Reference Building) are conducted to quantitatively assess the impacts of these customization options.

The remainder of this paper is organized as follows. In Section 2, we describe the typical configuration of the studied VAV system. We then propose a new MPC for this VAV system and explore different options when designing this MPC in Section 3. After that, we discuss how to evaluate the proposed MPC with the Reference Building in Section 4. In Section 5, we present the results from the MPC evaluation and provide insights on how to interpret those results. Conclusions can be found in Section 6.

2. Studied system

Figure 1 shows a schematic of a typical VAV system that serves several zones in a building. In this system, there are multiple air handling unit (AHU), terminals, a duct network, a chiller, and a boiler. An AHU generates cooling air and distributes the air among several zones via the duct network. The AHU consists of a fan, a cooling coil, a chilled water valve, and an outdoor damper. The fan circulates air between each zone and the AHU. The cooling coil removes heat from the air with cold water from the chiller. The chilled water valve adjusts the water flow rate through the cooling coil. The outdoor air damper adjusts the flow rate of the outdoor air. Each terminal serves one zone and adjusts the supply air flow rate via a damper and heats the supply air with a reheat coil, respectively. The reheat coil connects to the boiler, and there is a valve that modulates the flow rate of the hot water through the heating coil. The duct network distributes the supply air among terminals. It also collects the return air from all the zones and determines how much of the return air is circulated back to the AHU with one exhaust damper. The chiller and the boiler produce chilled water and hot water, respectively.

Figure 1. Typical configuration of a variable air volume system.

The operation of the VAV system is modulated by 6 local controllers. Specifically, controller #1 adjusts the position of the chilled water valve to maintain the temperature of the cooling air at or less than a predefined setpoint, namely the discharge air temperature setpoint. Controller #2 adjusts the position of the outdoor/exhaust/return air damper to maintain the temperature of the air entering the cooling coil at or less than a predefined setpoint, namely the mixed air temperature setpoint. Controller #3 adjusts the speed of the fan to maintain the static pressure in the duct network at a predefined setpoint. Controller #4 determines the desired supply air flow rate, namely the supply airflow setpoint, and the desired supply air temperature, namely the supply temperature setpoint, based on the zone temperature and thermostat settings. Controller #5 modulates the position of the damper in each terminal to maintain the supply airflow rate at the supply airflow setpoint. Controller #6 adjusts the position of the reheat coil to maintain the supply air temperature at the supply temperature setpoint.

3. Model predictive control design

We consider a procedure for designing MPC, as illustrated in Figure 2. This procedure consists of four steps. Step 1, Problem Statement, defines the control objectives and sets up the corresponding optimization problem to achieve those objectives. Steps 2 and 3 aim to solve the problem defined in Step 1. Specifically, Step 2, Optimization Reformation, modifies the objective function and/or adds additional constraints to help the optimization solver find the true optimum. Step 3, Constraints Relaxation, extends the search space by modifying the constraints from Steps 1 and 2. Step 4, Control Parametrization, assigns the values for the control parameters so that the MPC can be deployed. In the following, we elaborate how we execute the procedure for the VAV system and different options in Step 2, 3 and 4.

Figure 2. Procedure for designing model predictive control.

3.1. Problem statement

When operating the studied VAV system, there are two objectives:

• Reducing the energy consumption of the VAV system. More precisely, we aim to reduce the energy consumption of chillers, fans, and boilers.

• Improving the thermal comfort level of the zones served by the VAV system. We define the thermal comfort based on the zone temperature and keep the zone temperature within predefined limits.

Based on the control objectives, we define the following optimization problem: (1a) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t,\end{array}$(1a) (1b) $\begin{array}{rl}\mathrm{subject}\mathrm{to}& x^k=f(x^{k-1},u^k,d^k),\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(1b) (1c) $\begin{array}{rl}& x^0={\hat{x}}^0,\end{array}$(1c) (1d) $\begin{array}{rl}& x^k\in {\mathcal{X}}^k,\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(1d) (1e) $\begin{array}{rl}& u^k\in {\mathcal{U}}^k,\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(1e) (1f) $\begin{array}{rl}& u^k=u^{\mathrm{Mj}+1},\mathrm{\forall }k\in {\mathbb{N}}_{\mathrm{Mj}+2}^{M(j+1)},\\ & j\in {\mathbb{N}}_0^{\frac{N}{M}-1},\end{array}$(1f) where k and N are the index of the discrete time step and the number of time steps within one optimization period, respectively. $P^k$, $x^k$, $u^k$, and $d^k$ denote the power of the studied system, the system states, the control decision variables, and exogenous inputs, respectively, at time step k. ${\hat{x}}^k$ denotes the measured values of the system states at time step k. ${\mathcal{X}}^k$ and ${\mathcal{U}}^k$ are the set of accepted values for $x^k$ and $u^k$, respectively. M is the number of time steps when the values of $u^k$ are unchanged.

For the studied VAV system, $P^k$ is calculated by: (2) $P^k=P_{\mathrm{ch}}^k+P_{\mathrm{bo}}^k+\sum _{j=i}^n{P}_{\mathrm{fa},i}^k,$(2) where $P_{\mathrm{ch}}^k$, $P_{\mathrm{bo}}^k$, and $P_{\mathrm{fa},i}^k$ are the power of the chiller, the power of the boiler, and the fan power of AHU i, at time step k. n denotes the number of AHUs in the studied VAV system. $P_{\mathrm{ch}}^k$ is calculated by (3) $P_{\mathrm{ch}}^k=a_0+a_1{Q}_{\mathrm{co}}^k+a_2(Q_{\mathrm{co}}^k{)}^2+a_3(Q_{\mathrm{co}}^k{)}^3,$(3) where $a_0$, $a_1$, $a_2$, and $a_3$ are constants. $Q_{\mathrm{co}}^k$ is the cooling load at time step k and is calculated by (4) $\begin{array}{rl}{Q}_{\mathrm{co}}^k& =\frac{c_p}{\sigma }\sum _{i=1}^n\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k\\ & ({\zeta }_i^k{T}_o^k+(1-{\zeta }_i^k)\frac{\sum _{j=1,z_i}{\dot{m}}_{i,j}^k{T}_{i,j}^k}{\sum _{j=1,z_i}{\dot{m}}_{i,j}^k}-T_{s,i}^k),\end{array}$(4) where $c_p$ is the specific heat of air. σ is the sensible heat ratio. ${\dot{m}}_{i,j}^k$ and $T_{i,j}^k$ are the supply flow rate and the zone temperature of zone j, which is served by AHU i, respectively, at time step k. $z_i$ denotes the number of zones served by AHU i. ${\zeta }_i^k$ and $T_{s,i}^k$ are the outdoor air damper position and the supply air temperature of AHU i, respectively, at time step k. $T_o^k$ is the outdoor air temperature at time step k.

$P_{\mathrm{bo}}^t$ is calculated by (5) $P_{\mathrm{bo}}^k=b_0{c}_p\sum _{i=1}^n\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{s,i}^k),$(5) where $b_0$ is a constant. $T_{s,i,j}^k$ is the supply air temperature of zone j, which is served by AHU i, at time step k.

$P_{\mathrm{fa},i}^t$ is calculated by (6) $\begin{array}{rl}{P}_{\mathrm{fa},i}& =c_{0,i}+c_1\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k+c_{2,i}{\left(\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k\right)}^2\\ & +c_{3,i}{\left(\sum _{j=1,z_i}{\dot{m}}_{i,j}^k\right)}^3,\end{array}$(6) where $c_{0,i}$, $c_{1,i}$, $c_{2,i}$, and $c_{3,i}$ are constants.

Note that in Equations (4(4) $\begin{array}{rl}{Q}_{\mathrm{co}}^k& =\frac{c_p}{\sigma }\sum _{i=1}^n\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k\\ & ({\zeta }_i^k{T}_o^k+(1-{\zeta }_i^k)\frac{\sum _{j=1,z_i}{\dot{m}}_{i,j}^k{T}_{i,j}^k}{\sum _{j=1,z_i}{\dot{m}}_{i,j}^k}-T_{s,i}^k),\end{array}$(4) ), (5(5) $P_{\mathrm{bo}}^k=b_0{c}_p\sum _{i=1}^n\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{s,i}^k),$(5) ), and (6(6) $\begin{array}{rl}{P}_{\mathrm{fa},i}& =c_{0,i}+c_1\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k+c_{2,i}{\left(\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k\right)}^2\\ & +c_{3,i}{\left(\sum _{j=1,z_i}{\dot{m}}_{i,j}^k\right)}^3,\end{array}$(6) ), $T_o^k$ is one exogenous input, and $T_{i,j}^k$ is one system state variable, while $T_{s,i,j}^k$, ${\dot{m}}_{i,j}^k$, ${\gamma }_i^k$, and $T_{s,i}^k$ are control decision variables. In other words, $d^k$, $x^k$, and $u^k$ are defined as (7a) $\begin{array}{rl}& d^k=T_o^k,\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(7a) (7b) $\begin{array}{rl}& x^k=\{T_{i,j}^k|i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(7b) (7c) $\begin{array}{rl}& u^k=\{T_{s,i,j}^k,{\dot{m}}_{i,j}^k,{\zeta }_i^k,T_{s,i}^k|i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(7c) $T_{i,j}^k$ is calculated by (8) $T_{i,j}^k={\alpha }_{i,j}{T}_{i,j}^{k-1}+{\beta }_{i,j}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{i,j}^{k-1})+{\gamma }_{i,j}{T}_o^k+{\dot{q}}_{i,j},$(8) where ${\alpha }_{i,j}$, ${\beta }_{i,j}$, ${\gamma }_{i,j}$, and ${\dot{q}}_{i,j}$ are constants.

In addition, ${\mathcal{X}}^k$ and ${\mathcal{U}}^k$ are defined with (9a) $\begin{array}{rl}& {\mathcal{X}}^k=\{T_{i,j}^k|T_{i,j,l}^k\le T_{i,j}^k\le T_{i,j,h}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9a) (9b) $\begin{array}{rl}& {\mathcal{U}}^k=\{u_s|u_{s,\mathrm{low}}^k\le u_s\le u_{s,\mathrm{hig}}^k,u_s\in u^k\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9b) where $T_{i,j,l}^k$ and $T_{i,j,h}^k$ are the lower limit and the upper limit of $T_{i,j}^k$, respectively. $u_{s,\mathrm{low}}^k$ and $u_{s,\mathrm{hig}}^k$ are the lower limit and the upper limit of element s in $u^k$, respectively.

When implementing $u^k$ to the studied VAV system,

• $T_{s,i,j}^k$ and ${\dot{m}}_{i,j}^k$ are used to overwrite the output of Controller #4, the supply temperature setpoint and the supply airflow setpoint, respectively.

• ${\zeta }_i^k$ is used to overwrite the output of Controller #2.

• $T_{s,i}^k$ is used as the discharge air temperature setpoint and the mixed air temperature setpoint.

3.2. Optimization reformulation

We consider two modifications to the original formulation defined in Equation (1a(1a) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t,\end{array}$(1a) ). The first modification is to avoid simultaneous cooling and heating. For a VAV system, simultaneous cooling and heating usually occur when the supply airflow rate and the supply air temperature of the thermal zones are both set to be high. In such a case, the cooling power provided by the chiller is canceled out by the reheat coils. One can also explicitly eliminate simultaneous cooling and heating by changing the objective function (1a(1a) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t,\end{array}$(1a) ) to (10) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +\omega \sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k),\end{array}$(10) where ω is a weight factor. The second modification is to avoid the unrealistic operation of fans. Based on Equation (6(6) $\begin{array}{rl}{P}_{\mathrm{fa},i}& =c_{0,i}+c_1\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k+c_{2,i}{\left(\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k\right)}^2\\ & +c_{3,i}{\left(\sum _{j=1,z_i}{\dot{m}}_{i,j}^k\right)}^3,\end{array}$(6) ), $P_{\mathrm{fa},i}$ can be larger than 0 even if the flow rate is 0, which is not realistic. Usually, the coefficients in Equation (6(6) $\begin{array}{rl}{P}_{\mathrm{fa},i}& =c_{0,i}+c_1\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k+c_{2,i}{\left(\sum _{j=1}^{z_i}{\dot{m}}_{i,j}^k\right)}^2\\ & +c_{3,i}{\left(\sum _{j=1,z_i}{\dot{m}}_{i,j}^k\right)}^3,\end{array}$(6) ) are obtained through linear regressions, which might not guarantee that $c_{0,i}$ is equal to 0. One can avoid such unreasonable fan operations by adding the following constraint: (11) $c_{0,i}=0,\mathrm{\forall }i\in {\mathbb{N}}_1^n.$(11)

3.3. Constraint relaxation

The constraint relaxation is essentially to find a new problem formulation with a solution set that is the same as or similar to the original problem formulation but much easier to solve. In the following section, we elaborate on how we modify constraints and add penalty constraints for the proposed MPC.

In Equation (1a(1a) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t,\end{array}$(1a) ), there are two types of constraints. The first type is for the independent variables, such as the ones defined in Equation (1e(1e) $\begin{array}{rl}& u^k\in {\mathcal{U}}^k,\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(1e) ). Those constraints are usually easy to handle and thus are not considered in the constraint relaxation. The second type of constraint is for the dependent variables, such as the ones defined in Equation (1d(1d) $\begin{array}{rl}& x^k\in {\mathcal{X}}^k,\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(1d) ). Those constraints are typically difficult to handle and need to be modified. Specifically, we add a term, called a penalty function, to the objective function that consists of a penalty parameter multiplied by a measure of violation of the second constraint. There are two options for implementing this penalty function. In the first option, the objective function (1a(1a) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t,\end{array}$(1a) ) is changed to (12) $\underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t+\mathrm{\omega \xi }{\nu }^2,$(12) where ξ is a penalty coefficient and ν is a slack variable. ξ is used to adjust the trade-off between energy saving and thermal comfort; ν is calculated by the solvers of the optimization problem. Equation (9a(9a) $\begin{array}{rl}& {\mathcal{X}}^k=\{T_{i,j}^k|T_{i,j,l}^k\le T_{i,j}^k\le T_{i,j,h}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9a) ) is then changed to (13) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-\nu \le T_{i,j}^k\le T_{i,j,h}^k+\nu ,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(13) In the second option, the objective function (1a(1a) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t,\end{array}$(1a) ) is changed to (14) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +\omega \sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2,\end{array}$(14) where ${\nu }_{i,j}^t$ is a slack variable for zone j, which is served by AHU i at time t. Similar to ν, ${\nu }_{i,j}^t$ is also calculated by the solvers of the optimization problem. Equation (9a(9a) $\begin{array}{rl}& {\mathcal{X}}^k=\{T_{i,j}^k|T_{i,j,l}^k\le T_{i,j}^k\le T_{i,j,h}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9a) ) is then changed to (15) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-{\nu }_{i,j}^k\le T_{i,j}^k\le T_{i,j,h}^k+{\nu }_{i,j}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(15)

3.4. Control parametrization

In this study, MPC parameters refer to the prediction horizon and the control interval, which are defined as $N\mathrm{\Delta t}$ and $M\mathrm{\Delta t}$, respectively. The prediction horizon specifies how far we should look ahead, while the control interval measures how frequently we should update the control decision variables. When determining the values for MPC parameters, one needs to consider impacts from two aspects. On one hand, a longer prediction horizon and a smaller control interval might help the MPC to better achieve the control objectives. On the other hand, a longer prediction horizon might cause inaccurate predictions from the system models because of error accumulation (Huang, Lin, Chinde et al. 2021). Likewise, a smaller control interval might result in higher computational demand. Therefore, quantitatively evaluating the impacts of the MPC parameters is indispensable to identifying the appropriate values for those parameters. However, most of the previous study tends to select arbitrary values for MPC parameters. Based on the survey conducted by Serale et al. (2018), the value of the prediction horizon is usually selected to be larger than 1 h, while the control interval is usually selected from a range of 5 to 30 min.

4. Case study

This section details how the proposed MPC was evaluated over the Reference Building. We first introduce the Reference Building and how this building is modeled with EnergyPlus and Modelica. We then introduce how the evaluation is implemented and detail the settings of the evaluation. Finally, we discuss how the evaluation results are analyzed.

4.1. Reference building

When evaluating the proposed MPC, it is essential to include a representative operation condition of the VAV system. To fulfill this requirement, the Reference Building is employed. The Reference Building was developed by DOE and three of its national laboratories to represent the average characteristics of commercial building groups, categorized by building type and climate zone, in the United States. Specifically, the Reference Building for the large office buildings in Chicago is used in this study.

This higher-fidelity building model represents a VAV system that serves a typical high-rise building. This Reference Building has 12 floors with 5 zones per floor. In this building, the occupant density, the lighting power density, and the plug and process power density are $18.6{\mathrm{m}}^2/\mathrm{person}$, $10.8\mathrm{W}/{\mathrm{m}}^2$ and $10.8\mathrm{W}/{\mathrm{m}}^2$, respectively. The schedules of the internal load are based on ASHRAE (1989). The VAV system that serves this building contains 12 AHUs. Each AHU serves one floor and is connected to five VAV terminals.

The Reference Building is not a real building; thus, a simulation model of this building is used for evaluating the proposed MPC. In this study, we leverage a higher-fidelity building model from our previous study in Huang et al. (2018). When implementing this higher-fidelity building model, we employed EnergyPlus and Modelica Buildings Library (Wetter et al. 2014). More precisely, we modeled the building envelope and internal heat gain with EnergyPlus and approximated the VAV system and its associated control with Modelica Buildings Library. The data exchange between the EnergyPlus model and the Modelica model is realized via a functional mockup interface (Nouidui, Wetter, and Zuo 2014) and is illustrated in Figure 3. Note that the EnergyPlus model is derived from research associated with the Reference Building and has been reviewed extensively by building industry experts. The Modelica model of the physical systems was built based directly on existing Modelica libraries. As validations have already been performed for those libraries (Wetter et al. 2014), we do not perform a further valuation on this model in this study. On the other hand, for the Modelica model for the control systems, we have performed validations to make sure the systems can generate expected outputs when inputs are given. More details of the validation can be found in Huang et al. (2018) and Huang, Zuo, Vrabie et al. (2021).

Figure 3. Cosimulation between EnergyPlus and Modelica.

To facilitate the MPC evaluation, we implement a baseline control, as summarized in the following.

• Chiller plant and boiler plant,

  The operating number of chillers and boilers is controlled based on cooling load and heating load with state machines, respectively. The temperature of the chilled water leaving the chiller plant and the temperature of the hot water leaving the boiler plant are maintained based on predefined setpoints. More details can be found in Huang, Zuo, Vrabie et al. (2021).

• AHU,

  The supply air temperature setpoint is maintained as ${12.89}^{\circ }$C. The outdoor air damper position is controlled based on a fixed dry-bulb temperature economizer control (ASHRAE 2016) in which the threshold outdoor temperature is set at ${15.56}^{\circ }$C. Specifically, when the outdoor temperature exceeds ${15.56}^{\circ }$C, the outdoor air damper position is fixed as 10%; Otherwise, the outdoor air damper position is modulated to maintain the mixed air temperature as ${12.89}^{\circ }$C while the minimum damper position is set to be 10%. Note that the outdoor air damper position is closed during the unoccupied period. The fan speed is continuously modulated to maintain a constant static pressure during the occupied period and cycling on during the rest of the day.

• Terminals,

  The heating coil position and the supply airflow setpoint are controlled based on a dual-max control sequence (ASHRAE 2018). The supply airflow setpoint is then used by a feedback control to modulate the damper position.

4.2. Evaluation execution and settings

To facilitate the execution of the MPC test, we set up a cosimulation between the high-fidelity building model and the MPC. This cosimulation is based on an assumption that the CPU time cost by MPC is much lower than the control interval, that is, ${\delta }_i$ ¡ $\mathrm{\Delta }t$ ($i\in {\mathbb{N}}_1^{\mathrm{NM}}$). At the beginning of each control interval, the simulation pauses and sends the current measurement to trigger the MPC. The MPC then calculates the control actions based on this measurement and sends the control signals back to the simulation. After receiving the control signals, the simulation continues until it reaches the beginning of the next control interval. Note that among the decision variables, the outdoor damper position is implemented directly as the actuator status and the rest are implemented as setpoints for the local controller.

In this cosimulation, the MPC is implemented with JuMP (version 0.21.10) (Dunning, Huchette, and Lubin 2017), which is a domain-specific modeling language for mathematical optimization embedded in Julia (version 1.0.5) (Bezanson et al. 2017) and solved by IPOPT (version:0.7.0) (Wächter 2009). The cosimulation environment is implemented in the BOPTEST framework (Blum et al. 2021). We set the maximum iteration number of the MPC solver to be 1,000. If the iteration number exceeds the maximum iteration number, then the baseline control will be implemented. This cosimulation run on a DELL Latitude E7470 Laptop PC (Win10 OS, Intel^®CoreTM i7-6600, 16.0 GB RAM).

When executing the simulation, we consider two testing periods–one summer week (August 26–30) and one mild week (April 8–12)–and use the typical meteorological year (TMY) data (Hall et al. 1978) of Chicago, Illinois, United States, as the major input for those testing periods. When testing the MPC, we use the TMY data directly as the ideal weather forecast. In addition, for each testing period, one dedicated training period is used to generate training datasets for the system models. Figure 4 illustrates the outdoor air temperature during the testing period and the training period in the summer. To enhance data richness, we introduce time-varying thermostat setpoints for each zone to better excite the baseline control, as illustrated in Figure 5.

Figure 4. Outdoor air temperature during the testing period and the training period in the summer.

Figure 5. Excitation signals of the training period in the summer.

4.3. Performance metrics

We evaluate the performance of MPC from three aspects:

• how much the proposed MPC can reduce the energy consumption of the VAV systems;

• how well the proposed MPC can improve thermal comfort in the zones served by the VAV systems; and

• how fast the associated optimization problem can be solved.

When evaluating the MPC, we compare its performance to that of the baseline control with the following performance metrics:

• Electricity Consumption Saving ($\mathrm{\%}$). This metric characterizes the enhancement of the energy performance of the VAV system, which is defined as the percentage of saving in terms of electricity consumption; that is, (16) $\mathrm{\%}{C}_e=\frac{E_{\mathrm{base}}-E_{\mathrm{mpc}}}{E_{\mathrm{base}}},$(16) where $E_{\mathrm{base}}$ and $E_{\mathrm{mpc}}$ are the electricity consumption by the fans and the chiller with the baseline control and with MPC, respectively.

• Gas Consumption Saving ($\mathrm{\%}$). This metric characterizes the enhancement of the energy performance of the VAV system, which is defined as the percentage of saving in terms of gas consumption; that is, (17) $\mathrm{\%}{C}_g=\frac{G_{\mathrm{base}}-G_{\mathrm{mpc}}}{G_{\mathrm{base}}},$(17) where $G_{\mathrm{base}}$ and $G_{\mathrm{mpc}}$ are the gas consumption by the fans and the boiler with the baseline control and with MPC, respectively.

• Discomfort Reduction ($\mathrm{K}\cdot \min$). This metric characterizes the enhancement of building operations, which is defined as the reduction in the discomfort level; that is, (18) $\mathrm{DR}=\mathrm{Di}{s}_{\mathrm{base}}-\mathrm{Di}{s}_{\mathrm{mpc}},$(18) where $\mathrm{Di}{s}_{\mathrm{base}}$ and $\mathrm{Di}{s}_{\mathrm{mpc}}$ are the discomfort level without and with MPC, respectively. The discomfort level is calculated by (19) $\mathrm{Dis}=\sum _{i=1}^n\sum _{j=1}^{z_i}g(T_{i,j}^t)\frac{\mathrm{\Delta }t}{N},$(19) where the function $g()$ is defined as (20) $g(T_{i,j}^t)=\{\begin{array}{ll}{T}_{i,j}^t-T_{\mathrm{hig},i,j}^t,& \mathrm{if}{T}_{i,j}^t>T_{\mathrm{hig},i,j},\\ T_{\mathrm{low},i,j}^t-T_{i,j}^t,& \mathrm{if}{T}_{i,j}^t<T_{\mathrm{hig},i,j},\\ 0,& \mathrm{otherwise},\end{array}$(20) where $T_{\mathrm{low},i,j}^t$ and $T_{\mathrm{hig},i,j}^t$ are the lower and higher bounds for $T_{i,j}^t$.

• Computational Time (s). This metric characterizes the numerical performance of the MPC, which is defined as the sum of the computational time of the MPC runs; that is, (21) $\mathrm{CT}=\sum _{s=1}^M{t}_c^s,$(21) where $t_c^s$ is the CPU time used by MPC run s while M is the number of MPC runs during the testing period.

5. Results and discussion

5.1. MPC design assessment

To facilitate the assessment of different options in the optimization reformulation and the constraint relaxation, we consider four MPC designs, as illustrated in Table 1. Case 1 and case 2 use different methods for handling the thermal comfort constraints; Unlike case 1, case 3 doesn't consider simultaneous cooling and heating in the objective function; case 4 drops the constraint for avoiding the unreasonable fan operations in case 1.

Table 1. Model predictive control design.

Case	Objective function	Constraints
1	Equation (22(22) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +{\omega }_1\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k),\end{array}$(22) )	Equations (8(8) $T_{i,j}^k={\alpha }_{i,j}{T}_{i,j}^{k-1}+{\beta }_{i,j}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{i,j}^{k-1})+{\gamma }_{i,j}{T}_o^k+{\dot{q}}_{i,j},$(8) ), (1c(1c) $\begin{array}{rl}& x^0={\hat{x}}^0,\end{array}$(1c) ), (15(15) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-{\nu }_{i,j}^k\le T_{i,j}^k\le T_{i,j,h}^k+{\nu }_{i,j}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(15) ), (9b(9b) $\begin{array}{rl}& {\mathcal{U}}^k=\{u_s|u_{s,\mathrm{low}}^k\le u_s\le u_{s,\mathrm{hig}}^k,u_s\in u^k\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9b) ), (1f(1f) $\begin{array}{rl}& u^k=u^{\mathrm{Mj}+1},\mathrm{\forall }k\in {\mathbb{N}}_{\mathrm{Mj}+2}^{M(j+1)},\\ & j\in {\mathbb{N}}_0^{\frac{N}{M}-1},\end{array}$(1f) ), (11(11) $c_{0,i}=0,\mathrm{\forall }i\in {\mathbb{N}}_1^n.$(11) )
2	Equation (23(23) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t+{\omega }_1\xi {\nu }^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k).\end{array}$(23) )	Equations (8(8) $T_{i,j}^k={\alpha }_{i,j}{T}_{i,j}^{k-1}+{\beta }_{i,j}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{i,j}^{k-1})+{\gamma }_{i,j}{T}_o^k+{\dot{q}}_{i,j},$(8) ), (1c(1c) $\begin{array}{rl}& x^0={\hat{x}}^0,\end{array}$(1c) ), (13(13) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-\nu \le T_{i,j}^k\le T_{i,j,h}^k+\nu ,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(13) ), (9b(9b) $\begin{array}{rl}& {\mathcal{U}}^k=\{u_s|u_{s,\mathrm{low}}^k\le u_s\le u_{s,\mathrm{hig}}^k,u_s\in u^k\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9b) ), (1f(1f) $\begin{array}{rl}& u^k=u^{\mathrm{Mj}+1},\mathrm{\forall }k\in {\mathbb{N}}_{\mathrm{Mj}+2}^{M(j+1)},\\ & j\in {\mathbb{N}}_0^{\frac{N}{M}-1},\end{array}$(1f) ), (11(11) $c_{0,i}=0,\mathrm{\forall }i\in {\mathbb{N}}_1^n.$(11) )
3	Equation (24(24) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +\omega \sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2.\end{array}$(24) )	Equations (8(8) $T_{i,j}^k={\alpha }_{i,j}{T}_{i,j}^{k-1}+{\beta }_{i,j}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{i,j}^{k-1})+{\gamma }_{i,j}{T}_o^k+{\dot{q}}_{i,j},$(8) ), (1c(1c) $\begin{array}{rl}& x^0={\hat{x}}^0,\end{array}$(1c) ), (15(15) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-{\nu }_{i,j}^k\le T_{i,j}^k\le T_{i,j,h}^k+{\nu }_{i,j}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(15) ), (9b(9b) $\begin{array}{rl}& {\mathcal{U}}^k=\{u_s|u_{s,\mathrm{low}}^k\le u_s\le u_{s,\mathrm{hig}}^k,u_s\in u^k\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9b) ), (1f(1f) $\begin{array}{rl}& u^k=u^{\mathrm{Mj}+1},\mathrm{\forall }k\in {\mathbb{N}}_{\mathrm{Mj}+2}^{M(j+1)},\\ & j\in {\mathbb{N}}_0^{\frac{N}{M}-1},\end{array}$(1f) ), (11(11) $c_{0,i}=0,\mathrm{\forall }i\in {\mathbb{N}}_1^n.$(11) )
4	Equation (22(22) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +{\omega }_1\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k),\end{array}$(22) )	Equations (8(8) $T_{i,j}^k={\alpha }_{i,j}{T}_{i,j}^{k-1}+{\beta }_{i,j}{\dot{m}}_{i,j}^k(T_{s,i,j}^k-T_{i,j}^{k-1})+{\gamma }_{i,j}{T}_o^k+{\dot{q}}_{i,j},$(8) ), (1c(1c) $\begin{array}{rl}& x^0={\hat{x}}^0,\end{array}$(1c) ), (15(15) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-{\nu }_{i,j}^k\le T_{i,j}^k\le T_{i,j,h}^k+{\nu }_{i,j}^k,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(15) ), (9b(9b) $\begin{array}{rl}& {\mathcal{U}}^k=\{u_s|u_{s,\mathrm{low}}^k\le u_s\le u_{s,\mathrm{hig}}^k,u_s\in u^k\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9b) ), (1f(1f) $\begin{array}{rl}& u^k=u^{\mathrm{Mj}+1},\mathrm{\forall }k\in {\mathbb{N}}_{\mathrm{Mj}+2}^{M(j+1)},\\ & j\in {\mathbb{N}}_0^{\frac{N}{M}-1},\end{array}$(1f) )

(22) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +{\omega }_1\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k),\end{array}$(22) where ${\omega }_1$ and ${\omega }_2$ are weighting factors. (23) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t+{\omega }_1\xi {\nu }^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k).\end{array}$(23) (24) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +\omega \sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2.\end{array}$(24) Note that all the weight factors in the study are determined so that all the corresponding penalty terms are closer to $\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t$.

When evaluating the MPC designs in Table 1, we set the prediction horizon and the control interval to be 60 min and 5 min, respectively. For simplicity, we consider the first day of the mild week in the evaluation.

Table 2 summarizes the evaluation results:

• Equation (22(22) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +{\omega }_1\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k),\end{array}$(22) ) can better maintain the thermal comfort level than Equation (23(23) $\begin{array}{rl}\underset{u^1,\dots ,u^N}{\text{minimize}}& \sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t+{\omega }_1\xi {\nu }^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k).\end{array}$(23) ). In case 2, the MPC even increases the discomfort level compared with the baseline, while the discomfort level is reduced by 289.8 $\mathrm{K}\cdot min$ in case 1. This suggests that setting a dedicated, dynamic slack variable for each zone can better handle the thermal comfort constraint. As illustrated in Figure 6, in case 1, all the zone temperatures are within the comfort zones most of the time. In contrast, the zone temperatures in case 2 are beyond the comfort zones for most of the time, especially during the occupied period. This is due to the relaxation of the temperature constraints with Equation (13(13) $\begin{array}{rl}{\mathcal{X}}^k& =\{T_{i,j}^k|T_{i,j,l}^k-\nu \le T_{i,j}^k\le T_{i,j,h}^k+\nu ,i\in {\mathbb{N}}_1^n,j\in {\mathbb{N}}_1^{z_i}\},\\ & \mathrm{\forall }k\in {\mathbb{N}}_1^N.\end{array}$(13) ).

• Equation (22(22) $\begin{array}{rl}& \underset{u^1,\dots ,u^N}{\text{minimize}}\sum _{k=1}^N(P^k(x^k,u^k,d^k))\mathrm{\Delta }t\\ & +{\omega }_1\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}\xi ({\nu }_{i,j}^k{)}^2\\ & +{\omega }_2\sum _{k=1}^N\sum _{i=1}^n\sum _{j=1}^{z_i}(T_{s,i,j}^k-T_{s,i}^k),\end{array}$(22) ) reduces the chance for simultaneous heating and cooling. Compared with that in case 1, the electricity consumption and the gas consumption increase by 8% and 49%, respectively, in case 3. As shown in Figure 7, the supply flow rate setpoint conflicts with that of the supply temperature in case 3. The supply flow is sometimes set to be quite high, but the supply temperature is set to be high. In addition, there are times when supply flow is set to be quite high but the supply temperature is set to be higher. Those are typical ‘symptoms’ for simultaneous heating and cooling.

• Equation (11(11) $c_{0,i}=0,\mathrm{\forall }i\in {\mathbb{N}}_1^n.$(11) ) can significantly reduce electricity consumption. Compared with case 4, the MPC can reduce electricity consumption by 10% in case 1. As shown in Figure 8, the reduction in electricity consumption is due mainly to the reduction of the supply flow rate during the unoccupied period.

Figure 6. Zone temperatures in case 1 and case 2.

Figure 7. Supply flow rate and supply temperature of one perimeter zone in a typical floor.

Figure 8. Supply flow rate of the AHU that serves a typical floor.

Table 2. Evaluation results for model predictive control cases.

Case	Electricity consumption saving [%]	Gas consumption saving [%]	Discomfort reduction [$\mathrm{K}\cdot min$]	Computational time [s]
1	20.9	8.3	289.8	472.5
2	29.2	22.3	−3,932.4	553.1
3	12.1	−1.2	291.9	965.2
4	14.2	−54.4	235.4	581.1

5.2. Parameter sensitivity analysis

In this section, we evaluate how the MPC parameters affect control performance. We first study MPC performance under different prediction horizon lengths (the control interval is fixed as 5 min) for a mild day and a summer day, respectively. One can see the following from Table 3:

• For the mild day, the electricity consumption saving ratio and the gas consumption saving ratio don't change much when the prediction horizon increases. The electricity consumption saving ratio and the gas consumption saving ratio decrease when the prediction horizon increases for the summer day in general. It is interesting that the change of saving ratios by the prediction horizon is not monotonic.

• The discomfort reduction is not sensitive to the prediction horizon for both the mild day and the summer day. We did not see a tendency of change in discomfort reduction when the prediction horizon increased for the summer day. For the mild day, the discomfort reduction increases a little when the prediction horizon increases.

• The computational time increases when the prediction horizon increases for both the mild day and the summer day in general. With a larger prediction horizon, the number of decision variables increases, according to Equation (9b(9b) $\begin{array}{rl}& {\mathcal{U}}^k=\{u_s|u_{s,\mathrm{low}}^k\le u_s\le u_{s,\mathrm{hig}}^k,u_s\in u^k\},\mathrm{\forall }k\in {\mathbb{N}}_1^N,\end{array}$(9b) ), leading to a more complicated optimization problem. It is worth mentioning that the computational time does not increase in proportion to the prediction horizon. As shown in Figure 9, the distribution of computational time is quite uneven for MPC runs on the mild day. The outliers in Figure 9 explain the nonlinear increase of computational time by the prediction horizon.

Figure 9. Computational times of model predictive cotrol with different prediction horizons for the mild day.

Table 3. MPC evaluation results with different prediction horizons (control interval: 5 min, baseline electricity consumption: 9,980 kWh, baseline gas consumption: 53,1723 kJ).

Season	Prediction horizon $[min]$	Electricity consumption saving $[\mathrm{\%}]$	Gas consumption saving $[\mathrm{\%}]$	Discomfort reduction $[\mathrm{K}\cdot min]$	Computational time $[\mathrm{s}]$
summer	60	0.8	−3.2	277.1	454.0
	75	0.5	−6.4	269.2	663.0
	90	0.5	−5.7	278.1	790.7
	105	−1.9	−10.8	263.2	1,147.5
	120	−0.8	−4.8	278.7	1,122.0
mild	60	20.9	8.3	289.8	472.5
	75	20.6	8.7	288.6	633.5
	90	21.1	7.8	290.0	799.9
	105	21.5	7.8	290.6	1,026.3
	120	21.6	7.9	291.0	1,355.0

Table 4 illustrates the MPC performance under different control interval lengths (the prediction horizon is fixed at 60 min) for a mild day and a summer day, respectively. One can see the following:

• For the mild day, the electricity consumption saving ratio and the gas consumption saving ratio both decrease dramatically when the control interval increases for both the mild day and the summer day. This is because the MPC's ability to capture building dynamics decreases when the control interval increases.

• The discomfort reduction is decreased significantly when the control interval increases for both the mild day and the summer day. This means frequently updated control actions are necessary to maintain the desired thermal comfort level.

• The computational time decreases significantly when the control interval increases for both the mild day and the summer day. When the interval increases, the number of decision variables decreases, accelerating the optimization solving process.

Table 4. Model predictive control evaluation results with different control intervals (prediction horizon: 60 min).

Season	Control interval $[min]$	Electricity consumption saving $[\mathrm{\%}]$	Gas consumption saving $[\mathrm{\%}]$	Discomfort reduction $[\mathrm{K}\cdot min]$	Computational time $[\mathrm{s}]$
summer	5	0.9	−3.2	277.1	454.0
	15	−11.1	−8.9	258.8	139.5
	30	−33.4	−23.1	18.2	97.7
mild	5	20.9	8.3	289.8	472.5
	15	12.06	−1.1	215.0	205.4
	30	−5.9	−10.8	37.8	103.6

5.3. Detailed observation

In this section, we examine the details regarding the week-long simulations of the MPC operation to understand how the MPC achieves the control objectives. Note that in the week-long simulations, the control interval and the prediction horizon are set to be 5 min and 60 min, respectively, and the rest of the MPC implementation is set based on case 1. Table 5 summarizes the control performance of this MPC for the two typical weeks.

Table 5. Summary of model predictive control evaluation.

Season	Electricity consumption saving [%]	Gas consumption saving [%]	Discomfort reduction [$\mathrm{K}\cdot min$]	Computational time [s]
mild week	15.2	9.1	956.3	1,876.4
summer week	−1.1	−5.2	1,763.2	1,982.3

One can see the following from Table 5:

• The MPC results in 15% and 9% reductions in electricity consumption and gas consumption, respectively, in the mild week. However, in the summer week, this MPC does not result in any savings in energy consumption. This suggests that weather significantly affects the energy performance of the MPC. Note that the energy saving ratios for the entire mild week are lower than those of the first day of the mild week. This might be because the outdoor temperature during the second day, the third day, and the last day is generally higher than during the first day, as illustrated in Figure 10.

• This MPC can significantly improve the thermal comfort level for both the mild week and the summer week. Figure 11 shows the zone temperature under the MPC and the baseline during the summer week. We can see that improvements in comfort level usually occur when the operation mode of the VAV system changes from unoccupied to occupied. During the transition from unoccupied to occupied, the upper bound of the zone temperature decreases and the lower bound of the zone temperature increases.

• The computational time for the mild week is similar to that for the summer week. The average computational time for one MPC run is around 1.3 s, which is much lower than the control interval. This means that the assumption we made in Section 4.2 is valid.

Figure 10. Outdoor air temperature of the mild week.

Figure 11. Zone temperature during the mild week.

To understand how the reductions in electricity consumption and gas consumption are realized, we first plot the major operation variables of an AHU that serves a typical floor during the first day of the mild season. As shown in Figure 12:

• Compared to the baseline, the MPC can reduce the supply flow for most of the time, resulting in reduced fan power consumption. In addition, the MPC results in more hours of free cooling with the outdoor air, reducing the chiller consumption. This might partially explain why this MPC cannot reduce electricity consumption during the summer week: The outdoor temperature is too high for free cooling.

• The MPC turns on the fan 1 h earlier than the baseline control, causing additional electricity consumption during the occupied period. One hypothesis for this is that starting the fan earlier improves thermal comfort.

Figure 12. Operation of one air handling unit that serves the typical floor during the first day of the mild season.

To validate the hypothesis, we review the detailed results for individual VAV terminals. Figure 13 illustrates the operation of one VAV terminal that serves a core zone on the typical floor during the first day of the mild season. One can see that this MPC significantly reduces the amount of time that the zone temperature is above the upper limit in the earlier morning by starting the fan earlier.

Figure 13. Operation of one variable air volume terminal that serves a core zone in the typical floor during the first day of the mild season.

Starting the fan earlier can yield energy savings. Figure 14 shows the operation of one VAV terminal that serves a perimeter zone on the typical floor during the first day of the mild season. We can see that by starting the fan earlier, the MPC increases the zone temperature throughout the day. This avoids the need to provide heating and increase the flow rate during the heating mode, resulting in a reduction in gas and electricity consumption.

Figure 14. Operation of one air handling unit that serves a perimeter zone on typical floor during the first day of the mild season.

6. Conclusions

In this paper, we conduct a comprehensive evaluation of how optimization formulation, constraint relaxation, and parameter tuning affect MPC control performance with a high-fidelity building model. Based on the simulation results from different case studies, the following conclusions can be drawn:

• Including domain knowledge in the optimization reformulation might provide additional benefits. The testing results related to the fan coefficients and simultaneous heating and cooling reveal that the optimization algorithm might not be smart enough to avoid unreasonable optimization results.

• When a penalty term is used to handle the thermal comfort constraints, it is critical to consider each zone separately and make the slack variable dynamic. This might be because the accuracy of the zone temperature prediction varies among different zones and/or at different times of the optimization period.

• The electricity saving ratio, the gas saving ratio, and the discomfort reduction are sensitive to the control interval. The computational time can be significantly affected by both the prediction horizon and the control interval.

• The impact from the prediction horizon on the energy savings and thermal comfort varies by the seasons and can be non-monotonic.

• The MPC achieves better energy performance mainly via free-cooling and the "optimal start" of fans. This optimal start can reduce the energy cost of the VAV system during heating mode. It can also improve thermal comfort with a more smooth transition from unoccupied to occupied.

It is worth mentioning that, to ensure generality, this assessment doesn't consider specific needs from uncommon operation conditions. For example, in a very cold climate, the MPC design may need to include additional constraints to avoid a frozen cooling coil for the studied VAV system. In the future, we recommend replicating the assessment that is conducted in this paper, with individual MPC designs, to provide case-specific insights into real-world implementation.

This evaluation clearly emphasizes the need for further study on optimization formulation, constraint relaxation, and parameter tuning of the MPC. It will be beneficial to study how to systematically incorporate the domain knowledge into commonly used optimization formulations to better accumulate the needs for operating building systems. It will be also interesting to develop methods for automatizing the process for parameter tuning in a future study to support large-scale applications of MPC in building systems.

Acknowledgments

This work was supported by the US DOE Office of Energy Efficiency and Renewable Energy, Building Technologies Office. This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE).The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://www.energy.gov/doe-public-access-plan).

Data Availability Statement

Derived data supporting the findings of this study are available from the corresponding author on request.

Disclosure statement

No potential conflict of interest was reported by the author(s).