- Email: [email protected]
The algorithmic structure of a decision support system for a design of a district heating network
0305-0548190 $3.00+0.00 Copyright C 1990 PergamonPressplc
Compurers Opns Res. Vol. 17.No. 2,pp. 221-230,1990 rights reserved
Printedin Great Britain.All
THE ALGORITHMIC STRUCTURE OF A DECISION SYSTEM FOR A DESIGN OF A DISTRICT HEATING
SUPPORT NETWORK
A. SCIOMACHEN* and R. SOZZI? Department
of Information
Sciences, University of Milano, Via Cicognara 7, 20129 Milano, Italy
(Receioed June 1987;
revised
May
1989)
Scope and Purpose--Heating networks are gaining increasing importance as a cost effective and environmentally sensible way of distributing energy for the heating of residential and commercial buildings. The algorithms discussed in this paper, included in a decision support system (DSS), have been shown to be effective in a number of cases. The planning and the design of a network is a complex process: a DSS could reduce signilicantly the time and the cost involved.
Abstract-This paper is devoted to describe the algorithms implemented in a decision support system for the design of a district heating network. The design problem is formulated on a graph as the search for the minimum cost tree rooted at the source and spanning all the sinks. This search consists of two procedures: the first one finds, using greedy-like algorithms, an approximate solution which is subsequently improved in the second procedure, using local searches.
1. INTRODUCTION
The design of a district heating network is a complex process in which different factors like costs, reliability and the environmental impact must be taken into account. For this reason computer aided design techniques have been increasingly in use to assist the designer in assessing the trade-offs among all the factors involved and identifying the final design to be implemented. The first application of operations research techniques, at least to the authors’ knowledge, is in Refs [l] and [2] where the problem of deciding the operation of the production units and the dimensions of the pipe system is addressed, while the network topology is assumed to be given. The authors of this paper have been involved in a research activity, sponsored by the Energy Project of the Council for National Research, aimed at developing a decision support system (DSS) for the design of a district heating network. The first results obtained in this research are reported in Refs [3] and [4]. In this paper the cost is assumed as the performance index of the network: in this case the task of the designer is to link, at a minimum cost, a source to a number of demand nodes or sinks. This paper is devoted to illustrating the algorithmic procedures implemented in the DSS in order to compute the cost of a given network, to check whether it satisfies the constraints and to find the network of minimum cost. The main decision variables are the layout of the network and the capacity assignments to its arcs. The constraints are of two kinds. First, the features of the area to be served may restrict the number of possible layouts and second, the network is to be hydraulically feasible, i.e. the pressure must be within some range and this restricts the choice of pipe diameters. The main task of the DSS is to assist the designer in the formulation of the objective function and the constraints, to provide an optimization procedure and to assist the designer in controlling its execution and assessing the validity of the solutions obtained. The data typically available are the energy at the source and the requirements at the sinks, a description of the area *MS A. Sciomachen received the degree in information sciences from the University of Milano. She is Assistant Professor of Operation Research at the University of Milano, and has worked on space robotics at Tecnospazio S.P.A. and on manufacturing automation at IBM Yorktown Heights, New York. Her main interests are in combinatorial algorithms and modeling & simulation of manufacturing systems, particularly network design and Petri Net based models of manufacturing systems. tMr R. Sozzi received the degree in nuclear engineering from the Polytechnic Institute of Milano. He is working on applications of operations research to energy conservation and environmental assessment at CISE (Milano) and TIM (Roma). He was also Adjunt Professor of Operation Research at the University of Milano. His main interests are in combinatorial algorithms and modeling & optimization of environmental systems, particularly network design and optimal location of monitoring networks. 221
to be served by the network, i.e. where and at what cost the pipes can be placed, and the technical and economical data of the main network components (different diameters commercially available, pumping stations). The cost optimization problem, as described in Section 2, is formulated on a graph whose arcs represent the possible layouts and whose nodes are split inta a source, a set of sinks and a set of intermediate nodes. The search for the network of minimum cost, which is hydraulically feasible, is restricted to the trees of the graph rooted at the source and spanning all the sinks. This combinatorial problem is NP-hard [5]: actually, if there is only one available arc capacity, it is a particular case of the Steiner problem f6]. For this reason approximation algorithms are to be used to solve it. An interesting approximation algorithm to solve the Steiner problem has been recently suggested in [?I: the algorithm has a very good computational performance but cannot be applied to the proposed network design problem muse the cost of each arc in the heating network is expressed not only by a distance function but depends also on the arc capacity. The algorithm used in the DSS, described in Section 3.2, is based on another approach and consists of two procedures. The first one finds an approximation to the minimum cost network and the second improves the outcome of the first procedure using local searches. It is important to note that the input to the second procedure can also be provided by the designer on the basis of any given network design independently of the execution of the first procedure. Sections 3.1 and 3.3 are respectively devoted to the procedure for computing the cost of the network and to check the hydraulic constraint. The algorithmic structure of the DSS proposed in this paper has been implemented using the C language fg] on the IBM PC-RT using the AIX Operating System. The algo~thms and the data structures used have been carefully tailored to ensure a high level of interactivity and portability to a personal computing environment. This system has been used in the design of the district heating network in a medium sized town in northeastern Italy. The computational results reported in this paper are related to this application. One can observe that the DSS can assist in the design of other types of distribution networks. The model, introduced in the next section, is simple yet general to easily represent slurry pipelines [9}, gas or oil pipeline systems [lo], or drainage systems [l 11. 2. THE NETWORK
MODEL
The network design problem is modeled on a graph G = (I$ E) with the following speci~cations: tr= (sl tr, I) is the set of n -t- 1 nodes where s represents the source, U the set of sinks and I the set of intermediate nodes. At the intermediate nodes no energy is consumed or generated. By using these nodes, one can model some constraints of the area served by the network or the possibility for the pipes to join and fork at some points to obtain a more effective layout. In some cases the set I can be empty. To each node UE V is associated a thermal load r(o) such that r(u)=O,
DEI
r(u)>O,
YEU
r(s) = - 1 r(v), RY Le. it is assumed that the energy available at the source can meet all the requirements at the sinks. E c V x f/ is the set of m arcs connecting pairs of nodes of the graph. The set E is partitioned into two sets El and E2, where El represents the set of arcs which must be in the final layout of the heating network, and E2 is the set of the arcs among which a selection takes place. As an example of this problem, the design of a district heating network of a medium sized town in northeastern Italy is presented. A first analysis of the data yielded 21 sinks which should be connected to the source, for a global energy requirement of 11 x 10’ calorie/hr, that is the global flow of about 53 l/set, assuming that the flow temperature is 120°C out of the source and 60°C back to it. Table 1 lists the requirements of the 21 sinks. The possible network layouts are represented in Fig. 1. The resulting graph has 32 nodes and 40 arcs.
The algo~~c
structure of a DSS
223
Table I
Energyrequirements Sink
ww 6.54 1.62
1 L
3 4 5 6 7 8 9 10 11 12 13 14 I5 16 17 18 19 20 21 Total
0”::
1.77 1.84 1.84 0.74 12.98 0.20 1.00 0.44 0.59 0.44 0.47 0.53 Z:Z 1.92 15.20 2.61 S2.M
1
*
SOURCE
A
SINK
l
INTCRMEOIATE NOBE
100M -
Fig. 1
The topology of the network is assumed to be a tree G’ = (V’, E’), where E’ = El u E2’ and E2’ c E2, of the graph G with the root at the source and spanning all the sinks. This tree contains a subset f’ c 1. I’ can be 0 also when f # 0. Since trees are subgraphs without closed cycles, they are a reasonable model for the topology of a least cost ~st~bution network. To each arc (i, j)~ E is to be assigned a Sow fil, which satisfies the set cr of balance equations: 1. fil-0, @,j)BE 2. ~~~Li=~ffl=O,
ioI,k, joV
3. ~fku=~ful=M), u~U,k,jeV 4. cfsI=r(s) and fjs=O, VjeV.
(1)
Relations 2,3 and 4 express respectively the balance equations at the intermediate nodes, the sinks and the source. In order for the spanning tree to represent a network, i.e. to allow a feasible now, a diameter d, is to be assigned to each arc (i, j) E E’ which allows a flow fil satisfying the baIance equations
224
A. SCIOMACHEN and R. Sozzr
and satisfying a further constraint c2: max (Pi) < PN,
(2)
isV
where Pi is the pressure at node i, i.e. the pressure in the network is less than some nominal value P,,,.
3. THE
ALGORITHMIC
STRUCTURE
The algorithmic structure of the DSS can be summarized in the following steps, where G is the initial graph, G’ is the initial tree network and G* is the final solution. (1) Check whether G is connected. If G is not connected the problem has no feasible solution. (2) Check the existence and the feasibility of the initial solution G’. If G’ = 0 there is no initial tree network. If El # 0 and El has a cycle the problem has no feasible solution. Otherwise generate G’ by Procedure 3.2.1. If G’ # 0 and G’ has a cycle the problem has no feasible solution. (3) Computation of the cost of G’. Compute the cost of G’ according to Section 3.1 and 3.2.2. (4) Local search and optimization. Call Procedure 3.2.2. (5) Check the hydraulic feasibility. Compute the pressure along G’ according to Section 3.3. If G’ is hydraulically feasible G* = G’. Otherwise enter the diameters assignment policy chosen to compute the new diameters assignment according to Section 3.3. Note that steps 2-5 can be repeated more than one time with different starting solutions generated by Procedure 3.2.1 or provided by the user. 3.1. Computation of the cost The procedure for computing the cost is the following. For the available pipe diameters, the optimal flow and the corresponding cost of a unit length pipe diameter are computed. To this cost is added, for each arc of the graph, a possible additional unit cost related to the features of the area in which the arc is located. Pipelines are available in a finite selection of standard diameters. The unit investment cost of a pipeline inclusive of the installation is given for each available diameter. These costs must be amortized over some time period to convert total capital cost to an annual cost. Any amortization scheme can be used in the design method. The cost of the network is given by C=
C
C(dij)lijXij,
(LjkE’
where xii is a Boolean variable which is 1 if the arc (i, j) is chosen in the network and 0 otherwise, c(di,) is the cost of one unit length of pipe of diameter d,, and I, is the length of the acr (i, j)EE. The term c(d,,) includes (a) the investment cost, (b) operating costs (namely pumping costs, because the flow must be pumped along the network), and (c) the cost of heat losses proportional to the diameter of the pipe, the flow temperature, and the length of the period of operation. The value d, associated with each arc (i, j) is chosen in such a way that allows the flow fij at a minimum cost c(d,). For a given flow fij one computes the cost for each available diameter and chooses that diameter ensuring the minimum cost value under the constraint that the flow velocity along (i, j) does not exceed a value V,,,,,. This procedure performed for each f gives, for each diameter, a range Lfmin,f,,,.,] for which the value d,j is optimal. The proposed network design problem could be formulated as the determination of the minimum cost spanning tree subject to constraints (1) and (2). Table 2 reports the relevant data about commercially available diameters for the example considered in this paper. Column 1 is the diameter (cm); column 2 is the lower/upper bound of the optimal flow range (l/set) computed according to the above procedure (note that the lower bound is the same value as the upper bound of the closest smallest diameter), while column 3 gives the cost (S/m). Other parameters of the problem have been fixed. The technical life of the network is 20 yr, the capital amortization rate is 5 %, the electrical energy cost is 111 $/MWhr, the value of the heat losses is 23 %/IO’ calorie/hr, the efficiency of the pumps is 0.7 and the hours of operation per year are 2000 hr/year.
The algorithmic
structure
of a DSS
225
Table 2 Diameter (cm)
Min/max flow Wsec)
20 25 32 40 50 65 80 100 125 150 200 250 300
0.0, 0.3 0.3, 0.5 0.5, 0.7 0.7, 1.1 1.1, 1.8 1.8, 2.8 2.8, 4.5 4.5, 9.6 9.6, 15.3 15.3, 25.0 25.0, 56.1 56.1, 91.3 91.3, 150.6
Cost’ (S/m) 90 93 96 LOO 104 107 111 115 122 141 171 204 244
*The costs are computed on the basis of Italian currency. quoted as l$ = 1.350 Lira.
3.2. The network search procedure
The search for the minimum cost network is performed by two procedures. The first one, described in Section 3.2.1, is a greedy algorithm which provides a first approximation of the minimum cost network. This design is improved by the second procedure, described in Section 3.2.2, which is based on local search techniques. 3.2.1. The procedure for generating the starting solution. In the procedure used to generate the starting tree network, if the arcs are not provided by the designer, they can be chosen according to their lengths. In order to generate the initial solution for the local search technique, the algorithm performs the following two operations setting initially E’ = El. (1) Sort the arcs in E2 in nondecreasing order according to their lengths. (2) Add in turn each of the arcs of the sorted E2 to E’ if it does not cause a cycle in G’, until G’ is a tree that spans all the sinks. A better solution of the proposed problem can be obtained starting from (possibly) different initial solutions to be later improved by the second procedure. A natural way of obtaining different starting solutions is to replace the deterministic selection of each arc, performed in step 2, by a randomized selection on the arcs of E2 generated by using a discrete random variable; the corresponding probability mass function can be, for example, as shown in [3], a truncated geometric distribution over the integers 1,2, . . . , IE21, where (E21 is the number of arcs in E2. In the present case, the designer provided neither a starting tree network nor a selection of some arcs to be included in the final layout, i.e. in this case El = 0. The procedure summarized in Section 3.2.1 has therefore been executed with randomized arc selection, obtaining 3 different starting solutions. Figure 2 displays the 3 networks generated with three independent executions of the procedure. Network 1 has a length of 3885 m, a cost of $498,OOOand 31 nodes. The same values for network 2 and network 3 are respectively 3435 m, %405,OOO, 31 nodes and 4095 m, $51O,OOO and 32 nodes. 3.2.2. The local search procedure. For a complete analysis of the structure of local searches and their applications in combinatorial optimization the reader is referred to [12]. In this paper the authors outline the technique as it applies to the problem under consideration. This procedure uses an “elementary transformation” of a tree T, which consists in deleting one arc of T and adding to T one arc of G. Two trees Tl and T2 of G are said to be equivalent if one can be obtained from the other by this elementary transformation. Given a network G’ = (V’, E’), its neighborhood N(G’) is defined as the set of networks G” = (Y”, E”) such that G” is equivalent to G’. The network in N(G’) of minimum cost is termed locally optimal. In order to compute this locally optimal network, as in [lo], the neighborhood N(G’) is generated performing the following operations for each element of I/‘. (1) Let V” be the set of the nodes of Y’ not yet selected and set initially V” = V’. (2) If V” # 0 choose a node V,E I/” and let Y” = V’ - ui.
A. SCIOMACHEN
226
and
R. Sozv
Fig. 2 (3) Let A(ui) be the set of the b nonadjacent nodes closest to Uiand choose a node viva. (4) Add to E’ the arc (t+, u,) and let (e,,, ef2, . . . , erc), with (vi, uj) = e,,, be the c arcs in the cycle
thus determined. (5) Generate new trees G; = (E,, Vi) deleting in turn each arc e,, t = t,, t2, . . . , t,, not belonging to El, from E’ u (Vi, Uj). (6) Compute the pipe capacities and the flows of E; by the following procedure.
(6.1) Choose a node ulro Vi of degree 1. (6.2) Assign a flow to each arc (u,, u,,) incident to uk, i.e. for u,,~r(u~), where r(u,) is the list of all nodes in Vi adjacent to uk, such that fk,, = r(u,.). (6.3) For each arc considered in step 6.2 choose the pipe capacity as the smallest one such that dkh 2 fkh and compute the cost of the arcs such that C,, = c(dkh)lkh.
The algorithmic structure of a DSS
221
(6.4) Compute the requirement of each node U,,ET(Q) such that r(u,,) = r(u,,) + r(ur). (6.5) Remove, node ulrfrom V: and repeat steps 6.1 to 6.5 until Vi is empty. If c(G;) < c(G’) then retain Gi as the new current solution G’ and repeat steps l-6. Otherwise, remove ui from A(ui). If A(ui) is empty, repeat steps 2-6, otherwise, repeat steps 3-6. The output of this local search procedure is a network G* which is locally optimal. Its cost depends on the starting network and cannot be further improved by this elementary transformation. The computational cost of this procedure is given by the generation of 1VI * g WI w different networks, where
and w is the average number of arcs in the cycles generated in step 4. This local search technique, applied to each of the networks shown in network of length 3085 m, cost $408,000 and 30 nodes, depicted in Fig. each arc in the network depicted in Fig. 3(B). The local search has produced a further reduction in costs (8% less and a significant reduction in length (10% less than the same network). procedure took about 10 set of CPU time.
0 A
6
Fig. 3
Fig. 2, results in the same 3(A) with the diameter of than the initial network) The whole optimization
228
A.Scrom~c~e~ andR.Sozw
3.3. Hydraulic feasibility
In order to check for the hydraulic feasibility of the network, i.e. to satisfy constraint c2, one has to compute the maximum pressure along the network. If it is within a prescribed range the network is feasible. If not, the pipe diameter is increased/decreased where the pressure is too high/low and the network is checked for feasibility. In this paper the Darcy-Weisbach equation is used to compute the steady-state pressure along the arc (i, j). In particular, given the pressure P*at the node i, the pressure Pj at the node j is given by i pj = pi - dijlij
A,
&dij where I, = length of the arc (i, j), dij = diameter of the arc, Uij= and g = gravitational acceleration. Given the value of the Reynolds number
flOW
velocity,
6,
=
friction factor
Re -_----, dijrrij P where t is the specific weight and ,u is the cinematic viscosity, the value of the friction factor is given by the Colebrook relation: 6, = 0+25[lOg,,(E/376d<, 6,=[1.14-2lOg,o
+ S74/Re0’g)]2
E/dij]
if
Re>
if
Re<
10’
lo”,
where E is the relative roughness of the pipe. Note that the design technique proposed in this paper is independent of the particular flow equation used to compute the pressure. The graph depicting the behavior of the pressure along the critical path is called piezometric curve. The check of constraint (2) is performed by computing the piezometric function of the network
Critical
hydrouklc
Geodetic
lme
1000
Distance Fig. 4
CM 3
poth
229
The algorithmic structure of a DSS 70 _._,..,..
5
10
..,..........,,.......,,,.,.,.,......
15
20
25
30
35
40
45
50
Vertices Fig. 5. Execution times (for complete graphs).
and deciding from its values at the nodes whether the pressures are less than PN [13]. In particular, once the difference in pressure between the endpoints of each arc is known, one can compute the overall loss of load for all paths connecting the source to the sinks. The critical path of the network is that path for which the loss of load is largest. Once the loss of load along a path is known, one can compute at all nodes of the critical path the pressure of the fluid in both directions (from and to the source). If the largest pressure along this curve is greater than PN then the network is not feasible because it violates constraint cl. If the network is not feasible, the values of the diameters of the arcs incident to those nodes where the pressure is too large are increased until the constraint is satisfied. Alternatively, a decision about the diameters assignment is required from the designer of the district heating network under consideration. Figure 4 displays the piezometric function of the minimum cost network, shown in Fig. 3(A). The network thus obtained is hydraulically feasible since the maximum pressure along the network never exceeds the nominal given value PN = 16 kg/cm2. 4. CONCLUDING
REMARKS
The proposed DSS has a worst case computational complexity of O(m2), where m is the number of arcs of the graph, due mainly to Procedures 3.2.1 and 3.2.2. This theoretical result has been confirmed by the results of different tests, performed on graphs having lo-50 nodes, whose computational times are reported in Fig. 5. In order to assess the accuracy of the DSS, its solution for the problem under consideration has been compared with the exact solution obtained by an exhaustive approach. In particular, the algorithm proposed by Christofides in [6] for the generation of all the spanning trees of the graph G has been implemented and the cost of each tree has been evaluated using the procedure described in Section 3.1. Considering the computational time required by such an approach, it has been possible to obtain the exact solution of the problem only in the case of graphs with size of the order of that depicted in Fig. 1. These solutions coincide with those obtained by the DSS; in particular, it is possible to say that the considered district heating network problem has been solved by the DSS finding the optimal solution, reported in Fig. 3. Note that in this case the ratio between the solution times for the exact solution and the ones obtained by the heuristic procedure is of the order of 103; in fact, the exact solution has been obtained in about 10 hr, while the heuristic solution has been obtained in about 10 set on IBM PC-RT. Acknowledgernear-The
authors wish to thank the reviewer for his helpful suggestions.
REFERENCES 1. K. Fahlander and K. Svanberg. Optimization of a district heating system. In OR ‘78 (Edited by K. B. Haley). North-Holland, Amsterdam (1979). 2. K. Larsson and Y. El Mahgary, An energy model for overall optimization of a nuclear based district heating system. IAEA-AG-6215. 3. F. Archetti, A. Sciomachen and C. Vercellis, A randomized approximation algorithm for a network design problem. IAMI-CNR 84.11 (1984).
230
A.
%lOMACHeN
and R. Sozz~
4. F. Archctti, A. Sciomachen and C. Vercellis, Optimal design of a remote heating network. In Lecture Notes in Control & Information Sciences, No. 84. Springer, New York (1985). .5. M. R. Garey and D. S. Johnson, Computers ond Inteructability: A Guide fo the Theory of NP-Completeness. Freeman, San Francisco, Calif. (1979). 6. N. Christofides, Graph Theory. An Algorithmic Approach. Academic Press, New York (1975). 7. Y. F. Wu, P. Widmayer and C. K. Wong, A faster approximation algorithm for the Steiner problem in graphs. IBM Research Report RC 11041 (1985). 8. B. W. Kemingham and D. M. Ritchie, The C Programming Language. Prentice Hall, Englewood Cliffs, N.J. (1978). 9. A. U. Ulusoy and D. M. Miller, Pipeline network optimization: an application to slurry problems. Proceedings of the 8th IFIP Conference on Optimization Techniques, Wurzburg, F.R.G. (1977). 10. B. Rothfarb et al. Optimal design of offshore natural-gas pipeline systems. J. Opns Res. Sot. Am. 18, No. 6 (1970). 11. D. H. Lee, Low cost drainage network. Networks 6 (1976). 12. C. H. Papadimitriou and K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity. Prentice Hall, Englewood Cligs, N.J. (1983). 13. A. Sciomachen and R. Sozzi, Una procedura automatica per la progettazione di una rete di teleriscaldamento. Proceedings AIRO, “Ricerca Operativa ed Informatica”, Franc0 Angeli (1986).
Compurers Opns Res. Vol. 17.No. 2,pp. 221-230,1990 rights reserved
Printedin Great Britain.All
THE ALGORITHMIC STRUCTURE OF A DECISION SYSTEM FOR A DESIGN OF A DISTRICT HEATING
SUPPORT NETWORK
A. SCIOMACHEN* and R. SOZZI? Department
of Information
Sciences, University of Milano, Via Cicognara 7, 20129 Milano, Italy
(Receioed June 1987;
revised
May
1989)
Scope and Purpose--Heating networks are gaining increasing importance as a cost effective and environmentally sensible way of distributing energy for the heating of residential and commercial buildings. The algorithms discussed in this paper, included in a decision support system (DSS), have been shown to be effective in a number of cases. The planning and the design of a network is a complex process: a DSS could reduce signilicantly the time and the cost involved.
Abstract-This paper is devoted to describe the algorithms implemented in a decision support system for the design of a district heating network. The design problem is formulated on a graph as the search for the minimum cost tree rooted at the source and spanning all the sinks. This search consists of two procedures: the first one finds, using greedy-like algorithms, an approximate solution which is subsequently improved in the second procedure, using local searches.
1. INTRODUCTION
The design of a district heating network is a complex process in which different factors like costs, reliability and the environmental impact must be taken into account. For this reason computer aided design techniques have been increasingly in use to assist the designer in assessing the trade-offs among all the factors involved and identifying the final design to be implemented. The first application of operations research techniques, at least to the authors’ knowledge, is in Refs [l] and [2] where the problem of deciding the operation of the production units and the dimensions of the pipe system is addressed, while the network topology is assumed to be given. The authors of this paper have been involved in a research activity, sponsored by the Energy Project of the Council for National Research, aimed at developing a decision support system (DSS) for the design of a district heating network. The first results obtained in this research are reported in Refs [3] and [4]. In this paper the cost is assumed as the performance index of the network: in this case the task of the designer is to link, at a minimum cost, a source to a number of demand nodes or sinks. This paper is devoted to illustrating the algorithmic procedures implemented in the DSS in order to compute the cost of a given network, to check whether it satisfies the constraints and to find the network of minimum cost. The main decision variables are the layout of the network and the capacity assignments to its arcs. The constraints are of two kinds. First, the features of the area to be served may restrict the number of possible layouts and second, the network is to be hydraulically feasible, i.e. the pressure must be within some range and this restricts the choice of pipe diameters. The main task of the DSS is to assist the designer in the formulation of the objective function and the constraints, to provide an optimization procedure and to assist the designer in controlling its execution and assessing the validity of the solutions obtained. The data typically available are the energy at the source and the requirements at the sinks, a description of the area *MS A. Sciomachen received the degree in information sciences from the University of Milano. She is Assistant Professor of Operation Research at the University of Milano, and has worked on space robotics at Tecnospazio S.P.A. and on manufacturing automation at IBM Yorktown Heights, New York. Her main interests are in combinatorial algorithms and modeling & simulation of manufacturing systems, particularly network design and Petri Net based models of manufacturing systems. tMr R. Sozzi received the degree in nuclear engineering from the Polytechnic Institute of Milano. He is working on applications of operations research to energy conservation and environmental assessment at CISE (Milano) and TIM (Roma). He was also Adjunt Professor of Operation Research at the University of Milano. His main interests are in combinatorial algorithms and modeling & optimization of environmental systems, particularly network design and optimal location of monitoring networks. 221
to be served by the network, i.e. where and at what cost the pipes can be placed, and the technical and economical data of the main network components (different diameters commercially available, pumping stations). The cost optimization problem, as described in Section 2, is formulated on a graph whose arcs represent the possible layouts and whose nodes are split inta a source, a set of sinks and a set of intermediate nodes. The search for the network of minimum cost, which is hydraulically feasible, is restricted to the trees of the graph rooted at the source and spanning all the sinks. This combinatorial problem is NP-hard [5]: actually, if there is only one available arc capacity, it is a particular case of the Steiner problem f6]. For this reason approximation algorithms are to be used to solve it. An interesting approximation algorithm to solve the Steiner problem has been recently suggested in [?I: the algorithm has a very good computational performance but cannot be applied to the proposed network design problem muse the cost of each arc in the heating network is expressed not only by a distance function but depends also on the arc capacity. The algorithm used in the DSS, described in Section 3.2, is based on another approach and consists of two procedures. The first one finds an approximation to the minimum cost network and the second improves the outcome of the first procedure using local searches. It is important to note that the input to the second procedure can also be provided by the designer on the basis of any given network design independently of the execution of the first procedure. Sections 3.1 and 3.3 are respectively devoted to the procedure for computing the cost of the network and to check the hydraulic constraint. The algorithmic structure of the DSS proposed in this paper has been implemented using the C language fg] on the IBM PC-RT using the AIX Operating System. The algo~thms and the data structures used have been carefully tailored to ensure a high level of interactivity and portability to a personal computing environment. This system has been used in the design of the district heating network in a medium sized town in northeastern Italy. The computational results reported in this paper are related to this application. One can observe that the DSS can assist in the design of other types of distribution networks. The model, introduced in the next section, is simple yet general to easily represent slurry pipelines [9}, gas or oil pipeline systems [lo], or drainage systems [l 11. 2. THE NETWORK
MODEL
The network design problem is modeled on a graph G = (I$ E) with the following speci~cations: tr= (sl tr, I) is the set of n -t- 1 nodes where s represents the source, U the set of sinks and I the set of intermediate nodes. At the intermediate nodes no energy is consumed or generated. By using these nodes, one can model some constraints of the area served by the network or the possibility for the pipes to join and fork at some points to obtain a more effective layout. In some cases the set I can be empty. To each node UE V is associated a thermal load r(o) such that r(u)=O,
DEI
r(u)>O,
YEU
r(s) = - 1 r(v), RY Le. it is assumed that the energy available at the source can meet all the requirements at the sinks. E c V x f/ is the set of m arcs connecting pairs of nodes of the graph. The set E is partitioned into two sets El and E2, where El represents the set of arcs which must be in the final layout of the heating network, and E2 is the set of the arcs among which a selection takes place. As an example of this problem, the design of a district heating network of a medium sized town in northeastern Italy is presented. A first analysis of the data yielded 21 sinks which should be connected to the source, for a global energy requirement of 11 x 10’ calorie/hr, that is the global flow of about 53 l/set, assuming that the flow temperature is 120°C out of the source and 60°C back to it. Table 1 lists the requirements of the 21 sinks. The possible network layouts are represented in Fig. 1. The resulting graph has 32 nodes and 40 arcs.
The algo~~c
structure of a DSS
223
Table I
Energyrequirements Sink
ww 6.54 1.62
1 L
3 4 5 6 7 8 9 10 11 12 13 14 I5 16 17 18 19 20 21 Total
0”::
1.77 1.84 1.84 0.74 12.98 0.20 1.00 0.44 0.59 0.44 0.47 0.53 Z:Z 1.92 15.20 2.61 S2.M
1
*
SOURCE
A
SINK
l
INTCRMEOIATE NOBE
100M -
Fig. 1
The topology of the network is assumed to be a tree G’ = (V’, E’), where E’ = El u E2’ and E2’ c E2, of the graph G with the root at the source and spanning all the sinks. This tree contains a subset f’ c 1. I’ can be 0 also when f # 0. Since trees are subgraphs without closed cycles, they are a reasonable model for the topology of a least cost ~st~bution network. To each arc (i, j)~ E is to be assigned a Sow fil, which satisfies the set cr of balance equations: 1. fil-0, @,j)BE 2. ~~~Li=~ffl=O,
ioI,k, joV
3. ~fku=~ful=M), u~U,k,jeV 4. cfsI=r(s) and fjs=O, VjeV.
(1)
Relations 2,3 and 4 express respectively the balance equations at the intermediate nodes, the sinks and the source. In order for the spanning tree to represent a network, i.e. to allow a feasible now, a diameter d, is to be assigned to each arc (i, j) E E’ which allows a flow fil satisfying the baIance equations
224
A. SCIOMACHEN and R. Sozzr
and satisfying a further constraint c2: max (Pi) < PN,
(2)
isV
where Pi is the pressure at node i, i.e. the pressure in the network is less than some nominal value P,,,.
3. THE
ALGORITHMIC
STRUCTURE
The algorithmic structure of the DSS can be summarized in the following steps, where G is the initial graph, G’ is the initial tree network and G* is the final solution. (1) Check whether G is connected. If G is not connected the problem has no feasible solution. (2) Check the existence and the feasibility of the initial solution G’. If G’ = 0 there is no initial tree network. If El # 0 and El has a cycle the problem has no feasible solution. Otherwise generate G’ by Procedure 3.2.1. If G’ # 0 and G’ has a cycle the problem has no feasible solution. (3) Computation of the cost of G’. Compute the cost of G’ according to Section 3.1 and 3.2.2. (4) Local search and optimization. Call Procedure 3.2.2. (5) Check the hydraulic feasibility. Compute the pressure along G’ according to Section 3.3. If G’ is hydraulically feasible G* = G’. Otherwise enter the diameters assignment policy chosen to compute the new diameters assignment according to Section 3.3. Note that steps 2-5 can be repeated more than one time with different starting solutions generated by Procedure 3.2.1 or provided by the user. 3.1. Computation of the cost The procedure for computing the cost is the following. For the available pipe diameters, the optimal flow and the corresponding cost of a unit length pipe diameter are computed. To this cost is added, for each arc of the graph, a possible additional unit cost related to the features of the area in which the arc is located. Pipelines are available in a finite selection of standard diameters. The unit investment cost of a pipeline inclusive of the installation is given for each available diameter. These costs must be amortized over some time period to convert total capital cost to an annual cost. Any amortization scheme can be used in the design method. The cost of the network is given by C=
C
C(dij)lijXij,
(LjkE’
where xii is a Boolean variable which is 1 if the arc (i, j) is chosen in the network and 0 otherwise, c(di,) is the cost of one unit length of pipe of diameter d,, and I, is the length of the acr (i, j)EE. The term c(d,,) includes (a) the investment cost, (b) operating costs (namely pumping costs, because the flow must be pumped along the network), and (c) the cost of heat losses proportional to the diameter of the pipe, the flow temperature, and the length of the period of operation. The value d, associated with each arc (i, j) is chosen in such a way that allows the flow fij at a minimum cost c(d,). For a given flow fij one computes the cost for each available diameter and chooses that diameter ensuring the minimum cost value under the constraint that the flow velocity along (i, j) does not exceed a value V,,,,,. This procedure performed for each f gives, for each diameter, a range Lfmin,f,,,.,] for which the value d,j is optimal. The proposed network design problem could be formulated as the determination of the minimum cost spanning tree subject to constraints (1) and (2). Table 2 reports the relevant data about commercially available diameters for the example considered in this paper. Column 1 is the diameter (cm); column 2 is the lower/upper bound of the optimal flow range (l/set) computed according to the above procedure (note that the lower bound is the same value as the upper bound of the closest smallest diameter), while column 3 gives the cost (S/m). Other parameters of the problem have been fixed. The technical life of the network is 20 yr, the capital amortization rate is 5 %, the electrical energy cost is 111 $/MWhr, the value of the heat losses is 23 %/IO’ calorie/hr, the efficiency of the pumps is 0.7 and the hours of operation per year are 2000 hr/year.
The algorithmic
structure
of a DSS
225
Table 2 Diameter (cm)
Min/max flow Wsec)
20 25 32 40 50 65 80 100 125 150 200 250 300
0.0, 0.3 0.3, 0.5 0.5, 0.7 0.7, 1.1 1.1, 1.8 1.8, 2.8 2.8, 4.5 4.5, 9.6 9.6, 15.3 15.3, 25.0 25.0, 56.1 56.1, 91.3 91.3, 150.6
Cost’ (S/m) 90 93 96 LOO 104 107 111 115 122 141 171 204 244
*The costs are computed on the basis of Italian currency. quoted as l$ = 1.350 Lira.
3.2. The network search procedure
The search for the minimum cost network is performed by two procedures. The first one, described in Section 3.2.1, is a greedy algorithm which provides a first approximation of the minimum cost network. This design is improved by the second procedure, described in Section 3.2.2, which is based on local search techniques. 3.2.1. The procedure for generating the starting solution. In the procedure used to generate the starting tree network, if the arcs are not provided by the designer, they can be chosen according to their lengths. In order to generate the initial solution for the local search technique, the algorithm performs the following two operations setting initially E’ = El. (1) Sort the arcs in E2 in nondecreasing order according to their lengths. (2) Add in turn each of the arcs of the sorted E2 to E’ if it does not cause a cycle in G’, until G’ is a tree that spans all the sinks. A better solution of the proposed problem can be obtained starting from (possibly) different initial solutions to be later improved by the second procedure. A natural way of obtaining different starting solutions is to replace the deterministic selection of each arc, performed in step 2, by a randomized selection on the arcs of E2 generated by using a discrete random variable; the corresponding probability mass function can be, for example, as shown in [3], a truncated geometric distribution over the integers 1,2, . . . , IE21, where (E21 is the number of arcs in E2. In the present case, the designer provided neither a starting tree network nor a selection of some arcs to be included in the final layout, i.e. in this case El = 0. The procedure summarized in Section 3.2.1 has therefore been executed with randomized arc selection, obtaining 3 different starting solutions. Figure 2 displays the 3 networks generated with three independent executions of the procedure. Network 1 has a length of 3885 m, a cost of $498,OOOand 31 nodes. The same values for network 2 and network 3 are respectively 3435 m, %405,OOO, 31 nodes and 4095 m, $51O,OOO and 32 nodes. 3.2.2. The local search procedure. For a complete analysis of the structure of local searches and their applications in combinatorial optimization the reader is referred to [12]. In this paper the authors outline the technique as it applies to the problem under consideration. This procedure uses an “elementary transformation” of a tree T, which consists in deleting one arc of T and adding to T one arc of G. Two trees Tl and T2 of G are said to be equivalent if one can be obtained from the other by this elementary transformation. Given a network G’ = (V’, E’), its neighborhood N(G’) is defined as the set of networks G” = (Y”, E”) such that G” is equivalent to G’. The network in N(G’) of minimum cost is termed locally optimal. In order to compute this locally optimal network, as in [lo], the neighborhood N(G’) is generated performing the following operations for each element of I/‘. (1) Let V” be the set of the nodes of Y’ not yet selected and set initially V” = V’. (2) If V” # 0 choose a node V,E I/” and let Y” = V’ - ui.
A. SCIOMACHEN
226
and
R. Sozv
Fig. 2 (3) Let A(ui) be the set of the b nonadjacent nodes closest to Uiand choose a node viva. (4) Add to E’ the arc (t+, u,) and let (e,,, ef2, . . . , erc), with (vi, uj) = e,,, be the c arcs in the cycle
thus determined. (5) Generate new trees G; = (E,, Vi) deleting in turn each arc e,, t = t,, t2, . . . , t,, not belonging to El, from E’ u (Vi, Uj). (6) Compute the pipe capacities and the flows of E; by the following procedure.
(6.1) Choose a node ulro Vi of degree 1. (6.2) Assign a flow to each arc (u,, u,,) incident to uk, i.e. for u,,~r(u~), where r(u,) is the list of all nodes in Vi adjacent to uk, such that fk,, = r(u,.). (6.3) For each arc considered in step 6.2 choose the pipe capacity as the smallest one such that dkh 2 fkh and compute the cost of the arcs such that C,, = c(dkh)lkh.
The algorithmic structure of a DSS
221
(6.4) Compute the requirement of each node U,,ET(Q) such that r(u,,) = r(u,,) + r(ur). (6.5) Remove, node ulrfrom V: and repeat steps 6.1 to 6.5 until Vi is empty. If c(G;) < c(G’) then retain Gi as the new current solution G’ and repeat steps l-6. Otherwise, remove ui from A(ui). If A(ui) is empty, repeat steps 2-6, otherwise, repeat steps 3-6. The output of this local search procedure is a network G* which is locally optimal. Its cost depends on the starting network and cannot be further improved by this elementary transformation. The computational cost of this procedure is given by the generation of 1VI * g WI w different networks, where
and w is the average number of arcs in the cycles generated in step 4. This local search technique, applied to each of the networks shown in network of length 3085 m, cost $408,000 and 30 nodes, depicted in Fig. each arc in the network depicted in Fig. 3(B). The local search has produced a further reduction in costs (8% less and a significant reduction in length (10% less than the same network). procedure took about 10 set of CPU time.
0 A
6
Fig. 3
Fig. 2, results in the same 3(A) with the diameter of than the initial network) The whole optimization
228
A.Scrom~c~e~ andR.Sozw
3.3. Hydraulic feasibility
In order to check for the hydraulic feasibility of the network, i.e. to satisfy constraint c2, one has to compute the maximum pressure along the network. If it is within a prescribed range the network is feasible. If not, the pipe diameter is increased/decreased where the pressure is too high/low and the network is checked for feasibility. In this paper the Darcy-Weisbach equation is used to compute the steady-state pressure along the arc (i, j). In particular, given the pressure P*at the node i, the pressure Pj at the node j is given by i pj = pi - dijlij
A,
&dij where I, = length of the arc (i, j), dij = diameter of the arc, Uij= and g = gravitational acceleration. Given the value of the Reynolds number
flOW
velocity,
6,
=
friction factor
Re -_----, dijrrij P where t is the specific weight and ,u is the cinematic viscosity, the value of the friction factor is given by the Colebrook relation: 6, = 0+25[lOg,,(E/376d<, 6,=[1.14-2lOg,o
+ S74/Re0’g)]2
E/dij]
if
Re>
if
Re<
10’
lo”,
where E is the relative roughness of the pipe. Note that the design technique proposed in this paper is independent of the particular flow equation used to compute the pressure. The graph depicting the behavior of the pressure along the critical path is called piezometric curve. The check of constraint (2) is performed by computing the piezometric function of the network
Critical
hydrouklc
Geodetic
lme
1000
Distance Fig. 4
CM 3
poth
229
The algorithmic structure of a DSS 70 _._,..,..
5
10
..,..........,,.......,,,.,.,.,......
15
20
25
30
35
40
45
50
Vertices Fig. 5. Execution times (for complete graphs).
and deciding from its values at the nodes whether the pressures are less than PN [13]. In particular, once the difference in pressure between the endpoints of each arc is known, one can compute the overall loss of load for all paths connecting the source to the sinks. The critical path of the network is that path for which the loss of load is largest. Once the loss of load along a path is known, one can compute at all nodes of the critical path the pressure of the fluid in both directions (from and to the source). If the largest pressure along this curve is greater than PN then the network is not feasible because it violates constraint cl. If the network is not feasible, the values of the diameters of the arcs incident to those nodes where the pressure is too large are increased until the constraint is satisfied. Alternatively, a decision about the diameters assignment is required from the designer of the district heating network under consideration. Figure 4 displays the piezometric function of the minimum cost network, shown in Fig. 3(A). The network thus obtained is hydraulically feasible since the maximum pressure along the network never exceeds the nominal given value PN = 16 kg/cm2. 4. CONCLUDING
REMARKS
The proposed DSS has a worst case computational complexity of O(m2), where m is the number of arcs of the graph, due mainly to Procedures 3.2.1 and 3.2.2. This theoretical result has been confirmed by the results of different tests, performed on graphs having lo-50 nodes, whose computational times are reported in Fig. 5. In order to assess the accuracy of the DSS, its solution for the problem under consideration has been compared with the exact solution obtained by an exhaustive approach. In particular, the algorithm proposed by Christofides in [6] for the generation of all the spanning trees of the graph G has been implemented and the cost of each tree has been evaluated using the procedure described in Section 3.1. Considering the computational time required by such an approach, it has been possible to obtain the exact solution of the problem only in the case of graphs with size of the order of that depicted in Fig. 1. These solutions coincide with those obtained by the DSS; in particular, it is possible to say that the considered district heating network problem has been solved by the DSS finding the optimal solution, reported in Fig. 3. Note that in this case the ratio between the solution times for the exact solution and the ones obtained by the heuristic procedure is of the order of 103; in fact, the exact solution has been obtained in about 10 hr, while the heuristic solution has been obtained in about 10 set on IBM PC-RT. Acknowledgernear-The
authors wish to thank the reviewer for his helpful suggestions.
REFERENCES 1. K. Fahlander and K. Svanberg. Optimization of a district heating system. In OR ‘78 (Edited by K. B. Haley). North-Holland, Amsterdam (1979). 2. K. Larsson and Y. El Mahgary, An energy model for overall optimization of a nuclear based district heating system. IAEA-AG-6215. 3. F. Archetti, A. Sciomachen and C. Vercellis, A randomized approximation algorithm for a network design problem. IAMI-CNR 84.11 (1984).
230
A.
%lOMACHeN
and R. Sozz~
4. F. Archctti, A. Sciomachen and C. Vercellis, Optimal design of a remote heating network. In Lecture Notes in Control & Information Sciences, No. 84. Springer, New York (1985). .5. M. R. Garey and D. S. Johnson, Computers ond Inteructability: A Guide fo the Theory of NP-Completeness. Freeman, San Francisco, Calif. (1979). 6. N. Christofides, Graph Theory. An Algorithmic Approach. Academic Press, New York (1975). 7. Y. F. Wu, P. Widmayer and C. K. Wong, A faster approximation algorithm for the Steiner problem in graphs. IBM Research Report RC 11041 (1985). 8. B. W. Kemingham and D. M. Ritchie, The C Programming Language. Prentice Hall, Englewood Cliffs, N.J. (1978). 9. A. U. Ulusoy and D. M. Miller, Pipeline network optimization: an application to slurry problems. Proceedings of the 8th IFIP Conference on Optimization Techniques, Wurzburg, F.R.G. (1977). 10. B. Rothfarb et al. Optimal design of offshore natural-gas pipeline systems. J. Opns Res. Sot. Am. 18, No. 6 (1970). 11. D. H. Lee, Low cost drainage network. Networks 6 (1976). 12. C. H. Papadimitriou and K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity. Prentice Hall, Englewood Cligs, N.J. (1983). 13. A. Sciomachen and R. Sozzi, Una procedura automatica per la progettazione di una rete di teleriscaldamento. Proceedings AIRO, “Ricerca Operativa ed Informatica”, Franc0 Angeli (1986).