Optimization models for a single-plant District Cooling System

European Journal of Operational Research 247 (2015) 648–658

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Innovative Applications of O.R.

Optimization models for a single-plant District Cooling System Reem Khir∗, Mohamed Haouari Department of Mechanical and Industrial Engineering, College of Engineering, Qatar University, Doha-Qatar 2713, Qatar

a r t i c l e

i n f o

Article history: Received 7 September 2014 Accepted 30 May 2015 Available online 6 June 2015 Keywords: OR in energy Network design District Energy Systems Compact MIP models Reformulation Linearization Technique (RLT)

a b s t r a c t A District Cooling System (DCS) is an interconnected system encompassing a centralized chiller plant, a Thermal Energy Storage (TES) unit, a piping network, and clusters of consumers’ buildings. The main function of a DCS is to produce and deliver chilled water to satisfy the cooling demand of a scattered set of buildings. DCSs are recognized to be highly energy efficient, and therefore constitute an environment-friendly alternative to the traditional power-driven air conditioning systems being operated at individual buildings. In this paper, we investigate the optimal design and operation of a DCS so that the total investment and operational costs are minimized. This involves optimizing decisions related to chiller plant capacity, storage tank capacity, piping network size and layout, and quantities to be produced and stored during every period of time. To this end, mixed-integer programming (MIP) models, that explicitly capture the structural aspects as well as both pressure- and temperature-related requirements, are developed and tested. The results of computational experiments that demonstrate the practical effectiveness of the proposed models are also presented. © 2015 Elsevier B.V. and Association of European Operational Research Societies (EURO) within the International Federation of Operational Research Societies (IFORS). All rights reserved.

1. Introduction Driven by the anticipated increase in energy demand and knitted together with a rapid rise in CO2 emission levels, the notion of sustainable energy systems is gaining an ever-growing attention worldwide. One example of such systems is the so-called District Cooling System (DCS). The main function of a DCS is to mass-produce cooling requirements for a group of buildings while ensuring tailored provision to each of them, depending on their needs. Hence, DCS appears as an alternative to the power-driven air conditioning systems being operated at individual buildings. Worldwide, the current status quo of energy usage shows that at least 10 percent of electricity is used for cooling purposes (DHC/+ Technology Platform, 2009). This percentage is even much higher in some hot climatic countries such as the Gulf Cooperation Council (GCC) countries, where air conditioning accounts for 50 percent of its annual electricity consumption (Booz & Company, 2012). For that reason, DCSs are recently gaining a remarkable market position across the globe as they can reduce electricity consumption by 40 percent–50 percent comparing to conventional air conditioning systems (Arab Construction World, 2012). Currently, DCSs are adopted in the United States, Japan, Korea and many Western and Eastern European countries such as Aus∗

Corresponding author. Tel.: +97444034107. E-mail addresses: [email protected], [email protected] (R. Khir), [email protected] (M. Haouari). URL: http://faculty.qu.edu.qa/haouari (M. Haouari)

tria, Finland, France, Germany, Italy, Norway, Slovenia, and Sweden (Euro Heat & Power, 2013). Generally, the market of DCS has not matured enough as it emerged quite recently; yet, its adoption trend is on the rise. For instance, a 10-fold growth in its installed capacity was observed in Europe during the last decade (DHC/+ Technology Platform, 2009). Nonetheless, DCSs are recently gaining ground in the Middle East, particularly in the GCC countries due to their desert climate. For example, it was observed that 14 percent of GCC cooling demand was satisfied by means of DCSs in 2010 (Arab Construction World, 2012). Structurally, a DCS is an interconnected system encompassing a centralized chiller plant, a main distribution network and clusters of consumers’ buildings (see Fig. 1). The cooling effect, in form of chilled water, is produced in the plant and then distributed to individual customers through a piping network. In more advanced settings, DCS is complemented with a storage tank that is configured with the chiller plant. This optional installation increases the system’s overall flexibility as cooling requirements are not required to be met by chiller plant at all points of time; rather, excess cooling effect can be produced during non-peak hours and stored for future use, especially during peak hours (Powell, Cole, Ekarika, & Edgar, 2013). All together give a DCS the immense potential to unpack both environmental and economic benefits. Its consolidated production scheme allows tapping into the economy of scale; thus, making it a more energy efficient alternative when compared to building-specific cooling systems (Chow, Au, Yau, Cheng, Chan, & Fong, 2004). Moreover, its flexible operational conditions allows integrating various

http://dx.doi.org/10.1016/j.ejor.2015.05.083 0377-2217/© 2015 Elsevier B.V. and Association of European Operational Research Societies (EURO) within the International Federation of Operational Research Societies (IFORS). All rights reserved.

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

649

District served by a central chiller plant Central Chiller Plant

Chiller

Chiller

Chiller

Storage Tank

Fig. 1. Example of a DCS configuration.

cooling technologies that utilize non-carbon energy sources and environmental friendly cooling mediums such as treated sewage effluent and sea water. Hence, DCSs contribute to the reduction in pollution levels with considerable cut-off in energy-use costs (Hart & Rosen, 1996). Even though a DCS is considered to be a long-term sound alternative, the investment costs required to build its infrastructure may hinder its adoption in some areas. It was found that 60 percent of its investment cost is attributed to the distribution network (ASHRAE Inc., 2009). This suggests that the structural optimization of a DC network is paramount and well justified. Despite this highlighted importance, the current literature dealing with the structural optimization of DCSs is relatively scant. Indeed, and unlike the great attention given to District Heating Systems (DHSs) in the operations research literature (see Aringhieri, 2003; Henning, Amiri, & Holmgren, 2006; Ming, Fenglan, & Hairui, 2012; Pinson, Nielsen, Nielsen, Poulsen, & Madsen, 2009; Rastpour & Esfahani, 2010; Sakawa, Kato, & Ushiro, 2001; 2002; Uhlemair, Karschin, & Geldermann, 2014, to quote just a few), only few studies were devoted to the optimization of DCSs. Clearly, the special characteristics observed in a DCS suggest having a customized optimization model that essentially differentiates it from any other energy distribution system. For instance, the type of cooling medium imposes certain constraints on the piping configuration. This is due to the system’s thermal and hydraulic attributes such as temperature differential and pressure drop, which eventually have impact on both the design of piping network and the overall system performance. One of the earliest attempts in the area of DC networks was pioneered by Chan, Hanby, and Chow (2007) who aimed at finding a near optimal network topology by using a modified version of the genetic algorithm (GA; Chan et al., 2007). In their study, a sim-

plified version of the problem was addressed with the objective of minimizing both the fixed cost associated with building the piping configuration, and the pumping energy costs associated with running the plant operations. In parallel to this, relevant work was done by Soderman (2007) who developed an MIP model to find an optimal design of a DC network (Soderman, 2007). Proposed model sought to determine both structural and operational aspects of the problem so that both the annual operational and annualized investment costs are minimized. However, hydraulic- and temperature-related constraints were ignored. Likewise, Feng and Long (2008) directed their first effort to find a network layout that minimizes the total annual cost (Feng & Long, 2008). Interestingly, cooling energy losses in pipes were regarded as one of the cost parameters in their objective function. In their following study, Feng and Long (2010) utilized a standard GA to find an optimal network configuration while additionally incorporating the hydraulic stability constraints of the system (Feng & Long, 2010). In many cases, the structural optimization of a DCS was intended to reflect the optimization of one part of the system, being the piping configuration only without considering important network design aspects presented at both the upstream and downstream ends of the system. For instance, sizing the chiller plant is an important design parameter as the capacity produced by the plant drives the whole system. This is proved by the problem that was investigated by Soderman (2007) who considered finding an optimal operating capacity for both the chiller plant and the storage tank (Soderman, 2007). Moreover, the importance of ensuring a proper integration between the network and the connected customer substations was not evident in the reviewed published works. To stand out among the preceding few studies in the same area, we address the optimization of a DCS design while capturing

650

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

Fig. 2. Example of a daily cooling demand pattern.

structural constraints as well as both pressure- and temperaturerelated constraints. It was noticed that these technical constraints were overlooked in the majority of proposed models while they are crucially responsible for systems functionality and performance. In this paper, we make the following contributions: We develop MIP models for the optimal design of a DCS by specifying the chiller plant size, the storage tank size, the piping network size and layout, and the quantities produced and stored during each period of time while considering structural and technical constraints (including temperature- and pressure-related ones). The objective is to minimize the sum of investment and operations costs. • We present the results of computational experiments that demonstrate the practical usefulness of the proposed models, since solving large instances can be achieved within moderate CPU times. The remainder of this paper is organized as follows. Section 2 describes the formal problem settings and introduces some notation and assumptions used in building our models. Section 3 presents the MIP models for a DCS design problem. Section 4 shows the results of a computational study that aims at demonstrating the solvability and effectiveness of the proposed models. Finally, Section 5 highlights concluding remarks and future research directions. •

2. Structural optimization of a DCS We assume that a set of n customer buildings is to be connected to a single cooling source. The connection scheme is assumed to be indirect. It means that the customer cooling loads are to be supplied to the Energy Transfer Station (ETS) located at each customer building and not directly to the building air handling units. The cooling demand for each customer building is assumed to be deterministic, periodic (with a one-day periodicity), and stationary. For computational convenience, we assume that a day is divided into T time periods of equal duration (a period generally refers to one hour). We denote by τ the duration of every period t ∈ T. The demand for every customer j ∈ C during every period t ∈ T is denoted by dtj , and is expressed as a flow rate. Fig. 2 below illustrates an example of a daily cooling load pattern.

For modeling purposes, we assume that we are given an undirected graph G = (V, E ). The node set V includes three types of nodes: (i) a supply node (denoted by r), (ii) nodes representing customers’ substations (denoted by C), and (iii) facultative nodes representing pipe junctions (henceforth called Steiner nodes, denoted by S). These latter nodes are optional and may be included in the network to reduce the total investment cost. The edge set E includes all potential links. Given the input undirected graph G = (V, E ), we define a bidirected digraph by replacing each edge e = {i, j} ∈ E with two arcs (i, j) and (j, i) with opposite directions. Hence, we denote by A the set of arcs. The overall piping network is assumed to take the form of a treelike network. Hence, no loops/cycles are allowed. In other words, the network shall not accommodate any redundant pipes while ensuring connectivity. In our case, horizontal installation of pipes in underground trenches is assumed. The supply node represents a single facility that includes a chiller plant and a thermal storage tank. The location of this facility is known and given. However, we need to determine the capacity of both production and storage systems. For that reason, a set of possible capacities is given for each of the chiller plant and the storage tank based on market availability and demand needs. These sets are denoted by K and H, respectively. They are characterized by a fixed cost associated plant denotes the fixed cost assowith building and installing them. FCk ciated with every chiller capacity k (k ∈ K). Similarly, FCtank denotes h the fixed cost associated with every tank capacity h (h ∈ H). In addition to this, there are variable costs for producing and storing a unit prod of chilled water during every period t ∈ T and is denoted by VCt sto and VCt , respectively. An arc (i, j) ∈ A represents a possible direct connection between nodes i and j. If an arc is selected to be part of the main distribution network, then it should be assigned a pipe type for installation out of a set M of available pipe types. For every pipe of type m ∈ M, we are given: (i) λm : a unit investment cost covering all purchase and installation charges; (ii) dm : a pipe size (that is, the inner diameter m : minimum allowable velocity, and (iv) V m : maximum size); (iii) Vmin max allowable velocity. These two latter characteristics enforce lower and upper bounds on the allowable flow rate of chilled water in every pipe min and ϕ max , respectively. type. We denote these latter by ϕm m

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

Central Plant

12

pmax

p j=1

10

Customer2

Customer3

Customer4

Customer5

pj=5 pmin

p j=2

pj=3

8

Pmax

pj=4

Chiller Plant

Pressure Head (Bar)

Customer1

651

δPmin

6 4

P 2 min

0 0.0

2.0

4.0 Distance (km)

6.0

8.0

10.0

Fig. 3. Sample pressure gradient for a primary DC distribution network.

Part of the problem is to find an optimal layout while sizing pipes for installation so that the investment cost is minimized. The cost parameter presenting this part of the problem is the fixed investment cost denoted by Cimj , which is a function of both pipe’s unit cost λm and edge length lij (that is, Cimj = λm × li j ). 2.1. Hydraulic aspects Considering the hydraulic behavior of a DCS is vital to ensure the service continuity and the delivery of chilled water to all customers. Basically, a DCS is subject to two types of pressure changes; namely, pressure drop throughout the piping system, and differential pressure across the customers’ ETS. Both changes depend on the size and scheme of pumping system, along with the size and layout of the distribution network (that is, number, locations, loads and diversity of connected customers, type of connection schemes, and piping size and layout). Fig. 3 below illustrates typical pressure-related changes exhibited in a DCS. In the presented problem, a primary pumping scheme is assumed to be in use. Thus, the pumping head is located at the supply node and has the maximum pressure value in the system, denoted by pmax . Clearly and as the fluid moves down in the network, the pressure starts to drop due to fluid friction (major and minor losses). We denote by pj the pressure value of incoming chilled water at any node j ∈ V. In addition, we denote by Pimj the pressure drop along the edge that connects node i to node j when a pipe of size m is installed. More importantly and as stated earlier, another differential pressure is experienced across customer’s ETS and denoted by δ Pj . This drop is encountered while customer’s heat-exchangers extract needed energy before returning the chilled water through the return pipe. Such hydraulic-changes dictate maintaining certain values of pressure at every node so that the system’s operations remain feasible. Given that the pumping pressure is designed based on the pressure drop through the critical path (Chan et al., 2007), a maximum allowable pressure drop Pmax is defined to ensure delivery of chilled water to all customers. Likewise, a minimum allowable differential pressure δ Pmin is defined to ensure proper integration between the customer’s ETS and the main distribution network while ensuring the service delivery to the most remote customer (critical customer). Both parameters are considered in our model by imposing upper and

lower bounds on pressure values at every node in the network. To this end, two main assumptions were made to mathematically represent the pressure-related binding constraints. First, we assume that the pressure of chilled water decreases linearly as it moves downstream in a path. Secondly, we assume that the technical and structural behaviors of the supply and return lines are symmetric. As can be seen from Fig. 3, the pumping pressure depends on both the mean pres¯ and a safety factor/cushion of sure value along a path, denoted by P, which pressure cannot go below, denoted by πmin . In addition to this, the minimum pressure value in the system, denoted by pmin depends on two key parameters: (i) Pmax , and (ii) δ Pmin , as discussed earlier. Therefore, for any node j ∈ V: •

•

pj shall not exceed a maximum pressure value pmax , known as the pumping pressure. pj shall not drop below pmin which equals to the maximum value of (P¯ + 12 δ Pmin ) and ( pmax − Pmax ).

2.2. Thermal aspects While the DCS hydraulic profile is essential to ensure system’s functionality, service delivery, and compatibility with customer’s ETS, the system’s temperature differential is another key design parameter that has significant impact on the DCS overall performance and economics. Practically, low delta T syndrome is a very critical and common industry problem. Operating low delta T dictates a need to increase the amount of produced flow rates to meet customers’ requirements; thus, results in wasted energy, oversized capacities (in plant and distribution network) and consequently, increased capital and operational costs. This also contributes to experiencing increased heat gains in the system while the chilled water is being pumped and transported through the unnecessarily oversized pipes. Therefore, proper control of system’s thermal behavior is crucial to ensure that customers’ loads are met with the least possible wastes. Fig. 4 below illustrates typical temperature-related changes exhibited in a DCS. We denote by tj the incoming temperature of chilled water at every node j ∈ V. As explained, experiencing heat gains in the system is possible and therefore can result in increasing the chilled water temperature while being transported. This encountered temperature increase is denoted by Tim and is cumulative along the path that j

652

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

Central Plant

Customer1

Customer2

Customer3

Customer4

Customer5

RETURN PIPE GRADIENT

14

treturn

10 Chiller Plant

Temperature (Celsius)

12

Tmin

8 6 4

t j=3

t j=2

t j=1

tj

tmax t j=5

4

Tma T

tmin

SUPPLY PIPE GRADIENT

0.0

2.0

4.0 Distance (km)

8.0

6.0

10.0

Fig. 4. Sample temperature gradient for a primary DC distribution network.

connects the supply node with every customer node j ∈ C and vice versa. A maximum threshold for these heat gains is set depending on the contractual supply temperature as stipulated between the DC utility and the connected customers. Moreover, a minimum allowable temperature differential across customer’s ETS is introduced, and is denoted by δ Tmin . This is necessary to make sure that the system does not operate below-design delta T. Consequently, for any node j ∈ V: •

•

tj shall not drop below tmin , known as the supply temperature from the chiller plant. tj shall not exceed a maximum temperature value tmax , that is controlled by both the guaranteed maximum supply temperature as contracted with customers, and the system return temperature that is needed to meet δ T requirements of the DCS. Therefore, tmax equals to the minimum value of (T¯ − 12 δtmin ) and (tmin + Tmax ).

Our main focus is to develop and validate a comprehensive model that explicitly captures the practical features of the system, especially the ones that play vital role in determining the optimal network structure and operations as highlighted above. More precisely, the problem is to determine: (i) The capacity to be installed at the chiller plant. (ii) The capacity to be installed at the storage tank. (iii) The layout of the main distribution network as well as the sizing of each pipe in the network. (iv) The quantities of chilled water to be produced and stored at each period of time.

Table 1 Indices and sets. V: A: R: C: S: Vi− : Vi+ : M: K: H: T:

Set of all nodes. Set of all arcs (i, j) ∈ A. The source node (chiller plant node), R ⊂ V. Set of customer nodes, C ⊂ V. Set of Steiner nodes, S ⊂ V. Set of predecessor nodes to i. Set of successor nodes to i. Set of available pipe sizes. Set of available chiller plant sizes. Set of available tank sizes. Number of time periods of equal durations.

costs, capacities, and other parameters) are input data and to be collected from appropriate sources (e.g. market prices, service-specific metrics and load estimation). 3.1. Notation and decision variables The indices and sets, parameters, and decision variables, are displayed in Tables 1–3. 3.2. A mixed-integer nonlinear formulation

Minimize



FCplant yk + k

k∈K



FCstorage gh + h

The objective is to minimize the total investment and operational costs.

 t∈T

VCtprod Ft

t∈T

h∈H

+



VCtSto It +

 

cimj xm ij

(1)

(i, j )∈A m∈M

subject to:

3. Problem formulation In this section, a mixed integer linear programming (MILP) model is formulated for a DC design problem as described in Section 2. First, all index sets, notations and parameters which are used in the model are summarized. Second, an MIP model is formulated for a singlesource DC network. Then, the developed model is linearized using the RLT methodology. We assume that all relevant data (demand values,



yk = 1,

(2)

gh ≤ 1,

(3)

k∈K

 h∈H

Ft ≤

 k∈K

Qk yk ,

t ∈ T,

(4)

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

653

Table 2 Parameters. FCplant : k : FCtank h VCtprod : VCtsto : Cimj : Qk : Dh : τ: dtj : ϕmmin : ϕmmax : Timj : tmin : tmax : Pimj : pmin : pmax :

Fixed cost for building and installing chiller plant of capacity k ∈ K. Fixed cost for building and installing storage tank of capacity h ∈ H. Variable cost for producing a unit of chilled water during every period t ∈ T. Variable cost for storing a unit of chilled water during every period t ∈ T. Unit cost of installing a pipe m ∈ M on arc (i, j). Capacity associated with every chiller plant of size k ∈ K. Capacity associated with every storage tank of size h ∈ H. Duration of every period t ∈ T. Demand of every customer j ∈ C during every period t ∈ T. Minimum allowable flow rate of chilled water in every pipe m ∈ M. Maximum allowable flow rate of chilled water in every pipe m ∈ M. Temperature increase in chilled water being transported in pipe m ∈ M installed at arc connecting nodes i and j. Minimum allowable chilled water temperature at every node in the system. Maximum allowable chilled water temperature at every node in the system. pressure drop of chilled water being transported in pipe m ∈ M installed at arc connecting nodes i and j. Minimum allowable pressure value of chilled water at every node in the system. Maximum allowable pressure value of chilled water at every node in the system.

Table 3 Decision variables. Binary variable that takes value 1 if plant is built and installed at its kth capacity, and 0 otherwise, k ∈ K. Binary variable that takes value 1 if storage tank is built and installed at its hth capacity, and 0 otherwise, h ∈ H. Amount of chilled water produced by chiller plant during period t, t ∈ T, expressed as flow rate. Stock level in storage tank at the end of period t, t ∈ T. Binary variable that takes value 1 if arc (i, j) is selected, and 0 otherwise, (i, j) ∈ A. Binary variable that takes value 1 if pipe of type m is installed on arc (i, j), and 0 otherwise, m ∈ M, (i, j) ∈ A. Amount of flow rate in arc (i, j) during period t if pipe of type m is selected for installation on this arc, (i, j) ∈ A, m ∈ M, t ∈ T. Temperature of supplied chilled water at node j, j ∈ V. Pressure of supplied chilled water at node j, j ∈ V.

yk : gh : Ft : It : zij : : xm ij : fitm j tj : pj :

It ≤



Dh gh ,

t ∈ T,

(5)

p j z i j = pi z i j −

h∈H n 

dtj ,

t ∈ T,

(6)

pmin ≤ p j ≤ pmax

(7)

pmin

j=1

Io = IT ,



∀ j ∈ C,

zi j ≤ p j ≤ pmax

i∈V j−

(18)



(19)

zi j

∀ j ∈ S,

(20)

i∈V j−

zi j = 1,

∀ j ∈ C,

(8)

pr = pmax ,

zi j ≤ 1,

∀ j ∈ S,

(9)

Ft , It ≥ 0,

i∈V j−



Pimj xm ∀(i, j ) ∈ A, ij

m∈M

It−1 + τ Ft = It + τ





(21)

t ∈ T,

(22)

i∈V j−



∀(i, j ) ∈ A,

xm i j = zi j ,

(10)

yk , gh ∈ {0, 1},

k ∈ K, h ∈ H.

(23)

(11)

xm i j , zi j , ∈ {0, 1},

∀(i, j ) ∈ A, m ∈ M,

(24)

(12)

fitm j ,≥ 0

m∈M



fitm j −

i∈V j− m∈M



 

t f tm jk = d j ,

∀ j ∈ C, t ∈ T

k∈V j+ m∈M

fitm j =

i∈V j− m∈M

 

f tm jk ,

∀ j ∈ S, t ∈ T

k∈V j+ m∈M

tm max m ϕmmin xm ∀(i, j ) ∈ A, m ∈ M, t ∈ T i j ≤ f i j ≤ ϕm xi j ,

t j zi j = ti zi j +



t j, p j ≥ 0 (13)

Timj xm ∀(i, j ) ∈ A, ij

(14)

∀ j ∈ C,

(15)

m∈M

tmin ≤ t j ≤ tmax tmin



zi j ≤ t j ≤ tmax

i∈V j−

tr = tmin ,



zi j

∀ j ∈ S,

(16)

i∈V j−

(17)

∀(i, j ) ∈ A, t ∈ T, m ∈ M, ∀ j ∈ C ∪ S.

(25)

(26)

The objective function (1) minimizes the sum of the fixed costs of building chiller plant and storage tank, the variable production and storage costs and the fixed costs associated with purchasing and installing distribution pipes. Constraint (2) ensures that only one capacity is selected for chiller plant installation. Constraint (3) indicates that the installation of a storage tank is optional, and is assigned at most one capacity. Constraint (4) ensures that the produced amount of chilled water does not exceed the installed capacity of the chiller plant. Likewise, Constraint (5) ensures that the amount of stored chilled water does not exceed the installed storage tank capacity. Constraint (6) imposes the

654

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

balance constraints for the storage tank. Constraint (7) reflects modeling the steady state; hence, the transient state that occurs when we start operating the plant for the first time is not considered. Constraints (8) and (9) enforce the connectivity (as well as the non redundancy) of the network. They indicate that each customer node has exactly one incoming arc, and each Steiner node has at most one incoming arc, respectively. Since x is a binary variable, then by virtue of (8) and (9), z is necessarily a binary variable and the integrality of the z-variables is relaxed to z ≥ 0. Constraint (10) is introduced to ensure that only one pipe type is installed at any arc. Constraints (11) and (12) are the flow conservation conditions which must hold for all non-root nodes. Note that Constraint (11) holds only for customer nodes. At these nodes, part of the flow shall be delivered to satisfy cooling requirement, remaining flow (if any) shall be transported to the downstream node in the network. This is not applicable for Steiner nodes as these nodes are only pipe junction nodes with no cooling requirements, as expressed in Constraint (12). Constraint (13) represents flow capacity constraints. It enforces the logical relationship between the x- and the f-variables min , ϕ max ], other= 1, then fitm shall take any value in [ϕm (that is, if xm m ij j

= 0 ∀(i, j ) ∈ A, m ∈ M, t ∈ T ). wise fitm j The combination of Constraints (14–(17) ensure that supplied temperature of chilled water to all customer nodes and selected Steiner nodes is not exceeding the contractual/desired supply tem= 1 then we have perature. Indeed, we see that if zi j = 1 and xm ij m = 0 and we have t j = ti + Tim . On the other hand, if z = 0, then x i j j ij 0 = 0. Similarly, Constraints (18)–(21) guarantee that both pressure drop through pipes and across customer substation are not exceeded and within the allowable limits. Finally, Constraints (22)–(26) represent non-negativity, binary variables and integrality conditions.

It is possible to further enhance the model representation by appending the following valid inequalities. These constraints were initially proposed for the Steiner tree problem (see for e.g., Polzin and Daneshmand, 2001). They explicitly emphasize that the resulted graph is a directed, acyclic and connected graph.

z jk ≤



zi j ,

∀(i, j ) ∈ A, ∀ j ∈ S, ∀k ∈ V j+ ,

(27)

i∈V j−

zi j ≤



z jk ,

∀ j ∈ S.

⎡ ⎣

⎤



zi j = 1⎦ × t j ,

i∈V j−

⎡ ⎣



∀ j ∈ C,

(28)

(29)

k∈V j+

Constraint (27) is usually called the two-node subtour elimination constraint. It indicates that the resulting subgraph cannot include two arcs of opposite directions. Constraints (28) and (29) assert special connectivity rules on Steiner nodes due to their facultative nature. Both constraints are in the form of conditional constraints. Constraint (28) indicates that no outgoing arc from any Steiner node shall be included in the network if there is no incoming arc to this node. Similarly, Constraint (29) ensures that at least one outgoing arc shall be included if there is one incoming arc. 3.4. Model linearization As observed, the model is composed of equations (1)–(29), and contains a couple of nonlinear inequalities (namely, (14) and (18)). For that reason, we utilize the Reformulation-Linearization Technique (RLT) proposed by Sherali and Adams (1990) and Sherali, Adams, and Driscoll (1994) to derive an equivalent linear representation of the

(30)

⎤ zi j = 1⎦ × p j ,

∀ j ∈ C,

(31)

i∈V j−

2. Similarly, using (9), we construct the following inequalities:

⎡ ⎣

⎤



zi j ≤ 1⎦ × t j ,

i∈V j−

⎡ ⎣



∀ j ∈ S,

(32)

⎤ zi j ≤ 1⎦ × p j ,

∀ j ∈ S,

(33)

i∈V j−

3. In light of both (15) and (16), the following valid inequality can be derived.

⎤



⎣tmin

zi j ≤ t j ≤ tmax

i∈V j−



zi j ⎦ × zi j

∀ j ∈ S,

(34)

i∈V j−

4. Similarly, using Constraints (19) and (20), the following inequality is valid.

⎡

⎣ pmin

⎤



zi j ≤ p j ≤ pmax

i∈V j−

i∈V j−



1. Using (8), we construct the following equalities:

⎡

3.3. Valid inequalities

zi j + z ji ≤ 1,

model. The RLT is composed of two phases (i) reformulation phase and, (ii) linearization phase. In the reformulation phase, certain types of implied polynomial constraints, which include the aforementioned non-linear terms, are added to the formulation. This addition enables the utilization of new decision variables that replace the product terms in the original problem, as introduced during the linearization phase. A summary of the followed RLT procedures is presented in what follows. Reformulation phase: We append certain additional implied polynomial constraints to the problem. This is accomplished using the following steps:



zi j ⎦ × zi j

∀ j ∈ S,

(35)

i∈V j−

5. From Constraints (17) and (21); respectively, we construct the following valid equalities.

[tr = tmin ] × zr j ,

(36)

[pr = pmax ] × zr j ,

(37)

6. From (25), we construct the following valid inequalities:

[1 − zi j − z ji ≤ 0] × t j , [1 − zi j − z ji ≤ 0] × p j ,

∀(i, j ) ∈ A, ∀(i, j ) ∈ A,

(38) (39)

7. The product of Constraint (27) by the bound-factors in (15) and (19) yields to:

[1 − zi j − z ji ≤ 0] × (tmax − t j ), [1 − zi j − z ji ≤ 0] × ( pmax − p j ),

∀(i, j ) ∈ A, ∀(i, j ) ∈ A,

(40) (41)

8. Using (28), the following inequalities are valid:

⎡

⎣z jk ≤

⎤

 i∈V j−

zi j ⎦ × t j ,

∀ j ∈ S, ∀k ∈ V j+ ,

(42)

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

⎡

⎤

⎣z jk ≤



zi j ⎦ × p j , ∀ j ∈ S, ∀k ∈ V j+ ,

wi j + w ji ≥ p j − pmax + pmax zi j + (43)

9. From Constraint (29), we construct the following valid inequalities.

⎣

⎤





zi j ≤

i∈V j−

⎣

∀ j ∈ S,

(44)

k∈V j+

⎡ 

zi j ≤

i∈V j−

z jk ⎦ × p j ,

∀ j ∈ S,

(45)

k∈V j+

∀ j ∈ V, ∀(i, j ) ∈ A, ∀ j ∈ V, ∀(i, j ) ∈ A,



Timj xm ∀(i, j ) ∈ A, ij



Pimj xm ∀(i, j ) ∈ A, ij

(48)

(49)

After this, we replace each product term in (30)–(45) by its corresponding single variable from (46)–(49). Hence, forming below linearized constraints.



ui j +

Timj xm ∀ j ∈ C, ij = tj

i∈V j−

i∈V j− m∈M





wi j −

i∈V j−



i∈V j−

i∈V j−



i∈V j−

wi j −

i∈V j−

(50)

Pimj xm ∀ j ∈ C, ij = pj

(51)

Timj xm ∀ j ∈ S, ij ≤ tj

(52)

Pimj xm ∀ j ∈ S, ij ≤ pj

(53)

m∈M



ui j +

m∈M



ui j + u ji ≤ t j −

Timj xm ∀(i, j ) ∈ A, ij 

ui j + u ji ≥ t j − tmax + tmax zi j −

(54)

Timj xm ∀(i, j ) ∈ A, i j + tmax z ji

m∈M

(55)

u jk ≤

i∈V j−

 i∈V j−

ui j +

ui j +

 i∈V j−



Timj xm ∀ j ∈ S, ∀k ∈ V j+ , ij

(56)

m∈M

Timj xm ij ≤

i∈V j− m∈M



u jk

∀ j ∈ S,

(57)

k∈V j+

tmin zi j ≤ ui j ≤ tmax zi j wi j + w ji ≤ p j +



 m∈M

∀(i, j ) ∈ A, Pimj xm ∀(i, j ) ∈ A, ij

(61)

m∈M

Pimj xm ij ≤

i∈V j− m∈M

wr j = pmax zr j

∀ j ∈ Vr+ , ∀ j ∈ Vr+ ,



w jk j ∈ S,

(62)

k∈V j+

∀(i, j ) ∈ A,

(63) (64) (65)

Fundamentally, both system’s design and operations are considered in globally optimizing a DCS. In such event, the sizing of the plant together with the distribution network and its configuration are considered for optimization at once. However, and given the description in Section 2, one can notice that the problem can be decomposed into two related subproblems without compromising the system global optimality. Actually, the original problem includes two families of variables (as presented earlier in Table 3) •

•

(58) (59)

Family 1 includes plant design- and operation-related variables: y, g, F and I. Family 2 includes network design-related variables: z, x, f, t and p.

We observe that variables of Family 1 are used in Constraints (2)– (7), and (22)–(23) only. On the other hand, variables of Family 2 are used in the remaining constraints only. Hence, there are no coupling constraints for both families of variables (that is, there is no constraint that involves both families of variables). Consequently, the model can be decomposed into two independent sub-problems. These subproblems are:

•

m∈M



i∈V j−

pmin zi j ≤ wi j ≤ pmax zi j

•

i∈V j− m∈M



Pimj xm ∀ j ∈ S, k ∈ V j+ , ij



3.5. Model decomposition

m∈M



wi j −

(60)

(47)

m∈M

p j zi j = wi j −

wi j −

∀(i, j ) ∈ A,

Finally, the resulting model is an MILP model and is composed of (1)–(13), (15)–(17), (19)–(24), and (50)–(65).

Furthermore, we derive the following identities from (18) and (22) and p j xtm ; respectively. for linearizing t j xtm ij ij

t j zi j = ui j +

Pimj xm ij

(46)

and

w i j = pi z i j



ur j = tmin zr j

Linearization phase: We introduce two new RLT variables to linearize Constraints (14), (18), and (30)–(45). Toward this end, we define

ui j = ti zi j

 i∈V j−

⎤ 

w jk ≤

i∈V j−

z jk ⎦ × t j ,

 m∈M

+ pmax z ji

i∈V j−

⎡

655

Subproblem 1 addresses decisions related to the design and the operations of a chiller plant and a storage tank. Subproblem 2 addresses decisions related to the design of a main distribution network.

Even though these two subproblems are closely related, in the sense that both designs are essentially dependent on cooling demand and other technical parameters, they are considered to be independent given the described problem setting. This independence is explained by the following reasons. 1. The network design is only influenced by the quantity to be delivered to each customer at every period of time, and not impacted by the fact whether the distributed chilled water was just produced or was stored. 2. The sizing of chiller plant and storage tank depends on the total quantities demanded by customers and not on the network topology or its pipes’ sizes. As a matter of course, there are two determinants for sizing a chiller plant: (1) customers’ cooling load/demands (chilled water flow); and (2) system δ T. The latter is an important design consideration in a DCS, and any deviation from its designed value could result in significant efficiency and operational problems. In principle, the value of system δ T is a function of both supply and return temperature from and to the

656

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

plant (δ T = Treturn − Tsupply ); therefore, the return temperature is an equally important factor when optimizing systems capacity. Actually, low return temperature could adversely affect the available plant capacity. In our model, we catered for this by including δ T as a design parameter (input data) and not as a variable to be optimized. Given that the value of system δ T is controlled at customers’ ETS and not at the plant (through installing and operating proper heat exchangers, coils, valves etc.), imposing δ Tmin would necessarily control the value of return temperature to the plant; hence allows maintaining a desirable system performance level. This suggests that regardless how the solution network will look like; it has to adhere to temperature limits that are set in relation to a desired performance level.

Table 4 Characteristics of the problem classes. Problem

|V|

|C|

|S|

|T|

Sub-class

Peak cooling load (TR)

Class 1

15

10

5

4

20

10

6

40

20

8

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

10,000 20,000 40,000 80,000 100,000 10,000 20,000 40,000 80,000 100,000 10,0 00 20,000 40,000 80,000 100,000

Class 2

30

Class 3

60

In the sequel, both problems could be addressed and solved independently following the below MILP formulations. •

•

A model for optimal Plant Design and Operation (PDO): This model aims at optimizing Family 1 variables, by minimizing     plant storage prod yk + h∈H FCh gh + t∈T VCt Ft + t∈T VCtSto It , k∈K FCk subject to Constraints (2)–(7), and (22)–(23). A model for optimal Network Design (ND): This model aims at optimizing Family 2 variables along with the   , two additional RLT variables, by minimizing (i, j )∈A m∈M cimj xm ij subject to Constraints (8)–(13), (15)–(17), (19)–(21), (24)–(26), and (50)–(65).

4. Computational study In this section, we present the results of a computational study that aims at demonstrating the practical usefulness of the proposed models. The proposed MILP models were coded, and then implemented and tested on instances of various sizes with the aid of an optimization software package; namely, ILOG CPLEX 12.5. These numerical tests were guided by two objectives: • •

To test and validate the solvability of the proposed models. To examine their performance when generating optimal solutions and compare it against generating near optimal solutions by setting an optimality tolerance.

Given that our partner DC company refused to make real data public, hypothetical data were generated and then used to test the validity of the model while incorporating the full sense of a realistic case. More specifically, three classes of the problem were generated, in which each corresponds to a different network size. In our case, the size reflects the number of nodes that makes up a network; thus, generated problems cover small, medium and large scale networks. For each problem class, five instances were generated by varying the total cooling demands required during the peak hour/period. Hence, in so doing we generated a total of 15 instances. Table 4 below illustrates the key characteristics attributed to every problem class and subclass. Table 5 below gives an overview of the number of variables and constraints in models (PDO) and (ND) for each of the problem classes. 4.1. Performance of solving MILP to optimality The results that we have obtained after solving the developed instances optimally using IBM ILOG CPLEX 12.5 on an Intel(R) Core(TM) i7-3520M CPU @ 2.90 gigahertz with 4.00 gigabytes RAM are presented in what follows. The CPU time required to solve each of these instances is given in Table 6. Additionally, results are reported for each class as shown in Table 7. These reported results include mean, minimum and maximum CPU time required for solving each of the proposed models (computed over all solved instances in each class).

Table 5 Size of the problem instances. Problem

Class 1 Class 2 Class 3

Number of variables

Number of constraints

PDO

ND

PDO

ND

29 33 37

6755 22067 64092

15 21 27

6646 20744 47722

Table 6 CPU time for solving POD and ND models to optimality for every instance. Instance

PDO model CPU (seconds)

ND model CPU (minutes)

1.1 1.2 1.3 1.4 1.5 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 Average

1.74 0.87 1.12 0.92 0.99 1.1 1.09 0.92 0.98 1.12 1.03 1.06 0.96 0.93 1.07 1.06

0.14 0.29 0.26 0.16 0.15 6.83 21.87 16.90 23.24 15.18 14.49 42.09 142.95 378.88 750.31 94.25

As noticed, an optimal solution for the proposed PDO model can be generated within a very short time (within seconds). In addition to this, the proposed ND model can be solved optimally within 1 hour and 30 minutes on average. Moreover, it can be noticed that the largest instance, with 60 nodes, requires on average 4 hours and 25 minutes to be solved optimally. 4.2. Performance of solving ND model while setting an optimality tolerance In this section, we evaluate the performance of the heuristic strategy that amounts to solving the ND model while setting an optimality tolerance. In our experiments, we set the optimality tolerance to 1 percent, 2 percent and 5 percent. The results for each of these limits are displayed in Tables 8, 9 and 10, respectively. We see from this table that, if a 5 percent-optimality tolerance is set, then significant reduction of the CPU time is achieved (by 40 percent on average), while the proposed solution exhibit a reduced optimality gap that is often less than 1.5 percent.

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

657

Table 7 CPU time for solving POD and ND models to optimality for every class Class

|V|

Class 1 Class 2 Class 3

15 30 60

|C|

10 20 40

|S|

|T|

5 10 20

4 6 8

PDO model CPU (seconds)

ND model CPU (minutes)

Average

Minimum

Maximum

Average

Minimum

Maximum

1.13 1.04 1.01

0.87 0.92 0.93

1.74 1.12 1.07

0.20 16.80 265.75

0.14 6.83 14.49

0.29 23.24 750.31

Table 8 Performance of solving ND model by setting optimality tolerance of 1 percent.

Table 10 Performance of solving ND model by setting optimality tolerance of 5 percent.

Instances

GAP (in percent)

CPU time (seconds)

Percentage of reduction in CPU time

Instances

GAP (in percent)

CPU time (seconds)

Percentage of reduction in CPU time

1.1 1.2 1.3 1.4 1.5 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 Average

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.54 0.01 0.00 0.00 0.05

7.37 16.95 14.36 8.73 9.53 404.50 1275.46 972.21 1330.96 906.28 702.69 2017.39 4744.49 16,308.30 266,47.16 3691.09

14 3 6 7 10 1 3 4 5 0 19 20 45 28 41 12

1.1 1.2 1.3 1.4 1.5 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 Average

0.00 0.00 0.29 1.56 2.23 2.28 2.38 0.29 0.63 0.00 2.30 0.96 1.41 3.69 2.95 1.40

6.24 14.44 11.08 5.81 5.90 218.43 888.14 505.22 703.63 491.65 678.28 1171.82 4084.10 9419.89 26,085.58 2952.68

27 17 28 38 32 47 32 50 50 46 22 54 52 59 42 40

Table 9 Performance of solving ND model by setting optimality tolerance of 2 percent. Instances

GAP (in percent)

CPU time (seconds)

Percentage of reduction in CPU time

1.1 1.2 1.3 1.4 1.5 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 Average

0.00 0.00 0.29 0.00 1.14 0.57 1.60 0.29 0.00 0.00 0.20 0.96 0.01 1.14 0.82 0.47

7.49 15.82 13.53 8.33 8.36 387.40 1295.90 849.23 1195.97 554.61 667.59 1170.82 4689.20 15,355.00 25,620.90 3456.01

12 9 12 11 4 5 1 16 14 39 23 54 45 32 43 21

5. Conclusions and future research The inherent benefits of a DCS when adopted under right conditions and settings cannot be overstated. By being the least localized cooling option with the ability to bring in cost and environmental advantages, a DCS stands as a more attractive alternative when compared with the other conventional cooling schemes. To capitalize on such benefits, a proper planning is undoubtedly required. The significant initial investment cost required to build a DCS prompts having sound and more ideally, optimal, decisions to reap as much benefits as possible. In view of this, this research investigates MIP models that aid in making decisions related to chiller plant capacity, storage tank capacity, distribution network size and configuration, and quantities to be produced and stored during every period of time. Unprecedentedly, the proposed models were built to capture and reflect constraints related to the unique practical characteristics of DCSs

presented in both structural and technical sides. More specifically and unlike the previous efforts devoted in this area, models were built to include both thermal and hydraulic characteristics of the system. This was done by incorporating temperature- and pressure-related constraints. Given system’s complexity, the RLT was used to convert the resulted nonlinear formulation of the network design problem into an equivalent linear one. Computational experiments, carried out on networks with up to 60 nodes, demonstrated that both optimal and near-optimal solutions can be provided within a reasonable computing time. The presented research contributes to further strengthening the optimization efforts exhibited in the field of DCS design. More essentially, it is hoped that this work will not only help engineers to better design minimum cost DC networks, in a more efficient, effective, and operationally-feasible manner, but also further promotes DCS adoption in different part of the world. The optimization of a DCS is still in its infancy, and several research directions deserve to be investigated. Prospect works can build upon the current models by pursuing improvements or extensions, both scope and/or performance wise. This may include addressing trade-off between costs committed and Greenhouse Gas (GHG) emission levels, and studying the possible cost/GHG reductions benefiting both customers and cooling utilities. Additionally, a challenging generalization of the model that was investigated in this paper is the multiple-chiller plant variant. In this latter problem, we seek to find the optimal mix of cooling technologies, location and sizing of several chiller plants and storage tanks, and the allocation of customers to plants/tanks during every period of time in such a way that minimizes the overall costs. This multiple-chiller locationallocation problem offers the additional benefits of flexibility and reliability through tapping into the underutilized potential of renewable energy, reducing primary energy consumption; thus, offering lower GHG emission levels. We expect that the solution of this largescale complex network design problem will require the development of innovative exact and heuristic approaches. This is the topic of our on-going research.

658

R. Khir, M. Haouari / European Journal of Operational Research 247 (2015) 648–658

Acknowledgment The authors thank two anonymous referees and the Editor for their insightful comments that have helped improve the presentation of this paper. References Arab Construction World (2012). District cooling: A worthy opponent. ACW, 8, 10–12. Aringhieri, R. (2003). Optimal operations management and network planning of a district heating system with a combined heat and power plant. Annals of Operations Research, 120, 173–199. American Society of Heating, Refrigerating and Air-Conditioning Engineers (2009). ASHRAE handbook: Fundamentals. Atlanta, GA: ASHRAE. Booz and Company, (2012). Unlocking the potential of district cooling: The need for GCC governments to take actions. Technical Report. Booz and Co. Chan, A. L., Hanby, V. I., & Chow, T. T. (2007). Optimization of distribution piping network in district cooling system using genetic algorithm with local search. Energy Conversion and Management, 48, 2622–2629. Chow, T. T., Au, W. H., Yau, R., Cheng, V., Chan, A., & Fong, K. F. (2004). Applying districtcooling technology in hong kong. Applied Energy, 79, 275–289. DHC/+ Technology Platform, (2009). District heating and cooling: A vision towards 2020-2030-2050. www.dhcplus.eu/wp-content/uploads/2012/05/120529_ Vision_DHC_final.pdf. Euro Heat and Power, (2013). District heating and cooling- 2013 statistics. Available at http://www.euroheat.org/Statistics-69.aspx. Accessed 27.03.14. Feng, X., & Long, W. (2008). Applying single parent genetic algorithm to optimize piping network layout of district cooling system. In Proceedings of the Fourth IEEE international conference on natural computation (pp. 176–180). Feng, X., & Long, W. (2010). Optimal design of pipe network of district cooling system based on genetic algorithm. In Proceedings of the Sixth IEEE international conference on natural computation (pp. 2415–2418).

Hart, D. R., & Rosen, M. A. (1996). Environment and health benefits of district cooling using utility-based cogeneration on ontario, canada. Energy, 21, 1135–1146. Henning, D., Amiri, S., & Holmgren, k. (2006). Modelling and optimisation of electricity, steam and district heating production for a local swedish utility. European Journal of Operational Research, 175, 1224–1247. Ming, D., Fenglan, H., & Hairui, W. (2012). Energy supply network design optimization for distributed energy systems. Computers & Industrial Engineering, 63, 546–552. Pinson, P., Nielsen, T. S., Nielsen, H. A., Poulsen, N. K., & Madsen, H. (2009). Temperature prediction at critical points in district heating systems. European Journal of Operational Research, 194, 163–176. Polzin, T., & Daneshmand, S. V. (2011). A comparison of steiner tree relaxations. Discrete Applied Mathematics, 112, 241–261. Powell, K. M., Cole, W. J., Ekarika, U. F., & Edgar, F. T. (2013). Optimal chiller loading in a district cooling system with thermal storage. Energy, 50, 445–453. Rastpour, A., & Esfahani, M. S. (2010). Mathematical models for selection of optimal place and size of connections considering the time-value of money. European Journal of Operational Research, 200, 764–773. Sakawa, M., Kato, K., & Ushiro, S. (2001). Operation planning of district heating and cooling plants through genetic algorithms for nonlinear 0–1 programming. Computers & Mathematics with Applications, 42, 1365–1378. Sakawa, M., Kato, K., & Ushiro, S. (2002). Operational planning of district heating and cooling plants through genetic algorithms for mixed 0–1 linear programming. European Journal of Operational Research, 137, 677–687. Sherali, H. D., & Adams, W. P. (1990). A hierarchy of relaxations and convex hull characteristics for mixed-integer zero-one programming problems. Discrete Applied Mathematics, 52, 83–106. Sherali, H. D., Adams, W. P., & Driscoll, P. J. (1994). Exploiting special structures in constructing a hierarchy of relaxations for 0–1 mixed integer problems. Discrete Applied Mathematics, 46, 396–405. Soderman, J. (2007). Optimisation of structure and operation of district cooling networks in urban regions. Applied Thermal Engineering, 27, 2665–2676. Uhlemair, H., Karschin, I., & Geldermann, J. (2014). Optimizing the production and distribution system of bioenergy villages. International Journal of Production Economics, 147, 62–72.