Model Predictive Control of Central Chiller Plant With Thermal Energy Storage Via Dynamic Programming and Mixed-Integer Linear Programming

This work considers the optimal scheduling problem for a campus central plant equipped with a bank of multiple electrical chillers and a thermal energy storage (TES). Typically, the chillers are operated in ON/OFF modes to charge TES and supply chilled water to satisfy the campus cooling demands. A bilinear model is established to describe the system dynamics of the central plant. A model predictive control (MPC) problem is formulated to obtain optimal set-points to satisfy the campus cooling demands and minimize daily electricity cost. At each time step, the MPC problem is represented as a large-scale mixed-integer nonlinear programming problem. We propose a heuristic algorithm to obtain suboptimal solutions for it via dynamic programming (DP) and mixed integer linear programming (MILP). The system dynamics is linearized along the simulated trajectories of the system. The optimal TES operation profile is obtained by solving a DP problem at every horizon, and the optimal chiller operations are obtained by solving an MILP problem at every time step with a fixed TES operation profile. Simulation results show desired performance and computational tractability of the proposed algorithm. This work was motivated by the supervisory control need for a campus central plant. Plant operators have to decide a scheduling strategy to mix and match various chillers with a thermal energy storage to satisfy the campus cooling demands, while minimizing the operation cost. This work mathematically characterizes the system dynamics of a campus central plant and establishes a linear model to predict campus cooling load. It proposes a model predictive control (MPC) strategy to optimally schedule the campus central plant based on plant system dynamics and predicted campus cooling load. A heuristic algorithm is proposed to obtain suboptimal solutions for the MPC problem. The effectiveness and efficiency of the proposed approach are well demonstrated for the central plant at the University of California, Irvine.


I. INTRODUCTION
B UILDINGS are one of the primary consumers of energy in the United States. Buildings are responsible for 41% of primary energy consumption in 2010, compared to 28% by the industrial sector and 31% by the transportation sector [1]. Among all energy consumed by the buildings sector, residential buildings accounted for 54% and commercial buildings accounted for 46%. Among all building energy consumers, heating, ventilation, and air conditioning (HVAC) accounted for the largest share, which represented close to half of all energy consumption [1]. Energy management of HVAC systems becomes essentially important to improve the energy efficiency and reduce the energy cost of buildings.
A central chiller plant is a commonly used HVAC system for a campus with a large number of buildings [2]. The campus cooling loads can be huge and vary in a large range depending on the weather condition and number of buildings. A bank of multiple chillers is usually installed to meet various cooling demands with high operational flexibility. They may have different performance characteristics and are often operated in ON/OFF modes [3], that is, a chiller can be turned on with the maximum mass flow rate of chilled water supply, or be turned off with zero rate. Thermal energy storage (TES) can be used to shape the campus cooling load and reduce the peak one [2]. In the simplest case, TES is a large tank of chilled water. A typical campus central plant with a chiller bank and a TES is depicted in Fig. 1. The chiller bank can either supply chilled water directly to the campus for cooling, or store some of chilled water in TES for later use [4], [5]. Intuitively, TES is charged by the chiller bank at night when electricity rates and cooling demands are both low. During on-peak hours in daytime, it discharges the stored chilled water for campus cooling, thus reducing the peak load of chillers. A typical control system for the central plant has a two-level control structure: supervisory and local ones [6]- [8]. The goal of supervisory control is to optimize TES and chiller bank operations for minimizing energy consumption or electricity cost while satisfying campus cooling demands. The operation set-points are obtained by taking into account of system level characteristics and interactions among all subsystems. The local control modulates individual components, e.g., a relay, valve, or actuator, to track the set-points provided by the supervisory control. A closed-loop feedback control, e.g., PID control, is usually used to achieve the set-point tracking performance and guarantee the robust operation of components. The overall energy efficiency and control effectiveness of a central plant is determined by the performance and coordination of controls in both levels.
Well-designed control systems can help realize the full potential of both energy and cost savings for a central plant [8]- [11]. The focus of this paper is on the design of optimal supervisory control using a model predictive control (MPC) scheme. The main idea of MPC is to use the system model of a central plant to predict its future evolution and obtain its optimal control in a moving horizon basis. The objective is to find optimal set-points by using an MPC scheme to minimize the operational cost while satisfying campus cooling demands and system constraints. At each time step, a model-based optimization problem is solved over a finite horizon to obtain optimal set-points of operating the TES and chiller bank. More details on MPC theory can be found in [12].
For the supervisory control of a central chiller plant, MPC techniques have been investigated since the late 1980s [13], [14]. More early studies on modeling, prediction, and optimization for building HVAC systems can be found in [6], [7], [15], [16] and references therein. With the growing interests in optimal control, many new methods have recently been used to develop MPC schemes for the central chiller plant with TES [17]- [20]. For instance, the work [17] proposed a method by using dynamic programming to solve the optimal chiller sequencing problem. The potential energy savings of the proposed method was demonstrated with three example systems. But the chiller considered in [17] was not equipped with a TES. In [18], [20], the authors considered the MPC problem of a series of chillers equipped with a water tank-a similar setup as this paper. They validated the system model with experimental data and performed closed-loop experimental tests for the proposed control method. To simplify the optimization problem for MPC, they assumed that the mass flow rates of chilled water can be continuously adjusted for chillers, and the charging/discharging profile for TES is fixed or known in advance. Similar assumptions for chillers and TES have also been made.in [19].
However, in practice, the operation decisions for chillers are often ON/OFF scheduling sequences and TES operation profile varies depending on the actual cooling demands. The optimization/control frameworks proposed in previous studies do not directly apply for supervisory control of a central chiller plant with TES under practical considerations. To generate optimal ON/OFF decision sequences for chillers and charging/discharging decisions for TES, a mixed-integer nonlinear programming (MINLP) problem can be formulated and solved at each time step using an MPC scheme. The nonlinear constraints result from the system dynamics and integer decision variables are due to the operational constraints for chillers and TES. Since this MINLP problem involves considerable nonlinearity and a large number of integer variables, finding its global optimal solution is nontrivial and computationally expensive [21], [22]. In practice, it is of interest to develop efficient algorithms to approximately solve it.
In [23], we present a heuristic algorithm to search for the suboptimal solutions of the MINLP problem by mixed-integer linear programming (MILP). Its main idea is to explore the special structure of the original problem and solve a set of subproblems with reduced complexity. We formulate an MILP problem to find optimal chiller operations by linearizing the system dynamics with a fixed/given TES operation profile. Compared with [23], we newly propose a dynamic programming (DP)-based algorithm to find optimal TES operations. Now the optimal TES and chiller operations are obtained in two steps. In the first step, we generate an optimal TES operation profile for a fixed horizon via DP with system dynamics approximated by a finite state machine. In the second step, we fix the TES operation profile for a fixed future horizon, and solve for the optimal chiller operations via a reformulated MPC problem. The original MINLP problem is approximated by the reformulated MPC problem with the linearized system dynamics and fixed TES operation profile, which leads to a mixed integer linear program (MILP) at each time step. The nominal trajectory for linearization is obtained by simulating the system with the optimal predicted control trajectory calculated from the last time step MPC problem. The major advantage of the proposed algorithm is that both DP and MILP problems can be solved very efficiently with guaranteed performance. Moreover, the model mismatch due to the linearization can be compensated by the moving horizon scheme used in the MPC scheme. A description of the system model for a chiller bank and TES is provided in Section II. The MPC problem is formulated in Section III. The DP-based heuristic algorithm for obtaining the optimal TES operation profile is presented in Section IV. The MILP-based heuristic algorithm for obtaining optimal chiller operations is described in Section V. A case study is provided in Section VI, and conclusions are drawn in Section VII.

II. SYSTEM MODELING
The schematic diagram of a typical campus central plant is depicted in Fig. 1, with all pumps being omitted.

A. Nomenclature
We use to denote the derivative of variable with respect to time , and to denote a rate related input or constant, e.g., chilled water mass flow rate. Time-varying variables are indicated with appended.
Density of water (kg/m ).
Specific heat capacity of water (kJ/kg K).
Ambient humidity ratio.
Cross-layer area of TES tank (m ).
Total number of chillers.
Coefficient of Performance of the th chiller.
Chilled water mass flow rate of chiller (kg/s).
Chilled water supply temperature of chiller (K).
Chilled water return temperature of chiller (K).
Cooling amount provided by chiller (kW).
Chilled water circulation through TES (kg/s).
Top layer 'a' temperature of TES (K).
Bottom layer 'b' temperature of TES (K).
Chilled water mass flow rate of chiller bank (kg/s).
Chilled water supply temp. of chiller bank (K).
Chilled water return temperature of chiller bank (K).
Electricity consumed by chiller bank (kW).
Cooling amount provided by chiller bank (kW).
Cooling load of campus (kW).
Chilled water mass flow rate to campus (kg/s).
Chilled water supply temperature to campus (K).
Chilled water return temperature from campus (K).

B. Chiller Bank Model
A chiller bank is composed of electrical chillers operated only in ON/OFF modes. When a chiller is ON, it supplies chilled water at a fixed maximum mass flow rate and supply temperature; When OFF, it is shut down and supplies no chilled water. For , the set-point of the th chiller is a binary decision variable if chiller is ON if chiller is OFF. (1) The amount of cooling provided by chiller is (2) Note that both and are fixed constants, and varies according to the system dynamics.
The total electricity consumed by a chiller depends on several factors, e.g., supply water mass flow rate , generated cooling amount , and circulated water temperature difference . Note that is a fixed constant when chiller is on, and is a bilinear function of and . Thus the chiller electricity consumption mainly depends on temperature difference . Here, we simply use a linear regression model to estimate the total electricity consumption of a chiller as well as chilled water pump, condenser water pump, and cooling tower fans. Following [19], the electricity consumption of the th chiller is modeled as where coefficients and can be obtained by regression analysis of the actual measured chiller data. The Coefficient of Performance (COP) is employed as a metric to quantify the energy efficiency of a chiller as follows: The total chilled water mass flow rate supplied by the chiller bank is The total power consumed by the chiller bank is (6) The total cooling amount provided by the chiller bank is The supply and return temperatures of the chiller bank are denoted by and , respectively. All chillers are assumed to have the same and , i.e., for , and (8)

C. TES Model
We employ a stratified two-layer TES model developed in [19], where the state of TES is captured by the top-layer water temperature and the bottom-layer one . The model is further simplified by ignoring the time delays of the TES tank output. We also assume that the TES tank is always full, implying that the mass flow rate entering the tank is equal to the mass flow rate exiting the tank.
The TES tank operates in two modes. In the charging mode, chilled water from the chiller bank is supplied to both the TES and campus, as shown in Fig. 1, while in the discharging mode chilled water from the TES and chiller is supplied to the campus together. The mode is denoted by a binary variable When charging, i.e., , TES has the following dynamics: When discharging, i.e., , it has Parameters , , , and are time constant multipliers, and and are inter-layer overall heat transfer coefficients, which can be calibrated using the actually measured TES data.

D. Campus Load Model
The campus cooling load satisfies the following energy balance equation: (12) Note that , , and are taken as system variables, while is taken as a given exogenous input to the system from a linear regression model described in the Appendix.

E. Model Summary
By collecting system dynamic equations (1)- (12) and discretizing the system with the sampling time , we represent the system dynamics in a state-space form as where the state variable of the system is denoted by the intermediate variable is denoted by the control variable is denoted by and the input variable is denoted by Equation (13a) corresponds to the TES dynamic equations (10d), (10e), (11d), and (11e); (13b) corresponds to the TES and campus balance equations (10a)-(10c), (11a)-(11c), and (12); (13c) corresponds to the chiller balance equations (2)-(8).
Besides the system dynamics, we also impose the various constraints on system variables, i.e., there exist bounded sets , , and such that , , and . These constraints are either due to the equipment limitations, e.g., the allowed minimum/maximum mass flow rates, or the constraints to the operation decisions, e.g., the binary constraints for the chiller and TES operation decisions.

A. Problem Formulation
The objective of the optimization is to obtain a real-time optimal scheduling strategy for a chiller bank and TES to minimize the total electricity cost while satisfying the campus cooling demands and system constraints. We formulate an MPC problem to search for the optimal set points of the scheduling strategy. The main idea of MPC is to use the system model of the central plant to predict its future evolution and obtain its optimal control on a moving horizon basis.
At each time step , we are interested in solving the following finite-horizon optimization problem regarding the state and control actions at time steps , where is the prediction horizon of the MPC problem: (14) where denotes the state variable for the time step predicted at time step (similar notations also apply to other variables); , , and denote the predicted state, intermediate, and control sequences over the predicted horizon; electricity rate is assumed to be known for all time and denotes the electricity power consumed by the chiller bank; the cooling load is also the known input from the forecasting model described in the Appendix, i.e., ; denotes the actual state of the TES at time ; and denotes the terminal constraint set for the TES state, which is chosen to ensure stability and feasibility of the MPC problem [12], [18].
By assuming that the maximum campus cooling load never exceeds the maximum capacity of a chiller bank, there always exist feasible solutions for problem (14). Among all feasible solutions, let denote the optimal one for problem (14). At time step , only the first element of is implemented to the system, i.e., . At time step , the horizon of interests moves one-step forward to , and the optimization problem (14) is solved again by changing index from to . The initial point of the predicted state sequence is taken as , where denotes the newly measured state at time step . This procedure is repeatedly implemented by moving forward the horizon step by step, which yields a moving horizon control strategy.

B. Conventional Approaches Based on MINLP
For the MPC optimization problem (14), at each time-step , there are in total decision variables to be computed, of which there are binary on, and real on. We also note that the constraints due to (13) have nonconvex terms due to the bilinear system dynamics. Thus, the optimization problem (14) is a nonconvex MINLP problem, which is extremely difficult to solve [22].
There are many software packages that solve nonconvex MINLP problems to proven optimality, such as , , and [22]. But they require many computing resources and a large amount of computation time to achieve a certain global optimality. On the other hand, many convex-MINLP solvers can provide heuristic solutions to nonconvex MINLP in a more efficient way, e.g., , , and [24]. But they suffer from the high sensitivity to the initial starting point. All of these solvers can therefore be used off-line to compute the optimal set points but not suitable for real-time implementation of an MPC scheme.

C. Proposed Approach Based on DP and MILP
In this paper, we propose a heuristic algorithm to approach the nonconvex problem (14) based on DP and MILP. Our idea is to explore the special structure of the original problem and solve a set of subproblems with less complexity. The goal is to obtain suboptimal solutions to (14) in a time-efficient manner with certain performance guarantee. The proposed algorithm is described below.
• Generating optimal TES operation profiles using DP: We search the optimal TES and chiller operations in two steps.
In the first step, we generate an optimal TES operation profile for a fixed future horizon by using DP. The system dynamics is approximated by a finite state machine (FSM) to reduce the computational complexity. Then, we fix the TES operation profile for the fixed future horizon, and solve for the optimal chiller operations by using a reformulated MPC problem.
• Linearizing system dynamics: The system dynamics is linearized along the predicted trajectory of state and control, where the trajectory is generated by simulating the nonlinear system model. The linearized dynamics can well approximate nonlinear one around the operating points. • Obtaining optimal chiller operations by solving an MILP problem: In the second step, we reformulate the MPC problem with the linearized system dynamics and fixed TES operation profile, which leads to an MILP problem at each time step. It can be solved very efficiently with certain performance guarantee. • Generating nominal trajectories for linearization: By simulating the nonlinear system model with the predicted control sequence obtained from MILP, we obtain the predicted trajectories of state and control for the next horizon. They are then used as nominal values in the next time step to linearize the system dynamics. The conceptual framework of the proposed approach is shown in Fig. 2 with its details to be provided next.

IV. DP-BASED ALGORITHM FOR OPTIMAL TES OPERATION PROFILE
Here, we describe a DP-based algorithm that generates an optimal TES operation profile . In Section V, will be fixed in the MPC problem (14) for a fixed future horizon to solve for the optimal chiller operations. The basic idea is to quantize the state variable pair ( ) and then approximate the differential equation-based TES model with an FSM. Then we can obtain a sub-optimal TES operation profile with drastically lower computational complexity.

A. Approximating TES Dynamics by FSM
An FSM consists of a set of states including an initial state, an input alphabet, and a transition function that maps input symbols and current states to the next ones [25]. We use the following notation to define it: , where denotes the set of states, denotes the input alphabet, also termed as "actions," denotes the transition function of the state dynamics, and denotes the initial state. If is the current state and is the current action signal at the th time step, then the next state . To develop a simple FSM model to approximate the chiller and TES system dynamics, we impose the following assumption only in this section.
Assumption 1: Both chilled water supply temperature and return one are constants. In reality, Assumption 1 is approximately satisfied with small fluctuations on and . This assumption helps develop a simple FSM model to approximate chiller and TES dynamics.
The set of states in FSM is defined as the collection of discretized states of TES: . Recall that a TES state is denoted by , where is the top layer water temperature and is the bottom layer water temperature. By partitioning the interval uniformly into levels, we obtain the discretized version of and , denoted as and , respectively. We add a further constraint that to be consistent with the reality. After the discretization, the cardinality of set is . The action of FSM at time step is taken as , a -by-1 vector containing the ON/OFF set-points for chillers. The action set is the -dimensional binary vector space with the cardinality . The transition function is a map that takes each pair to some state . This map is obtained by approximating the chiller and TES dynamics given by (10) - (11). The procedure to obtain is described as follows: At time step , given the action , variables , , and are completely obtained via (2) -(7) under Assumption 1. Note that the campus cooling load is a given/known exogenous input predicted from a linear regression model (Appendix A). The TES operating mode is then obtained as if if .
Given the current FSM state , we then compute variables , , , and . If , i.e., TES is being charged, we compute these variables as follows:

If
, i.e., TES is being discharged, we compute The next FSM state is obtained by evolving discretized TES dynamics (10d)-(10e) or (11d)-(11e) one time-step forward with the above computed variables. while satisfying all constraints are marked as "feasible" ones. The most cost-effective action at state is found and recorded in . The corresponding cost value is computed and stored in . Based on the outputs generated from Algorithm 1, the optimal TES operation profile can be easily derived. This is achieved by traversing the optimal state-action pairs forward in time, starting from initial state . This procedure is summarized in Algorithm 2. 13: end while 14: return .
One major advantage of a DP-based approach to choose a TES operation profile is that, unlike in [18], [19] where the TES operation profiles are fixed as a priori, they are decided by an optimization algorithm by taking the predicted campus cooling load into consideration. The complexity of the DP-based algorithm depends on the number of FSM quantization levels and the number of chillers , i.e., the algorithmic complexity is of the order . Note that the number of chillers and the length of horizon are fixed here. With more quantization levels, we can have a better approximation of TES dynamics. By varying the number of quantization levels , we can make a tradeoff between its complexity and accuracy to approximate the TES dynamics.

V. DECIDING OPTIMAL CHILLER OPERATIONS
Here, we present an MILP-based method to optimally operate the chiller bank in real time. Its main idea is to reformulate an MPC problem with linearized system dynamics and a fixed TES operation profile. The latter is generated via the DP algorithm. This formulation leads to an MILP problem, which can be solved very efficiently at every time step.

A. Linearizing System Dynamics
It appears that the major difficulties of solving a nonconvex MINLP problem are due to their non-convexity, which prevents the implementation of efficient algorithms [22]. For problem (14), nonconvexity mainly comes from the bilinear terms in the dynamics of the chiller bank and TES. Note that all bilinear terms appearing in (13c) are in the form of a binary valued variable and a real valued function. We apply Property 1 to convert the bilinear equation (13c) equivalently to a set of linear inequalities

B. Formulating an MILP Problem
Next, we define a simplified MPC optimization problem with a linearized model. The TES operation profile is fixed as a known input to reduce the computational complexity. We fix the TES operation profile as , where the TES operation profile is generated every time steps using DP as described in Section IV. At time step , if , we set and generate via Algorithm 1; otherwise, if , we use generated before. In the simplified MPC formulation, is then obtained at each time step as At each time step , we choose a nominal trajectory of state and control . The linearized system dynamics for the predicted horizon is obtained along the nominal trajectory of the state and control. Then, we consider the following optimization problem as an approximation to problem (14): (15) Problem (15) is an MILP problem, with binary variables and real on. Due to its convexity, it is much easier to solve than the nonconvex MINLP problem in general. Many software packages can efficiently solve an MILP problem with certain performance guarantee, e.g., , , , and [24], [26].

C. Generating Nominal Trajectories for Linearization
Recall that, to solve problem (15) at each time step , we require a nominal trajectory of the state and control for the predicted horizon . In this paper, we generate it through the simulation of the nonlinear system model (13) with optimal predicted control computed from problem (15). Let denote the optimal solution to (15) at time step . For the next predicted horizon , we heuristically choose the following control sequence as the nominal trajectory of control: Then, we apply the control sequence to simulate the nonlinear system model (13) with the predicted load input and initial state . The simulated state sequences and are taken as the nominal trajectory of states for the predicted horizon . Note that the initial trajectory of the state and control at time step 0 can be generated with any control sequence based on the experience or other heuristics.
The MILP-based algorithm is described in Algorithm 3 to compute suboptimal system control set-points of the chiller bank. Note that the TES operation profile is determined via Algorithm 1 for every horizon.

VI. CASE STUDY
Here, we present a case study to demonstrate MPC of a central chiller plant in summer time at the University of California, Irvine (UCI). We first use the measured data from the UCI campus to calibrate the parameters of the system model described in Section II. Then, the system model is built up in the MATLAB environment. An MPC strategy is compared with a baseline strategy through MATLAB simulations.

A. System Setup
The UCI central plant is depicted in Fig. 1. A bank of chillers supplies chilled water to satisfy the campus cooling demand. The supply water temperature of all chillers is the same and maintained at a constant value 40 F. 1 The return water temperature is the same for all chillers, and varies according to the system dynamics. For the MPC strategy considered below, we constrain the return water temperature in the range of F to ensure all chillers to work with high efficiency. They can only work in an ON/OFF mode, either with the maximum or zero mass flow rate. The maximum ones are fixed constant values and summarized in Table I. The parameters of the chiller electricity model (3) are calibrated through the regression analysis of the available measured data, which are also summarized in Table I. We calculate COP of each chiller using model (4) with parameters given in Table I. The calculated COP is depicted in Fig. 3 versus return water temperature in the range of interests. Intuitively, we observe that the seventh chiller has the highest efficiency and the fifth chiller has the lowest one in most range of return water temperature.
Chilled water supplied by the chiller bank can also be stored in a water tank TES for shifting the campus cooling load. The cross-layer area of TES is 564.95 m with a height of 30.48 m. At the initial time, TES is assumed to be fully discharged with the uniform water temperature of 60 F. The parameters of the TES model in (10) and (11) are calibrated through the regression analysis of the available measured TES data. The calibrated time-constant multipliers are , , , and . The calibrated inter-layer overall heat transfer coefficients are 0.0283 kW/(m K and 0.0197 kW/(m K . The campus cooling load is modeled and predicted using an ARX (autoregressive exogenous) model, where the exogenous inputs are taken as ambient temperature and humidity ratio (see the Appendix for more information). The parameters used in cooling load model (17) are trained using the actually measured campus load and weather data in the month of June 2009. The values of trained parameters are omitted here due to the space limitation, while the training results are depicted in Fig. 4. We observe that the actual cooling loads are well approximated by the cooling load prediction model (17). In the MPC problem, we employ the trained cooling load model (17) to predict the future cooling load based on historical campus cooling load data and weather forecast data. The electricity rate schedule for summer hours are obtained from Pacific Gas and Electric Company and summarized in Table II, for simplicity only the Time-Of-Use (TOU) energy charge is used in the cost function and the demand charge part is not included.

B. Two Scheduling Strategies
Here, we consider two strategies to schedule the operation of TES and chiller bank. One is a greedy-search-based heuristic approach which serves as a baseline strategy, and the other is the proposed DP and MILP-based MPC strategy. To simplify the comparison study, in both strategies, the TES charging/discharging operation profile is generated by using Algorithms 1 and 2 with horizon and desired terminal state 58 F 58 F . The system dynamics are approximated by an FSM quantized with levels. With the campus cooling This suggests that TES should be charged for the first 9 hours and switched to a discharge mode thereafter. 1) Baseline: This strategy generates suboptimal chiller and TES operations based on a greedy search algorithm. All possible chiller ON/OFF modes are first ordered according to their efficiency. Then, under a given TES charging/discharging profile, an initial feasible solution is obtained which requires more than necessary chillers staged ON. Starting from this initial solution, a greedy algorithm is applied to search for better chiller operation sequences that most drastically reduce the TOU-based cost function, until some given constraints are violated. This strategy has been tested with historical operation data and is shown to achieve 30% saving over the plant operator's heuristic operation in [19].
2) MPC: For the MPC strategy, we fist consider the MINLP approach to solve the MPC problem (14) with 192 binary variables and 692 real variables at every time step. Due to the large 1 1   number of variables and many nonlinear equations, the conventional nonconvex MINLP solvers (e.g., ) can not provide any feasible solution or does not converge in a reasonable time frame (up to 5 h). The convex-MINLP solvers (e.g., ) can always heuristically generate some feasible solutions in around 20 min but always keep running without stopping up to 5 h (this means it is still searching for better solutions). For a comparison purpose, we record those solutions obtained from software by stopping it when the computation time reaches one hour.
Then, we consider the proposed approach for solving the MPC problem based on DP and MILP. As discussed above, the optimal TES operation profile is generated using DP. By fixing as a known input, we formulate an MILP-based MPC problem (15) to solve for optimal chiller operations. The MILP-based heuristical Algorithm 3 is employed to search for the optimal set points of the system. The sampling time 1 h is used to discretize the system model. The prediction horizon of the MPC scheme is taken as 24 h. At each hour, the software package is employed to solve the MILP problem (15) with performance gap . For the MILP problem considered in this paper, it usually takes less than 10 minutes only to achieve the performance requirement. The initial nominal trajectory used by Algorithm 3 is obtained by simulating the system with the baseline strategy. By comparing the solutions obtained using two approaches, it is found that the solutions obtained through DP and MILP are always better than the solutions obtained by convex-MINLP approach. In Section VI-C, we will compare the baseline and MPC strategies with the solution computed via DP and MILP approaches.

C. Comparison of Two Strategies
The performance of both strategies is compared through MATLAB simulations with historical weather data for a specific day in June 2010. The total simulation time is 24 h. The predicted campus cooling load profile is depicted in Fig. 5(a) overlaid with the electricity rate for summer hours. We duplicate the campus cooling load profile for the next day to simulate the MPC scheme.
The controls computed from two strategies are applied to the simulation model of the UCI central plant. The chilled water temperature through the campus and the chiller bank are depicted in Fig. 6(a) and (b), respectively. The top and bottom temperatures of TES are depicted in Fig. 6(c). We observe that chilled water temperatures in the MPC strategy have larger variations than those in the baseline strategy. We also observe that for the MPC strategy the top and bottom layers of TES reach the almost same temperature at the end of day. But for the baseline strategy the bottom layer has much lower temperature than  For both strategies, the hourly cooling amount provided by individual chillers and TES are illustrated in Fig. 7, where the negative cooling amount means that TES is in a charging mode. We observe that in both strategies chiller 5 is never operated, and chillers 1, 2, and 7 are operated in most of time. This makes sense since chiller 5 has the lowest COP while chillers 1, 2, and 7 have almost highest COP in the range of interests for the chilled water return temperature as shown in Fig. 3. However, two strategies are different in the way to charge/discharge TES. We observe that TES charges more cooling amount in the baseline strategy than in MPC during the off-peak hours. However, TES discharges less cooling in the baseline strategy than in MPC during the on-peak hours. The overcharge in off-peak hours and insufficiently discharge in on-peak hours are the main causes of energy inefficiency under the baseline strategy. The energy and economical inefficiency can also be observed from electricity consumption and cost plot in Fig. 5(b) and (c). We summarize the total electricity consumption and cost in Table III. From it, the MPC-based strategy performs better than the baseline by around 9.70%, in terms of total electricity consumption, and by around 10.84%, in terms of total electricity cost.

1) Horizon Length:
The horizon length has an important impact on the performance of solutions and the computation time of algorithms. For the case study presented in this paper, we choose the horizon length as 24 h for the illustration of simulation results. We have also tested different horizon values for a comparison purpose. Generally speaking, by choosing a longer horizon length (e.g., multiplies of 24 h), we obtain optimal solutions that take into consideration of the forecast information and system dynamics in a long time range. But the computation time for obtaining the solution increases a great deal as the length of horizon increases. On the other hand, by choosing a shorter horizon length, we obtain the solutions with much less computation time. But the solutions may be only optimal for a short time range and the performance of the MPC strategy can be largely degraded. In practice, one can try a couple of horizon lengths and choose the one that satisfies the performance and computation time requirements.
2) Uncertainty Analysis: The simulation results presented here are obtained under the assumption that the parameters are accurate and chiller system models can perfectly represent the real system. However, in reality, there usually exist uncertainties in parameters and models. In this section, a simple Monte Carlo uncertainty analysis is carried out to demonstrate how the modeling uncertainty would affect the potential energy saving of proposed MPC strategies [27]. Here, we consider the uncertainties for some major parameters of the chiller system, namely, (chilled water supply temperature), and , , , and , and (fitted parameters of the TES model). We assume the considered parameters all follow the normal distribu-tions with means equal to their nominal values and standard deviations equal to 5% of their nominal values. A sample matrix is generated based on random sampling of the considered parameters, where denotes the number of parameters and denotes the sample size. Here we have 7 parameters and we take as the sample size. Then we repeatedly apply baseline and MPC strategies using the randomly generated parameter samples. The probabilistic distributions of energy savings (electricity consumption or cost) can be obtained via the numerical approximations. For the case study considered here, the standard deviations of energy savings are less than 2.3% compared with around 10% energy savings for nominal values. The results show that the fluctuation of energy savings caused by parameter uncertainties is relatively small. Thus, the proposed MPC strategy is not very sensitive to the modeling parameters specified here.

VII. CONCLUSION
This work formulates an MPC problem to obtain optimal set-points for supervisory control of a central chiller plant with thermal energy storage (TES). It proposes an efficient heuristic algorithm to solve the MPC problem based on dynamic programming (DP) and mixed integer linear programming (MILP). The optimal TES operation profile is obtained for a fixed horizon by solving a DP problem, and then optimal chiller operations are obtained at each time step by solving an MILP problem with a fixed TES operation profile. Simulation results show that both electricity consumption and cost can be significantly reduced for the campus central plant at the University of California, Irvine. Future research includes more benchmark studies with other control methods [28], [29] and using the optimization methods successfully used in other fields to solve this problem [30]- [49].

APPENDIX CAMPUS COOLING LOAD PREDICTION
For the MPC problem, we need to know the predicted campus cooling load for the next hours, where is the prediction horizon. In this paper, we employ a linear ARX model to predict the campus cooling loads [50]. Mathematically, at each time step , we predict the campus cooling loads for the next time steps. For , we have (17) where , , and denote the predicted cooling loads, ambient temperature, and ambient humidity ratio for time step at the time step , respectively. The parameters , , and can be obtained using historically measured campus cooling load and weather data. Note that here we only consider the ambient temperature and humidity ratio as the inputs to the ARX model, which have most significant correlation to the campus cooling loads. To further improve the prediction accuracy, one could take more variables as the inputs to the ARX model, e.g., ambient solar radiation and campus occupancy information.
The past values of campus cooling loads are available on an hourly basis from the campus plant management system. At time step , we assume the knowledge of the actual values of up to time step : , for all . The past information of ambient temperature and humidity ratio are available from on-site measurement. The predicted weather information (ambient temperature and ambient humidity ratio for next few days) can be obtained from National Digital Database via NOAA's National Weather Service. At time step , we assume the knowledge of the actual values and up to time step , and the predicted values of and for future time steps: .