Cooling in Singapore

Apart from the aforementioned existing district cooling system [1], a recent study on Singapore [7] focused on comparing surface and subsurface cooling in tropical countries using an EnergyPlus model of a representative building. It concludes that an open-loop groundwater cooling system might be a better option for tropical cities like Singapore compared to conventional systems. The viability of district cooling in the CBD area of Singapore has already been assessed in a working paper [8] by Singapore Power Ltd.

District cooling

A thermal ice storage system and its influence on the cost-optimal design of district cooling distribution network is investigated [9] using non-linear optimization in two exemplary networks. It finds that ice storage can help reduce the required pipe diameters. A TRNSYS model is built to investigate a similar question [10] for a case study in Hong Kong; storage is not found to be an economically attractive option there.

A genetic algorithm is used [11] to assess the optimum shares of five building types to form a cooling load time series that is most suitable for district cooling. This approach was later embedded in an optimization scheme for a development project in South East Kowloon, Hong Kong. Another genetic algorithm was used by Feng and Long [12] to determine and compare piping layouts.

Khir [13] presents a mixed integer linear programming (MILP) model for designing district cooling system. It explicitly models flow rates, temperatures and pressures in a graph of arcs and nodes. Linearization of the physical constraints leads to a large amount of auxiliary variables. Both plant design and operation and network design can be optimized. However, to achieve a tractable problem size, those two subproblems are solved separately. However, this separation requires fixing the supply task. The paper also highlights the relationship with Steiner tree problems.

Södermann [14] presents a similar MILP model for designing district cooling systems in urban areas. It uses a linearized cost function for sizing equipment (plants, pipes, storages). A small set of representative load periods (8 in the example study) approximate the annual cooling load duration curve.

District heating

Due to the basic principle, district cooling and heating are comparable types of supply. In both types, with few central generation plants heat or cold is generated which is mostly transferred to water. This tempered water is transported through pipes to customers buildings. Having transferred the heat or cold to the buildings, the cooled or heated water returns to the generation plants where the process starts from the beginning. The biggest difference between both are the technologies used in generation plants, which is not the focus of this article.

Jamsek [15] presents a linear programming (LP) approach for determining cost-minimal pipe topology and capacity configuration for a given structure and cost data. A more conceptual, but worthwhile discussion of district heating and cooling systems is given by Rezaie [16].

Udomsri [17] presents and compares the different options to combine centralized and decentralized heating and cooling. It points towards solutions that use decentralized, thermally driven chillers to produce cooling on-site. This configuration can have advantages in regions with both heating and cooling demands throughout the year.

Multi-commodity district energy system

A detailed model on combined energy system planning for electric, heating and cooling demand is presented by Chinese from 2008 [18]. It presents a MILP model. A case study in Udine, Italy, is performed with 7 nodes and 41 individually weighted time steps to represent the seasonal demand characteristics of a full year.

A recent paper [19] investigates the case of combined district energy system for electricity, heating and cooling using a mixed integer linear optimization problem. Both spatial (7 nodes) and temporal (2 days in 48 time steps) resolution are low in the presented case study, but the non-linear performance of gas turbines is represented by a set of linearized constraints.

A broader review of computer models for renewable energy system analysis and optimization is provided by Connolly et al. [20]. Weber presents in her thesis [21] a comprehensive methodology to plan cost-optimal polygeneration energy systems.

Sets

A district or city is represented as a graph of vertices and arcs. This graph is derived from the street network. It should be derived in such a way to include all considerable locations for network pipes. The spatial resolution can be chosen as fine as required. Here, it is on the level of building blocks, i.e. a single street segment between two intersections.

Let V be the set of vertices v_i, corresponding to connecting or terminating points of the graph. Set A of arcs then comprises ordered tuples of vertices a_ij=(v_i,v_j) with i≠j. A is symmetric, that means either both or none of the pair a_ij and a_ji are elements of A. For readability, the subscript ∘_i is used to denote any parameter or variable that is defined over vertices v_i∈V, while ∘_ij is used to denote a quantity defined over arcs a_ij∈A.

The set V₀⊆V defines so-called source vertices, which represent locations of possible cooling sources. Usually, the number of source vertices is low (≤10) compared to the size of the graph.

The neighbor sets N_i of a vertex v_i are defined by the indices of all vertices that are connected to it by an arc: N_i={ k∣a_ki∈A }.

Time is represented by a set T of discrete time steps t, which represent a small number of representative operational situations. The minimum viable number of time steps is two: one for peak demand with a duration of one to several hours, a second for the average annual load with a duration of the remaining year. A more comprehensive choice are three to four time steps: The third can refer to the all-year base load, while the fourth can be used to represent a common intermediate load level. Additional time steps need to be introduced if redundancy requirements are to be stated. Refer to the discussion of the parameter availability below for more details.

Parameters

The numerical given facts for this model are shown in Table 2. They are grouped by their defining domains. First are vertex parameters, then arc parameters, then constant or global parameters, and finally time-dependent parameters.

Vertices have a single parameter: their maximum power capacity IEq6_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[Q^{\text {max}}_{i}\] \end{document}. It is the thermal output power given in kW for that location. All non-source vertices have this parameter set to 0.

Arcs have three main attributes: their length l_ij (m), the thermal peak demand d_ij (kW) of their adjacent buildings and g_ij (binary), which states whether a pipe already exists in this arc.

As secondary parameter, IEq7_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[C^{\text {max}}_{ij}\] \end{document} indicates the maximum thermal power capacity (kW) in an arc, which can be derived from the maximum available pipe diameter. For arcs with existing pipes (g_ij=1), this value should be set according to the diameter of that existing pipe. As both arcs a_ij and a_ji refer to the same street segment, parameter values for both members are always identical.

Global parameters are technical and economic parameters that refer to the district cooling system as a whole.

Economic parameters are all costs and revenues. The investment costs are split into a fixed part c_fix and a variable part c_var. The fixed part, given in S$/m, captures all costs that are not dependent on the capacity of the pipe to be built, mainly earth works. The variable part, given in S$/(kW m), refers to the capacity-dependent cost component of building a pipe. Both values must be tailored to the study area to compensate for differences in cost structure and availability of different pipe sizes. The same is true for the operation & maintenance (OM) cost parameter c_om in S$/m. In contrast to the fixed investment cost term, it is also to be paid for existing pipes. The cost for providing cooling IEq8_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[c^{\text {cool}}_{v}\] \end{document} here is dependent on the vertex to allow for representing cheap surface water cooling in contrast to more expensive re-cooling options.

The letter c and a subscript denote economic parameters. These are investment costs for building the pipe network c_fix and c_var, maintenance c_om, costs of providing IEq9_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[c^{\text {cool}}_{i}\] \end{document} cooling, and revenue for delivering c_rev cooling to consumers. While c_fix contains the size-independent costs (mainly earth works), c_var contains costs that are dependent on the thermal capacity (diameter) of the pipe.

Time-dependent parameters s_t and w_t represent scaling factor (dimensionless) and weight or duration (h) of a time step t. A value of s_t=1 refers to peak demand, while smaller values correspond to moments of partial load. Together, these two parameters encode a discretized annual load curve.

Redundancy requirements can be stated in the model by setting the binary availability parameter y_it to value 0 in certain time steps for 1,2,… source vertices. If one time step for each foreseen failure configuration is introduced, a full n−1,n−2,… failure safety (against cooling source outages) can be enforced in the produced solution. Clarification: in the conventional time steps discussed above in paragraph time step set, the value y_it=1 is to be set for all source vertices v_i∈V₀.

Variables

The main optimization task is to find values for the binary decision variable ξ_ij. If its value is one, a pipe is built in the street segment corresponding to the arc a_ij. For each time step, the actual use of a given pipe is decided by setting the binary pipe usage decision variable ψ_ijt. If its value is one, the pipe in arc a_ij is used in direction from vertex v_i to vertex v_j. Consequently, there must be a power flow IEq10_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\pi ^{\text {in}}_{ijt}\] \end{document} into the pipe. Simultaneously, a value ξ_ij=1 requires that the demand d_ij of this arc has to be satisfied at all times. The power flow variable at the other end of the pipe is called IEq11_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\pi ^{\text {out}}_{ijt}\] \end{document}. Its value is determined by the ingoing power flow minus losses and demand.

The non-negative variable ρ_it represents thermal power output from a source vertex v_i∈V₀ at time t in kW. It is limited in value by the source vertex capacity parameter IEq12_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[Q_{i}^{\text {max}}\] \end{document} and possibly by the availability parameter y_it. Table 3 summarizes all model variables. The following line shows their domains, on which they are defined:

Equations

Equations fall into two categories: The first is the cost function whose value is to be minimized. The second are constraints that codify all the previously discussed rules in mathematical form. By that, they define the region of feasible solutions, under which the solver selects a (close to) cost optimal solution.

The cost function value ζ is the sum of costs for providing cooling, minus revenue for that cooling. Network costs occur for building the network (annualized investment) and for operation & maintenance. The following derived parameters are introduced to shorten the definition of the cost function whose definition is given below: IEq16_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[k^{\text {fix}}_{ij}\] \end{document} (S$/a) and IEq17_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[k^{\text {var}}_{ij}\] \end{document} (S$/(kW a)) cover network costs (investment, O&M). IEq18_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[k^{\text {cool}}_{ij}\] \end{document} (S$/kWh) represents cooling generation costs, while IEq19_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[r^{\text {cool}}_{ij}\] \end{document} (S$/h) represents revenue for provided cooling. With these, the cost function is equal to

The first summand forms a piece-wise linear function that is depicted in Fig. 2. The plot also depicts that the linear approximation possibly over- and underestimates the cost of either small or large capacities. Depending on the parameterization, one can either force an optimistic (by deliberate underestimation) or cautious (by deliberate overestimation) planning decision. The linear approximation exhibits unsteady jump from zero to the fixed cost uc_fix+c_om for any positive value of IEq20_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\bar {\pi }_{ij}\] \end{document} for any non-existing (g_ij=0) pipe. The link between the binary decision variable ξ_ij and the continuous pipe capacity variable IEq21_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\bar {\pi }_{ij}\] \end{document} is enforced through Eqs. (8), (12), and (14) below.

These are the definitions for the four derived parameters used in the cost equation above. The division by value 2, common to all three arc parameters, compensates for the double-representation of one street segment {v_i,v_j} as a pair two directed arcs a_ij, a_ji. More about this issue can be found at the explanation of the symmetry constraints (15) and (16) below. The cooling cost parameter IEq23_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[k_{i}^{\text {cool}}\] \end{document} is divided by the concurrence effect parameter b to remove the load reduction effect that is introduced in Eq. 10. In other words, the power flow that needs to be satisfied from a source vertex is lower by b than the energy flow that needs to be delivered. The division by b removes that difference in terms of costs.

The constraints formalize all physical laws and technical rules that need to be satisfied by a feasible solution. The first constraint concerns energy conservation in vertices v_i, with respect to its neighbors N_i:

This inequality is thus a relaxed version of the law of energy conservation. Relaxed, because it allows for power to vanish at any vertex. As power does not come for free, any solution returned by the solver will satisfy this constraint to equality. For all non-source vertices (i.e. ρ_it=0), the difference between outgoing and incoming power must be smaller than or equal to zero. In source vertices, a positive difference may remain when it is met by an equal amount of input power ρ_it into the network at that vertex.

Demand satisfaction is the main arc constraint. It creates the logical connection between the discrete pipe usage decision ψ_ijt and the demand satisfaction. Like the cost equation, it relies on two derived parameters x_ij and y_ijt. This constraint is also graphically explained in Fig. 3.

Parameter x_ij refers to losses w^var that are proportional to the amount of power flow into the arc a_ij. Parameter y_ijt refers to fixed thermal losses w^fix that only depend on the pipe length, not on the power flow, and the time-dependent demand d_ij·s_t that is further reduced by the concurrence effect b and the connect quota parameter q. In notation:

Pipe capacity is a technical constraint that limits the power flow through a pipe by the built capacity. This is accomplished by first limiting the ingoing power flow into a pipe IEq26_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\pi _{ijt}^{\text {in}}\] \end{document} by the built pipe capacity IEq27_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\bar {\pi }_{ij}\] \end{document}. The second constraint is required to set the pipe usage decision variable ψ_ijt. Its value is forced to 1 if a flow in direction from v_i to v_j at time step t is required. The reason for this seemingly artificial constraint is explained in the following constraint.

Unidirectionality of the power flow IEq28_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\pi _{ijt}^{\text {in}}\] \end{document} is required, as otherwise the solver would happily use a pipe’s capacity in both directions simultaneously, which would not be physically possible. Therefore, the pipe usage decision variable ψ_ijt may only have the value 1 in one direction at any given time step t. This way, the direction of flow may change between time steps, but any given pipe may only be used in one direction at a given time.

Build capacity limits the pipe capacity variables to the maximum available diameter, whose capacity in kW is represented by parameter IEq29_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[C_{ij}^{\text {max}}\] \end{document}. At the same time, this constraint sets the value of the building decision variable ξ_ij, if a non-zero pipe capacity IEq30_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\bar {\pi }_{ij}\] \end{document} is needed.

Symmetry of building decision are two constraints that enforce that the mathematically independent arcs a_ij and a_ji must have identical variable values for their pipe capacity IEq31_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\bar {\pi }_{ij}\] \end{document} and building decision ξ_ij. This way, a pipe can then be used in both directions at different time steps.

Use if built is the last piece in the puzzle to link the building decision ξ_ij to the pipe usage decision ψ_ijt. Up to this point, nothing in the formulation requires that the cooling demand of a customer is satisfied at all times. This constraint enforces exactly that, by requiring that the sum of ψ_ijt+ψ_jit is greater or equal to 1, if a pipe is built along that arc. Otherwise, the condition should not be enforced. This is done by calculating the expression (ξ_ij+ξ_ji)/2. It is equal to 1 if ξ_ij=1, as Eq. (15) requires ξ_ij=ξ_ji. Otherwise, its value is equal to 0. In other words: if an arc is has a pipe, its demand must always be satisfied by a power flow from any of its two sides.

Source vertices is the last constraint. It limits the source vertex power flow ρ_it by its maximum allowed capacity IEq32_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[Q_{i}^{\text {max}}\] \end{document}. The availability parameter y_itusually has value 1, except for special time steps in which the source vertex v_i is made unavailable by setting y_it=0:

Equations 2 to (18) together define a linear mixed-integer program, whose optimal solution is the desired vector of pipe capacities IEq33_TeX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} \[\boldsymbol {\bar {\pi }}^{\star }\] \end{document}, accompanied by optimal cooling power flows.

Base case

The base case is designed to show today’s demand situation. In this scenario, the size and layout of a profitable district cooling network is to be determined.

In the base case, a revenue of 0.14 S$/kWh is assumed. A sensitivity analysis with reduced revenues in steps of 0.02 S$/kWh is also conducted to determine how sensitive the optimal network size is on the price for cooling.

Growing demand

In order to assess how the total cooling costs (generation and distribution) are affected by a possible demand growth in the study area, additional load is introduced in the edge (54, 114) in the south-eastern corner. It’s value is changed from 0MW to 200MW with steps of 50 MW. As this step creates more load than the existing cooling stations could satisfy while satisfying the redundancy constraint, all cooling stations’ capacity is increased by 50%.

Base case

The resulting network for the base case is shown in Fig. 7. It shows a district cooling network that satisfies almost all (1240MW of 1241MW) cooling demand. In other words: under the presented cost assumptions and a cooling revenue of 0.14 S$/kWh, district cooling is an economically viable strategy. A sensitivity check with a reduced cooling revenue of 0.10 S$/kWh still leads to an almost full saturation of 1237 MW supplied demand. The unsatisfied demands are mainly short, isolated edges with low demand (≤3MW) that cannot offset the network investment cost. Only when revenue is further reduced, between 0.08 and 0.06 S$/kWh, the optimal network drops from satisfying 1200 MW to 0 MW. This value range thus can be taken as a lower bound for profitability of a district cooling system (under the given cost assumptions).

Figures 8 and 9 compare the cooling power flow in two of the 18 time steps considered for designing the network capacities in Fig. 7. The first situation shows the cost-optimal cooling power flow when all cooling stations are operational. The three cheap cooling stations (57, 110, 154) are delivering their full capacity of 150 MW each, while several of the other stations only run in partial load. The second situation in Fig. 9 shows how the redundancy requirement leads to the consideration of very different flow situations. Here, cooling station 13 now runs at full capacity and satisfies the middle of the study area.

Growing demand

The resulting cooling cost for gradually increasing the demand in edge (54, 114) for each case is shown in Fig. 10. The cooling cost here is calculated by dividing total costs for generation and network by the amount of provided cooling (in kWh). Cost for cooling grows slightly – along with demand – from approximately 0.045 S$/kWh to 0.048 S$/kWh in the case of additional 200MW peak cooling demand. This increase is caused by a higher utilization of the more expensive 0.07 S$/kWh cooling stations. The specific network costs actually exhibit a minor reduction as demand grows, but this effect is negligibly small here.