- Email: [email protected]
A mixed-integer linear programming model for real-time cost optimization of building heating, ventilation, and air conditioning equipment
Accepted Manuscript Title: A Mixed-Integer Linear Programming Model for Real-Time Cost Optimization of Building Heating, Ventilation, and Air Conditioning Equipment Author: Michael J. Risbeck Christos T. Maravelias James B. Rawlings Robert D. Turney PII: DOI: Reference:
S0378-7788(17)30647-3 http://dx.doi.org/doi:10.1016/j.enbuild.2017.02.053 ENB 7414
To appear in:
ENB
Received date: Revised date: Accepted date:
23-3-2016 7-10-2016 23-2-2017
Please cite this article as: Michael J. Risbeck, Christos T. Maravelias, James B. Rawlings, Robert D. Turney, A Mixed-Integer Linear Programming Model for Real-Time Cost Optimization of Building Heating, Ventilation, and Air Conditioning Equipment, (2017), http://dx.doi.org/10.1016/j.enbuild.2017.02.053 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A Mixed-Integer Linear Programming Model for Real-Time Cost
cr
Conditioning Equipment
ip t
Optimization of Building Heating, Ventilation, and Air
an
October 7, 2016
us
Michael J. Risbeck,a,∗ Christos T. Maravelias,b James B. Rawlings,c Robert D. Turneyd
Abstract
M
In this paper, we present a framework for the formulation and solution of mixed-integer linear programming (MILP) models for operational planning of HVAC systems in commercial buildings. We introduce the general concepts of generators (e.g., chillers, boilers, cooling towers) and resources (e.g., electricity, chilled water), which allow
d
us to model a wide range of central plants. By using discrete variables for on/off states of central plant equipment and
te
continuous variables for equipment load, storage tank usage, and building temperature trajectory, all critical-to-cost decisions are taken into account. In addition, we consider both time-varying use charges and peak demand charges as well as building models of varying complexity. Because equipment models are often nonlinear, piecewise-linear
Ac ce p
approximations are used, which can be made arbitrarily accurate and enable real-time solution to the resulting optimization problems. We also employ a symmetry-free formulation to enhance the solution of the MILP model. Such features lead to improved performance compared to approaches employing heuristics or optimization without discrete variables. We demonstrate optimization for a small-scale cooling system and show favorable scaling of solution time with respect to number of units and temperature zones. In addition, we show simultaneous optimization of heating and cooling loops when a forecast of resource demand is available. Finally, we demonstrate closed-loop operation of our proposed optimization scheme. In all cases examined, a solution with an optimality gap below 1% can be obtained within 5 minutes, and thus the proposed optimization can be solved in real time, allowing for rapid response to changing weather, price, or demand forecasts. Keywords: online optimization; central plants; cost reduction; thermal energy storage; scheduling ∗Corresponding Author a University of Wisconsin–Madison. b University of Wisconsin–Madison. c University of Wisconsin–Madison. d Johnson Controls, Milwaukee, WI.
Page 1 of 38
1
Introduction
2
Heating, ventilation, and air conditioning (HVAC) equipment consumes a significant amount of energy within com-
3
mercial buildings. Such systems use (often simple) heuristic methods for scheduling and control [1], which result in
4
lower system efficiency. In addition, time-varying occupancy and ambient conditions lead to nonuniform electricity
5
usage [2–4], which causes utility companies to implement alternative pricing structures that strive to reduce peak elec-
6
tricity demand by consumers [5] because electricity production is most efficient when demand is relatively constant.
7
Such pricing structures include time-of-use charges (assessed using time-varying electricity rates) and peak demand
8
charges (assessed based on a customer’s one-time maximum rate of electricity consumption). With thermal energy
9
storage (TES) that can temporally shift electricity use, informed customers can leverage these pricing structures to
10
reduce total cost. This storage comes in the form of both active TES, such as purpose-built water or ice storage tanks,
11
and passive TES, which takes advantage of the natural thermal capacitance of buildings [6]. A combination of active
12
and passive storage techniques is often the most efficient solution [7, 8], but using these techniques requires a nontrivial
13
decision-making process [9].
an
us
cr
ip t
1
For these systems, a key challenge in planning using math programming methods is that there are both continuous
15
decision variables for equipment setpoints, utility consumption, storage usage, etc., and discrete decision variables
16
to specify whether process equipment is switched on or off. This complexity can be captured by mixed-integer pro-
17
gramming (MIP) models, which have been employed to make other building efficiency and sustainability decisions,
18
including selection of renewable and non-renewable energy sources [10–12], optimizing building retrofits based on
19
life-cycle analysis [13], and scheduling for “smart” appliances connected to a microgrid [14]. MIP models are most
20
computationally tractible when they have only linear constraints and objective (mixed-integer linear programming,
21
MILP), but more general mixed-integer nonlinear programming (MINLP) models have also been developed. Such
22
models can consider process energy balances in greater detail [15–18], but they are often limited to systems with few
23
pieces of equipment or few degrees of freedom.
Ac ce p
te
d
M
14
24
Because large-scale MIP models are computationally demanding, decomposition strategies are often employed. A
25
typical architecture is to is to make equipment selection choices in a small subproblem, either through near-optimal
26
heuristics [19–21], or via solving single-period scheduling problem [22–27]. Results of subproblems can then be used
27
to construct aggregate system curves [28], which simplify optimization at the higher level by eliminating all the dis-
28
crete, on/off choices for equipment. For asymmetric models of TES (i.e., models with different dynamics depending
29
on whether the system is charging or discharging), discrete choices must also be made for TES operation. Meth-
30
ods such as backward reach set computation [9] and dynamic programming [29] have been used to precompute the
31
charge/discharge profile, which also serves to simplify optimization at other levels. However, an important drawback
32
of such decomposition schemes is that equipment on/off states are often not considered over a long horizon, which
2
Page 2 of 38
means the effects of switching are not completely considered. Thus, units may be switched on and off rapidly in
34
response to small variations in system load, which must be addressed via further ad-hoc adjustments to the top-level
35
optimization problem. In addition, most approaches in the literature consider only chilled water production in a single
36
chiller, whereas the central plants we wish to address contain multiple units of different types and must produce hot
37
and chilled water concurrently.
ip t
33
In some cases, it may be desirable to determine demand by optimizing the temperature trajectory within buildings
39
using a model instead of simply following a fixed demand forecast. Building models can be either derived from first
40
principles [30] or identified from data [31]. Such models often split buildings into multiple interacting temperature
41
zones with separate control inputs, requiring the temperature trajectory of each zone to be considered. While buildings
42
often exhibit nonlinear dynamics, simplified linear building models can be sufficient for control [32]. The most
43
important feature of these models is that they relate occupant comfort to zone heating or cooling. Constraints on
44
occupant comfort are often implemented in terms of temperature, although more complicated metrics are possible
45
[33]. Studies have shown that such model-based control can out-perform conventional methods in terms of cost and
46
energy savings [34]. However, many existing approaches do not consider the complexities of the central plant (such
47
as on/off decisions or active storage) by assuming a constant coefficient of performance for cooling equipment. Thus,
48
a complete treatment must account for both buildings and the central plant.
M
an
us
cr
38
Accordingly, the goal of this paper is to propose an optimization-based framework that treats in detail both the
50
building and central plant subsystems. We propose a mathematical programming model that considers multiple pieces
51
of central plant equipment (for both heating and cooling) and multiple temperature zones over a finite prediction
52
horizon. In [35], an MILP formulation for the central plant subsystem was presented. This formulation approximates
53
critical-to-cost nonlinearities as piecewise-linear, which allows powerful MILP software to be used for optimization.
54
Improved performance compared to formulations without discrete decision variables was then demonstrated. In this
55
work, we extend the previous formulation to include a building model that allows optimal demand to be determined
56
online. In addition, we use a symmetry-free reformulation for identical pieces of equipment (not present in [35]) so
57
that explicit symmetry-breaking constraints are no longer necessary. These augmentations provide wider applicability
58
and more favorable scaling (with respect to horizon length and to number of central plant units) compared to the
59
previous formulation and also to alternative nonlinear formulations.
Ac ce p
te
d
49
60
The remainder of this paper is structured as follows. In Section 2, we present the problem statement and complete
61
MILP formulation (subsuming the formulation from [35]). In Section 3, we discuss the equipment models required
62
by the formulation. In Section 4, we present a series of simulations to validate the assumptions made by the model.
63
In particular, we discuss solution times for the combined building/central plant formulation; examine the cost savings
64
for a combined heating/cooling plant; address the handling of variable supply temperatures by benchmarking against
65
a more detailed (but computationally challenging) nonlinear formulation; and, demonstrate closed-loop operation as 3
Page 3 of 38
66
would be required in a real implementation. Finally, Section 5 summarizes the results of the study.
67
2
68
2.1
69
The overall objective is to maintain comfort of building occupants at the lowest economic cost. Depending on climate,
70
the focus may be on producing chilled water for building cooling in chillers or on producing hot water for building
71
heating in boilers. Nevertheless, hot and chilled water can be used for purposes unrelated to the HVAC system, and
72
thus simultaneous production of both resources must be considered. In these cases, heat-recovery chillers (HRCs) can
73
produce both hot and chilled water at higher efficiencies than separate units, as waste heat is rejected to a process hot
74
water stream instead of to the ambient. Thus, low-cost operation requires use of all types of equipment. An additional
75
challenge is that building temperature is often nonuniform with multiple interacting regions, each with their own
76
control inputs. Thus, a complete treatment of the problem requires simultaneous optimization of hot water production,
77
chilled water production, and building temperature trajectory.
ip t
cr
us
an
M
79
Problem Statement
The overall problem is as follows: given forecasts of hourly utility prices, ambient conditions and building occupancy, decide how to operate the HVAC system at lowest economic cost via the following decisions:
d
78
Optimization Problem
• What is the temperature profile of each zone?
81
• When is storage charging or discharging?
82
• Is each piece of equipment running?
83
• How much chilled and hot water is being producing?
Ac ce p
te
80
84
While making these decisions, temperature comfort constants must be satisfied, equipment minimum and maximum
85
capacities must be respected, and switching of equipment must be limited.
86
For simplicity, we assume that the supply and return temperatures of the chilled and hot water are known ahead
87
of time and fixed. While these decisions could be made within the optimization, it does not appear that variable
88
temperatures result in significant cost savings (See Section 4.3 for more details). In addition, cooling water supply
89
temperature is treated as a fixed parameter, although its value need not be constant throughout the prediction horizon.
90
Finally, we assume that building temperature dynamics can be sufficiently approximated using a linear model, or that
91
a forecast of building hot and chilled water demand is available.
4
Page 4 of 38
ip t
cr
Figure 1: Diagram of abstract optimization model. Within the load, primary demand is represented by the red, blue, and green curves at the top, and secondary demand is represented by the orange curve at the bottom.
2.2
Mathematical Formulation
93
Instead of formulating separate, individual constraints and variables for each type of equipment, we use an abstract
94
representation of the central plant in terms of “resources” and “generators.” This abstraction is similar to the state-task
95
network [36] and resource-task network [37] models that have been employed to great success for chemical production
96
scheduling [38–40]. Within the abstract model, a generator j ∈ J represents any piece of equipment in the central
97
plant (e.g., chillers, HRCs, pumps), while a resource k ∈ K represents any material flow in or out of the plant (e.g.,
98
chilled water, electricity). The prediction horizon is split into discrete intervals of equal length ∆ that are indexed by
99
t ∈ T. This overall abstract system is illustrated in Figure 1. In the abstract model, the system purchases resources
100
from the market, uses generators to convert these to other resources, and finally delivers these resources to the load,
101
optionally employing storage along the way. Once at the load, the resources are used to control the temperature
102
dynamics of zones i ∈ I or to serve auxiliary processes.
Ac ce p
te
d
M
an
us
92
103
We consider two different types of demand: (1) “primary demand,” which is is a decision variable and is used
104
to affect zone temperature (e.g., chilled water sent to air handlers); and, (2) “secondary demand,” which is a fixed
105
parameter that accounts for resource use by auxiliary processes (e.g., hot water for domestic use or chilled water for
106
specific internal cooling needs). Note that in different problem instances, demand may be entirely primary or entirely
107
secondary.
108
In the following sections, we use lower-case Roman letters for decision variables and subscripts, Greek letters for
109
fixed parameters, and bold upper-case Roman letters for sets. A complete listing of all symbols is given in Appendix A.
110
2.2.1
111
For a resource k and time t, primary demand in zone i is modeled as a variable gikt , while secondary demand φkt is
112
assumed to be known throughout the prediction period. Satisfaction of both types of demand is ensured through the
Demand Constraint
5
Page 5 of 38
113
following constraint: X
qjkt + ykt + pkt + bkt ≥ φkt +
j∈J
X
gikt ,
k ∈ K, t ∈ T.
(1)
i∈I
Thus, total demand is met via production qjkt , net storage withdrawal ykt , and direct purchase pkt . To soften this
115
constraint, the slack variable bkt ≤ φkt can be used, albeit at a heavy penalty in the objective function. By softening
116
the secondary demand constraint, the optimizer can still find a feasible solution even if secondary demand cannot
117
possibly be met (e.g., due to unit breakdowns). Note that primary demand cannot be ignored here; instead, temperature
118
comfort constraints are softened as described later.
cr
Stored resource inventory skt evolves based on its charge or discharge rate ykt according to the following dynam-
119
ics: skt = σk sk(t−1) + ykt ,
us
120
ip t
114
k ∈ K, t ∈ T,
(2)
where the term σk ≤ 1 is included to model storage inefficiencies, e.g., addition of heat to a chilled water tank from the
122
ambient. While TES dynamics are most accurately described by more complex effectiveness models with temperature
123
dependence [16], simpler state-of-charge models such as (2) have been shown to provide sufficient accuracy [6].
124
2.2.2
125
The true production/consumption relationships in the generators are often nonlinear, which requires enforcing con-
126
straints of the form
M
an
121
te
d
Equipment Models
Ξjt (ujt , qjk1 t , qjk2 t , . . . , qjkK t ) = 0,
that is, a nonlinear relationship between the production or consumption qjkt of the various resources, with additional
128
dependence on the on/off state uj of the generator. Note that these functions may be time-varying due to dependence
129
on external time-varying parameters (e.g., ambient temperature or humidity). For computational tractability, we wish
130
to solve the optimization as an MILP, and thus, we must use piecewise-linear approximations of the true nonlinear
131
model.
Ac ce p
127
132
For each type of generator j, we partition the (bounded) domain of Ξjt into a set Mj of polyhedrons, and we also
133
build the set Nmj of extreme points for each polyhedron m. For each extreme point n in each region m, a continuous
134
variable zjmnt is declared to choose the weighting of the given point. With these definitions, generator production and
135
consumption rates are calculated by
qjkt =
X
X
ζjknt zjmnt ,
j ∈ J, k ∈ K, t ∈ T,
(3)
m∈Mj n∈Njm
136
giving qjkt as a convex combination of the operating points stored in the parameter ζjknt .
6
Page 6 of 38
137
138
To decide how many generators of a given type are on, we use integer variables ujt , and to choose from among the available subdomains, we declare integer variables vjmt . Then, the following constraints are used: X
vjmt = ujt ,
(4)
j ∈ J, t ∈ T,
X
j ∈ J, m ∈ Mj , t ∈ T.
zjmnt = vjmt ,
n∈Njm
ip t
m∈Mj
(5)
Equation (4) chooses exactly one subdomain for each generator of type j that is on, and (5) distributes weight to the
140
chosen subdomains’ extreme points. Implicit in these constraints is that whenever ujt = 0, all corresponding qjkt = 0,
141
i.e., that when all generators are off, there is no production or consumption. By using integer-valued variables ujt and
142
vjmt for identical generators, we avoid computationally costly symmetry, which would be present if multiple sets
143
binary variables were used. Details of how the piecewise linearizations are constructed are presented in Appendix B.
144
The functional forms Ξjt of the true nonlinear models are presented in Section C.5 in Appendix C.
145
2.2.3
146
With integer variables to count the number of type-j generators that are on, rapid switching can be restricted via
147
constraints. First, the number of switching events is calculated using
an
us
cr
139
d
M
Switching Constraints
te
− ujt − uj(t−1) = u+ jt − ujt ,
j ∈ J, t ∈ T,
(6)
− where u+ jt and ujt count the number of switch-on and switch-off events for the given period. With these values, dwell
149
times are enforced as follows:
Ac ce p
148
ujt ≥
δj+ −1
X
u+ j(t−τ ) ,
j ∈ J, t ∈ T,
(7)
τ =0
µj − ujt ≥
δj− −1
X
u− j(t−τ ) ,
j ∈ J, t ∈ T,
(8)
τ =0
150
in which µj is the total number of generators of type j. Equation (7) requires that the number of “on” generators is
151
at least equal to the number of “switch on” events in the last δj+ time periods, while (8) enforces the corresponding
152
condition for “off” generators and “switch off” events.
153
2.2.4
154
Within large buildings, there are typically multiple independent (but interacting) temperature zones. Letting fit denote
155
the temperature in zone i at the end of time period t (i.e., at absolute time t∆, with ∆ the discrete timestep), and
Building Model
7
Page 7 of 38
156
declaring gikt to be the flow of resource k (e.g., cold water or hot water) into zone i during period t. With these
157
variables, temperature evolves according to ! fit =
X
λii0 fi(t−1) +
i0 ∈I
X
ωii0 k gi0 kt
+ θit ,
i ∈ I, t ∈ T.
(9)
ip t
k∈K
In this constraint, the coefficients λii0 and ωii0 k describe a linear time-invariant model for temperature evolution, while
159
θit accounts for time-varying disturbances. Note that while the energy flow into zone i does not instantaneously affect
160
the temperature of zone i0 , there is fill-in when the model is converted from continuous to discrete time, and thus ωii0 k
161
requires the double index ii0 . Details of this model can be found in Section 3.2.
cr
158
Since the objective is to maintain temperature bounds at minimum cost, constraints need to be introduced to enforce
163
temperature bounds. To avoid infeasibility, these constraints are softened by adding nonnegative slack variables fit+
164
and fit− as follows:
an
us
162
U fit ≤ ψit + fit+ ,
(10)
i ∈ I, t ∈ T,
(11)
M
L fit ≥ ψit − fit− ,
i ∈ I, t ∈ T,
L U with ψit and ψit the lower and upper comfort bounds. This constraint softening is to allow primary demand gikt
166
to be reduced without making the solution infeasible. Together with the slack bkt for secondary demand, softened
167
constraints ensure that the model has a feasible solution regardless of the demand and ambient temperature forecasts.
168
2.2.5
169
The objective function minimizes total resource cost ctot k , with additional penalty terms of soft constraint slack vari-
170
ables bkt , fit+ , and fit− as follows:
te
d
165
Ac ce p
Objective Function
min
X
k∈K
ctot k
+
! X
βkt bkt
+
t∈T
XX
+ − − χ+ it fit + χit fit .
(12)
t∈T i∈I
171
− Typically, the coefficients βkt , χ+ it and χit are be chosen quite large to ensure that the corresponding soft constraints
172
are satisfied whenever possible.
173
Economic cost ctot k considers both time-of-use and demand charges, calculated according to max max ctot pk + k = ρk
X
ρkt pkt ,
∀k ∈ K,
(13)
t∈T
pkt ≤ pmax , k
∀k ∈ K, t ∈ T.
(14)
8
Page 8 of 38
Here, ρkt and ρmax are time-of-use and demand charge prices (assumed known), while pkt and pmax k kt are the decision
175
variables representing how much and when a given resource is purchased. If peak demand charge is considered
176
only for certain times throughout the day (e.g., designated “peak hours”), constraint (14) need only be included for
177
the appropriate subset of periods. Multiple demand charges can also be considered by introducing additional pmax k
178
variables and constraints (14) for each demand “window.” For demand charges assessed over times longer than the
179
prediction horizon, pmax can take a lower bound Πmax equal to the prior maximum. Additionally, the weight ρmax of k k k
180
the demand charge can be increased throughout the period to reflect its greater importance near the end of the billing
181
cycle.
182
3
183
The optimization problem with constraints (1) through (11), (13), and (14), minimizing the objective function (12) is
184
referred to as MMILP . In the following subsections, we discuss the models used by the formulation and the generality
185
afforded by the abstract representation.
186
3.1
187
To ensure that the optimization problem correctly calculates total operating cost, accurate models of the generators’
188
production/consumption relationships are needed. Although these models are typically nonlinear, by replacing them
189
with piecewise-linear approximations as described in Appendix B, there are no restrictions on the algebraic form of
190
the models. However, a decision still must be made about what terms should be included as variables versus what
191
terms can be taken as parameters.
us
cr
ip t
174
M
an
Discussion
Ac ce p
te
d
Equipment Models
192
Numerous empirical models are available for chillers [41], often relating coefficient of performance (COP) as a
193
function of cooling load and one or more stream temperatures. For example, the Gordon-Ng model uses cooling
194
load, chilled water supply temperature, and cooling water supply temperature as independent variables [3]. For cost
195
optimization, the most relevant variable is cooling load, as it appears in the demand satisfaction constraint. To allow
196
for further optimization, supply temperature can also be considered as a variable (indirectly via volumetric flow; see
197
Section 4.3). Cooling water supply temperature, however, is mostly determined by ambient wet bulb temperature, and
198
it is thus treated as a (time-varying) parameter. Other models, including the Modified DOE-2 Model [42], empirical
199
polynomial fits, [6, 9, 43], and black-box lookup tables [15], use similar decision variables and can treated in the same
200
manner. It is also common to assume a constant COP [21, 29, 30, 44], which eliminates the need to consider any
201
stream temperatures. Thus, chiller models typically require at most two independent variables.
202
Because auxiliary equipment (e.g., pumps and cooling towers) comprises only a small part of total resource cost
203
and is often inflexible (i.e., has few or no degrees of freedom once other variables have been chosen), it is often 9
Page 9 of 38
ignored in optimization [42]. However, because auxiliary units affect peak demand charges, improved performance
205
(in terms of reduced cost) can be achieved by considering them even in approximate forms. For pumps, empirical
206
models of electricity consumption as a function of volumetric flow have been constructed from historical data [3, 27].
207
Cooling towers can be explicitly described by an effectiveness model [45, 46]. Alternatively, because cooling tower
208
performance depends on ambient temperature and humidity, these conditions can either be included in the chiller
209
power equation explicitly [15], or they can be used to estimate the optimal cooling water temperature in an empirical
210
function and then pass this value to chiller models [42]. In any case, choosing nonconvex or multivariable models for
211
auxiliary equipment can significantly slow optimization, and so care must be taken to balance model accuracy with
212
computational tractability.
213
3.2
214
To avoid unnecessary computational complexity, we use a linear dynamic model for temperature evolution in the
215
buildings. This model allows us to account for multiple interacting temperature zones and determine the optimal
216
temperature profile for each zone. Without any control action, the temperature trajectory is determined completely
217
by interactions with the external environment and among zones. Using an index i ∈ I for individual zones, a simple
218
model is
us
cr
ip t
204
M
an
Temperature Model
X dfi = −κoi (fi − f o ) − κii0 (fi − fi0 ) + gio , dt 0
d
αi
(15)
i∈I
in which fi is the temperature of zone i, αi is the thermal capacitance of zone i, f o is the external temperature, and
220
gio is the direct heat addition to zone i due, e.g., to solar radiation or heat production within the zone. In general, gio
221
is a function of ambient conditions and zone occupancy, but we assume a forecast is available. The coupling among
222
the zones and ambient is determined by the interaction coefficients κii0 and κoi . To provide control over temperature
223
evolution, a term gi is added to represent direct heat addition to (or rejection from) zone i by the controller.
Ac ce p
te
219
224
To use this model within the MILP formulation, we must convert from the continous-time ODE (15) to a discrete-
225
time difference equation of the form (9). Assuming all quantities except fi are constant between times t and t + ∆,
226
(15) is discretized using standard formulas to give
fi (t + ∆) =
X
(λii0 fi (t) + ωii0 gi0 ) + θi f o (t).
(16)
i0 ∈I 227
with the λ, ω, and θ parameters related to various matrix exponentials involving the zone masses and the interaction
228
coefficients κ. With minor modification, the ambient temperature f o can also be allowed to vary within each timestep,
229
e.g., using a first-order hold instead of a zero-order hold. In either case, the result is a linear difference equation, and so
230
it can be embedded within an MILP optimization problem assuming forecasts for f o and g o are available. Switching
10
Page 10 of 38
Return to Chillers
Hot Section
Return from Building Tank Discharging
Tank Charging
Supply to Building
ip t
Supply from Chillers
cr
Cold Section
us
Figure 2: Diagram of stratified tank model. Tank is treated as two separate volumes of water that interact with each other and the ambient. Only two of the four streams are active at a given time depending on whether the tank is charging or discharging. back to discrete-time index notation, t is treated as an integer, with xt denoting the (constant) value of quantity x
232
between times (t − 1)∆ and t∆ for the predetermined period length ∆. After lumping θi fi (t) into θit , we arrive at the
233
form given in (9). Note that we must add an extra index k to ωii0 so that the sign of gklt can be adjusted appropriately
234
for different resources. For example gklt is nonnegative for both hot and chilled water, but for chilled water, a negative
235
sign must be added to ensure that temperature within the zone is reduced.
236
3.3
237
A typical model for chilled water storage tanks is a single tank split into a hot and a cold section [47], as illustrated
238
in Figure 2. Whenever volume is withdrawn from one section, it is added to the opposite section, and thus total tank
239
volume remains constant. Cooling capacity is stored within the cold section of the tank, and the useful cooling amount
240
is determined by the total enthalpy therein. To arrive at the approximate model (2), a more accurate nonlinear tank
241
model is linearized at a nominal operating point, and at that operating point, we make the following assumptions:
Ac ce p
te
Storage Model
d
M
an
231
242
• The temperature of the cold section is equal to the temperature of the supply from chillers.
243
• The temperature of the hot section is equal to the temperature of the building return stream.
244
• Heat loss from the tank to the ambient is negligible compared to heat exchange between the sections.
245
A detailed derivation can be found in Appendix D.
246
Two potential issues with this model are that it overestimates the decay of enthalpy within the tank and that it does
247
not independently model the “cold volume” within the tank. If the simplified model (21) predicts significantly lower
248
efficiency for the storage tank, then the controller will be overly conservative with respect to utilization of storage, but
249
this effect could be attenuated by adjusting the decay coefficient σ. To model volume, the one-state linear model could
250
be replaced by a two-state linear model without significantly affecting solution times, but such an addition would 11
Page 11 of 38
be unlikely to affect optimal behavior. We note that with closed-loop operation and adequate state estimation, the
252
plant/model mismatch for the storage tanks is not a significant issue. Finally, the nonlinear stratified tank model could
253
be approximated more closely using a piecewise linearization similar to the nonlinear equipment models, but these
254
additions would increase solution times. Furthermore, because the stratified tank model is itself an approximation,
255
state estimation and closed-loop operation would still be recommended.
256
3.4
257
One of the key assumptions of this approach is that the equipment and building models are reasonably accurate.
258
While a detailed procedure for model fitting and maintenance is beyond the scope of this work, we recommend some
259
general procedures. For equipment models, empirical or semi-empirical models should be used so that they can be
260
periodically updated (e.g., every week) to account for the most recent historical data. For example, the unknown
261
coefficients for each chiller model in [41] can be found via linear regression, and the piecewise-linear approximation
262
is then updated by evaluating the new model at a finite number of operating points (see Appendix B). Because this
263
process is computationally cheap, it can be performed as often as necessary. However, while optimization time is
264
insensitive to small changes in the equipment models, more extreme changes (such as altering the convexity of an
265
equipment model) can significantly slow optimization. Thus, care must be taken to avoid major updates to the models,
266
for example by using nonlinear regression and penalizing large deviation from the old model. For the buildings, grey-
267
box models such as RC networks have been shown to be reasonably accurate [48] and also amenable for use with
268
standard system identification techniques. Because these models represent dynamics (and not just static relationships
269
like the equipment models), they require online use of a state estimator such as a Kalman Filter as in [49]. Note that
270
the building model may also need to be periodically updated to account for seasonal changes.
271
3.5
272
To more finely tune optimal behavior, constraints or fictitious costs can be added to MMILP . To reduce equipment
273
switching, a penalty coefficient for u+ it can be added to the objective function. Thus, the optimal solution would tolerate
274
higher economic cost if it can reduce the number of switches. To account for time-varying resource consumption by
275
the building but not by HVAC equipment, a forecast term can be added on the right-hand side of (14) so that the
276
demand charge is calculated using overall building-wide resource consumption. Such an addition could afford more
277
flexibility, as there may already be an intrinsic peak not created by HVAC equipment. Other HVAC-derived resource
278
consumption in zones (e.g., electricity consumption in zone air handlers) is also readily modeled by modifying the
279
gikt sum in the demand constraint (1). For recursive feasibility, a periodic constraint sk0 = skT (i.e., that the initial
280
and finale charges must match) could be added for the TES system. Similar periodicity can be enforced on zone
ip t
251
Ac ce p
te
d
M
an
us
cr
Model Validation
Further Considerations
12
Page 12 of 38
281
temperatures flt to ensure recursive feasibility with respect to zone comfort constraints. Pre-scheduled maintenance
282
or unit breakdowns are treated easily by setting appropriate ujt to zero. By using an abstract representation of the central plant, MMILP also accommodates additional resources and gen-
284
erators without reformulation. For example, with a suitable model for electricity production as a function of gas
285
consumption, a gas turbine could be considered simply by adding the corresponding index j and parameters ζiknt .
286
Absorption chillers can also be treated by adding steam (or other heat source) as a new resource k and augmenting the
287
appropriate models. Thus, this formulation is applicable to a wide variety of systems.
288
3.6
289
One of the key requirements of this formulation is that central-plant equipment can be treated as operating at steady
290
state within each sample time. Although there will certainly be transient startup and shutdown dynamics, we assume
291
that they resolve sufficiently quickly so that the static equipment models describe the actual equipment operating
292
characteristics throughout a large portion of the sample. Because discrete variables are present for equipment on/off
293
states, single-period deviations (for example, a flat increase in electricity consumption during the startup period)
294
could be considered in the cost calculation, but incorporation of detailed nonlinear dynamics would require significant
295
extension of this framework. In addition, for systems with a large number of temperature zones (greater than 50),
296
solving the complete problem in real time becomes more challenging, as the large number of continuous variables
297
needed to model the zones adds significant overhead to the MILP optimization. In such cases, it is necessary to
298
decompose the original problem into subproblems, for example a building temperature subproblem using a single
299
aggregate model of the central plant (e.g., constant COP) followed by a more detailed subproblem of just the central
300
plant. We will not discuss these ideas in this work, but we do note that the various subproblems can readily be
301
formulated as special cases of the full formulation described here, and thus the same considerations are relevant.
302
Finally, significant transport delays in delivering resources to the building are not explicitly considered, as they are
303
assumed to be short with respect to the sample time. Higher-order building models could be used to approximate
304
transportation time, although direct treatment may be necessary for significant delays.
305
4
306
In each of the following examples, price and demand profiles have been provided by Johnson Controls using data
307
from real buildings. Chillers and HRCs were modeled using Gordon-Ng models [41], while the boiler was assumed to
308
operate at constant efficiency. Pumps are modeled by empirical polynomial fits, and cooling towers are modeled with
309
a simplified explicit model [46]. All models are approximated as piecewise-linear, and coefficients for all equipment
310
models are provided by Johnson Controls. All timesteps ∆ are 1 h. Optimization was performed using Gurobi 6 on a
cr
ip t
283
Ac ce p
te
d
M
an
us
Limitations
Simulation Results
13
Page 13 of 38
Zone 1
g4o
Ambient
Zone 3
T1
Zone 2
fo
T3
T2
Zone 4
T4
Chiller 1 g1
1 2 Pumps
Chiller 2
Central Plant
Chiller 3
g3
us
Storage
ip t
fo
cr
g2o
Chiller 4
an
Figure 3: Example of combined system with building and central plant. desktop computer with 8 GB RAM and a 2.66 GHz × 4 CPU.
312
4.1
313
In this example, we model both the building and the central plant, with the equipment and temperature zones as shown
314
in Figure 3. For the building subsystem, we consider four temperature zones, two of which are “air” zones, and two
315
of which are “mass” zones. The air zones 1 and 3 must be maintained within comfort bounds, and they couple only
316
with the direct cooling input and the corresponding mass zones (2 and 4, respectively). The mass zones interact with
317
the ambient and are coupled to each other, but there are no comfort bounds for these zones. In effect, we have a
318
second-order linear model for temperature evolution within the two air zones, but in principle arbitrary order linear
319
models for each zone could be considered via appropriate choice of parameters. For the central plant subsystem, we
320
consider four conventional chillers and two pumps. There is also a chilled water storage tank available for use. Thus,
321
chilled water is produced in chillers and pumped either to the air zones to provide immediate cooling or to the storage
322
tank for later use.
M
311
Ac ce p
te
d
Building Cooling
323
Within the optimization, we use a time-varying forecasts for secondary demand, electricity price, and ambient
324
temperature (all of which are parameters). These profiles are shown in Figure 4, from which the daily cycle is clear.
325
In particular, we note that electricity price is highest when ambient temperature and cooling demand are highest, and
326
so there is significant potential for cost savings via load shifting. Solving MMILP with a 48 h horizon leads to a model
327
with 2.3k constraints and 4.4k variables, of which 1.4k are integer (after presolve to reduce the problem size). After
328
30 s, the optimality gap was 0.81%, and after 60 s, the gap was 0.17%. Due to nonconvexity in the pump models, it is
329
difficult to completely close the optimality gap, and thus the solver was stopped after 60 s. However, since equipment 14
Page 14 of 38
5.0 2.5 0.0
0
12
24
36
0
12
24
36
0
12
48
ip t
100 80 60
28 26 24 22 20 18 16
24
48
us
20
cr
40
36
an
Secondary Demand (MW)
7.5
Ambient Temp. (◦ C)
Electricity Price ($/MW)
10.0
48
Time (h)
M
Figure 4: Electricity price, auxiliary cooling demand, and ambient temperature for example system. and building models are only approximate, there is no need to solve to optimality in a practical implementation due
331
to feedback. The cooling and temperature dynamics for this solution are shown in Figure 5. We see that the optimal
332
profile is to charge the storage tank overnight (as shown by the large negative values for the “Storage” curve). Note
333
that “Demand” refers to total cooling demand, i.e., the sum of secondary demand from Figure 4 and primary demand
334
delivered to the temperature zones as determined by the optimization. For the temperature zones, there is no pre-
335
cooling, as zone temperature is left to float until it reaches the upper comfort bound of 22 ◦ C (recall that Zones 2 and
336
4 do not have comfort bounds, as they represent the surrounding solid mass for the air zones 1 and 4). Equipment
337
utilization shown as a Gantt chart in Figure 6. We see that equipment is switched off during peak hours to reduce
338
cost. However, due to the demand charge, it is more expensive to run all of the equipment at full capacity. Indeed, the
339
optimal peak loading would be very difficult to determine using heuristics alone, which illustrates the superiority of
340
optimization-based approaches.
341
4.1.1
342
Large commercial buildings have many temperature zones, and furthermore, it may become important to consider
343
longer horizons. Thus, we examine the scaling of our formulation with respect to building size by solving models of
344
various sizes with a 7 day horizon. Sizes and solution times for each optimization problem are shown in Table 1.
345
Note that changing the number of zones does not affect the number of discrete variables. Despite the slowdown as the
346
number of zones increases, performance is sufficiently fast for online use. We see that in all instances, a solution no
Ac ce p
te
d
330
Scaling with Number of Zones
15
Page 15 of 38
30 27 24 21 18 15
Demand Production Storage Unmet 0
12
24
36
48 Ambient Zone 1 Zone 2 Zone 3 Zone 4
0
12
24
36
ip t
Chilled Water (MW)
Temperature (◦ C)
60 40 20 0 −20 −40
48
cr
Time (h)
Chiller 1 Chiller 2 Chiller 3 Chiller 4 Pump 1 Pump 2 0
6
12
18
an
us
Figure 5: System and zone temperature dynamics for building cooling optimization.
24
30
36
42
48
M
Time (h)
Figure 6: Equipment utilization for building cooling optimization. Boxes indicate when a unit is on, while lines show equipment loading between minimum and maximum capacity. more than 1% of the best possible solution was obtained within just over 6 minutes. Recalling from that all models are
348
only approximate, a 1% optimality gap is sufficient for online use, in which optimization must be completed roughly
349
every 15 minutes.
te
d
347
Often, airside variables often do not appreciably change after the root node if sufficient active storage is available,
351
which allows a near-optimal forecast of demand to be obtained from even highly suboptimal solutions to the overall
352
problem. If such a forecast of primary demand is available, then the optimization problem need not include the building
353
model at all, and all demand terms are lumped into the secondary demand parameter φkt . To illustrate the benefits of
354
omitting the airside model, we consider the 25 zone example from before, sum the gklt variables, add them to φkt , and
355
repeat the optimization without the building model (9) through (11). As shown in the “No Building” entry in Table 1,
Ac ce p
350
Table 1: Problem sizes for 7 day horizon and varying number of zones. Sizes are as reported by Gurobi following pre-solve reductions. “1% Gap Time” is the wall-clock time needed to obtain a solution that is no more than 1% worse than the optimal solution. The “No Building” case simply meets a forecast of the optimal 25-zone demand. Size
Variables (k)
Integers (k)
Constraints (k)
1% Gap Time (s)
2 Zones 8 Zones 15 Zones 25 Zones No Building
6.2 12.3 18.7 29.2 4.4
1.8 1.8 1.8 1.8 1.8
3.8 7.8 11.9 19.0 2.5
4.0 26.0 106.3 365.5 2.0
16
Page 16 of 38
Table 2: Problem sizes for 7 day horizon and varying central plant size. Sizes are as reported by Gurobi following pre-solve reductions. “1% Gap Time” is the wall-clock time needed to obtain a solution that is no more than 1% worse than the optimal solution. Variables (k)
Integers (k)
Constraints (k)
1% Gap Time (s)
Nominal 2× 4× 8×
4.4 4.4 4.4 4.4
1.8 1.8 1.8 1.8
2.5 2.5 2.5 2.5
2.0 2.5 3.0 1.5
ip t
Size
the solution time for this instance is only 2 s. This result indicates that decomposition is be highly advantageous when
357
the number of zones is large.
358
4.1.2
359
Finally, we wish to address how the amount of equipment affects problem scaling. For these tests, we use the waterside-
360
only problem from the previous section and increase the amount of equipment present. Note that demand and storage
361
capacity are also scaled accordingly. Instance sizes and solution times are shown in Table 2. Here we see that scaling
362
with number of identical units is almost insignificant because on/off states are modeled using general integer variables
363
and not separate binary variables for each piece of equipment. Thus, adding additional units requires changing only
364
the upper bounds µj on ujt , which does not affect the number of variables or constraints, as seen in Table 2. Note
365
that if additional distinct types of equipment (e.g., chillers with different capacities) were added, then scaling would
366
be less favorable.
367
4.2
368
In central plants serving campuses or large buildings, simultaneous production of chilled and hot water is often re-
369
quired, e.g., due to secondary demand. Because of demand charges and heat-recovery chillers, these two process are
370
codependent and thus require simultaneous treatment. Below, sample optimization results are presented for the system
371
in Figure 7, which contains chilled and hot water loops (interacting via HRCs), both with active storage tanks and
372
pumps. The chillers also require cooling water produced in cooling towers, which are not shown in the diagram.
us
cr
356
te
d
M
an
Scaling with Number of Generators
Ac ce p
Combined Heating/Cooling Systems
373
In the optimization, a prediction horizon of 48 h is used. All generators are given on/off dwell times of 2 h, which
374
requires that after any switching event, that unit may not be switched again for at least 2 h. Forecasts for electricity
375
and gas price are given in Figure 8 along with heating and cooling demand forecasts. Here, the main periodicity
376
is in cooling demand and electricity prices, with a smaller effect in heating demand. For the system in Figure 7, a
377
solution with an optimality gap of 1.8% was found after 10 s of CPU time. Optimization was stopped after 300 s with
378
an optimality gap of 0.6%. Chilled and hot water production are shown in Figure 9, and equipment usage is given
379
in Figure 10. In this case, we see that production is decreased when electricity is expensive, although the effect is 17
Page 17 of 38
ip t cr us
12
an
Hot Water
M
0
24
Electricity
te
100 75 50 25 0
Chilled Water
36
48
36
48
Gas
d
120 90 60 30 0
0
12
24 Time (h)
Ac ce p
Price ($/MW)
Demand (MW)
Figure 7: Diagram of central plant for combined heating/cooling system. Storage tanks are available for both hot water (red) and chilled water (blue). Cooling towers (not shown) are connected to the chillers via a cooling water loop.
Figure 8: Resource demand and price forecasts. The central plant must supply time-varying demand hot and chilled water demand subject to time-varying electricity price and constant gas price.
380
nowhere near as pronounced as in the previous section. This difference is due to the demand charge and the reduced
381
flexibility of having to meet a fixed forecast of resource demand. A breakdown for electricity costs is provided in
382
Table 3. From this table, we see that chillers and HRCs are by far the largest contribution to total cost, indicating that
383
auxiliary units need not be modeled to a high degree of accuracy.
384
4.2.1
385
Because discrete decision variables are the main cause of computational difficulty, we compare MMILP to an alter-
386
native formulation that requires only continuous decision variables. By ignoring minimum capacities and switching
387
constraints, the piecewise-linear approximations of nonlinear generator models are replaced with single linear approx-
388
imations. These new models are quite inaccurate, but they no longer require the discrete variables vjmt associated
Comparison to Simpler Models
18
Page 18 of 38
Demand Production Storage Unmet
0
12
24
36
48 Demand Production Storage Unmet
20 10 0 −10
0
12
24
36
48
ip t
30
cr
Chilled Water (MW) Hot Water (MW)
60 45 30 15 0 −15 −30
Time (h)
us
Figure 9: Production of chilled and hot water for MMILP . Negative values for “Storage” indicate that storage tanks are being charged, while positive values denote discharge.
M
an
Chiller 1 Chiller 2 Chiller 3 Chiller 4 HRC 1 HRC 2 Boiler 1 Cold Pump 1 Cold Pump 2 Hot Pump 1 Hot Pump 2 Tower 1 Tower 2 6
12
18
d
0
24
30
36
42
48
Time (h)
Ac ce p
te
Figure 10: Equipment utilization for MMILP . Colored bars indicate that a unit is on, while solid lines show loading fraction between minimum and maximum capacity. Note that a first-in first-out queueing system was applied to disaggregate integer variables ujt into individual on/off states.
389
with piecewise-linearization, and indeed the on/off variables ujt can be ignored, as a unit’s on/off state is determined
390
by whether any qjt 6= 0. We refer to this scheme as MLP . Using the same parameters given in Figure 8, MLP is solved
391
as a linear program, and the equipment utilization plot in Figure 11 is obtained. In particular there is much more
392
equipment switching for MMILP compared to MMILP , as MLP cannot directly constrain equipment on/off switching.
393
Note that there are other solutions with the same total cost due to non-uniqueness of this formulation.
394
Table 4 compares costs for MMILP and MLP . Costs as predicted by the (piecewise-) linear models are shown in
395
the “Model” columns To provide a fair comparison, all electricity costs are recalculated using the underlying nonlinear
396
equipment models and shown in the “Actual” columns. In this example, we see that the fully linear models of MLP
397
are reasonably accurate due to units operating near capacity (which is where the linear approximations were made).
398
However, even in this favorable case, MMILP is noticeably more accurate (as demonstrated by the small difference
399
between “Model” and “Actual” values of electricity consumption). Furthermore, when the optimal policy requires
400
variable operation, there may be significant inaccuracy imposed by the linear models of MLP relative to the more
19
Page 19 of 38
Table 3: Electricity charges by equipment type for MMILP . Note that “Cost” column accounts for time-varying use charges, and so percentages do not exactly match with “Usage.” Usage, MWh
Cost, 1000$
168.8 106.7 24.8 6.0 1.6
(54.83%) (34.64%) (8.06%) (1.95%) (0.52%)
11.14 6.93 1.62 0.39 0.11
(55.16%) (34.33%) (8.05%) (1.94%) (0.53%)
Total
307.9
(100.00%)
20.19
(100.00%)
cr
Chillers HRCs Cold Pumps Hot Pumps Cooling Towers
ip t
Equipment
0
6
12
18
an
us
Chiller 1 Chiller 2 Chiller 3 Chiller 4 HRC 1 HRC 2 Boiler 1 Cold Pump 1 Cold Pump 2 Hot Pump 1 Hot Pump 2 Tower 1 Tower 2
24
30
36
42
48
M
Time (h)
401
accurate piecewise-linear models of MMILP .
d
Figure 11: Equipment utilization for MLP . Equipment models were approximated using a single linear term, and integer decision variables were not used.
In addition to greater model accuracy, MMILP also out-performs MLP in terms of unit switching. As seen in
403
Figure 11, there times when a unit is switched on for one period and then immediately switched off. Such behavior
404
cannot occur with MMILP because 2 h dwell times are enforced directly. Implicit in the problem formulation is that
405
equipment dynamics are irrelevant for the given timesteps, but without being able to enforce dwell times, solution
406
quality of MLP may deteriorate significantly when applied to the true system due to the burden (both in terms of
407
resource consumption and equipment wear) of extra switching. Although various ad-hoc additions could be made
408
to deter rapid switching, the discrete on/off nature of equipment scheduling decisions cannot be fully captured by
409
continuous variables. Thus, while MLP is much easier to solve due to its lack of discrete decision variables, it
410
performs noticeably worse than MMILP due to its significant limitations.
411
4.2.2
412
Currently, HVAC central plants with storage tanks are operated using heuristics instead of any sort of online opti-
413
mization [4, 9, 16, 47]. As will be demonstrated, heuristic methods can lead to significant cost increases if heuristics
414
fail to fully account for pricing structure. One such heuristic is to run equipment at full capacity over night to fully
415
charge storage tanks, and then keep equipment off during peak hours. To simulate this heuristic, we embed these
Ac ce p
te
402
Comparison to Heuristic Methods
20
Page 20 of 38
Table 4: Electricity cost comparison for MMILP , MLP , and heuristic approaches. “Change” row is total cost relative to Formulation MMILP Optimal MLP Heuristic MMILP Charges Optimal MMILP (1000$) Model Actual Model Actual Model Actual 20.21 2.24 22.45
20.19 2.24 22.43
20.32 2.21 22.53
21.29 2.63 23.93
21.05 3.59 24.63
21.02 3.59 24.61
Change
−
−
+0.37%
+6.68%
+9.73%
+9.72%
0
6
12
18
an
us
cr
Chiller 1 Chiller 2 Chiller 3 Chiller 4 HRC 1 HRC 2 Boiler 1 Cold Pump 1 Cold Pump 2 Hot Pump 1 Hot Pump 2 Tower 1 Tower 2
ip t
Usage Demand Total
24
30
36
42
48
M
Time (h)
Figure 12: Equipment utilization with heuristic rules enforced in the model. All chillers must be fully on from 1am to 7am and fully off from 2pm to 6pm. HRCs are left free. considerations as hard constraints in the optimization problem. Electricity cost is shown in Table 4, and equipment
417
usage is shown in Figure 12. In this case, the boiler must be activated to provide sufficient hot water, as production
418
is lost due to the heat recovery chillers being shut off in the afternoon. From these results, we see that the heuristic
419
solution does not perform well, even though a majority of its decisions were made optimally. Had the entire schedule
420
been determined by similar rules, performance would deteriorate even further. The main problem with such heuristics
421
is that the presence of peak demand charges requires choosing an upper bound for total electricity that accounts for
422
the tradeoff between demand and use charges. Although more sophisticated heuristics could better account for peak
423
demand charges, it is unlikely that performance would approach optimization-based methods.
424
4.3
425
The load Q supplied by a chiller can be varied by adjusting both the volumetric flow V through the chiller and the
426
supply temperature TCHWS . In the previous examples, a degree of freedom has been eliminated by assuming a fixed
427
TCHWS . While fixed supply temperature is a common practice in existing buildings, there is potential for cost savings
428
by allowing both load and supply temperature to vary independently. The main reason to allow TCHWS to vary is
429
that there is a tradeoff between chiller electricity usage and pump electricity usage that depends on TCHWS . Thus, to
430
deliver a specified load Q with the lowest electricity consumption potentially requires a different supply temperature
Ac ce p
te
d
416
Variable Stream Temperature
21
Page 21 of 38
431
for each possible on/off configuration of pumps and chillers. Based on common forms of chiller models [41], the most
432
relevant set of variables is (Q, TCHWS ), although if the temperature of water entering the chiller is known, then any
433
two of {Q, TCHWS , V } can be chosen. The optimal supply temperature is given by some relationship TCHWS = f (Q, u), in which Q is a vector with the
435
loads of all chillers, and u is a binary vector with the on/off states of all chillers and pumps. Even if we could find this
436
function f via offline optimization, there are too many degrees of freedom to embed a piecewise linear approximation
437
while maintaining fast solution times. Thus, a better alternative is to pick the optimal TCHWS online, which can be
438
considered by linearizing the chiller models in two variables instead of one. However, this addition can significantly
439
increase the computational burden, and so it is important to know what cost savings, if any, can be obtained.
440
4.3.1
441
We first wish to determine whether the piecewise-linear approximations made by MMILP are sufficiently accurate to
442
determine the optimal supply temperature. To assess this accuracy, we compare MMILP to a more detailed alternate
443
formulation MMINLP that uses the full nonlinear equipment models and includes energy balances to accurately cal-
444
culate the temperatures of the cooling water streams (which are assumed constant in MMILP ). Full details of MMILP
445
can be found in Appendix C. Due to its computational complexity, MMILP is intractable when multiple time periods
446
are considered, and so we perform a series of 1-period optimizations to calculate the optimal equipment configuration
447
and supply temperature for varying total system loads. Optimal supply temperature using the piecewise-linear and
448
nonlinear models are shown in Figure 13. From this plot we see that the accuracy of the piecewise-linear approx-
449
imation is quite good with the only noticeable differences are near switch points where MMINLP and MMILP select
450
different numbers of pumps. Indeed, the maximum temperature difference between the two models is below 3 ◦ C, and
451
at loads above 75% capacity where low cost is particularly important, the temperature difference is well below 1 ◦ C.
452
Thus, even with the simplifications made by MMILP , it is still possible to accurately determine the optimal supply
453
temperature.
454
4.3.2
455
To assess the effects of supply temperature on total electricity cost, we also calculate the optimal electricity consump-
456
tion (as predicted by MMILP and MMINLP ) as a function of varying total load. To provide a fair comparison, we obtain
457
the loads for each piece of equipment using optimization and calculate electricity consumption using the full nonlinear
458
models. These values are shown in Figure 14. We see once again that there is good agreement between MMILP and
459
MMINLP , indicating that the piecewise-linear approximations are quite accurate. We also observe that there is only
460
small savings to be gained by allowing variable water temperature. While the maximum individual savings is 4.1%
us
cr
ip t
434
Ac ce p
te
d
M
an
Accuracy of MMILP
Electricity Cost Savings
22
Page 22 of 38
11
Return MMILP Supply MMINLP Supply
10 9 8 7 6 5 0
10
20
30
40
50
cr
Total Load (MW)
ip t
Water Temperature (◦ C)
12
us
Figure 13: Optimal supply temperature using the MMILP and MMINLP models. The MMILP model uses piecewiselinear models, while the MMINLP model uses nonlinear models and includes cooling water energy balances. Both models have variable supply temperature. at a relatively low total load of 4 MW, when the plant is operating above half capacity, savings are significantly re-
462
duced. Average savings (weighted by total load) are 1.2%. Because of the increased computational burden, online
463
optimization of supply temperature may not worthwhile. However, a strength of MMILP is that this tradeoff can be
464
evaluated.
465
4.4
466
We conclude with an example of the system being run in closed loop. For simplicity of discussion, we consider a
467
system with only chillers. As described in Section 3.3, the model used for the storage tanks is only approximate,
468
and so we provide state feedback by giving the optimizer an accurate estimate of tank enthalpy (i.e., initial condition
469
sk0 ) at each time point. Figure 15 shows electricity price and equipment utilization over a 1 week period. In this
470
case, the electricity price and cooling demand were chosen to be periodic, and thus the operating profile is periodic as
471
well (excepting for small start and end effects). Note that each optimization has a a 24 prediction horizon with a 1 h
472
timestep.
M
an
461
Ac ce p
te
d
Closed-loop Operation
473
For comparison, we also simulate the system in open loop, i.e., without providing the actual charge of the storage
474
tank. Enthalpy, temperature, and volume profiles for the storage tank are shown in Figure 16. In this figure, only the
475
first 72 h of operation are shown, after which open-loop operation leads to tank overflow. From this simulation, we
476
see that open-loop operation leads to instability due to model mismatch. Because this model overestimates decay of
477
energy, open-loop operation leads to over-charging the tank, and eventually leads to a tank completely full of cold
478
water. By contrast, when the system is operated in closed-loop, the optimizer is aware of the true current charge of the
479
tank, and thus does not over-fill. Furthermore, the one-step prediction for tank enthalpy using the simplified model is
480
sufficiently accurate, as is the estimated temperature rise when the tank is non-empty. This agreement indicates that
481
the simplified model is accurate enough for closed-loop operation. 23
Page 23 of 38
6 4 2 0
0
10
0
10
20
30
40
20
30
50
cr
1.02
Relative Load (-)
ip t
Variable MMILP Constant MMILP Variable MMINLP Constant MMINLP
8
us
Electric Load (MW)
10
1.00 0.98
40
50
an
0.96
Thermal Load (MW)
d te
0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 0.0
Ac ce p
Electricity Price ($/MWh)
M
Figure 14: Comparison of optimal electricity consumption. Legend refers to whether chilled water supply temperature is variable or constant, and whether the approximate MMILP or full nonlinear MMINLP was used. Bottom plot is normalized to Constant MMILP .
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
Chilled Water (MW)
45 30
Demand Production Storage Unmet
15
0
−15
−30 0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
Chiller 1 Chiller 2 Chiller 3 Chiller 4
0.0
Time (days)
Figure 15: Electricity price, production, and equipment utilization for closed-loop simulation.
24
Page 24 of 38
0
CL Model
CL Actual
12
24
36
48
60
72
0
12
24
36
48
60
72
0
12
24
36
48
cr
Volume m3
8 6 4 2
Time (h)
us
Temp. Rise (◦ C)
10
0
ip t
8 7 6 5 4 3 2 1 0
OL Actual
Enthalpy (MJ)
OL Model 0 −50 −100 −150 −200 −250 −300 −350
60
72
an
Figure 16: Dynamics of storage tank for closed-loop and open-loop operation. More negative values of enthalpy indicate more stored cooling capacity. Temperature rise is calculated relative to the starting tank temperatures.
5
Conclusions
483
In this work, we have proposed a framework for the formulation of mathematical programming models for the op-
484
timization of equipment usage in a commercial HVAC systems. This framework accounts for time-varying utility
485
prices, demand charges, occupant comfort constraints, and on/off equipment dwell times. Using a series of reformu-
486
lations, including piecewise linearizations of unit models and an aggregate formulation for equipment switching, we
487
have developed a mixed-integer linear programming model, MMILP , that can be solved quickly and leads to solutions
488
with lower cost and less equipment switching than a similar scheme without discrete variables. In its most general
489
form, MMILP includes a linear model of building dynamics that is used to estimate cooling or heating loads required
490
to induce a given temperature trajectory in the building zones. Thus, resource demand can be determined within the
491
optimization, and use of passive TES is implicit. Alternatively, if a forecast of resource demand is available, a building
492
model need not be included at all, and MMILP can determine optimal usage of central plant equipment and active TES.
493
This formulation requires static models for resource consumption as a function of load for central plant equipment, as
494
well as dynamic models for storage tank charge and (possibly) building temperature evolution, but these models can
495
be obtained from manufacturers’ or historical data and can be updated as additional operating data is obtained.
Ac ce p
te
d
M
482
496
With respect to computational tractability, we have demonstrated that the formulation scales favorably with number
497
of temperature zone and number of units in the central plant. For a 7 day horizon with 1 h timestep, quality suboptimal
498
solutions can be found for a 25 zone building within 400 s, and with a demand forecast instead of the building model,
499
such solutions are found after only 10 s. For a central plant with both heating and cooling equipment, we show 6.7%
25
Page 25 of 38
improvement over a simpler formulation without discrete variables and 9.7% improvement over a simulated heuristic.
501
Thus, MILP-based optimization is the best choice for equipment scheduling in central plants. In addition, we have
502
shown that while variable supply temperature can be considered by the formulation, it does not significantly decrease
503
total cost, with a maximum reduction of 4.1% at low total loads and an average reduction of 1.2% across the entire
504
operating range. Therefore, a fixed supply temperature assumption is suitable for online use, but variable temperature
505
can be considered as well with additional computational effort. Finally, we have shown that while the simplified first-
506
order storage tank model is not accurate for long-term open-loop predictions, performance does not deteriorate in a
507
closed-loop implementation with measurement feedback. Based on its flexibility and modest computational time, our
508
model is suitable for use in new or existing buildings, and it can achieve significant cost reduction relative to heuristic
509
scheduling methods.
510
6
511
Funding, equipment models, and sample data provided by Johnson Controls, Inc.
512
Appendix A
513
A.1
514
The following sets are used in the abstract model description.
M
an
Acknowledgments
Ac ce p
temperature zones
te
d
Notation
Sets
i∈I 515
516
us
cr
ip t
500
j∈J
type of generators
k∈K
resources/utilities
t∈T
time periods
The following sets are used for piecewise linearization.
517
m∈M
interpolation regions for equipment models
n∈N
interpolation points for equipment models
Mj ⊆ M
interpolation regions for equipment j
Njm ⊆ N
interpolation points contained in region m for equipment j
518
A.2
Parameters
519
The following parameters are given based on available forecasts or desired system performance.
26
Page 26 of 38
523
demand charge
βkt
penalty for unmet secondary demand
Ξjt ( · )
nonlinear model function
L U ψit , ψit
lower and upper temperature comfort bounds for zones
Πmax k
previous maximum resource demand
ip t
ρmax k
cr
forecast resource prices
Number of units of each type available for use
ζjknt
interpolation points for piecewise-linear models
σk
fractional retention of stored resource
δj+ , δj−
minimum on/off dwell time for generators
Πk
Upper bound for resource purchase
Υk
Bound for single-period storage charge/discharge
Σk
Storage capacity
an
µj
us
The following parameters are calculated from equipment models.
M
522
ρkt
The following parameters are calculated from building models. αi
thermal capacitance of zone
κoi , κii0
heat transfer coefficients with ambient and between zones
λii0 , ωii0 k
coefficients in discretized model for building temperature evolution
θit
time-varying temperature disturbance for zones
d
521
forecast secondary resource demand
te
520
φkt
Ac ce p
524
− χ+ it , χit
penalties for comfort violation in zones
Γik
upper bound on resource flow to zones
525
Decision Variables
526
The following decision variables are primarily for the waterside subsystem. ujt ∈ {0, . . . , µj }
integer number of generators on during a given time
− u+ jt , ujt
number of generators switched on/off at the beginning of a time period
qjkt
production of resources in equipment
pkt ∈ [0, Πk ]
amount of purchased resources
pmax ∈ [Πmax , ∞) k k
maximum single-period purchase of resource k
ykt ∈ [−Υk , Υk ]
charge (< 0) or discharge (> 0) of stored resources
skt ∈ [0, Σk ]
stored inventory at the end of a time period
bkt ∈ [0, φkt ]
unmet secondary demand
527
27
Page 27 of 38
528
The following variables are added for the piecewise-linear approximations of equipment models. vjmt ∈ {0, . . . , µj }
number of units operating in each interpolation region
zjmnt ∈ [0, µj ]
weighting of interpolation points in each region
529
zone temperature of zone at end of a time period
gikt ∈ [0, Γik ]
primary demand of resources by zones
fit+ , fit−
positive and negative comfort violation for zones
532
Appendix B
533
The general relationship we are trying to capture is
Piecewise Linearization
x ∈ X ⊂ Rn ,
with X and Y compact. We assume that X can be partitioned into a finite number of polyhedrons {P1 , . . . , PM }, with SM X = i=1 Pi and int(Pi ) ∩ Pj = ∅ ∀i, j 6= i (that is, the polyhedrons only intersect on their boundaries).
M
535
y ∈ Y = f (X) ⊂ Rm
an
y = f (x),
534
ip t
fit
cr
531
The following decision variables are used when the building airside model is included.
us
530
Since we are modeling an equality relationship, we can choose either a relaxation of f or an approximation of f .
537
To relax f , we create a set-valued mapping F from X to Y with the requirement that f (x) ∈ F (x) ⊂ Y and use
538
the constraint y ∈ F (x) The advantage here is that any solution to the exact nonlinear model is a valid solution to
539
the relaxed model. However, there is now an additional degree of freedom that is not present in the actual nonlinear
540
formulation, and thus the optimal solution using the relaxed model is often infeasible to the true model. Alternatively,
541
we may simply approximate f by creating a piecewise-linear function f˜ such that |f˜(x) − f (x)| ≤ ε ∀x ∈ X with ε
542
some predetermined tolerance. The advantage here is that f˜ is actually a function and does not afford the system any
543
additional flexibility. However, if nonlinear constraints are tight, the approximate model may be infeasible even when
544
the true model was feasible. In our case, the variables corresponding to y are generally associated with production
545
amounts and do not have any tight feasibility requirements; thus, we use approximations of f .
Ac ce p
te
d
536
546
For a single dependent variable (i.e., y ∈ I ⊂ R), algorithms exist to find efficient piecewise linearizations for
547
convex and nonconvex functions [50, 51], and for two variables, triangulations can be used [52, 53]. We could make
548
use of these techniques by splitting f into component functions f1 , . . . , fm and use a different partitioning of X for
549
each, but this greatly increases the size of the resulting MILP. Thus, we require that a single partitioning be used for
550
all components of f , and so many existing techniques cannot be applied. Instead, we simply use a regular partitioning
551
of X into simplexes and use a secant approximation on each. This strategy ensures continuity of the resulting ap-
552
proximation f˜, and it also means f˜(x) = f (x) at extreme points of the Pj . Although such an approximation would
553
likely not minimize error throughout the entire domain, this feature is desirable because the optimal solutions to the 28
Page 28 of 38
554
MILP will often lie on these extreme points. To choose the final partitioning, we simply increase the resolution until
555
the maximum approximation error is below some prespecified value. Though more sophisticated techniques could be
556
developed to create more accurate approximations with fewer subdomains, the proposed scheme is sufficient for the
557
cases we consider.
559
The nonlinear relationships we must model are between the qjkt variables for a given j and t. These equations
ip t
558
often take the form
(17)
cr
(qjk1 t , . . . , qjk2 ,t ) = Ξjt (qjk3 t , . . . , qjk4 t ),
in which we explicitly identify which resources are the independent variables x and which are the dependent variables
561
y. Without loss of generality, we may reparameterize our functions so that all of the qjkt are outputs, i.e., we have
562
only relationships of the form
us
560
(18)
an
ˆ jt (z), (qj1t , qj2t , . . . , qjKt ) = Ξ
for a new dummy parameter z. For example, for f defined on q ∈ [Qmin , Qmax ], we can take z ∈ [0, 1] and add the
564
ˆ We then define the set N to index the mapping q = Qmin + z(Qmax − Qmin ) to Ξ, yielding the reparameterized Ξ.
565
extreme points zn in z-space and set
M
563
(ζjk1 nt , . . . , ζjkK nt ) = Ξjt (zn )
(19)
and leave z as the free optimization variable with q calculated according to (3). By using this general formulation, it
567
does not matter what the specific dependent and independent variables are for each generator.
568
Appendix C
569
To validate the piecewise-linearization approach, we can compare the optimal solution returned from the model MMILP
570
to results obtained from using full nonlinear equipment models and energy balances, which is denoted MMINLP . This
571
system is diagrammed in Figure 17. For simplicity, only one of each piece of equipment is shown. Variables and
572
parameters are described in the sections below. Functional forms of equipment models are shown in Section C.5.
573
C.1
574
Each chiller j has the following variables:
Ac ce p
te
d
566
Nonlinear Equipment Selection Formulation
Chillers
575
• Binary uCH j for its on/off state
576
CH • Continuous QCH for its cooling and electrical loads j and Wj
577
• Continuous VjCHW and VjCHW for its volumetric flows of chilled and cooling water 29
Page 29 of 38
Cooling Tower mA l WlCT
Tamb φamb
QCT l
T CWR
ip t
mCW l T CWS
QCH j
T CHWR
Pump
VjCHW Chiller
WkP
TjCHWS
us
WjCH
cr
VjCW
• Continuous TjCHWS for the supply temperature of chilled water
M
578
an
Figure 17: Diagram of chiller/cooling tower/pump system.
Parameters are c and ρ, the heat capacity and density of water, and T CHWR for the chilled water return temperature.
580
Equations are as follows:
d
579
te
CH CHWS WjCH = uCH , T CWS ) j fchiller (Qj , Tj CH CH CH CH uCH j Qmin ≤ Qj ≤ uj Qmax
Ac ce p
CW CW CW uCH ≤ uCH j Vmin ≤ Vj j Vmax
CHW QCH (T CHWR − TjCHWS ) j = ρcVj
581
The first equation calculates electric load using the nonlinear model fchiller (see (20) in Section C.5). Note that T CWS
582
is the supply temperature of cooling water determined by cooling tower operation. The next two equations enforce
583
minimum and maximum capacities. The final equation calculates the volumetric flow of chilled water using an energy
584
balance.
585
C.2
586
Each pump k has the following variables:
Pumps
587
• Binary uPk for its on/off state
588
• Continuous VlP and VlP for its volumetric flow of chilled water
30
Page 30 of 38
589
Equations are as follows: WkP ≥ uPk fpump (VkP )
ip t
P P uPk Vmin ≤ VkP ≤ uCH j Vmax
The first equation gives electricity consumption as a function of volumetric flow, and the second equation enforces
591
minimum and maximum capacities.
592
C.3
593
Each cooling tower l has the following variables:
cr
590
us
Cooling Towers
• Binary uCH j for its on/off state
595
CT • Continuous QCT for its cooling and electrical loads l and Wl
596
• Continuous mCW and mA l for its mass flows of cooling water and air l
M
597
an
594
Equations are as follows:
d
CW A CWR QCT ) l = ftower (ml , ml , T
te
WlCT = ffan (mA l )
Ac ce p
CT CT CT CT uCT l Qmin ≤ Ql ≤ ul Qmax
598
The first equation calculates cooling duty as a function of mass flows and cooling water return temperature. Note that
599
because of the model form, we need not calculate individual values for T CWS . The second equation uses the tower fan
600
model to give electric load, while the third equation enforces minimum and maximum capacities.
31
Page 31 of 38
601
C.4
Balances
602
Finally, various material balances must be enforced to ensure that sufficient equipment is switched on. These are as
603
follows: Total QCH j ≥Q
ip t
X j
X
VjCHW =
X
j
VjCW =
X
j
X
mCW l
l
CH QCH = j + Wj
j
X
CWR mCW − T CWS ) l c(T
l
W Total =
X
cr
k
X
WjCH +
X
X
WlCT
an
j
WkP +
us
ρ
VkP
k
l
The first balance ensures that enough chilled water is produced to meed the prespecified cooling duty QTotal . The
605
second balance ensures that enough pumping is available for chilled water. The third balance manages cooling water
606
flow among the chillers and cooling towers. The fourth balance is used to calculate the aggregate cooling water supply
607
temperature T CWS , which affects the electricity consumption in the chillers. Finally, the fifth balance calculates the
608
total electricity consumption W Total , which serves as the objective to be minimized.
609
C.5
610
For modeling of chillers, we use the semi-empirical Gordon-Ng model, [41] defined below:
Ac ce p
Equipment Models
te
d
M
604
fchiller (Q, T CHWS , T CWS ) :=
T CWS T CHWS −Q Q + a1 T CHWS + a2 1 − CWS T T CHWS − a3 Q
(20)
611
The parameters a1 , a2 , and a3 are obtained via regression with measured data, and they are related to entropy changes
612
within the process.
613
For pumps, the following black-box empirical model is used:
fpump (V ) := b1 ln (1 + b2 V ) + b3 V + b4
614
615
Here, b1 through b4 are regression coefficients. Finally, for cooling towers, a simplified effectiveness model [46] is used for calculating cooling duty, while for fan
32
Page 32 of 38
v+ , Hhot /Vhot
Hhot Vhot
Tank Discharging
v+ , h+
v− , Hcold /Vcold
cr
Hcold Vcold
ip t
Tank Charging
v− , h−
electricity consumption, a simple cubic fit is used [20]. ftower (mCW , mA , T CWR ) :=
an
616
us
Figure 18: Diagram of stratified tank model. Tank is treated as two separate volumes of water that interact with each other and the ambient. Only two of the four streams are active at a given time depending on whether the tank is charging or discharging.
c1 (mCW )c3 c3 (T CWR − Tamb ) mCW 1 + c2 m A
M
ffan (mA ) := κ(mA )3 c1 , c2 , c3 , and κ are obtained via regression.
618
Appendix D
619
The storage model (2) is an approximation of a stratified tank model [47]. Within this section, we use v+ to denote
620
the volumetric flow for when the tank is charging, and v− for when the tank is being discharged. We assume the tank
621
volume is held constant, and thus model it as consisting of two stratified regions, each with separate volume V and
622
enthalpy H. This is illustrated in Figure 18.
te
d
617
Ac ce p
Tank Model Derivation
623
Within the tank, here is energy exchange between the two layers of the tank, and there is also external exchange
624
between each layer and the ambient. Often the tank is assumed to be well-insulated, which means the only enthalpic
625
change to the tank is due to heat exchange between the two layers (i.e., k = 0 in the model). However, we retain the
33
Page 33 of 38
ambient term for now and arrive at the following mass and energy balances.
627
= −v+ + v− = v+ − v− Hhot Hhot Hcold Hhot =− v+ + h − v− − K − −k − hamb Vhot Vhot Vcold Vhot Hcold Hhot Hcold Hcold = h+ v+ − v− + K − −k − hamb Vcold Vhot Vcold Vcold
cr
dVhot dt dVcold dt dHhot dt dHcold dt
ip t
626
in which parameters h+ and h− are the (per-volume) supply and return enthalpies.
0 0 0 0 0 0 To simplify the dynamics, we choose an operating point (Vhot , Vcold , Hhot , Hcold ) such that Hhot /Vhot = h− and
629
0 0 Hcold /Vcold = h+ (i.e., that the extra “hot” water at the top of the tank is the same temperature as the “spent” chilled
630
water returning from the building). Linearizing about this point yields the following energy balances:
an
us
628
M
Vhot Vcold dHhot Hhot Hcold ≈ −h− v+ + h− v− − (K + k)h− 1 + 0 − 0 + Kh+ 1 + 0 − 0 + khamb dt Hhot Vhot Hcold Vcold dHcold Hhot Vhot Vcold Hcold ≈ h+ v+ − h+ v− + Kh− 1 + 0 − 0 − (K + k)h+ 1 + 0 − 0 + khamb dt Hhot Vhot Hcold Vcold Then, noting that we are free to choose a reference point for the enthalpy terms, we set h− = 0 and substitute
632
0 0 h+ = Hcold /Vcold in the second term, which leaves the remaining cold enthalpy balance
te
d
631
633
Ac ce p
H0 dHcold = h+ (v+ − v− ) − (K + k) cold 0 dt Vcold K ≈ h+ (v+ − v− ) − 0 Hcold Vcold
1+
Hcold Vcold − 0 0 Hcold Vcold
+ khamb
0 assuming Vcold ≈ Vcold and k ≈ 0. Making the change of variables s := Hcold and y := h+ (v+ − v− ), we arrive at
ds = −˜ σs + y dt
(21)
634
0 with σ ˜ := K/Vcold . After discretization, (21) becomes the storage model (2) used within the optimization problem.
635
References
636
[1] M. Fadzli Haniff, H. Selamat, R. Yusof, S. Buyamin, F. Sham Ismail, Review of HVAC scheduling techniques
637
for buildings towards energy-efficient and cost-effective operations, Renewable and Sustainable Energy Reviews
638
27 (2013) 94–103.
34
Page 34 of 38
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
ip t
cr
645
storage, Journal of Process Control 24 (8) (2014) 1292–1300.
[5] M. H. Albadi, E. F. El-Saadany, Demand response in electricity markets: An overview, in: 2007 IEEE Power Engineering Society General Meeting, IEEE, 2007, pp. 1–5.
us
644
[4] C. R. Touretzky, M. Baldea, Integrating scheduling and control for economic MPC of buildings with energy
[6] G. P. Henze, Energy and cost minimal control of active and passive building thermal storage inventory, Journal of Solar Energy Engineering 127 (3) (2005) 343–351.
an
643
through energy optimization, Processes 1 (3) (2013) 312–329.
[7] M. Kintner-Meyer, A. F. Emery, Optimal control of an HVAC system using cold storage and building thermal capacitance, Energy and Buildings 23 (1) (1995) 19–31.
M
642
[3] K. Kapoor, K. M. Powell, W. J. Cole, J. S. Kim, T. F. Edgar, Improved large-scale process cooling operation
[8] G. P. Henze, C. Felsmann, G. Knabe, Evaluation of optimal control for active and passive building thermal storage, International Journal of Thermal Sciences 43 (2) (2004) 173–183.
d
641
energy planning, Energy and Buildings 58 (2013) 281–291.
[9] M. Behl, T. X. Nghiem, R. Mangharam, Green scheduling for energy-efficient operation of multiple chiller plants, Real-Time Systems Symposium (RTSS), 2012 IEEE 33rd (2012) 195–204.
te
640
[2] L. Duanmu, Z. Wang, Z. J. Zhai, X. Li, A simplified method to predict hourly building cooling load for urban
[10] H. Dagdougui, R. Minciardi, A. Ouammi, M. Robba, R. Sacile, Modeling and optimization of a hybrid system
Ac ce p
639
for the energy supply of a green building, Energy Conversion and Management 64 (2012) 351–363. [11] A. Ashouri, S. S. Fux, M. J. Benz, L. Guzzella, Optimal design and operation of building services using mixedinteger linear programming techniques, Energy 59 (2013) 365–376. [12] M. Trifkovic, W. A. Marvin, P. Daoutidis, M. Sheikhzadeh, Dynamic real-time optimization and control of a hybrid energy system, AIChE Journal 60 (7) (2014) 2546–2556. [13] S.-I. Gustafsson, Mixed integer linear programming and building retrofits, Energy and Buildings 28 (2) (1998) 191–196. [14] D. Zhang, N. Shah, L. G. Papageorgiou, Efficient energy consumption and operation management in a smart building with microgrid, Energy Conversion and Management 74 (2013) 209–222. [15] Y. Ma, F. Borrelli, B. Hencey, B. Coffey, S. Bengea, P. Haves, Model predictive control for the operation of building cooling systems, IEEE Control Systems Technology 20 (3) (2012) 796–803. 35
Page 35 of 38
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
ip t
cr
673
Energy and Buildings 110 (2016) 94–107.
[19] J. E. Braun, A near-optimal control strategy for cool storage systems with dynamic electric rates, HVAC&R Research 13 (4) (2007) 557–580.
us
672
[18] C. R. Touretzky, M. Baldea, A hierarchical scheduling and control strategy for thermal energy storage systems,
[20] J. E. Braun, G. T. Diderrich, Near-optimal control of cooling towers for chilled-water systems, ASHRAE Transactions 96 (CONF-9006117).
an
671
I: Simulation environment, Energy Conversion and Management 64 (2012) 499–508.
[21] J. E. Braun, A general control algorithm for cooling towers in cooling plants with electric and/or gas-driven chillers., HVAC&R Research 13 (4) (2007) 581–598.
M
670
[17] A. Hajiah, M. Krarti, Optimal control of building storage systems using both ice storage and thermal mass—part
[22] M. Ali, V. Vukovic, M. H. Sahir, G. Fontanella, Energy analysis of chilled water system configurations using simulation-based optimization, Energy and Buildings 59 (2013) 111–122.
d
669
building cooling system, Applied Energy 111 (2013) 1032–1045.
[23] A. Aravelli, S. S. Rao, Energy optimization in chiller plants: A novel formulation and solution using a hybrid optimization technique, Engineering Optimization 45 (10) (2013) 1187–1203.
te
668
[16] J. A. Candanedo, V. R. Dehkordi, M. Stylianou, Model-based predictive control of an ice storage device in a
[24] A. J. Ardakani, F. F. Ardakani, S. H. Hosseinian, A novel approach for optimal chiller loading using particle
Ac ce p
667
swarm optimization, Energy and Buildings 40 (12) (2008) 2177–2187. [25] W.-S. Lee, Y.-T. Chen, Y. Kao, Optimal chiller loading by differential evolution algorithm for reducing energy consumption, Energy and Buildings 43 (23) (2011) 599–604. [26] Z. W. Geem, Solution quality improvement in chiller loading optimization, Applied Thermal Engineering 31 (10) (2011) 1848–1851.
[27] L. Lu, W. Cai, Y. C. Soh, L. Xie, S. Li, HVAC system optimization—condenser water loop, Energy Conversion and Management 45 (4) (2004) 613–630. [28] K. M. Powell, W. J. Cole, U. F. Ekarika, T. F. Edgar, Optimal chiller loading in a district cooling system with thermal energy storage, Energy 50 (2013) 445–453. [29] G. P. Henze, B. Biffar, D. Kohn, M. P. Becker, Optimal design and operation of a thermal storage system for a chilled water plant serving pharmaceutical buildings, Energy and Buildings 40 (6) (2008) 1004–1019. 36
Page 36 of 38
699
700
701
702
703
704
electricity pricing, Energy and Buildings 60 (2013) 199–209.
ip t
698
[31] M. Avci, M. Erkoc, A. Rahmani, S. Asfour, Model predictive HVAC load control in buildings using real-time
[32] J. Ma, J. Qin, T. Salsbury, P. Xu, Demand reduction in building energy systems based on economic model predictive control, Chemical Engineering Science 67 (1) (2012) 92–100.
cr
697
in: American Control Conference (ACC), 2012, 2012, pp. 3669–3674.
[33] R. Z. Freire, G. H. Oliveira, N. Mendes, Predictive controllers for thermal comfort optimization and energy savings, Energy and Buildings 40 (7) (2008) 1353–1365.
us
696
[30] D. I. Mendoza-Serrano, D. J. Chmielewski, Controller and system design for HVAC with thermal energy storage,
[34] S. C. Bengea, A. D. Kelman, F. Borrelli, R. Taylor, S. Narayanan, Implementation of model predictive control for an HVAC system in a mid-size commercial building, HVAC&R Research 20 (1) (2014) 121–135.
an
695
[35] M. J. Risbeck, C. T. Maravelias, J. B. Rawlings, R. D. Turney, Cost optimization of combined building heat-
706
ing/cooling equipment via mixed-integer linear programming, in: American Control Conference, Chicago, IL,
707
2015, pp. 1689–1694.
710
711
d
MILP formulation, Computers & Chemical Engineering 17 (2) (1993) 211–227. [37] C. C. Pantelides, Unified frameworks for optimal process planning and scheduling, in: Proceedings on the second
te
709
[36] E. Kondili, C. Pantelides, R. Sargent, A general algorithm for short-term scheduling of batch operations—I.
conference on foundations of computer aided operations, 1994, pp. 253–274.
Ac ce p
708
M
705
712
[38] C. A. M´endez, J. Cerd´a, I. E. Grossmann, I. Harjunkoski, M. Fahl, State-of-the-art review of optimization meth-
713
ods for short-term scheduling of batch processes, Computers & Chemical Engineering 30 (6-7) (2006) 913–946.
714
[39] C. T. Maravelias, General framework and modeling approach classification for chemical production scheduling,
715
AIChE Journal 58 (6) (2012) 1812–1828.
716
[40] I. Harjunkoski, C. T. Maravelias, P. Bongers, P. M. Castro, S. Engell, I. E. Grossmann, J. Hooker, C. M´endez,
717
G. Sand, J. Wassick, Scope for industrial applications of production scheduling models and solution methods,
718
Computers & Chemical Engineering 62 (2014) 161–193.
719
720
721
722
[41] T.-S. Lee, K.-Y. Liao, W.-C. Lu, Evaluation of the suitability of empirically-based models for predicting energy performance of centrifugal water chillers with variable chilled water flow, Applied Energy 93 (2012) 583–595. [42] W. J. Cole, T. F. Edgar, A. Novoselac, Use of model predictive control to enhance the flexibility of thermal energy storage cooling systems, in: American Control Conference (ACC), 2012, 2012, pp. 2788–2793. 37
Page 37 of 38
726
727
728
729
730
[44] D. I. Mendoza-Serrano, D. J. Chmielewski, HVAC control using infinite-horizon economic MPC, in: Decision and Control (CDC), 2012 IEEE 51st Annual Conference on, 2012, pp. 6963–6968.
ip t
725
Energy and Buildings 33 (4) (2001) 371–378.
[45] Z. Ma, S. Wang, X. Xu, F. Xiao, A supervisory control strategy for building cooling water systems for practical and real time applications, Energy Conversion and Management 49 (8) (2008) 2324–2336.
cr
724
[43] B. C. Ahn, J. W. Mitchell, Optimal control development for chilled water plants using a quadratic representation,
[46] G.-Y. Jin, W.-J. Cai, L. Lu, E. L. Lee, A. Chiang, A simplified modeling of mechanical cooling tower for control and optimization of HVAC systems, Energy Conversion and Management 48 (2) (2007) 355–365.
us
723
[47] Y. Ma, F. Borrelli, B. Hencey, A. Packard, S. Bortoff, Model predictive control of thermal energy storage in build-
732
ing cooling systems, in: Decision and Control, 2009 held jointly with the 2009 28th Chinese Control Conference.
733
CDC/CCC 2009. Proceedings of the 48th IEEE Conference on, 2009, pp. 392–397.
737
738
739
740
741
742
743
744
745
M
ˇ acˇ ekov´a, Building modeling as a crucial part [49] S. Pr´ıvara, J. Cigler, Z. V´anˇ a, F. Oldewurtel, C. Sagerschnig, E. Z´ for building predictive control, Energy and Buildings 56 (2013) 8–22.
d
736
Energy and Buildings 108 (2015) 454–462.
[50] A. Imamoto, B. Tang, Optimal piecewise linear approximation of convex functions, in: Proceedings of the world
te
735
[48] K. J. Kircher, K. Max Zhang, On the lumped capacitance approximation accuracy in rc network building models,
congress on engineering and computer science, 2008, pp. 1191–1194.
Ac ce p
734
an
731
[51] C. Frenzen, T. Sasao, J. T. Butler, On the number of segments needed in a piecewise linear approximation, Journal of Computational and Applied Mathematics 234 (2) (2010) 437–446. [52] S. Rippa, Adaptive approximation by piecewise linear polynomials on triangulations of subsets of scattered data, SIAM Journal on Scientific and Statistical Computing 13 (5) (1992) 1123–1141. [53] A. Toriello, J. P. Vielma, Fitting piecewise linear continuous functions, European Journal of Operational Research 219 (1) (2012) 86–95.
38
Page 38 of 38
S0378-7788(17)30647-3 http://dx.doi.org/doi:10.1016/j.enbuild.2017.02.053 ENB 7414
To appear in:
ENB
Received date: Revised date: Accepted date:
23-3-2016 7-10-2016 23-2-2017
Please cite this article as: Michael J. Risbeck, Christos T. Maravelias, James B. Rawlings, Robert D. Turney, A Mixed-Integer Linear Programming Model for Real-Time Cost Optimization of Building Heating, Ventilation, and Air Conditioning Equipment, (2017), http://dx.doi.org/10.1016/j.enbuild.2017.02.053 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A Mixed-Integer Linear Programming Model for Real-Time Cost
cr
Conditioning Equipment
ip t
Optimization of Building Heating, Ventilation, and Air
an
October 7, 2016
us
Michael J. Risbeck,a,∗ Christos T. Maravelias,b James B. Rawlings,c Robert D. Turneyd
Abstract
M
In this paper, we present a framework for the formulation and solution of mixed-integer linear programming (MILP) models for operational planning of HVAC systems in commercial buildings. We introduce the general concepts of generators (e.g., chillers, boilers, cooling towers) and resources (e.g., electricity, chilled water), which allow
d
us to model a wide range of central plants. By using discrete variables for on/off states of central plant equipment and
te
continuous variables for equipment load, storage tank usage, and building temperature trajectory, all critical-to-cost decisions are taken into account. In addition, we consider both time-varying use charges and peak demand charges as well as building models of varying complexity. Because equipment models are often nonlinear, piecewise-linear
Ac ce p
approximations are used, which can be made arbitrarily accurate and enable real-time solution to the resulting optimization problems. We also employ a symmetry-free formulation to enhance the solution of the MILP model. Such features lead to improved performance compared to approaches employing heuristics or optimization without discrete variables. We demonstrate optimization for a small-scale cooling system and show favorable scaling of solution time with respect to number of units and temperature zones. In addition, we show simultaneous optimization of heating and cooling loops when a forecast of resource demand is available. Finally, we demonstrate closed-loop operation of our proposed optimization scheme. In all cases examined, a solution with an optimality gap below 1% can be obtained within 5 minutes, and thus the proposed optimization can be solved in real time, allowing for rapid response to changing weather, price, or demand forecasts. Keywords: online optimization; central plants; cost reduction; thermal energy storage; scheduling ∗Corresponding Author a University of Wisconsin–Madison.
Page 1 of 38
1
Introduction
2
Heating, ventilation, and air conditioning (HVAC) equipment consumes a significant amount of energy within com-
3
mercial buildings. Such systems use (often simple) heuristic methods for scheduling and control [1], which result in
4
lower system efficiency. In addition, time-varying occupancy and ambient conditions lead to nonuniform electricity
5
usage [2–4], which causes utility companies to implement alternative pricing structures that strive to reduce peak elec-
6
tricity demand by consumers [5] because electricity production is most efficient when demand is relatively constant.
7
Such pricing structures include time-of-use charges (assessed using time-varying electricity rates) and peak demand
8
charges (assessed based on a customer’s one-time maximum rate of electricity consumption). With thermal energy
9
storage (TES) that can temporally shift electricity use, informed customers can leverage these pricing structures to
10
reduce total cost. This storage comes in the form of both active TES, such as purpose-built water or ice storage tanks,
11
and passive TES, which takes advantage of the natural thermal capacitance of buildings [6]. A combination of active
12
and passive storage techniques is often the most efficient solution [7, 8], but using these techniques requires a nontrivial
13
decision-making process [9].
an
us
cr
ip t
1
For these systems, a key challenge in planning using math programming methods is that there are both continuous
15
decision variables for equipment setpoints, utility consumption, storage usage, etc., and discrete decision variables
16
to specify whether process equipment is switched on or off. This complexity can be captured by mixed-integer pro-
17
gramming (MIP) models, which have been employed to make other building efficiency and sustainability decisions,
18
including selection of renewable and non-renewable energy sources [10–12], optimizing building retrofits based on
19
life-cycle analysis [13], and scheduling for “smart” appliances connected to a microgrid [14]. MIP models are most
20
computationally tractible when they have only linear constraints and objective (mixed-integer linear programming,
21
MILP), but more general mixed-integer nonlinear programming (MINLP) models have also been developed. Such
22
models can consider process energy balances in greater detail [15–18], but they are often limited to systems with few
23
pieces of equipment or few degrees of freedom.
Ac ce p
te
d
M
14
24
Because large-scale MIP models are computationally demanding, decomposition strategies are often employed. A
25
typical architecture is to is to make equipment selection choices in a small subproblem, either through near-optimal
26
heuristics [19–21], or via solving single-period scheduling problem [22–27]. Results of subproblems can then be used
27
to construct aggregate system curves [28], which simplify optimization at the higher level by eliminating all the dis-
28
crete, on/off choices for equipment. For asymmetric models of TES (i.e., models with different dynamics depending
29
on whether the system is charging or discharging), discrete choices must also be made for TES operation. Meth-
30
ods such as backward reach set computation [9] and dynamic programming [29] have been used to precompute the
31
charge/discharge profile, which also serves to simplify optimization at other levels. However, an important drawback
32
of such decomposition schemes is that equipment on/off states are often not considered over a long horizon, which
2
Page 2 of 38
means the effects of switching are not completely considered. Thus, units may be switched on and off rapidly in
34
response to small variations in system load, which must be addressed via further ad-hoc adjustments to the top-level
35
optimization problem. In addition, most approaches in the literature consider only chilled water production in a single
36
chiller, whereas the central plants we wish to address contain multiple units of different types and must produce hot
37
and chilled water concurrently.
ip t
33
In some cases, it may be desirable to determine demand by optimizing the temperature trajectory within buildings
39
using a model instead of simply following a fixed demand forecast. Building models can be either derived from first
40
principles [30] or identified from data [31]. Such models often split buildings into multiple interacting temperature
41
zones with separate control inputs, requiring the temperature trajectory of each zone to be considered. While buildings
42
often exhibit nonlinear dynamics, simplified linear building models can be sufficient for control [32]. The most
43
important feature of these models is that they relate occupant comfort to zone heating or cooling. Constraints on
44
occupant comfort are often implemented in terms of temperature, although more complicated metrics are possible
45
[33]. Studies have shown that such model-based control can out-perform conventional methods in terms of cost and
46
energy savings [34]. However, many existing approaches do not consider the complexities of the central plant (such
47
as on/off decisions or active storage) by assuming a constant coefficient of performance for cooling equipment. Thus,
48
a complete treatment must account for both buildings and the central plant.
M
an
us
cr
38
Accordingly, the goal of this paper is to propose an optimization-based framework that treats in detail both the
50
building and central plant subsystems. We propose a mathematical programming model that considers multiple pieces
51
of central plant equipment (for both heating and cooling) and multiple temperature zones over a finite prediction
52
horizon. In [35], an MILP formulation for the central plant subsystem was presented. This formulation approximates
53
critical-to-cost nonlinearities as piecewise-linear, which allows powerful MILP software to be used for optimization.
54
Improved performance compared to formulations without discrete decision variables was then demonstrated. In this
55
work, we extend the previous formulation to include a building model that allows optimal demand to be determined
56
online. In addition, we use a symmetry-free reformulation for identical pieces of equipment (not present in [35]) so
57
that explicit symmetry-breaking constraints are no longer necessary. These augmentations provide wider applicability
58
and more favorable scaling (with respect to horizon length and to number of central plant units) compared to the
59
previous formulation and also to alternative nonlinear formulations.
Ac ce p
te
d
49
60
The remainder of this paper is structured as follows. In Section 2, we present the problem statement and complete
61
MILP formulation (subsuming the formulation from [35]). In Section 3, we discuss the equipment models required
62
by the formulation. In Section 4, we present a series of simulations to validate the assumptions made by the model.
63
In particular, we discuss solution times for the combined building/central plant formulation; examine the cost savings
64
for a combined heating/cooling plant; address the handling of variable supply temperatures by benchmarking against
65
a more detailed (but computationally challenging) nonlinear formulation; and, demonstrate closed-loop operation as 3
Page 3 of 38
66
would be required in a real implementation. Finally, Section 5 summarizes the results of the study.
67
2
68
2.1
69
The overall objective is to maintain comfort of building occupants at the lowest economic cost. Depending on climate,
70
the focus may be on producing chilled water for building cooling in chillers or on producing hot water for building
71
heating in boilers. Nevertheless, hot and chilled water can be used for purposes unrelated to the HVAC system, and
72
thus simultaneous production of both resources must be considered. In these cases, heat-recovery chillers (HRCs) can
73
produce both hot and chilled water at higher efficiencies than separate units, as waste heat is rejected to a process hot
74
water stream instead of to the ambient. Thus, low-cost operation requires use of all types of equipment. An additional
75
challenge is that building temperature is often nonuniform with multiple interacting regions, each with their own
76
control inputs. Thus, a complete treatment of the problem requires simultaneous optimization of hot water production,
77
chilled water production, and building temperature trajectory.
ip t
cr
us
an
M
79
Problem Statement
The overall problem is as follows: given forecasts of hourly utility prices, ambient conditions and building occupancy, decide how to operate the HVAC system at lowest economic cost via the following decisions:
d
78
Optimization Problem
• What is the temperature profile of each zone?
81
• When is storage charging or discharging?
82
• Is each piece of equipment running?
83
• How much chilled and hot water is being producing?
Ac ce p
te
80
84
While making these decisions, temperature comfort constants must be satisfied, equipment minimum and maximum
85
capacities must be respected, and switching of equipment must be limited.
86
For simplicity, we assume that the supply and return temperatures of the chilled and hot water are known ahead
87
of time and fixed. While these decisions could be made within the optimization, it does not appear that variable
88
temperatures result in significant cost savings (See Section 4.3 for more details). In addition, cooling water supply
89
temperature is treated as a fixed parameter, although its value need not be constant throughout the prediction horizon.
90
Finally, we assume that building temperature dynamics can be sufficiently approximated using a linear model, or that
91
a forecast of building hot and chilled water demand is available.
4
Page 4 of 38
ip t
cr
Figure 1: Diagram of abstract optimization model. Within the load, primary demand is represented by the red, blue, and green curves at the top, and secondary demand is represented by the orange curve at the bottom.
2.2
Mathematical Formulation
93
Instead of formulating separate, individual constraints and variables for each type of equipment, we use an abstract
94
representation of the central plant in terms of “resources” and “generators.” This abstraction is similar to the state-task
95
network [36] and resource-task network [37] models that have been employed to great success for chemical production
96
scheduling [38–40]. Within the abstract model, a generator j ∈ J represents any piece of equipment in the central
97
plant (e.g., chillers, HRCs, pumps), while a resource k ∈ K represents any material flow in or out of the plant (e.g.,
98
chilled water, electricity). The prediction horizon is split into discrete intervals of equal length ∆ that are indexed by
99
t ∈ T. This overall abstract system is illustrated in Figure 1. In the abstract model, the system purchases resources
100
from the market, uses generators to convert these to other resources, and finally delivers these resources to the load,
101
optionally employing storage along the way. Once at the load, the resources are used to control the temperature
102
dynamics of zones i ∈ I or to serve auxiliary processes.
Ac ce p
te
d
M
an
us
92
103
We consider two different types of demand: (1) “primary demand,” which is is a decision variable and is used
104
to affect zone temperature (e.g., chilled water sent to air handlers); and, (2) “secondary demand,” which is a fixed
105
parameter that accounts for resource use by auxiliary processes (e.g., hot water for domestic use or chilled water for
106
specific internal cooling needs). Note that in different problem instances, demand may be entirely primary or entirely
107
secondary.
108
In the following sections, we use lower-case Roman letters for decision variables and subscripts, Greek letters for
109
fixed parameters, and bold upper-case Roman letters for sets. A complete listing of all symbols is given in Appendix A.
110
2.2.1
111
For a resource k and time t, primary demand in zone i is modeled as a variable gikt , while secondary demand φkt is
112
assumed to be known throughout the prediction period. Satisfaction of both types of demand is ensured through the
Demand Constraint
5
Page 5 of 38
113
following constraint: X
qjkt + ykt + pkt + bkt ≥ φkt +
j∈J
X
gikt ,
k ∈ K, t ∈ T.
(1)
i∈I
Thus, total demand is met via production qjkt , net storage withdrawal ykt , and direct purchase pkt . To soften this
115
constraint, the slack variable bkt ≤ φkt can be used, albeit at a heavy penalty in the objective function. By softening
116
the secondary demand constraint, the optimizer can still find a feasible solution even if secondary demand cannot
117
possibly be met (e.g., due to unit breakdowns). Note that primary demand cannot be ignored here; instead, temperature
118
comfort constraints are softened as described later.
cr
Stored resource inventory skt evolves based on its charge or discharge rate ykt according to the following dynam-
119
ics: skt = σk sk(t−1) + ykt ,
us
120
ip t
114
k ∈ K, t ∈ T,
(2)
where the term σk ≤ 1 is included to model storage inefficiencies, e.g., addition of heat to a chilled water tank from the
122
ambient. While TES dynamics are most accurately described by more complex effectiveness models with temperature
123
dependence [16], simpler state-of-charge models such as (2) have been shown to provide sufficient accuracy [6].
124
2.2.2
125
The true production/consumption relationships in the generators are often nonlinear, which requires enforcing con-
126
straints of the form
M
an
121
te
d
Equipment Models
Ξjt (ujt , qjk1 t , qjk2 t , . . . , qjkK t ) = 0,
that is, a nonlinear relationship between the production or consumption qjkt of the various resources, with additional
128
dependence on the on/off state uj of the generator. Note that these functions may be time-varying due to dependence
129
on external time-varying parameters (e.g., ambient temperature or humidity). For computational tractability, we wish
130
to solve the optimization as an MILP, and thus, we must use piecewise-linear approximations of the true nonlinear
131
model.
Ac ce p
127
132
For each type of generator j, we partition the (bounded) domain of Ξjt into a set Mj of polyhedrons, and we also
133
build the set Nmj of extreme points for each polyhedron m. For each extreme point n in each region m, a continuous
134
variable zjmnt is declared to choose the weighting of the given point. With these definitions, generator production and
135
consumption rates are calculated by
qjkt =
X
X
ζjknt zjmnt ,
j ∈ J, k ∈ K, t ∈ T,
(3)
m∈Mj n∈Njm
136
giving qjkt as a convex combination of the operating points stored in the parameter ζjknt .
6
Page 6 of 38
137
138
To decide how many generators of a given type are on, we use integer variables ujt , and to choose from among the available subdomains, we declare integer variables vjmt . Then, the following constraints are used: X
vjmt = ujt ,
(4)
j ∈ J, t ∈ T,
X
j ∈ J, m ∈ Mj , t ∈ T.
zjmnt = vjmt ,
n∈Njm
ip t
m∈Mj
(5)
Equation (4) chooses exactly one subdomain for each generator of type j that is on, and (5) distributes weight to the
140
chosen subdomains’ extreme points. Implicit in these constraints is that whenever ujt = 0, all corresponding qjkt = 0,
141
i.e., that when all generators are off, there is no production or consumption. By using integer-valued variables ujt and
142
vjmt for identical generators, we avoid computationally costly symmetry, which would be present if multiple sets
143
binary variables were used. Details of how the piecewise linearizations are constructed are presented in Appendix B.
144
The functional forms Ξjt of the true nonlinear models are presented in Section C.5 in Appendix C.
145
2.2.3
146
With integer variables to count the number of type-j generators that are on, rapid switching can be restricted via
147
constraints. First, the number of switching events is calculated using
an
us
cr
139
d
M
Switching Constraints
te
− ujt − uj(t−1) = u+ jt − ujt ,
j ∈ J, t ∈ T,
(6)
− where u+ jt and ujt count the number of switch-on and switch-off events for the given period. With these values, dwell
149
times are enforced as follows:
Ac ce p
148
ujt ≥
δj+ −1
X
u+ j(t−τ ) ,
j ∈ J, t ∈ T,
(7)
τ =0
µj − ujt ≥
δj− −1
X
u− j(t−τ ) ,
j ∈ J, t ∈ T,
(8)
τ =0
150
in which µj is the total number of generators of type j. Equation (7) requires that the number of “on” generators is
151
at least equal to the number of “switch on” events in the last δj+ time periods, while (8) enforces the corresponding
152
condition for “off” generators and “switch off” events.
153
2.2.4
154
Within large buildings, there are typically multiple independent (but interacting) temperature zones. Letting fit denote
155
the temperature in zone i at the end of time period t (i.e., at absolute time t∆, with ∆ the discrete timestep), and
Building Model
7
Page 7 of 38
156
declaring gikt to be the flow of resource k (e.g., cold water or hot water) into zone i during period t. With these
157
variables, temperature evolves according to ! fit =
X
λii0 fi(t−1) +
i0 ∈I
X
ωii0 k gi0 kt
+ θit ,
i ∈ I, t ∈ T.
(9)
ip t
k∈K
In this constraint, the coefficients λii0 and ωii0 k describe a linear time-invariant model for temperature evolution, while
159
θit accounts for time-varying disturbances. Note that while the energy flow into zone i does not instantaneously affect
160
the temperature of zone i0 , there is fill-in when the model is converted from continuous to discrete time, and thus ωii0 k
161
requires the double index ii0 . Details of this model can be found in Section 3.2.
cr
158
Since the objective is to maintain temperature bounds at minimum cost, constraints need to be introduced to enforce
163
temperature bounds. To avoid infeasibility, these constraints are softened by adding nonnegative slack variables fit+
164
and fit− as follows:
an
us
162
U fit ≤ ψit + fit+ ,
(10)
i ∈ I, t ∈ T,
(11)
M
L fit ≥ ψit − fit− ,
i ∈ I, t ∈ T,
L U with ψit and ψit the lower and upper comfort bounds. This constraint softening is to allow primary demand gikt
166
to be reduced without making the solution infeasible. Together with the slack bkt for secondary demand, softened
167
constraints ensure that the model has a feasible solution regardless of the demand and ambient temperature forecasts.
168
2.2.5
169
The objective function minimizes total resource cost ctot k , with additional penalty terms of soft constraint slack vari-
170
ables bkt , fit+ , and fit− as follows:
te
d
165
Ac ce p
Objective Function
min
X
k∈K
ctot k
+
! X
βkt bkt
+
t∈T
XX
+ − − χ+ it fit + χit fit .
(12)
t∈T i∈I
171
− Typically, the coefficients βkt , χ+ it and χit are be chosen quite large to ensure that the corresponding soft constraints
172
are satisfied whenever possible.
173
Economic cost ctot k considers both time-of-use and demand charges, calculated according to max max ctot pk + k = ρk
X
ρkt pkt ,
∀k ∈ K,
(13)
t∈T
pkt ≤ pmax , k
∀k ∈ K, t ∈ T.
(14)
8
Page 8 of 38
Here, ρkt and ρmax are time-of-use and demand charge prices (assumed known), while pkt and pmax k kt are the decision
175
variables representing how much and when a given resource is purchased. If peak demand charge is considered
176
only for certain times throughout the day (e.g., designated “peak hours”), constraint (14) need only be included for
177
the appropriate subset of periods. Multiple demand charges can also be considered by introducing additional pmax k
178
variables and constraints (14) for each demand “window.” For demand charges assessed over times longer than the
179
prediction horizon, pmax can take a lower bound Πmax equal to the prior maximum. Additionally, the weight ρmax of k k k
180
the demand charge can be increased throughout the period to reflect its greater importance near the end of the billing
181
cycle.
182
3
183
The optimization problem with constraints (1) through (11), (13), and (14), minimizing the objective function (12) is
184
referred to as MMILP . In the following subsections, we discuss the models used by the formulation and the generality
185
afforded by the abstract representation.
186
3.1
187
To ensure that the optimization problem correctly calculates total operating cost, accurate models of the generators’
188
production/consumption relationships are needed. Although these models are typically nonlinear, by replacing them
189
with piecewise-linear approximations as described in Appendix B, there are no restrictions on the algebraic form of
190
the models. However, a decision still must be made about what terms should be included as variables versus what
191
terms can be taken as parameters.
us
cr
ip t
174
M
an
Discussion
Ac ce p
te
d
Equipment Models
192
Numerous empirical models are available for chillers [41], often relating coefficient of performance (COP) as a
193
function of cooling load and one or more stream temperatures. For example, the Gordon-Ng model uses cooling
194
load, chilled water supply temperature, and cooling water supply temperature as independent variables [3]. For cost
195
optimization, the most relevant variable is cooling load, as it appears in the demand satisfaction constraint. To allow
196
for further optimization, supply temperature can also be considered as a variable (indirectly via volumetric flow; see
197
Section 4.3). Cooling water supply temperature, however, is mostly determined by ambient wet bulb temperature, and
198
it is thus treated as a (time-varying) parameter. Other models, including the Modified DOE-2 Model [42], empirical
199
polynomial fits, [6, 9, 43], and black-box lookup tables [15], use similar decision variables and can treated in the same
200
manner. It is also common to assume a constant COP [21, 29, 30, 44], which eliminates the need to consider any
201
stream temperatures. Thus, chiller models typically require at most two independent variables.
202
Because auxiliary equipment (e.g., pumps and cooling towers) comprises only a small part of total resource cost
203
and is often inflexible (i.e., has few or no degrees of freedom once other variables have been chosen), it is often 9
Page 9 of 38
ignored in optimization [42]. However, because auxiliary units affect peak demand charges, improved performance
205
(in terms of reduced cost) can be achieved by considering them even in approximate forms. For pumps, empirical
206
models of electricity consumption as a function of volumetric flow have been constructed from historical data [3, 27].
207
Cooling towers can be explicitly described by an effectiveness model [45, 46]. Alternatively, because cooling tower
208
performance depends on ambient temperature and humidity, these conditions can either be included in the chiller
209
power equation explicitly [15], or they can be used to estimate the optimal cooling water temperature in an empirical
210
function and then pass this value to chiller models [42]. In any case, choosing nonconvex or multivariable models for
211
auxiliary equipment can significantly slow optimization, and so care must be taken to balance model accuracy with
212
computational tractability.
213
3.2
214
To avoid unnecessary computational complexity, we use a linear dynamic model for temperature evolution in the
215
buildings. This model allows us to account for multiple interacting temperature zones and determine the optimal
216
temperature profile for each zone. Without any control action, the temperature trajectory is determined completely
217
by interactions with the external environment and among zones. Using an index i ∈ I for individual zones, a simple
218
model is
us
cr
ip t
204
M
an
Temperature Model
X dfi = −κoi (fi − f o ) − κii0 (fi − fi0 ) + gio , dt 0
d
αi
(15)
i∈I
in which fi is the temperature of zone i, αi is the thermal capacitance of zone i, f o is the external temperature, and
220
gio is the direct heat addition to zone i due, e.g., to solar radiation or heat production within the zone. In general, gio
221
is a function of ambient conditions and zone occupancy, but we assume a forecast is available. The coupling among
222
the zones and ambient is determined by the interaction coefficients κii0 and κoi . To provide control over temperature
223
evolution, a term gi is added to represent direct heat addition to (or rejection from) zone i by the controller.
Ac ce p
te
219
224
To use this model within the MILP formulation, we must convert from the continous-time ODE (15) to a discrete-
225
time difference equation of the form (9). Assuming all quantities except fi are constant between times t and t + ∆,
226
(15) is discretized using standard formulas to give
fi (t + ∆) =
X
(λii0 fi (t) + ωii0 gi0 ) + θi f o (t).
(16)
i0 ∈I 227
with the λ, ω, and θ parameters related to various matrix exponentials involving the zone masses and the interaction
228
coefficients κ. With minor modification, the ambient temperature f o can also be allowed to vary within each timestep,
229
e.g., using a first-order hold instead of a zero-order hold. In either case, the result is a linear difference equation, and so
230
it can be embedded within an MILP optimization problem assuming forecasts for f o and g o are available. Switching
10
Page 10 of 38
Return to Chillers
Hot Section
Return from Building Tank Discharging
Tank Charging
Supply to Building
ip t
Supply from Chillers
cr
Cold Section
us
Figure 2: Diagram of stratified tank model. Tank is treated as two separate volumes of water that interact with each other and the ambient. Only two of the four streams are active at a given time depending on whether the tank is charging or discharging. back to discrete-time index notation, t is treated as an integer, with xt denoting the (constant) value of quantity x
232
between times (t − 1)∆ and t∆ for the predetermined period length ∆. After lumping θi fi (t) into θit , we arrive at the
233
form given in (9). Note that we must add an extra index k to ωii0 so that the sign of gklt can be adjusted appropriately
234
for different resources. For example gklt is nonnegative for both hot and chilled water, but for chilled water, a negative
235
sign must be added to ensure that temperature within the zone is reduced.
236
3.3
237
A typical model for chilled water storage tanks is a single tank split into a hot and a cold section [47], as illustrated
238
in Figure 2. Whenever volume is withdrawn from one section, it is added to the opposite section, and thus total tank
239
volume remains constant. Cooling capacity is stored within the cold section of the tank, and the useful cooling amount
240
is determined by the total enthalpy therein. To arrive at the approximate model (2), a more accurate nonlinear tank
241
model is linearized at a nominal operating point, and at that operating point, we make the following assumptions:
Ac ce p
te
Storage Model
d
M
an
231
242
• The temperature of the cold section is equal to the temperature of the supply from chillers.
243
• The temperature of the hot section is equal to the temperature of the building return stream.
244
• Heat loss from the tank to the ambient is negligible compared to heat exchange between the sections.
245
A detailed derivation can be found in Appendix D.
246
Two potential issues with this model are that it overestimates the decay of enthalpy within the tank and that it does
247
not independently model the “cold volume” within the tank. If the simplified model (21) predicts significantly lower
248
efficiency for the storage tank, then the controller will be overly conservative with respect to utilization of storage, but
249
this effect could be attenuated by adjusting the decay coefficient σ. To model volume, the one-state linear model could
250
be replaced by a two-state linear model without significantly affecting solution times, but such an addition would 11
Page 11 of 38
be unlikely to affect optimal behavior. We note that with closed-loop operation and adequate state estimation, the
252
plant/model mismatch for the storage tanks is not a significant issue. Finally, the nonlinear stratified tank model could
253
be approximated more closely using a piecewise linearization similar to the nonlinear equipment models, but these
254
additions would increase solution times. Furthermore, because the stratified tank model is itself an approximation,
255
state estimation and closed-loop operation would still be recommended.
256
3.4
257
One of the key assumptions of this approach is that the equipment and building models are reasonably accurate.
258
While a detailed procedure for model fitting and maintenance is beyond the scope of this work, we recommend some
259
general procedures. For equipment models, empirical or semi-empirical models should be used so that they can be
260
periodically updated (e.g., every week) to account for the most recent historical data. For example, the unknown
261
coefficients for each chiller model in [41] can be found via linear regression, and the piecewise-linear approximation
262
is then updated by evaluating the new model at a finite number of operating points (see Appendix B). Because this
263
process is computationally cheap, it can be performed as often as necessary. However, while optimization time is
264
insensitive to small changes in the equipment models, more extreme changes (such as altering the convexity of an
265
equipment model) can significantly slow optimization. Thus, care must be taken to avoid major updates to the models,
266
for example by using nonlinear regression and penalizing large deviation from the old model. For the buildings, grey-
267
box models such as RC networks have been shown to be reasonably accurate [48] and also amenable for use with
268
standard system identification techniques. Because these models represent dynamics (and not just static relationships
269
like the equipment models), they require online use of a state estimator such as a Kalman Filter as in [49]. Note that
270
the building model may also need to be periodically updated to account for seasonal changes.
271
3.5
272
To more finely tune optimal behavior, constraints or fictitious costs can be added to MMILP . To reduce equipment
273
switching, a penalty coefficient for u+ it can be added to the objective function. Thus, the optimal solution would tolerate
274
higher economic cost if it can reduce the number of switches. To account for time-varying resource consumption by
275
the building but not by HVAC equipment, a forecast term can be added on the right-hand side of (14) so that the
276
demand charge is calculated using overall building-wide resource consumption. Such an addition could afford more
277
flexibility, as there may already be an intrinsic peak not created by HVAC equipment. Other HVAC-derived resource
278
consumption in zones (e.g., electricity consumption in zone air handlers) is also readily modeled by modifying the
279
gikt sum in the demand constraint (1). For recursive feasibility, a periodic constraint sk0 = skT (i.e., that the initial
280
and finale charges must match) could be added for the TES system. Similar periodicity can be enforced on zone
ip t
251
Ac ce p
te
d
M
an
us
cr
Model Validation
Further Considerations
12
Page 12 of 38
281
temperatures flt to ensure recursive feasibility with respect to zone comfort constraints. Pre-scheduled maintenance
282
or unit breakdowns are treated easily by setting appropriate ujt to zero. By using an abstract representation of the central plant, MMILP also accommodates additional resources and gen-
284
erators without reformulation. For example, with a suitable model for electricity production as a function of gas
285
consumption, a gas turbine could be considered simply by adding the corresponding index j and parameters ζiknt .
286
Absorption chillers can also be treated by adding steam (or other heat source) as a new resource k and augmenting the
287
appropriate models. Thus, this formulation is applicable to a wide variety of systems.
288
3.6
289
One of the key requirements of this formulation is that central-plant equipment can be treated as operating at steady
290
state within each sample time. Although there will certainly be transient startup and shutdown dynamics, we assume
291
that they resolve sufficiently quickly so that the static equipment models describe the actual equipment operating
292
characteristics throughout a large portion of the sample. Because discrete variables are present for equipment on/off
293
states, single-period deviations (for example, a flat increase in electricity consumption during the startup period)
294
could be considered in the cost calculation, but incorporation of detailed nonlinear dynamics would require significant
295
extension of this framework. In addition, for systems with a large number of temperature zones (greater than 50),
296
solving the complete problem in real time becomes more challenging, as the large number of continuous variables
297
needed to model the zones adds significant overhead to the MILP optimization. In such cases, it is necessary to
298
decompose the original problem into subproblems, for example a building temperature subproblem using a single
299
aggregate model of the central plant (e.g., constant COP) followed by a more detailed subproblem of just the central
300
plant. We will not discuss these ideas in this work, but we do note that the various subproblems can readily be
301
formulated as special cases of the full formulation described here, and thus the same considerations are relevant.
302
Finally, significant transport delays in delivering resources to the building are not explicitly considered, as they are
303
assumed to be short with respect to the sample time. Higher-order building models could be used to approximate
304
transportation time, although direct treatment may be necessary for significant delays.
305
4
306
In each of the following examples, price and demand profiles have been provided by Johnson Controls using data
307
from real buildings. Chillers and HRCs were modeled using Gordon-Ng models [41], while the boiler was assumed to
308
operate at constant efficiency. Pumps are modeled by empirical polynomial fits, and cooling towers are modeled with
309
a simplified explicit model [46]. All models are approximated as piecewise-linear, and coefficients for all equipment
310
models are provided by Johnson Controls. All timesteps ∆ are 1 h. Optimization was performed using Gurobi 6 on a
cr
ip t
283
Ac ce p
te
d
M
an
us
Limitations
Simulation Results
13
Page 13 of 38
Zone 1
g4o
Ambient
Zone 3
T1
Zone 2
fo
T3
T2
Zone 4
T4
Chiller 1 g1
1 2 Pumps
Chiller 2
Central Plant
Chiller 3
g3
us
Storage
ip t
fo
cr
g2o
Chiller 4
an
Figure 3: Example of combined system with building and central plant. desktop computer with 8 GB RAM and a 2.66 GHz × 4 CPU.
312
4.1
313
In this example, we model both the building and the central plant, with the equipment and temperature zones as shown
314
in Figure 3. For the building subsystem, we consider four temperature zones, two of which are “air” zones, and two
315
of which are “mass” zones. The air zones 1 and 3 must be maintained within comfort bounds, and they couple only
316
with the direct cooling input and the corresponding mass zones (2 and 4, respectively). The mass zones interact with
317
the ambient and are coupled to each other, but there are no comfort bounds for these zones. In effect, we have a
318
second-order linear model for temperature evolution within the two air zones, but in principle arbitrary order linear
319
models for each zone could be considered via appropriate choice of parameters. For the central plant subsystem, we
320
consider four conventional chillers and two pumps. There is also a chilled water storage tank available for use. Thus,
321
chilled water is produced in chillers and pumped either to the air zones to provide immediate cooling or to the storage
322
tank for later use.
M
311
Ac ce p
te
d
Building Cooling
323
Within the optimization, we use a time-varying forecasts for secondary demand, electricity price, and ambient
324
temperature (all of which are parameters). These profiles are shown in Figure 4, from which the daily cycle is clear.
325
In particular, we note that electricity price is highest when ambient temperature and cooling demand are highest, and
326
so there is significant potential for cost savings via load shifting. Solving MMILP with a 48 h horizon leads to a model
327
with 2.3k constraints and 4.4k variables, of which 1.4k are integer (after presolve to reduce the problem size). After
328
30 s, the optimality gap was 0.81%, and after 60 s, the gap was 0.17%. Due to nonconvexity in the pump models, it is
329
difficult to completely close the optimality gap, and thus the solver was stopped after 60 s. However, since equipment 14
Page 14 of 38
5.0 2.5 0.0
0
12
24
36
0
12
24
36
0
12
48
ip t
100 80 60
28 26 24 22 20 18 16
24
48
us
20
cr
40
36
an
Secondary Demand (MW)
7.5
Ambient Temp. (◦ C)
Electricity Price ($/MW)
10.0
48
Time (h)
M
Figure 4: Electricity price, auxiliary cooling demand, and ambient temperature for example system. and building models are only approximate, there is no need to solve to optimality in a practical implementation due
331
to feedback. The cooling and temperature dynamics for this solution are shown in Figure 5. We see that the optimal
332
profile is to charge the storage tank overnight (as shown by the large negative values for the “Storage” curve). Note
333
that “Demand” refers to total cooling demand, i.e., the sum of secondary demand from Figure 4 and primary demand
334
delivered to the temperature zones as determined by the optimization. For the temperature zones, there is no pre-
335
cooling, as zone temperature is left to float until it reaches the upper comfort bound of 22 ◦ C (recall that Zones 2 and
336
4 do not have comfort bounds, as they represent the surrounding solid mass for the air zones 1 and 4). Equipment
337
utilization shown as a Gantt chart in Figure 6. We see that equipment is switched off during peak hours to reduce
338
cost. However, due to the demand charge, it is more expensive to run all of the equipment at full capacity. Indeed, the
339
optimal peak loading would be very difficult to determine using heuristics alone, which illustrates the superiority of
340
optimization-based approaches.
341
4.1.1
342
Large commercial buildings have many temperature zones, and furthermore, it may become important to consider
343
longer horizons. Thus, we examine the scaling of our formulation with respect to building size by solving models of
344
various sizes with a 7 day horizon. Sizes and solution times for each optimization problem are shown in Table 1.
345
Note that changing the number of zones does not affect the number of discrete variables. Despite the slowdown as the
346
number of zones increases, performance is sufficiently fast for online use. We see that in all instances, a solution no
Ac ce p
te
d
330
Scaling with Number of Zones
15
Page 15 of 38
30 27 24 21 18 15
Demand Production Storage Unmet 0
12
24
36
48 Ambient Zone 1 Zone 2 Zone 3 Zone 4
0
12
24
36
ip t
Chilled Water (MW)
Temperature (◦ C)
60 40 20 0 −20 −40
48
cr
Time (h)
Chiller 1 Chiller 2 Chiller 3 Chiller 4 Pump 1 Pump 2 0
6
12
18
an
us
Figure 5: System and zone temperature dynamics for building cooling optimization.
24
30
36
42
48
M
Time (h)
Figure 6: Equipment utilization for building cooling optimization. Boxes indicate when a unit is on, while lines show equipment loading between minimum and maximum capacity. more than 1% of the best possible solution was obtained within just over 6 minutes. Recalling from that all models are
348
only approximate, a 1% optimality gap is sufficient for online use, in which optimization must be completed roughly
349
every 15 minutes.
te
d
347
Often, airside variables often do not appreciably change after the root node if sufficient active storage is available,
351
which allows a near-optimal forecast of demand to be obtained from even highly suboptimal solutions to the overall
352
problem. If such a forecast of primary demand is available, then the optimization problem need not include the building
353
model at all, and all demand terms are lumped into the secondary demand parameter φkt . To illustrate the benefits of
354
omitting the airside model, we consider the 25 zone example from before, sum the gklt variables, add them to φkt , and
355
repeat the optimization without the building model (9) through (11). As shown in the “No Building” entry in Table 1,
Ac ce p
350
Table 1: Problem sizes for 7 day horizon and varying number of zones. Sizes are as reported by Gurobi following pre-solve reductions. “1% Gap Time” is the wall-clock time needed to obtain a solution that is no more than 1% worse than the optimal solution. The “No Building” case simply meets a forecast of the optimal 25-zone demand. Size
Variables (k)
Integers (k)
Constraints (k)
1% Gap Time (s)
2 Zones 8 Zones 15 Zones 25 Zones No Building
6.2 12.3 18.7 29.2 4.4
1.8 1.8 1.8 1.8 1.8
3.8 7.8 11.9 19.0 2.5
4.0 26.0 106.3 365.5 2.0
16
Page 16 of 38
Table 2: Problem sizes for 7 day horizon and varying central plant size. Sizes are as reported by Gurobi following pre-solve reductions. “1% Gap Time” is the wall-clock time needed to obtain a solution that is no more than 1% worse than the optimal solution. Variables (k)
Integers (k)
Constraints (k)
1% Gap Time (s)
Nominal 2× 4× 8×
4.4 4.4 4.4 4.4
1.8 1.8 1.8 1.8
2.5 2.5 2.5 2.5
2.0 2.5 3.0 1.5
ip t
Size
the solution time for this instance is only 2 s. This result indicates that decomposition is be highly advantageous when
357
the number of zones is large.
358
4.1.2
359
Finally, we wish to address how the amount of equipment affects problem scaling. For these tests, we use the waterside-
360
only problem from the previous section and increase the amount of equipment present. Note that demand and storage
361
capacity are also scaled accordingly. Instance sizes and solution times are shown in Table 2. Here we see that scaling
362
with number of identical units is almost insignificant because on/off states are modeled using general integer variables
363
and not separate binary variables for each piece of equipment. Thus, adding additional units requires changing only
364
the upper bounds µj on ujt , which does not affect the number of variables or constraints, as seen in Table 2. Note
365
that if additional distinct types of equipment (e.g., chillers with different capacities) were added, then scaling would
366
be less favorable.
367
4.2
368
In central plants serving campuses or large buildings, simultaneous production of chilled and hot water is often re-
369
quired, e.g., due to secondary demand. Because of demand charges and heat-recovery chillers, these two process are
370
codependent and thus require simultaneous treatment. Below, sample optimization results are presented for the system
371
in Figure 7, which contains chilled and hot water loops (interacting via HRCs), both with active storage tanks and
372
pumps. The chillers also require cooling water produced in cooling towers, which are not shown in the diagram.
us
cr
356
te
d
M
an
Scaling with Number of Generators
Ac ce p
Combined Heating/Cooling Systems
373
In the optimization, a prediction horizon of 48 h is used. All generators are given on/off dwell times of 2 h, which
374
requires that after any switching event, that unit may not be switched again for at least 2 h. Forecasts for electricity
375
and gas price are given in Figure 8 along with heating and cooling demand forecasts. Here, the main periodicity
376
is in cooling demand and electricity prices, with a smaller effect in heating demand. For the system in Figure 7, a
377
solution with an optimality gap of 1.8% was found after 10 s of CPU time. Optimization was stopped after 300 s with
378
an optimality gap of 0.6%. Chilled and hot water production are shown in Figure 9, and equipment usage is given
379
in Figure 10. In this case, we see that production is decreased when electricity is expensive, although the effect is 17
Page 17 of 38
ip t cr us
12
an
Hot Water
M
0
24
Electricity
te
100 75 50 25 0
Chilled Water
36
48
36
48
Gas
d
120 90 60 30 0
0
12
24 Time (h)
Ac ce p
Price ($/MW)
Demand (MW)
Figure 7: Diagram of central plant for combined heating/cooling system. Storage tanks are available for both hot water (red) and chilled water (blue). Cooling towers (not shown) are connected to the chillers via a cooling water loop.
Figure 8: Resource demand and price forecasts. The central plant must supply time-varying demand hot and chilled water demand subject to time-varying electricity price and constant gas price.
380
nowhere near as pronounced as in the previous section. This difference is due to the demand charge and the reduced
381
flexibility of having to meet a fixed forecast of resource demand. A breakdown for electricity costs is provided in
382
Table 3. From this table, we see that chillers and HRCs are by far the largest contribution to total cost, indicating that
383
auxiliary units need not be modeled to a high degree of accuracy.
384
4.2.1
385
Because discrete decision variables are the main cause of computational difficulty, we compare MMILP to an alter-
386
native formulation that requires only continuous decision variables. By ignoring minimum capacities and switching
387
constraints, the piecewise-linear approximations of nonlinear generator models are replaced with single linear approx-
388
imations. These new models are quite inaccurate, but they no longer require the discrete variables vjmt associated
Comparison to Simpler Models
18
Page 18 of 38
Demand Production Storage Unmet
0
12
24
36
48 Demand Production Storage Unmet
20 10 0 −10
0
12
24
36
48
ip t
30
cr
Chilled Water (MW) Hot Water (MW)
60 45 30 15 0 −15 −30
Time (h)
us
Figure 9: Production of chilled and hot water for MMILP . Negative values for “Storage” indicate that storage tanks are being charged, while positive values denote discharge.
M
an
Chiller 1 Chiller 2 Chiller 3 Chiller 4 HRC 1 HRC 2 Boiler 1 Cold Pump 1 Cold Pump 2 Hot Pump 1 Hot Pump 2 Tower 1 Tower 2 6
12
18
d
0
24
30
36
42
48
Time (h)
Ac ce p
te
Figure 10: Equipment utilization for MMILP . Colored bars indicate that a unit is on, while solid lines show loading fraction between minimum and maximum capacity. Note that a first-in first-out queueing system was applied to disaggregate integer variables ujt into individual on/off states.
389
with piecewise-linearization, and indeed the on/off variables ujt can be ignored, as a unit’s on/off state is determined
390
by whether any qjt 6= 0. We refer to this scheme as MLP . Using the same parameters given in Figure 8, MLP is solved
391
as a linear program, and the equipment utilization plot in Figure 11 is obtained. In particular there is much more
392
equipment switching for MMILP compared to MMILP , as MLP cannot directly constrain equipment on/off switching.
393
Note that there are other solutions with the same total cost due to non-uniqueness of this formulation.
394
Table 4 compares costs for MMILP and MLP . Costs as predicted by the (piecewise-) linear models are shown in
395
the “Model” columns To provide a fair comparison, all electricity costs are recalculated using the underlying nonlinear
396
equipment models and shown in the “Actual” columns. In this example, we see that the fully linear models of MLP
397
are reasonably accurate due to units operating near capacity (which is where the linear approximations were made).
398
However, even in this favorable case, MMILP is noticeably more accurate (as demonstrated by the small difference
399
between “Model” and “Actual” values of electricity consumption). Furthermore, when the optimal policy requires
400
variable operation, there may be significant inaccuracy imposed by the linear models of MLP relative to the more
19
Page 19 of 38
Table 3: Electricity charges by equipment type for MMILP . Note that “Cost” column accounts for time-varying use charges, and so percentages do not exactly match with “Usage.” Usage, MWh
Cost, 1000$
168.8 106.7 24.8 6.0 1.6
(54.83%) (34.64%) (8.06%) (1.95%) (0.52%)
11.14 6.93 1.62 0.39 0.11
(55.16%) (34.33%) (8.05%) (1.94%) (0.53%)
Total
307.9
(100.00%)
20.19
(100.00%)
cr
Chillers HRCs Cold Pumps Hot Pumps Cooling Towers
ip t
Equipment
0
6
12
18
an
us
Chiller 1 Chiller 2 Chiller 3 Chiller 4 HRC 1 HRC 2 Boiler 1 Cold Pump 1 Cold Pump 2 Hot Pump 1 Hot Pump 2 Tower 1 Tower 2
24
30
36
42
48
M
Time (h)
401
accurate piecewise-linear models of MMILP .
d
Figure 11: Equipment utilization for MLP . Equipment models were approximated using a single linear term, and integer decision variables were not used.
In addition to greater model accuracy, MMILP also out-performs MLP in terms of unit switching. As seen in
403
Figure 11, there times when a unit is switched on for one period and then immediately switched off. Such behavior
404
cannot occur with MMILP because 2 h dwell times are enforced directly. Implicit in the problem formulation is that
405
equipment dynamics are irrelevant for the given timesteps, but without being able to enforce dwell times, solution
406
quality of MLP may deteriorate significantly when applied to the true system due to the burden (both in terms of
407
resource consumption and equipment wear) of extra switching. Although various ad-hoc additions could be made
408
to deter rapid switching, the discrete on/off nature of equipment scheduling decisions cannot be fully captured by
409
continuous variables. Thus, while MLP is much easier to solve due to its lack of discrete decision variables, it
410
performs noticeably worse than MMILP due to its significant limitations.
411
4.2.2
412
Currently, HVAC central plants with storage tanks are operated using heuristics instead of any sort of online opti-
413
mization [4, 9, 16, 47]. As will be demonstrated, heuristic methods can lead to significant cost increases if heuristics
414
fail to fully account for pricing structure. One such heuristic is to run equipment at full capacity over night to fully
415
charge storage tanks, and then keep equipment off during peak hours. To simulate this heuristic, we embed these
Ac ce p
te
402
Comparison to Heuristic Methods
20
Page 20 of 38
Table 4: Electricity cost comparison for MMILP , MLP , and heuristic approaches. “Change” row is total cost relative to Formulation MMILP Optimal MLP Heuristic MMILP Charges Optimal MMILP (1000$) Model Actual Model Actual Model Actual 20.21 2.24 22.45
20.19 2.24 22.43
20.32 2.21 22.53
21.29 2.63 23.93
21.05 3.59 24.63
21.02 3.59 24.61
Change
−
−
+0.37%
+6.68%
+9.73%
+9.72%
0
6
12
18
an
us
cr
Chiller 1 Chiller 2 Chiller 3 Chiller 4 HRC 1 HRC 2 Boiler 1 Cold Pump 1 Cold Pump 2 Hot Pump 1 Hot Pump 2 Tower 1 Tower 2
ip t
Usage Demand Total
24
30
36
42
48
M
Time (h)
Figure 12: Equipment utilization with heuristic rules enforced in the model. All chillers must be fully on from 1am to 7am and fully off from 2pm to 6pm. HRCs are left free. considerations as hard constraints in the optimization problem. Electricity cost is shown in Table 4, and equipment
417
usage is shown in Figure 12. In this case, the boiler must be activated to provide sufficient hot water, as production
418
is lost due to the heat recovery chillers being shut off in the afternoon. From these results, we see that the heuristic
419
solution does not perform well, even though a majority of its decisions were made optimally. Had the entire schedule
420
been determined by similar rules, performance would deteriorate even further. The main problem with such heuristics
421
is that the presence of peak demand charges requires choosing an upper bound for total electricity that accounts for
422
the tradeoff between demand and use charges. Although more sophisticated heuristics could better account for peak
423
demand charges, it is unlikely that performance would approach optimization-based methods.
424
4.3
425
The load Q supplied by a chiller can be varied by adjusting both the volumetric flow V through the chiller and the
426
supply temperature TCHWS . In the previous examples, a degree of freedom has been eliminated by assuming a fixed
427
TCHWS . While fixed supply temperature is a common practice in existing buildings, there is potential for cost savings
428
by allowing both load and supply temperature to vary independently. The main reason to allow TCHWS to vary is
429
that there is a tradeoff between chiller electricity usage and pump electricity usage that depends on TCHWS . Thus, to
430
deliver a specified load Q with the lowest electricity consumption potentially requires a different supply temperature
Ac ce p
te
d
416
Variable Stream Temperature
21
Page 21 of 38
431
for each possible on/off configuration of pumps and chillers. Based on common forms of chiller models [41], the most
432
relevant set of variables is (Q, TCHWS ), although if the temperature of water entering the chiller is known, then any
433
two of {Q, TCHWS , V } can be chosen. The optimal supply temperature is given by some relationship TCHWS = f (Q, u), in which Q is a vector with the
435
loads of all chillers, and u is a binary vector with the on/off states of all chillers and pumps. Even if we could find this
436
function f via offline optimization, there are too many degrees of freedom to embed a piecewise linear approximation
437
while maintaining fast solution times. Thus, a better alternative is to pick the optimal TCHWS online, which can be
438
considered by linearizing the chiller models in two variables instead of one. However, this addition can significantly
439
increase the computational burden, and so it is important to know what cost savings, if any, can be obtained.
440
4.3.1
441
We first wish to determine whether the piecewise-linear approximations made by MMILP are sufficiently accurate to
442
determine the optimal supply temperature. To assess this accuracy, we compare MMILP to a more detailed alternate
443
formulation MMINLP that uses the full nonlinear equipment models and includes energy balances to accurately cal-
444
culate the temperatures of the cooling water streams (which are assumed constant in MMILP ). Full details of MMILP
445
can be found in Appendix C. Due to its computational complexity, MMILP is intractable when multiple time periods
446
are considered, and so we perform a series of 1-period optimizations to calculate the optimal equipment configuration
447
and supply temperature for varying total system loads. Optimal supply temperature using the piecewise-linear and
448
nonlinear models are shown in Figure 13. From this plot we see that the accuracy of the piecewise-linear approx-
449
imation is quite good with the only noticeable differences are near switch points where MMINLP and MMILP select
450
different numbers of pumps. Indeed, the maximum temperature difference between the two models is below 3 ◦ C, and
451
at loads above 75% capacity where low cost is particularly important, the temperature difference is well below 1 ◦ C.
452
Thus, even with the simplifications made by MMILP , it is still possible to accurately determine the optimal supply
453
temperature.
454
4.3.2
455
To assess the effects of supply temperature on total electricity cost, we also calculate the optimal electricity consump-
456
tion (as predicted by MMILP and MMINLP ) as a function of varying total load. To provide a fair comparison, we obtain
457
the loads for each piece of equipment using optimization and calculate electricity consumption using the full nonlinear
458
models. These values are shown in Figure 14. We see once again that there is good agreement between MMILP and
459
MMINLP , indicating that the piecewise-linear approximations are quite accurate. We also observe that there is only
460
small savings to be gained by allowing variable water temperature. While the maximum individual savings is 4.1%
us
cr
ip t
434
Ac ce p
te
d
M
an
Accuracy of MMILP
Electricity Cost Savings
22
Page 22 of 38
11
Return MMILP Supply MMINLP Supply
10 9 8 7 6 5 0
10
20
30
40
50
cr
Total Load (MW)
ip t
Water Temperature (◦ C)
12
us
Figure 13: Optimal supply temperature using the MMILP and MMINLP models. The MMILP model uses piecewiselinear models, while the MMINLP model uses nonlinear models and includes cooling water energy balances. Both models have variable supply temperature. at a relatively low total load of 4 MW, when the plant is operating above half capacity, savings are significantly re-
462
duced. Average savings (weighted by total load) are 1.2%. Because of the increased computational burden, online
463
optimization of supply temperature may not worthwhile. However, a strength of MMILP is that this tradeoff can be
464
evaluated.
465
4.4
466
We conclude with an example of the system being run in closed loop. For simplicity of discussion, we consider a
467
system with only chillers. As described in Section 3.3, the model used for the storage tanks is only approximate,
468
and so we provide state feedback by giving the optimizer an accurate estimate of tank enthalpy (i.e., initial condition
469
sk0 ) at each time point. Figure 15 shows electricity price and equipment utilization over a 1 week period. In this
470
case, the electricity price and cooling demand were chosen to be periodic, and thus the operating profile is periodic as
471
well (excepting for small start and end effects). Note that each optimization has a a 24 prediction horizon with a 1 h
472
timestep.
M
an
461
Ac ce p
te
d
Closed-loop Operation
473
For comparison, we also simulate the system in open loop, i.e., without providing the actual charge of the storage
474
tank. Enthalpy, temperature, and volume profiles for the storage tank are shown in Figure 16. In this figure, only the
475
first 72 h of operation are shown, after which open-loop operation leads to tank overflow. From this simulation, we
476
see that open-loop operation leads to instability due to model mismatch. Because this model overestimates decay of
477
energy, open-loop operation leads to over-charging the tank, and eventually leads to a tank completely full of cold
478
water. By contrast, when the system is operated in closed-loop, the optimizer is aware of the true current charge of the
479
tank, and thus does not over-fill. Furthermore, the one-step prediction for tank enthalpy using the simplified model is
480
sufficiently accurate, as is the estimated temperature rise when the tank is non-empty. This agreement indicates that
481
the simplified model is accurate enough for closed-loop operation. 23
Page 23 of 38
6 4 2 0
0
10
0
10
20
30
40
20
30
50
cr
1.02
Relative Load (-)
ip t
Variable MMILP Constant MMILP Variable MMINLP Constant MMINLP
8
us
Electric Load (MW)
10
1.00 0.98
40
50
an
0.96
Thermal Load (MW)
d te
0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 0.0
Ac ce p
Electricity Price ($/MWh)
M
Figure 14: Comparison of optimal electricity consumption. Legend refers to whether chilled water supply temperature is variable or constant, and whether the approximate MMILP or full nonlinear MMINLP was used. Bottom plot is normalized to Constant MMILP .
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
Chilled Water (MW)
45 30
Demand Production Storage Unmet
15
0
−15
−30 0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
7.0
Chiller 1 Chiller 2 Chiller 3 Chiller 4
0.0
Time (days)
Figure 15: Electricity price, production, and equipment utilization for closed-loop simulation.
24
Page 24 of 38
0
CL Model
CL Actual
12
24
36
48
60
72
0
12
24
36
48
60
72
0
12
24
36
48
cr
Volume m3
8 6 4 2
Time (h)
us
Temp. Rise (◦ C)
10
0
ip t
8 7 6 5 4 3 2 1 0
OL Actual
Enthalpy (MJ)
OL Model 0 −50 −100 −150 −200 −250 −300 −350
60
72
an
Figure 16: Dynamics of storage tank for closed-loop and open-loop operation. More negative values of enthalpy indicate more stored cooling capacity. Temperature rise is calculated relative to the starting tank temperatures.
5
Conclusions
483
In this work, we have proposed a framework for the formulation of mathematical programming models for the op-
484
timization of equipment usage in a commercial HVAC systems. This framework accounts for time-varying utility
485
prices, demand charges, occupant comfort constraints, and on/off equipment dwell times. Using a series of reformu-
486
lations, including piecewise linearizations of unit models and an aggregate formulation for equipment switching, we
487
have developed a mixed-integer linear programming model, MMILP , that can be solved quickly and leads to solutions
488
with lower cost and less equipment switching than a similar scheme without discrete variables. In its most general
489
form, MMILP includes a linear model of building dynamics that is used to estimate cooling or heating loads required
490
to induce a given temperature trajectory in the building zones. Thus, resource demand can be determined within the
491
optimization, and use of passive TES is implicit. Alternatively, if a forecast of resource demand is available, a building
492
model need not be included at all, and MMILP can determine optimal usage of central plant equipment and active TES.
493
This formulation requires static models for resource consumption as a function of load for central plant equipment, as
494
well as dynamic models for storage tank charge and (possibly) building temperature evolution, but these models can
495
be obtained from manufacturers’ or historical data and can be updated as additional operating data is obtained.
Ac ce p
te
d
M
482
496
With respect to computational tractability, we have demonstrated that the formulation scales favorably with number
497
of temperature zone and number of units in the central plant. For a 7 day horizon with 1 h timestep, quality suboptimal
498
solutions can be found for a 25 zone building within 400 s, and with a demand forecast instead of the building model,
499
such solutions are found after only 10 s. For a central plant with both heating and cooling equipment, we show 6.7%
25
Page 25 of 38
improvement over a simpler formulation without discrete variables and 9.7% improvement over a simulated heuristic.
501
Thus, MILP-based optimization is the best choice for equipment scheduling in central plants. In addition, we have
502
shown that while variable supply temperature can be considered by the formulation, it does not significantly decrease
503
total cost, with a maximum reduction of 4.1% at low total loads and an average reduction of 1.2% across the entire
504
operating range. Therefore, a fixed supply temperature assumption is suitable for online use, but variable temperature
505
can be considered as well with additional computational effort. Finally, we have shown that while the simplified first-
506
order storage tank model is not accurate for long-term open-loop predictions, performance does not deteriorate in a
507
closed-loop implementation with measurement feedback. Based on its flexibility and modest computational time, our
508
model is suitable for use in new or existing buildings, and it can achieve significant cost reduction relative to heuristic
509
scheduling methods.
510
6
511
Funding, equipment models, and sample data provided by Johnson Controls, Inc.
512
Appendix A
513
A.1
514
The following sets are used in the abstract model description.
M
an
Acknowledgments
Ac ce p
temperature zones
te
d
Notation
Sets
i∈I 515
516
us
cr
ip t
500
j∈J
type of generators
k∈K
resources/utilities
t∈T
time periods
The following sets are used for piecewise linearization.
517
m∈M
interpolation regions for equipment models
n∈N
interpolation points for equipment models
Mj ⊆ M
interpolation regions for equipment j
Njm ⊆ N
interpolation points contained in region m for equipment j
518
A.2
Parameters
519
The following parameters are given based on available forecasts or desired system performance.
26
Page 26 of 38
523
demand charge
βkt
penalty for unmet secondary demand
Ξjt ( · )
nonlinear model function
L U ψit , ψit
lower and upper temperature comfort bounds for zones
Πmax k
previous maximum resource demand
ip t
ρmax k
cr
forecast resource prices
Number of units of each type available for use
ζjknt
interpolation points for piecewise-linear models
σk
fractional retention of stored resource
δj+ , δj−
minimum on/off dwell time for generators
Πk
Upper bound for resource purchase
Υk
Bound for single-period storage charge/discharge
Σk
Storage capacity
an
µj
us
The following parameters are calculated from equipment models.
M
522
ρkt
The following parameters are calculated from building models. αi
thermal capacitance of zone
κoi , κii0
heat transfer coefficients with ambient and between zones
λii0 , ωii0 k
coefficients in discretized model for building temperature evolution
θit
time-varying temperature disturbance for zones
d
521
forecast secondary resource demand
te
520
φkt
Ac ce p
524
− χ+ it , χit
penalties for comfort violation in zones
Γik
upper bound on resource flow to zones
525
Decision Variables
526
The following decision variables are primarily for the waterside subsystem. ujt ∈ {0, . . . , µj }
integer number of generators on during a given time
− u+ jt , ujt
number of generators switched on/off at the beginning of a time period
qjkt
production of resources in equipment
pkt ∈ [0, Πk ]
amount of purchased resources
pmax ∈ [Πmax , ∞) k k
maximum single-period purchase of resource k
ykt ∈ [−Υk , Υk ]
charge (< 0) or discharge (> 0) of stored resources
skt ∈ [0, Σk ]
stored inventory at the end of a time period
bkt ∈ [0, φkt ]
unmet secondary demand
527
27
Page 27 of 38
528
The following variables are added for the piecewise-linear approximations of equipment models. vjmt ∈ {0, . . . , µj }
number of units operating in each interpolation region
zjmnt ∈ [0, µj ]
weighting of interpolation points in each region
529
zone temperature of zone at end of a time period
gikt ∈ [0, Γik ]
primary demand of resources by zones
fit+ , fit−
positive and negative comfort violation for zones
532
Appendix B
533
The general relationship we are trying to capture is
Piecewise Linearization
x ∈ X ⊂ Rn ,
with X and Y compact. We assume that X can be partitioned into a finite number of polyhedrons {P1 , . . . , PM }, with SM X = i=1 Pi and int(Pi ) ∩ Pj = ∅ ∀i, j 6= i (that is, the polyhedrons only intersect on their boundaries).
M
535
y ∈ Y = f (X) ⊂ Rm
an
y = f (x),
534
ip t
fit
cr
531
The following decision variables are used when the building airside model is included.
us
530
Since we are modeling an equality relationship, we can choose either a relaxation of f or an approximation of f .
537
To relax f , we create a set-valued mapping F from X to Y with the requirement that f (x) ∈ F (x) ⊂ Y and use
538
the constraint y ∈ F (x) The advantage here is that any solution to the exact nonlinear model is a valid solution to
539
the relaxed model. However, there is now an additional degree of freedom that is not present in the actual nonlinear
540
formulation, and thus the optimal solution using the relaxed model is often infeasible to the true model. Alternatively,
541
we may simply approximate f by creating a piecewise-linear function f˜ such that |f˜(x) − f (x)| ≤ ε ∀x ∈ X with ε
542
some predetermined tolerance. The advantage here is that f˜ is actually a function and does not afford the system any
543
additional flexibility. However, if nonlinear constraints are tight, the approximate model may be infeasible even when
544
the true model was feasible. In our case, the variables corresponding to y are generally associated with production
545
amounts and do not have any tight feasibility requirements; thus, we use approximations of f .
Ac ce p
te
d
536
546
For a single dependent variable (i.e., y ∈ I ⊂ R), algorithms exist to find efficient piecewise linearizations for
547
convex and nonconvex functions [50, 51], and for two variables, triangulations can be used [52, 53]. We could make
548
use of these techniques by splitting f into component functions f1 , . . . , fm and use a different partitioning of X for
549
each, but this greatly increases the size of the resulting MILP. Thus, we require that a single partitioning be used for
550
all components of f , and so many existing techniques cannot be applied. Instead, we simply use a regular partitioning
551
of X into simplexes and use a secant approximation on each. This strategy ensures continuity of the resulting ap-
552
proximation f˜, and it also means f˜(x) = f (x) at extreme points of the Pj . Although such an approximation would
553
likely not minimize error throughout the entire domain, this feature is desirable because the optimal solutions to the 28
Page 28 of 38
554
MILP will often lie on these extreme points. To choose the final partitioning, we simply increase the resolution until
555
the maximum approximation error is below some prespecified value. Though more sophisticated techniques could be
556
developed to create more accurate approximations with fewer subdomains, the proposed scheme is sufficient for the
557
cases we consider.
559
The nonlinear relationships we must model are between the qjkt variables for a given j and t. These equations
ip t
558
often take the form
(17)
cr
(qjk1 t , . . . , qjk2 ,t ) = Ξjt (qjk3 t , . . . , qjk4 t ),
in which we explicitly identify which resources are the independent variables x and which are the dependent variables
561
y. Without loss of generality, we may reparameterize our functions so that all of the qjkt are outputs, i.e., we have
562
only relationships of the form
us
560
(18)
an
ˆ jt (z), (qj1t , qj2t , . . . , qjKt ) = Ξ
for a new dummy parameter z. For example, for f defined on q ∈ [Qmin , Qmax ], we can take z ∈ [0, 1] and add the
564
ˆ We then define the set N to index the mapping q = Qmin + z(Qmax − Qmin ) to Ξ, yielding the reparameterized Ξ.
565
extreme points zn in z-space and set
M
563
(ζjk1 nt , . . . , ζjkK nt ) = Ξjt (zn )
(19)
and leave z as the free optimization variable with q calculated according to (3). By using this general formulation, it
567
does not matter what the specific dependent and independent variables are for each generator.
568
Appendix C
569
To validate the piecewise-linearization approach, we can compare the optimal solution returned from the model MMILP
570
to results obtained from using full nonlinear equipment models and energy balances, which is denoted MMINLP . This
571
system is diagrammed in Figure 17. For simplicity, only one of each piece of equipment is shown. Variables and
572
parameters are described in the sections below. Functional forms of equipment models are shown in Section C.5.
573
C.1
574
Each chiller j has the following variables:
Ac ce p
te
d
566
Nonlinear Equipment Selection Formulation
Chillers
575
• Binary uCH j for its on/off state
576
CH • Continuous QCH for its cooling and electrical loads j and Wj
577
• Continuous VjCHW and VjCHW for its volumetric flows of chilled and cooling water 29
Page 29 of 38
Cooling Tower mA l WlCT
Tamb φamb
QCT l
T CWR
ip t
mCW l T CWS
QCH j
T CHWR
Pump
VjCHW Chiller
WkP
TjCHWS
us
WjCH
cr
VjCW
• Continuous TjCHWS for the supply temperature of chilled water
M
578
an
Figure 17: Diagram of chiller/cooling tower/pump system.
Parameters are c and ρ, the heat capacity and density of water, and T CHWR for the chilled water return temperature.
580
Equations are as follows:
d
579
te
CH CHWS WjCH = uCH , T CWS ) j fchiller (Qj , Tj CH CH CH CH uCH j Qmin ≤ Qj ≤ uj Qmax
Ac ce p
CW CW CW uCH ≤ uCH j Vmin ≤ Vj j Vmax
CHW QCH (T CHWR − TjCHWS ) j = ρcVj
581
The first equation calculates electric load using the nonlinear model fchiller (see (20) in Section C.5). Note that T CWS
582
is the supply temperature of cooling water determined by cooling tower operation. The next two equations enforce
583
minimum and maximum capacities. The final equation calculates the volumetric flow of chilled water using an energy
584
balance.
585
C.2
586
Each pump k has the following variables:
Pumps
587
• Binary uPk for its on/off state
588
• Continuous VlP and VlP for its volumetric flow of chilled water
30
Page 30 of 38
589
Equations are as follows: WkP ≥ uPk fpump (VkP )
ip t
P P uPk Vmin ≤ VkP ≤ uCH j Vmax
The first equation gives electricity consumption as a function of volumetric flow, and the second equation enforces
591
minimum and maximum capacities.
592
C.3
593
Each cooling tower l has the following variables:
cr
590
us
Cooling Towers
• Binary uCH j for its on/off state
595
CT • Continuous QCT for its cooling and electrical loads l and Wl
596
• Continuous mCW and mA l for its mass flows of cooling water and air l
M
597
an
594
Equations are as follows:
d
CW A CWR QCT ) l = ftower (ml , ml , T
te
WlCT = ffan (mA l )
Ac ce p
CT CT CT CT uCT l Qmin ≤ Ql ≤ ul Qmax
598
The first equation calculates cooling duty as a function of mass flows and cooling water return temperature. Note that
599
because of the model form, we need not calculate individual values for T CWS . The second equation uses the tower fan
600
model to give electric load, while the third equation enforces minimum and maximum capacities.
31
Page 31 of 38
601
C.4
Balances
602
Finally, various material balances must be enforced to ensure that sufficient equipment is switched on. These are as
603
follows: Total QCH j ≥Q
ip t
X j
X
VjCHW =
X
j
VjCW =
X
j
X
mCW l
l
CH QCH = j + Wj
j
X
CWR mCW − T CWS ) l c(T
l
W Total =
X
cr
k
X
WjCH +
X
X
WlCT
an
j
WkP +
us
ρ
VkP
k
l
The first balance ensures that enough chilled water is produced to meed the prespecified cooling duty QTotal . The
605
second balance ensures that enough pumping is available for chilled water. The third balance manages cooling water
606
flow among the chillers and cooling towers. The fourth balance is used to calculate the aggregate cooling water supply
607
temperature T CWS , which affects the electricity consumption in the chillers. Finally, the fifth balance calculates the
608
total electricity consumption W Total , which serves as the objective to be minimized.
609
C.5
610
For modeling of chillers, we use the semi-empirical Gordon-Ng model, [41] defined below:
Ac ce p
Equipment Models
te
d
M
604
fchiller (Q, T CHWS , T CWS ) :=
T CWS T CHWS −Q Q + a1 T CHWS + a2 1 − CWS T T CHWS − a3 Q
(20)
611
The parameters a1 , a2 , and a3 are obtained via regression with measured data, and they are related to entropy changes
612
within the process.
613
For pumps, the following black-box empirical model is used:
fpump (V ) := b1 ln (1 + b2 V ) + b3 V + b4
614
615
Here, b1 through b4 are regression coefficients. Finally, for cooling towers, a simplified effectiveness model [46] is used for calculating cooling duty, while for fan
32
Page 32 of 38
v+ , Hhot /Vhot
Hhot Vhot
Tank Discharging
v+ , h+
v− , Hcold /Vcold
cr
Hcold Vcold
ip t
Tank Charging
v− , h−
electricity consumption, a simple cubic fit is used [20]. ftower (mCW , mA , T CWR ) :=
an
616
us
Figure 18: Diagram of stratified tank model. Tank is treated as two separate volumes of water that interact with each other and the ambient. Only two of the four streams are active at a given time depending on whether the tank is charging or discharging.
c1 (mCW )c3 c3 (T CWR − Tamb ) mCW 1 + c2 m A
M
ffan (mA ) := κ(mA )3 c1 , c2 , c3 , and κ are obtained via regression.
618
Appendix D
619
The storage model (2) is an approximation of a stratified tank model [47]. Within this section, we use v+ to denote
620
the volumetric flow for when the tank is charging, and v− for when the tank is being discharged. We assume the tank
621
volume is held constant, and thus model it as consisting of two stratified regions, each with separate volume V and
622
enthalpy H. This is illustrated in Figure 18.
te
d
617
Ac ce p
Tank Model Derivation
623
Within the tank, here is energy exchange between the two layers of the tank, and there is also external exchange
624
between each layer and the ambient. Often the tank is assumed to be well-insulated, which means the only enthalpic
625
change to the tank is due to heat exchange between the two layers (i.e., k = 0 in the model). However, we retain the
33
Page 33 of 38
ambient term for now and arrive at the following mass and energy balances.
627
= −v+ + v− = v+ − v− Hhot Hhot Hcold Hhot =− v+ + h − v− − K − −k − hamb Vhot Vhot Vcold Vhot Hcold Hhot Hcold Hcold = h+ v+ − v− + K − −k − hamb Vcold Vhot Vcold Vcold
cr
dVhot dt dVcold dt dHhot dt dHcold dt
ip t
626
in which parameters h+ and h− are the (per-volume) supply and return enthalpies.
0 0 0 0 0 0 To simplify the dynamics, we choose an operating point (Vhot , Vcold , Hhot , Hcold ) such that Hhot /Vhot = h− and
629
0 0 Hcold /Vcold = h+ (i.e., that the extra “hot” water at the top of the tank is the same temperature as the “spent” chilled
630
water returning from the building). Linearizing about this point yields the following energy balances:
an
us
628
M
Vhot Vcold dHhot Hhot Hcold ≈ −h− v+ + h− v− − (K + k)h− 1 + 0 − 0 + Kh+ 1 + 0 − 0 + khamb dt Hhot Vhot Hcold Vcold dHcold Hhot Vhot Vcold Hcold ≈ h+ v+ − h+ v− + Kh− 1 + 0 − 0 − (K + k)h+ 1 + 0 − 0 + khamb dt Hhot Vhot Hcold Vcold Then, noting that we are free to choose a reference point for the enthalpy terms, we set h− = 0 and substitute
632
0 0 h+ = Hcold /Vcold in the second term, which leaves the remaining cold enthalpy balance
te
d
631
633
Ac ce p
H0 dHcold = h+ (v+ − v− ) − (K + k) cold 0 dt Vcold K ≈ h+ (v+ − v− ) − 0 Hcold Vcold
1+
Hcold Vcold − 0 0 Hcold Vcold
+ khamb
0 assuming Vcold ≈ Vcold and k ≈ 0. Making the change of variables s := Hcold and y := h+ (v+ − v− ), we arrive at
ds = −˜ σs + y dt
(21)
634
0 with σ ˜ := K/Vcold . After discretization, (21) becomes the storage model (2) used within the optimization problem.
635
References
636
[1] M. Fadzli Haniff, H. Selamat, R. Yusof, S. Buyamin, F. Sham Ismail, Review of HVAC scheduling techniques
637
for buildings towards energy-efficient and cost-effective operations, Renewable and Sustainable Energy Reviews
638
27 (2013) 94–103.
34
Page 34 of 38
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
ip t
cr
645
storage, Journal of Process Control 24 (8) (2014) 1292–1300.
[5] M. H. Albadi, E. F. El-Saadany, Demand response in electricity markets: An overview, in: 2007 IEEE Power Engineering Society General Meeting, IEEE, 2007, pp. 1–5.
us
644
[4] C. R. Touretzky, M. Baldea, Integrating scheduling and control for economic MPC of buildings with energy
[6] G. P. Henze, Energy and cost minimal control of active and passive building thermal storage inventory, Journal of Solar Energy Engineering 127 (3) (2005) 343–351.
an
643
through energy optimization, Processes 1 (3) (2013) 312–329.
[7] M. Kintner-Meyer, A. F. Emery, Optimal control of an HVAC system using cold storage and building thermal capacitance, Energy and Buildings 23 (1) (1995) 19–31.
M
642
[3] K. Kapoor, K. M. Powell, W. J. Cole, J. S. Kim, T. F. Edgar, Improved large-scale process cooling operation
[8] G. P. Henze, C. Felsmann, G. Knabe, Evaluation of optimal control for active and passive building thermal storage, International Journal of Thermal Sciences 43 (2) (2004) 173–183.
d
641
energy planning, Energy and Buildings 58 (2013) 281–291.
[9] M. Behl, T. X. Nghiem, R. Mangharam, Green scheduling for energy-efficient operation of multiple chiller plants, Real-Time Systems Symposium (RTSS), 2012 IEEE 33rd (2012) 195–204.
te
640
[2] L. Duanmu, Z. Wang, Z. J. Zhai, X. Li, A simplified method to predict hourly building cooling load for urban
[10] H. Dagdougui, R. Minciardi, A. Ouammi, M. Robba, R. Sacile, Modeling and optimization of a hybrid system
Ac ce p
639
for the energy supply of a green building, Energy Conversion and Management 64 (2012) 351–363. [11] A. Ashouri, S. S. Fux, M. J. Benz, L. Guzzella, Optimal design and operation of building services using mixedinteger linear programming techniques, Energy 59 (2013) 365–376. [12] M. Trifkovic, W. A. Marvin, P. Daoutidis, M. Sheikhzadeh, Dynamic real-time optimization and control of a hybrid energy system, AIChE Journal 60 (7) (2014) 2546–2556. [13] S.-I. Gustafsson, Mixed integer linear programming and building retrofits, Energy and Buildings 28 (2) (1998) 191–196. [14] D. Zhang, N. Shah, L. G. Papageorgiou, Efficient energy consumption and operation management in a smart building with microgrid, Energy Conversion and Management 74 (2013) 209–222. [15] Y. Ma, F. Borrelli, B. Hencey, B. Coffey, S. Bengea, P. Haves, Model predictive control for the operation of building cooling systems, IEEE Control Systems Technology 20 (3) (2012) 796–803. 35
Page 35 of 38
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
ip t
cr
673
Energy and Buildings 110 (2016) 94–107.
[19] J. E. Braun, A near-optimal control strategy for cool storage systems with dynamic electric rates, HVAC&R Research 13 (4) (2007) 557–580.
us
672
[18] C. R. Touretzky, M. Baldea, A hierarchical scheduling and control strategy for thermal energy storage systems,
[20] J. E. Braun, G. T. Diderrich, Near-optimal control of cooling towers for chilled-water systems, ASHRAE Transactions 96 (CONF-9006117).
an
671
I: Simulation environment, Energy Conversion and Management 64 (2012) 499–508.
[21] J. E. Braun, A general control algorithm for cooling towers in cooling plants with electric and/or gas-driven chillers., HVAC&R Research 13 (4) (2007) 581–598.
M
670
[17] A. Hajiah, M. Krarti, Optimal control of building storage systems using both ice storage and thermal mass—part
[22] M. Ali, V. Vukovic, M. H. Sahir, G. Fontanella, Energy analysis of chilled water system configurations using simulation-based optimization, Energy and Buildings 59 (2013) 111–122.
d
669
building cooling system, Applied Energy 111 (2013) 1032–1045.
[23] A. Aravelli, S. S. Rao, Energy optimization in chiller plants: A novel formulation and solution using a hybrid optimization technique, Engineering Optimization 45 (10) (2013) 1187–1203.
te
668
[16] J. A. Candanedo, V. R. Dehkordi, M. Stylianou, Model-based predictive control of an ice storage device in a
[24] A. J. Ardakani, F. F. Ardakani, S. H. Hosseinian, A novel approach for optimal chiller loading using particle
Ac ce p
667
swarm optimization, Energy and Buildings 40 (12) (2008) 2177–2187. [25] W.-S. Lee, Y.-T. Chen, Y. Kao, Optimal chiller loading by differential evolution algorithm for reducing energy consumption, Energy and Buildings 43 (23) (2011) 599–604. [26] Z. W. Geem, Solution quality improvement in chiller loading optimization, Applied Thermal Engineering 31 (10) (2011) 1848–1851.
[27] L. Lu, W. Cai, Y. C. Soh, L. Xie, S. Li, HVAC system optimization—condenser water loop, Energy Conversion and Management 45 (4) (2004) 613–630. [28] K. M. Powell, W. J. Cole, U. F. Ekarika, T. F. Edgar, Optimal chiller loading in a district cooling system with thermal energy storage, Energy 50 (2013) 445–453. [29] G. P. Henze, B. Biffar, D. Kohn, M. P. Becker, Optimal design and operation of a thermal storage system for a chilled water plant serving pharmaceutical buildings, Energy and Buildings 40 (6) (2008) 1004–1019. 36
Page 36 of 38
699
700
701
702
703
704
electricity pricing, Energy and Buildings 60 (2013) 199–209.
ip t
698
[31] M. Avci, M. Erkoc, A. Rahmani, S. Asfour, Model predictive HVAC load control in buildings using real-time
[32] J. Ma, J. Qin, T. Salsbury, P. Xu, Demand reduction in building energy systems based on economic model predictive control, Chemical Engineering Science 67 (1) (2012) 92–100.
cr
697
in: American Control Conference (ACC), 2012, 2012, pp. 3669–3674.
[33] R. Z. Freire, G. H. Oliveira, N. Mendes, Predictive controllers for thermal comfort optimization and energy savings, Energy and Buildings 40 (7) (2008) 1353–1365.
us
696
[30] D. I. Mendoza-Serrano, D. J. Chmielewski, Controller and system design for HVAC with thermal energy storage,
[34] S. C. Bengea, A. D. Kelman, F. Borrelli, R. Taylor, S. Narayanan, Implementation of model predictive control for an HVAC system in a mid-size commercial building, HVAC&R Research 20 (1) (2014) 121–135.
an
695
[35] M. J. Risbeck, C. T. Maravelias, J. B. Rawlings, R. D. Turney, Cost optimization of combined building heat-
706
ing/cooling equipment via mixed-integer linear programming, in: American Control Conference, Chicago, IL,
707
2015, pp. 1689–1694.
710
711
d
MILP formulation, Computers & Chemical Engineering 17 (2) (1993) 211–227. [37] C. C. Pantelides, Unified frameworks for optimal process planning and scheduling, in: Proceedings on the second
te
709
[36] E. Kondili, C. Pantelides, R. Sargent, A general algorithm for short-term scheduling of batch operations—I.
conference on foundations of computer aided operations, 1994, pp. 253–274.
Ac ce p
708
M
705
712
[38] C. A. M´endez, J. Cerd´a, I. E. Grossmann, I. Harjunkoski, M. Fahl, State-of-the-art review of optimization meth-
713
ods for short-term scheduling of batch processes, Computers & Chemical Engineering 30 (6-7) (2006) 913–946.
714
[39] C. T. Maravelias, General framework and modeling approach classification for chemical production scheduling,
715
AIChE Journal 58 (6) (2012) 1812–1828.
716
[40] I. Harjunkoski, C. T. Maravelias, P. Bongers, P. M. Castro, S. Engell, I. E. Grossmann, J. Hooker, C. M´endez,
717
G. Sand, J. Wassick, Scope for industrial applications of production scheduling models and solution methods,
718
Computers & Chemical Engineering 62 (2014) 161–193.
719
720
721
722
[41] T.-S. Lee, K.-Y. Liao, W.-C. Lu, Evaluation of the suitability of empirically-based models for predicting energy performance of centrifugal water chillers with variable chilled water flow, Applied Energy 93 (2012) 583–595. [42] W. J. Cole, T. F. Edgar, A. Novoselac, Use of model predictive control to enhance the flexibility of thermal energy storage cooling systems, in: American Control Conference (ACC), 2012, 2012, pp. 2788–2793. 37
Page 37 of 38
726
727
728
729
730
[44] D. I. Mendoza-Serrano, D. J. Chmielewski, HVAC control using infinite-horizon economic MPC, in: Decision and Control (CDC), 2012 IEEE 51st Annual Conference on, 2012, pp. 6963–6968.
ip t
725
Energy and Buildings 33 (4) (2001) 371–378.
[45] Z. Ma, S. Wang, X. Xu, F. Xiao, A supervisory control strategy for building cooling water systems for practical and real time applications, Energy Conversion and Management 49 (8) (2008) 2324–2336.
cr
724
[43] B. C. Ahn, J. W. Mitchell, Optimal control development for chilled water plants using a quadratic representation,
[46] G.-Y. Jin, W.-J. Cai, L. Lu, E. L. Lee, A. Chiang, A simplified modeling of mechanical cooling tower for control and optimization of HVAC systems, Energy Conversion and Management 48 (2) (2007) 355–365.
us
723
[47] Y. Ma, F. Borrelli, B. Hencey, A. Packard, S. Bortoff, Model predictive control of thermal energy storage in build-
732
ing cooling systems, in: Decision and Control, 2009 held jointly with the 2009 28th Chinese Control Conference.
733
CDC/CCC 2009. Proceedings of the 48th IEEE Conference on, 2009, pp. 392–397.
737
738
739
740
741
742
743
744
745
M
ˇ acˇ ekov´a, Building modeling as a crucial part [49] S. Pr´ıvara, J. Cigler, Z. V´anˇ a, F. Oldewurtel, C. Sagerschnig, E. Z´ for building predictive control, Energy and Buildings 56 (2013) 8–22.
d
736
Energy and Buildings 108 (2015) 454–462.
[50] A. Imamoto, B. Tang, Optimal piecewise linear approximation of convex functions, in: Proceedings of the world
te
735
[48] K. J. Kircher, K. Max Zhang, On the lumped capacitance approximation accuracy in rc network building models,
congress on engineering and computer science, 2008, pp. 1191–1194.
Ac ce p
734
an
731
[51] C. Frenzen, T. Sasao, J. T. Butler, On the number of segments needed in a piecewise linear approximation, Journal of Computational and Applied Mathematics 234 (2) (2010) 437–446. [52] S. Rippa, Adaptive approximation by piecewise linear polynomials on triangulations of subsets of scattered data, SIAM Journal on Scientific and Statistical Computing 13 (5) (1992) 1123–1141. [53] A. Toriello, J. P. Vielma, Fitting piecewise linear continuous functions, European Journal of Operational Research 219 (1) (2012) 86–95.
38
Page 38 of 38