Evolutionary Algorithms for Optimal Operation of Gas Networks by Andy Jongsuk Lee A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate and Postdoctoral Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) January 2018 © Andy Jongsuk Lee, 2018 ii Abstract Natural gas is one of the cleanest fossil fuels and most reliable energy source readily available for everyday use. Due to changes in weather conditions and in economic growth and market conditions, demand for natural gas by municipal, as well as commercial and industrial consumers, has greatly increased, and delivering gas to these varied consumers has become more complex. The problem of optimally operating gas transmission networks is generally formulated to minimize the total supply cost of the network while satisfying the demand of consumers at different delivery points, at minimal guaranteed pressures. If necessary, compressors are installed and operated to supply more energy to the network and increase gas pressures, so as to compensate for pressure losses in the pipe network. The nonlinear relationships between system discharges and headlosses due to friction and mechanical devices in pipes form a complex set of nonlinear constraints in optimization models. In order to solve this complex problem, researchers have used approaches that use linear approximation in classical optimization techniques. Recently, heuristic algorithms have been applied, as these techniques can be used to find approximate solutions to complex optimization problems more quickly than classical optimization methods, and can often obtain an acceptable solution when classical methods fail. This thesis compares the performance of four heuristic algorithms for optimizing operations of a gas transmission network, namely, the: Genetic Algorithm (GA), Differential Evolution Algorithm (DE), Artificial Bee Colony Algorithm (ABC), and Covariance Matrix Adaptation Evolution Strategy (CMA-ES). The algorithms are compared in terms of their computational burden, and the quality of the obtained solutions with respect to practical management concerns. Among the four algorithms, CMA-ES is found to achieve the global optimum, compared with that generated with iii classical optimization, most consistent results in solutions that meet pressure and flow requirements, and conserve compressor power. iv Lay Summary Gas suppliers aim to efficiently transport the gas to the end users and satisfy the pressure and demand requirements, often increases in usage by residents, commercial businesses, and industries. Transporting gas from extraction sites to the customers is challenging due to long transporting routes and the need to meet pressure targets set by contracts. When gas travels along great distances, the friction of the pipe walls causes significant energy losses. Such losses in energy can lead to the need to acquire more gas in order to meet the pressure targets of the end users. Gas suppliers may save on operating costs of their transmission systems by applying a set of reliable problem-solving tools that are known generally as optimization algorithms. This thesis evaluates the performance of four of these algorithms, for finding optimal operating conditions of gas transmission systems in the face of constant demands for distributed end users. v Preface The thesis is the original intellectual product of the author, J. Lee. The development of the simulation and optimization models of the case study in Chapter 4 and the analysis of the case study results was undertaken by the author. The development of the simulation and optimization codes used to analyze the case study was undertaken by N. Moosavian and the author. vi Table of Contents Abstract .................................................................................................................................................. ii Lay Summary ........................................................................................................................................ iv Preface .................................................................................................................................................... v Table of Contents .................................................................................................................................. vi List of Tables....................................................................................................................................... viii List of Figures ....................................................................................................................................... ix List of Symbols .................................................................................................................................... xii Acknowledgements ............................................................................................................................ xvii Dedication ......................................................................................................................................... xviii Chapter 1: Introduction .......................................................................................................................... 1 Chapter 2: Literature Review ................................................................................................................. 6 2.1 Gas network problem definition ................................................................................................... 6 2.1.1 Optimal operation and design ............................................................................................... 7 2.1.2 Model differences between the operation and design problems ........................................... 9 2.2 Simulation and optimization solutions of pipe networks ........................................................... 10 2.2.1 Hydraulic simulation ........................................................................................................... 10 2.2.2 Pipe network optimization .................................................................................................. 13 2.2.2.1 Optimization of pipe networks ..................................................................................... 13 2.2.2.2 Classical optimization of pipe networks ...................................................................... 13 2.2.2.3 Heuristic optimization of pipe networks ...................................................................... 15 Chapter 3: Optimal Operations Models for Gas Networks and Methods of Analysis ......................... 18 3.1 Decision model for gas provision ............................................................................................... 18 3.1.1 Constraints ........................................................................................................................... 19 3.1.2 Formulation of gas network problem .................................................................................. 21 3.1.2.1 Fuel consumption minimization ................................................................................... 21 3.1.2.2. Gas purchase minimization ......................................................................................... 22 3.1.2.2 Least cost gas pipe network design .............................................................................. 23 3.2 Heuristic algorithms for solving decision models of optimal operations of gas networks ......... 23 3.2.1 Genetic Algorithm (GA) ..................................................................................................... 23 3.2.1.1 Coding of strings .......................................................................................................... 25 3.2.1.2 Populations of strings ................................................................................................... 26 3.2.1.3 Reproduction ................................................................................................................ 26 vii 3.2.1.3.1 Fitness reproduction or selection ........................................................................... 27 3.2.1.3.2 Crossover ............................................................................................................... 27 3.2.1.3.3 Mutation ................................................................................................................ 28 3.2.1.4 Stopping criteria ........................................................................................................... 29 3.2.2 Differential Evolution (DE) Algorithm ............................................................................... 29 3.2.2.1 Mutation ....................................................................................................................... 30 3.2.2.2 Crossover ...................................................................................................................... 30 3.2.2.3 Selection ....................................................................................................................... 31 3.2.2.4 Stopping criteria ........................................................................................................... 31 3.2.3 Artificial Bee Colony (ABC) Algorithm ............................................................................. 32 3.2.3.1 Algorithm structure ...................................................................................................... 33 3.2.4 Covariance Matrix Adaptation Evolution Strategy (CMA-ES) .......................................... 35 3.2.4.1 Initialize the parameters ............................................................................................... 36 3.2.4.2 Selection and recombination ........................................................................................ 37 3.2.4.3 Step size control ........................................................................................................... 38 3.2.4.4 Covariance matrix adaption ......................................................................................... 39 3.2.4.5 Stopping criteria ........................................................................................................... 40 3.3 Hydraulic simulators for assessing the performance of gas networks ....................................... 40 Chapter 4: Case Study: Belgium Gas Network .................................................................................... 48 4.1 Belgium gas network operations optimization model and evolutionary algorithm implementation ................................................................................................................................. 50 4.2 Optimal operations of the Belgium gas transmission network ................................................... 51 Chapter 5: Conclusions and Future Work .......................................................................................... 102 References .......................................................................................................................................... 105 viii List of Tables Table 1. Binary Code Example of GA ................................................................................................. 47 Table 2. Belgian Daily Maximum High-Calorific Gas Demand by Location in 1989 ........................ 95 Table 3. Allowable Limits of Flow and Pressure at the Nodes ............................................................ 96 Table 4. Pipe Layout Description ......................................................................................................... 97 Table 5. Parameter Values for Evolutionary Algorithms ..................................................................... 98 Table 6. Pressure Violations Resulting from Application of CMA-ES for the Stopping Criterion of 50,000 Function Evaluations ................................................................................................................ 99 Table 7. Pressure Violations Resulting from Application of GA for the Stopping Criteria of the 20,000 and 50,000 Function Evaluations ........................................................................................... 100 Table 8. Pressure Violations Resulting from Application of DE for the Stopping Criteria of the 20,000 and 50,000 Function Evaluations ........................................................................................... 101 ix List of Figures Figure 1. Gas Transmission Network for the British Columbia, Canada. (from Fortis BC, https://www.fortisbc.com/About/ServiceAreas/Pages/default.aspx) ..................................................... 5 Figure 2. Network Schematic ............................................................................................................... 42 Figure 3. Simple GA Flowchart ........................................................................................................... 43 Figure 4. Differential Evolution Flowchart .......................................................................................... 44 Figure 5. Artificial Bee Colony Algorithm Flowchart ......................................................................... 45 Figure 6. CMA-ES Flowchart .............................................................................................................. 46 Figure 7. Belgium Gas Transmission Network (Drawn on Google Earth Image and modified from DeWolf and Smeers (2000)) PLEASE ADD POPPEL ........................................................................ 55 Figure 8. The Total Cost, million $, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ........................................................................................................................................... 56 Figure 9. Pressure at Node 1 (Zeebrugge), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 57 Figure 10. Pressure at Node 2, (Dudzele), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 58 Figure 11. Pressure at Node 3 (Bruges), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 59 Figure 12. Pressure at Node 4 (Zomergem), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 60 Figure 13. Pressure at Node 5 (Loenhout), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 61 Figure 14. Pressure at Node 6 (Antwerp), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 62 Figure 15. Pressure at Node 7 (Ghent), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 63 Figure 16. Pressure at Node 8 (Voeren), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 64 x Figure 17. Pressure at Node 9 (Berneau), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 65 Figure 18. Pressure at Node 10 (Liège), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 66 Figure 19. Pressure at Node 11 (Warnant), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 67 Figure 20. Pressure at Node 12 (Namur), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 68 Figure 21. Pressure at Node 13 (Anderlues), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 69 Figure 22. Pressure at Node 14 (Péronnes-lez-Binche), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ........................................................................................................... 70 Figure 23. Pressure at Node 15 (Mons), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 71 Figure 24. Pressure at Node 16 (Blaregnies), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 72 Figure 25. Pressure at Node 17 (Wanze), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 73 Figure 26. Pressure at Node 18 (Sinsin), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 74 Figure 27. Pressure at Node 19 (Arlon), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 75 Figure 28. Pressure at Node 20 (Pétange), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 76 Figure 29. Supply at Node 1 (Zeebrugge), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 77 Figure 30. Supply at Node 2 (Dudzele), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 78 Figure 31. Demand at Node 3 (Bruges), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 79 xi Figure 32. Demand at Node 5 (Loenhout), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 80 Figure 33. Demand at Node 6 (Antwerp), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 81 Figure 34. Demand at Node 7 (Ghent), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 82 Figure 35. Supply at Node 8 (Voeren), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 83 Figure 36. Demand at Node 10 (Liège), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 84 Figure 37. Demand at Node 12 (Warnant), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 85 Figure 38. Supply at Node 13 (Namur), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 86 Figure 39. Supply at Node 14 (Péronnes-lez-Binche), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ......................................................................................... 87 Figure 40. Demand at Node 15 (Mons), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 88 Figure 41. Demand at Node 16 (Blaregnies), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 89 Figure 42. Demand at Node 19 (Arlon), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 90 Figure 43. Demand at Node 20 (Pétange), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ................................................................................................................ 91 Figure 44. Power for Compressor at Pipe 9, kW, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 92 Figure 45. Power for Compressor at Pipe 19, kW, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 93 Figure 46. Power for Compressor at Pipe 22, kW, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations ............................................................................................................................ 94 xii List of Symbols Gas Pipe Network Notation i= Node indicator j= Node indicator K= Set of pipes N= Set of nodes TK= Total number of pipes TN= Total number of nodes CS= Set of compressor stations Pi= Gas pressure at node i (bar) Pj= Gas pressure at node j (bar) Pi min= Minimum allowable gas pressure at node i (bar) Pi max= Maximum allowable gas pressure at node i (bar) Pj max= Maximum allowable gas pressure at node j (bar) Si= Net gas supply at node i (106 m3/day) Sj= Net gas supply at node j (106 m3/day) Di= -Si, demand at node i (106 m3/day) Dj= -Sj, demand at node j (106 m3/day) Qi in= Flow into the node i (106 m3/day) Qi out= Flow out of node i (106 m3/day) qi= Flow available at node i (106 m3/day) Qk= Flow in the kth pipeline (106 m3/day) Fk= Pipe characteristic coefficient of the kth pipeline (106 m3/day/bar) HPcsk= Power required for compressor in the kth pipeline (kW) λcsk= Compressor characteristic coefficient in the kth pipeline (Unitless) Rcsk= Compressor characteristic coefficient in the kth pipeline (Unitless) Bcsk= Compressor characteristic coefficient in the kth pipeline (Unitless) Pj max= Maximum allowable gas pressure at the exit end of a pipe node j (bar) Ci = Purchase cost of gas at supply node i ($/106 m3/day) xiii C(Diamk)= Specified unit cost of pipe with a corresponding diameter, Diamk, per linear length of the kth pipeline ($/mm/km) Lk= Length of the kth pipeline (km) CPF(Diamk)= Specified unit cost of a pipe fitting in the kth pipeline with a corresponding diameter, Diamk ($/mm) CCS(Diamk)= Specified unit cost of a compressor station in the kth pipeline with a corresponding diameter, Diamk ($/mm) CV(Diamk)= Specified unit cost of a valve in the kth pipeline with a corresponding diameter, Diamk ($/mm) Tgas= Gas temperature in a pipe, 281.15 (°K) δ= Density ratio of the gas relative to air at Standard Temperature and Pressure (Unitless) Z= Gas compressibility factor (-) ε= Pipe wall roughness of a carbon steel pipe (mm) Related to the Genetic Algorithm (GA) U= Solution vector T= Solution vector U’= Solution vector after crossover T’= Solution vector after crossover ul= The value of the lth decision variable in a given solution vector Related to the Differential Evolutionary Algorithm (DE) E= Decision Variable NP= Number of decision variables in a solution vector, used in mutation g= Current generation indicator for DE and CMAES g+1= Next generation indicator for DE and CMAES vm, g+1= A mutant variable of E in generation g+1, which is a combination of Er1,g, Er2,g, and Er3,g Er1,g= Randomly selected first variable value of E within a solution vector at generation g Er2,g= Randomly selected second variable value of E within a solution vector at generation g Er3,g= Randomly selected third variable value of E within a solution vector at generation g R= A mutation factor that is a real number that is greater than 0 and less than 2 xiv CO= Number of decision variables in a solution vector, used in crossover uz,g= The zth trial solution vector at generation g uz,g+1= The zth trial solution vector at generation g+1 umz, g+1= mth variable value for the zth trial solution vector for generation g+1 randb(s)= Randomly generated number ranging from 0 to 1 CR= Crossover threshold value ranging from 0 to 1 Related to the Artificial Bee Colony Algorithm (ABC) b= Indicator of a food source SN= Number of food sources e= Indicator of optimization parameter OP= Number of optimization parameters h= Indicator of employed bees BN= Number of employed bees Cy= Indicator of a cycle of the search process Cymax= Number of repeated cycles in the search process Prb= The probability of attracting onlooker bees to a given food source fb = Fitness value of bth food source posbe= Current food source position or current value of the decision vector newposbe= New food source position or a new value of the decision vector ɸbe= A random number that enhances the random search, is between -1 and 1 Related to the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) g= Current generation indicator for DE and CMAES g+1= Next generation indicator for DE and CMAES xa(g+1) = ath offspring of generation g+1 M(g)= Mean value of the distribution of decision vectors generated in the search that takes place in generation g M(g+1) =Mean value of the distribution of decision vectors generated in the search that takes place in generation g+1 xv σ(g)= Standard deviation of the distribution of decision vectors generated in the search that takes place in generation g σ(g+1) = Standard deviation of the distribution of decision vectors generated in the search that takes place in generation g+1 C(g)= Covariance matrix of the distribution of decision vectors generated in the search that takes place in generation g C(g+1) = Covariance matrix of the distribution of decision vectors generated in the search that takes place in generation g+1 ps= Population size or number of offspring in a given generation, which is greater than or equal to 2 Cm= Learning rate coefficient of the covariance matrix for updating the mean value of the distribution from one generation to the next µ= Selected best individuals in a given generation µeff= Variance effective selection mass for the mean, using the method of weighted intermediate recombination, where 1 < µeff < µ wτ= Positive weight coefficients of the method of weighted average xτa(g+1) = τth-best offspring from the initial sample space Cσ= Learning rate coefficient for the accumulation of the step size control, < 1 Dσ= Damping parameter for the step size updates N(0,I)= Expected value of a multivariate normal distribution with zero mean and unity covariance matrix Pσ= Length of an evolution path Pc(g+1) = Evolution path at generation g+1 Pc= Evolution path at generation g Cc= Learning rate coefficient for the accumulation of the rank-one update of the covariance matrix C(g+1) = Covariance matrix at generation g+1 C(g)= Covariance matrix at generation g Cone= Learning rate coefficient for the rank one update of the covariance matrix Cµ= Learning rate coefficient for the rank -µ update of the covariance matrix Related to the Global Gradient Algorithm (GGA) A10 = Matrix representing the relationships between source nodes and pipes linkage A11= Diagonal matrix representing the relationships between nodes containing equipment including fittings, compressors, valves, and pumps and specified pipelines xvi A12= Matrix of nodes and pipes linkage (the topology of the network) where, {-1=if the flow of the pipe leaves the node0=if the pipe is not connected to the node+1=if the flow of the pipe enters the node; A21= Transpose of A12 matrix DIAG-1= Inverse of a diagonal matrix of flow equation in pipe or pump Q= Unknown pipe flows at each node (106 m3/day) Qq= gth iteration of flow matrix Qq+1= Next iteration flow head matrix H= Unknown pressure head at a given node (bar) H0= Known pressure head at a source node (bar) Hg+1=Next iteration pressure head matrix Du= Known demand flow at each node (106 m3/day) xvii Acknowledgements I would like to thank my supervisor Dr. Barbara Lence, for her guidance and tremendous input throughout the process. I would like to thank Naser Moosavian, Ph.D. student, for his valuable support throughout the process. I would like to thank my second reader Dr. Ziad Shawwash for his time and comments regarding the review of the thesis. I would also like to thank my wife, Angela Jihye Kim for all her love and encouragement during my two years of studies at University of British Columbia. Lastly, I would like to thank my colleagues for their support. xviii Dedication To my wife, Angela Jihye Kim, my parents in law (Gregorius and Marianna), my grandmother Catharine (Suk kyu Lee) and my Edmontonian family (Jerome, Theresa, and Lillian). 1 Chapter 1: Introduction Natural gas is one of the cleanest fossil fuels and most reliable energy sources readily available for everyday use. It is used for cooking, heating, fuel for vehicles, electricity generation, and producing materials and chemicals. It is also known to be one of the most versatile forms of energy compared with other fossil fuels. It is colorless, odorless and tasteless, is made of compounds of hydrogen and methane, or simply methane, and is extracted from both offshore and onshore sources, i.e., reservoirs, such as those in deeper sections of the ocean bed. When natural gas is extracted from the ocean bed, it is usually in liquid form. Separation processes are employed at the underground storage fields or by distribution companies to convert the liquid gas to a gaseous state that is considered to be consumer-grade natural gas. It is then transported to consumers, often through transmission pipe networks that connect sources and land delivery points. Due to changes in weather conditions and in economic growth and market conditions, demand for natural gas by municipal, as well as commercial and industrial consumers have greatly increased, and delivering gas to these varied consumers has become more complex. This thesis compares the performance of four heuristic algorithms for optimizing operations of a gas transmission network. The algorithms are compared in terms of their computational burden, and the quality of the obtained solutions with respect to practical management concerns and reliability in the face of uncertainty. A gas pipe network is comprised of two different types of systems: transmission networks that transmit high-pressure gas from supplies to regional distributors; and distribution networks that distribute low-pressure gas to regional demand sites. Transmission networks differ from distribution networks in terms of several aspects. Relative to distribution networks, transmission systems: use larger pipes; employ more pressure valves, compressors, and nozzles; are branched in topology; and operate at higher pressures. They are usually more complex than distribution networks due to their 2 use of complex mechanical equipment and use larger pipes to efficiently carry a sizeable amount of pressurized gas over long distances. They are often laid out as branch networks. An example transmission network for British Columbia, Canada, is shown in Figure 1. This transmission network connects northern Alberta, the source and exporting region, to the southern British Columbia, the location of the demand sites in the region, and is long and branched. It is similar to the Belgian Gas Network, which is the case study for this thesis. Like these systems, most transmission networks are extensively long in distance, usually cross provincial, state, and national boundaries, and connect to few demand hubs, such as at major border crossings and large cities. In contrast, a distribution network usually consists of small diameter pipes with no mechanical equipment and is operated at low and medium pressures. The complexity of distribution networks is very similar to water distribution networks that have no chlorine booster stations, both of which have been extensively analyzed and optimized. A distribution network is often laid out as a looped network or a combination of loops and branches. Compared with Figure 1, the distribution network is denser or clustered, has many short distance pipes, contains several pipe crossings, and many intra-connections between pipes. These pipelines are usually installed along roads, freeways, and highways and thus have a layout that is similar to a road map of the region. While many investigations have been conducted regarding the application of systems analysis techniques for optimizing the design and operation of distribution networks, relatively fewer have been conducted for transmission networks. However, optimization of transmission networks has shown the greatest promise for exacting cost savings because these networks employ equipment that is considerably more expensive than those used in distribution networks, in terms of both capital and operating costs. 3 As the usage of natural gas as a source of energy is increased, the benefits of optimizing the handling and processing of this substance should also increase. Consumers in inland regions rely on imported natural gas from other countries, in particular, stand to benefit from the efficient allocation of expenditures for purchasing natural gas and for operating pipe transmission networks more efficiently. If such regions have pre-existing transmission networks, they may be facing the decision of whether to replace or expand their existing pipe network or to implement a strategy that improves its operation. The former is considered a design, or capacity expansion, systems problem, and the latter is considered a systems operation problem. The problem of optimally operating gas transmission networks is generally formulated to minimize the total supply cost of the network while satisfying the demand of consumers at different delivery points, at a minimal guaranteed pressure. If necessary, compressors are installed and operated to supply more gas to the network and increase gas pressures, so as to compensate for pressure losses caused by the pipe hydraulic friction. Compressor stations are installed along the pipelines and are powered by gas taken directly from the network gas flow. The nonlinear relationships between system discharges and headlosses due to friction and mechanical devices in pipes form a complex set of nonlinear constraints. These constraints can be found in any pipe network, however, the option of using compressors to compensate for headlosses in gas networks adds variable demands to the system, changes the nature of the optimization problem, and makes it more difficult to solve. The relationship between pressure increases due to compressor operations, and resulting pipe flows, is a nonlinear function, which creates a nonlinear and nonconvex feasible solution region. In order to solve this complex problem, researchers have used approaches that link linear approximation with classical optimization techniques, such as linear, dynamic, nonlinear, and successive-linear programming. Recently, heuristic algorithms have been applied, as these techniques can be used to find approximate solutions to complex optimization problems more quickly than 4 classical optimization methods, and can often obtain an acceptable solution when classical methods fail to provide sufficient solutions. In this thesis, four different heuristic algorithms are investigated for determining the optimal operation strategy for a gas transmission network. These algorithms are the Genetic Algorithm (GA), Differential Evolution Algorithm (DE), Artificial Bee Colony Algorithm (ABC), and Covariance Matrix Adaptation Evolution Strategy (CMA-ES). The algorithms are applied to the Belgian Gas Network and compared in terms of the quality of the solutions obtained, and the necessary computational effort. The research questions addressed in this thesis are: 1) Among the heuristic algorithms most widely used for the design of networks, how efficient are these algorithms for optimizing gas transmission operations? 2) Do these heuristics provide consistent optimal designs, or do their solutions differ from one another? 3) What insights may be gained from comparing their solutions in terms of cost, reliability, and robustness? The next chapter provides a synthesis of the literature related to applications of optimization to gas transmission problems. Next, the models which are typically used to optimize network design and operations are presented along with a brief overview of the GA, DE, ABC, and CMA-ES optimization heuristics, and the Global Gradient Algorithm (Todini and Pilati, 1987) for conducting hydraulic analyses, which are employed in this project. In Chapter 4, the case study based on the Belgian Gas Network is presented, and the results of the application of the heuristics for optimizing the operation of this network are presented. Finally, in Chapter 5, conclusions and suggestions for future work in the application of optimization heuristics for gas network operations are presented. 5 Figure 1. Gas Transmission Network for the British Columbia, Canada. (from BC Oil and Gas Commission, http://www.bcogc.ca/node/11470/download) 6 Chapter 2: Literature Review The field of gas network optimization includes the application of simulation and decision analytic tools for designing the layout and operation strategies for existing and proposed networks. This literature review focuses on studies undertaken to: 1) formulate the gas network operation and design optimization models; 2) simulate the hydraulic conditions resulting from proposed or existing gas network system designs; 3) implement classical and heuristic optimization algorithms, that link the formulated optimization models with hydraulic analyses of the given gas network, to determine the optimal system design or operating strategy; and 4) evaluate the performance of an array of methods for undertaking items 2-3. 2.1 Gas network problem definition Similar to other networks of conduits for distributing incompressible fluids, gas pipe network designers and operators face concerns related to optimizing the distribution of materials, meeting demands while satisfying the physical limitations and other contractual obligations of the network, and reducing energy losses in the pipes, machines, and fittings of the gas network during its operation. Proper design and operation approaches may save a great deal of money for the company that is distributing and selling the gas to consumers. The gas network optimization problem consists of two main categories: operation and design. Depending on the needs of the users, and the available hydraulic information for the system, the gas network optimization problem may focus on only one of these categories or both of them simultaneously. 7 2.1.1 Optimal operation and design The optimal operation of a gas network is generally undertaken in a manner which maximizes revenue to the gas company, while also maximizing efficiencies and minimizing operating costs of the gas network. For example, the common objectives for the operation of a given network which may be desired solely, or in combination, are: minimization of operating costs, minimization of fuel costs, maximization of line pack, and minimization of carbon dioxide emissions. The goal of the operating cost minimization problem is to find the lowest total operating costs, often including the cost of energy lost in the pipes and valves, and compressors. Fuel cost minimization is directly related to compressor operations in a network. Compressors are essential to the operation of the network for restoring pressure levels in the face of pressure losses due to friction and in order to maintain desired pressures at delivery points. They boost gas pressures in a network and in turn use some fuel to run. According to Tabkhi (2007), the operating cost of compressors can be as much as 25% to 50 % of the total company’s operating budget. Therefore, minimizing fuel costs of compressors can make a significant impact on the operating cost of a gas company. Wu et al. (2000) and Rios-Mercado et al. (2002) propose and investigate a model reduction technique that minimizes the fuel consumption of compressor stations in steady-state transmission networks. The method defines and maximizes the line pack and thereby significantly reduces the dimensions of the optimization model without disrupting the core problem structure. That is, it maximizes the volume of gas in the pipe network between the compressor that is providing a given discharge pressure and the delivery point that is providing gas at a given pressure. The length of pipe between the compressor and the delivery point can represent a temporary storage unit when the gas service to the customer is inactive, and the technique of maximizing the line pack, or volume of gas, in this length of pipe can provide a reserve of fuel to withstand peak hour usage and to manage supply and demands efficiently. Carbon dioxide emission minimization may be a key approach to consider as the world switches to more environmentally friendly operations that reduce the causes of global warming. The network may be 8 located in a country having strict laws limiting the emission of carbon dioxide. In these cases, companies may face disincentives to emit carbon dioxide, such significant legal pressures, carbon taxes, and/or hefty industrial penalties, and thus be encouraged to make carbon dioxide emission minimization a primary objective. Optimal design of a gas network may involve selection of the optimal configuration, or topology, and pipe sizes of the system before construction, but more often than not, it is employed for redesigning an existing gas network for capacity expansion due to increases in demand, or optimizing the size and composition of a pipe that replaces an existing damaged or end-of-life pipe. Capacity expansion may be accommodated by redesigning the topology of a network where selected pipe sizes may be increased, or by selecting pipes for resurfacing for improved performance, or modifying the number and locations of compressors. Olorunniwo (1982) formulate the gas network capacity expansion problem as a minimization of investment planning outcomes. He summarized the design and capacity expansion of natural gas transmission networks as having four main elements: 1) optimal layout of the network, 2) optimal operation of the static system, 3) determination of the steady-state flow and pressure characteristics, and 4) optimal dynamic design which strategically expands the network over time with given steady-state flow and pressure information. The optimal design of the gas network has also been formulated as a maximization of the net present value of the network while minimizing the total investment and maintenance costs. This model has been used to identify the characteristics of compressor stations such as the optimal location and type of compressors, and optimal pipe dimensions, length and topology to be installed. Babu et al (2003), develop a nonlinear optimization model to minimize the total investment cost of transmission networks, considered the yearly operating and maintenance costs of compressors, and solve the model using the DE Algorithm. From an asset management perspective, end-of-service-life pipes or broken pipes need to be replaced with new pipes, and the replacement strategy generally considers pipe conditions, 9 installation dates, and environmental and operating conditions. Replacing one or more pipes can affect the entire system and its operation, due to the hydraulic properties of the new pipe configuration. In addition, selecting an optimal size of the pipe to install also requires a great deal of effort. The most common goal of design optimization of a gas network is to minimize investment and maintenance costs while meeting all operating service requirements and constraints. Models that optimize the network design and operation simultaneously have been developed and these are generally multi-objective and complex. The goal of the design and operation problem of a gas network is to minimize both investment cost and operating costs of a network simultaneously, while also accounting for the tradeoff that exists when optimizing each of these costs. Babonneau et al. (2012) formulated the optimal operation and design model of the natural gas system as a joint objective problem where the objective minimizes the investment cost and energy cost at the same time. Uster and Dilaveroglu (2014) developed a model for on optimizing a new natural gas transmission network or expanding an existing network, using an integrated large mixed integer nonlinear optimization model that selects design and operation variables that determine the location and size of pipelines in the network, compressor stations and their capacities, timings of these installations in a multi-period planning horizon, natural gas purchase and steady-state flow decisions, for each period, that minimize the total investment and operating costs. 2.1.2 Model differences between the operation and design problems The most significant difference between the design and operation problems is the formulation of the objective function. The variables that are optimized for design are generally the pipe diameter and configuration, and those that are optimized for operation are the volume or flow quantities and pressures at the nodes of the network. The design is usually focused on the location and size of compressors in the pipelines and minimizing the investment cost which can be viewed as a cost of 10 installation and maintenance of facilities and instrumentation. The costs of installation and maintenance are generally expressed with functions that represent fixed and variable costs. On the other hand, the operating costs are generally related to the energy loss in pipes and the fuel cost of compressors which are variable costs. The operating problem focuses on flows in pipes to satisfy demand requirements without changing the existing topology of the network, while the design problem considers bettering the topology of a network. 2.2 Simulation and optimization solutions of pipe networks 2.2.1 Hydraulic simulation In order to optimize a pipe network system, it is essential to understand the hydraulic behavior of the pipe networks. Hydraulic simulation models are typically used to estimate nodal pressures and pipe flows, considering such input as head pressure in reservoirs, nodal demands, and characteristics of pipes, pumps, and valves. These simulation models solve for the coupled linear continuity equations at nodes and nonlinear energy equations along pipes using one of five formulations of these equations. These are the: 1) ΔQ equations, for which the unknowns are the flows in pipes and these are corrected based on Taylor Series Expansions of the energy equations for the pipes in the loops, in terms of the flows in the pipes; 2) ΔH equations, for which the unknowns are nodal pressures and these are corrected using equations based on Taylor Series Expansions of the energy equations for the pipes in the loops, in terms of nodal headloss; 3) H equations, for which the unknowns are nodal heads; 4) Q equations, for which the unknowns are pipe flows; and 5) Q-H equations, for which the unknowns are both nodal pressures and pipe flows. The Hardy-Cross (1936) approach is the first systematic iterative solution procedure and is based on the ΔQ, or loop-flow correction, equations. For each successive loop in the network, flow values for each pipe are assigned so that continuity is satisfied, and the net energy losses are 11 evaluated for this flow distribution. The flow rates are then adjusted incrementally to balance the energy losses in the loop so that the net energy loss is approximately zero. The ΔQ equations estimate the incremental changes in flow and are based on the ratio of the first-order Taylor Series Expansion of the energy equations for the pipes in the loop, and the second-order Taylor Series Expansion of the energy equations for these pipes. This process is repeated until the net energy loss in the loop is approximately zero. The method then proceeds to each subsequent adjacent loop in the network, where the flow in one pipe member of the loop is pre-defined, at the beginning of the loop analysis, based on the analysis for the previous loop. Once the net energy losses in each loop are driven to approximately zero, the entire network is evaluated in terms of whether continuity is satisfied at each node, and net energy losses are approximately zero for each loop. Cornish (1939) applies the same principal but rearranges the problem into nodal head correction equations, or ΔH equations, in a method which is called a method of mass balance. The ΔH equations are similar to the ΔQ equation, in that they are based on the ratio of the first-order Taylor Series Expansion of the energy equations for the pipes in the loop in terms of nodal headloss, and the second-order Taylor Series Expansion of the energy equations for these pipes. These two approaches are both known as the Hardy-Cross Method and are commonly used for analysis of looped networks. Because both approaches rely on iterative sets of assumptions about flow and flow direction, and estimate corrections for one loop at a time, they require a high computational effort, making the method inefficient for analyzing large sized networks. While modified Hardy Cross Methods that significantly reduce the number of iterations and improve the convergence rate have been developed, these have limited applicability to the very large sized network. The Node Equations formulation is developed by Martin and Peters (1963). They initially allocate flows for one loop of the network, based on continuity of flow at the nodes, and apply the Newton-Raphson Method to solve the set of simultaneous nonlinear equations that describe the headloss in pipes in a loop, the H equations. They find an initial solution of the H equations and then 12 correct the nodal heads in successive iterations so that the changes in headloss values in pipes, from one iteration to the next, are minimized. The resultant flow is then assigned to the given pipe, and the method proceeds to subsequent adjacent loops in the network, where the flow in one pipe member of the loop is pre-defined, at the beginning of the Newton-Raphson application, based on the analysis for the previous loop. The disadvantage of using this method is its slow convergence, especially for large sized networks. Similar to solving the H equations, Wood and Charles (1972) develop the Q equation method which defines headlosses in pipes in terms of flow. They initially allocate flows for one loop of the network, based on continuity of flow at the nodes, and apply the Newton-Raphson Method to solve the set of simultaneous nonlinear equations that describe the headlosses around loops, in terms of flow. They also solve these equations by linearizing the nonlinear relationships between the pipe resistance coefficients and the headloss in a pipe. The method proceeds to subsequent adjacent loops in the network, where the flow in one pipe member of the loop is pre-defined, at the beginning of the Newton-Raphson application, based on the analysis for the previous loop. Todini and Pilati (1987) introduce the Global Gradient Algorithm (GGA) in which they employ the Q-H equations, of continuity of flow at the nodes and conservation of energy in the loops, for the entire network, and solve these equations simultaneously with a Global Gradient search method. The Q-H equations are described in a matrix format that relates nodes and connecting pipes. The GGA is initialized with a set of trial flow solutions for the entire network, for which continuity holds, and searches for improvements in these flow allocations, that result in the net headloss in all pipe loops being driven as close to zero as possible. The algorithm searches for improvement from each of the trial solutions in parallel. Being a global search algorithm, the GGA both explores dissimilar solutions and exploits solutions that indicate improvement. Recently, Giustolisi et al. (2009) developed an improved version of GGA, which is called Enhanced GGA and allows for 13 correcting headloss error due to representation of uniformly distributed demands as nodal concentrated withdrawals. 2.2.2 Pipe network optimization 2.2.2.1 Optimization of pipe networks Once the optimal network operation and/or design problem has been formulated, an optimization solution algorithm must be applied. Two types of optimization algorithms, classical and heuristic optimization approaches, have been investigated in their application to water supply networks, gas networks and road networks. However, the most widely applied algorithms to optimization of gas networks are classical optimization methods that use an exhaustive search from one trial solution, and differential calculus to determine the globally optimal solution. Examples of classical optimization techniques include Dynamic, Linear, Nonlinear, and Mixed Integer Programming. In contrast, heuristic optimization uses classical optimization techniques, along with stochastic elements that identify good trial solutions to exploit and good areas of the solution space to explore, oftentimes based on behaviors that occur in nature, to solve for the local optimal solutions and select the best among these. Examples of heuristic optimization techniques include GA, DE, Simulated Annealing, Particle Swarm Optimization, ABC and Ant Colony Optimization. 2.2.2.2 Classical optimization of pipe networks Classical optimization solutions of pipe network design first appeared in the water distribution system literature in the early 1960’s. Shamir (1964) develops a benchmark example network and demonstrates the application of gradient like search technique in the design of a water distribution system for this network, where the decision variables are the design pipe diameters. For a design optimization of a pipe network, the objective function is considered to be a function of the 14 loading condition (i.e., demand on the system) and operating costs as a function of the energy loss in the network. Shamir and Howard (1968) reformulate the optimization model and conduct a sensitivity analysis using the Jacobian of the flow solution. The gas pipe network optimization problem was formulated as a Dynamic Programming model by Larson and Wong (1968) and since then, Dynamic Programming has been one of the most popular methods applied in the optimization of gas pipeline operations (see, e.g., Larson and Wong., 1968; Cheesman, 1971; Martch and McCall, 1972; Olorunniwo, 1981; Oloruniwo and Jensen, 1982;). Dynamic Programming has an advantage of reaching a guaranteed global optimum and can easily be applied in reformulations of complex nonlinear problems. The disadvantages of Dynamic Programming, however, include its high computational burden and its limitations to being only applicable for non-cyclic networks. Larson and Wong (1968), use Dynamic Programming to determine optimal pipe network solutions for systems operating under both transient and steady-state flows. For the transient problem, they develop an example of a single compressor in a long pipeline where the optimization problem is to minimize the compressor energy required to satisfy a known time-varying demand and pressure requirements. For the steady-state problem, they develop an example of a long series of many connecting compressors and pipelines where the optimization problem is to minimize the total compressor energy required to meet flow restrictions, pressure requirements, and constraints regarding compressor operation. They formulate the steady-state problem as a multi-staged decision and solve it using Dynamic Programming. Martch and McCall (1972) use the same technique as Larson and Wong (1968) to minimize the annual cost of owning and operating a given network. Olorunniwo (1981) and Jensen (1982) were the first researchers to apply gas pipe network optimization to an existing natural gas network, and develop the optimal design and capacity expansion planning models to identify: “1) The maximum number of compressor stations that will ever be required during the time horizon, 2) The optimal locations of these compressor stations, 3) The initial construction dates of the stations, 4) The optimal schedule of expansions of the compressors, 5) The optimal diameter sizes of the main 15 pipes for each arc of the network, 6) The minimum recommended thickness of these main pipes, 7) The optimal diameter sizes and thickness of the required loop pipes on each arc of the network, 8) The timing of construction of these loop pipes, 9) The fraction of the length of the existing pipe to be looped at any time in the planning horizon, and 10) The operating pressure in the system.” They consider two sets of constraints, these are: “1) Demands for gas are met at all times within the time horizon, and 2) Pressure and other engineering specifications are never violated.” O’Neill et al. (1979) introduce a benchmark design problem which has a transmission pipeline with one compressor problem and solves this problem with Successive Linear Programming. The focus of the study is to minimize the cost of operation with linear and nonlinear constraints. De Wolf and Smeers (1996) introduce a Piecewise Linear Programming approach and demonstrate it with the benchmark design problem developed by O’Neill et al. (1979). De Wolf and Smeers (2000) propose a two-step procedure for the piecewise linear approximation of the solution of the operations problem for the Belgian Gas Network. The problem is solved with Mixed Integer Nonlinear Programming. Given newer approximation methods, a comparison study was also undertaken. They conclude that Piecewise Linear Programming converges much faster than the other approximation methods examined, which is successive Linear Programming. Test example networks with 24, 34, and 60 arcs are analyzed for the effectiveness of piecewise linear approximations of the nonlinear constraints linked with an extension of the Simplex Algorithm. 2.2.2.3 Heuristic optimization of pipe networks In recent decades, many heuristic optimization techniques have been used for water distribution network design and operation optimization problems. Simpson et al. (1993) are the first to apply GAs to these problems. GAs are based on the mechanics of natural genetics, whereby successive generations of a population are conceived, born, and raised until they are ready to reproduce, Holland (1975). Since most water resources researchers were applying classical 16 optimization techniques at that time, heuristic optimization applications to the water pipe networks were implemented very gradually, though are widely used today (Maier et al. 2015). On the other hand, GAs were first applied, along with rule learning, in engineering practice in the field of gas network optimization, by Goldberg (1983), a Ph.D. thesis written under the supervision of Holand. Goldberg (1983) focuses on gas dispatching since determining the market demand for gas and adjusting the operations of a pipeline system to meet demands is very important for gas industries. He also applies GAs to Larson and Wong’s (1968) problem formulation for his case study. However, application of this and other heuristic algorithms did not find their way into the gas network optimization literature until Babu et al. (2003) successfully applied the DE Algorithm to the optimal design of a gas transmission network with a complex nonlinear optimization model. They compare their final solution(s) with solutions using Nonlinear Programming and a Branch and Bound Mixed Integer Nonlinear Programming Algorithm. Using DE, they successfully identify an optimal solution with less computational time compared with the other methods. As mentioned in Section 3.1.2.1, Rios-Mercado et al. (2003) propose using the Tabu Search Algorithm for minimizing fuel consumption of compressor stations in a network. Surprisingly, the Tabu Search Algorithm was developed much earlier than GAs though it only came into use at a much later date. de Mélo Duarte et al. (2006) propose using a Tabu Search Algorithm for the optimization of gas distribution networks to find the least cost pipe design while satisfying the minimum pressure constraints at each node. They compare their results with those identified with a GA, from Surry et al. (1995), and two other versions of the Tabu Search Algorithms, from Cunha and Ribeiro (2004). Chebouba et al. (2009) analyze the optimal operation of a steady flow gas pipeline using the Ant Colony Optimization Algorithm. They focus on identifying the number of turbo compressors and the discharge pressure for each compressor and compare their results with those identified with Dynamic Programming. El-Mahdy et al. (2010) use a GA with binary coding to optimize the pipe sizes in a gas network. A penalty factor is introduced to allow solutions that violate soft constraints, such as the minimum 17 permissible pressure constraints, to remain in the population during the solution process, while the GA is guided to convergence and a minimum network cost is identified. 18 Chapter 3: Optimal Operations Models for Gas Networks and Methods of Analysis 3.1 Decision model for gas provision The formulation of a gas network problem addressed in this study is viewed from an operating management perspective and focused on the optimal operation of the gas network that maximizes the net benefits to one or more supply companies. Gas supply, demand, and sale prices, as well as pipe and mechanical device conditions, must be considered. The network is controlled by the gas supplying company, whom, in turn, purchases certain quantities of gas from several sources in order to satisfy the demand and required pressures at consumer points located throughout the network. The network consists of several supply points where gas can be purchased or input to the system and several demands points where gas flows out of the system. Possible intermediate points can exist within the system at which compressors may be located, the network may be serviced, potential supply points may be located, and gas may be rerouted within the network. The gas transmission network operations model consists of arcs, or pipelines, and nodes. A pipeline is described as a connection between two nodes, denoted as nodes i and j. The set of nodes in a network is denoted by N and the set of pipelines is denoted by K. If a pipeline links to a node in which a compressor is located, the compressor may be activated in order to satisfy the minimum contracted pressures at downstream supply and demand nodes, should the contracted pressures be violated. The set of compressors in a network is denoted by CS. There are two state variables associated with each node, and these are Pi which represents the gas pressure at node i, and Si which represents the net gas supply at node i. When the value of Si is positive, a supply of gas is entering node i, and when it is negative, gas is demanded, or withdrawn at node i. A demand at node i is also denoted as Di, and is defined to be equal to -Si. A flow schematic is shown in Figure 2. 19 3.1.1 Constraints The constraints of the gas network are similar to those in water distribution. Flow continuity at nodes and energy conservation in pipeline loops are the most common constraints in pipe networks. However, there are more constraints associated with nodes and pipelines in a gas distribution network than in a water distribution network, including the contracted flow rate and pressure requirements at each node, and the power limits and flow requirements of compressors. For each node, flow continuity must be satisfied, which is formulated as: ∑ Qi in- ∑ Qi out=qi, ∀ (i)∈ N (1) where qi= flow available at node i (106 meters3/day); N= a set of nodes; and Qi in and Qi out = flow into and out of node i (106 meters3/day), respectively. If qi is a negative value, demand exists at node i and Di is positive, and if qi is a positive value, supply exists at node i and Si is positive. At a supply node i, the gas inflow must satisfy the minimum mass flow rate that is required, or contracted, by consumers. The contract is set to specify an average daily amount of gas to be supplied to each consumer, and usually, a range of the required gas quantities is specified. Similar to the demand requirements at demand nodes, the gas amount, and pressure required at, or exiting, supply nodes is also specified by contract. For continuity to hold, the gas supplied to any node must be greater than or equal to the amount demanded. Similar to the demand requirements at demand nodes, the amount and pressure of the gas exiting supply nodes can also be specified as a range of minimum and maximum allowable values. For each node, the supply quantity per day must be satisfied, such that: Si min ≤ Si ≤ Si max, ∀ (i) ∈ N (2) 20 where Si= net gas supplied at node i (106 meters3/day); Si min = minimum allowable net gas to be supplied at node i (106 meters3/day), set by contract; and Si max = maximum allowable net gas to be supplied at node i, (106 meters3/day), set by contract. For each node, the gas pressure requirements must be satisfied, such that: Pi min ≤ Pi ≤ Pi max, ∀ (i) ∈ N (3) where, Pi= gas pressure of a node i (bar); Pi min = minimum allowable gas pressure at node i (bar) set by contract; and Pi max= maximum allowable gas pressure at node i (bar) set by contract. The energy relationship in a pipe can be described as a function of gas flow in a pipeline and the pressure at both ends of the pipeline, and is formulated as: Qk2=Fk2(Pi2-Pj2), ∀ (i,j)∈ K (4) where Fk= coefficient that describes pipe characteristics such as the length, diameter, and roughness of the contacting surface of the kth pipeline, and the gas composition (106 meters3/day/bar); Qk= flow in the kth pipeline, note that Qk are unrestricted in sign (106 meters3/day); Pi= gas pressure at node i (bar) at one end of the k th pipeline; Pj= gas pressure at node j (bar) at the other end of the kth pipeline. According to Goldberg (1989), for a compressor, the pressure and flow relationship can be obtained by the adiabatic compression of an ideal gas and is formulated as. HPcskQk=λcsk* (PiPj)Rcsk-Bcsk, ∀ (i,j)∈ K (5) where HPcsk= power required for compressor (horsepower); Qk = flow in the kth pipeline (106 meters3/day); λcsk, Rcsk, and Bcsk are characteristic coefficients that describe the compressor in the kth pipeline (unitless). 21 For a compressor, there is an upper limit on the allowable pressure at the exit of the compressor, and this is formulated as: Pj ≤ Pj max, ∀ (j)∈ N (6) where, Pj= gas pressure at the exit end of the pipe (bar); Pj max=the maximum allowable gas pressure at the exit end of the pipe (bar). 3.1.2 Formulation of gas network problem The most widely accepted gas network problems optimize network operations either to minimize fuel consumption of the network, or to minimize gas quantities required to meet the demand at specified network nodes, given the supply and pressure constraints and compressor operation limitations described in Equations 1 through 6 in Section 3.1.1. The choice of objective influences the: formulation of the objective function and constraint set; effectiveness of the algorithms used to solve the models; and assessment of the performance of the solution for addressing the needs of the model user. The following section summarizes these objective functions and their corresponding constraints for the operation of the gas network. 3.1.2.1 Fuel consumption minimization According to Luongo et al. (1989), the fuel and operating cost of running the compressors represents between 25% and 50% of the total company’s operating budget, thus the overall operating cost of a company is highly dependent on the operating procedures for its compressors. The operating cost of each compressor is represented by the fuel consumed, and the network problem is formulated to minimize compressor operations that, in turn, consume fuel. The total fuel consumption of the 22 compressors is optimized by minimizing the power consumed by the compressor while satisfying specified flows and minimum pressure at each node. This can be written in mathematical form as: min z, z= ∑ HPcskTKk=1(7) where HPcsk= power required for compressor (horsepower); TK= Total number of pipes. The objective function of the problem is the sum of the fuel costs over all of the compressors in the network. In addition to fuel consumption minimization, a multi-objective optimization problem can be formulated by introducing an environmental variable that relates to the resulting environmental impact of compressor operations. Efficient compressor processes can produce less green-house gas emissions, and lead to lower environmental impacts. However, the increased cost of these compressors, and the reduction in network performance through unattractive changes in pressures may need to be traded-off against the environmental benefits of such operations. 3.1.2.2. Gas purchase minimization The minimization of gas purchase cost while satisfying all system demands and other operating constraints in a network most often has an objective function that is as follows: min z, z= ∑ Ci*SiTNi=1(8) where Ci= purchase cost of gas at the supply nodes ($/106 meters3); Si=required quantity of gas at the supply nodes (106 meters3/day) 23 3.1.2.2 Least cost gas pipe network design Least cost gas pipe network design is a minimization of the investments for the network layout. It is used in cases where the operation of the network is to be improved by either reinforcing existing pipelines or network equipment such as fittings, compressors, valves or adding new pipes to accommodate the system to a higher standard of operating tolerances. To minimize the pipe network design layout under the prescribed operating constraints described in Equations 1 through 6 in Section 3.1.1., the cost of a new pipe with a specified diameter and length and reinforcing parts. min z, z= ∑ C(Diamk)*Lk++CPF(Diamk)+CCS(Diamk)+CV(Diamk)TKk=1(9) where, C(Diamk)= specified unit cost of a pipe categorized by its diameter per linear length ($/km); Diamk= diameter of pipeline k (mm); Lk= length of pipeline k (km); CPF(Diamk)= specified unit cost of a pipe fitting in the kth pipeline with a corresponding diameter ($/mm); CCS(Diamk)= specified unit cost of a compressor station in the kth pipeline with a corresponding diameter ($/mm); CV(Diamk)= specified unit cost of a valve in the kth pipeline with a corresponding diameter ($/mm). 3.2 Heuristic algorithms for solving decision models of optimal operations of gas networks 3.2.1 Genetic Algorithm (GA) GAs, introduced by Holland (1975), are based on the mechanics of natural genetics, whereby successive generations of a population are conceived, born, and raised until they are ready to reproduce. Once new generations are “born” and older generations “die away” the algorithm codes the values of the decision variables of a problem, in the form of a genetic string, or chromosome, that represents a single population member, or solution vector. While several coding mechanisms exist, all 24 of them translate the solution vector into a series of real numbers. A large number of solutions, or strings, are randomly generated, and the natural process of a Darwinian survival of the fittest is mimicked for these strings. In GAs, the fitness of a genetic string is defined as some non-negative value that makes the objective function better. The fittest strings are allowed to reproduce, exchanging genetic material, i.e., portions of their strings, with other fit strings. Thus, the population of strings in subsequent generations is comprised of portions of the fittest members of previous generations. The algorithm works iteratively, successively generating and testing each new generation of strings. GAs use a stochastic approach for determining the breeding (or crossover) rate, as well as the likelihood, and location, of possible mutations. The search within the population is continually improving by focusing on areas of the region of solutions, or fitness landscape, with a strong likelihood of high fitness. The initial population of strings is randomly generated throughout the fitness landscape, and the algorithm proceeds to explore the entire landscape, climbing many peaks of the landscape in parallel, thus reducing the probability of being trapped in any one local optimum. The algorithm is computationally simple and powerful for undertaking searches of complex fitness landscapes. There are no restrictions on the variety of problem type to which GAs may be applied. Goldberg (1989) summarizes the difference between GAs and other search methods as follows: “GAs rely on variable sets that are coded, not on the variables themselves; Unlike many nonlinear optimization approaches that search the feasible region from a single point, GAs search the feasible region from an array of points; GAs require only the objective function (or fitness) information, not a gradient, derivative, or other auxiliary data; and Rather than using deterministic transition rules, GAs rely on probabilistic transition rules.” While many variations of GAs have been introduced, this thesis uses their simplest form, as presented in Holland (1975). The general steps of this simple GA are described in the flowchart in Figure 3. 25 Figure 3 shows the simplest GA flowchart in which an iteration begins after the initialization of population of strings. Each iteration contains reproductive operators such as selection, crossover, and mutation. Once reproductive operators are successfully implemented, this simple GA evaluates the fitness of results. This process is continued in successive iterations until the stopping criterion is met. The stopping criteria may be a specified degree of convergence of the solution or a predetermined number of generations or function evaluations. 3.2.1.1 Coding of strings GAs require the values of the decision variables of the optimization problem to be coded into strings, or solution vectors, of finite length. The most common coding method, Binary Coding, assigns a series of 0’s or 1’s to describe a variable value. Goldberg (1989) provides a simple example of Binary Coding, whereby there is only one decision variable, XGA, and the optimization problem is to maximize XGA2 subject to the bounds 0 ≤ XGA ≤ 31. The value of variable XGA is defined by the string GAS, where GAS is comprised of five strings, GAS1 through GAS5, any of which either equals 0 or 1. To code the string for XGA, Goldberg uses the following transformation of XGA: XGA=GAS1∙24+GAS2∙23+GAS3∙22+GAS4∙21+GAS5∙20 (10) At the extremes, the string (1,1,1,1,1) yields XGA = 31, and the string (0,0,0,0,0) yields XGA = 0, and therefore the transformation above can cover the whole range of values that XGA can take. Goldberg shows that by a random process, it is possible to construct multiple strings, or solution vectors, which would constitute a population like those in Table 1. Here, the first column provides the string values, the second column, the resulting value for XGA, the third column, the objective function value, or fitness value, and the fourth column, the distribution of the fitness values for this population. Because GAs work directly with the underlying code, they are not dependent upon continuity of the variable space or the existence of a derivative for the fitness function. Of course, finite string 26 coding may require more creative mappings than is shown in this example, such as those employing Mapped Fixed Point, Floating Point, and Gray codes, or a variety of hybrid coding schemes. In this thesis, simple Binary Coding is used to solve the gas transmission network operations problem. 3.2.1.2 Populations of strings The first step in implementing a GA is to set the size of the population and to randomly generate the starting population of strings. The size of the population should be sufficiently large to generate a reasonable representation of the fitness landscape, but small enough to limit the required computation time for each generation. Goldberg (1989) suggests a typical population size of 50 to 100 to provide a rich exploration of the decision space, that, while it may not guarantee fast convergence, will identify high-performance strings quickly. De Jong (1975) undertook a parametric analysis of the simple GAs and deduced that the population size should be moderate (50 to 100) to prevent genetic drift. 3.2.1.3 Reproduction Once the initial population of strings is generated, GAs perform a series of simple operations on the population to generate a new population, the next generation. These operations can be divided into two types: reproductive plans and genetic operators. Reproductive plans determine the maximum number of copies of any given string that may survive to the next generation, and genetic operators determine how the existing population of strings may be combined and modified to create the next generation of strings. For a simple GA, three genetic operators are typically used, the fitness reproduction or selection, crossover, and mutation operators. 27 3.2.1.3.1 Fitness reproduction or selection The fitness reproduction operator determines whether a given string is copied into the new generation, based on the string’s fitness. During reproduction, more highly fit strings are allowed to produce higher numbers of offspring in the mating pool. To select the strings with high fitness, Goldberg (1989) proposes that the ranking method or tournament selection be used. Under the ranking method, the probability of reproduction during a given cycle is proportional to the relative fitness of the individual string, i.e., the fitness of the string divided by the average fitness of the population. Under tournament selection, all possible pairs of strings are compared in terms of their fitness. Some of these high-fitness strings are copied into the next generation, and some of them are randomly paired for breeding, i.e., to create one or more offspring comprised of their mixed gene information. The breeding process is undertaken by the crossover operator. In this thesis, tournament selection is used as the selection operator. 3.2.1.3.2 Crossover Consider the solution vector represented as U, where: U = u1u2u3…un and ul represents the value of the lth decision variable in this solution. The crossover operator exchanges the genetic information between two selected high-fitness strings to produce one or more offspring. Crossover determines the location in the strings where the genetic material is exchanged. There are three types of crossover mechanisms, single point, multipoint, and uniform crossover. Under single point crossover, genetic material is exchanged at one randomly selected location. For example, consider two strings U and T that are to mate: U = u1u2u3u4u5 28 T = t1t2t3t4t5 Suppose the location where single point crossover is to occur is immediately after the second bit of the string. The offspring strings for this mating scenario are: U’ = u1u2t3t4t5 T’= t1t2u3u4u5 Under multi-point crossover, genetic material is exchanged at multiple randomly selected locations, and under uniform crossover, genetic information may be exchanged at an unlimited number of locations in the parent strings of genetic material and each location of which is selected at random. With the probability of 0.5, each genetic information from two first strings (parents) is assigned to two new strings (offspring). For example, consider two strings U and T from above: Before: U = u1u2u3u4u5 T = t1t2t3t4t5 After: U’ = u1t2u3u4t5 T’= t1u2t3t4u5 3.2.1.3.3 Mutation Under mutation, each offspring may be randomly changed at a random location in its string. For instance, in a string coded with Binary Coding, a random “1” could be changed to a “0.” Mutation results in somewhat of a random walk through the fitness landscape and acts as a guard against being trapped in a local optimum. While it has been shown that fitness values are not highly 29 sensitive to the rates of crossover and mutation, Goldberg (1989) recommends the use of a high crossover probability, and low mutation probability, on the order of 0.1, and 0.001, respectively. 3.2.1.4 Stopping criteria As shown in Figure 3, once the fitness is determined for the population of strings for a given generation, the stopping criteria is evaluated. In general, GAs stop when they reach a maximum number of generations or when the fitness functions for the next population of strings cannot be further improved. 3.2.2 Differential Evolution (DE) Algorithm Developed by Price and Storn (1997), DE, like GAs, is a simple population-based evolutionary algorithm. It is sufficiently powerful for optimizing nonlinear and non-differentiable continuous fitness functions and has been demonstrated to be faster and more robust than other search algorithms for solving network problems. It is designed to fulfill the following requirements (from Price and Storn, 1997): “Applicability to non-differentiable, nonlinear, and multimodal objective functions; Parallel search capabilities to cope with computationally intensive objective functions; Easily implemented, i.e., having few control variables to guide the minimization, and variables that are robust and easy to identify; and Having attractive convergence properties, i.e., repeatable and consistent convergence to the global optimum over consecutive independent trials.” The general steps of the DE are described in the flowchart in Figure 4. Figure 4 describes a similar process structure as in GAs, and provided in Figure 3, except that the order of the reproduction operators is reversed. DE initiates the reproduction with mutation, then crossover, and 30 then selection. Stopping criteria are similar to those in GAs, as the algorithm iterates until convergence of the solution or a predetermined number of evaluation trials are reached. DE undertakes a stochastic parallel direct search and operates on a set of E-dimensional solution vectors, E, comprising the population for each generation. In order for the DE algorithm to work, the initial random population must cover the entire fitness landscape. Similar to GAs, the vectors of decision variable values, or solution vectors, are coded into a finite string. 3.2.2.1 Mutation DE generates new solution vectors from the initial population differently than GAs. DE begins with the mutation operator which generates the new variable values for each decision variable by adding the weighted difference between the variable values for two randomly selection population members to a third variable value for a third randomly selected population member. For each decision variable E = 1, … , NP, this can be written in a mathematical form as: vm, g+1= Er1, g+R*(Er2, g-Er3, g) (11) Where NP is a number of members in the population; r1, r2, and r3 are random indexes (1,2, …, NP); vm, g+1 is a mutant variable of E in generation g+1, which is a combination of Er1,g, Er2,g, and Er3,g; Er1, g, Er2, g, and Er3, g represent mutually different values for variable E, for three solution vectors that are randomly selected from the initial population, and R = a real number that is greater than 0 and less than 2. 3.2.2.2 Crossover Crossover is undertaken by creating trial vectors using a specified crossover threshold, CR. Define 𝑢𝑧,g to be the E-dimensional set of decision variables for solution vector u, at generation g, where: 31 uz,g=(u1,g,u2,g,…,uE,g) (12) For each z = 1, …CO variable in a solution vector, random numbers, randb(s), are generated and compared with CR. For each of s = 1, …E variables in us,g , if randb(s) is less than or equal to CR, the variable value of the mutation is assigned to this variable for the trial vector uz,g+1. If randb(s) is greater than CR, the variable value remains unchanged in the trial vector. Mathematically, for each solution vector, z, and each generation, g: umz,g+1= {vmz,g+1 if(randb(s)≤CR) or s=rnbr(z)umz,g if(randb(s)>CR) and s≠rnbr(z) , m=1,2,…., E (13) Where 𝑢𝑚𝑧,𝑔+1 is the mth variable value for the zth trial vector for generation g+1, rnbr(z) is a randomly chosen index which ensures umz, g+1 obtains at least one parameter from vmz, g+1, and randb(s) and CR range from 0 to 1. 3.2.2.3 Selection Finally, the selection process is undertaken to select the elite vectors to survive to the next generation. If the trial vector, 𝑢𝑧,𝑔+1 , generated with mutation and crossover has a better fitness function compared with the initial solution vector, uz,g, then the trial vector survives to the next generation. If not, the initial solution vector, uz,g, survives. 3.2.2.4 Stopping criteria As shown in Figure 4, once the fitness is determined for the population of strings for a given generation, the stopping criteria is evaluated. As with the GA, DE stops when it reaches a maximum number of generations or when the fitness functions in the next population of strings cannot be further improved. 32 3.2.3 Artificial Bee Colony (ABC) Algorithm The ABC Algorithm, developed by Karaboga (2005), is a technique for solving combinatorial, multidimensional and multimodal optimization problems. The algorithm is motivated by the intelligent behavior of honey bees and is one of a number of algorithms that mimic swarm intelligence. Other examples of swarm intelligence-based optimization are ant colony, immune system, and particle swarm, i.e., the social behavior of birds flocking, optimization. ABC uses two fundamental concepts of swarm intelligence, the self-organization ability of the swarm and the division of labour within it. In nature, the main food source of bees is nectar, which is located in flowers. There are three main components involved in bee nectar foraging, the food source, employed forager bees, and unemployed forager bees, and two behavior modes: recruitment to a nectar source and abandonment of a source. The value of the food source may be dependent on how close the food is to the bees, how much food source is available and how easy it is to extract the nectar from the food source. In the ABC algorithm, there are three types of bees, employed forager, unemployed onlooker, and employed scout bees. Potential optimization solutions represent food sources and the quality of the solution represents the value of the food source or its nectar. The employed forager bees find the food sources, i.e., find potential solutions, and determine their values. The values of the food sources are then shared with the unemployed onlooker bees, with a certain probability, in what is akin to a dance in the nest. The onlooker bees wait in the nest and absorb the food source information shared by the employed forager bees. By observing the employed forager bees’ dance, they determine the value of any one food source and whether they are attracted to the food source. The longer the duration of any bee’s dance, the greater probability an onlooker bee being drawn to that food source. This exchange of information is the key process of the collective knowledge of the hive. 33 Initially, there are equal amounts of employed forager and unemployed onlooker bees and no scout bees. The nectar, or information about one source or solution, is gathered by each employed bee, and this information is shared with onlooker bees. Onlooker bees produce new solutions, evaluate and select the best among these, and determine which solutions will be explored further, based on the fitness value of each solution, relative to the total fitness of all solutions. The employed bees then go to these identified sources and continue to search in the neighborhood of the source, in adaptively reduced step lengths. If, after searching in the neighborhood of the solution for a specified number of trials, no improved solutions are identified, the given solution is abandoned and the employed bee becomes a scout bee and searches elsewhere in the solution space. The attributes of the ABC algorithm, identified by Karaboga (2005), include: “Positive feedback: As the nectar amount of the food sources increases, the number of onlookers visiting the food sources also increases. Negative feedback: The employed bees do not continue to search around poor food sources. Fluctuations: The scouts carry out a random search process for discovering new food sources. Multiple interactions: Bees share their information about food sources with their nestmates on the dance area.” 3.2.3.1 Algorithm structure The key structural elements of the ABC algorithm are: The artificial bee colony has three different types of bees; employed, onlooker, and scout bees. Initially, the number of employed forager bees is equal to the number of onlooker bees and there are no scouts. Every food source attracts only one employed bee. This means the number of food sources being examined and operated on is the same as the number of employed bees at any time. Initially, the algorithm generates a randomly distributed population of positions of food sources or solutions. Each position of the food source, posbeh is determined. Where b= 1, …, SN (number of food sources); e= 1, …, OP (number of optimization parameters); and h= 1, …, BN (number of employed bees). The initial population of the 34 positions is subjected to the repeated cycles, Cy= 1, …,Cy max, of the search processes of the employed bees, the onlooker bees, and scout bees. Onlookers are attracted to the food source with the probability proportional to the food source fitness value, i.e., higher-valued food sources attract more onlooker bees. The probability of attracting onlooker bees, Prb, associated with the food source b, is determined as follow: Prb= fb∑ fnSNn=1(14) Where 𝑓𝑏 is the fitness value of the bth solution, and SN is the number of food sources or solutions at the current iteration, which is equal to the number of employed bees (BN). After a given food source has been evaluated and shared with onlooker bees, the employed bee selects the new food source based on the most recent one selected, according to the following expression: newposbe=posbe+ϕbe(posbe-poshe) (15) Where newposbe= new food source position; posbe= the current food source position; 𝜙𝑏𝑒 =a random number between -1 and 1. The value of h cannot be same value of b; poshe= randomly chosen food source position If the value of a food source for an employed bee does not improve relative to that for the last food source found by that bee, the employed bee becomes a scout bee. The selection of a scout bee is controlled by a parameter called “limit”. The limit is the value of the number of trials for releasing a food source by the employed bee if the solution is not improved by the predetermined limit on the number of trials. The scout bees search for new food sources randomly throughout the solution space. In general, the algorithm can be described in the following steps: the employed bees locate the food source, determine the value of the food source, and share this information with onlooker 35 bees; depending on the value of the food source, onlooker bees are attracted to the food source and the employed bees are invited to search further around the attractive food source in adaptively finer step sizes; if no higher-valued food sources are identified by an employed bee, scout bees are created and these scouts randomly search for better nearby food sources. The detailed flow chart of the algorithm is shown in Figure 5. 3.2.4 Covariance Matrix Adaptation Evolution Strategy (CMA-ES) The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) an evolutionary algorithm developed by Hansen (2001) that may be used to solve difficult nonlinear, nonconvex optimization problems, in continuous domains. According to Hansen (2001), key features of the CMA-ES that allows it to overcome problems associated with other evolutionary algorithms, are: “The covariance matrix adapts the search distribution to improve performance when decision problems are poorly scaled or have highly non-separable objective functions. In the CMA-ES, the population size can be freely chosen, because the learning rates that are inherent in the covariance matrix equation prevent degeneration of the population into subspaces, even for small population sizes. The step size control feature of the CMA-ES prevents the population from converging prematurely, though the search may still resolve at a local optimum. “ Like other evolutionary algorithms, CMA-ES undergoes iterations or generations, and at each generation, solution vectors are selected and recombined to produce the next generation. Unlike other evolutionary algorithms, CMA-ES relies on a parameterized multivariate normal distribution model to generate an initial population and to spans the feasible region. In the initialization step, the mean, variance, covariance matrix and parameters such as a number of offspring, recombination weights, and several coefficients for the learning rate of the covariance matrix, are selected. Next, the initial population is generated using the initial parameter values. Then, the solution vectors are modified using a weighted average, based on the weighted intermediate recombination method, to generate the next generation of offspring solutions and the new mean of these are identified. At each subsequent 36 generation, the algorithm parameters, such as the mean, variance, and covariance matrix of the population, and the step size for trial exploitation are updated and the process continues until the stopping criteria have been met. The updated variance of the population distribution, the covariance matrix, is considered to be an important property because it limits the premature convergence of the algorithm. Over a number of search iterations, information is accumulated and this facilitates the ability to adapt the covariance matrix, even when using small populations. The adaptation of the step size allows for independent timescale control. 3.2.4.1 Initialize the parameters The initial values of the parameters include values for the offspring population size, parent population size, initial mean, variance, and step size, and several coefficients of covariance learning rates. At the first iteration, the covariance matrix is set to the identity matrix. The population of solutions consists of a set of vectors that are generated by sampling a multivariate normal distribution (N) around the mean (M) for a given generation (iteration). This can be written as follows: xa(g+1)=N (M(g),(σ(g))2,C(g)) for a =1, 2, …, ps (16) Where xa(g+1)represents the ath offspring of generation (g +1), M(g)is the mean value of the search distribution at generation g , (𝜎(𝑔)) is the standard deviation of the search distribution at generation g, C(g) is the covariance matrix of the search distribution at generation g, and ps is the population size or number of offspring, which is required to be greater than or equal to 2. Each generation requires that its parameters, which are stated above, be updated after the selection and recombination process. 37 3.2.4.2 Selection and recombination Similar to other evolutionary algorithms, the selection and recombination process occurs at each iteration. To perform selection and recombination, the new mean M(g+1) of the search distribution needs to be updated and the mean vector is then mutated to generate the offspring. First, the vectors of solutions are rank ordered in decreasing order of their objective functions, and the μ selected best individuals are identified, where μ = ½ a. The method of updating the mean uses a weighted average of these μ parents. This can be written in mathematical form as follows: M(g+1)=M(g)+Cm ∑ wτ (xτa(g+1)-M(g))µτ=1(17) With condition, of ∑ wτµτ=1=1, w1≥ w2≥… wμ>0 (18) where, wτ= a positive weight coefficient; Cm = a learning rate coefficient, which is set to 1 most of the time but is always less than or equal to 1; xτ:a(g+1)= τ th best offspring from the initial sample space; M(g)= mean of the search distribution for the current generation. The weighted vector can be either a uniform vector or a linearly varying vector. The uniform vector method equally weights all the selected samples while the linear vector method allocates more weight towards the fittest solutions than the less fit solutions. According to Hansen (2001), the individual weights may be found using implementation of weighted intermediate recombination and denotes μeff , as the variance effective selection mass, where 1< μeff< μ, and is expressed as: μeff=1∑ wτ2µτ=1(19) 38 3.2.4.3 Step size control The step size needs to be updated in each iteration because the covariance matrix can only increase or decrease the scale in one direction for each evaluation. Hansen (2001) describes the justification for the step size control mechanism, as a complement to the covariance matrix, as: “1. The optimal overall step length cannot be well approximated by the covariance matrix, in particular, if μeff is chosen to be larger than the value of one; and 2. The largest reliable learning rate for the covariance matrix update is too slow to achieve competitive change rates for the overall step length.” Thus, as proposed by Hansen (2001), the step size is independent of the covariance matrix. To update the step size, cumulative step length adaption (CSA) is used. σ(g+1)= σg × exp (Cσdσ(‖Pσ‖N(0,I)-1)) (20) Where Cσ= learning rate coefficient for the accumulation of the step size control and is < 1; dσ= damping parameter for the step size updates; N(0,I) = expected value of a multivariate normal distribution with zero mean and unity covariance matrix; 𝑃𝜎 = length of an evolution path that can be described as a sum of consecutive steps. The formulation of an evolution path can be written as follows: Pc(g+1)=(1-Cc)Pc(g)+√Cc(2-Cc)μeff M(g+1)-M(g)σ(g)(21) Where, Cc = learning rate coefficient for the accumulation of the rank-one update of the covariance matrix; Pc = evolutionary path for covariance matrix and step size. According to Hansen (2001), the length of an evolution path is exploited due to the following reasons: “ Whenever the evolution path is short, single steps cancel each other out, and the steps are anti-correlated. Thus, step size must be reduced so that the steps do not cancel each other out. 39 Whenever the evolution path is long, the single steps point in similar directions, and the steps are correlated. Thus, step size must be increased so that the feasible region can be covered by fewer, but longer, steps in the same direction. Steps are approximately perpendicular in expectation, and thus are uncorrelated. “ 3.2.4.4 Covariance matrix adaption The covariance matrix determines the second order model of the objective function of the distribution. The process uses samples from the previous iteration with a new mean at generation g+1 and a covariance matrix of the model distribution multiplied by the step size around the mean. The new generation covariance matrix is constructed using a combination of two rank update methods and they are the rank -µ-update and rank-one-update. These two rank updates provide a reliable estimation of a new generation of the covariance matrix whether the previous generation of covariance matrix was driven from a single solution or a small population set. Rank- µ - update uses the entire population and rank-one update uses the evolution path to provide the information of correlations between generations. The covariance matrix describes the dependencies between the variables in the distribution. For a multivariate normal distribution, the covariance matrix is a positive definite matrix. The following mathematical expression describes how to determine the covariance matrix of the new generation. C(g+1)= (1-Cone-Cμ ∑ wy) C(g)+ConePc(g+1)Pc(g+1)T+Cμ ∑ wτ(x1:a(g+1)-M(g)σ(g)λτ=1) (x1:a(g+1)-M(g)σ(g))T(22) where, Cone= learning rate coefficient for the rank-one update of the covariance matrix; Cµ= learning rate coefficient for the rank-µ update of the covariance matrix, and 1- Cone = Cµ; and ∑wy = similar weighting value as ∑ wτλτ=1 . 40 3.2.4.5 Stopping criteria The execution of the algorithm is iterative and continues until one of three stopping criteria is met. CMA-ES can run for a predetermined number of iterations, it can stop when the algorithm does not make any improvement in the optimal solution, or it can stop at a predetermined objective function value. The detailed flow chart of the algorithm is shown in Figure 6. 3.3 Hydraulic simulators for assessing the performance of gas networks In order to verify the solution of an optimization model for pipe network operations, one generally checks the solutions with a hydraulic simulator of the distribution of pressure and flows throughout the network due to continuity and energy considerations. In water distribution analyses, such simulators are also used in post-optimization analyses, such as cut-set analyses, which indicate the robustness of the network design in the face of mechanical failures. In 1987, Todini and Pilati introduced the Global Gradient Algorithm (GGA) which is a hybrid simulator that uses optimization to implement the Newton Raphson Method for minimizing errors in predicting pressure and flows values at network nodes. The GGA is designed for general looped water distribution network problems, however, a partial formulation can be used for branch networks. According to Todini and Pilati (1987), the general formulation of the GGA is as follows: {Hg+1=(A21DIAG-1A12)-1*(A21Qg+Du-A21DIAG-1A11Qg-A21DIAG-1A10H0)Qg+1=Qg-DIAG-1 (A11Qg+A12Hg+1+A10H0)(23) Where, A10= a matrix representing the relationships between source nodes and pipes linkage; 41 A11= a diagonal matrix representing the relationships between nodes containing equipment including fittings, compressors, valves, and pumps and specified pipelines; A12= a matrix of nodes and pipes linkage (the topology of the network) where, {-1=if the flow of the pipe leaves the node0=if the pipe is not connected to the node+1=if the flow of the pipe enters the node; A21= transpose of A12; DIAG-1 = inverse of a diagonal matrix of flow equation in pipe or pump; Q= unknown pipe flows at each node; Qg= gth iteration of the flow matrix; Qg+1= g th +1 next iteration of the flow matrix; H= unknown pressure head at a given node; H0= known pressure head at a source node; Hg+1= next iteration pressure head matrix; and Du= known demand flow at each node The GGA requires at least one node with a known pressure head in order to solve the hydraulic simulation using the Newton Raphson Algorithm. The GGA undergoes an iterative process to update all matrices as it moves toward convergence of the Newton Raphson Method, while an iteration process is undertaken, the entire system of equations updates pipe flows and the node pressures. 42 Figure 2. Network Schematic 43 Figure 3. Simple GA Flowchart (Created based on Goldberg (1989)) 44 Figure 4. Differential Evolution Flowchart (Created based on Price and Storn (1997)) 45 Figure 5. Artificial Bee Colony Algorithm Flowchart (Created based on Karaboga (2005)) 46 Figure 6. CMA-ES Flowchart (Created based on Hansen (2001)) 47 Table 1. Binary Code Example of GA String Associated XGA Fitness value (XGA2) % of total (0,0,1,1,1) (0,1,1,0,1) (1,0,1,0,1) (0,1,0,0,1) 7 13 21 9 49 169 441 81 6.6 22.8 59.6 11 740 100 48 Chapter 4: Case Study: Belgium Gas Network The performance of the four evolutionary algorithms, for identifying the optimal operation strategy for gas networks, is demonstrated in this Chapter for a case study based on the primary pipeline elements for delivering high calorific gas in the Belgium Gas Transmission Network. Figure 7 provides a simplified schematic of this network, where the types of the high- and low-calorific gas supplied by a given branch of the network, and the location of compressors and storage reservoirs, are identified. The case study was developed by DeWolf and Smeers (2000) for analyzing the application of the extended Simplex Algorithm, a classical optimization technique, to the gas network operations problem. As shown in Figure 7, Belgium has no domestic gas resources and relies solely on importing natural gas from several countries including: the Netherlands, Norway, and Algeria. Dutch gas enters the system in Belgium at Poppel, Algerian gas enters at Zeebrugge (Node 1), and Norwegian gas enters via the Netherlands at Voeren (Node 8). The gas is transmitted through the network and delivered to France and Luxemburg, exiting the system at Blaregnies (Node 16) and Pétange (Node 20), in Belgium, respectively. The network transmits two types of gas: high-calorific and low-calorific gas. The high-calorific gas contains10,000 kilocalories per cubic meter and is sourced from Algeria and Norway. The low-calorific gas contains 8,400 kilocalories per cubic meter and is sourced from the Netherlands. For simplicity, this study focuses on the former gas type, because it is the major source of gas for the region, and requires compression at several nodes. Assuming that the branches for the two types are separated and thus operated in parallel, and that substitutions between gas types are not allowed, this application could be generalized and extended to analyze the transmission of both types of gas with little additional reformulation effort. 49 The following tables describe the demand for high-calorific gas at the nodes in the network (Table 2), the allowed limits for flow and pressure at the nodes (Table 3), and the pipeline characteristics of Belgium Gas Network (Table 4), as provided by Dewolf and Smeers (2000). In Table 3, each node is assigned to each location of selected Belgium cities. Each city has a minimum allowable net gas to be supplied (Lower Limit of Supply) and a maximum allowable net gas to be supplied (Upper Limit of Supply). A negative value of supply means ‘demand’ and positive supply means ‘supply’. If the lower limit of flow is 0 then there is no demand in that node and if the lower and upper limit of flow is 0 then the node is only used for bypassing. Negative infinity means there is no minimum allowable net gas to be supplied. Purchasing cost of gas is listed in Table 3. Algerian gas and Norwegian gas are purchased from Node 1 (Zeebrugge) and 8 (Voeren), respectively, and is distributed through the pipe transmission line. Some of the gas is stored in a storage facility and the gas can be taken from the storage facility when it is needed. As shown in Figure 7, there are 4 storage facilities for gas in the Belgium gas network and they are located at Nodes 2 (Dudzele), 5 (Loenhout), 13 (Anderlues), and 14 (Péronnes-lez-Binche). Storage nodes have an implicit price when the gas is drawn from the storage. Withdrawals from Node 2 (Dudzele) and 5 (Loenhout) are priced based on Algerian prices and Node 13 (Anderlues) and 14 (Péronnes-lez-Binche) are priced based on Norwegian prices. The Fk value as defined by Equation 4 in Chapter 3, is the coefficient that describes pipe characteristics, and is calculated based on formula developed by Stoner (1969) which it is derived from the Waymouth (1912) equation as: Fk=96.074 830 10-15 Diamk5ZTLkδ*[2 log (3.7Diam𝑘ε) ]2(24) where Lk = length of the pipe (km); Diamk = Diameter of the pipe (mm); Tgas= 281.15, gas temperature (°K); δ = density of the gas relative to air at Standard Temperature and Pressure (IUPAC, 50 1982) (=0.6106, dimensionless); Z= gas compressibility factor (0.8, dimensionless); and ε= wall roughness of a carbon steel pipe (=0.05 mm). As summarized in Tables 2-4, the network consists of 20 nodes, 24 pipes where five of these are paired pipes in parallel, and 3 compressor stations. Each pipe has a different Fk coefficient value due to its length and diameter, though the material of the pipes is constant. Each node has its own upper and lower limits for pressure and flow, regardless of whether they are demand or supply nodes. 4.1 Belgium gas network operations optimization model and evolutionary algorithm implementation The objective function of the model to optimize the operation of the Belgium Gas Network, as proposed by DeWolf and Smeers (2000), is a gas purchase minimization model, with the objective function of Equation 8 from Chapter 3. The constraints are continuity, supply quantity daily limits, pressure limits, and energy equations which are described by Equations 1-4, and 6 from Chapter 3. Note that the constraint in Equation 6 is formulated through the use of penalty functions for the GA and DE Algorithms. The parameters for implementing the evolutionary algorithms defined in Equation 11 through 22 are provided in Table 5. Each algorithm is implemented to solve for Equations 1-4, 6 and 8 in ten trials. Results of the model solutions are obtained for two stopping criteria: one for which 20,000 function evaluations are undertaken, and one for which 50,000 function evaluations are undertaken. All computations are implemented in the MATLAB programming language environment with an Intel® Core™ 2Duo CPU P8700 @ 2.53 GHz and 4.00 GB RAM. 51 4.2 Optimal operations of the Belgium gas transmission network Figure 8 shows the box-whisker plots of the optimal total cost of gas supplied during operation of the Belgium Gas Transmission Network for delivery of high-calorific gas, for the ten trials of algorithm implementation, with stopping criteria of 20,000 (Figure 8a) and 50,000 (Figure 8b) function evaluations. The solid bars represent the range of the ten minimum cost values between the 25th and 75th percentile values, the values marked with an “x” and a “-” correspond to the mean and median minimum cost values, respectively. The extreme values of the plots are the highest and lowest minimum cost determined. The optimal cost determined by DeWolf and Smeers (2000) is also shown by the green dashed-line. The results for the CMA-ES for the stopping criterion of 20,000 function evaluations are not included because the CMA-ES did not converge for this number of evaluations. For the stopping criterion of 50,000 function evaluations, the CMA-ES converges to the minimum cost obtained by DeWolf and Smeers (2000) for every one of the ten trials. The lowest of the minimum cost values for gas for the Belgian Gas Network that is determined using all four evolutionary algorithms is 91.06 million $. This value is obtained in applications of the GA, and the DE, ABC, and CMA-ES Algorithm, under the 20,000 and 50,000 function evaluation stopping criterion, respectively. The value of the minimal gas purchase cost is same as that obtained by Dewolf and Smeers (2000) using the extended simplex algorithm. The mean minimum cost for each algorithm varies, as shown in Figure 8. For the 20,000 function evaluations stopping criterion, the DE produces the lowest mean minimum cost value of 91.48 million $, and the narrowest range of minimum cost values among the algorithms. The ABC and GA algorithms have the second and the third smallest mean minimum cost values of 92.95 million $ and 93.32 million $, respectively, for these stopping criteria. For the 50,000 function evaluations stopping criteria, the CMA-ES had the best performance among the four algorithms with the mean minimum cost value of 91.06 million $ and a minimal variation among all ten trails. The mean value for the DE algorithm, in this case, increases slightly compared with the 20,000 function evaluations stopping criterion, to a 52 value of 91.80 million $, with a corresponding minimum value of 91.48 million $. The reason that the DE performance decreases is because of a penalty function associated with constraints in Equation 6, the pressure limit. Unlike to DE, the performance of the GA and ABC algorithms improves as the number of function evaluations is increased to 50,000. The mean minimum value of the ABC and GA algorithms again have the second and the third smallest mean minimum cost values of 91.80 million $ and 92.20 million $, respectively. While the CMA-ES performs well in terms of cost, the 2nd and 7th trials of CMA-ES results in pressure violations at the Node 8, Voeren, and Node 17, Wanze, as shown in Table 6. These pressure violations the CMA-ES are relatively small and the percent violation is generally acceptable. The solutions found with the ABC Algorithm does not violate any pressure or flow requirements for both stopping criteria, while those found with the GA and DE Algorithms violate pressure requirements for both stopping criteria. As shown in Table 7, the GA has the worst performance among the algorithms in terms of pressure violations. Among the ten trials, 15 nodes experience pressure violations, though these violations are not at the same nodes for all ten trials. Among these violations, for the upper pressure limit and the lower pressure limit, the maximum violation relative to the limit is 64.13% and 13.39%, respectively. As shown in Table 8, among the ten trials of the implementation of the DE algorithm, the pressure limits at five nodes are violated, though these violations are not at the same nodes for all ten trials. Among these violations, for the upper-pressure limit and the lower pressure limit, the maximum violation relative to the limit is 2.87% and 0.25%, respectively. Note that for Node 8, Voeren, and 9 (Berneau) for the 20,000 function evaluations stopping criterion, all trials of the DE result in the upper-pressure limit being violated and for the 50,000 function evaluations stopping criterion, the pressure limits at these nodes are violated for nine trials. In general, Tables 6-8 shows that the algorithms have difficulty converging to solutions that meet the pressure requirements at Nodes 8 (Voeren), 9 (Berneau), 16 (Blaregnies), 17 (Wanze), and 53 20 (Pétange). Nodes 8 (Voeren), 9 (Berneau), and 17 (Wanze) are particularly problematic, whereas the solutions found with the evolutionary algorithms are somewhat more robust for meeting the pressure targets for Nodes 16 (Blaregnies) and 20 (Pétange). Figures 9-28 are the box-whisker plots for the resulting pressures (bar) at Nodes 1 to 20, respectively. As shown in these figures, there are some similarities among nodes in the network. For example, for each algorithm, the pressure response of Nodes 1 (Zeebrugge), 2 (Dudzele), 3 (Bruges), 6 (Antwerp) and 7 (Ghent) are nearly identical for all four algorithms. Node 4 (Zomergem) is only slightly different from these nodes as well, though the resulting pressure decreases because the lower pressure limit is relaxed. Node 5 (Loenhout) has the most freedom in terms of the pressure range, and as shown in the plot has more variance because it is also a storage node. For each algorithm, the resultant pressures for Nodes 8-16, are very similar, where the pressure value may decrease as the nodes change, but the range of resulting pressures remains the same. For Nodes 17 (Wanze) and 18 (Sinsin), the resulting pressures under application of each algorithm tend to stay near the upper pressure limit. At Node 19 (Arlon) pressure is roughly midway between two pressure limits and at Node 20 (Pétange) the pressure is near to the lower pressure limit. For Nodes 8-20, DE, ABC, and CMA-ES have a similar range of variation for both stopping criteria, while for Nodes 8-17, the performance range for the application of GA improves when the stopping criterion is changed from 20,000 to 50,000 function evaluations. There is no improvement in the performance of GA across stopping criteria for Nodes 18-20. The DE algorithm decreases in performance across stopping criteria at Nodes 1-4, while all other algorithms improve. If the pressure results are compared with those found by DeWolf and Smeer (2000), the ABC Algorithm performs most closely to their results, followed by the CMA-ES, DE and GA Algorithms, in that order. CMA-ES shows a strong performance for Nodes 8 to 20 with very small variation in pressure levels, while the ABC Algorithm ranks second in terms of the consistency of the variance in pressures throughout all nodes, followed by the DE and GA Algorithm, in that order. 54 Figures 29-43 are the box-whiskers plots for the available flow to meet demands, or to be supplied to, Nodes 1-20. In terms of providing flows for meeting demands and supply contracts, all algorithms manage to satisfy flow restrictions, and, similar to the algorithm employed by DeWolf and Smeer (2000), result in flow levels that are close to the upper flow limit. Compressor power usage, in kW, was also calculated and represented as box-whisker plots. Figures 44-46 are the plots of the required compressor power (in kW) at the three compressor locations, Pipe 9, 19, and 22. Based on a comparison of these figures, the compressor located in Pipe 9 consumes a wide range of power in contrast to the other compressors, in order to compensate the pressure loss and meet the pressure requirement at Node 14 (Péronnes-lez-Binche). The CMA-ES Algorithm results in solutions that utilize significantly more power for this compressor compared with the compressors in Pipes 19 and 22. Among the four algorithms, the GA solution results in the widest range of power variation for each compressor. The DE Algorithm results in solutions that utilize a minimum range of compressor power for all three compressors. Perhaps, the reason that the DE Algorithm could not meet the pressure demands as well as other algorithms is due to the fact that its solutions results a minimum usage of compressors. In general, the application of the ABC and CMA-ES Algorithms result in solutions that meet pressure requirements and demand and supply restrictions, and conserve compressor power. Overall, the best performance is exhibited by the CMA-ES Algorithm, followed, in order of effectiveness, by the ABC, DE and GA Algorithms, respectively.55 Figure 7. Belgium Gas Transmission Network (Drawn on Google Earth Image and modified from DeWolf and Smeers (2000)) 56 Figure 8. The Total Cost, million $, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 57 Figure 9. Pressure at Node 1 (Zeebrugge), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 58 Figure 10. Pressure at Node 2, (Dudzele), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 59 Figure 11. Pressure at Node 3 (Bruges), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 60 Figure 12. Pressure at Node 4 (Zomergem), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 61 Figure 13. Pressure at Node 5 (Loenhout), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 62 Figure 14. Pressure at Node 6 (Antwerp), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 63 Figure 15. Pressure at Node 7 (Ghent), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 64 Figure 16. Pressure at Node 8 (Voeren), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 65 Figure 17. Pressure at Node 9 (Berneau), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 66 Figure 18. Pressure at Node 10 (Liège), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 67 Figure 19. Pressure at Node 11 (Warnant), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 68 Figure 20. Pressure at Node 12 (Namur), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 69 Figure 21. Pressure at Node 13 (Anderlues), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 70 Figure 22. Pressure at Node 14 (Péronnes-lez-Binche), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 71 Figure 23. Pressure at Node 15 (Mons), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 72 Figure 24. Pressure at Node 16 (Blaregnies), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 73 Figure 25. Pressure at Node 17 (Wanze), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 74 Figure 26. Pressure at Node 18 (Sinsin), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 75 Figure 27. Pressure at Node 19 (Arlon), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 76 Figure 28. Pressure at Node 20 (Pétange), bar, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 77 Figure 29. Supply at Node 1 (Zeebrugge), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 78 Figure 30. Supply at Node 2 (Dudzele), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 79 Figure 31. Demand at Node 3 (Bruges), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 80 Figure 32. Demand at Node 5 (Loenhout), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 81 Figure 33. Demand at Node 6 (Antwerp), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 82 Figure 34. Demand at Node 7 (Ghent), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 83 Figure 35. Supply at Node 8 (Voeren), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 84 Figure 36. Demand at Node 10 (Liège), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 85 Figure 37. Demand at Node 12 (Warnant), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 86 Figure 38. Supply at Node 13 (Namur), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 87 Figure 39. Supply at Node 14 (Péronnes-lez-Binche), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 88 Figure 40. Demand at Node 15 (Mons), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 89 Figure 41. Demand at Node 16 (Blaregnies), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 90 Figure 42. Demand at Node 19 (Arlon), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 91 Figure 43. Demand at Node 20 (Pétange), 106 m3/day, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 92 Figure 44. Power for Compressor at Pipe 9, kW, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 93 Figure 45. Power for Compressor at Pipe 19, kW, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations 94 Figure 46. Power for Compressor at Pipe 22, kW, for the Stopping Criteria of a) 20,000 and b) 50,000 Function Evaluations95 Table 2. Belgian Daily Maximum High-Calorific Gas Demand by Location in 1989 Node Location Demand (106 m3/day) 3 Brugge 3.918 6 Antwerp 4.034 7 Ghent 5.256 10 Liège 6.365 12 Namur 2.12 15 Mons 6.848 16 Blaregnies (To France) 15.616 19 Arlon 0.222 20 Pétange (To Luxemburg) 1.919 Total 46.298 96 Table 3. Allowable Limits of Flow and Pressure at the Nodes Node Location Lower Limit of Supply (106 m3/day) Upper Limit of Supply (106 m3/day) Lower Limit of Pressure (bar) Upper Limit of Pressure (bar) Purchasing Cost of Gas ($/MBTU) 1 Zeebrugge 8.870 11.594 0 77 2.28 (Algeria) 2 Dudzele 0 8.4 0 77 2.28 (Storage) 3 Bruges -∞ -3.918 30 80 0 4 Zomergem 0 0 0 80 0 5 Loenhout 0 4.8 0 77 2.28 (Storage) 6 Antwerp -∞ -4.034 30 80 0 7 Ghent -∞ -5.256 30 80 0 8 Voeren 20.344 22.012 50 66.2 1.68 (Norway) 9 Berneau 0 0 0 66.2 0 10 Liège -∞ -6.365 30 66.2 0 11 Warnant 0 0 0 66.2 0 12 Namur -∞ -2.120 0 66.2 0 13 Anderlues 0 1.2 0 66.2 1.68 (Storage) 14 Péronnes-lez-Binche 0 0.96 0 66.2 1.68 (Storage) 15 Mons -∞ -6.848 0 66.2 0 16 Blaregnies -∞ -15.616 50 66.2 0 17 Wanze 0 0 0 66.2 0 18 Sinsin 0 0 0 63 0 19 Arlon -∞ -0.222 0 66.2 0 20 Pétange -∞ -1.919 25 66.2 0 97 Table 4. Pipe Layout Description Pipe From To Diameter (mm) Length (km) Fk2 Value Compressor Present? 1 Zeebrugge Dudzele 890 4 9.07027 No 2 Zeebrugge Dudzele 890 4 9.07027 No 3 Dudzele Bruges 890 6 6.04685 No 4 Dudzele Bruges 890 6 6.04685 No 5 Bruges Zomergem 890 26 1.39543 No 6 Loenhour Antwerp 890 43 0.100256 No 7 Antwerp Ghent 590.1 29 0.148655 No 8 Ghent Zomergem 590.1 19 0.226895 No 9 Zomergem Péronnes-lez-Binche 590.1 55 0.659656 Yes 10 Voeren Berneau 590.1 5 7.25622 No 11 Voeren Berneau 890 5 0.108033 No 12 Berneau Liège 890 20 1.81405 No 13 Berneau Liège 395.5 20 0.0270084 No 14 Liège Warnant 890 25 1.45124 No 15 Liège Warnant 395.5 25 0.0216067 No 16 Warnard Namur 890 42 0.863836 No 17 Namur Anderlues 890 40 0.907027 No 18 Anderlues Péronnes-lez-Binche 890 5 7.25622 No 19 Péronnes-lez-Binche Mons 890 10 3.62811 Yes 20 Mons Blaregnies 890 25 1.45124 No 21 Warnant Wanze 395.5 10.5 0.0514445 No 22 Wanze Sinsin 315.5 26 0.00641977 Yes 23 Sinsin Arlon 315.5 98 0.00170320 No 24 Arlon Pétange 315.5 6 0.0278190 No 98 Table 5. Parameter Values for Evolutionary Algorithms Genetic Algorithm (GA) Differential Evolution (DE) Artificial Bee Colony (ABC) Covariance Matrix Adaptation Evolution Strategy (CMA-ES) Population size 100 Population size 100 Bee colony size (employed bees and onlooker bees) 100 Step-size 0.5 Mutation rate 0.01 Lower bound of scaling factor 0.2 Ratio of employed bees /onlooker bees 50/50 Population size 100 Crossover rate 0.8 Upper bound of scaling factor 0.8 Number of food sources 50 Number of parents 50 Crossover probability 0.3 Trial limit for abandoning a food source given no improvement 100 99 Table 6. Pressure Violations Resulting from Application of CMA-ES for the Stopping Criterion of 50,000 Function Evaluations Algorithm Node Function evaluations Maximum violation Mean violation of each function evaluation CMA-ES 8, Voeren 50000 0.004% (relative to upper limit) 0.004% 17, Wanze 50000 0.002% (relative to upper limit) 0.001% 100 Table 7. Pressure Violations Resulting from Application of GA for the Stopping Criteria of the 20,000 and 50,000 Function Evaluations Algorithm Node Function evaluations Maximum violation Mean violation of each function evaluation GA 6, Antwerp 20000 0.22% (relative to lower limit) 0.22% 50000 0.25% (relative to lower limit) 0.25% 7, Ghent 20000 0.39% (relative to lower limit) 0.39% 50000 0.29% (relative to lower limit) 0.29% 8, Voeren 20000 36.94% (relative to upper limit) 23.74% 50000 6.32% (relative to upper limit) 4.52% 9, Berneau 20000 36.50% (relative to upper limit) 23.24% 50000 5.75% (relative to upper limit) 3.94% 10, Liège 20000 34.72% (relative to upper limit) 28.96% 50000 3.43% (relative to upper limit) 3.43% 11, Warnant 20000 33.60% (relative to upper limit) 27.78% 50000 1.95% (relative to upper limit) 1.95% 12, Namur 20000 31.85% (relative to upper limit) 25.92% 13, Anderlues 20000 30.67% (relative to upper limit) 24.67% 14, Péronnes-lez-Binche 20000 30.49% (relative to upper limit) 24.48% 15, Mons 20000 29.26% (relative to upper limit) 23.19% 16, Blaregnies 20000 27.77% (relative to both limits) 13.14% 50000 7.67% (relative to lower limit) 3.84% 17, Wanze 20000 61.63% (relative to upper limit) 28.75% 50000 23.75% (relative to upper limit) 11.96% 18, Sinsin 20000 64.13% (relative to upper limit) 32.81% 50000 22.93% (relative to upper limit) 9.60% 19, Arlon 20000 33.76% (relative to upper limit) 16.63% 20, Pétange 20000 32.58% (relative to upper limit) 15.27% 101 Table 8. Pressure Violations Resulting from Application of DE for the Stopping Criteria of the 20,000 and 50,000 Function Evaluations Algorithm Node Function evaluations Maximum Violation Mean violation of each function evaluation DE 8, Voeren 20000 2.87% (relative to upper limit) 2.07% 50000 2.68% (relative to upper limit) 2.14% 9, Berneau 20000 2.29% (relative to upper limit) 1.50% 50000 2.13% (relative to upper limit) 1.60% 16, Blaregnies 50000 0.25% (relative to lower limit) 0.25% 17, Wanze 50000 1.49% (relative to upper limit) 1.49% 20, Pétange 20000 0.07% (relative to lower limit) 0.07% 50000 0.06% (relative to lower limit) 0.05% 102 Chapter 5: Conclusions and Future Work The thesis presents an analysis of the use of evolutionary algorithms for determining the optimal operating strategy of gas transmission networks that minimize the total operating cost of the gas network and identify the optimal demand and supply flows. In comparison with water distribution networks, gas transmission networks, are complex due to their use of compressor equipment and the large lengths and diameters of their pipelines. The optimal operation strategies for pre-existing gas networks are difficult to determine with classic optimization algorithms, given the limitations associated with their existing hydraulic conditions and compressor characteristics. Four evolutionary algorithms, the GA, DE, ABC, and CMA-ES Algorithms are evaluated for their performance in determining the optimal operation strategy of the Belgium Gas Network. The hydraulics of the system are analyzed with the Global Gradient Algorithm. Among the four algorithms, CMA-ES is found to achieve the global optimum, compared with that generated with classical optimization, most consistently. CMA-ES performs the best among the four algorithms because CMA-ES results in solutions that meet pressure requirements and demand and supply restrictions, and conserve compressor power. The ABC Algorithm performance ranks second-best with good consistency and an ability to satisfy all constraints on optimal operation such as continuity, energy equation, demand and supply limits, and compressor constraints. The DE and GA Algorithms exhibit the least reliable results for the given case study due to constraint violations. The DE Algorithm results in solutions that utilize a minimum range of compressor power for all three compressors, however, because of the minimal usage of the compressors, the DE Algorithm cannot meet the pressure targets at several nodes. The GA Algorithm results in the widest 103 range of power for each compressor and fails to meet the pressure targets at many nodes. The GA Algorithm could not utilize the compressors efficiently to meet all pressure targets. Compressors play a key role in the gas network for meeting pressure requirements. When the CMA-ES is applied, the compressors are placed in the network in a way that boosts the pressure of the gas to meet pressure targets but also at the same time, the compressors run efficiently to not consume much gas for their own use. Pressure constraints appear to be problematic since the GA and DE Algorithms have trouble meeting these targets even though they activate use of the compressors. Water distribution networks have similar aspects as gas networks, such as pipes that haul resources, booster stations to boost pressure of the system, and the layout of the network. Given this similarity, the optimization model of the gas network can be used for water distribution networks after adjustments are made to model parameters. Algorithms that are used in this thesis, have also been used for water distribution systems and perform analogously in water distribution problems. This thesis and its findings for gas networks will benefit the gas industry, water industry, and other industries facing network problems. By determining the optimal operation of a given network, the industries can save money in operations. With the use of CMA-ES to determine the optimal operating pressures and demand and supply flows, the industries can easily analyze the network and still achieve reliable optimal solutions. Further work, to advance the findings of this thesis, should include a sensitivity analysis of the demand and hydraulic conditions, and optimization of operations of the entire Belgium Gas Network including the low-calorific gas sources and demand sites, the addition of other machinery, such as valves, and analysis of the most current layout of the network. The sensitivity analysis may be undertaken to determine which node most significantly influences the optimal operating solution, and how sensitive is the optimal solution to variations in demands and pressure requirements at the nodes. The Belgium Gas Network analyzed in this thesis is in 1989 layout. The network has grown 104 and become more complicated since then. The most recent layout of the Belgium Gas Network contains more machinery, and therefore should be investigated in future work. An uncertainty analysis of the effect of uncertainty of the demand of gas users on operations in the network should also be conducted. The uncertain parameters can be described with known statistical distributions or with random distributions, depending on the availability of demand information. A reliability analysis which examines both of two failure modes, performance failure and structural failure, may also provide insights into the robustness of the optimal operating solution. 105 References An, S. (2004). Natural Gas and Electricity Optimal Power Flow. Doctoral dissertation, Oklahoma State University, Stillwater, OK. Andre, J., Bonnans, F., & Cornibert, L. (2009). Optimization of capacity expansion planning for gas transportation networks. European Journal of Operational Research, 197(3), 1019-1027. doi:10.1016/j.ejor.2007.12.045 Babonneau, F., Nesterov, Y., & Vial, J.-P. (2012). Design and Operations of Gas Transmission Networks. Operations Research, 60(1), 34-47. doi:10.1287/opre.1110.1001 Babu, B., Angira, R., Chakole, P. G., & Syed Mubeen, J. (2003). Optimal design of gas transmission network using differential evolution. Department of Chemical Engineering, Birla Institute of Technology & Science, Pilani-333031 (Rajasthan) India. Babu, B., Chakole, P. G., & Syed Mubeen, J. (2005). Differential evolution strategy for optimal design of gas transmission network. Multidiscipline Modeling in Materials and Structures, 1(4), 315-328. Bi, W., Dandy, G. C., & Maier, H. R. (2015). Improved genetic algorithm optimization of water distribution system design by incorporating domain knowledge. Environmental Modelling & Software, 69, 370-381. Bargmann, D., Ebbers, M., Heinecke, N., Koch, T., Kühl, V., Pelzer, A.,Spreckelsen, & K. (2015). State-of-the-art in evaluating gas network capacities. In T., Koch Editor (eds.), Evaluating Gas Network Capacities. Society for Industrial and Applied Mathematics. doi:10.1137/1.9781611973693.ch4, 65-84. Berardi, L., Giustolisi, O., & Todini, E. (2009). Enhanced global gradient algorithm: A general formulation. In World Environmental and Water Resources Congress 2009: Great Rivers, Kansas city, Missouri, USA. 17-21 May (pp. 1-10). Borraz-Sánchez, C., & Haugland, D. (2011). Minimizing fuel cost in gas transmission networks by dynamic programming and adaptive discretization. Computers & Industrial Engineering, 61(2), 364-372. doi:10.1016/j.cie.2010.07.012 Borraz-Sanchez, C., & Rios-Mercado, R. Z. (2004). A non-sequential dynamic programming approach for naturalgas network optimization. WSEAS Transactions on Systems, 3(4), 1384-1389. Borraz-Sánchez, C., & Ríos-Mercado, R. Z. (2009). Improving the operation of pipeline systems on cyclic structures by tabu search. Computers & Chemical Engineering, 33(1), 58-64. doi:10.1016/j.compchemeng.2008.07.009 Botros, K., Sennhauser, D., Jungowski, K., Poissant, G., Golshan, H., & Stoffregen, J. (2004). Multi-objective optimization of large pipeline networks using genetic algorithm. Paper presented at the 2004 International Pipeline Conference, Calgary, Alberta, Canada, 4-8 October, 2004 (pp. 2005-2015). 106 Brkić, D. (2009). An improvement of Hardy Cross method applied on looped spatial natural gas distribution networks. Applied Energy, 86(7-8), 1290-1300. doi:10.1016/j.apenergy.2008.10.005 Chebouba, A. (2015). Multi objective optimization of line pack management of gas pipeline system. Journal of Physics: Conference Series, 574, 012114. doi:10.1088/1742-6596/574/1/012114 Chebouba, A., Yalaoui, F., Smati, A., Amodeo, L., Younsi, K., & Tairi, A. (2009). Optimization of natural gas pipeline transportation using ant colony optimization. Computers & Operations Research, 36(6), 1916-1923. doi:10.1016/j.cor.2008.06.005 Cheesman, A. P. (1971). Understanding origin of pressure is a key to better well planning. Oil and Gas Journal, 69(46), 146. Coelho, P. M., & Pinho, C. (2007). Considerations about equations for steady state flow in natural gas pipelines. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 29(3), 262-273. Cornish, R. J. (1939). The analysis of flow in networks of pipes.(Includes plates and indeces). Journal of the Institution of Civil Engineers, 13(2), 147-154. Correa-Posada, C., & Sánchez-Martín, P. (2014). Gas network optimization: A comparison of piecewise linear models. Optimization Online. Available: http://www.optimization-online.org/DB_FILE/2014/10/4580.pdf. Cross, H. (1936). Analysis of Flow in Networks of Conduits or Conductors. University of Illinois at Urbana Champaign, College of Engineering. Engineering Experiment Station. Bulletin No.286, 34(22). De Jong, K. A. (1975). An Analysis of the Behavior of a Class of Genetic Adaptive Systems. Doctoral dissertation, University of Michigan, Ann Arbot, MI. (Order No. 7609381). Available from ProQuest Dissertations & Theses Global. (302788400). de Mélo Duarte, H., Goldbarg, E., & Goldbarg, M.(2006). A tabu search algorithm for optimization of gas distribution Networks. In: Gottlieb J., Raidl G.R. (eds) Evolutionary Computation in Combinatorial Optimization. EvoCOP, Springer, Berlin, Heidelberg, Vol 3906. De Wolf, D., & Smeers, Y. (1996). Optimal dimensioning of pipe networks with application to gas transmission networks. Operations Research, 44(4), 596-608. De Wolf, D., & Smeers, Y. (2000). The gas transmission problem solved by an extension of the simplex algorithm. Management Science, 46(11), 1454-1465. Demissie, A., & Zhu, W. (2015). A survey on gas pipelines operation and design optimization. IIE Annual Conference. Proceedings, 734-742. Retrieved fromhttp://ezproxy.library.ubc.ca/login?url=https://search-proquest-com.ezproxy.library.ubc.ca/docview/1791989125?accountid=14656 Djebedjian, B., El-Naggar, M., & Shahin, I. (2011). Optimal design of gas distribution network: a case study. Mansoura Engineering Journal, 36(3), 35-51. 107 Djebedjian, B., Shahin, I., & El-Naggar, M. (2008). Gas distribution network optimization by genetic algorithm. Paper presented at the Ninth International Congress of Fluid Dynamics & Propulsion, Alexandria, Egypt. Domschke, P., Geißler, B., Kolb, O., Lang, J., Martin, A., & Morsi, A. (2011). Combination of Nonlinear and Linear Optimization of Transient Gas Networks. INFORMS Journal on Computing, 23(4), 605-617. doi:10.1287/ijoc.1100.0429 El-Mahdy, O. F. M., Ahmed, M. E. H., & Metwalli, S. (2010). Computer aided optimization of natural gas pipe networks using genetic algorithm. Applied Soft Computing, 10(4), 1141-1150. doi:10.1016/j.asoc.2010.05.010 Fahimnia, B., Molaei, R., & Ebrahimi, M. (2008). Genetic algorithm optimization of fuel consumption in compressor stations. WSEAS Transactions on Systems and Control, 3(1), 1-10. Farzaneh-Gord, M., & Rahbari, H. R. (2016). Unsteady natural gas flow within pipeline network, an analytical approach. Journal of Natural Gas Science and Engineering, 28, 397-409. doi:10.1016/j.jngse.2015.12.017 Flores-Villarreal, H. J., & Ríos-Mercado, R. Z. (2003). Computational experience with a GRG method for minimizing fuel consumption on cyclic natural gas networks. Computational methods in circuits and systems applications, 90-94. Forouzanfar, M., Doustmohammadi, A., Menhaj, M. B., & Hasanzadeh, S. (2010). Modeling and estimation of the natural gas consumption for residential and commercial sectors in Iran. Applied Energy, 87(1), 268-274. doi:10.1016/j.apenergy.2009.07.008 Fügenschuh, A., Geißler, B., Gollmer, R., Morsi, A., Pfetsch, M. E., Rövekamp, J., . . . Steinbach, M. C. (2015). Physical and technical fundamentals of gas networks. by T. Koch, B. Hiller, ME Pfetsch, and L. Schewe. SIAM-MOS series on Optimization. SIAM, 17-44. Gagganapalli, S. R. (2015). Implementation and Evaluation of CMA-ES Algorithm, Master thesis, North Dakota State University, Fargo, ND. Available from NDSU Library. (http://hdl.handle.net/10365/25507) Giustolisi, O., & Berardi, L. (2011). Water distribution network calibration using enhanced GGA and topological analysis. Journal of hydroinformatics, 13(4), 621-641. Goldberg, D. E. (1983). Computer-aided Gas Pipeline Operation using Genetic Algorithms and Rule Learning. Doctoral dissertation, University of Michigan, Ann Arbor, MI. Goldberg, D. E. (1987). Computer-aided pipeline operation using genetic algorithms and rule learning. Part I: Genetic algorithms in pipeline optimization. Engineering with Computers, 3(1), 35-45. Goldberg, D. E., & Kuo, C. H. (1987). Genetic algorithms in pipeline optimization. Journal of Computing in Civil Engineering, 1(2), 128-141. 108 Gunes, E. F. (2013). Optimal Design of a Gas Transmission Network: A Case Study of the Turkish Natural Gas Pipeline Network System. Doctoral dissertation, Iowa State University, Ames, IA. Haddad, J., & Behbahani, R. (2013). Optimization of a natural Gas transmission system. International Journal of Computer Applications, 66(11). Hansen, C. T., Madsen, K., & Nielsen, H. B. (1991). Optimization of pipe networks. Mathematical Programming, 52(1-3), 45-58. Hansen, N. (2016). The CMA evolution strategy: A tutorial. arXiv preprint arXiv:1604.00772. Hansen, N., & Ostermeier, A. (2001). Completely derandomized self-adaptation in evolution strategies. Evolutionary computation, 9(2), 159-195. Harik, G. R., Lobo, F. G., & Goldberg, D. E. (1999). The compact genetic algorithm. IEEE transactions on evolutionary computation, 3(4), 287-297. Herrán-González, A., De La Cruz, J. M., De Andrés-Toro, B., & Risco-Martín, J. L. (2009). Modeling and simulation of a gas distribution pipeline network. Applied Mathematical Modelling, 33(3), 1584-1600. doi:10.1016/j.apm.2008.02.012 Holland, J. H. (1992). Adaptation in Natural and Artificial Systems. An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence. MIT Press, Cambridge, MA, USA, pp 221. Kabirian, A., & Hemmati, M. R. (2007). A strategic planning model for natural gas transmission networks. Energy Policy, 35(11), 5656-5670. doi:10.1016/j.enpol.2007.05.022 Karaboga, D. (2005). An Idea Based on Honey Bee Swarm for Numerical Optimization. Technical report-tr06, Erciyes university, Kayseri, Turkey. Karaboga, D., & Basturk, B. (2007). A powerful and efficient algorithm for numerical function optimization: artificial bee colony (ABC) algorithm. Journal of Global Optimization, 39(3), 459-471. doi:10.1007/s10898-007-9149-x Karaboga, D., & Basturk, B. (2008). On the performance of artificial bee colony (ABC) algorithm. Applied Soft Computing, 8(1), 687-697. doi:10.1016/j.asoc.2007.05.007 Khedr, A., & Tolson, B. (2015). Comparing optimization techniques with an engineering judgment approach to WDN design. Journal of Water Resources Planning and Management, 142(5), C4015014. Li, C., Jia, W., Yang, Y., & Wu, X. (2010). Analysis of operation optimization of gas pipelines based on adaptive fenetic algorithm. Paper presented at the Information Engineering and Electronic Commerce (IEEC), 2010 2nd International Symposium, Ternopil, Ukraine, 23-25 July, 2010 1-6. Li, C., Jia, W., Yang, Y., & Wu, X. (2011). Adaptive Genetic Algorithm for Steady-State Operation Optimization in Natural Gas Networks. Journal of Software, 6(3). doi:10.4304/jsw.6.3.452-459 109 Li, X., Armagan, E., Tomasgard, A., & Barton, P. I. (2011). Stochastic pooling problem for natural gas production network design and operation under uncertainty. AIChE Journal, 57(8), 2120-2135. doi:10.1002/aic.12419 Liu, E., Li, C., & Yang, Y. (2014). Optimal energy consumption analysis of natural gas pipeline. ScientificWorldJournal, 2014, 506138. doi:10.1155/2014/506138 Luongo, C. A., Gilmour, B. J., & Schroeder, D. W. (1989, April). Optimization in natural gas transmission networks: A tool to improve operational efficiency. In the Proceedings of the Third SIAM Conference on Optimization, Boston, MA, USA, April 1989. Mahlke, D., Martin, A., & Moritz, S. (2010). A mixed integer approach for time-dependent gas network optimization. Optimization Methods and Software, 25(4), 625-644. doi:10.1080/10556780903270886 Mak, T. W., Van Hentenryck, P., Zlotnik, A., Hijazi, H., & Bent, R. (2016). Efficient dynamic compressor optimization in natural gas transmission systems. Paper presented at the American Control Conference (ACC), Boston, MA, USA, 6-8, July, 2016. Martch, H. B., & McCall, N. J. (1972). Optimization of the design and operation of natural gas pipeline systems. Paper presented at the Fall Meeting of the Society of Petroleum Engineers of AIME. , San Antonio, Texas, USA, 8-11 October, 1972. Martin, A., Möller, M., & Moritz, S. (2005). Mixed Integer Models for the Stationary Case of Gas Network Optimization. Mathematical Programming, 105(2-3), 563-582. doi:10.1007/s10107-005-0665-5 Martin, D. W., & Peters, G. (1963). The application of Newton’s method to network analysis by digital computer. Journal of the institute of Water Engineers, 17(2), 115-129. Miller, B. L., & Goldberg, D. E. (1995). Genetic algorithms, tournament selection, and the effects of noise. Complex Systems, 9(3), 193-212. Molaei, R., Ebrahimi, M., Sadeghian, S., & Fahimnia, B. (2007). Genetic Algorithm-Based Optimization of Fuel Consumption in Network Compressor Stations. Paper presented at the WSEAS International Conference on Applied and Theoretical Mechanics, Spain, 14-16, December, 2007. 3(1). O'Neill, R. P., Williard, M., Wilkins, B., & Pike, R. (1979). A mathematical programming model for allocation of natural gas. Operations Research, 27(5), 857-873. Olorunniwo, F. O. (1981). A Methodology for the Optimal Design and Capacity Expansion Planning of Natural Gas Transmission Networks. Doctoral dissertation, University of Texas, Austin, TX Olorunniwo, F. O., & Jensen, P. A. (1982). Optimal capacity expansion policy for natural gas transmission networks—a decomposition approach. Engineering optimization, 6(1), 13-30. Olorunniwo, F. O., & Jensen, P. A. (1982). Natural gas transmission system optimization. Transp. Eng. J. ASCE;(United States), 108. 110 Osiadacz A.J., Bell D.J. (1990) steady-state optimization of large gas networks. In: Manley, J., McKee, S., Owens, D. (eds) Proceedings of the Third European Conference on Mathematics in Industry. European Consortium for Mathematics in Industry, vol 5. Springer, Dordrecht. Pambour, K. A., Bolado-Lavin, R., & Dijkema, G. P. J. (2016). An integrated transient model for simulating the operation of natural gas transport systems. Journal of Natural Gas Science and Engineering, 28, 672-690. doi:10.1016/j.jngse.2015.11.036 Percell, P. B., & Ryan, M. J. (1987). Steady state optimization of gas pipeline network operation. Paper presented at the PSIG Annual Meeting, Tulsa, Oklahoma, USA, 22-23 October, 1987. Puranik, Y., Kilinç, M., Sahinidis, N. V., Li, T., Gopalakrishnan, A., Besancon, B., & Roba, T. (2016). Global optimization of an industrial gas network operation. AIChE Journal, 62(9), 3215-3224. doi:10.1002/aic.15344 Radcliffe, N. J., & Surry, P. D. (1995). Fundamental limitations on search algorithms: Evolutionary computing in perspective. In Computer Science Today (pp. 275-291). Springer, Berlin, Heidelberg. Ríos-Mercado, R. Z., & Borraz-Sánchez, C. (2015). Optimization problems in natural gas transportation systems: A state-of-the-art review. Applied Energy, 147, 536-555. doi:10.1016/j.apenergy.2015.03.017 Ríos-Mercado, R. Z., Kim, S., & Boyd, E. A. (2006). Efficient operation of natural gas transmission systems: A network-based heuristic for cyclic structures. Computers & Operations Research, 33(8), 2323-2351. doi:10.1016/j.cor.2005.02.003 Ríos-Mercado, R. Z., Wu, S., Scott, L. R., & Boyd, E. A. (2002). A reduction technique for natural gas transmission network optimization problems. Annals of Operations Research, 117(1-4), 217-234. Rodriguez, G. H., Pibouleau, L. G., Pantel, C. A., & Domenech, S. (2010). Optimization of gas transmission networks under energetic and environmental considerations. International Journal of Chemical Reactor Engineering, 8(1). Ruan, Y., Liu, Q., Zhou, W., Batty, B., Gao, W., Ren, J., & Watanabe, T. (2009). A procedure to design the mainline system in natural gas networks. Applied Mathematical Modelling, 33(7), 3040-3051. doi:10.1016/j.apm.2008.10.008 Schewe, L., Koch, T., Martin, A., & Pfetsch, M. E. (2015). Mathematical optimization for evaluating gas network capacities. In T., Koch Editor (eds.), EVALUATING GAS NETWORK CAPACITIES (pp. 87-102). Society for Industrial and Applied Mathematics. 10.1137/1.9781611973693.ch5 Schmidt, M., Steinbach, M. C., & Willert, B. M. (2015). High detail stationary optimization models for gas networks: validation and results. Optimization and Engineering, 17(2), 437-472. doi:10.1007/s11081-015-9300-3 Schreider, S., Plummer, J., McInnes, D., & Miller, B. (2014). Sensitivity analysis of gas supply optimization models. Annals of Operations Research, 226(1), 565-588. doi:10.1007/s10479-014-1709-0 111 Shahin, E. I. A. S. (2012). Optimization of Gas Pipeline Networks using Genetic Algorithms, Doctoral dissertation, Mansoura University, Mansoura, Egypt. Shamir, U. (1964). Minimum Cost Design of Water Distribution Networks, Doctoral dissertation, Dep. of Civil Engineering, MIT, Cambridge, MA. Shamir, U. Y., & Howard, C. D. (1968). Water distribution systems analysis. Journal of the Hydraulics Division, 94(1), 219-234. Shamir, U. (1974). Optimal design and operation of water distribution systems. Water resources research, 10(1), 27-36. Simpson, A., Murphy, L., & Dandy, G. (1993). Pipe network optimisation using genetic algorithms. Paper presented at American Society of Civil Engineers Specialty Conference of Water Resources Planning and Management Division, Seattle, Washington, USA, 1-5, May, 1993. Singh, R. R., & Nain, P. (2012). Optimization of Natural Gas Pipeline Design and Its Total Cost Using GA. International Journal of Scientific and Research Publications, 2(8). Storn, R., & Price, K. (1997). Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11(4), 341-359. Suribabu, C. R. (2010). Differential evolution algorithm for optimal design of water distribution networks. Journal of Hydroinformatics, 12(1), 66. doi:10.2166/hydro.2010.014 Szoplik, J. (2012). The Gas Transportation in a Pipeline Network, Advances in Natural Gas Technology. In: Dr. Hamid Al-Megren (Ed.) ADVANCES IN NATURAL GAS TECHNOLOGY, InTech, Croatia, (pp. 339-358). Tabkhi, F. (2007). Optimization of Gas Transmission Networks, Doctoral dissertation, Institut National Polytechnique de Toulouse, Toulouse, France. Talatahari, S., Mohaggeg, H., Najafi, K., & Manafzadeh, A. (2014). Solving parameter identification of nonlinear problems by Artificial Bee Colony Algorithm. Mathematical Problems in Engineering, 2014, 1-6. doi:10.1155/2014/479197 Todini, E., & Pilati, S. (1987). A gradient method for the analysis of pipe networks, International Conference on Computer Applications for Water Supply and Distribution 1987, Leicester Polytechnic, U.K., 8-10 September, 1987. Todini, E. (2008). On the convergence properties of the different pipe network algorithms. Paper presented at Eighth Annual Water Distribution Systems Analysis Symposium (WDSA). Cincinnati, Ohio, United States, 27-30, August, 2008. Todinov, M. T. (2013). Flow Networks : Analysis and Optimization of Repairable Flow Networks, Networks with Disturbed Flows, Static Flow Networks and Reliability Networks, Elsevier Science, Oxford, U.K., 2013. pp. 320. Tolson, B. A., Maier, H. R., Simpson, A. R., & Lence, B. J. (2004). Genetic algorithms for reliability-based optimization of water distribution systems. Journal of Water Resources Planning and Management, 130(1), 63-72. 112 Üster, H., & Dilaveroğlu, Ş. (2014). Optimization for design and operation of natural gas transmission networks. Applied Energy, 133, 56-69. doi:10.1016/j.apenergy.2014.06.042 Vasan, A., & Simonovic, S. P. (2010). Optimization of water distribution network design using differential evolution. Journal of Water Resources Planning and Management, 136(2), 279-287. Vondráček, J., Pelikán, E., Konár, O., Čermáková, J., Eben, K., Malý, M., & Brabec, M. (2008). A statistical model for the estimation of natural gas consumption. Applied Energy, 85(5), 362-370. doi:10.1016/j.apenergy.2007.07.004 Woldeyohannes, A. D., & Majid, M. A. A. (2011). Simulation model for natural gas transmission pipeline network system. Simulation Modelling Practice and Theory, 19(1), 196-212. doi:10.1016/j.simpat.2010.06.006 Wong, P., & Larson, R. (1968). Optimization of natural-gas pipeline systems via dynamic programming. IEEE Transactions on Automatic Control, 13(5), 475-481. Wong, P. J., & Larson, R. E. (1968). Optimization of tree-structured natural-gas transmission networks. Journal of Mathematical Analysis and Applications, 24(3), 613-626. Wood, D. J., & Charles, C. O. (1972). Hydraulic network analysis using linear theory. Journal of the Hydraulics division, 98(7), 1157-1170. Wu, S., Rios-Mercado, R. Z., Boyd, E. A., & Scott, L. R. (2000). Model relaxations for the fuel cost minimization of steady-state gas pipeline networks. Mathematical and Computer Modelling, 31(2-3), 197-220. You, Y., Li, L., & Gao, L. Z. (2011). Analysis of Hydraulic Conditions in Urban Natural Gas Transmission and Distribution Network. Applied Mechanics and Materials, 121-126, 2426-2430. doi:10.4028/www.scientific.net/AMM.121-126.2426 Zavala, V. M. (2014). Stochastic optimal control model for natural gas networks. Computers & Chemical Engineering, 64, 103-113. doi:10.1016/j.compchemeng.2014.02.002 Zheng, F., Zecchin, A. C., & Simpson, A. R. (2013). Self-Adaptive Differential Evolution Algorithm Applied to Water Distribution System Optimization. Journal of Computing in Civil Engineering, 27(2), 148-158. doi:10.1061/(asce)cp.1943-5487.0000208
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Evolutionary algorithms for optimal operation of gas...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Evolutionary algorithms for optimal operation of gas networks Lee, Andy Jongsuk 2018
pdf
Page Metadata
Item Metadata
Title | Evolutionary algorithms for optimal operation of gas networks |
Creator |
Lee, Andy Jongsuk |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Natural gas is one of the cleanest fossil fuels and most reliable energy source readily available for everyday use. Due to changes in weather conditions and in economic growth and market conditions, demand for natural gas by municipal, as well as commercial and industrial consumers, has greatly increased, and delivering gas to these varied consumers has become more complex. The problem of optimally operating gas transmission networks is generally formulated to minimize the total supply cost of the network while satisfying the demand of consumers at different delivery points, at minimal guaranteed pressures. If necessary, compressors are installed and operated to supply more energy to the network and increase gas pressures, so as to compensate for pressure losses in the pipe network. The nonlinear relationships between system discharges and headlosses due to friction and mechanical devices in pipes form a complex set of nonlinear constraints in optimization models. In order to solve this complex problem, researchers have used approaches that use linear approximation in classical optimization techniques. Recently, heuristic algorithms have been applied, as these techniques can be used to find approximate solutions to complex optimization problems more quickly than classical optimization methods, and can often obtain an acceptable solution when classical methods fail. This thesis compares the performance of four heuristic algorithms for optimizing operations of a gas transmission network, namely, the: Genetic Algorithm (GA), Differential Evolution Algorithm (DE), Artificial Bee Colony Algorithm (ABC), and Covariance Matrix Adaptation Evolution Strategy (CMA-ES). The algorithms are compared in terms of their computational burden, and the quality of the obtained solutions with respect to practical management concerns. Among the four algorithms, CMA-ES is found to achieve the global optimum, compared with that generated with classical optimization, most consistent results in solutions that meet pressure and flow requirements, and conserve compressor power. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-01-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0363065 |
URI | http://hdl.handle.net/2429/64385 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_february_lee_jongsuk.pdf [ 3.88MB ]
- Metadata
- JSON: 24-1.0363065.json
- JSON-LD: 24-1.0363065-ld.json
- RDF/XML (Pretty): 24-1.0363065-rdf.xml
- RDF/JSON: 24-1.0363065-rdf.json
- Turtle: 24-1.0363065-turtle.txt
- N-Triples: 24-1.0363065-rdf-ntriples.txt
- Original Record: 24-1.0363065-source.json
- Full Text
- 24-1.0363065-fulltext.txt
- Citation
- 24-1.0363065.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0363065/manifest