A Hierarchical Decomposition Approach for the Optimal Design of a
District Cooling System
Bingqian Liu
1 a
, C
ˆ
ome Bissuel
1 b
, Franc¸ois Courtot
2
, C
´
eline Gicquel
3 c
and Dominique Quadri
3
1
EDF R&D, Chatou, France
2
EDF R&D China, China
3
Universit
´
e Paris Saclay, LRI, France
Keywords:
Local Energy System Design, Mixed-integer Linear Programming, Branch & Bound, Hierarchical
Decomposition.
Abstract:
A district cooling system is a centralized cooling supply system providing air conditioning to a set of buildings
located in the same district. Designing and sizing such a system is very complex, as both the initial construction
cost and the operation cost of the cooling system during its entire life must be considered. We first propose a
modeling approach aiming at formulating this combinatorial optimization problem as a mixed-integer linear
program of tractable size. We then extend a previously published hierarchical decomposition technique in
order to find the optimal solution in an efficient way. Finally, we provide preliminary computational results
based on a real-life case study located in China.
1 INTRODUCTION
A district cooling system (DCS) is a centralized cool-
ing supply system. It consumes electricity to cool
down water and distributes it through an underground
pipe network to the buildings in the district to provide
them with air conditioning. DCSs usually are highly
energy-efficient cooling systems. Thus, according to
the Electrical and Mechanical Services Department
of Hong Kong EMSD (2020), using DCSs instead
of traditional air-cooled air-conditioning systems re-
sults in energy savings of 35%. DCSs also compare
well with individual water-cooled air-conditioning
systems using cooling towers as the energy savings
may be up to 20%. Furthermore, this lower energy
consumption leads to lower greenhouse gas emis-
sions, which helps reduce the environmental impact
of the system.
Designing a DCS involves choosing the type and
number of chillers to be installed as well as the ice
storage capacity. These decisions should take into
account not only the construction costs, but also the
operation costs of the system during its whole life-
time. Computing these operation costs is a chal-
a
https://orcid.org/0000-0001-7493-4277
b
https://orcid.org/0000-0002-5430-3168
c
https://orcid.org/0000-0002-2719-7443
lenging problem. Namely, the demand for cooling
power is highly variable and features a daily, weekly
and yearly seasonality together with random varia-
tions. Moreover, due to technical reasons owing to
the chillers, these operations costs are not at all pro-
portional to the produced cooling power. Thus, in or-
der to accurately estimate them, a detailed schedule
describing, on a hourly basis, the on/off status and the
load allocation of each chiller should be built for an
horizon spanning a whole year. Furthermore, the de-
ployment of a district cooling system is usually not
a one-shot decision but rather a process in which in-
vestment decisions are made step by step, following
the development of the district and the upward trend
of the average demand over the years. This implies
that a multi-phase strategic deployment plan should
be built.
This optimization problem can be formulated as a
mixed-integer program. However, its resolution poses
several difficulties. The first one comes from the non-
linearity of the chillers’ performance curves. These
performance curves give, for each chiller, the amount
of electricity consumed as a function of the amount
of produced cooling power and thus play a key role
in the estimation of the system operation costs. Sec-
ond, the need to simultaneously build a multi-year
phasing plan and a detailed operational schedule for
each day of the planning horizon leads to the formu-
Liu, B., Bissuel, C., Courtot, F., Gicquel, C. and Quadri, D.
A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System.
DOI: 10.5220/0010220403170328
In Proceedings of the 10th International Conference on Operations Research and Enterprise Systems (ICORES 2021), pages 317-328
ISBN: 978-989-758-485-5
Copyright
c
2021 by SCITEPRESS Science and Technology Publications, Lda. All rights reserved
317
lation of a huge mathematical program, which can-
not be solved directly by current mathematical pro-
gramming solvers. Finally, the use of classical de-
composition methods, such as the Benders’ decom-
position approach, is not straightforward as it would
imply sub-problems involving binary and/or integer
decision variables. The mixed-integer program mod-
eling the problem is thus computationally intractable
as such.
However, solution approaches exploiting the nat-
ural hierarchy between decisions relative to the sys-
tem design and decisions relative to the daily op-
eration schedules have been investigated in several
works. Weber et al. (2007) thus studied the optimal
design of a multi-energy system. A two-level opti-
mization method is implemented. The master opti-
mization level explores the set of possible system de-
signs with the use of an evolutionary algorithm. For
each considered system design, the slave optimiza-
tion level calculates the optimal cost by linear pro-
gramming. This approach does not provide a guaran-
teed optimal solution and the formulation of the slave
sub-problems as linear programs does not allow to ac-
curately model the way the system operates in prac-
tice. Iyer and Grossmann (1998) also use a bi-level
method to optimize the choice and sizes of equipment
for utility systems. At the design optimization level,
the problem aims at fixing the system infrastructure.
All binary decision variables relative to the opera-
tion schedules are removed from this master problem.
Each time a potential infrastructure is found, the op-
eration optimization level solves a set of single-period
scheduling sub-problems taking the current system
infrastructure as input data. Design cuts are used to
tighten the gap between the solutions found at both
levels. Yokoyama et al. (2015) consider local energy
supply systems and propose a customized Branch &
Cut algorithm exploiting the hierarchical relationship
between the design and operation decision variables
of the mathematical program. The upper level prob-
lem corresponds to the initial optimization problem
in which all operational integer and binary variables
are kept but relaxed to be continuous. Each time an
integer feasible solution is found at the upper level,
a sequence of single-period independent operation
sub-problems is solved to check the feasibility and
value of the current design solution. Note that Iyer
and Grossmann (1998) and Yokoyama et al. (2015)
both consider a single-phase variant of the problem,
i.e. a variant in which all design decisions are one-
shot decisions. Moreover, their design infrastruc-
ture did not allow short-term intra-day energy stor-
age. This allows them to consider single-period (i.e.
one-hour) scheduling sub-problems rather than multi-
period (i.e. 24-hour) scheduling sub-problems, which
significantly decreases the size of these sub-problems.
In the present work, we propose a solution ap-
proach for the optimal design, over a multi-phase in-
vestment horizon, of a local district cooling system
in which intra-day ice storage is allowed. This ap-
proach relies on three key elements. First, we seek
to reduce the size of the initial optimization problem.
We thus consider a deployment plan involving a lim-
ited number of phases, some of which spanning sev-
eral years. Moreover, we use the clustering approach
described in Zatti et al. (2019) to select a small set
of typical days to represent with the smallest pos-
sible loss of accuracy the various conditions under
which the system will be operated. Second, we build
a piecewise linear approximation of the performance
curves of each chiller and propose a way to exploit
their convexity to reduce the size of the formulation
of the operation scheduling sub-problems. This re-
sults in the formulation of a large-size mixed-integer
linear program (MILP). Thirdly, we develop an ex-
act solution algorithm based on the hierarchical de-
composition method recently proposed by Yokoyama
et al. (2015) to solve this MILP.
The remaining of the paper is organized as fol-
lows. Section 2 gives a detailed description of the
problem under study. In Section 3, the approximation
of the nonlinear performance curves and the cluster-
ing method used to identify typical days are first pre-
sented. The formulation of the optimization problem
as an MILP is then provided. Section 4 introduces the
hierarchical decomposition algorithm used to solve
this MILP. In Section 5, the results of our preliminary
computational experiments based on a real-life case
study in China are reported.
2 PROBLEM PRESENTATION
Before the actual construction of a cooling system,
a thorough and reliable system design is necessary.
The system design consists in selecting the appropri-
ate chillers and in sizing the ice storage capacity to
be installed. The main objective is to design a sys-
tem which will be able to meet the clients’ cooling
demand at all time while leading to the lowest long-
term investments and operation costs.
2.1 Resources
A chiller is a machine that removes heat from a liq-
uid by using a variety of techniques such as vapor-
compression. We consider here electric-powered
chillers which are used to cool water. These chillers
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
318
can be classified into two main categories depend-
ing on their functionality. Standard chillers (de-
noted by STDC) only produce cooling power to sat-
isfy the instantaneous demand of the customers. Ice
chillers (denoted by ICEC) have two distinct operat-
ing modes: they either produce cooling power, but
usually with a lower efficiency than the one of a
STDC, or they produce ice. This ice can be stored
for a few hours in an ice storage tank and be melted
afterwards to provide cooling power.
Within each category of chillers (STDC or ICEC),
there are chillers with different production capacity
levels. Each level corresponds to a predefined produc-
tion range, i.e. to a minimum and maximum cooling
power (or ice) it can provide per hour when turned
on, and to a set of performance curves. Namely, the
energy consumption of a chiller is a function of the
produced cooling power and the ambient tempera-
ture. We thus have performance curves representing
the relation between the electric power consumed by a
chiller and the cooling power (or ice) produced under
different ambient temperatures. Note that these per-
formance curves are usually not linear. In the present
paper, we will focus on the case of non-linear convex
performance curves. The more general non-convex
case is left for future work.
Another key resource in the system is the ice stor-
age tank. This tank is linked to all the ice chillers of
the system, can store the produced ice for a few hours
and release it afterwards to produce cooling power. It
has a maximal storage capacity, which can be chosen
within a predefined range.
See Figure 1 for a schematic representation of the
studied DCS.
Figure 1: Overview of a DCS.
2.2 Cooling Demand
As highlighted in the introduction, the cooling de-
mand to be satisfied by the local cooling system is
highly variable. First, the demand for cooling power
varies throughout the day and is usually much higher
at daytime than at night. There are also weekly and
yearly variations: the demand pattern of a weekday
thus significantly differs from the one observed dur-
ing the week-end and the total daily demand varies
during the year, in particular with the summer and
winter seasons. Finally, during the first operating
years of the system, the demand displays a general
upward trend as new clients join the district cooling
system. After this initial period, the demand usually
gets more or less stable for the rest of the system life-
time.
Parts of these demand variations display a sea-
sonal pattern and are thus to some extent predictable.
However, the demand is also subject to the presence
of extreme weather conditions which are difficult to
anticipate. Days with such effects should be consid-
ered separately from others. This translates in partic-
ular into the existence of extreme days, i.e. days in
which the total daily or hourly demand is exception-
ally low or high as compared to its usual value.
A key requirement is that the system should be
able to satisfy the demand for cooling power at all
time, whatever the hour of the day, the day in week or
the season. In particular, it should be able to satisfy
all the demand, even if it is exceptionally low or high.
2.3 Electricity Supply
The chillers are powered by electricity bought from
an external utility provider.
Similar to the cooling power demand, the elec-
tricity price displays daily, weekly and yearly varia-
tions. In particular, within a day, three types of peri-
ods, termed peak, flat and valley periods, can be dis-
tinguished. They correspond respectively to the high-
est, intermediate or lowest prices. Peak periods are
usually at noon and in the evening, the valley peri-
ods at night and early in the morning and the rest of
the day corresponds to flat periods. These price varia-
tions can be exploited to reduce the total energy cost,
e.g. by producing ice at night when the demand for
cooling power is low and the electricity rather cheap,
storing it for a few hours and releasing ice to produce
cooling power at daytime when the demand is high
and the electricity more expensive.
Moreover, in many cases, the contract with the
electricity provider includes a maximum allowed in-
stantaneous power consumption from the grid. This
upper limit is set contractually at the beginning of
each year. The corresponding cost, called the con-
tract fee, is proportional to the subscribed maximum
power. This contract fee differs from the electricity
consumption cost. It namely buys the permission of
A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System
319
consuming electric power and sets an upper limit (ex-
pressed in kW) to this instantaneous consumption. In
contrast, the electricity consumption cost depends on
the total amount of consumed electric energy which
is expressed in kWh.
2.4 Costs
The total cost of the system comprises two main parts:
the design cost and the operation cost.
The design cost is the sum of the purchase, instal-
lation and maintenance costs of the chillers and ice
storage tank, and of the annual contract fee. Once the
cooling system design is chosen, this cost will be de-
termined.
The operation cost corresponds to the cost of the
electricity consumed when the system operates to sat-
isfy the customers’ cooling demand. Due to the non-
linearity of the chillers’ performance curves, this cost
is not proportional to the amount of produced cooling
power. In order to accurately estimate it, we need to
build a detailed production schedule for each day of
the planning horizon using a fine time discretization.
The objective of this optimization problem is to
find a system design which can minimize the sum of
the design and operation costs. Therefore, the best
system design as well as the most cost-efficient daily
optimal operation strategy corresponding to the cho-
sen system should be searched.
3 PROBLEM MODELING
3.1 Selection of Typical Days
The huge size of the optimization problem, which
combines both long-term design and detailed (usually
with an hourly timestep) production scheduling deci-
sions, makes it computationally intractable. A pos-
sible way of getting over this difficulty consists in
selecting, within the available data, a subset of typ-
ical days and extreme days which will represent the
various conditions under which the system will be
operated. These days should be carefully chosen as
they will have a strong impact on the selection of the
chillers and sizing of the ice storage tank.
This can be done by solving a clustering problem
such as the k-medoid problem. In our case, this prob-
lem consists in partitioning the set of days belonging
to the available historical time series into groups or
clusters and in choosing one member (i.e. one day)
in each cluster to represent it so as to minimize the
total Euclidian distance between each day in the time
series and its representative. In this work, we use the
extension of the k-medoid problem proposed by Zatti
et al. (2019) to select typical days for optimizing the
design of our system.
3.2 Piecewise Linear Approximation of
the Performance Curves
The performance curves of a given type of chiller are
supplied by its manufacturer. They give, for a dis-
crete set of values of the ambient temperature, the
amount of electric energy consumed as a function of
the output cooling power. Their non-linearity is a ma-
jor source of difficulty for the resolution of our opti-
mization problem. We thus propose to build a piece-
wise linear approximation of each curve.
This approximation is built using the following
procedure. We first fix a predefined number of break-
points in the piecewise linear approximation to be
built. We then determine the coordinates of these
breakpoints by heuristically solving a small non-
linear optimization problem aiming at minimizing the
distance between the approximate performance curve
and the actual one.
Then, for each hour of each selected day, we use
historical data about the average ambient temperature
to determine the approximate performance curve that
should be used to compute the energy consumption of
each type of chiller during each time period.
3.3 Notation
The problem modeling relies on a three-level time dis-
cretization. The forecast lifetime of the system is first
divided into a set of phases or investment periods,
each one typically spanning one or several years: we
assume that design decisions such as the installation
of a new chiller can be only made at the beginning
of a new phase. Let Φ be the number of considered
phases. Within each considered phase, the various
conditions under which the system will be operated
will be represented by a number of preselected typ-
ical days and extreme days. We denote by D
φ
the
set of typical days and extreme days used to repre-
sent the various daily demand patterns during phase
φ {1, ...,Φ}. Each selected day d of phase φ has a
weight w
φ,d
corresponding to the number of days it
represents, i.e. to the number of days of the original
historical time series which were assigned to the clus-
ter it belongs to. Finally, in order to describe the intra-
day variations of the cooling demand and electricity
price, each selected day is divided into 24 one-hour
periods. For the sake of readability, in what follows,
we use the letter t to represent the time period (φ,d,h)
corresponding to phase φ, selected day d and hour h.
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
320
Let Dem
t
(resp. EP
t
) represent the cooling power de-
mand (resp. the electricity price) at time period t.
There are different types of chillers that may be
installed in the system. Each type of chiller m = (p, l)
can be described by its category p {ST DC, ICEC}
(i.e. standard or ice chiller), which defines the list of
commodities C
p
{COLD, ICE} it can produce, and
its production capacity level l {1, ..., L
p
}. Let M =
{(p, l)|p {ST DC, ICEC}, l {1,...,L
p
}} be the set
of all chiller types. We denote by P
min
m,c
and P
max
m,c
the
minimum and maximum output power of a chiller of
type m = (p, l) producing commodity c C
p
. As ex-
plained above, the performance curve of a chiller of
type m producing commodity c at time period t is ap-
proximately represented by a piecewise linear func-
tion comprising B
t,m,c
breakpoints. Let a
t,m,c,b
and
o
t,m,c,b
be the abscissa and ordinate of breakpoint b
of this piecewise linear function.
With respect to the design of the system, some
restrictions have to be taken into account. There is
namely an upper bound SD
max
m
to the total number of
chillers of type m that may be installed in the sys-
tem. The total ice storage capacity built in the system
should also stay below a maximum allowed capacity
StoC
max
. Moreover, in terms of design costs, we de-
note by FC
φ,m
the fixed cost of investing in a chiller
of type m at phase φ. This fixed cost comprises the
construction or capital expenditure cost at phase φ and
the total maintenance cost of the chiller over its whole
lifetime (from phase φ to phase Φ ). Regarding the
ice storage capacity, the installation cost is broadly
proportional to the installed capacity. We denote by
LC
φ
the cost of building one unit of ice storage ca-
pacity at phase φ. Finally, the unit subscription cost
for the maximum instantaneous power is assumed de-
noted by SC
φ
. Note that all the design costs FC
φ,m
,
LC
φ
and SC
φ
are assumed uniform along each deploy-
ment phase.
3.4 Variables
In order to model the problem as a mixed-integer lin-
ear program, we introduce two sets of decision vari-
ables.
The first set of decision variables corresponds to
long-term design decisions which will determine the
general structure of the system, together with the
multi-year phasing plan. Thus, SD
φ,m
represents the
integer number of chillers of type m to be installed at
the beginning of phase φ, StoC
φ
is a continuous vari-
able representing the ice storage capacity built at the
beginning of phase φ and C
φ
is the maximum allowed
electric power consumption contracted with the elec-
tricity provider for phase φ. These variables will be
referred to as design decision variables in what fol-
lows.
The second set of decision variables are used to
build the operational schedule for each selected day
d D
φ
of each phase φ. In the present work, we will
focus on the special case in which all performance
curves of the chillers are convex. This assumption
allows us to use aggregate performance curves (see
Appendix for more detail about this) and to build the
schedule while considering, in each period t, aggre-
gate scheduling variables, i.e. variables pertaining
to the set of chillers of identical type producing the
same commodity c during t, instead of disaggregate
scheduling variables, i.e. variables pertaining to the
status and output of each individual chiller installed
in the system in period t.
We thus introduce S
t,m,c
the integer number of
chillers of type m producing commodity c during pe-
riod t, P
t,m,c
the total amount of commodity c pro-
duced by the chillers of type m during period t and
Q
t,m,c
the total electric consumption of the chillers of
type m producing commodity c in period t. More-
over, in order to monitor the ice storage and release,
we define STO
t
the amount of ice stored in the tank
at the beginning of period t and R
t
the amount of ice
released during t. These decisions will be referred to
as operation decision variables in what follows.
3.5 Objective Function
We seek to minimize the sum of the design and op-
eration costs of the system over its whole lifetime.
The objective function of the mathematical program
is thus given by:
min
Φ
φ=1
α
φ
"
mM
FC
φ,m
SD
φ,m
+ LC
φ
StoC
φ
+SC
φ
C
φ
+
dD
φ
w
φ,d
23
h=0
mM
cC
p
EP
φ,d,h
Q
φ,d,h,m,c
#
(1)
where α
φ
is the actualization rate for the costs in-
curred in phase φ.
3.6 Constraints
Similar to the decision variables, the constraints of the
mathematical model can be classified into two groups:
design constraints and operation constraints.
3.6.1 Design Constraints
Design constraints are constraints involving only de-
sign decision variables. In the present case, for each
A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System
321
type of chiller m M , we have the following con-
straint which states that the total number of chillers of
type m included in the DCS should be less than the
maximum number of chillers of this type allowed.
Φ
ϕ=1
SD
ϕ,m
SD
max
m
(2)
Similarly, Constraint (3) makes sure that the to-
tal ice storage capacity of the system stays below the
maximum allowed value.
Φ
ϕ=1
StoC
ϕ
StoC
max
(3)
3.6.2 Operation Constraints
Operations constraints are constraints involving op-
eration decision variables, together with design vari-
ables in some cases. For each time period t, we have
the following set of constraints.
Demand Satisfaction. The demand for cooling
power must be satisfied at all time. This can be done
either by directly using the cooling power produced
by the currently turned on chillers and by releasing
some ice from the ice storage tank.
mM
P
t, m,COLD
+ R
t
= Dem
t
(4)
Chillers. Regarding the chillers, two sets of con-
straints should be introduced in the formulation.
First, for each type of chiller m M , the total
number of operating (i.e. turned on) chillers should
be less than the number of chillers currently installed
in the DCS:
cC
p
S
t,m,c
φ
ϕ=1
SD
ϕ,m
(5)
Second, for each type of chiller m M and each
commodity c C
p
it can produce, the total amount
produced by the turned on chillers should stay within
the allowed production range:
P
t,m,c
P
max
m,c
S
t,m,c
(6)
P
t,m,c
P
min
m,c
S
t,m,c
(7)
Electric Consumption. As shown in Appendix,
when the chillers’ performance curve are convex, the
aggregate electric consumption of the set of chillers
corresponding to a given type m = (p,l) M and pro-
ducing commodity c C
p
in t belongs to the epigraph
of a piecewise linear function defined by breakpoints
b = 1...B
t,m,c
. We thus have, for each m = (p,l) M ,
each c C
p
and each b {1...B
t,m,c
1}, the follow-
ing inequality:
Q
t,m,c
s
t,m,c,b
P
t, m, c
+ c
t,m,c,b
S
t,m,c
(8)
where s
t,m,c,b
=
o
t,m,c,b+1
o
t,m,c,b
a
t,m,c,b+1
a
t,m,c,b
is the slope of the
line segment between breakpoints b and b + 1 and
c
t,m,c,b
= o
t,m,c,b
s
t,m,c,b
a
t,m,c,b
the corresponding
constant value. Note that Constraints (8) accurately
compute the electric consumption of a set of chillers
only if the corresponding performance curve is con-
vex. The more general case of non-convex perfor-
mance curve is left for future work.
Moreover, the total amount of electric energy con-
sumed in period t is limited by an upper bound, C
φ
×
1hour, which represents the maximum instantaneous
power C
φ
contracted with the electricity provider
times the duration of the period (one hour):
mM
cC
p
Q
t,m,c
C
φ
(9)
Ice Storage. Regarding the ice storage, three con-
straints should be considered in each time period.
First, the amount of ice stored in the tank should
not exceed the current storage capacity.
STO
t
φ
ϕ=1
StoC
ϕ
(10)
Second, the amount of ice released during the pe-
riod should be less than the amount stored in the tank
at the beginning of the period.
R
t
STO
t
(11)
Third, the evolution of the ice inventory stored in
the tank should comply with inventory balance equa-
tions. We first consider time periods (φ, d,h) corre-
sponding to h {0...22} as the last hour of the day
h = 23 requires a special treatment. For each period
(φ,d, h) such that h {0...22}, we have:
STO
φ,d,h
+
L
ICEC
l=1
P
φ,d,h,ICEC,l,ICE
R
φ,d,h
= STO
φ,d,h+1
(12)
Constraints (12) state that the amount of ice stored
at the beginning of hour h +1, STO
φ,d,h+1
, is equal to
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
322
the amount of ice already stored at the beginning of
hour h, STO
φ,d,h
, plus the total amount of ice pro-
duced by the ice chillers of various capacity levels
during h minus the ice melted during h. Note that the
loss of ice stored in the tank during a day is assumed
to be negligible.
Moreover, for each time period (φ, d, h) corre-
sponding to h = 23, i.e. to the last hour of the se-
lected day, we have the following inventory balance
equation:
STO
φ,d,23
+
L
ICEC
l=1
P
φ,d,23,ICEC,l,ICE
R
φ,d,23
= STO
φ,d,0
(13)
Namely, in practice, the entering inventory of a
given day in the scheduling horizon is imposed by
the leaving inventory of the previous day. In our
case, we do not consider each individual day of the
scheduling horizon but rather a number of represen-
tative days which will not necessarily occur succes-
sively in practice. We thus impose that the leav-
ing ice inventory of a selected day d, computed as
STO
φ,d,23
+
L
ICEC
l=1
P
φ,d,23,ICEC,l,ICE
R
φ,d,23
should
be equal to the entering inventory of the same selected
day d, STO
φ,d,0
. This might be understood as the fact
that the selected day d will be cyclically repeated w
φ,d
times in the simplified scheduling horizon used in our
optimization problem for phase φ.
4 SOLUTION APPROACH
The mathematical program (1)-(13) formulated in
Section 3 displays a particular structure. We namely
have a set of independent sub-problems. Each of them
corresponds to optimizing the detailed schedule for a
single selected day d of a single phase φ and involves
operation variables and operation constraints relative
only to the corresponding day (φ, d). These indepen-
dent sub-problems are linked together by the design
variables.
A Benders decomposition approach might seem
appropriate for a problem displaying such a structure
in which coupling variables link together a set of inde-
pendent sub-problems. However, its application is not
straightforward here as each sub-problem involves in-
teger variables (namely variables S
t,m,c
) and Benders
decomposition algorithms rely on the strong duality
theory to generate Benders cuts. We thus investigate
another hierarchical decomposition approach recently
proposed by Yokoyama et al. (2015).
4.1 Hierarchical Structure of the
Problem
To better highlight the hierarchical structure of the
problem and more easily explain the decomposition
approach, we will use a compact formulation of the
problem.
Vector x represents in a synthetic way the de-
sign variables relative to all phases, i.e. x =
(SD
1
,StoC
1
,C
1
,...,SD
φ
,StoC
φ
,C
φ
). Vector y
k
rep-
resents all the continuous operation variables relative
to a given selected day k = (φ, d) whereas z
k
stands
for all the binary or integer variables relative to day k.
With this notation, Problem (1)-(13) can be for-
mulated as follows.
minZ = f
0
(x) +
kK
f
k
(y
k
,z
k
) (14)
h(x) 0 (15)
g
k
(x,y
k
,z
k
) 0 k K (16)
x Z
ν
(17)
y
k
R
µ
k K (18)
z
k
Z
λ
k K (19)
In the objective function (14), the term f
0
(x)
corresponds to the design cost whereas the term
f
k
(y
k
,z
k
) computes the operation cost for each se-
lected day k in the set K = {(φ,d), φ = 1...Φ,d
D
φ
}. Constraints (15) correspond to the design con-
straints. Constraints (16) represent in a concise man-
ner all the operations constraints relative to day k
K . Constraints (17)-(19) give the definition domain
of each variable vector. Problem (14)-(19) will be re-
ferred to as the Complete Model (CM) in what fol-
lows.
Note how all design variables are defined as in-
teger variables in (17). This restriction is added to
the problem as it is a necessary condition for the use
of the hierarchical decomposition algorithm proposed
by Yokoyama et al. (2015).
This hierarchical decomposition algorithm uses as
a starting point a relaxation of (CM) in which the de-
sign variables x are kept integer or binary whereas
the binary or integer operation variables z are relaxed.
This gives the following semi-relaxed or upper level
problem denoted by (SRM).
A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System
323
minZ
SRM
= f
0
(x) +
kK
f
k
(y
k
,
˜
z
k
) (20)
h(x) 0 (21)
g
k
(x,y
k
,
˜
z
k
) 0 k K (22)
x Z
ν
(23)
y
k
R
µ
k K (24)
˜
z
k
R
λ
k K (25)
Problem (SRM) thus involves the same number of
variables and constraints as the initial problem (14)-
(19). However, the number of binary and integer vari-
ables is drastically reduced, which should ease its res-
olution.
Furthermore, the hierarchical decomposition algo-
rithm relies on the key observation that in Problem
(CM), when the design of the system is determined,
the problem can be decomposed into a set of small
independent operation sub-problems. Let Problem
OM
k
(x
]
) be the Operation Model relative to the op-
timization of the schedule of day k K , given a fixed
design of the system described by vector x
]
. It is for-
mulated as:
minZ
k
(x
]
) = f
k
(y
k
,z
k
) (26)
g
k
(x
]
,y
k
,z
k
) 0 (27)
y
k
R
µ
(28)
z
k
Z
λ
(29)
4.2 Decomposition Algorithm
The decomposition algorithm proposed by
Yokoyama et al. (2015) exploits the hierarchical
structure described above. Thus, at the upper level,
a relaxed version of the initial problem, i.e. problem
SRM, is solved by a Branch & Cut algorithm. Each
time a potential incumbent solution (x
]
,
˜
y,
˜
z) is found
during this tree search, the corresponding values of
the design variables x
]
are used as input data to solve
a series of independent scheduling sub-problems,
namely problems OM
k
(x
]
),k K . This gives an
accurate estimation of the feasibility and value of the
potential design solution x
]
at the operation level.
If x
]
is found to be feasible and less expensive than
the current incumbent solution, it is accepted as the
new incumbent solution. Otherwise, it is rejected.
When all the branches are searched in the upper level
Branch & Bound search tree, the current incumbent
solution gives the optimal solution of the original
problem. This hierarchical decomposition approach,
provided it converges within the allotted computation
time, thus guarantees the optimality of the found
solution.
More precisely, the algorithm comprises the fol-
lowing main steps. See also Figure 2 for a flow chart
of the decomposition algorithm.
Step 1. A Branch & Bound search tree is carried out
within the feasible space of problem SRM: see
1
on
the flow chart. The branching is done on the design
variables x until there is no open node left, in which
case we end the calculation, or an integer feasible
solution of SRM, denoted by (x
]
,
˜
y,
˜
z), is obtained,
in which case we go to Step 2. The value
˜
Z(x
]
) =
f
0
(x
]
) +
kK
f
k
(
˜
y
k
,
˜
z
k
) provides a lower bound of
the actual cost of the design solution x
]
.
Step 2. Before accepting x
]
as a new incumbent so-
lution for the upper level problem, we check its feasi-
bility and compute its actual cost Z(x
]
).
- We first initialize the value of Z(x
]
) as Z(x
]
)
˜
Z(x
]
).
Then, for each k K :
- We solve Sub-problem OM
k
(x
]
) using a standard
Branch & Cut algorithm: see
2
on the flow chart.
- If OM
k
(x
]
) is unfeasible, x
]
cannot be feasible for
Problem (CM). We stop and go to Step 4.
- Otherwise, we record the optimal integer solution of
OM
k
(x
]
), (y
k
,z
k
), and its optimal cost f
k
(y
k
,z
k
).
- We update the current estimation of the actual cost
of the design solution x
]
by computing Z(x
]
)
Z(x
]
) + f
k
(y
k
,z
k
) f
k
(
˜
y
k
,
˜
z
k
).
- If Z(x
]
) is larger than the incumbent value, x
]
cannot
be an optimal solution of Problem (CM). We stop and
go to Step 4: see
3
on the flow chart.
- If all days in K have been considered, we go to Step
3. Otherwise, we go on with the next day in K .
Step 3. We replace the incumbent solution by (x
]
,y,z)
and the incumbent value by Z(x
]
) and go to Step 4:
see
4
on the flow chart.
Step 4. We reject the current node in Problem SRM to
prevent the solver from taking the semi-relaxed objec-
tive value of the design solution x
]
,
˜
Z(x
]
), as a valid
upper bound to be used for the rest of the Branch &
Bound search at the upper level. This allows us to
guarantee that the cutoff value taken into account to
close nodes in the branching process is the actual cost
Z(x
]
) of the current incumbent design solution x
]
, see
Yokoyama et al. (2015). Then go back to Step 1: see
5
on the flow chart.
Finally, we noted in our preliminary experiments,
that while running the hierarchical algorithm de-
scribed above, we often had to solve multiple times
the same operations problem OM
k
(.). Namely, two
different design solutions x
1
and x
2
may e.g. be sim-
ilar for the first phases and differ only for the last
phases. It may also happen that two different design
solutions x
1
and x
2
differ with respect to the decisions
relative to the first phases of the horizon but give the
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
324
Figure 2: Flow chart of hierarchical decomposition algorithm.
same system infrastructure for the last phases. In both
cases, operations problems OM
k
(x
1
) and OM
k
(x
2
)
will be equivalent for all the selected days k corre-
sponding to the phases in which x
1
and x
2
provide the
same system infrastructure.
In order to avoid this useless computational ef-
fort, we thus modify the implementation of the hi-
erarchical decomposition algorithm. When the op-
eration cost of a newly encountered system infras-
tructure needs to be evaluated for a given phase,
the operation sub-problems relative to this phase are
solved and the corresponding optimal operation cost
k st. dD
φ
f
k
(y
k
,z
k
) is recorded in memory. Then,
over the course of the algorithm, each time the op-
eration cost of a design solution corresponding to the
same system infrastructure in the same phase needs to
be evaluated, we simply use the recorded value with-
out recomputing it from scratch.
5 PRELIMINARY NUMERICAL
RESULTS
In order to assess the performance of the proposed
solution approach, we consider a real-life case study
corresponding to a DCS under construction in China.
The expected lifetime of the cooling system is
30 years. The total annual cooling demand is an-
ticipated to increase in the first three years and stay
stable afterwards. We thus consider Φ = 3 invest-
ment phases. Phases 1 and 2 correspond to the first
two years whereas Phase 3 corresponds to the last 28
years. To show the increase in demand, we compare
in Table 1 the total yearly demand, which is the sum
of demand of one year, and the maximum hourly de-
mand for each phase. New chillers and storage ca-
pacity should be installed at the beginning of the first
three years to guarantee that the cooling supply will
meet the increasing demand. The hourly cooling de-
mand value is predicted by combining historical data
on the cooling consumption in the area and forecasts
on the future number of clients which will connect to
the DCS.
There are L
ST DC
= 3 types of standard chillers and
L
ICEC
= 2 types of ice chillers available. Their maxi-
mum output capacity P
max
p,l,c
and installation cost FC
φ,m
(in millions of CNY) are shown in Table 3. Note that
chillers of type (ST DC,1) and (ST DC, 2) have the
same maximum output capacity P
max
p,l,COLD
but chillers
of type (ST DC, 2) are less efficient and less expen-
sive than the ones of type (ST DC,1). For each type of
chiller and each corresponding commodity, the min-
imum output power P
min
p,l,c
is equal to 0.10P
max
p,l,c
. Fig-
ures 3 and 4 display the performance curves of the
available standard and ice chillers at an ambient tem-
perature of 30
C. We use B
t,m,c
= 4 breakpoints to
build the piecewise linear approximation of the per-
formance curve of the chillers of type m producing
commodity c at time period t. The coordinates are
determined by solving a small non-linear optimiza-
tion problem thanks to a heuristic method belonging
to the numpy package in Python. Finally, for each
type of chiller m, the maximum number of chillers
that can be installed, SD
max
m
is set to 10.
The unit cost of installing ice storage capacity LC
φ
is 222.16CNY per kWh and the unit subscription cost
for the maximum instantaneous power allowed, SC
φ
,
is 276CNY per kW. Since the installed storage capac-
ity and the subscribed contract power are discretized,
a unit of storage capacity is 6000kWh and a unit of
contract power is 3000kW. The installed storage ca-
pacity should be lower than StoC
max
= 200GWh and
the contract power should be no more than 30GW.
A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System
325
The discount rate α
φ
is 8%.
The electricity price features a daily seasonality
but no weekly nor yearly variations. Figure 5 shows
the unit price of electricity as a function of the hour
of the day. As the electricity price does not change
with the day in the week or the season, the typical
days and extreme days are selected so as to represent
as best as possible the variations in the cooling power
demand. For each phase, 30 typical days are selected
using the approach proposed by Zatti et al. (2019) and
4 extreme days are identified: the day with the highest
hourly demand, the one with the highest total demand,
the one with the lowest non-zero hourly demand and
the one with the lowest non-zero total demand. There-
fore, there are a total of 34 × 3 = 102 operation sub-
problems.
This results in the formulation of a large-size
MILP. The model involves 64381 variables, among
which 22053 are integer and 600 are binary, as well as
101422 constraints. This MILP is solved using the hi-
erarchical decomposition algorithm presented in Sec-
tion 4. This algorithm is implemented in Python using
the mathematical programming solver CPLEX12.8.
All the problems are solved using a machine with a
Intel Xeon 2.90GHz processor and 16GB RAM.
Table 4 describes the optimal solution in terms of
the number of chillers and the storage capacity built
in each phase, together with the contract power. Fig-
ure 6 shows the operation of chillers and storage in
one of the typical days. The periods in which the ice
production of (ICEC,2) seems to take negative value
correspond to periods in which ice is produced and
stored in the tank.
Furthermore, we carried out additional experi-
ments in order to get a preliminary assessment of the
computational efficiency of the hierarchical decom-
position algorithm and to evaluate the impact of the
selection of representative days on the infrastructure
design. We thus created 4 additional instances based
on our case study by varying the number of selected
days per phase from 6 to 34. The corresponding re-
sults are provided in Table 2. The first line in the table
indicates the number of selected days for each phase.
Lines 2 to 5 provide indications on the size of the
MILP to be solved: we note that the number of vari-
ables and constraints in the problem increases linearly
with the number of selected days. Lines 6 to 9 corre-
spond to the results obtained while directly solving
Problem (CM) with CPLEX12.8 whereas Lines 10 to
13 correspond to the results obtained with the hierar-
chical decomposition algorithm. For each algorithm,
we report the computation time (which was limited
to a maximum of 2 hours), the total and design cost
of the best solution found within the available time
Table 1: Phases and Demand.
Phase 1 2 3
Length (yr) 1 1 28
Total Yearly Demand (GWh) 54.7 250.0 276.8
Max Hourly Demand (MWh) 13.6 62.2 93.8
limit and the remaining gap (i.e. the relative differ-
ence between the values of the best solution and the
best lower found after 2 hours of calculation).
These results first show that the hierarchical de-
composition algorithm is significantly more efficient
than the standard Branch & Cut algorithm embedded
in CPLEX12.8. Namely, with this algorithm, we were
able to solve to optimality the 5 considered instances,
and this with a computation time divided by a factor
of at least 2.25.
Moreover, results from Table 2 also show that the
system infrastructure found by solving the MILP de-
pends on the number of selected days. This can be
seen among others by the fact that the investment cost
varies with the number of selected days, which indi-
cates that the structure of the system and the deploy-
ment plan vary with this number. We note however
that, when using the hierarchical decomposition ap-
proach, both the total cost and the investment cost
become stable when more than 26 selected days are
considered for each phase. Additional computational
experiments are needed to determine the exact value
of the minimum number of selected days per phase
above which the system design does not change any
more.
Figure 3: Performance curves of STDC at 30
C.
Figure 4: Performance curves of ICEC at 30
C.
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
326
Table 2: Comparison between complete model and hierarchical decomposition.
Instance Number of Days 6 10 18 26 34
Size Number of Variables 11461 19021 34141 49216 64381
Number of Integer Variables 3909 6501 11685 16869 22053
Number of Binary Variables 180 240 360 480 600
Number of Constraints 17938 29854 53710 77566 101422
CPLEX Computation time 1726s >7200s >7200s >7200s >7200s
Objective Value 5.50 × 10
8
5.70 × 10
8
5.75 × 10
8
5.75 × 10
8
5.75 × 10
8
Investment Cost 2.45 × 10
8
2.35 × 10
8
2.36 × 10
8
2.36 × 10
8
2.37 × 10
8
Relative Gap 0.00% 0.02% 0.15% 0.34% 0.42%
Hier. Dec. Computation time 220s 922s 942s 1435s 1831s
Objective Value 5.50 × 10
8
5.70 × 10
8
5.74 × 10
8
5.75 × 10
8
5.75 × 10
8
Investment Cost 2.45 × 10
8
2.35 × 10
8
2.35 × 10
8
2.35 × 10
8
2.35 × 10
8
Relative Gap 0.00% 0.00% 0.00% 0.00% 0.00%
Table 3: Available chillers.
Type P
max
p,l,COLD
P
max
p,l,ICE
FC
m
(STDC,1) 8791 - 13
(STDC,2) 8791 - 12
(STDC,3) 5000 - 7.7
(ICEC,1) 8087 5626 13
(ICEC,2) 5000 3478 8.3
Table 4: Optimal system design and phasing obtained with
30 typical days and 4 extreme days at the operation level.
Phase 1 2 3
(STDC,1) 0 0 0
(STDC,2) 1 5 3
(STDC,3) 0 1 0
(ICEC,1) 0 0 0
(ICEC,2) 1 0 0
Storage Capacity (MWh) 24 0 0
Contract Power (MW) 3 9 15
Figure 5: Electricity Price of a Day.
6 CONCLUSION
We presented a modeling and solving approach for
the optimal design of a district cooling system involv-
ing intra-day energy storage over a multi-phase hori-
zon. The modeling approach relies on a clustering
algorithm to identify a subset of typical days and on a
Figure 6: Operation Result of a Typical Day in Phase 3.
piecewise linear approximation of the chillers’ perfor-
mance curves. It results in the formulation of a large-
size mixed-integer linear program. An improved hier-
archical decomposition method is then implemented
to optimally solve this MILP. This decomposition
method exploits the hierarchy between upper-level in-
frastructure design decision and lower-level operation
scheduling decisions. Our preliminary computational
results carried out on a real-life case study located
in China show that the proposed hierarchical decom-
position algorithm significantly outperforms a mathe-
matical programming solver at providing optimal so-
lutions of the MILP. Moreover, our results also show
that, provided a minimum number of typical days are
taken into account to estimate the operation costs, the
final system infrastructure and the deployment phas-
ing do not change with the subset of selected typical
days, which is an important point to gain the trust of
the decision makers.
There are several possible directions for future re-
search suggested by the present work. First, addi-
tional computational experiments are needed to evalu-
ate the impact of the use of a limited number of repre-
sentative days and of the piecewise linear approxima-
tion of the chillers’ performance curves on the solu-
tion provided by the proposed approach. Second, on
a longer term perspective, it would be interesting to
extend the proposed approach to consider non-convex
A Hierarchical Decomposition Approach for the Optimal Design of a District Cooling System
327
performance curves for the chillers and to study more
general local energy systems, in particular systems si-
multaneously providing a district with several sources
of energy (electricity, heat, cold, ...).
REFERENCES
EMSD (2020). Benefits of dcs. http:
//aiweb.techfak.uni-bielefeld.de/content/
bworld-robot-control-software/. [September 9,
2020].
Iyer, R. and Grossmann, I. E. (1998). Synthesis and op-
erational planning of utility systems for multiperiod
operation. Computers & Chemical Engineering, 22(7-
8):979–993.
Weber, C., Mar
´
echal, F., and Favrat, D. (2007). Design and
optimization of district energy systems. In Computer
aided chemical engineering, volume 24, pages 1127–
1132. Elsevier.
Yokoyama, R., Shinano, Y., Taniguchi, S., Ohkura, M.,
and Wakui, T. (2015). Optimization of energy supply
systems by MILP branch and bound method in con-
sideration of hierarchical relationship between design
and operation. Energy conversion and management,
92:92–104.
Zatti, M., Gabba, M., Freschini, M., Rossi, M., Gambarotta,
A., Morini, M., and Martelli, E. (2019). k-milp: A
novel clustering approach to select typical and ex-
treme days for multi-energy systems design optimiza-
tion. Energy, 181:1051–1063.
APPENDIX
Convex Performance Curves
When the performance curves of the chillers are con-
vex, we have the following property.
Lemma 1. For a set of identical chillers producing
the same commodity, the optimal load allocation con-
sists in equally distributing the total output power be-
tween the chillers.
Proof. The proof is done by contradiction. Suppose
we have two identical chillers, producing a total cool-
ing power of P with chiller 1 producing P
1
and chiller
2 producing P
2
> P
1
. Let π : P 7→ Q = π(P) be
the convex performance curve of these two chillers.
The total amount of electricity consumed by the two
chillers producing P is Q = π(P
1
) +π(P
2
).
We show that this load allocation is not optimal,
i.e. that it is possible to reduce the total amount of
consumed electricity. Namely, let δP be a small varia-
tion in the output. By decreasing the output of chiller
2 by δP, we can obtain a decrease in the electricity
consumed by this chiller of π
0
(P
2
)δP, where π
0
is the
derivative function of π. In order to still be able to
provide a total output of P, we increase the output of
chiller 1 by δP, which leads to an increase in its elec-
tricity consumption of π
0
(P
1
)δP.
By convexity of function f , π
0
(P
1
) < π
0
(P
2
).
Hence the total amount of electricity consumed with
the load allocation (P
1
+ δP, P
2
δP) is smaller than
the one consumed with the load allocation (P
1
,P
2
).
The result follows.
Let us now focus on the case in which the perfor-
mance curve π is convex and piecewise linear. When
multiple identical chillers are simultaneously produc-
ing the same commodity, the relation providing the
total amount of consumed electricity as a function of
the total amount of output power can be plotted as an
aggregate performance curve. We have the following
property:
Lemma 2. Let π be the piecewise linear and convex
performance curve of a given type of chiller. π in-
volves B breakpoints. Let (a
b
,o
b
) be the abscissa and
ordinate of breakpoint b.
The aggregate performance curve Π
S
of S identi-
cal chillers of this type is also piecewise linear and
convex. It involves B breakpoints whose coordinates
are given by (Sa
b
,So
b
).
Proof. Let us consider the case where the total out-
put of the S chillers is P [SP
min
;SP
max
] where
[P
min
,P
max
] is the production range of a single chiller.
By Lemma 1, the optimal load allocation consists
in requiring each chiller γ = 1...S to produce the same
output P
γ
=
P
S
. Let b be the index of the breakpoint
of function f such that P
γ
[a
b
,a
b+1
]. The energy
consumed by each chiller γ is thus given by: Q
γ
=
s
b
P
γ
+ c
b
where s
b
and c
b
are the slope and constant
value of the b
th
line segment of π. The total energy
consumed by the S chillers is thus equal to Q = s
b
P +
c
b
S. This equality holds for any value of P such that
P
S
[a
b
,a
b+1
], i.e. any value of P [Sa
b
,Sa
b+1
]. This
means that Π
S
is linear over the segment [Sa
b
,Sa
b+1
],
with a slope equal to s
b
and a constant value of c
b
S.
By generalizing this result to all possible values
of the total output P, we have that Π
S
is a piecewise
linear function involving B breakpoints of coordinates
(Sa
b
,So
b
). Moreover the slope of Π
S
on its b
th
seg-
ment is s
b
. As π is convex, we have s
b
s
b+1
,b =
1...B. Π
S
is thus convex.
ICORES 2021 - 10th International Conference on Operations Research and Enterprise Systems
328