Efficient management of HVAC systems through coordinated operation of parallel chiller units: An economic predictive control approach

Energy & Buildings 304 (2024) 113879 Available online 3 January 2024 0378-7788/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Energy & Buildings journal homepage: www.elsevier.com/locate/enbuild Efficient management of HVAC systems through coordinated operation of parallel chiller units: An economic predictive control approach J.A. Borja-Conde∗, J.M. Nadales, J.G. Ordonez, F. Fele, D. Limon Department of Systems Engineering and Automation, University of Seville, Seville, Spain A R T I C L E I N F O A B S T R A C T Keywords: Heating ventilation and air-conditioning (HVAC) Chiller sequencing Energy efficiency Building management Economic predictive control Embedded systems Gradient methods In developed countries, air conditioning systems have become major contributors to energy consumption in buildings. Cooling installations made up of independent chiller units connected in parallel pose a challenge in finding the most energy-efficient management strategy. This article proposes a novel economic model predictive control strategy to optimize the operation of multiple-chiller HVAC systems according to an economic cost comprising energy consumption as well as a thermal comfort index. Provided that the gradient of the economic cost function can be calculated or estimated, the proposed controller only entails the solution of a single quadratic programming (QP) problem at each sampling period, reducing computational requirements and thus facilitating the deployment on commercial embedded HVAC management platforms. The proposed controller improves economic performance. As our numerical analysis shows, optimal sequencing is achieved for the case of dualchiller plants. Moreover, the controller adapts to possible variations in the economic criteria, enabling the system to respond to changes in electricity price and user preferences. A realistic case study, using a high fidelity model of a building simulated in TRNSYS, demonstrates the effectiveness of the proposed methodology. In particular, the proposed controller demonstrates to be more efficient than state-of-the-art QP-based economic predictive controllers in achieving a reduction of the energy consumption up to 5.19%, which is in line with the targeted reduction of 7.9% by 2030 in the EU Green Deal’s ‘Fit for 55’ package. 1. Introduction Recent research studies reveal that residential and commercial buildings alone are responsible for up to 40% of the total energy consumption in developed countries, a phenomenon accentuated in recent years due to the increasing demand for building services and comfort levels [1,2]. Furthermore, it is estimated that almost 50% of the amount of energy consumed in buildings is due to heating, ventilation, and air conditioning (HVAC). This has driven significant efforts in recent years to improve the efficiency of this type of systems without compromising comfort and indoor air quality [3]. The importance of this issue is reflected in international standards. For example, the European Union has included the need to increase energy efficiency in buildings in new plans for the EU’s climate transition [4], where heating and cooling are presented as one of the key factors [5]. In the particular case of cooling and heating systems, one of the open challenges—in addition to an adequate architecture design of the installation—is the development of control strategies that optimize the cost of operating multiple-chiller units that supply cooling power to * Corresponding author. E-mail addresses: [email protected] (J.A. Borja-Conde), [email protected] (J.M. Nadales), [email protected] (J.G. Ordonez), ff[email protected] (F. Fele), [email protected] (D. Limon). the installation. This leads to the widely known chiller sequencing control problem (CSCP), which aims at determining the optimal switching thresholds for the different units while guaranteeing that sufficient cooling capacity is provided to the system. This problem, presented as a key challenge in [6], has awakened great interest among the research community during the last few years. Many are the examples of different control approaches that have been proposed to solve the CSCP [7]. For instance, in [8]a non-cooperative game-based formulation is proposed in order to simplify the typically complex centralized CSCP. In [9]a stochastic control algorithm is proposed to deal with the uncertainties associated with the measurement of the cooling load and the estimation of the maximum cooling capacity of the chiller. In [10], a model-free optimization strategy is proposed to maximize the energy efficiency of chilled-water plants with parallel chillers. The strategy utilizes extremum seeking control with penalty terms and adjusts control inputs based on total power consumption. In [11], a collaborative neurodynamic optimization approach to distributed chiller loading is presented, which considers non-convex power consumption functions. https://doi.org/10.1016/j.enbuild.2023.113879 Received 10 August 2023; Received in revised form 6 December 2023; Accepted 28 December 2023 Energy & Buildings 304 (2024) 113879 2 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. Nomenclature: variables and coefficients Primary system  𝑇𝑤𝑟 Temperature of return water to all chiller units....... ◦C 𝑇𝑤𝑠,𝑖 Chiller 𝑖output water temperature................... ◦C  𝑇𝑤𝑠 Temperature of supply water to the secondary system ◦C 𝑇𝑤𝑠 Vector of chiller units’ output water temperatures.... ◦C 𝑞𝑖Water flow of chiller unit 𝑖........................ kg⋅s-1 𝑞𝑇Total recirculated water flow...................... kg⋅s-1 𝜓𝑖Partial load ratio of chiller unit 𝑖...................... % 𝜓Vector of partial load ratios of all chiller units......... % 𝐶𝐸,𝑖 Electrical consumption of chiller unit 𝑖............... kW 𝜅𝑖Maximum cooling capacity available of chiller unit 𝑖.kW 𝑄𝐵,𝑖 Cooling supply by chiller unit 𝑖...................... kW 𝑄𝐵Total cooling supply................................. kW 𝜂𝑖Energy efficiency ratio of chiller unit 𝑖 𝜂 Mean energy efficiency ratio of operating chiller units Secondary system 𝑇𝑗Air temperature of room 𝑗............................ ◦C 𝑇Vector of air temperatures of all rooms............... ◦C 𝑇𝑟𝑒𝑓 𝑗Reference air temperature of room 𝑗.................. ◦C 𝑞𝑓,𝑗 Water flow of fan coil 𝑗........................... kg⋅s-1 𝑇𝑓,𝑗 Fan coil 𝑗output water temperature.................. ◦C 𝑄𝑗Heat absorbed from room 𝑗.......................... kW 𝑣𝑗Opening of the 3-way valve of fan coil 𝑗............... % 𝑣Vector of valve openings of all fan coils............... % 𝐶𝑇,𝑗 Thermal comfort cost of room 𝑗..................... kW External 𝑇𝑜𝑎 Outdoor air dry bulb temperature.................... ◦C 𝑅𝑡ℎ Total horizontal radiation........................ 𝑊⋅m-2 Set of external disturbance parameters Parameters 𝐶𝐸Nominal electrical consumption of chiller unit....... kW 𝜅Nominal cooling capacity of chiller unit............. kW 𝑐0,…,𝑐5Coefficients of function, 𝑓1(𝑇𝑤𝑠,𝑖, 𝑇𝑜𝑎), of the cooling capacity 𝑝0,…,𝑝5Coefficients of functions, 𝑓2(𝑇𝑤𝑠,𝑖, 𝑇𝑜𝑎), of the electrical consumption 𝑘1,…,𝑘3Coefficients of functions, 𝑓3(𝜓), of the electrical consumption 𝑐𝑝Specific heat capacity of water.............. kJ⋅(◦C⋅kg)-1 𝜌Density of water.................................. kg⋅m-3 Economic weighting parameters 𝜇Electricity per-unit cost.......................... $⋅kW-1 𝜆𝑗Thermal discomfort cost of room 𝑗................ $⋅◦C-2 Furthermore, within the range of different control approaches to solve CSCP, it is worth highlighting the use of predictive control strategies due to their ability to deal with constraints while optimizing some user-defined performance index. The reader is referred to [15] and [16], where extensive and detailed reviews of the state of the art on model predictive control (MPC) in HVAC systems are conducted, concluding that MPC outperforms other approaches such as PID and different heuristic control strategies in terms of energy consumption and thermal comfort. This positions MPC as a promising solution for further research. Many examples can be cited that address CSCP: in [17], an efficient heuristic algorithm is proposed to solve the MPC problem based on dynamic programming and mixed integer linear programming (MILP) for obtaining the optimal management strategy of chiller units; in [18], an MPC strategy is employed to anticipate future cooling demand and schedule the operation of a multiple-chiller system where performance maps are obtained through a novel data-driven approach using neural networks; [19] proposes a hybrid predictive control to reduce the switching frequency based on an auto-regressive (ARX) cooling load prediction model embedded into the chiller sequencing control logic; in [20], a comparison between energetic and economical optimization criteria is made on a cooling plant based on air-cooled water chillers and a cold storage system; in [6], an economic MPC based on a MILP formulation is presented for large-scale HVAC systems with multiple chiller units; in [21], the optimal decision sequences for chillers are generated using a mixed-integer nonlinear programming problem. Distributed MPC (DMPC) has also received wide consideration for HVAC control [16]. For instance, [22] proposes a two-layer scheme, where global information is leveraged to optimize overall energy efficiency at the upper layer, while the bottom layer only uses local measures for tracking. In [23]a DMPC approach is considered to optimize a tradeoff between energy cost and thermal comfort of the occupants; however, coupling effects and thermal disturbances are neglected. This aspect is instead addressed by [24,25]and in [26]: in particular, the latter proposes a convex relaxation of the thermal couplings of adjacent zones to allow for the distributed computation of the optimal control sequence. Energy consumption is of primary importance when optimizing the operation of the HVAC system: still, thermal comfort should also be taken into account for a fair evaluation of the performance. While arguably subjective, several standard formulations exist which relate thermal comfort to different factors such as air temperature, mean radiant temperature, relative humidity, clothing insulation and even metabolic rate [27,28]. In [29], a Pareto multi-objective approach to minimize energy cost and thermal discomfort is applied considering different discomfort indicators to compare Fanger’s model [30]with the ASHRAE 55 adaptive model. However, the evaluation of these comfort models is complex and computationally demanding [31], prompting the use of simplified approaches instead. A common adaptation is to constrain the room temperature within a suitable range, as done in [32–34]. In practice, however, it is desirable to dynamically update the target room temperature online [35]. In [36], a temperature range constraint is also used, but this can be adjusted according to personal preferences. Alternative approaches in the literature rely on a simplified comfort model. In [37], the authors present a multi-objective MPC for commercial building HVAC. They show that the relationship between temperature and comfort (based on ASHRAE formula) is well approximated by a quadratic function, whereby an equivalent comfort index can derived as the square error between the temperature of each room and some target temperature setpoint; such an approach is also adopted in [38]. Other works proposed the use of linearized comfort models [39,40]. In this context, it is also worth of mention the use of personalized climatization devices [41,42]. These provide a means to finely adjust the comfort of each individual in a given environment, by directly influencing the thermal conditions in a circumscribed space around them. Their use, in combination with centralized HVAC systems, can allow significant increases in the building energy efficiency. In [43], control algorithms and machine learning heuristics are used to improve energy efficiency and individual comfort compared to traditional HVAC systems. The work of [44] presents a comfort system combining personal thermal climatization technologies to regulate individual comfort. Another interesting avenue of research is on the use of variable water flow HVAC systems, with the objective of reducing the electrical consumption of the pumps, and operating the chiller units more efficiently. For Energy & Buildings 304 (2024) 113879 3 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. Fig. 1. Comparison between a standard two-layer RTO+MPC scheme and the proposed single-layer gradient-based scheme. The proposed approach avoids the typical drawbacks of the RTO layer [12]: Replacing the NLP-RTO problem by its first-order Taylor approximation—which can be estimated in real time by means of online measurements [13]—in the tracking MPC cost function reduces the control problem to a quadratic programming (QP) formulation, thus facilitating its implementation [14]. Please refer to Sections 3and 4for a detailed description of each approach and its notation. instance, in [45] optimal control strategies for variable speed pumps is presented, leading to energy savings of 12-32% in the chilled-water pumps’ energy demand; in [46]the optimization of the pump operation is shown to afford up to 51.11% reduction in the pump’s demand. In [47], an MPC approach is presented where the water flow rates are decision variables. The energy-saving potential of variable water flow approaches can be increased by combining them with split-range control in the fancoils, i.e., using both water flow and fan speed as optimization variables [48,49]. While these approaches are particularly attractive in new installations, the benefits may not be as pronounced in retrofit projects [50]. In this work, we focus on economic MPC to optimally determine the management of a system with multiple chiller units in a constantflow cooling installation to meet the cooling demand while optimizing economic performance criteria compromising energy costs and user comfort. It is a common practice to address the CSCP by defining fixed cooling load thresholds to determine the number of chillers in operation [51–53]. However, the optimal choice for the latter actually depends on the operating conditions, therefore such an approach is inherently inefficient. Thus, in this work we take a different approach, and seek dynamic chiller sequencing based on an ad-hoc definition of the economic MPC cost function. Economic optimization of a dynamical process is typically approached in the literature using a two-layer control structure [54]. In the first high-level layer, a real-time optimizer (RTO) determines the optimal operating point by minimizing an economic cost function using a high-fidelity model of the system. This typically involves solving an online nonlinear programming (NLP) problem. In the second layer, these optimal setpoints are used as references by a tracking MPC [55], which drives the system towards the corresponding equilibrium point. As it has been demonstrated in [56], real-time optimization of multiple chiller cooling systems outperforms weekly and daily control strategies. Furthermore, due to the varying time scales at which each layer operates and the utilization of diverse models, there is a possibility that the desired setpoint provided by the RTO may prove unattainable or exhibit sluggish responsiveness to potential economic changes. In addition, the need for an accurate thermal building model adds complexity to this solution [57], and the implementation of an NLP problem restricts its applicability on target hardware platforms. 1.1. Contribution In this work we propose a novel QP-based economic predictive controller for optimal management of multiple-chiller HVAC systems according to an economic cost that takes into account energy consumption as well as a user comfort index. Inspired by [12], the proposed controller integrates economic real-time optimization into the tracking predictive controller by means of the gradient of the energy consumption curve which, within the scope of the manuscript, is assumed to be known (see the comparison diagram depicted in Fig. 1). The resulting control strategy features the following properties: (i) The controller only requires the solution of a single QP problem at each sampling interval. This enables its implementation on commercial industrial platforms for HVAC systems [14]. (ii) The proposed controller reduces the economic cost associated with the operation of the plant; furthermore, we show through numeri- Energy & Buildings 304 (2024) 113879 4 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. Fig. 2. Scheme of the cooling system with 𝑀-chiller units of constant-flow and the secondary system with 𝑁rooms. cal analysis that optimal sequencing is achieved in the case of two parallel chillers without requiring a mixed-integer approach. (iii) The proposed control strategy automatically adapts to possible changes in the economic cost function derived, for instance, from variations in the energy price or in the users’ preferences. To the best of the author’s knowledge, there are no predictive control approaches in the literature featuring aforementioned properties for multiple-chiller HVAC plants. The only comparable approach we have found is [58], which employs a QP-based economic MPC to operate a single chiller plant considering a linear economic cost function. In the case study presented in Section 6, these properties are demonstrated and our approach is compared to that proposed in [58], suggesting a significant improvement in energy efficiency. In this study, the central chilled water system is assumed to consist of a primary-only configuration with constant water flow rate and two chillers operating in parallel. 3-way valves are used for controlling the chilled water flow through the fan coils. If the system to be operated included more chillers or additional heating sources, such as boilers, the proposed MPC-based methodology could be also used by extending appropriately the prediction model and its manipulable inputs. The rest of the paper is structured as follows. In Section 2the cooling production installation to be controlled is presented and the economic characterization of the chiller units is detailed. In Section 3 the optimal management of the chiller units is justified by analyzing consumption curves with respect to cooling production. In Section 4 the proposed single-layer gradient-based control strategy is described. Section 5introduces the case study and provides a detailed description of the building model implemented in TRNSYS. In Section 6the main results are analyzed and benchmarked against a baseline strategy. Finally, Section 7summarizes the main conclusions of this work. 2. Centralized cooling system description A general scheme of the system under study is shown in Fig. 2. A primary system composed of 𝑀chiller units CH1, CH2,..., CH𝑀refrigerates water in a constant-flow cooling installation with the purpose of regulating the temperature of the secondary system—a building with a total of 𝑁rooms. The temperature of each room is denoted as 𝑇𝑗, 𝑗=1, … , 𝑁. The constant total water flow 𝑞𝑇has temperature  𝑇𝑤𝑟 before being split and distributed across the chiller units, where the individual flows are denoted by 𝑞𝑖, 𝑖 =1, … , 𝑀. Each of these are characterized by a temperature 𝑇𝑤𝑠,𝑖 at the outlet of the chiller. All output flows are then mixed again to produce the same constant flow 𝑞𝑇at a new temperature  𝑇𝑤𝑠, which depends on  𝑇𝑤𝑟 and the cooling effect produced by the chiller unit. The output flow of the chiller units is distributed among the different rooms—the secondary system shown in Fig. 2, where 𝑞𝑓,𝑗 denotes the flow corresponding to room 𝑗. In each room, heat transfer 𝑄𝑗occurs through a fan coil unit 𝐹𝑗(positive 𝑄𝑗denotes heat absorbed from the room). We consider all fan coils operating at constant fan speed. The value of 𝑄𝑗is modulated by the aperture 𝑣𝑗of a three-way valve which determines the portion of 𝑞𝑓,𝑗 that circulates through coil 𝑗. The latter flow exits the coil at temperature 𝑇𝑓,𝑗. Finally, the output flows of all rooms are mixed together, resulting in the recirculating flow to the chiller units. The objective of the proposed control approach is to control the temperature of each room 𝑇𝑗to maximize thermal comfort while optimizing energy consumption. This is accomplished by manipulating the cooling effect produced by each chiller unit and adjusting the aperture of the three-way valves. 2.1. Model of the chiller unit Modeling cooling machines is generally a challenging task due to their internal complexity. Furthermore, in the case of commercial chiller units, their core architecture and configuration are confidential business designs, and manufacturers provide only limited information on the system performance. For a chiller unit 𝑖, the available data to describe its model typically includes yield curves based on performance parameters. These parameters usually include the chiller capacity (CAP) 𝜅𝑖, which represents the maximum cooling power available to be applied to the water flow; the electrical power consumed by the machine 𝐶𝐸,𝑖; and the energy efficiency ratio (EER) 𝜂𝑖, i.e., the ratio between the actual cooling production 𝑄𝐵,𝑖 and the power consumption, which is commonly used as the standard measure of the system’s performance [59]. Therefore, in this work the chiller units are modeled using the curves that describe the behavior of these parameters [60]. These parameters are expressed through standard curves, which can be identified using operation data provided by manufacturer or obtained through experimental measures [61,62]. The CAP is described as a function of the external outdoor air temperature 𝑇𝑜𝑎 and the chiller supply water flow temperature 𝑇𝑤𝑠,𝑖. Using these variables, the CAP of each chiller unit 𝜅𝑖is calculated as 𝜅𝑖=𝜅𝑖⋅𝑓1,𝑖(𝑇𝑜𝑎,𝑇𝑤𝑠,𝑖),(1) where 𝜅𝑖is the nominal capacity of the chiller unit 𝑖and is a constant parameter provided by the system manufacturer. The function 𝑓1,𝑖 is defined as 𝑓1,𝑖(𝑇𝑜𝑎,𝑇𝑤𝑠,𝑖)=𝑐0,𝑖 +𝑐1,𝑖𝑇𝑤𝑠,𝑖 +𝑐2,𝑖𝑇2 𝑤𝑠,𝑖 +𝑐3,𝑖𝑇𝑜𝑎 +𝑐4,𝑖𝑇2 𝑜𝑎 +𝑐5,𝑖𝑇𝑜𝑎𝑇𝑤𝑠,𝑖 , (2) Energy & Buildings 304 (2024) 113879 5 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. where 𝑐0,𝑖, … , 𝑐5,𝑖 are the estimated coefficients. The partial load ratio (PLR) 𝜓𝑖∈[0, 1] of a chiller unit 𝑖is defined as the actual cooling effect 𝑄𝐵,𝑖 produced by the chiller divided by the maximum cooling effect available, and can be expressed as 𝜓𝑖=𝑄𝐵,𝑖 𝜅𝑖 =𝑐𝑝⋅𝜌⋅𝑞𝑖⋅( 𝑇𝑤𝑟 −𝑇𝑤𝑠,𝑖) 𝜅𝑖 ,(3) where 𝑐𝑝and 𝜌are the specific heat and density of water, and 𝑞𝑖is the input flow. Knowing the value of 𝑇𝑤𝑠,𝑖 and 𝜓𝑖at which the chiller unit is operating, the electrical consumption 𝐶𝐸,𝑖 can be calculated as 𝐶𝐸,𝑖 =𝐶𝐸,𝑖 ⋅𝑓2,𝑖(𝑇𝑜𝑎,𝑇𝑤𝑠,𝑖)⋅𝑓3,𝑖(𝜓𝑖).(4) This expression depends on the nominal consumption of the chiller 𝐶𝐸,𝑖, which is provided by the manufacturer, and on two functions whose expressions are 𝑓2,𝑖(𝑇𝑜𝑎,𝑇𝑤𝑠,𝑖)= 𝑝0,𝑖 +𝑝1,𝑖𝑇𝑤𝑠,𝑖 +𝑝2,𝑖𝑇2 𝑤𝑠,𝑖 +𝑝3,𝑖𝑇𝑜𝑎 +𝑝4,𝑖𝑇2 𝑜𝑎 +𝑝5,𝑖𝑇𝑜𝑎𝑇𝑤𝑠,𝑖 ,(5) 𝑓3,𝑖(𝜓)=𝑘1,𝑖𝜓𝑖−𝑘2,𝑖𝜓2 𝑖+𝑘3,𝑖𝜓3 𝑖,(6) where 𝑝0,𝑖, … , 𝑝5,𝑖 and 𝑘1,𝑖, … , 𝑘3,𝑖 are the estimated coefficients. Finally, the EER is calculated as the effective cooling power divided by the machine electrical consumption, i.e., 𝜂𝑖=𝑄𝐵,𝑖 𝐶𝐸,𝑖 .(7) In typical commercial chillers, the operator is allowed to set the outlet flow temperature, which is used as reference by some embedded control loop. Here, we assume that the units are equipped with low level controllers of the water temperature 𝑇𝑤𝑠,𝑖 that result in overdamped dynamics of the closed loop system. On the other hand, the steady-state characteristic curve of the chiller is a nonlinear function of 𝑇𝑤𝑠,𝑖,  𝑇𝑤𝑟 and 𝑇𝑜𝑎. Taking this into account, we propose a Hammerstein model [63]to describe the dynamics of the chiller as a first order system with a time constant 𝜏𝑖and unitary DC-gain (expressing the over-damped behavior) in series with the steady-state characteristic function, considering  𝑇𝑤𝑟 and 𝑇𝑜𝑎 as exogenous inputs. Remark 1. Motivated by the fact that energy performance of chillers is usually characterized with respect to the PLR and to simplify the decision making problem, we directly consider the PLR of each unit as a decision variable for the proposed controller; this is a common practice in the literature on HVAC systems control (see, e.g., [64]and references therein). Given the calculated PLR and the measures of the inlet and outlet temperature of each unit, the setpoint of the outlet temperature can be obtained from (1)-(3). 2.2. Multiple-chiller cooling system For a cooling system composed of 𝑀chiller units connected in parallel, from the mass and energy balance, the mixed supply flow temperature  𝑇𝑤𝑠 can be calculated as  𝑇𝑤𝑠 = 𝑀 ∑ 𝑖=1 (𝑇𝑤𝑠,𝑖 ⋅ 𝑞𝑖 𝑞𝑇),(8) and the total cooling power supply 𝑄𝐵can be calculated as 𝑄𝐵=𝑐𝑝⋅𝜌⋅𝑞𝑇⋅( 𝑇𝑤𝑟 − 𝑇𝑤𝑠).(9) The total electrical consumption 𝐶𝐸is calculated as the sum of the electrical consumption of each chiller unit as 𝐶𝐸= 𝑀 ∑ 𝑖=1 𝐶𝐸,𝑖 .(10) Fig. 3. Input/output model of the secondary system in Fig. 2. Then, the average EER of the whole system 𝜂 is calculated as 𝜂 =𝑄𝐵 𝐶𝐸 .(11) 2.3. Model of the secondary system The total cooling power supply, 𝑄𝐵, can also be defined as the total cooling demand of the building, since we assume that the energy losses in the water distribution system are negligible. Therefore, 𝑄𝐵can also be calculated as the sum of the cooling power given to each room 𝑄𝑗as 𝑄𝐵= 𝑁 ∑ 𝑗=1 𝑄𝑗,(12) where 𝑄𝑗can be calculated as 𝑄𝑗=𝑐𝑝⋅𝜌⋅𝑣𝑗⋅𝑞𝑓,𝑗 ⋅(𝑇𝑓,𝑗 − 𝑇𝑤𝑠).(13) The flows out of the individual fan coils 𝑞𝑓,𝑗 are mixed again obtaining a total flow 𝑞𝑇at temperature  𝑇𝑤𝑟 which is eventually recirculated into the chiller units. Measures of the output temperature of the air temperature of the room 𝑇𝑗are assumed to be available. See Fig. 2. The overall secondary system is then a nonlinear dynamic system with  𝑇𝑤𝑠 and the vector 𝑣as inputs; and the temperatures of the rooms 𝑇and the return temperature to the chiller units  𝑇𝑤𝑟 as outputs. The set of exogenous disturbances comprises occupancy, solar radiation and outdoor temperature, among others. This is depicted in Fig. 3, where Σrepresents the nonlinear dynamical system, relating the inputs and outputs: (𝑇, 𝑇𝑤𝑟)=Σ(𝑣,  𝑇𝑤𝑠,).(14) It is worth emphasizing that it is generally hard to characterize Σin (14), due to the system’s complexity and nonlinear dynamics. Hence, in this work we use a linear time invariant (LTI) model of the secondary system obtained from identification. 3. Economic optimization In this work we aim at minimizing the total demand of energy from a series of chiller units connected in parallel while complying with certain thermal comfort criteria. Many different research studies can be found in the literature addressing thermal comfort models in indoor environments. The reader is referred to [65]where a detailed review of different personal comfort models is provided. On the other hand, different standards specify comfort levels in indoor environments (see [27] or [28]for the case of building located in Spain). Based on this, a function 𝐶𝑇,𝑗, which accounts for thermal comfort of room 𝑗, can be defined. Let 𝑓𝑒𝑐𝑜 ∶ℝ1+2𝑀+2𝑁→ℝbe the economic cost function to be optimized given by the following expression: 𝑓𝑒𝑐𝑜(𝑇𝑜𝑎,𝑇𝑤𝑠,𝜓,𝑇)=𝑓𝑒𝑐𝑜,𝐶𝐸(𝑇𝑜𝑎,𝑇𝑤𝑠,𝜓)+𝑓𝑒𝑐𝑜,𝐶𝑇(𝑇),(15) with Energy & Buildings 304 (2024) 113879 6 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. Fig. 4. EER of the commercial chiller unit used in the case study (Daikin EWAD740, 𝜅1=𝜅2=11kW and constant flow rate 𝑞1=𝑞2=0.5m3⋅h-1), as a function of 𝑄𝐵and  𝑇𝑤𝑠, when 𝑇𝑜𝑎 =27.5◦C. The figure evidences a monotonic relation between  𝑇𝑤𝑠 and 𝜂, indicating that higher output water temperature values correspond to greater efficiency. Additionally, the data reveals that operating with a single unit is more efficient when the cooling supply is lower, whereas using two units becomes more advantageous as the cooling supply increases. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) 𝑓𝑒𝑐𝑜,𝐶𝐸(𝑇𝑜𝑎,𝑇𝑤𝑠,𝜓)=𝜇⋅ 𝑀 ∑ 𝑖=1 𝐶𝐸,𝑖(𝑇𝑜𝑎,𝑇𝑤𝑠,𝑖,𝜓𝑖),(16) and 𝑓𝑒𝑐𝑜,𝐶𝑇(𝑇)=⋅ 𝑁 ∑ 𝑗=1 𝜆𝑗⋅𝐶𝑇,𝑗(𝑇𝑗),(17) where 𝐶𝐸,𝑖 is the electrical consumption of each chiller unit according to (4), and 𝐶𝑇,𝑗 is the thermal discomfort cost; 𝜆𝑗and 𝜇are the economic weighting parameters. Note that a different value of 𝜆𝑗 could be applied in order to prioritize certain rooms if required. In this work, we assume that the same priority is given to all rooms, i.e., 𝜆𝑗=𝜆, 𝑗=1, … , 𝑁. Note that the thermal discomfort, as motivated in the literature review, can be accounted by using the functions of the standards focused on comfort evaluation [27,28], or by using approximated functions around desired temperature setpoints—such as quadratic [37,38]or linearized forms [39,40]. In our case study, described in Section 5, the approximated quadratic form is employed. Then, an economic optimization problem, where the decision variables are 𝜓={𝜓1, … , 𝜓𝑀}and 𝑣 ={𝑣1, … , 𝑣𝑁}, is expressed as 𝜓†∗,𝑣†∗∈arg min 𝜓,𝑣 𝑓𝑒𝑐𝑜(𝑇𝑜𝑎,𝑇𝑤𝑠,𝜓,𝑇) s.t. (1), (3), (8), (14), (18) where 𝜓†∗and 𝑣†∗minimize the value of the cost function. Given these considerations, the performance curves of the chiller units can be analyzed to determine the expected optimal behavior of the controlled system. Fig. 4illustrates the EER of the system when operating with either one or two units at the same point of operation. It can be noticed that as  𝑇𝑤𝑠 increases, 𝜂 also improves. Therefore, the controller should strive to maintain a higher  𝑇𝑤𝑠. Additionally, it is observed that there is a range in the lower cooling supply where operating with a single chiller is more efficient, while employing two chiller units becomes advantageous in the higher cooling supply. It is important to note that this practical phenomenon—extensively examined in Fig. 5 focusing on electrical consumption instead of EER—arises when all co- Fig. 5. Electrical consumption of the scaled system used in the case study (Daikin EWAD740, 𝜅1=𝜅2=11kW). Blue shaded area corresponds to chiller unit CH1being prevalent over CH2(i.e., 𝜓1>𝜓 2), up to CH2off (𝜓2=0 in bright blue). Red shaded area indicates instead where the cooling load is shared (almost) equally over the chiller units (i.e., 𝜓1≈𝜓2)—bright red edge represents 𝜓1=𝜓2. It can be appreciated that there is a range in the lower heat demand where the consumption of the cooling system is minimized operating only one chiller. However, as the demand increases, the operation with paired chiller units is more efficient. We emphasize that there is no trivial solution for this problem, since the operation point at which this transition occurs depends on the operation conditions (and is automatically found by the controller); note that in the figure the supply water temperature and flow rate, as well as the external temperature, are fixed to  𝑇𝑤𝑠 =7◦C, 𝑇𝑜𝑎 =27.5◦Cand 𝑞𝑇=1m3⋅h-1for illustrative purposes. efficients in (6)are positive, which is commonly found in commercial chiller units. This underscores the importance of independent operation for each chiller unit in a multiple-chiller system. Our proposed approach effectively addresses these issues and offers a solution to optimize the system’s performance, as presented in Section 6. 4. Gradient-based single-layer economic MPC Typically, the primary goal of achieving economic optimization in a control system is accomplished through the use of an upper level RTO, which computes (18). This is responsible for determining the optimal operation point, which is then passed down to a real-time control layer (e.g., an MPC controller). The latter layer calculates the control actions required for the plant to attain the optimum, considering a simplified dynamic model of the plant and constraints on the state and control variables. However, as motivated in Section 1, this approach leads to different issues. Therefore, we adopt the single-layer control approach proposed in [12](see the comparison in Fig. 1). By adhering to this strategy, the functionality of the RTO is achieved through the incorporation of the first-order Taylor approximation of 𝑓𝑒𝑐𝑜 into the cost function of the tracking MPC. The MPC requires a model of the plant, which we identify as a LTI model. The manipulable inputs of this model are the PLR of each chiller unit 𝜓𝑖, and the opening of each room’s fan coil valve 𝑣𝑗. The outputs are the temperatures of each of the 𝑁rooms in the building, 𝑇𝑗, and the temperature of the mixed supply water flow of the chiller units,  𝑇𝑤𝑠. The linear model employed by the MPC is thus in the form 𝑥𝑘+1 =𝐴𝑥𝑘+𝐵𝑢𝑘, 𝑦𝑘=𝐶𝑥𝑘+𝑑𝑘,(19) where 𝑘is the discrete-time index, 𝑥𝑘∈ℝ𝑛𝑥is the state of the system, 𝑑𝑘∈ℝ𝑛𝑦is the external disturbance (included to achieve an Energy & Buildings 304 (2024) 113879 7 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. offset free control), 𝑦𝑘∈ℝ𝑛𝑦is the measured output of the system (𝑇𝑗, 𝑗=1, … , 𝑁and  𝑇𝑤𝑠), 𝑢𝑘∈ℝ𝑛𝑢is the manipulable input of the system (𝑣𝑗, 𝑗=1, … , 𝑁and 𝜓𝑖, 𝑖 =1, … , 𝑀), 𝑛𝑥, 𝑛𝑦=𝑁+1 and 𝑛𝑢=𝑁+𝑀are the number of states, measured outputs and manipulable inputs, respectively, and 𝐴, 𝐵and 𝐶are the matrices that describe the model obtained through identification. To this end, the system is excited by each of the different control actions—𝑣𝑗, 𝑗=1, … , 𝑁and 𝜓𝑖, 𝑖 =1, … , 𝑀, and a linear regression procedure is applied to obtain linear models of the stationary behavior of the system. Some preprocessing might be required in order to remove negligible states [66]. Note that a compromise between accuracy and overfitting must be achieved by adequately selecting the dimensions of the identified model. If the state of the system is not accessible, an observer is designed to estimate the state of the model together with the output disturbance. To this end, a linear observer with gains 𝐿𝑥for the state and 𝐿𝑑for the disturbance can be employed.  𝑑𝑘= 𝑑𝑘−1 +𝐿𝑑(𝑦𝑘−𝐶𝑥𝑘− 𝑑𝑘−1)(20a) 𝑥𝑘+1 =𝐴𝑥+𝐵𝑢𝑘+𝐿𝑥(𝑦𝑘−𝐶𝑥𝑘− 𝑑𝑘)(20b) where 𝑥𝑘and  𝑑𝑘denote the estimated state and disturbance at time 𝑘, respectively. Then, the proposed single-layer gradient-based MPC is formulated through the following optimization problem: min 𝒙,𝒖,𝑥𝑠𝑠,𝑢𝑠𝑠 𝑁𝑝−1 ∑ 𝑘=0 (|𝑥𝑘−𝑥𝑠𝑠|2 𝑄+|𝑢𝑘−𝑢𝑠𝑠|2 𝑅) +|𝑥𝑁−𝑥𝑠𝑠|2 𝑃 +∇ 𝑦𝑓𝑒𝑐𝑜(𝑢𝑘−1,𝑦𝑘−1)⋅(𝑦𝑠𝑠 −𝑦𝑘−1) +∇ 𝑢𝑓𝑒𝑐𝑜(𝑢𝑘−1,𝑦𝑘−1)⋅(𝑢𝑠𝑠 −𝑢𝑘−1) +|𝑢𝑠𝑠 −𝑢𝑘−1|2 𝑆(21a) s.t. 𝑥0=𝑥0(21b) 𝑥𝑘+1 =𝐴𝑥𝑘+𝐵𝑢 𝑘(21c) 𝑦𝑘=𝐶𝑥 𝑘+ 𝑑0(21d) 𝑥𝑠𝑠 =𝐴𝑥𝑠𝑠 +𝐵𝑢 𝑠𝑠 (21e) 𝑦𝑠𝑠 =𝐶𝑥𝑠𝑠 + 𝑑0(21f) 𝑦𝑘,𝑦𝑠𝑠 ∈,𝑢 𝑘,𝑢𝑠𝑠 ∈(21g) for 𝑘=0,…,𝑁 𝑝−1, where 𝑁𝑝is the prediction horizon of the finite optimization problem, 𝑥𝑠𝑠 ∈ℝ𝑛𝑥is the target artificial steady state of the system, 𝑢𝑠𝑠 ∈ℝ𝑛𝑢 and 𝑦𝑠𝑠 ∈ℝ𝑛𝑦are the corresponding artificial steady input and output when 𝑥𝑠𝑠 is reached. Note that this is a reachable equilibrium point of the system. The MPC calculates this point in such a way that the solution is feasible while minimizing the discrepancy between the system trajectory in the prediction horizon and the artificial steady state. The cost function (21a)is composed of the regulation terms weighted by matrices 𝑄 ≻0 ∈ℝ𝑛𝑥×𝑛𝑥and 𝑅 ≻0 ∈ℝ𝑛𝑥×𝑛𝑢, the terminal cost penalty weighted by 𝑃≻0 ∈ℝ𝑛𝑥×𝑛𝑥to ensure stability [67], and the first-order Taylor approximation of the economic cost function 𝑓𝑒𝑐𝑜. (21b)sets the initial state; (21c)and (21d)are the state transition and output model of the system; (21e)and (21f) defines the stationary state of the system; and (21g) defines the constraints on the system variables. Note that the set must be appropriately defined to ensure that the feasibility of (21) is maintained in case of uncertainty. The proposed controller ensures nominal asymptotic stability to the optimal equilibrium point even when the economic cost changes. Moreover, the controller exhibits inherent robustness from model mismatches whenever the optimization problem remains feasible [12]. Table 1 Physico-thermal building materials properties. Parameter Ext. Wall Int. Wall Floor Roof Window Thickness (m) 0.171 0.100 1.106 0.141 0.017 Density (kg⋅m-3)852.61400 121.3148.5626.8 Capac. (J⋅(kg⋅K)-1)1139 1000 852.6848.1188 𝑈 (W⋅(m2⋅K)-1) Air-air 0.512 2.550 0.039 0.318 3.0 Surf-surf 0.556 2.550 0.040 0.334 3.0 5. Case study The example case under study reflects a typical office building in Madrid (Spain), equipped with a HVAC system designed in accordance with the Spanish standards and regulations. This is a centralized waterbased system with two chiller units. There are five air-conditioned rooms equipped with independent fan coil units. In addition, the building model includes a corridor which connects the rooms and contributes to the thermal coupling between them (see Fig. 7). Temperature of this corridor is considered to be uncontrollable. To simulate the building’s thermal dynamics we rely on the specialized software TRNSYS [68,69]. In particular, by using the TESS Component Library, TRNSYS allows to replicate the transient behavior of HVAC systems coupled with building thermal loads [70]. The TRNSYS model scheme is shown in Fig. 6. 5.1. Building layout and dimensions The building is divided into six different zones: five offices and a corridor. See Figs. 7and 8. The dimensions of Zone 1 are 8 ×4 ×3.5m, the dimensions of Zones 2 to 5 are 4 ×4 ×3.5m, and the dimensions of the corridor are 12 ×3.5 ×3.5m (last dimension indicates ceiling height). Windows are 1.75 ×1.2m. The north facade has two windows in Zone 1 and one in Zone 2, the south facade has one window in each zone, the west facade has one window in the corridor. We note that the thermal properties of windows differ from those of walls, making them a primary source of various disturbances, such as solar radiation. With this information, TRNSYS3D—a plugin for SketchUp [71]that allows the user to draw multizone buildings and import the geometry (including building self-shading and internal view factors for radiation exchange) into the TRNSYS Building environment (TRNBuild)” [72]—is used to develop the building model. 5.2. Building materials Construction materials are defined according to the “High Mass Case Study” of ANSI/ASHRAE Standard 140-2004 [73].1The main physicothermal properties are specified in Table 1. Hemispherical infrared emittance of ordinary uncoated glass is 90%. Zone capacitance, because of air, furniture, equipment, etc., is ten times air capacitance inside the zone. These data are incorporated in the TRNSYS building model generated with TRNSYS3D. 5.3. Internal gains and ventilation TRNSYS allows to set internal gains by setting occupancy, lighting and working equipment schedules, which are specified in Table 2for five weekdays (Monday to Friday). This table also specifies the standard used to consider each type of gain. 1To address numerical convergence issues, the floor construction in this case study includes a thick insulating layer to effectively decouple the building’s floor from the ground, as per [73]; we note that this results in a very low heat transfer coefficient. Energy & Buildings 304 (2024) 113879 8 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. Fig. 6. TRNSYS model scheme. Controller is a Type155 block which applies a control law specified in MATLAB. Chiller 1 and Chiller 2 are equation blocks modeling the chiller. FC 1 to 5model 3-way valves and volumetric flow based two-pipe fan coil (based on equation and Type928 blocks). Valve is an equation block which models a water diverter at the inlet of the chillers. Two Collector components (based on equation blocks) simulate water mixers at the outlet of the chillers (left) and fan coils (right). Building is a Type56 block which models the dynamics of multiple thermal zones in a building. Weather is a Type15-6 block which applies weather conditions based on given input weather data. Fig. 7. Building layout. Outside air is delivered directly to the space through the fan coil unit, which includes a mixer to blend it with the room’s recirculated air. The air exchange rate is set in compliance with current Spanish regulations [74]. In particular, for office buildings as the one considered here, a fresh air inflow of 12.5 l⋅(s⋅person)-1is required by law. Therefore, the outside air flow for Zone 1 is set to 50 l⋅s-1, for Zone 2 to 87.5 l⋅s-1, and for Zones 3, 4 and 5 to 12.5 l⋅s-1. 5.4. HVAC system A water pump (with power rating of 0.16 kW) circulates water in the central chilled water system; this power load remains constant while the system is operational. For the remaining elements of the HVAC system, Table 2 Building internal gains schedule. Zone Occ. Occ. Sched. Activity No. PCs3Lighting4 Z1 (4 p.) 27-15:00 Office127-18:00 29-18:00 2 Z2 (7 p.) 713-14:00 Meeting1113-14:00 Z3 (1 p.) 17-15:00 Office117-15:00 Z4 (1 p.) 19-18:00 Office119-14:00 15-18:00 Z5 (1 p.) 19-14:00 15-18:00 Office119-18:00 Corridor 37-7:05 Light work207-18:05 49-9:05 18-18:05 515-15:05 713-13:05 914-14:05 1ISO 7730:2005 [28]: seated, light work, typing. 2ISO 7730:2005: standing, light work or working slowly. 3TRNSYS documentation: 230 W with color monitor. 4TRNSYS documentation: related floor area equals the total room’s, 10 W/m2, 40% fluorescent tube. the following two classes of cooling equipment are integrated in the TRNSYS model: 5.4.1. Chiller units Two identical chiller units, independently operated, constitute the cold production system. The chiller model reflects a scaled version of the Daikin EWAD740, with nominal capacity of 11 kW. These are characterized as described in Section 2.1 and modeled in TRNSYS by using equation blocks. In particular, (1), (2), (3)and (8)are relative to the output water temperature, while (4), (5), (6)and (10) express the energy consumption. The values of the coefficients are specified in Ta- Energy & Buildings 304 (2024) 113879 9 J.A. Borja-Conde, J.M. Nadales, J.G. Ordonez et al. Fig. 8. Building dimensions. Table 3 Chiller unit coefficients. Coef. f1 Value Coef. f2 Value Coef. f3 Value 𝑐0(-) 9.1e-1𝑝0(-) 4.8e-1–– 𝑐1(◦C-1)4.1e-2𝑝1(◦C-1)8.6e-3𝑘1(-) 2.10 𝑐2(◦C-2)5.3e-4𝑝2(◦C-2)5.2e-4𝑘2(-) 1.78 𝑐3(◦C-3)2.9e-3𝑝3(◦C-1)7.3e-3𝑘3(-) 1.46 𝑐4(◦C-4)-1.6e-4𝑝4(◦C-2)1.6e-4–– 𝑐5(◦C-5)-4.9e-4𝑝5(◦C-2)7.7e-5–– Table 4 Fan coil units parameters. Zone Fan coil model (Carrier) Nominal air flow rate Fan input power Nominal sensible capacity Nominal latent capacity (l⋅s-1)(W) (kW) (kW) 142GWE 012 203 45 5.9 4.8 242GWE 010 167 50 4.7 3.9 342GWE 008 86 20 4 3.3 442GWE 004 125 35 2.4 2.3 542GWE 008 86 20 4 3.3 Fig. 9. Outdoor air temperature and total horizontal solar radiation. ble 3; these were obtained by least squares estimation and validated on operational data provided by the manufacturer through private communication. 5.4.2. Fan coils Fan coil units are sized according to building cooling demand of each zone, considering maximum loads of sensible and latent heat, in accordance with the Spanish standards and regulations [74]. The parameters of these units are specified in Table 4. Fan coil units are modeled in TRNSYS by Type928, that models a volumetric flow based two-pipe fan coil [70]. 5.5. Weather conditions The weather conditions have been obtained from the weather database available in TRNSYS for the city of Madrid (Spain). Together with the Type15 component (weather data reading and processing) [76]. The experiments of the case study simulate typical summer weather. The system is subject to realistic time-varying weather conditions, such as solar radiation and outdoor temperature, among others (see [77]). Fig. 9shows the time profile of the main weather variables considered in our simulations. 5.6. Controller The control algorithm detailed in Section 4is implemented in MATLAB using the Yalmip tool [78]and the Gurobi solver [79], while communication with TRNSYS is established through Type155 component [76]. We chose a sampling time of 4 minutes, given the time constant of the system is in the order of a few hours. The prediction horizon is 1hour and 40 minutes, a trade-off between computing time and performance. 𝑄, 𝑅, and 𝑆are unitary matrices. The terminal cost weight matrix, 𝑃, is calculated solving the Ricatti equation. The weighting parameters for thermal comfort, 𝜆, and for electricity cost, 𝜇, are set as detailed in Sections 6.1 and 6.2. In this case study, we assume that desired reference temperatures are externally given as setpoints—which can be calculated, e.g., following given comfort standards [27,28], which are included in the objective function of the MPC as 𝐶𝑇,𝑗 ∶=|𝑇𝑗−𝑇𝑟𝑒𝑓 𝑗|2.(22) Then, thermal comfort is maximized by minimizing the square error between the temperature of each room and the given comfort temperature setpoint. We note that, for the case of Fanger’s model, this is tantamount to considering a convex approximation of Fanger’s discomfort curve in the objective function of our EMPC scheme. This can be seen, for instance, in [37, Fig. 5], which shows a comparison of discomfort indices based on ASHRAE formula and their quadratic approximation; and in [80, Fig. 5], showing that Fanger’s curve is well approximated by a quadratic function for a deviation from the optimum temperature within ±2 ◦C. Note that in this case it follows from (22)that 𝑓𝑒𝑐𝑜,𝐶𝑇is quadratic, therefore this term can be directly included in (21a)instead of its Taylor approximation. This is possible because the result has still a QP formulation. Note also that, since the gradient of 𝑓𝑒𝑐𝑜,𝐶𝐸(𝑢𝑘−1, 𝑇𝑤𝑠)with respect to 𝑇𝑤𝑠 is approximately zero, and the increment of 𝑇𝑤𝑠 at each sampling