An Exact Algorithm for the Optimal Chiller Loading Problem and Its Application to the Optimal Chiller Sequencing Problem

: The optimal management of multiple chiller systems calls for the solution of the so-called optimal chiller loading (OCL) problem. Due to the interplay of continuous and logical constraints, OCL is an NP-hard problem, so that a variety of heuristic algorithms have been proposed in the literature. Herein, an algorithm for its exact solution, named X-OCL, is developed under the assumption that the chillers’ power consumption curves are quadratic. The proposed method hinges on a decomposition of the solution space so that the overall OCL problem is decomposed to a set of equality constrained quadratic programming problems that can be solved in closed form. By applying the new X-OCL solver to well known case studies, we assess and compare the performances of several literature algorithms, highlighting also some errors in the published results. Moreover, X-OCL is used to design a greedy optimal chiller sequencing (OCS) solver, called X-OCS. The X-OCS is tested on two literature benchmarks and on the model of the heating, ventilation and air-conditioning (HVAC) system of a semiconductor plant, over a two-year period. The performances of X-OCS are remarkably close to the theoretical optimal performance.


Introduction
Heating, ventilation and air-conditioning (HVAC) systems typically utilize a large percentage of the total building energy consumption, amounting to approximately 25-30% for dwellings [1] and up to more than 50% in industries relying on clean-rooms [2]. Between 40% and 60% of this requirement is due to chillers, which are therefore responsible for a significant fraction of total energy use. Such numbers indicate that the efficiency of HVAC systems is closely dependent on the efficiency of the chiller unit. Since a multiple chiller system typically employs machines, there are usually several combinations of chillers' part loads that are able to satisfy the load demand. The problem of determining the load fraction that each chiller has to deliver in order to minimize the system power consumption is known as optimal chiller loading (OCL) problem. In the last decade, several methods have been proposed to solve the OCL problem [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Chang et al. [3] assume that the chillers' efficiency is a quadratic function of the partialization and maximize the sum of chillers' coefficients of performance (COP) using the Lagrangian method. Of course, maximizing the sum of the efficiencies is not equivalent to the "canonical" OCL problem, where system power consumption is instead minimized. In a second paper, Chang [4] has addressed the canonical OCL problem again by the Langrangian multipliers methods, assuming that all the chiller power consumption curves are convex and cubic in the part load ratio. The poor convergence properties of the lambda iteration at low cooling load demand were obviated by a suitable gradient method.
Nevertheless, it is well-known that the chillers' power consumption could be a concave function, which makes OCL an NP-hard mixed-integer problem and poses a further challenge to the search of the optimal solution as global convergence of iterative methods cannot be guaranteed. This has motivated the development of several heuristic methods which do not guarantee optimality, but are effective to obtain a fair solution in a viable execution time.
Many of these algorithms have been tested on the Hsinchu benchmark, a widely used case study consisting of a six-chiller system installed in a semiconductor factory located in the Hsinchu Scientic Garden (Taiwan). Hence, it is therefore natural to use the performances achieved on this benchmark to compare and rank alternative solution methods. A selection consisting of the six best algorithms is reported in Table 1. In order of publication, they are IFA (2013), DCSA (2014), GAMS (2015), TLBO (2017), EIWO (2018), DCEDA (2020). Two recent algorithms, ICA (2020) and WGA (2020), claim to have improved over the existing literature. However, an inspection of the chillers' cooling capacities used in these two papers shows that they are different from those of the standard benchmark so that the claim of improvement is void. It is worth noting that, since an exact solver of the Hsinchu benchmark is not available, the final word on the optimality of solutions found by heuristic algorithms still needs to be had.
In the practical management of a real chiller system, solving the OCL is not sufficient, because there exist further dynamic constraints, namely minimum up-time and down-time requirements on chillers' operation. When these constraints are accounted for, the power consumption minimization problem goes under the name of optimal chiller sequencing (OCS). Again, this is a problem that is hardly tractable without resorting to some heuristics. In particular, the knowledge of all future cooling loads is required, which raises the problem of forecasting it with reasonable accuracy, a task that can be successfully addressed only on a finite prediction horizon. In the literature, OCS solvers resulting from the combination of heuristic OCL algorithms with dynamic programming schemes have been proposed [8,[23][24][25].
The present paper addresses both the OCL and OCS problems. Concerning the former one, two main issues are investigated. First of all, we derive an exact algorithm when the chillers' power consumption is a quadratic function of the partial load ratio (PLR), as it happens for the Hsinchu benchmark. Our X-OCL algorithm allows one to have the final word on the existing heuristic methods, highlighting also some erroneous results reported in the literature. The second issue has to do with the practical applicability of X-OCL to real-world plants. In particular, the execution time is compared with a state-of-art mixed integer solver and the adequacy of the quadratic power consumption model is discussed.
Concerning the OCS problem, we exploit the X-OCL algorithm to derive a lower bound on the minimum power consumption achievable by any OCS solver. Second, a greedy OCS algorithm leveraging on X-OCL is proposed and compared with [23]. Finally, the lower bound is used to quantitatively assess the degree of suboptimality ensuing from the lack of preview implicit in the greedy approach. Literature benchmarks, i.e., the Hsinchu one and two OCS benchmarks are used for test and comparison. Moreover, an extensive dataset collected during two years in a semiconductor fab is used to build and demonstrate a comprehensive solution of both OCL and OCS, including the data-based estimation of the chillers' power consumption models. In particular, the potential energy saving with respect with the current HVAC energy management is assessed.

The Optimal Chiller Loading Problem
Assuming n chillers operated in parallel, let Q i , i = 1, . . . , n, denote the cooling power delivered by the i-th chiller and P i = P i (Q i ), i = 1, . . . , n the associated power consumption. For a prescribed overall cooling load demand Q load , the goal of the OCL problem is finding the cooling powers Q i , i = 1, . . . , n that the n chillers have to deliver in order to minimize the system total energy consumption P tot : For each chiller, PLR (part load ratio) denotes the cooling load fraction, given by where Q 100%,i is the maximum power supplied under full capacity operation. The vector of all PLRs is denoted by When the i-th chiller is turned on, it should not operate under a minimum PLR, denoted by PLR min,i , 1. For the subsequent derivations, it is convenient to introduce a binary variable δ i that indicates the status of the i-th chiller and a real-valued variable x i , such that: The power consumption P i of the i-th chiller is assumed to depend mainly on PLR i and the condenser inlet water temperature T i , that is P i = P i (PLR i , T i ).
The consumption surface P i (PLR i , T i ) is either obtained from laboratory experiments or field data collected during operation.
Note that the problem of finding the optimal part load ratios PLR * can be formulated as a mixed-integer nonlinear program (MINLP): In the above problem, two types of constraints are present: the cooling demand constraint, Equation (1b) and a set of operational constraints Equation (1c-e) regarding the admissible operating regions of the chillers.
In the following, it is implicitly assumed that the cooling demand constraint is such that the admissible solutions set is nonempty. Then, given that the constraints define a closed set of admissible solution, the cost function admits a minimum.

Quadratic Power Consumption Model
With the exception of particular structures, mixed-integer programming problems are classified as NP-hard, which means that in the worst case, the solution time grows at least exponentially with the problem size. Although its combinatorial nature might suggest the use of heuristics, we will show that a significant subclass of industrial OCL problems, characterized by a quadratic power consumption model quadratic in PLR i , may still be successfully attacked by a carefully designed exact method. It is worth asking whether the quadratic assumption of the power consumption is more or less realistic. We will discuss this point later. Even if the consumption is not quadratic, a good quadratic approximation is likely computable. Moreover, in the case in which it is not easy to workout a good quadratic approximation, it can still be addressed by means piece-wise quadratic model, as suggested in Section 6. Assumption 1. The power consumption P i of the i-th chiller obeys the following model: where β p,i are the model parameters and f (·) is a suitable function. Moreover, it is assumed that β 2,i = 0.
When the OCL problem is solved in a given time slot, the condenser inlet water temperature T i can be assumed to be known. Then, for a given T i , in the interval [PLR min , 1], the consumption surface is a quadratic function of PLR i alone: with a i = β 0,i + f (T i ), c i = β 1,i , and q i = β 2,i . In view of the quadratic nature of the cost function, the system total energy consumption can be expressed in matrix form as follows: Moreover, the equality constraint (1b) can be rewritten as represent the standard OCL problem with quadratic power consumption, that so far has been attacked by either mixed integer programming solvers or a variety of heuristic algorithms. Differently from these approaches, below we derive an algorithm that computes the exact solution.

Partition of the Solution Space
In view of (4), the admissible set F for PLR i is where each subset F σ is associated with one of the following operating conditions: switched off (σ = 0), minimum part load (σ = 1), maximum part load (σ = 2), intermediate part load (σ = 3). In the sequel, σ i ∈ {0, 1, 2, 3} will denote the state of the i-th chiller.
In order to satisfy the operational constraints of a multiple chiller system, the solution of the OCL problem must be searched within the cartesian product of the chillers' admissible sets, i.e., PLR ∈ S, S = F n Since the admissible set F of a single chiller can be partitioned in four subsets, the overall admissible set S can be partitioned in 4 n subsets S j , j = 1, . . . , 4 n , each of which is in a one-to-one correspondence with the n−digit multichiller code formed by the state codes σ i , i = 1, . . . , n, of the n chillers.
To make an example, consider the case of n = 3 chillers. Then, the possible 4 n = 64 subsets S j are associated to the multi-chiller codes as follows: The n elements of the set S j are in a one-to-one correspondence with the chiller number i, 1 ≤ i ≤ n. Given the chiller number i, the notation S j [i] will denote the operating condition, either a point or a range, of the i-th chiller. For instance, within S 8 , whose multi-chiller code is s 8 = [0, 1,3], the three chillers i = 1, 2, 3 will operate at the following conditions:

Divide and Conquer Strategy
For j = 1, . . . , 4 n , let QP(j) indicate the OCL problem (4) restricted to the subset S j .
The 4 n problems QP(j) can be partitioned in two subsets, C andC: , ∀i}: all the partial load ratios are fixed so that S j has cardinality one; -C, when there is at least one chiller operating at intermediate part load (∃i : s ji = 3).
Consider, for example, S 7 , whose multi-chiller code is s 7 = [0, 1, 2]: chiller #1 is switched off, chiller #2 is operating at minimum part load, and chiller #3 is operating at maximum part load. Within this subset, no optimization is actually needed because all the chillers' PLRs are fixed and S 7 = 0 PLR min 1 T has cardinality one. Therefore, only a feasibility check is required: if constraint (8b) is satisfied, then S 7 is the optimal solution of QP(7). Otherwise, QP(7) does not admit a solution. It is easy to see that the number of elements of the subset C is 3 n .
In this case, QP(8) is a constrained quadratic programming problem in the unknown x 3 . As will be shown later, for j ∈C, the QP(j) problems enjoy a remarkable property: their optimal solution, if it exists, is a critical point and no more than one critical point exists. In the following, PLR * (j) will denote: the optimal solution of QP(j), if j ∈ C; 2.
the unique feasible critical point of QP(j), if j ∈C.
Let A denote the set of integers j s.t. PLR * (j) exists. The associated value of the cost function will be denoted by Then, the key idea is to exploit the partition of the solution space by a two-step procedure: 1.
In the next subsection, it is shown how to reduce the inequality-constrained problems associated with j ∈C to equality-constrained quadratic problems (EQP), for which a closed form solution is available.

Reduction to Equality-Constrained Problems
For a given j ∈C, we denote by V j = {i ∈ {1, 2, . . . , n} | s ji = 3} the set of chillers operating at intermediate part loads, i.e., between PLR min and 1. We also let κ(j) denote the cardinality of V j . It is then possible to rewrite QP(j) as a reduced-order quadratic problem in the κ(j) ≤ n unknowns is the cooling load that must be supplied by the chillers operating at intermediate part loads. Now, we associate to each QP(j) the corresponding equality-constrained quadratic problem EQP(j), that is obtained by removing the inequality constraints (10c).

Problem 2. EQP(j)
It is convenient to rewrite EQP(j) in matrix form. For this purpose, we introduce a selection matrix M(V j ) ∈ R κ(j)×n that selects κ(j) elements out of n.
Note that, being integers, the elements of V j admit an obvious ordering. Then, In short, M(V j ) will be denoted by M j . The reduced-order matrices are thus given by: Then, EQP(j) can be restated as: By applying the first-order Karush-Kuhn-Tucker (KKT) necessary condition to the EQP(j) problem, the following linear system is obtained: wherex * identifies a critical point, either maximum, minimum or saddle, and λ * ∈ R is the associated Lagrange multiplier. In order to guarantee the existence of the solution, a technical assumption is introduced.
Note that Assumption 2 is immediately satisfied if q i > 0, ∀i, that is when the quadratic power consumption curves (3) are all convex, although this is not necessary. Theorem 1. Under Assumption 1, EQP(j) admits a unique critical point x * , given by Proof of Theorem 1. In view of Assumption 1, det(Q) = 0, because q i = 0, ∀i. Moreover, Assumption 2 guarantees thatẼQ −1ẼT = 0. Then, it is immediate to see that the KKT condition (15) admits (16)- (17) as a unique solution.
Given the critical pointx * , it is easy to obtain a critical point for QP(j), as well. For this purpose, it is convenient to introduce an auxiliary vectorx ∈ R n×1 , such that: The candidate critical point for QP(j) is thus given by: The keystone of the solution procedure is the connection between the critical points of QP(j) and EQP(j), as stated in the following theorem.
Proof of Theorem 2. Sufficiency. Assume that the critical point PLR * for EQP(j) belongs to S j .
Observe that EQP(j) has less contraints than QP(j). Therefore, if a critical point for EQP(j) is feasible for QP(j), it is ipso facto a critical point for QP(j), as well. Necessity. Assume that QP(j) admits a critical point, say PLR * (j). Given that EQP(j) and QP(j) differ only for strict inequality constraints, any critical point for QP(j) is critical also for EQP(j). Since EQP(j) admits at most one critical point, necessity is proven.

Summary of the X-OCL Algorithm
We are now in a position to summarize the steps of the proposed algorithm, hereafter named X-OCL, for the exact solution of the OCL problem.
The partition of the solution set in the 4 n subsets S j , j = 1, . . . 4 n allows one to divide the MINLP problem into 4 n sub-problems QP(j) (8). For each subset S j , two situations can occur: 1.
all the chillers part load ratios are fixed; 2.
there is at least one chiller working at intermediate part load.
In the former case, only a feasibility check is required to decide whether the PLR configuration is to be kept as a candidate solution. In the latter case, the solution of the QP(j) problem is reduced to the solution of an EQP(j) problem that admits a unique critical point, easily computable in closed form. As stated in Theorem 2, if the critical point of EQP(j) belongs to S j , it coincides with the critical point of QP(j). Otherwise, QP(j) does not admit a critical point.
Once all the QP(j) have been processed, the corresponding set of critical points {PLR * (j), j = 1, . . . , 4 n } include the optimal solution of the overall OCL problem, which can be found just by comparing the associated system power consumptions P * tot (j). A pseudo-code summary is reported in Algorithm 1.
Algorithm 1: X-OCL. 1: Input: a, c, q ∈ R n×1 , s ∈ R 4 n ×n , Q load ∈ R 2: Output: PLR * , P * { feasibility check} 8: all the chillers part load ratios are fixed 9: x ← S j 10: if PLR * (j) ∈ S j then 25: Summarizing the difference between OCL and OCS is that the former one is a static problem, whereas the second one is intrinsically a dynamic problem. An OCS problem stripped of its dynamic constraints could be solved just as a sequence of OCL problems.

Hsinchu Cooling Plant Model
The Hsinchu chiller system, originally described in [7], has become a widely used benchmark for the testing and comparison of OCL algorithms [6,15,[17][18][19]. The case study involves six chillers installed in a semiconductor factory located in Hsinchu Scientific Garden (Taiwan) with a 7620 kW total cooling capacity. Quadratic models of the chillers' energy consumption were obtained and validated from data collected every 5 min over a 5-month period [7]. The benchmark problem assumes that the condenser inlet water temperature is 24.5 • C. The coefficients of the six chillers' energy consumption models (3) are reported in [7] (Table 3) and the corresponding P-PLR curves are displayed in Figure 1. It is asked to solve the OCL problem for five different cooling loads, ranging from 70% to 90% of the system total cooling capacity. It is also required that the partial load ratio of each chiller never goes below 0.3. According to our notation, the following parameter settings are used: a i , c i , and q i from [7]; -Q Load = 90%, 85%, 80%, 75%, and 70% of the chillers' maximum capacity (∑ n i=1 Q nom ); -

OCL Benchmark: Results
The comparison between the X-OCL optimal solution and those obtained by the six literature methods in Table 1 sheds light on some errors in the published results.
For the three highest cooling load demands (Q Load = {6858, 6477, 6096} [kW]), the X-OCL optimal solution coincides with the common solution provided by the six algorithms, thus confirming that they had reached the optimum. In the remaining two cases (Q Load = 5717 and 5334 [kW]), the solution computed by X-OCL coincides with those of GAMS, EIWO and DCEDA, which, however, are apparently outperformed by IFA and DCSA. The worst performance is that of TLBO.
However, concerning the solution provided by IFA and DCSA, four wrong values of power consumption (identified by asterisks in Table 1) were published in [15,17]: -P tot,IFA at load 75% and 70%; -P tot,DCSA at load 75% and 70%.
In fact, Table 1 reports P tot,IFA = P tot,DCSA ≈ 3507. 3 [kW] at 70% cooling load, which is inconsistent with the PLR i reported in the same table. Such inconsistency is easily verified by plugging the PLR i 's into the chillers' power consumption Equation (3) in order to obtain the individual chiller power consumptions of the P i [kW] column (the papers [15,17] describing algorithms IFA and DCSA do not report these individual consumptions). Apparently, the source of the error is the mechanical use of the chiller's quadratic power consumption model out of its operational range, that is, in correspondence with the null partial load ratio. Obviously, the associated power consumption is null as well, but if the quadratic model has a negative constant term, as is the case for chillers 3-5, the formula will return a negative power consumption as if a turned-off chiller could generate free power.
In Table 2, the published chillers' power consumption P p,i and the corresponding corrected values P c,i are reported for both IFA and DCSA. The values highlighted in red are the unfeasible chillers' power consumption, of which the correction is reported in column P c,i .
Once the errors have been corrected, IFA and DCSA are no more optimal at 75% and 70% loads.  Table 3.
In conclusion, on the Hsinchu benchmark, GAMS, EIWO and DCEDA prove to be the best heuristic algorithms, as their solutions coincide with the optimal ones computed by X-OCL, for all considered loads.

The Optimal Chiller Sequencing Problem
The cooling load demand of a building can be subject to significant variations during the day. Consequently, solving the OCL problem in each time step t (for example, 20 min), just ignoring minimum up/down time constraints on the chillers, could lead to frequent switchings (chiller startups and shutdowns). In order to preserve chillers from excessive mechanical stress and increase their operating life, each machine should not be switched off before a minimal up-time is reached. Analogously, it should not be switched on too quickly. To comply with these requirements, minimum up/down-time constraints must be enforced in the formulation of the so-called optimal chiller sequencing (OCS) problem. In its full formulation, OCS is a dynamical problem, because, in order to minimize the cumulative power consumption, the current decision should also take into account future constraints. As a consequence, the solution approaches proposed in the literature range from dynamic programming to heuristic methods were designed to alleviate the complexity of the problem.

A Lower Bound to the OCS Problem
Any solution to the OCS problem must face some level of approximation. Even when dynamic programming is used, there is the necessity of forecasting future loads, which introduces a suboptimality margin with respect to the ideal solution based on perfect knowledge of the future load profile. The availability of an easy-to-compute limit of performance against which the results of heuristic methods can be benchmarked is therefore of interest. In order to derive such a limit, one can consider the relaxed OCS problem (R-OCS), that is, an OCS problem without up-and down-time constraints. The relaxed OCS problem boils down to a sequence of independent OCL problems to be solved at each step in correspondence with the associated load. The availability of an exact OCL solver, such as X-OCS, makes it possible to compute the exact solution of the relaxed OCS as well. Notably, in view of the independence of the OCL problems, what matters is not the load sequence but the load distribution, so that the R-OCS bound could be easily derived based on statistical distributions reflecting different production and weather scenarios.
Such a bound can be used to quantitatively assess the existing margin of improvement for a given heuristic OCS solver. In fact, if the achieved power consumption is close enough to the bound, there is no scope for the search of further improvements. Along this direction, in the following section, the X-OCL solver is used to derive a greedy OCS algorithm, whose performance is then assessed against the R-OCS bound.

X-OCS, a Greedy Approach to OCS
X-OCS is a greedy algorithm that reduces OCS to a sequence of OCL problems, solvable through X-OCL. The approach is greedy because at each time step, future constraints are ignored, and the optimal OCL solution, compatible with the current minimum up/down time constraint, is searched for. In the mathematical form, the greedy OCS problem can be written as a mixed-integer quadratic problem with linear constraints: where It is easy to observe that the greedy OCS problem is an OCL problem with the two additional constraints (19e-f) which force some chillers to be online/offline depending on their previous states At each time step t, two situations can occur: all the chillers' states δ i (t) are free, i.e., all the chillers have been online/offline for more time steps than those prescribed by MUT/MDT; -there is at least one chiller, say the i-th one, whose state δ i (t) is constrained to be online or offline (δ i (t) = 1 or δ i (t) = 0) by the MUT or MDT.
Concerning the first case, the minimum up/down-time constraints are not active, therefore the step of the greedy OCS boils down to an OCL problem and its optimal solution can be found by the X-OCL algorithm. In the second case, instead, at each time t, the optimal solution is found by applying the X-OCL algorithm to a suitable subset of the OCL solution space S. At each time step t, the solution space of the greedy OCS is therefore given by: are the sets of chillers which must be on and off, respectively.
The idea is to exploit the partitions of the solution space P j by the typical two-step procedure of X-OCL:
Herein, A(t), which denotes that the set of integers j s.t. PLR * (j) exists, is a function of t because the load changes with time.

OCS Simulated Example
In this section, the performance of the X-OCS method is assessed on two case studies taken from the literature [23]. Case study 1 involves a hotel in Taipei with two 450 refrigeration tons (RT) chillers and two 1000 RT chillers, while case study 2, again from the Hsinchu Science Industrial District, features nine 1250 RT chillers. In both cases, chillers are described by their COP-PLR curves, expressed by a second-order polynomial model: where α i , β i and γ i are the chiller's coefficients, reported in Table A1.
The first step was the identification of quadratic power consumption models (3) using data sampled from the COP-PLR curves. Further details regarding the identification procedure are reported in Appendix A. The coefficients of the quadratic P-PLR curves in Figure 2 are reported in Table A2. The aim of both benchmarks is to compute the sequence of chillers' partializations over one day, assuming 20-min stages, so as to minimize the cumulative power consumption, while satisfying the cooling demand constraint at each stage. The load demand profiles are displayed in Figure 3.
The parameters' settings were as follows: - Recall that MUT i = 3 means that the i-th chiller must be on at least 3 consecutive steps before being turned off. Likewise, MDT i = 1 indicates that the i-th chiller, once turned off, must remain off at least 1 step.
The results obtained via X-OCS were compared with those obtained by dynamic programming, as reported in [23]. Where optimization was carried out under the ideal condition that all future loads are known in advance.
For sake of comparability with DP, the chillers power consumption associated with X-OCS was evaluated by plugging the PLRs computed by X-OCS into the original benchmark's COP model [23] and not via the approximate quadratic power consumption curves used by X-OCS.

OCS Benchmark: Results
The results are shown in Tables 4 and 5. In case study 1, the power consumption obtained by the X-OCS method was, at each stage, lower than or equal to that obtained by the DP method. From stage 1 to stage 29, the DP and X-OCS methods selected the same chillers. Since the MUT and MDT constraints were not active, the solutions coincided with the R-OCS ones. At stage 30, the R-OCS solution had chillers 4 and 3 switched on. However, the MUT constraint forced the DP and X-OCS to leave chillers 1 and 2 switched on. The same goes for stage 31. At stage 32, 44 and 59 the X-OCS method performed better than the DP one. Although the chillers were not constrained by minimum up/down time limits, apparently the DP method, as implemented in [23], could not find the global minimum. Notably, the greedy X-OCS method achieves a total electric power consumption (46,495.91 [kW]) which is less than 0.1 % greater than the minimum achievable bound R-OCS (64,432.56 [kW]), meaning that there is no practical margin of improvement left).
In the second case study, the cooling load demand varied slowly over time, so that the MUT and MDT constraints were never active. Therefore, the X-OCS, notwithstanding its greedy nature, attains the best achievable performance bound R-OCS. For the majority of the cooling loads, DP and X-OCS gave the same results, the only exceptions being at stages 22, 48-62, 68-72, where X-OCS performed slightly better. For both the case studies, the cumulative daily power consumptions obtained by DP (marked by asterisks) had been reported incorrectly in [23]: As a matter of fact, the daily power consumption values reported in the paper did not match with the sums of the power consumptions at each step, which we used for the correction.
So far, the performances of alternative methods have been compared on OCL and OCS benchmark problems whose quadratic consumption models were taken from the literature. Moreover, a limited number of loads were considered.
In this section, the feasibility of HVAC efficient management based on the exact solution of the OCL problem is validated against a real-world scenario that includes the estimation of the chillers' consumption models from the field data.

Field Data
The study was based on experimental data from a large HVAC system installed in a semiconductor fab located in Austria. The chiller unit is composed by five water-cooled centrifugal chillers, see Table 6, connected in parallel, in a constant primary flow chilled water system. Quarter-hourly data were recorded at different working conditions over a period of almost two years, from February 2017 to January 2019. Collected data include: temperatures, PLR, power consumption.
The time series of the cooling load demand is shown in Figure 4. The chillers were subject to the following operating constraints: -PLR min = 0.2; -MUT i = 4 (one hour), ∀i; -MDT i = 2 (half an hour), ∀i.

Chiller Energy Consumption Models
The chillers' power consumption models were estimated using the evaporator cooling capacity Q evap [kW] and the condenser inlet water temperature T [ • C], as covariates, and the compressor power consumption data P [kW], as target.
For each chiller, the dataset was randomly partitioned in two subsets: 70% for training and 30% for testing, respectively. The parameters β p,i of the model (2) were estimated via least squares fitting of the training data, discarding data with PLR < PLR min = 0.2. The values of the estimated parameters β p,i are reported in Table 6, together with the percentage coefficient of variation, defined as CV% = 100 × SE(β p,i )/|β p,i |. The cooling capacity Q nom are the nominal values provided by the manufacturer.
In Figure 5, the fitted surfaces P(PLR, T) are displayed against the validation data for each of the five chillers. It is seen that the the quadratic model, in spite of a few outliers that are not uncommon when data are collected in industrial frameworks, predicts the consumptions at different operating conditions well, as also confirmed by the goodness-of-fit (GOF) plots in Figure 6.   Figure 7 shows that T, i.e., the condenser inlet water temperature depicted in Figure 8, is mainly concentrated in a narrow range centered around the setpoint, namely 21.5 • C. Following what is usually done in the literature benchmarks, one could neglect temperature variations around the set point and solve the OCL and OCS problems using the chillers consumption curves at 21.5 • C, displayed in Figure 9. However, the inspection of Figure 5 shows that for chiller 4, the power consumption is significantly affected by the temperature, especially in summer. Therefore, differently from other literature studies, OCL and OCS solutions, we computed this based on the complete model P(PLR, T).

Real HVAC System: Assessment of Potential Savings
The field data were used to perform a retrospective analysis of the efficiency of the HVAC system management. More precisely, starting from the historical decisions and the associated cumulative power consumption, two comparisons were performed for the 2-year OCS problem. First of all, the lower bound R-OCS on the best achievable consumption was computed in order to quantitatively assess the potential improvement margin. Being a lower bound, R-OCS may be overly optimistic, so that it is important to evaluate the performance that can be obtained in practice. This was done by running the X-OCS solver, whose energy consumption could then be compared with the (ideal) R-OCS bound and the historically recorded power consumption.
The cumulative energy consumption recorded during the 2-year monitoring was 1.758 × 10 6 [kWh]. This figure can be compared with the R-OCS lower bound, equal to 1.600 × 10 6 [kWh]. This means that the potential margin of improvement is not larger than 8.97%.
When the X-OCS algorithm was applied, the cumulative energy power consumption was 1.601 × 10 6 [kWh]. As a matter of fact, for this HVAC system, the loss of performance due to the suboptimality of the greedy algorithm is definitely negligible (it is less than 0.1%). On average, the power saving is 0.157 × 10 6 /(2 × 365 × 24) = 8.96 [kW]. In Figure 10, the difference between the consumption achieved by X-OCS and the lower bound is plotted on a weekly basis. The black and blue lines overlap almost perfectly, so that they cannot be distinguishable.

Execution Time
The X-OCL and X-OCS algorithms, coded in Matlab R , were executed on a standard laptop (Intel(R) i7-7500U dual-core with hyperthreading, RAM 16GB, 2.7 GHz). While no explicit parallelization of the algorithm was implemented, the solution of the 4 n QP problems was formulated as a unique algebraic computation using sparse matrices. This means that the algorithm may have benefited from some optimization automatically enforced by the Matlab R compiler. In order to compare the numerical efficiency, the Hsinchu benchmark was also solved using the CPLEX solver for costrained mixed integer problems as implemented in the GAMS environment, used by [6]. The results concerning the computational times are summarized in Table 7.
Concerning the Hsinchu OCL benchmark, the computational cost of X-OCL was negligible (0.11 s) and, more importantly, definitely lower than the time (1.991 s) spent by GAMS.
Concerning the three OCS benchmarks, in all cases, the X-OCS computation time was larger than that for R-OCS, because the greedy algorithm performs the additional task of looking for solutions that satisfy the up-/down-time constraints.
In the OCS examples, the larger number of loads makes it possible to assess the average computation time per load. For a small number of chillers, the average computation time per load of X-OCS is remarkably small: 0.368/72 = 5 × 10 −3 s for OCS benchmark 1 (4 chillers) and 573.96/68, 110 = 8.4 × 10 −3 s for the OCS on field data (5 chillers). In the latter case, X-OCS took less than 6 min to solve the OCS problem over 2-year data with quarter-hour sampling.
As expected, in view of the exponential growth of the number of QP problems, the maximal average computation time per load is found in correspondence with the OCS benchmark 2, where 9 chillers are present. Nevertheless, the average time per load amounts to 228.29/72 = 3.17 s, which is totally affordable.
Finally, a simulated benchmark was set up in order to better assess the average computation time per load as a function of the number of chillers. In particular, nine OCL simulated experiments, characterized by a variable number of chillers, ranging from two to ten, were considered. The ten-chiller system was composed by the six Hsinchu chillers, whose models are reported in [7] (Table 3), with the first four replicated. For each experiment, the OCL problem was solved for five different cooling loads ranging from 70% to 90% of the chiller's total cooling capacity. In Figure 11, the average computation time per load is displayed as a function of the number of chillers. The plot highlights the exponential growth but, at the same time, it shows that, even for a large HVAC system made of 10 chillers, the computation time per load (10 s) is all but prohibitive.  Figure 11. Simulated X-OCL experiment: average computation time per load.

Discussion
The main purpose of the paper was the derivation of an exact algorithm for the solution to the OCL problem. This goal was successfully completed by a decomposition approach that exploits a suitable partition of the solution space. The availability of an exact method has been immediately exploited along two directions. First, it has become possible to give a definitive assessment on the performances of different literature methods that had been applied to some consolidated benchmark problems. By the way, the comparison with the exact solution revealed that some unrealistic performances had been declared in the literature, due to erroneous extrapolations of the power consumption curve. The exact method has been exploited also in relation with the optimal chiller sequencing problem. For a given cooling demand profile, if the dynamic up/down-time constraints are neglected, a sequence of OCL problems can be exactly solved to yield a lower bound limit, called R-OCS, to the performance of any feasible solution.
Concerning the practical applicability of the exact solution in place of heuristic approaches discuss in the literature, two main objections may be raised: (i) the explosion of the computational cost with the number of chillers, (ii) the need to make overly restrictive assumptions on the shape of the power consumption curves.
On the first side, it is indeed true that the OCL problem as formulated in (4) is NP-hard, as confirmed by Figure 11, where the exponential growth of the computations time is apparent. At the same time, the figure shows that, even for a medium/large-sized chiller system, e.g., 6-10 chillers, the computation time per load is less than 2 s. This suggests that, even for chiller systems used in large semiconductor factories, the cost of the exact solution is all but prohibitive. The numerical efficiency is a direct consequence of the tiny number of computations required to solve the elementary EQP problems, see (12). Moreover, the partition strategy underlying X-OCL implies that it is totally parallelizable into 4 n threads, a feature that has not been explicitly exploited and that could further speed up the solution.
The second objection to the general applicability of the exact solution has to do with the quadratic assumption made on the power consumption curve. Even if the majority of benchmarks satisfy this assumption, there is no stringent reason to rule out other functional descriptions. Nevertheless, when confronted with real data (see the field data benchmark), we found that a quadratic power consumption fitted well the recorded data, see Figure 6. Moreover, even when a single quadratic function is not adequate, it is still possible to switch to a piece-wise quadratic description. In that case, it would be immediate to generalize the exact algorithm by just increasing the number of partitions S j , in such a way that the problem still boils down to a set of easy-to-solve EQP problems.
In view of its numerical efficiency, the application of X-OCS to the solution of the OCS problem also appears very promising. Indeed, in the OCS benchmarks taken from the literature and in the field data OCS benchmark, the performance of X-OCS is very close to the lower bound, implying that there is no scope for the use of more sophisticated algorithms.

Conflicts of Interest:
The authors declare no conflict of interest.

3.
Using the training set made of input-output paired samples {PLR i (k), P i (k)}, k = 1, . . . N estimate the parameters of the quadratic model via ordinary least squares (OLS).
The above procedure was applied to the case studies presented in Section 4.3. For both the case studies, 50 data points were randomly sampled (uniformly in [0.5,1]) from the COP-PLR curves reported in Table A1 to obtain the training datasets represented in the left panels of Figures A1 and A2 as red dots. The estimated P-kW curves are shown in the right panels of Figures A1 and A2, and their parameters are reported in Table A2.