Optimal scheduling of chiller plant with thermal energy storage using mixed integer linear programming

In this paper, we consider 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 the TES and supply chilled water to the campus. A bilinear model is established to describe the system dynamics. A model predictive control (MPC) problem is formulated to obtain optimal set-points to satisfy the campus cooling demands and minimize daily electricity costs. At each time step, the MPC problem is represented as a large-scale mixed integer nonlinear programming (MINLP) problem. We propose a heuristic algorithm to search for suboptimal solutions to the MINLP problem based on mixed integer linear programming (MILP), where the system dynamics is linearized along the simulated trajectories of the system. Simulation results show good performance and computational tractability of the proposed algorithm.


I. INTRODUCTION
For a campus with a large number of buildings, a central chiller plant is commonly used to serve its cooling loads.Depending on the size and number of buildings, the campus cooling loads can be huge and vary in a large range.A bank of multiple chillers is usually installed at the central plant to meet various cooling demands with high operational flexibility.These chillers can have different performance characteristics and they are often operated in ON/OFF modes [1].Plant operator needs to decide a scheduling strategy to mix and match various chiller combinations to meet the campus cooling demands, while considering the need to reduce the energy consumptions and operating costs.
Thermal Energy Storage (TES) can be used to shape the campus cooling loads.In the simplest case, a TES is a large tank for the storage 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.TES is charged by chiller bank at night when electricity rates and cooling demands are both low.During daytime hours, TES discharges the stored chilled water for campus cooling, thus reducing peak loads of chillers.
Well-designed control system can help realizing the full potential of both energy and cost savings for campus central Kun   plant.The goal of control here is to minimize energy consumption or operation cost while satisfying campus cooling demands and system operational constraints.The focus of this paper is on the design of high-level supervisory control, using model predictive control (MPC) scheme, to determine the optimal set-points for scheduling the chiller bank with TES.Several studies have identified promising saving potentials through the MPC scheme for the central chiller plant [2]- [7].Most of these studies assumed that the mass flow rates of the chilled water can be continuously adjusted for chillers.However, in many practical retro-fit cases, chiller operation decisions are ON/OFF scheduling sequences.For such problems the optimization/control framework proposed in above works would not directly apply.To generate optimal ON/OFF sequences for chillers, a mixed integer nonlinear programming (MINLP) problem needs to be solved at each time step for the MPC scheme.This MINLP problem involves considerable nonlinearity and a large number of integer variables, therefore finding a global optimal solution is nontrivial and computationally expensive.
In this paper, we propose a MILP-based heuristic algorithm to search for the suboptimal solutions of MINLP problem.The main idea is to formulate a MILP problem to approximate the MINLP problem by linearizing the system dynamics along the nominal trajectory of the system.Here the nominal trajectory 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 the MILP problem can be solved very efficiently with guaranteed performance.Moreover, the model mismatch due to the linearization can be compensated by the moving horizon scheme of MPC.We demonstrate the proposed algorithm through a case study for the campus central plant at University of California, Irvine.

II. SYSTEM MODELING
The schematic diagram of a campus central plant is depicted in Fig. 1, with all pumps omitted.In this section, we describe simplified system models for the real-time optimization purpose.

A. Nomenclature
The notations used in this paper are summarized below: i Index of chiller T amb Ambient temperature (K) c pw Specific heat capacity of water (kJ/kg

B. Chiller bank model
The chiller bank is composed of n c 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 a fixed supply temperature; When a chiller is OFF, it is shut down and supplies no chilled water.For i = 1, . . ., n c , the operating set-point of the ith chiller is given by a binary decision variable The amount of cooling provided by the ith chiller is The total electricity consumed by a chiller depends on several factors, e.g., supply water mass flow rates ṁchw,i , generated cooling amount Qchw,i , circulated water temperature differentiation (T chwr,i − T chws,i ).Note that ṁchw,i is a fixed constant when chiller is power on, and Qchw,i is a bilinear function of ṁchw,i and (T chwr,i − T chws,i ).Thus the chiller electricity consumption mainly depends on temperature differentiation (T chwr,i − T chws,i ).Following [6], the electricity consumption W chw,i of the ith chiller is simply modeled as: where the coefficients a i and b i can be obtained by regression analysis of the measured chiller data.The Coefficient of Performance (COP) is employed as a metric to quantify the energy efficiency of a chiller The total chilled water mass flow rate supplied by the chiller bank is The total power consumed by the chiller bank is The supply temperature and return temperature of the chiller bank are denoted by T chws and T chwr , respectively.All chillers are assumed to have the same supply temperature T chws and the same return temperature T chwr , i.e., for i = 1, . . ., n c ,

C. TES model
We employ a stratified two layer TES model developed in [6], where the state of TES is captured by the top layer water temperature T a and the bottom layer water temperature T b .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 is operated in two modes: in charging mode, chilled water from chiller bank is supplied to both the TES and the campus; in discharging mode, chilled water from the TES and the chiller is supplied to the campus together.The TES operating mode is denoted by a binary variable: In charging mode (i.e., σ = 1), the dynamics of TES is In discharging mode (i.e., σ = 0), the dynamics of TES is Parameters f ac , f bc , f ad , f bd are time constant multipliers, and U c ,U d 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: Here ṁcam , T cams , and T camr are taken as system variables, while Qcam is taken as a given exogenous input to the system.Note that, Qcam can be predicted based on the historical load data and the weather forecast of ambient temperature, humidity, and solar radiation [3].

E. Model summary
By collecting system dynamic equations ( 1) to ( 11) and discretizing the system with the sampling time ∆t, we represent the system dynamics in the state-space form: where the state variable of the system is denoted by the intermediate variable of the system is denoted by the control variable of the system is denoted by and the input variable to the system is denoted by The equation (12a) corresponds to TES dynamic equations (9d), (9e), (10d), and (10e); The equation (12b) corresponds to TES and campus balance equations (9a) to (9c), (10a) to (10c), and (11); The equation (12c) corresponds to chiller balance equations ( 2) to (7).Besides system dynamics, we also impose various constraints for system variables, i.e., x ∈ X, z ∈ Z, and u ∈ U for some bounded sets X, Z, and U.These constraints are due to equipment limitations, or desirable operating conditions of main components.

III. MODEL PREDICTIVE CONTROL PROBLEM
The objective of the optimization is to obtain a real-time optimal scheduling strategy for chiller bank and TES to minimize the total electricity costs, and satisfy the campus cooling demands and system constraints.We formulate a Model Predictive Control (MPC) problem to search for the optimal set points of the scheduling strategy.At each time step k, we are interested in solving the following finitehorizon optimization problem regarding the state and control actions at time steps {k, . . ., k + N − 1}, where N is the prediction horizon of the MPC problem: min (13) where x(l|k) denotes the state variable for the time step l predicted at time step k (similar notations also apply for other variables); the sets denote the predicted state, intermediate, and control sequences over the predicted horizon; electricity rate c(l|k) is assumed to be known for all time and z 1 = W chw denotes the electricity power consumed by the chiller bank; the cooling load d(l|k) is a given input that can be obtained from weather forecast; x(k) denotes the actual state of the TES at time k; and X f denotes the terminal constraint set for the TES state variables, which is chosen to ensure stability and feasibility of the MPC problem [3], [8].
By assuming that the maximum campus cooling load never exceeds the maximum capacity of the chiller bank, we know that there always exist feasible solutions to the problem (13).Among all feasible solutions, let (X * (k), Z * (k), U * (k)) denote the optimal solution to the problem (13).At time step k, only the first element of U * (k) is implemented to the system, i.e., u(k) = u * (k|k).At time step k + 1, the horizon of interests moves one-step forward to {k + 1, . . ., k + N}, and the optimization problem (13) is solved again by changing index from k to k + 1.The initial point of the predicted state sequence is taken as x(k + 1|k + 1) = x(k + 1), where x(k + 1) denotes the newly measured state at time step k + 1.This procedure is repeatedly implemented by moving forward the horizon step-by-step, which yields a moving horizon control strategy or model predictive control strategy.
For the MPC optimization problem (13), at each timestep k, there are a total of N * (4n c + 8) decision variables to be computed, of which there are N * (n c + 1) binary variables, and N * (3n c + 7) real variables.We also note that the constraints due to the system dynamics ( 12) have non-convex terms due the bilinear system dynamics.Thus, the optimization problem (13) is a non-convex mixed-integer nonlinear programming (MINLP) problem, which is quite difficult to solve in both theory and practice [9].There are many software packages that can solve non-convex MINLP problem to proven optimality, such as BARON, α − BB, LINDO − Global, and Couenne [9].But these software packages require a lot of computing resources and a large amount of computation time to achieve a certain global optimality.These software packages can therefore be used off-line to compute the optimal set points, but they are not suitable for real-time implementation of the MPC scheme.

IV. MILP-BASED HEURISTIC ALGORITHM
In this paper, we propose a heuristic algorithm to solve the non-convex MINLP problem (13) by further exploring its special structure.Our goal is to obtain sub-optimal solutions to (13) in a time efficient way with performance guarantee.

A. Deciding TES operation mode profile
TES is the major component in the system to shift the campus cooling load from the peak hours to the off-peak hours.Intuitively, we should charge the TES during off-peak hours when electricity rate is low, and discharge the TES during the peak hours or partial-peak when electricity rate is high.In order to reduce the complexity of optimization [3], [6], we also assume a fixed operation profile, denoted by σ f ix , for charging/discharging the TES.

B. Linearizing system dynamics
It appears that the major difficulties of solving a nonconvex MINLP problem are due to the non-convexities of the problem, which prevents the implementation of efficient algorithms [9].For problem (13), non-convexities mainly come from bilinear terms in chiller and TES dynamics.Let x(k), z(k), and ū(k) denote the nominal values of the state and control at the kth time step.We approximate the bilinear dynamic equation (12a) by a set of linear dynamic equations where the matrices A 1 (k), A 2 (k), B(k), E 0 (k) are obtained through linearization of (12a) around nominal values x(k), z(k), ū(k).Similarly, we approximate bilinear balance equation (12b) by a set of linear balance equations To remove the bilinear terms introduced by the product of a binary valued variable and a real valued function, we employ an equivalence modeling trick described below [10]: (P1): Let ξ (ω) ∈ R be a real valued scalar function defined on a bounded set Ω such that m ≤ ξ (ω) ≤ M for any ω ∈ Ω, κ ∈ {0, 1} be a binary valued scalar variable, µ ∈ R be a real valued scalar variable.For any ω ∈ Ω, we have the following equivalence: We apply Property (P1) to convert the bilinear equation (12c) equivalently to a set of linear inequalities

C. Formulating an MILP problem
In this section, we define a simplified MPC optimization problem with the linearized model described in Section IV-B.We fix the TES operation mode profile as σ f ix .At each time step k, we choose a nominal trajectory of state and control The linearized dynamics for the predicted horizon {k, . . ., k + N − 1} is obtained along the nominal trajectory of state and control.Then we consider the following finite-horizon optimization problem as an approximation to the problem (13): min The optimization problem ( 14) is a mixed integer linear programming (MILP) problem, with N * n c binary variables and N * (3n c + 7) real variables.Due to the convexity of the problem, the MILP problem is much easier to solve than the non-convex MINLP problem in general.There exist many software packages that can efficiently solve an MILP problem with certain performance bound, e.g., SCIP, GLPK, lpsolve, and CBC [11], [12].Here the performance bound provides a quantitative estimation of the goodness of suboptimal solutions compared to the global optimal solution.

V. CASE STUDY
In this section, we present a case study to demonstrate the model predictive control of the central chiller plant in summer time at the University of California, Irvine (UCI).

A. System setup
System diagram of the UCI central plant is depicted in Fig. 1.A bank of n c = 7 chillers supplies chilled water to satisfy the campus cooling demand.The supply water  The chilled water supplied by the chiller bank can also be stored in a water tank TES for shifting the campus cooling loads.The cross-layer area of the TES is A = 564.95m 2 with a height of h = 30.48m.At the initial time, the TES is assumed to be fully discharged with the uniform water temperature of 60 • F. kkThe calibrated time-constant multipliers of the TES model ( 9) and ( 10) are f ac = 1.740, f bc = 2.222, f ad = 4.138, and f bd = 1.599.The calibrated inter-layer overall heat transfer coefficients are U c = 0.0283kW/(m 2 • K) and U d = 0.0197kW/(m 2 • K).

B. Two scheduling strategies
In this section, we consider two strategies for scheduling chiller bank and the TES.One is a greedy-search-based heuristic approach which serves as a baseline strategy, and the other is the proposed MILP-based MPC strategy.In both strategies, the TES is operated with a fixed operation profile δ f ix (k) =1 for 0 ≤ k ≤ 9 and δ f ix (k) = 0 for 10 ≤ k ≤ 23.
1) Baseline strategy: This strategy generates sub-optimal chiller and TES operation sequences based on a greedy search algorithm.All possible 2 7 chiller ON/OFF modes are first ordered according to their efficiency.Then for a given TES charging and 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 mostly reduce the Time-Of-Use (TOU) based cost 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 [6].
2) MILP-based MPC strategy: This strategy implements the MPC scheme by considering the MILP formulation (14).The MILP-based heuristical algorithm described in Section IV is employed to search for the optimal set points of the system.The sampling time is taken as ∆t = 1 hours to discretize the system model.The prediction horizon of the MPC scheme is taken as N = 24 hours.At each hour, the software package SCIP is employed to solve the MILP problem (14) with performance gap ε = 5%.For the problem considered in this paper, it usually takes less than 10 minutes to achieve the performance gap.The initial nominal trajectory is obtained by simulating the system with the baseline strategy.

C. Comparison of two strategies
The performance of two strategies are compared through MATLAB simulations for a specific day in June, 2010.The predicted campus cooling load profile is depicted in Fig. 2 (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 control set points computed from the two strategies are applied to the simulation model of the UCI central plant.For both strategies, the hourly cooling amount provided by individual chillers and the TES are illustrated in Fig. 4, where the negative cooling amount means that the TES is in charging mode.We observe that in both strategies chiller 5 is never operated, and chillers 1, 2, 7 are operated in most of time.This makes sense since chiller 5 has the lowest COP while chillers 1, 2, 7 have almost highest COP in the range of interests for the chilled water return temperature (see Fig. 3 for COP plot).But two strategies are different in the way to charge or discharge the TES.We observe that TES charges more cooling amount in baseline strategy than in MPC strategy during off-peak hours.TES discharges less cooling in baseline strategy than in MPC strategy during onpeak hours.The overcharge in off-peak hours and insufficient discharge in on-peak hours are the main causes of energy inefficiency using baseline strategy.The hourly electricity consumption and hourly electricity cost of two strategies are depicted in Fig. 2 (b) and (c), respectively.We observe that MPC strategy outperforms the baseline in almost all time.

VI. CONCLUSIONS
In this paper, we formulated a model predictive control (MPC) approach for the ON/OFF scheduling problem of a central chiller plant with thermal energy storage.We proposed an efficient heuristic algorithm to search for a suboptimal scheduling strategies based on mixed integer linear programming.Simulation results showed that electricity costs are largely reduced using MILP-based MPC strategy.

Fig. 1 .
Fig. 1.A systematic diagram of a typical campus central plant.The chilled water is supplied by a bank of n c chillers, and circulates through campus and Thermal Energy Storage (TES).More details can be found in Section II.

Fig. 4 .
Fig. 4. Two strategies scheduling results: hourly cooling amount generated by individual chillers and TES.