General methodology combining engineering optimization of primary HVAC & R plants with decision analysis methods--Part I: deterministic analysis.

This paper is the first of a two-part sequence that proposes a general methodology for dynamic scheduling and optimal control of complex primary HVAC & R plants, which combines engineering analyses within a practical decision analysis framework by modeling risk attitudes of the operator. The methodology involves a computationally efficient, deterministic engineering optimization phase for scheduling and controlling primary systems over the planning horizon, followed by a systematic and comprehensive stochastic sensitivity and decision analysis phase, where various sources of uncertainties are evaluated along with alternative non-optimal but risk-averse operating strategies. This paper describes the deterministic component of the analysis methodology, which essentially involves the development of response surface models for different combinations of system configurations to be used for static optimization and then using them in conjunction with the modified Dijkstra's algorithm for dynamic scheduling and optimal control under different operating conditions and pricing signals. The proposed methodology is illustrated for a semi-real hybrid cooling plant operated under two different pricing schemes: real-time pricing and time-of-use with electricity demand. We feel that the general methodology framework proposed sacrifices very little in accuracy while being much more efficient computationally than the more complicated optimization methods proposed in the general literature. Moreover, this approach is suitable for online implementation, and it is also general enough to be relevant to other energy systems.

INTRODUCTION

A literature review on supervisory optimal control applied to the operation of complex HVAC & R systems including cooling plants and BCHP plants has been done by Jiang (2005). There are several studies that have adopted mixed integer linear programming (MILP) techniques to this general problem (for example, Dotzauer [1997, 2003] and Yokoyama et al. [2002]). However, the structural design problem has often been treated by considering only a single-period operation (Papoulias and Grossmann 1983) or a multi-period one with a small number of periods (Horii et. al. 1987). Some approaches based on meta-heuristics, such as simulated annealing (SA) and genetic algorithms (GA), have also been proposed (for example, Sakamoto et al. [1999], Curti et al. [2000], and Yin and Wong [2001]). However, these approaches are said to have limitations in the determination of values of search parameters, the judgment of optimality, and the requirement of extensive computation times (Yokoyama et al. 2002).

Several levels of optimal control schemes have been proposed for existing cooling plant operation. These can be grouped broadly as follows:

1. Cookbook solutions, which are simple rules and guidelines for operators to follow (Hydeman 2002).

2. Heuristic control schemes, widely used in the current building control profession, that are developed based on local optimization, system model simplification, estimation, and experience. The 2003 ASHRAE Handbook--HVAC Applications (ASHRAE 2003) describes in detail such control heuristics for operating HVAC systems and components. Further, control heuristics could be used as a starting point in an optimization scheme. In addition, a heuristic type of suboptimal control is often desirable for online implementation purposes.

3. Rigorous optimization algorithms that follow the strict definition of optimal control by proposing optimization algorithms to minimize the objective function (which is often the cost). The approach proposed in this paper falls in this category.

Because of the variety of energy sources used in complex HVAC & R systems, the interdependency between sources, and the variation of technical and economic conditions with time, e.g., change of load, deterioration of equipment, change of fuel and electricity prices, etc., the planning of plant day-to-day operation and evaluation of alternative performance options is not simple. Much of the difficulty is mainly due to the following reasons:

* The objective functions and models are usually nonlinear functions that may contain both discrete (for example, equipment on/off status) and continuous variables--locating the global optimum is not guaranteed. Further, one may have to deal with multiple objectives, which make the problem even more complicated.

* The possible number of independent or decision variables for the problem is large, with a large set of diverse constraints, therefore presenting the engineer with the difficult, if not impossible, task of determining the best operating strategy. Further, if the problem is a multi-period dynamic problem (i.e., involving several stages), with the number of binary scheduling variables increasing with the number of periods, the conventional solution algorithm, which combines the branch and bound method with the simplex method, may require computation times that are not practical (Yokoyama et al. 2002).

Despite recent advances in computer power and the development of better optimization algorithms, only a few are used in industry. What is more remarkable is that most complex HVAC & R plants are still scheduled by humans in a heuristic manner without the aid of computer supporting tools. One possible reason for this often voiced by professionals is the lack of consideration of how to combine pure engineering solutions with individual risk attitudes of how system operators weigh risk over predicted outcome. It is, in essence, this aspect that is addressed by this research.

OBJECTIVE AND SCOPE

The primary objective of this paper is to propose a general and computationally efficient methodology for minimizing the operating costs, including both energy costs and demand costs, of complex HVAC & R plants over the planning horizon, which is taken as 12 hours. The operating cost would include electricity usage cost, gas usage cost, and equipment start-up cost. A true optimization would require the simultaneous optimization of all cost components under the pre-specified thermal load and well-defined performance characteristics and maintenance costs of equipment. An even finer level of analysis would be to consider the reliability associated with different equipment, since the large equipment could be of different vintage and level of degradation.

Utility costs differ with different pricing signals, resulting in different formulations of the optimization cost function, which, in turn, may require different optimization techniques. Only two cases are considered: (a) real-time pricing, which has no demand charge involving energy (electricity and gas) over a certain time period along with start-up cost, and (b) TOU (time of use) with demand, which is more complex since the cost function includes gas and electricity cost (energy cost + demand cost) and start-up cost. The optimization function should explicitly consider start-up cost caused both by additional energy consumption and increased demand. In order to minimize the demand charge, equipment must be operated so that situations that cause large spikes in power consumption (due to having to accelerate components such as fans and motors up to their design speeds) do not occur during periods of peak power.

EQUIPMENT MODELS

Selecting a performance model is an important and essential first step in optimizing the operation of any engineering system. The main components in a cooling plant include chillers, cooling towers, fans, and pumps.

Gordon-Ng (GN) Chiller Model

The semi-empirical GN chiller model (Gordon and Ng 2000) predicts the dependence of chiller COP (defined as the ratio of chiller thermal cooling capacity divided by the electrical power consumed by the compressor) with certain key (and easily measurable) parameters such as the fluid (water or refrigerant) return temperature from the condenser, fluid temperature leaving the evaporator (or the chilled-water supply temperature to the building), and the thermal cooling capacity of the evaporator. A detailed evaluation consisting of over 50 chillers of all types (one-stage, two-stage centrifugal with inlet-guide and VSD, screw, scroll, and reciprocating) and sizes has been conducted by Jiang and Reddy (2003). It was found that the fundamental GN formulation for all types of vapor compression chillers is excellent in terms of its predictive ability.

The following equation is the GN fundamental model for vapor compression chillers:

([1/COP] + 1)[[T.sub.cho]/[T.sub.cdi]] - 1 = [a.sub.1][[T.sub.cho]/[Q.sub.ch]] + [a.sub.2][([T.sub.cdi] - [T.sub.cho])/[[T.sub.cdi][Q.sub.ch]]] + [a.sub.3][[(1/COP + 1)[Q.sub.ch]]/[T.sub.cdi]] (1)

where [a.sub.1], [a.sub.2], and [a.sub.3] are regression coefficients, [T.sub.cho] is chilled-water outlet temperature (K), [T.sub.cdi] is condenser water inlet temperature (K), and [Q.sub.ch] is chiller load (kW).

GN models for single-stage absorption chillers are also valid and have been shown by Jiang and Reddy (2003) to be accurate (with coefficient of variation [CV] about 6%-8%) for steam and hot water two-stage absorption systems.

([[[T.sub.gni] - [T.sub.cdi]]/[[T.sub.gni]COP]] - [[[T.sub.gni] - [T.sub.cho]]/[T.sub.cho]])[Q.sub.ch] = [b.sub.0] + [b.sub.1][[T.sub.cdi]/[T.sub.gni]] (2)

where [b.sub.0] and [b.sub.1] are regression coefficients and [T.sub.gni] is generator inlet temperature (K).

Effectiveness Cooling Tower Model

The effectiveness-NTU model concept, originally proposed for sensible heat exchangers, was modified by Braun (1988) and Braun et al. (1989) to model performance of cooling towers by utilizing the assumption of a linearized air saturation enthalpy. The following general correlation for NTU in terms of the flow rates is used with estimates of the coefficients c and n identified from measurements at different air flow rates [dot.m.sub.a], (with water flow rate [dot.m.sub.w] being, in most cases, constant):

NTU = c([dot.m.sub.w]/[dot.m.sub.a])[.sup.1 + n] (3)

Using the standard expression for effectiveness of a counterflow cooling tower, outlet water temperature is determined from an energy balance on the cooling tower,

[T.sub.w,o] = [T.sub.ref] + [[[dot.m.sub.w,j]([T.sub.w,i] - [T.sub.ref])[C.sub.pw] - [dot.m.sub.a]([h.sub.a,o] - [h.sub.a,i])]/[dot.m.sub.w,o][C.sub.pw]]. (4)

Fan and Pump Model

A third-order polynomial model representing the relationship between fan (or pump) power [P.sub.fan] and [m.sub.a] used in the HVAC 2 Toolkit (Brandemuehl 1993) is adopted in this study. Performance at off-rated conditions is calculated from the rated performance using the part-load ratio model:

[P.sub.fan](t) = FMP{[e.sub.0] + [e.sub.1][PLR(t)] + [e.sub.2][PLR(t)][.sup.2] + [e.sub.3][PLR(t)][.sup.3]} (5)

where

FMP = fan motor power at rated condition, kW

[e.sub.0], [e.sub.1], [e.sub.2], and [e.sub.3] = fan performance coefficients

PLR(t) = fan part-load ratio defined as the ratio of total air flow to fan capacity

In constant-flow systems, [P.sub.pump] is essentially constant. However, for variable-flow pumps, [P.sub.pump] is a function of the building loads or, more specifically, of the fluid flow rate [m.sub.w]. Phelan et al. (1997) studied the predictive ability of linear and quadratic models between [P.sub.pump] and [m.sub.w] and concluded that quadratic models are superior to linear models. A curve similar to that introduced for fan power (Equation 5) can also be used to model the relationship between [P.sub.pump] and [m.sub.w].

GENERAL METHODOLOGY FOR OPTIMAL OPERATION

Formulation of the Objective Function

The objective function is the utility cost function, which is different for different pricing signals.

A real-time pricing case has no demand charge and is the sum of the operating cost of equipment under steady-state operation and the cost of additional energy use during equipment start-up. The objective function can be written as

[F.sub.RTP] = Min{[T.summation over (t = 1)][T.summation over (k = 1)][[R.sub.t,k][P.sub.t,k] + (1 - [u.sub.t - 1,k])[u.sub.t,k]S[C.sub.t,k]]}, (6)

subject to

[h.sub.i]([x.sub.1], [x.sub.2],..., [x.sub.K]) = [b.sub.i], i = 1, 2,..., m, (7)

[g.sub.j]([x.sub.1], [x.sub.2],..., [x.sub.K]) < [c.sub.j], j = 1, 2,..., n, and (8)

[x.sub.k.sup.l] [less than or equal to] [x.sub.k] [less than or equal to] [x.sub.k.sup.u], k = 1, 2,..., K, (9)

where t and k are indices for time intervals and equipment, respectively, and T and K denote the total number of time intervals and equipment, respectively. Equations 7 through 9 represent equality constraints, inequality constraints, and boundary or range constraints, respectively; [x.sub.1], [x.sub.2],..., [x.sub.K] are control variables and P and SC are functions of control variable x.

The TOU with electricity demand case is more complex since the cost function includes steady-state costs of gas and electricity (energy cost + demand cost) as well as start-up cost. The objective function, subjected to similar equality and inequality constraints as Equations 7 through 9, is expressed as

[F.sub.TOU] = Min{[T.summation over (t = 1)][K.summation over (k = 1)][([R.sub.k][P.sub.k]) + (1 - [u.sub.t - 1,k])[u.sub.t,k]S[C.sub.t,k]] + [R.sub.de][T.max.[t = 1]][K.summation over (k = 1)][P.sub.ele,t,k]}. (10)

In Equation 10, the first term on the right-hand side of the equation is the steady-state energy cost, the second term is the additional energy cost due to equipment start-up, and the last term is the demand cost. Note that this equation applies to the demand-setting day of the month when the additional demand charges are incurred.

Estimation of Expected Energy Consumption Based on Static Optimization Case

The optimization problem for operating a plant that meets a pre-specified thermal load consists of a two-level hierarchical structure because of the following two different types of decision variables. The first, or higher-level unit commitment problem, involves discrete control variables that are not continuously adjustable but are discrete, such as the number of operating chillers, cooling tower cells, condenser water pumps, chilled-water pumps, and the relative speeds for multi-speed fans. The second, or lower-level economic dispatch problem, involves control variables that need to be controlled continuously. Independent continuous control variables might include the chilled-water temperature setpoints, relative water flow rates to the chillers, cooling tower cells, and the speeds for variable-speed fans or pumps and so on. Therefore, the optimization of plant operation should consider both selection of which equipment to use and how to operate it.

Description of Static Optimization Case. The static optimization case involves optimizing the operating cost for each time step, i.e., each hour. The cost components include only steady-state hourly energy costs for electricity and gas. So the quantity to be minimized, [F.sub.ss], is the total cost of energy consumption, summed over all components that are operating. The energy consumption [P.sub.k] for each of the k components is a function of the component's characteristics and is dependent on the controlled variables as given by a set of output equations. Energy usage for each component has an associated cost rate [R.sub.k], which can be time-dependent (e.g., time-of-day electrical rates). A set of optimization variables, [x.sub.k], is sought, which minimizes the objective function [F.sub.s] over an hour with respect to the independent continuous and discrete control variables.

[F.sub.s] = Min{[K.summation over (k = 1)][R.sub.k][P.sub.k]} (11)

Generation of Plant Operating Modes. There are a number of feasible operating modes (i.e., potential ways of combining various equipment to meet the demand at current weather conditions). Note that a feasible operating mode is one that satisfies the constraints, while the whole set of operating modes consists of all the possible combinations of controllable discrete variable states for each individual plant component without imposing any constraints. The approach adopted here is to generate tables (in the form of matrices) for each type of equipment, where rows represent different state combinations of that type of equipment. Heuristic constraints are applied to these tables to generate feasible combinations for each type of equipment automatically, after which the feasible combinations for the whole plant can be generated. The algorithm to automatically generate the operating modes for the entire plant is illustrated in the following example.

Consider a small cooling plant that has two parallel chillers of the same type and size, two one-cell cooling towers of the same type, and two chilled-water pumps and two condenser water pumps of the same type. Here, fans of the cooling towers are assumed to be two-speed, and all the pumps are taken to be fixed-speed. Each chilled-water pump is dedicated to a chiller and each condenser water pump is dedicated to a cooling tower.

First, we generate all the possible operating modes for the two chillers that are represented by the matrix ComC[H.sub.t]. Since the status of a chiller can only be OFF/ON, this can be represented by a binary value (0/1). So the maximum number of chiller combinations is [2.sup.N] (where N is the number of chillers). For this simple case, N = 2 and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The first row represents the case when both chillers are off, while the last row represents the case when both chillers are on. Some combinations can be eliminated based on physical constraints. For example, either row 2 or 3 can be deleted because the chillers are of the same type and size. Subsequently, the additional constraint that the total chiller load should be larger than the building load can be used to identify feasible operating modes. For example, if the building load is such that one chiller alone can satisfy it, the feasible operating modes would include operating one chiller as well as both chillers. The corresponding matrix is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Similarly, cooling tower operating modes can be generated and represented by the matrix ComC[T.sub.t]. The operating modes for fans are not just ON and OFF; they also have low-speed and high-speed status. Hence, fan status needs to be represented by (0/1/2). Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The number of possible fan combinations for one simple cooling system is [3.sup.N]. Again, the physical constraint that all fans are similar can be used to reduce the cooling towers' operating modes to yield the feasible matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Chilled-water pumps' operating modes represented by ComPME and condenser water pumps denoted by ComPMC need not be considered separately because they are both dedicated to chillers and cooling towers, respectively. Finally, the feasible operating modes for the whole cooling plant COM can be generated by combining each of the feasible operating modes for individual components, i.e., by combinatorial operation of rows of the matrices ComCH and ComCT,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where ComCH(n) (n = 1, 2,...) represents the nth row of matrix ComCH. Thus, in expanded form:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Determination of Optimal Control Setpoints of Equipment. Next, the sequential quadratic programming (SQP) algorithm (Fletcher 2001) is used to determine the optimal operating setpoints (for example, chiller load fraction, cooling tower fan speed) that result in minimal energy cost for each mode. By performing this analysis for all feasible equipment combinations, the equipment operational control settings that result in minimal energy cost for each equipment combination are determined, from which the best combination of equipment to meet the current load can be selected.

Response Surface Analysis Based on Optimized Operation for Each Operating Mode. Although it is possible to use the procedure described above to determine the optimal control strategy to meet a specified cooling load at a given time, it may be impractical from a computational viewpoint to use it to predict the optimal path over a given time horizon (i.e., the dynamic case). Hence, there is a distinct advantage in simplifying the process of identifying the optimal control setting (as was done by Koeppel et al. [1995] for a double-effect water/Li-Br absorption chiller). Instead, we propose to use experimental design techniques in general, and a response surface model in particular, as a computationally efficient alternative strategy. One performs a number of simulation runs beforehand for different operating conditions (a key feature is the factorial design method to determine these operating conditions), identifies the optimal control operation for each condition, and regresses these optimal control strategies using a polynomial model. This model is then used as a replacement or proxy for the numerical simulation model, and all inferences related to optimization/uncertainty analysis, requiring several thousands of simulations for the original model, are derived from this fitted model. This would allow reasonable computational time for estimating expected energy consumption when building loads are met with various combinations of large equipment, (here, large equipment means the equipment that consumes much more energy than the rest and has a start-up penalty; in a cooling plant, the large equipment would be chillers). These estimates could be used to identify when it is desirable to turn large equipment on or off. The response surface models need to be generated only once for a specific plant and need to be updated only when equipment is replaced or added or when its performance degrades appreciably.

In most response surface model problems, the form of the relationship between the response and the independent variables is unknown. So an important step is to define a suitable approximation for the true functional relationship between response and the set of independent variables. Usually, a polynomial model in some region of the independent variables is employed. Quadratic models are usually sufficient for most industrial engineering applications, though higher orders may be necessary in certain instances. Of course, it is unlikely that a polynomial model will be a reasonable approximation of the true functional relationship over the entire space of the independent variables, but for a relatively small region, they usually work quite well (Montgomery 1991).

Least-squares estimation of the model coefficients is most effective if proper experimental designs are used to collect the data. Central composite design (CCD) is probably the most widely used experimental design for fitting a second-order response surface (NIST 2005). A CCD contains an imbedded factorial or fractional factorial design with center points that is augmented with a group of axial points that allow estimation of curvature. The factorial or "cube" portion and center points may serve as a preliminary stage where one can fit a first-order (linear) model but still provide evidence regarding the importance of a second-order contribution or curvature. Center point runs (which are essentially repeats of the center point) are included to provide a measure of process stability and inherent variability and to provide a check for curvature (NIST 2005). For a four-factor experiment design, the CCD typically generates 30 data points, including 16 factorial points and 8 axial points. For a more elaborate discussion, refer to any appropriate statistical text (for example, Montgomery [1991] or Jiang [2005]).

To apply the response surface method to our problem, the optimization algorithm described above is applied to the plant system over a wide range of conditions that are selected by CCD for each combination. From the detailed optimization procedures for the static case, we can identify four important forcing input variables: cooling loads ([Q.sub.ch]), wet-bulb temperature ([T.sub.wb]), electricity rate ([R.sub.ele]), and gas rate ([R.sub.gas]). The diurnal variation of the first two variables can be obtained from past measured data. Models for electricity use or gas use can be calculated following a relation such as

{[P.sub.ele], [P.sub.gas]} = f([Q.sub.ch], [T.sub.wb], [R.sub.ele], [R.sub.gas]). (12)

The optimal setpoints are determined for each combination under different conditions of [Q.sub.ch], [T.sub.wb], [R.sub.ele], and [R.sub.gas]. A general polynomial regression model is fitted with parsimony being achieved by using stepwise least-squares regression at a significance level of 0.05.

How well the model fits the data can be ascertained from two statistical measures--the coefficient of determination adjusted for degrees of freedom (Adj-[R.sup.2]) and the coefficient of variation of the root mean square error (CV-RMSE). This, however, is but the preliminary step (a necessary but not sufficient criterion) in evaluating the predictive ability of a model (this type of evaluation is referred to as internal predictive ability or simulation error). A model with high internal predictive ability may not necessarily be robust enough to guarantee accurate predictions under different sets of operating conditions. This is a major issue in black box models but less so in physical models (Reddy and Andersen 2002). The well-accepted approach to evaluate the external predictive ability of such models is to use a data set other than the one used for model identification (or model training) to determine the predictive accuracy of the model. Therefore, an additional data set is created by running the optimization algorithm introduced previously under conditions other than those used in the experimental design conditions. This extra data set can then be used to evaluate the prediction capability of the fitted response models.

Minimization of Plant Energy Cost over a Period

The expected energy consumption of each possible operating mode during the current time interval (i.e., for each hour) can be determined from the response surface regression model. The next step is to find the optimal operating strategy over a certain time horizon that consists of several hourly periods. Hence, current loads and future loads need to be considered simultaneously. When demand charges are to be included, the optimization problem is to schedule the available equipment so that it meets the load as efficiently as possible while avoiding start-up spikes in energy consumption. For each of the anticipated loads in the planning horizon, any proposed methodology should determine which combination of equipment can be used to meet the load and how the equipment should be operated to minimize the energy cost.

A general strategy proposed by Olson and others (Olson 1988; Olson and Liebman 1990; Olson et al. 1993, 1994), called dynamic chiller sequencing (DCS), is adapted and used in our study (Jiang 2005). The method called Dijkstra's algorithm (Dijkstra 1959) is used to solve the equipment scheduling problem based on the shortest path algorithm, which in this case is the lowest cost path algorithm. It is one of the most computationally efficient algorithms for solving single-source shortest path problems. For example, its time complexity is proportional to the square of the number of time intervals as compared to the cube for the well-known discrete dynamic programming algorithm. The algorithm basically involves three steps: (1) construct the cost network over the planning time horizon for different feasible equipment combinations, (2) search for the optimal scheduling strategy, and (3) determine the optimal setpoints for practical implementation based on the optimal scheduling strategy found from step 2. In this study, certain realistic assumptions have been made: (1) that the energy consumption as a result of chiller start-up is 20% greater than the steady-state condition (Olson 1988) and (2) that it lasts for 15 minutes for a vapor compression chiller and one hour for an absorption chiller. These values can, of course, be modified as necessary for the specific equipment. Furthermore, accurate values for start-up costs can be obtained given power consumption measurements during the equipment start-up period. The reader can refer to the published papers by Olson et al. (1993, 1994) for a complete description of how the Dijkstra algorithm can be applied to HVAC & R systems.

To summarize, the engineering deterministic methodology proposed in this paper involves performing equipment selection and load allocation as well as setpoint determination simultaneously through a combination of nonlinear optimization and a solution of a shortest path problem. Optimization of the scheduling strategy and control variables under different rate structures is greatly accelerated by using response surface models to determine the energy cost for different equipment configurations under different operating conditions and price signals.

CASE STUDY: HYBRID COOLING PLANTS

As an illustration, our dynamic methodology is applied to a hybrid cooling plant so as to determine its optimal operating strategies under different utility rate structures. Deterministic optimization results are discussed, and the effectiveness of the methodology is demonstrated. The effects of various sources of uncertainty are discussed in the companion paper (Jiang et al. 2007).

Description of Hybrid Cooling Plant

An actual cooling plant serves as a model for this case study. The cooling plant, shown in Figure 1, consists of three chillers of equal rated capacity (one direct-fired absorption chiller and two centrifugal chillers), three variable-speed cooling towers, and six fixed-speed pumps (three for the evaporator water loop and three for the condenser water loop). The three chillers of 2110 kW (600 ton) rated cooling capacity each are parallel to each other, and so are the cooling towers. Usually the size of the cooling tower used for an absorption chiller is larger than that for a centrifugal chiller of the same cooling capacity, so the larger cooling tower is dedicated to the absorption chiller. The pumps are dedicated to the cooling towers and chillers as shown in Figure 1. Table 1 summarizes the specifications of the various components of the system assumed in this case study, and Table 2 lists the coefficients used in the equipment models, data source, and regression goodness-of-fit results.

[FIGURE 1 OMITTED]

Classification of Variables

The forcing variables (or exogenous variables) are measurable quantities that cannot be controlled but that affect the component outputs and/or costs, such as the total building thermal load ([Q.sub.ch]) and ambient wet-bulb temperature ([T.sub.wb]). The next step is to identify both the control variables for which the optimization is to be performed and the forcing variables that affect the system performance over time. With variable-speed fans and dedicated pumps, the only significant discrete control variables are the number of operating chillers and cooling towers. The independent continuous control variables considered include the chilled-water supply temperatures ([T.sub.cho]) and the speeds for the variable-speed cooling tower fans. The optimal control variables change over time in response to the forcing variables. To optimize this model for a particular chiller/cooling tower combination, the decision variables (optimization variables) are the cooling load of each chiller and the water temperatures leaving each cooling tower. A classification of the relevant variables of the system is given in Table 3.

Formulation of Constraints and Range of Variations

Both equality and inequality constraints influence the optimization of the cooling plant system. Different constraints are chosen based on performance characteristics of equipment and energy balance relationships of energy flow for different systems. The most obvious constraint is for the cooling plant system to meet the required cooling load at any time. The simplest type of inequality constraint is to place bounds on control variables. For example, lower and upper limits are necessary for the chilled-water set temperature in order to avoid freezing in the evaporator and to provide adequate dehumidification for the zone. For the system being studied, the assumed constraints are listed in Table 4.

Response Surface Analysis Based on Optimization Results of Static Case

As described earlier, predictions of energy consumption when the total cooling load is met with various combinations of chillers are conveniently done using the response surface models obtained by regressing the simulation results of the static optimization under different system operating modes.

Experimental Design and Response Surface Model Fitting. For a hybrid cooling plant, the four important forcing input variables are cooling loads ([Q.sub.ch]), wet-bulb temperature ([T.sub.wb]), electricity rate ([R.sub.ele]), and gas rate ([R.sub.gas]). The feasible configurations (or groups) for running the plant during each hour to meet the building load can be determined. Table 5 shows all the feasible configurations (or groups) for this case study. For example, group G1 corresponds to the case when only chiller CH3 and cooling tower CT3 are operating.

Response Surface Model Fitting. The CCD is adopted to generate the data set to fit models for each group using response surface methods. (The interested reader can refer to Jiang [2005]) for details. Under each operating condition generated by CCD, the SQP optimization algorithm is applied to the hybrid cooling plant system so as to find the optimal control setpoints for each configuration that minimize the plant operating cost. Recall that a cooling tower is assumed to be dedicated to a chiller. Under the same chiller configuration, only those cooling tower configurations that have the least total energy consumption are used for regression. Note that the power consumption of the whole plant as well as the power consumption of an individual chiller (the dedicated cooling tower and pumps are also included) are regressed separately.

Regression results of the linear, second-order, and third-order polynomial models are given in Table 6, with the model order finally selected shown in boldface type. For groups G1, G2, and G3, the quadratic model captures the variability in the data quite well, with adjusted [R.sup.2] values being almost 100% and CV-RMSE being less than 0.1%. However, for G4 and G5, the quadratic models are poor; hence, higher-order models are necessary. An incomplete third-order polynomial form was finally found to capture the relationship between the inputs and output to an acceptable CV-RMSE level (less than 3%) (Jiang 2005).

In order to evaluate the external prediction ability of the fitted-response surface models, numerous additional cases (in the range of 50-100) have been simulated for each group and used to evaluate the predictive accuracy of the models. The analysis results are summarized in Table 7. It is clear that the CV-RMSE values of both internal and external predictions are very good (less than 3%) for G1, G2, and G3, indicating that their response surface models are robust with sound and reliable predictive ability. For G4 and G5, CV-RMSEs of external prediction for whole-plant energy use are less than 5%, which is an acceptable level. Though some of the individual equipment models are poor (for example, CV-RMSE values for G4, [P.sub.ele(abs)] = 11.7%, and G5, [P.sub.ele(cen)] = 7.4%), the corresponding plant models are acceptable (CV-RMSE values for G4, [P.sub.ele] = 4.7% and [P.sub.gas] = 1.2%). These are the models used in determining the minimal operating cost of the plant.

Optimization of the Hybrid Cooling Plant Operation under Real-Time Pricing

The methodology is applied to the hybrid cooling plant for minimizing the operation cost over a specified time period under RTP, with the planning horizon chosen to be 12 hours. This is a practical choice since this corresponds to the occupied period of the day for most commercial buildings and campuses. Since there is no demand cost under this electricity rate structure, Dijkstra's algorithm can be directly applied to solve the optimization problem.

Typical building cooling load profiles and ambient wet-bulb temperature profiles for a hot day and a mild day in Philadelphia were determined from TMY climatic data (NREL 1995), as shown in Figures 2 and 3. Cooling load profiles and rate structures used for simulation are artificial but realistic and would promote load-shifting and use of an absorption chiller during periods of high electricity rate. Three chillers are needed to satisfy the building cooling load during the peak of the hot day; however, only two chillers are needed to meet building cooling load during the peak period of the mild day. Figure 4 shows typical RTP electricity rate profiles (adapted from Henze and Krarti [1999]). For the gas rate, we chose two values, $0.4/therm and $1.2/therm, to study its effect on the optimization results. For convenience, we denote six conditions as listed in Table 7. After the optimal chiller sequencing path is determined, the SQP optimization program is used to find the optimal control setpoints.

The cost difference between using detailed static optimization and using the dynamic scheduling algorithm with response surface modeling approach is shown in Table 8. The results indicate that the least-cost path algorithm is very accurate, with the maximal error being 1.8% under RTP-6. This evaluation, therefore, serves as a validation of the proposed dynamic scheduling algorithm.

Figures 5, 6, and 7 show the chiller scheduling paths for diurnal conditions RTP-1, RTP-3, and RTP-6, with the bolded dash-dot line representing the least-cost path (or chiller operating strategy) and the dashed lines representing the various feasible scheduling strategies over the planning horizon. The labeled numbers beside lines denote the chiller feasible operating modes (groups). The following observations can be made:

* Under the RTP rate structure, least-cost paths are equally least-energy paths because no demand charge is considered.

* As expected, during those hours with low electricity rates (hours 8:00-12:00 and 18:00-20:00 for RTP-3, RTP-4, RTP-6), operating the centrifugal chiller is always preferred, even for a medium electricity rate (RTP-C) and high gas rate.

* An absorption chiller is preferred under the high electricity rate, even during a mild day (for example, from Figure 6, during the hours 14:00-16:00).

* On a hot day, there are not many options for scheduling chillers since all three chillers must be operated to satisfy the building cooling load; therefore, the relative benefit of optimal scheduling and control is not significant. However, on a mild day, with many more options available, the optimal operating strategy is more difficult to determine and the savings due to optimal control are more significant.

Optimization of the Hybrid Cooling Plant under TOU with Electricity Demand

Because of the demand cost component, the modified Dijkstra's algorithm is used to determine the least-cost path for operating the hybrid cooling plant over 12 hours. Factors studied under TOU with demand are building cooling load, wet-bulb temperature, electricity and gas rates, and electricity demand rate. We assume the same cooling load profiles and wet-bulb temperature profiles for typical hot and mild days in Philadelphia, as before. Usually, a TOU rate structure has two time periods--on-peak and off-peak hours. Both on- and off-peak hours can be different for different electricity rate profiles. Although this methodology can be applied to TOU with both on- and off-peak periods without any extra modification, we have assumed the planning horizon (from 9:00 to 20:00) to be on-peak hours, with electricity rates assumed to be $0.1/kWh and $0.15/kWh, and demand rates to be $5/kW and $20/kW. The gas rates assumed are $0.4/therm and $1.2/therm, representative of low and high values. Note that the 16 conditions summarized in Table 9 include various combinations of hot and mild days and the electricity rates, electricity demand rates, and gas rates listed previously. Results of the various optimizations performed can be found in Jiang (2005), while some salient observations follow.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

* Due to the inclusion of the demand cost component, the least-cost operating strategy is not the same as the least-energy strategy (see Figures 8-10). We note that during those hours that peak is likely to be hit, the absorption chiller replaces the centrifugal chiller, although running an absorption chiller consumes more energy. Also from Figure 8, under diurnal condition TOU-1, we note that the second centrifugal chiller is turned on at 13:00 in order to avoid the demand charge at 14:00, although running only two chillers could meet the cooling load.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

* Only when the electricity rates are so high that running the absorption chiller incurs lower energy cost than running the centrifugal chiller plus the avoided electricity demand charge is the least-cost path identical to the least-energy cost path. From Figure 10, under diurnal condition TOU-13 with electricity and gas rates being $0.15/kWh and $0.4/therm, from 12:00 to 16:00, running the absorption chiller incurs lower energy cost and the least-cost path overlaps with the least-energy path.

* Operating the absorption chiller is preferred whenever the demand charge occurs (Figures 8-10). Even with low electricity rate and low demand rate during a mild day (Figure 9), it is preferable to operate the absorption chiller in order to avoid the demand charge, although running it has higher energy cost. Obviously, with increased electricity rate (Figure 10), running the absorption chiller is preferred since it can avoid the demand charge as well as lead to lower energy cost.

* We find that different demand rates ($5/kW or $20/kW) yield similar optimal scheduling strategies in most cases (for example, TOU-1 and TOU-3, TOU-5 and TOU-7 of Table 9). This is because in most cases, even the lowest demand rate chosen ($5/kW) is still not low enough to make the demand charge lower than the extra energy cost (compared to the least-energy path). The only exception among these cases is, under diurnal conditions, TOU-2 and TOU-4, with the same electricity rate ($0.05/kW) and gas rate ($1.2/therm). When the demand rate is $20/kW, the absorption chiller starts replacing the centrifugal chiller at 11:00, which is one hour earlier than the $5/kW case. Here, the reason is that the demand charge incurred by running "0 1 1" at 12:00 is less than the energy cost difference between "0 1 1" and "1 0 1" and vice versa for the $20/kW case. Note, however, that these observations are valid only for the demand-setting day of the month.

* Similar to the RTP electricity rate structure, the gas rate affects the operation strategy. For example, because of the increased gas rate, the time to turn on the absorption chiller is postponed by one hour under diurnal condition TOU-2 as compared to that under diurnal condition TOU-1.

* During those hours with low cooling load or low electricity rate, it is preferable to run the centrifugal chiller because no demand charge would occur. For example, the centrifugal chiller is always preferred in the morning (9:00-11:00) and late afternoon (17:00-20:00) on a mild day.

* Due to the inclusion of demand cost, dynamic optimal scheduling on a hot day can still yield very significant benefits, although the alternative scheduling strategies on a hot day are less than those on a mild day.

A final evaluation was done in order to gauge the benefit of selecting an optimal dynamic strategy versus the simpler strategy of performing static optimization at each hour. The results, summarized in Figure 11, show that the cost savings by adopting our optimal dynamic operating strategy can be as much as 70% over a 12-hour planning horizon, although on average it is about 37.5% for the cases studied. It must be pointed out that these results represent the "worst case" performance of the static optimization case versus the least-cost operating strategy, since demand charges are applied only to a single peak-setting day in the month. However, the methodology can be applied to determining an optimal operating strategy over the whole month.

SUMMARY

This paper presents a general and computationally efficient methodology for optimal scheduling and control of primary HVAC & R plants involving numerous equipment based on different objective functions. The novelties of the method are that it is optimal in the number of static optimizations needed to develop a response surface model and also that it proposes a computationally efficient dynamic scheduling algorithm.

The proposed methodology involves two stages. The first stage relates to the static optimization case, where one determines the optimal cost of operating different combinations of equipment. The first step can be further split into four subproblems: (1) generate feasible operating modes under different operating conditions--select which equipment to run and the relative speeds for multi-speed fans or pumps; (2) generate, using experimental design techniques, a finite set of plant operating conditions including cooling loads, wet-bulb temperatures, and utility rates; (3) determine the optimal values of the continuous control variables using the minimization solver at each plant operating mode; and (4) perform response surface analysis on the minimal operation cost versus forcing variables (for example, cooling load, wet-bulb temperature, and utility rates) for each plant operating mode.

[FIGURE 11 OMITTED]

The second stage involves dynamic optimal scheduling using the computationally efficient Dijkstra algorithm, which was originally modified to handle multiple vapor compression chillers and which we further adapted to hybrid chiller plants. The algorithm is able to handle energy cost as well as costs due to electric demand peaks. A semi-real hybrid cooling case study is presented to illustrate the entire methodology and demonstrate how meaningful trends and practical heuristics can be identified for the specific plant under different price signals. The companion paper (Jiang et al. 2007) extends the deterministic component of the engineering analysis presented here to address the effect of various sources of uncertainty and how these could be combined within a decision analysis framework that includes risk attitudes of the operator.

ACKNOWLEDGMENTS

We thank Dr. Itzhak Maor for providing us with detailed information of an actual hybrid cooling plant, which served as the basis of the case study in this research, and also for useful advice and practical insights throughout this work. Insightful comments from Drs. P. Gurian, J. Wen, S.V. Smith, and J.R. Weggel are also acknowledged. The paper also benefited greatly from the detailed comments of the anonymous reviewers.

NOMENCLATURE 

COP = coefficient of performance 

CV-RMSE = coefficient of variation of the root mean square error

[C.sub.pw] = constant pressure specific heat of liquid water, kJ/kg x K

F = objective function

FMP = fan-motor power at rated conditions, kW

[h.sub.a] = enthalpy of moist air per mass of dry air, kJ/kg

K = total number of equipment or components

k = equipment or component index

m = index of equality constraints

[m.sub.a] = mass flow rate of dry air, kg/s

[m.sub.w] = mass flow rate of water, kg/s

NTU = number of transfer units of a heat exchanger 

n = index of inequality constraints, index in Equation 3

[n.sub.f] = number of factorial runs

P = energy consumption per unit time interval, kWh/h or therm/h

PLR = part-load ratio

[Q.sub.ch] = chiller cooling load, assumed equal to the building cooling demand, kW

R = unit cost of energy, $/kWh or $/therm

RTP = real-time pricing for electricity

[R.sub.de] = electricity demand rate, $/kW

SC = start-up cost, $

T = total number of time intervals

TOU = time-of-use rate for electricity pricing

t = time interval index

[T.sub.cdi] = condenser water inlet temperature, [degrees]C

[T.sub.cho] = chilled-water outlet temperature, [degrees]C

[T.sub.chs] = supply water temperature to building after mixing, [degrees]C

[T.sub.cto] = cooling tower water outlet temperature, [degrees]C

[T.sub.gni] = generator inlet temperature for the absorption chiller, [degrees]C

[T.sub.ref] = reference temperature for zero enthalpy of liquid water, [degrees]C

[T.sub.w] = water temperature, [degrees]C

[T.sub.wb] = ambient wet-bulb temperature, [degrees]C

u = chiller operating status; for example, at any time interval the status of, say, three chillers can be [0 0 1], which would signify the first two chillers are off and the third one is on

x = control variables

Subscripts

i = inlet

o = outlet

ele = electricity

f = fan

gas = gas

r = return

s = supply

ss = steady state

REFERENCES

ARI. 2000. Standard 560-2000, Absorption Water Chilling and Water Heating Packages. Arlington, VA: Air Conditioning and Refrigeration Institute.

ARI. 2003. Standard 550/590-2003, Performance Rating of Water-Chilling Packages Using the Vapor Compression Cycle. Arlington, VA: Air Conditioning and Refrigeration Institute.

ASHRAE. 2003. 2003 ASHRAE Handbook--HVAC Applications. Atlanta: American Society of Heating, Refrigerating and Air-Conditioning Engineers.

Brandemuehl, M. 1993. HVAC2 Toolkit, Algorithms and Subroutines for Secondary HVAC System Energy Calculations. Atlanta: American Society of Heating, Refrigerating and Air-Conditioning Engineers.

Braun, J.E. 1988. Methodology for the design and control of central cooling plants. PhD dissertation, University of Wisconsin, Madison.

Braun, J.E., S.A. Klein, and J.W. Mitchell. 1989. Effectiveness models for cooling towers and cooling coils. ASHRAE Transactions 95(2):164-74.

Curti, V., D. Favrat, and M.R. von Spakovsky. 2000. An environomic approach for the modeling and optimization of a district heating network based on centralized and decentralized heat pumps, cogeneration and/or gas furnace. Part II: Application. Int. J. Therm. Sci. 39:731-41.

Dijkstra, E.W. 1959. A note on two problems in connection with graphs. Numerische Mathematik 1:269-71.

Dotzauer, E. 1997. Algorithms for short-term production-planning of cogeneration plants. PhD thesis, Linkoping University, Linkoping, Sweden.

Dotzauer, E. 2003. Experiences in mid-term planning of district heating systems. Energy 28:1545-55.

Fletcher, R. 2001. Practical Methods of Optimization. Chichester, New York: John Wiley and Sons.

Gordon, J.M., and K.C. Ng. 2000. Cool Thermodynamics. Cambridge, UK: Cambridge International Science Publishing.

Henze, G.P., and M. Krarti. 1999. The impact of forecasting uncertainty on the performance of a predictive optimal controller for thermal energy storage systems. ASHRAE Transactions 105(2):553-61.

Horii, S., K. Ito, P.S. Pak, and Y. Suzuki. 1987. Optimal planning of gas turbine cogeneration plants based on mixed-integer linear programming. Int. J. Energy Research 11(4):507-18.

Hydeman, M. 2002. Optimizing the design and control of chilled water plants. Professional Developent Seminar (PDS), ASHRAE Learning Institute.

Jiang, W. 2005. Framework combining static optimization, dynamic scheduling and decision analysis applicable to complex primary HVAC & R systems. PhD dissertation, Drexel University, Philadelphia, PA.

Jiang, W., and T.A. Reddy. 2003. Re-evaluation of the Gordon-Ng performance models for water-cooled chillers. ASHRAE Transactions 109(2):272-87.

Jiang, W., T.A. Reddy, and P. Gurian. 2007. Combining engineering optimization of primary HVAC & R plants with decision analysis methods, Part II: Uncertainty and decision analysis. HVAC & R Research 13(1).

Koeppel, E.A., S.A. Klein, J.W. Mitchell, and B.A. Flake. 1995. Optimal supervisory control of an absorption chiller system. HVAC & R Research 1(4):325-42.

Montgomery, D.C. 1991. Design and Analysis of Experiments, 3d ed. New York: John Wiley & Sons.

NIST. 2005. Engineering Statistics Handbook, chapter 5. Gaithersburg, MD: National Institute of Standards and Technology.

NREL. 1995. TMY2S, Typical Meteorological Years Derived from the 1961-1990 National Solar Radiation Database. Golden, CO: National Renewable Energy Laboratory.

Olson, R.T. 1988. Optimal allocation of building cooling loads to chilled water equipment. PhD dissertation, University of Illinois at Urbana-Champaign, Urbana, IL.

Olson, R.T., and J.S. Liebman. 1990. Optimization of a chilled water plant using sequential quadratic programming. Engineering Optimization 15:171-91.

Olson, R.T., C.O. Pedersen, and J.S. Liebman. 1993. A dynamic procedure for the optimal sequencing of plant equipment, Part I: Algorithm fundamentals. Engineering Optimization 21(1):63-78.

Olson, R.T., C.O. Pedersen, and J.S. Liebman. 1994. A dynamic procedure for the optimal sequencing of plant equipment, Part I: Validation and sensitivity analysis. Engineering Optimization 22:163-83.

Papoulias, S.A., and I.E. Grossmann. 1983. A structural optimization approach in process synthesis. Part II: Heat recovery networks. Computers and Chemical Engineering 7:707-21.

Phelan, J., M.J. Brandemuehl, and M. Krarti. 1997. In-situ Performance testing of fans and pumps for energy analysis (RP-827). ASHRAE Transactions 103(1):318-32.

Reddy, T.A., and K.K. Andersen. 2002. An evaluation of classical steady-state off-line linear parameter estimation methods applied to chiller performance data. HVAC & R Research 8(1):101-24.

Sakamoto Y., A. Nagaiwa, S. Kobayasi, and T. Shinozaki. 1999. An optimization method of district heating and cooling plant operation based on genetic algorithm. ASHRAE Transactions 105(2):104-15.

Yin, S., and W. Wong. 2001. Hybrid simulated annealing/genetic algorithm to short-term hydro-thermal scheduling with multiple thermal plants. Electrical Power and Energy Systems 23:565-75.

Yokoyama, R., Y. Hasegawa, and K. Ito. 2002. A MILP decomposition approach to large scale optimization in structural design of energy supply systems. Energy Conversion and Management 43:771-90.

Wei Jiang, PhD

Associate Member ASHRAE

T. Agami Reddy, PhD, PE

Member ASHRAE

Received January 13, 2006; accepted April 28, 2006

Wei Jiang is a research engineer at Pacific Northwest National Laboratory, Richland, WA, and T. Agami Reddy is a professor in the Department of Civil, Architectural and Environmental Engineering, Drexel University, Philadelphia, PA.

Table 1. Summary of Component Specifications for the Hybrid Cooling
Plant

Component                          Number  Rated Specifications

Chillers          Direct-fired     CH1     Capacity: 2110 kW (600 ton)
                  absorption               Burner fan electricity
                    chiller                  consumption: 40 kW
                  Centrifugal      CH2     Capacity: 2110 kW (600 tons)
                    chiller        CH3
Cooling towers    For absorption           Cell number: 1
  with variable-    chiller                Fan speed: 384 rpm
  speed fan and                            Airflow: 131 kg/s
  fixed-speed                                (231,900 cfm)
  condenser                        CT1     Motor: 45 kW (60 hp)
  pumps                            CDP1    Flow rate: 136 kg/s
                                             (2160 gpm)
                                           Pressure drop: 38.1 m
                                             (125 ft)
                                           Motor: 67 kW (90 bhp)
                                           Pump efficiency: 75%
                  For centrifugal          Variable speed
                    chiller                Fan speed: 342 rpm
                                   CT2     Airflow: 108 kg/s
                                             (190,100 cfm)
                                   CT3     Motor: 30 kW (40 hp)
                                   CDP2    Flow rate: 113 kg/s
                                             (1800 gpm)
                                   CDP3    Pressure drop: 38.1 m
                                             (125 ft)
                                           Motor: 56 kW (75 bhp)
                                           Pump efficiency: 75%
Pumps with        Evaporator       CHP1    Flow rate: 91 kg/s (1440 gpm)
  fixed speed       water loop     CHP2    Pressure drop: 12.2 m (40 ft)
                                   CHP3    Motor: 15 kW (20 bhp)
                                           Pump efficiency: 75%

Table 2. Pertinent Information Relating to Cooling Plant Component
Models (The Pump Models are Not Shown Since They are Fixed Speed)

                                                    Regression Results
                          Model                   Adj-[R.sup.2]  CV-RMSE
Component Model           Equation  Coefficients  (%)            (%)

Chiller  Absorption       2         [b.sub.0] =   95.6            4.6
         (CH1)                        14342.443
                                    [b.sub.1] =
                                      -18093.925
         Centrifugal      1         [a.sub.1] =   94.3            3.2
         (CH2, CH3)                   0.226
                                    [a.sub.2] =
                                      -281.218
                                    [a.sub.3] =
                                      0.00689
Cooling  For absorption   3         [c.sub.1] =   74.8            7.7
  tower  (CT1)                        1.325
                                    [n.sub.1] =
                                      -0.426
         For centrifugal  3         [c.sub.2] =   72.5           10.1
         (CT2, CT3)                   1.281
                                    [n.sub.2] =
                                      -0.573
         Fan              5         [e.sub.0] =   99.8            2.2
         (CT1, CT2, CT3)              0.03247
                                    [e.sub.1] =
                                      -0.07882
                                    [e.sub.2] =
                                      1.1708
                                    [e.sub.3] =
                                      -0.1184

Table 3. Different Types of Variables Influencing Cooling Plant System
Operation

Forcing variables  [Q.sub.ch]         Building cooling load
  (exogenous       [T.sub.wb]         Wet-bulb temperature
  variables)       [R.sub.k]          Electricity/gas rate
Optimization       [Q.sub.ch,k]       Cooling load of each chiller
  variables          ([Q.sub.ch,1],
                     [Q.sub.ch,2],
                     [Q.sub.ch,3])
                   [T.sub.cto,k]      Outlet water temperature of each
                     ([T.sub.cto,1],    cooling tower
                     [T.sub.cto,2],
                     [T.sub.cto,3])
Control variables  [T.sub.cho,k]      Chilled-water setpoint of each
                     ([T.sub.cho,1],    chiller
                     [T.sub.cho,2],
                     [T.sub.cho,3])
                   [m.sub.a,k]        Cooling tower fan speed/air mass
                     ([m.sub.a,1],      flow rate of Each cooling tower
                     [m.sub.a,2],
                     [m.sub.a3])
Fixed variables    [m.sub.chw,k]      Chilled-water flow rate of each
                     ([m.sub.chw,1],    chiller
                     [m.sub.chw,2],
                     [m.sub.chw,3])
                   [m.sub.ctw,k]     Condenser water flow rate of each
                     ([m.sub.ctw,1]    chiller
                     [m.sub.ctw,2],
                     [m.sub.ctw,3])
                   [m.sub.chw]       Total chilled-water flow rate
                   [T.sub.chs]       Supply water temperature to
                                       building

Table 4. Constraints and Range of Different Variable Variation of the
Hybrid Cooling Plant System

Constraint  Description                 Range of Variation

 1          Constraint on               [T.sub.chs] = 6.7[degrees]C
              chilled-water               (44[degrees]F) (a)
              supply temperature
              after mixing
 2          Range of ambient            15.6[degrees]C (60[degrees]F)
              wet-bulb temperature        [less than or equal to]
                                          [T.sub.wb]
                                          [less than or equal to]
                                          29.4[degrees]C
                                          (85[degrees]F) (a)
 3          Range of evaporator         4.4[degrees]C (40[degrees]F)
              outlet water                [less than or equal to]
              temperature                 [T.sub.cho,k]
                                          [less than or equal to]
                                          10[degrees]C
                                          (50[degrees]F) (a)
 4          Range of condenser          18.3[degrees]C (65[degrees]F)
              outlet water                [less than or equal to]
              temperature                 [T.sub.cdo]
                                          [less than or equal to]
                                          37.8[degrees]C
                                          (100[degrees]F) (a)
 5          Range of condenser inlet    18.3[degrees]C (65[degrees]F)
              water temperature           [less than or equal to]
                                          [T.sub.cdi]
                                          [less than or equal to]
                                          29.4[degrees]C
                                          (85[degrees]F) (a)
 6          Range of airflow mass of    10% x [m.sub.a_rated,k]
              fan                         [less than or equal to]
                                          [m.sub.a,k]
                                          [less than or equal to]
                                          [m.sub.a_rated,k]
 7          Constraint on the chiller   [Q.sub.ch] =
              load                        [summation over k][Q.sub.ch,k]
 8          Due to common header for    [T.sub.chi] = [T.sub.chr]
              the evaporator inlet
              water temperature
 9          Due to common header for    [T.sub.cti,k] = [T.sub.cdo] =
              the cooling tower inlet     [summation over k]
              water                       ([T.sub.cdo,k] x
                                          [m.sub.ctw,k])/[m.sub.ctw]
10          Due to common sump for the  [T.sub.cdi] = [T.sub.cto] =
              condenser water             [summation over k]
                                          ([T.sub.cdo,k] x
                                          [m.sub.ctw,k])/[m.sub.ctw]
11          Due to common header for    [summation over k]([T.sub.cho,k]
              the evaporator outlet       x [m.sub.chw,k])/[m.sub.chw]
              water temperature
12          Energy balance on chiller   [Q.sub.ch] + P = [Q.sub.cd]
13          Constraint on minimum       0.15 x [Q.sub.rated,k]
              chiller capacity            [less than or equal to]
                                          [Q.sub.ch,k]
14          Constraint on maximum       [P.sub.ele,ch,k]
              chiller power draw          [less than or equal to]
                                          [P.sub.ele_rated,ch,k]

a. Suggested by ARI standards (ARI 2000, 2003).

Table 5. Feasible Configurations or Groups of the Various Plant
Components According to Their On/Off Status

           Status of Components
Group  CH1  CH2  CH3  CT1  CT2  CT3

G1     0    0    1    0    0    1
G2     1    0    0    1    0    0
G3     0    1    1    0    1    1
G4     1    0    1    1    0    1
G5     1    1    1    1    1    1

Note: 0 = on, 1 = off.

Table 6. Summary of Response Surface Models for Predicting Plant and
Individual Chiller Energy Use for All Five Groups--Significant Model
Coefficients Are Identified Using Forward Step-Wise OLS Regression

Group
(Number                        Third Order/Second Order/Linear
of Data  Energy                                     CV-RMSE
Points)  Consumption       RMSE            Internal (%)   External (%)

G1       [P.sub.ele]       --/0.28/5.1     --/0.11/2.0    --/0.14/--
(81)     [P.sub.gas]       --/--/--        --/--/--       --/--/--
G2       [P.sub.ele]       --/1.0/2.3      --/0.6/1.7     --/2.9/--
(81)     [P.sub.gas]       --/0.05/0.5     --/0.0/1.1     --/1.1/--
G3       [P.sub.ele]       --/0.6/10.1     --/0.1/1.9     --/1.1/--
(108)    [P.sub.gas]       --/--/--        --/--/--       --/--/--
G4       [P.sub.ele]        4.3/28.9/41.2  1.4/9.1/13.6    4.7/--/--
(108)    [P.sub.gas]       11/12.9/51.5    0.4/0.5/2.0     1.2/--/--
         [P.sub.ele(abs)]   0.4/0.4/1.1    0.3/0.3/0.8    11.7/--/--
         [P.sub.ele(cen)]   4.3/28.3/40.9  2.4/15.1/21.9   7.4/--/--
G5       [P.sub.ele]        5.9/25.5/27.7  1.1/4.5/4.9     4.2/--/--
(126)    [P.sub.gas]       19.1/63.0/71.7  0.7/2.4/2.8     2.3/--/--
         [P.sub.ele(abs)]   0.2/0.44/1.5   0.2/0.3/1.2     2.4/--/--
         [P.sub.ele(cen)]   5.3/39.2/47.8  1.2/9.0/10.9    5.0/--/--

Group
(Number                    Third Order/Second
of Data  Energy            Order/Linear
Points)  Consumption       Adj-[R.sup.2] (%)

G1       [P.sub.ele]       --/100.0/99.3
(81)     [P.sub.gas]       --/--/--
G2       [P.sub.ele]       --/92.7/77.7
(81)     [P.sub.gas]       --/100.0/99.3
G3       [P.sub.ele]       --/100.0/99.3
(108)    [P.sub.gas]       --/--/--
G4       [P.sub.ele]       98.5/72.1/43.2
(108)    [P.sub.gas]       99.5/99.3/89.6
         [P.sub.ele(abs)]  99.3/98.8/91.4
         [P.sub.ele(cen)]  98.2/71.7/40.9
G5       [P.sub.ele]       99.7/94.1/93.0
(126)    [P.sub.gas]       97.7/77.8/71.2
         [P.sub.ele(abs)]  97.7/98.8/86.4
         [P.sub.ele(cen)]  99.7/96.0/8.6

Note: The table shows results of RMSE, CV, and Adjusted [R.sup.2] for
third-order, second-order, and linear polynomial models. The model order
finally selected for the subsequent analysis is shown in boldface type.

Table 7. Diurnal Conditions Studied under RTP Rate Structure

Diurnal    Cooling Load*  Wet-Bulb                           Gas Rate
Condition  Profile        Temperature Profile*  RTP Profile  ($/therm)

RTP-1      hot day        hot day               RTP-A        0.4
RTP-2      hot day        hot day               RTP-A        1.2
RTP-3      mild day       mild day              RTP-B        0.4
RTP-4      mild day       mild day              RTP-B        1.2
RTP-5      hot day        hot day               RTP-C        0.4
RTP-6      hot day        hot day               RTP-C        1.2

*See Figures 3-5.

Table 8. Comparison of Static Optimization (Using SQP Algorithm) vs.
Least-Cost Path Algorithm under Different RTP Signals

                   Cost of Least-Cost Strategy  Cost Difference (a)
Diurnal Condition  ($)                          (%)

RTP-1              3055.8                       -0.33
RTP-2              3689.2                        0.035
RTP-3               663.3                       -0.20
RTP-4               691.0                       -0.027
RTP-5              1814.4                       -0.29
RTP-6              2267.7                        1.8

a. Cost difference between using detailed static optimization and the
dynamic scheduling algorithm with response surface models. The values
are close to zero, which demonstrates the accuracy of using the
computationally efficient response surface modeling approach compared to
the detailed static optimization.

Table 9. Diurnal Conditions Studied under TOU Rate Structure

                         Wet-Bulb     Electricity  Demand
Diurnal    Cooling Load  Temperature  Rate         Rate    Gas Rate
Condition  Profile*      Profile*     ($/kWh)      ($/kW)  ($/therm)

TOU-1      hot day       hot day      0.1           5      0.4
TOU-2      hot day       hot day      0.1           5      1.2
TOU-3      hot day       hot day      0.1          20      0.4
TOU-4      hot day       hot day      0.1          20      1.2
TOU-5      hot day       hot day      0.15          5      0.4
TOU-6      hot day       hot day      0.15          5      1.2
TOU-7      hot day       hot day      0.15         20      0.4
TOU-8      hot day       hot day      0.15         20      1.2
TOU-9      mild day      mild day     0.1           5      0.4
TOU-10     mild day      mild day     0.1           5      1.2
TOU-11     mild day      mild day     0.1          20      0.4
TOU-12     mild day      mild day     0.1          20      1.2
TOU-13     mild day      mild day     0.15          5      0.4
TOU-14     mild day      mild day     0.15          5      1.2
TOU-15     mild day      mild day     0.15         20      0.4
TOU-16     mild day      mild day     0.15         20      1.2

*See Figures 3-5.