02 Final

SCHOOL OF INDUSTRIAL AND INFORMATION ENGINEERING
Master of Science in Automation and Control Engineering
Object-Oriented Modelling and Simulation

of the Italian Gas Network
Thesis advisor: M.Sc. Thesis of:

Prof. Francesco Casella Simone Bosotti
Matr. 905728
Research advisor:
Prof. Emanuele Martelli
Academic Year 2019-2020

Abstract
The push towards the decarbonization is not reducing the importance of the natural
gas in the Italian and European fuel mix, because it is the fossil fuel with the lowest
environmental impact. For this reason, the optimization of the Italian transport
network is crucial not only to guarantee a reliable energy source, but also to reduce
its economic and environmental cost.
This thesis lays the foundations for a long-term collaboration between Politecnico
di Milano and Snam, the Italian transmission system operator.
In this work, a Modelica package that contains models for the simulation of the
network components, such as pipes, compressor stations, and control valves, has
been developed. This language allows to build a set of components that, thanks to
the object-oriented approach, can be easily employed in different simulations and
with various purposes.
Afterwards, a Python script has been realized. It collects the files provided by Snam
and parses them, storing in a graph the information about the network compon-
ents and configuration, the machines statuses and set-points, and the pressures and
flow-rates at the nodes. Then, the script can simplify the graph, keeping only the
fundamental components and merging the series of pipes with similar geometrical
features. Finally, it can generate the Modelica model that represents the network.
The result is a toolchain, based only on open-source software, that is able to gen-
erate a fully customizable Modelica model and the relative map for the network
visualization. This can be considered as the first step of the wider process that will
lead to the development of an optimization suite based on dynamic optimization
algorithms, e.g. MILP.
i
ii
Estratto
La spinta verso la decarbonizzazione non sta diminuendo l’importanza del gas natu-
rale nel mix energetico italiano ed europeo, dato che questo è il combustibile fossile
con il minor impatto ambientale. Per questa ragione, l’ottimizzazione della rete di
trasporto italiana è cruciale non solo per garantire una fonte di energia affidabile,
ma anche per ridurre il suo costo economico e ambientale.
Questa tesi getta le basi per una collaborazione a lungo termine tra il Politecnico di
Milano e Snam, il gestore del sistema di trasporto italiano.
In questo lavoro è stato sviluppato un pacchetto Modelica che contiene i modelli

per la simulazione dei componenti della rete, come tubi, stazioni di compressione e
valvole di controllo. Questo linguaggio permette di costruire una serie di compo-
nenti che, grazie all’approccio object-oriented, possono essere facilmente impiegati
in diverse simulazioni e con vari scopi.
Inoltre, è stato sviluppato uno script Python che raccoglie i file forniti da Snam e
li tratta salvando in un grafo le informazioni sui componenti e sulla configurazione
della rete, sugli stati e set-point delle macchine, sulle pressioni e portate nei nodi.
Quindi, questo script è in grado di semplificare questo grafo, conservando solo i suoi
componenti fondamentali e unendo le serie di tubi con caratteristiche geometriche
simili. Infine, può generare il modello Modelica che rappresenta la rete.
Il risultato è un pacchetto di strumenti, basati solo su software open-source, che è in

grado di generare un modello Modelica completamente personalizzabile e la relativa
mappa per la rappresentazione della rete. Questo può essere considerato come il
primo passo per un processo più grande che porterà allo sviluppo di una suite di
ottimizzazione basata su algoritmi di ottimizzazione dinamica, ad esempio MILP.
iii
iv
Ringraziamenti
Il primo ringraziamento va ai proff. Casella e Martelli per avermi guidato in questo

lavoro e per essere sempre stati disponibili per proficui momenti di confronto.
A Snam va la mia gratitudine per avermi dato l’opportunità di contribuire a questo

interessante progetto.
Un grazie anche a Paolo e Matteo, miei compagni di viaggio, per aver bilanciato (e
sopportato) la mia impulsività nell’approccio ai problemi. L’unica nota stonata che
rimarrà sarà il non esserci potuti incontrare di persona più spesso.
Ai miei genitori, Elena e Maurizio, va tutta la mia riconoscenza per avermi sempre
sostenuto, credendo nelle mie capacità e nelle mie aspirazioni. Non potrò mai rin-
graziarvi abbastanza per aver supportato il percorso che mi ha portato fino a qui,
anche quando nemmeno io sapevo dove fossi diretto.
A mia sorella Chiara va tutta la mia ammirazione per aver sopportato la trasforma-
zione della nostra camera in sala video-conferenze.
Ringrazio in modo speciale la mia fidanzata Federica per essere stata la fondamentale
stampella del mio inglese zoppicante durante la stesura di questa tesi, ma soprattutto
per aver sempre ascoltato tutti i miei dubbi e le mie paure e per avermi sempre
spronato a dare il massimo.
Infine ringrazio tutti gli amici, universitari e non, che mi hanno accompagnato in
tutti questi anni.
v
vi
Contents
1 Introduction 1
1.1 Description of a natural gas transport network . . . . . . . . . . . . . 2
1.2 Thesis goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Simulation program . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Network data analysis and simplification . . . . . . . . . . . . . . . . 7
2 Modelling of the Gas Network Components 9

2.1 Modelica language . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 GasNetworks package . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Methodological introduction . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Gas model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.5 Connectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.6 Pipe models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.6.1 Pipe transfer function . . . . . . . . . . . . . . . . . . . . . . 19
2.6.2 Spatial discretization . . . . . . . . . . . . . . . . . . . . . . . 23
2.6.3 Modelica models for the pipes . . . . . . . . . . . . . . . . . . 27
2.7 Compressor model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.8 Gas turbine model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.9 Gas Pumping Unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.10 Bypass Valve for Compressor Station . . . . . . . . . . . . . . . . . . 42
2.11 Compressor Station . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.11.1 Ideal control for compressor station . . . . . . . . . . . . . . . 44
2.12 Control valve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.12.1 Control valve with mass flow-rate and outlet pressure set-points 47
2.12.2 Control valve with outlet pressure set-point . . . . . . . . . . 48
2.13 Other models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.13.1 Node . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.13.2 Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
vii
viii CONTENTS
2.13.3 System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.13.4 System initializer . . . . . . . . . . . . . . . . . . . . . . . . . 51
3 Modelling of the Italian Gas Network 53

3.1 Source files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2 Python libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3 Data elaboration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.3.1 Data import and cleaning . . . . . . . . . . . . . . . . . . . . 63
3.3.2 Data merge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.4 Graph handling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.4.1 Graph generation . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.4.2 Graph simplification . . . . . . . . . . . . . . . . . . . . . . . 71
3.5 Map export . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
3.5.1 Graphical representation of nodes . . . . . . . . . . . . . . . . 88
3.5.2 Graphical representation of edges . . . . . . . . . . . . . . . . 90
3.5.3 Pop-ups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.6 Modelica model generation . . . . . . . . . . . . . . . . . . . . . . . . 97
3.6.1 writeModelicaModel function . . . . . . . . . . . . . . . . . . . 97
3.6.2 generateModelicaStationModels function . . . . . . . . . . . . 99
4 Model validation 105

4.1 Case-study 1: Central-Southern Italy . . . . . . . . . . . . . . . . . . 106
4.1.1 Description of the scenario . . . . . . . . . . . . . . . . . . . . 106
4.1.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.2 Case-study 2: Malborghetto station start transient . . . . . . . . . . 114
4.2.1 Description of the scenario . . . . . . . . . . . . . . . . . . . . 114
4.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5 Conclusions 119
5.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
List of Figures
1.1 Natural gas price evolution in the 2014-2020 time interval. . . . . . . 1

1.2 Worldwide natural gas consumption scenarios for the time interval
2018-2040. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Structure of a transport network for natural gas. . . . . . . . . . . . . 3
1.4 Consumptions and pollutant emissions caused by the transport net-
work managed by Snam Rete Gas. . . . . . . . . . . . . . . . . . . . 4
1.5 Snam network map, early 2020. . . . . . . . . . . . . . . . . . . . . . 6
2.1 Plot containing the temperature of the natural gas at the inlet of the
machines of Malborghetto, Gallese, Tarsia, and Enna. . . . . . . . . . 14
2.2 Linear interpolation of Z values obtained with ASPEN (T = 15 ◦ C)
for the mean composition. . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Modelica code of the PneumaticPort connector. . . . . . . . . . . . . 16
2.4 Modelica code of the MechPowerPort connector. . . . . . . . . . . . . 17
2.5 Modelica code of the PowerMeteringPort connector. . . . . . . . . . . 17
2.6 The equivalent RCI circuit of the infinitesimal building block for a pipe. 21
2.7 Two-port linear network described by the matrix ABCD. . . . . . . . 21
2.8 The RCR pipe building block . . . . . . . . . . . . . . . . . . . . . . 24
2.9 Example of RCR pipe, split in four finite volumes. . . . . . . . . . . . 24
2.10 The CRC pipe building block . . . . . . . . . . . . . . . . . . . . . . 25
2.11 Icon for pipes in GasNetworks Modelica package . . . . . . . . . . . . 27
2.12 Modelica code for the mass balance equation in an RCR pipe. . . . . 28
2.13 Modelica code for the momentum balance equation in a CRC pipe. . 28
2.14 Modelica code of the function that pairs the consumption positions
with the pipe finite volumes. . . . . . . . . . . . . . . . . . . . . . . . 30
2.15 Comparison between the values of ∆his obtained from the polynomial
and from ASPEN simulations. . . . . . . . . . . . . . . . . . . . . . . 33
2.16 Comparison between the values of Tout,is obtained from the polyno-
mial and from ASPEN simulations. . . . . . . . . . . . . . . . . . . . 34
ix
x LIST OF FIGURES
2.17 Example of map of a centrifugal compressor. . . . . . . . . . . . . . . 35

2.18 Icon for the CompressorPolynomials model. . . . . . . . . . . . . . . 35
2.19 Map of 23 MW gas turbine. . . . . . . . . . . . . . . . . . . . . . . . 38
2.20 Icon for the GasTurbine model. . . . . . . . . . . . . . . . . . . . . . 40
2.21 GasPumpingUnit graphical representations. . . . . . . . . . . . . . . 41
2.22 Icon for the Bypass model. . . . . . . . . . . . . . . . . . . . . . . . . 42
2.23 Modelica code of the Bypass model. . . . . . . . . . . . . . . . . . . . 43
2.24 BaseCompressorStation graphical representations. . . . . . . . . . . . 44
2.25 Graphical representation of the curvilinear coordinates used in the
ideal control. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.26 Modelica code for the ideal control equations. . . . . . . . . . . . . . 45
2.27 Code that performs the control of the bypass valve. . . . . . . . . . . 46
2.28 Modelica code of the model for the control valve with set-points for
mass flow-rate and outlet pressure. . . . . . . . . . . . . . . . . . . . 48
2.29 Modelica code of the model for the control valve with outlet pressure
set-point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.30 Icon for Node model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.31 Modelica code of the Node model. . . . . . . . . . . . . . . . . . . . . 50
2.32 Icons for the flow-rate source, the time-varying pressure source, and
the constant pressure source. . . . . . . . . . . . . . . . . . . . . . . . 50
2.33 Modelica code of the FlowSource model. . . . . . . . . . . . . . . . . 51
2.34 Icon for SystemInitializer model. . . . . . . . . . . . . . . . . . . . . 52
2.35 Modelica code for the SystemInitializer model. . . . . . . . . . . . . . 52
3.1 UML diagram describing the content of the XML file. . . . . . . . . . 54

3.2 Example of the definition of a node in the XML file. . . . . . . . . . . 54
3.3 Example of the definition of an element in the XML file. . . . . . . . 54
3.4 UML diagram describing the content of the “Mercato” CSV file. . . . 55
3.5 Lines extracted from the “Mercato” CSV file. . . . . . . . . . . . . . 55
3.6 Two lines extracted from an “NSI” file. . . . . . . . . . . . . . . . . . 56
3.7 UML diagram describing the content of the “NSI” file. . . . . . . . . 56
3.8 Example of rows extracted from the Aggregate file. . . . . . . . . . . 57
3.9 UML diagram describing the content of the “Aggregate” file. . . . . . 57
3.10 UML diagram describing the content of the REQ file. . . . . . . . . . 57
3.11 Examples of rows from the REQ file. . . . . . . . . . . . . . . . . . . 58
3.12 Block diagram representing the procedure that elaborates the data
acquired, leading to the geographical map and the Modelica model. . 60
3.13 Examples of pandas instructions used in the program. . . . . . . . . . 60
LIST OF FIGURES xi
3.14 Examples of BeautifulSoup instructions used in the program. . . . . . 61

3.15 Examples of the application of the merge algorithm with different
methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.16 Example of the use of pandas.merge in the program. . . . . . . . . . . 62
3.17 Diagram of the connections between items in different source files
used when merging data about nodes. . . . . . . . . . . . . . . . . . . 65
3.18 Diagram of the connections between items in different source files
used when merging data about edges. . . . . . . . . . . . . . . . . . . 66
3.19 Python code that performs the merge operation for both nodes and
edges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.20 Excerpt of the program code that creates the nodes in the graph and
saves into them their attributes, read from the dataframe. . . . . . . 69
3.21 UML diagram of the graph. . . . . . . . . . . . . . . . . . . . . . . . 71
3.22 Detail of the map representing the original graph around the Enna
pumping station. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.23 Examples of different placements of control valves in the network. . . 75
3.24 Flowchart representing the sub-functions contained in edge removal
and their exit conditions. . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.25 Code that performs the removal of every edge that has one end not
connected to another element. . . . . . . . . . . . . . . . . . . . . . . 77
3.26 Code of the function int edges removal. . . . . . . . . . . . . . . . . . 77
3.27 int edge simplification, initial condition. . . . . . . . . . . . . . . . . . 78
3.28 int edge simplification, wrong procedure - step 1. . . . . . . . . . . . . 79
3.31 int edge simplification, correct procedure - step 1. . . . . . . . . . . . 80
3.35 Problematic configuration for int edge simplification function. . . . . 81
3.36 safe add edge example. . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.37 Example of valve configuration that can be solved only using slowEle-
mentRemoval. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.38 Comparison between the original, detailed map and the simplified one. 85
3.39 Flowchart of the external loop around the function followPipe. . . . . 85
3.40 Flowchart of the internal structure of the function followPipe. . . . . 86
3.41 Representations of the network. . . . . . . . . . . . . . . . . . . . . . 89
xii LIST OF FIGURES
3.42 A real example of the representation of nodes in the map. . . . . . . . 89

3.43 Code that manages the creation of nodes on the map. . . . . . . . . . 91
3.44 Part of the map containing every category of edges of the network. . . 92
3.45 Code that generates the circles highlighting active machines on the
map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
3.46 Example of a node pop-up collapsed and expanded. . . . . . . . . . . 94
3.47 Example of a pop-up referred to an edge. In this case, it describes
the attributes of a control valve. . . . . . . . . . . . . . . . . . . . . . 94
3.48 Example of a compressor station pop-up: the basic information page
on the left, the advanced data on the right. . . . . . . . . . . . . . . . 96
3.49 Example of a pop-up containing the information about a composite
pipe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
3.50 Examples of the code that generates the map. . . . . . . . . . . . . . 97
3.51 Comparison between the output obtained using a TimeTable and an
IntegerTable with the same parameter table. . . . . . . . . . . . . . . 100
3.52 Example of tpMarket variable taken as input by importTimeSeries-
Market. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.53 Example of the conversion performed by the function convertArray-
ForModelica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
3.54 Example of the behaviour of the convertTimeSeriesForModelica func-
tion with active relaxationTransient option. . . . . . . . . . . . . . . 103
4.1 Map of the Central-Southern Italy network. . . . . . . . . . . . . . . 106

4.2 Evolution of the mass flow-rates at the injection points in the network.107
4.3 Comparison between the SCADA measurements and the simulated
values for the inlet pressure of the Enna station. . . . . . . . . . . . . 108
values for the pressures at Mazara terminal outlet, Gela terminal,
and Enna inlet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
values for the pressures at Tarsia station outlet, Melizzano station
inlet, and Gallese station inlet. . . . . . . . . . . . . . . . . . . . . . . 110
values for the mass flow-rate in Melizzano station. . . . . . . . . . . . 111
values for the dimensionless rotational speed of the active units in
Enna, Tarsia, and Gallese stations. . . . . . . . . . . . . . . . . . . . 112
LIST OF FIGURES xiii

values for the fuel mass flow-rate of the active units in Enna, Tarsia,
and Gallese stations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.9 Map of the North-Eastern Italy network. . . . . . . . . . . . . . . . . 115
4.10 Pressure profile of the nodes 5394 and 5404. . . . . . . . . . . . . . . 115
values for the pressure at Malborghetto station outlet, at Node 21349,
and at Istrana station. . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4.12 Behaviour of the simulated mass flow-rates that flows in Malborghetto
compressor station, in its machine and in its bypass valve. . . . . . . 117
xiv LIST OF FIGURES
List of Tables
2.1 Calculations for CH4 (90%) and C2 H6 (10%) mix using GERG-2008 . 13
2.2 Mean composition that has been considered in ASPEN simulations
to obtain the interpolating polynomial for Z . . . . . . . . . . . . . . 15
2.3 Coefficients evaluated for the considered mixture at T̄ = 15 ◦ C. . . . . 15
2.4 Values of the sound velocity for the considered composition of natural
gas, evaluated with REFPROP. . . . . . . . . . . . . . . . . . . . . . 19
2.5 τRI time constants and parameters used for the evaluation for a set
of pipes extracted from the real network. . . . . . . . . . . . . . . . . 23
2.6 τRC,n time constants for a set of pipes extracted from the real network. 27
2.7 Reference values for the compressors map. . . . . . . . . . . . . . . . 31
2.8 Example of coefficients for the interpolating polynomial for ∆his . . . 33
2.9 Example of coefficients for the interpolating polynomial for Tout,is . . 33
2.10 Reference values for gas turbine maps. . . . . . . . . . . . . . . . . . 38
2.11 Scheme of the possible different values for the w at the outlet. . . . . 47
2.12 Scheme of the possible different values for the pout . . . . . . . . . . . 47
3.1 Colours of the segments representing the different edge categories. . . 90
4.1 Numerical comparison between the simulation and the SCADA meas-
ures for the station of Enna. . . . . . . . . . . . . . . . . . . . . . . . 108
ures for pressures at Mazara del Vallo terminal, Gela terminal, and
Enna station inlet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
ures for pressures at Tarsia station outlet, Melizzano station inlet, and
Gallese station inlet. . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
ures for the mass flow-rate in Melizzano station. . . . . . . . . . . . . 111
xv
xvi LIST OF TABLES
4.5 Toolchain execution times and simulation parameters in the Case-

Study 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
ures for the dimensionless rotational speed of the active units in Enna,
Tarsia, and Gallese stations. . . . . . . . . . . . . . . . . . . . . . . . 112
ures for the fuel mass flow-rate of the active units in Enna, Tarsia,
and Gallese stations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.8 Toolchain execution times and simulation parameters in the Case-
Study 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Chapter 1
Introduction
Natural gas represents the fossil fuel with the lowest environmental impact because
it produces a low quantity of CO2 per unit energy, thanks to the highest ratio
between the number of C-H bonds and number of carbon atoms. Natural gas also
has a modest sulfur content, and it does not contain solid residue. Furthermore,
the natural gas price declined in the last decade (Fig. 1.1) [21] and, thanks to the
spread of the liquefied natural gas, the difficulties of long-range distribution have
been overcome.
Figure 1.1: Natural gas price evolution in the 2014-2020 time interval.
Source: [21]
The International Energy Agency (IEA) scenarios, developed before the economic
crisis related to the COVID-19 pandemic, reported an increase in the natural gas
consumption in the next ten years, even in the most conservative case (Fig. 1.2).
1
2 CHAPTER 1. INTRODUCTION
Figure 1.2: Worldwide natural gas consumption scenarios for the time interval 2018-
2040.
Source: [14]
Despite the precariousness related to the demand contraction due to the current
crisis, it is legit to foresee a consumption increase, even in the first phase of the
energy transition towards decarbonization [13, 15]. Natural gas is widely used in
the industrial and residential sector, and it is strategic in electricity production
because it is employed in high flexibility and efficiency plants to compensate for the
volatile production of renewable sources.
Moreover, the network, used for natural gas, is operated to transport biomethane
and, in some trials, a mixture containing hydrogen [34].
1.1 Description of a natural gas transport net-

work
The function of a national network for the transport of Natural Gas (NG) is to take
the fluid from the nodes of import, extraction, or regasification and to bring it to the
redelivery points. These can be represented by important and delimited users, e.g.
power stations, or, otherwise, connections to local networks operated by different
societies (Fig. 1.3) [35].
The movement of the gas in the network is guaranteed by the compression power
generated by specific stations, that, consuming a small fraction of the fuel pumped
in pipes, provide the necessary amount of energy to the fluid to overcome the static
head and the pressure losses caused by friction.
In the context of European regulations, companies that manage the NG transport

1.1. DESCRIPTION OF A NATURAL GAS TRANSPORT NETWORK 3
Figure 1.3: Structure of a transport network for natural gas.

Source: [35]
networks are also responsible for the dispatching of the gas. For this reason, they
have the role of Transmission System Operator (TSO). The dispatching process
consists of the balancing of the quantities of gas injected, extracted and collected
in the network, with the aim of respecting the operational pressure limits of the
infrastructure [7].
While in electric grids the power balance must be respected within a very short time
span of about 15 seconds, in gas networks, thanks to the storage capacity provided
by reservoirs and even by the pipes themselves, the balancing can be performed
on time scales of days. However, the European regulation states that the usage of
pipe storage, called linepack, must be “consistent with the economic and efficient
operation of the transmission network ” [7].
The current environmental policies in Europe promote the reduction of greenhouse

gases and pollutants emissions through the efficiency optimization of the existent
processes and the increase of the power quota produced by renewable sources. These
indications are critical for the energy sector, and, as a consequence, the network
operators must apply important measures to fulfil them. In this context, software
that simulates and optimizes the operation of gas networks, and in particular their
compression stations, can be useful, and it can reduce consumptions, lowering both
the emissions of pollutants and the costs. In addition, the future plans for the
development of the network are considering the adoption of compressors moved by
electrical motors. This approach can not only reduce the quantity of burned gas,
Figure 1.4: Consumptions and pollutant emissions caused by the transport network
managed by Snam Rete Gas.
Source: [36]
but it can also be used as an active resource for the balancing of the electrical grid
[9].
1.2 Thesis goals

The work of this thesis is a part of a wider project that sees the collaboration of
Snam, the Italian network TSO, and Politecnico di Milano. The final goal of this
industrial activity is the development of software that simulates and optimizes the
operation of the national network for natural gas transport. By doing so, it is
possible to reduce the operational costs and the pollutants emissions caused by the
machines, respecting the technical and environmental limits and guaranteeing the
gas flow [36].
This work, that is the result of this thesis and the one of two students belonging to
the Department of Energy in Politecnico di Milano, is focused on the realization of
the first model for dynamic simulation. The principal steps are:
Data analysis on the material provided by Snam about the topology of the
network and the injections and extractions flows of natural gas;
Creation of the models of the network components, e.g. pipes and gas pumping
units, using the Modelica language;
Development of a script, programmed in Python, that reconstructs the network

1.2. THESIS GOALS 5
from the data, then derives a simplified version of it for the simulation and
optimization suites;
Simulation and validation of the model.
In Figure 1.5, the map of the Italian transport network is represented. Its commercial
management is entrusted to Snam Rete Gas, and it is split into the National Network
(NN) and the Regional Network (RN).
The National Network is composed of pipes and stations sized and tested consider-
ing the requirements set by the imports, the regasification terminals, the national
production plants of natural gas and biomethane, and the storages. Its function is
to transfer substantial quantities of gas from the inlet nodes to the consumption
areas. The National Network comprehends five entry points for the foreign import
lines, three connections to regasification plants, and thirteen compression stations.
In addition, the NN is connected to twelve storage sites, which, in the scope of this
work, will be represented in the model as injection/extraction points, avoiding the
explicit modelling of their compression units.
The Regional Network is composed of the remaining part of pipelines managed by

Snam Rete Gas not included in the NN and the relative plants. Its purpose is to
move the natural gas in smaller regions, bringing it from the collection points on
the NN to the final redelivery points. These pipes are characterized by a diameter
that becomes smaller with the distance from the National Network [33]. The mass
flow-rate withdrawn by customers connected to the RN is, in the model, measured
by means of an aggregate value that is placed in correspondence to the connection
between the National and Regional distribution system.
In the last months of 2020, the Trans-Adriatic Pipeline (TAP), the new pipeline
that carries the natural gas from Azerbaijan to southern Italy crossing Greece and
Albania, became operational, but it is not represented on the map in Figure 1.5.
The information provided by Snam for this work contains almost the entire Na-
tional Network, integrated with a part of the Regional one, localized in Toscana.
It contains:
6188 nodes, where 1094 are redelivery points and 28 entry points;
1840 pipes;
2251 gate valves;

Figure 1.5: Snam network map, early 2020.

Source: Snam
15 compressor stations;
68 metering stations;
64 control valves;
29 check valves;
2589 joints.
In order to have a model suitable for the dynamic simulation, the system must be
transformed into a simplified one, which retains the parts that are relevant for the
simulation. In this way, the number of equations to solve is also lowered. This
procedure is performed in the Python script and it is described in Section 3.4.2. In
the end, the simplified network contains only pipes, compressor stations and control
valves.
1.3 Simulation program

The modelling of the transport network is crucial in this thesis. For this reason,
the work of the first months was focused on the development of the models of the
components. To obtain simulations representative of the entire network, it is not
1.4. NETWORK DATA ANALYSIS AND SIMPLIFICATION 7
possible to overlook a precise comprehension and representation of the phenomena

taking place in every element.
To achieve this goal effectively, the Modelica programming language has been em-
ployed. In fact, Modelica is equation-based and object-oriented. This means that
the equations written in the models are effortlessly readable because there is a strict
correspondence between the language jargon and the mathematical symbols. Be-
sides, the object-oriented approach guarantees that the base components can be
simply connected to each other. Additionally, modularity will be useful when the
code developed for this thesis will be employed to realize an optimizer based on
MILP or MINLP algorithms. The advantages of Modelica are explained in detail in
Section 2.1.
1.4 Network data analysis and simplification

The transport network managed by Snam is a complex system, composed of thou-
sands of different components, and its configuration changes frequently. In fact,
compression stations can be switched on or off at any moment and the state of a
group of valves can be changed by the operators, leading to a different pipe layout.
The files, provided by Snam, that describe the system and that are used as sources
for this work, are:
A file that contains the static description of the elements and their connections.
This is exported from SIMONE, a software that is currently used by Snam to
simulate the network.
A file that describes the states of every controlled element, e.g. valves and
stations. This is the log of the Supervisory Control And Data Acquisition
(SCADA) system of the network.
A file that reports the mass flow-rates in correspondence with import nodes
and redelivery points. This is the register where the flow-rate sensors store
their data.
A table that matches components with their statuses and measures, also ex-
ported from SIMONE.
The description of the network components and their associations with statuses
and measures is static. In fact, this information can vary only after a physical
modification of the items, e.g. a new pipe is installed, or one sensor is added in a
certain location. On the contrary, the operational data are constantly acquired, and
the files contain time series for every considered attribute.
For these reasons, the approach followed in this thesis is to create a fully automatic
workflow, that acquires the files, then handles them, and finally generates a Modelica
model that is stand-alone and can directly be simulated. For this purpose, Python
has been chosen as a programming language. It has a wide availability of open-
source libraries, especially in scientific fields, and, between them, there is also one
that supports the direct integration of the Modelica simulation part.
The used libraries, the applied procedures, and other details are explained in Chapter 3.
The simulation that results from the source files automatic elaboration has been
validated with the real data contained in the SCADA logs. The network has been
generated and simplified using a snapshot of the configuration in a certain time,
and then it has been simulated for six hours. Then, the evolution of some interest-
ing quantities has been compared with the real trend measured by the sensors of
the physical network. The procedures, the initial assumptions and the results are
presented in Chapter 4.
The result of this work is a complete toolchain that, starting from the source files
provided by Snam, is able to create customizable Modelica models representing the
network. These can be generated for different purposes: either a more detailed
model for the simulation, or a simpler one for the MILP optimization. Besides, in
literature, there are no other examples of gas networks of this size and complexity
modelled using the Modelica language.
Chapter 2
Modelling of the Gas Network

Components
In this work, the modelling of the gas network is of crucial importance. Even if
the single components are described with a few equations each, the complexity of
the system is given by the numerousness of parts and required connections. Indeed,
the simplified system consists of 13 compressor stations, composed of three to five
groups of compressor and turbine each, about 160 pipes, characterized by different
lengths, diameters, and inclinations, and thirty control valves, with their set-points
of pressure or mass flow-rate.
To master this variety, the best approach is to build modular models that can be
easily connected but that are written in a way that is completely independent from
the components they will be connected to. This idea is the core of the object-oriented
paradigm, and, in the modelling area, it is perfectly declined by Modelica.
2.1 Modelica language

Modelica [3] is a non-proprietary, object-oriented, equation-based language developed
since 1996 by the non-profit Modelica Association. This language allows modelling
the dynamic behaviour of physical systems, using algebraic and differential equations
to describe the component’s models, and to create multi-domain complex systems.
The association maintains also a set of libraries for electrical, mechanical, fluid, and
control systems, in addition to the open-source Modelica Standard Library.
The language and the libraries are used in several simulation environments, such
9
10 CHAPTER 2. MODELLING OF THE GAS NETWORK COMPONENTS
as Dymola, Amesim, and MapleSim. To develop the project of this thesis, Open-
Modelica, an open-source simulation environment supported by the Open Source
Modelica Consortium, has been chosen [27].
The main advantages of this language are: [5]
A-causal, declarative modelling: the model of every component is composed

of a set of algebraic and differential equations that describes how this object
behaves and not how the equations are solved. The boundary conditions are
not necessarily defined as inputs or outputs: the model is the same, irrespective
of what is connected to.
Code transparency: equations can be written as they appear in textbooks, the

program manages the procedure to solve them. This eases off the code writing
and verification.
Encapsulation: every model has rigorously defined connectors that allows to

bound it to any other compatible component, regardless of their internal de-
tails.
Inheritance: models can be defined with a hierarchical approach. One basic

model, that contains all the common variables and equations, can be extended
with specific features.
Multi-physics modelling: being Modelica an object-oriented language, it is

straightforward to combine models that belong to different physical domains.
Re-usability: the previous ideas are strong incentives toward reuse of models.
Since every model is well-defined, regardless the environment it is used in, it
is possible to place it in other projects without any additional effort.
2.2 GasNetworks package

One of the main objectives of this thesis is to build a Modelica package that can
be used to simulate the behaviour of gas networks. A package is a specialized class
that may only contain declarations of classes and constants. These classes can be
models, functions, connectors, types, or other sub-packages [2].
The hierarchical structure allowed by Modelica packages eases the organization of

their contents, for example dividing basic components from tests or ancillary func-
tions from types definitions.
2.3. METHODOLOGICAL INTRODUCTION 11
The GasNetworks package contains nine sub-packages:
Interfaces: the connectors used to describe the ports of the components;
Sources: models that generate pneumatic quantities and that are used as
boundary conditions;
BaseClasses: components that contain only the parts that are in common
between similar models. They are then extended to create complete compon-
ents;
Components: models of items, used to create other models;
Test: models that represent simple situations in which the components are
tested;
Media: models that contain the state equations of the natural gas, and the
correlated parameters;
Types: package where the variable types used in this package are defined;
Choices: contains the definitions of the enumerations used to let the user
choose a parameter within a list;
Utilities: contains the custom functions created for this package.
In this thesis, the content of Interfaces, Sources, BaseClasses, Components, Media,

and Utilities is analysed.
2.3 Methodological introduction

In transport and distribution networks for natural gas, laws state that price-related
measurements must be referred to standard conditions. This guarantees that the
quantities are evaluated independently of the pressure and temperature in a specific
location and time instant [1]. In particular, the standard cubic meter (Sm3 ) is
employed: it defines the gas quantity contained in a cubic meter of the same gas
at p0 = 101325 Pa and T0 = 288.15 K. Using the ideal gas model (2.1), it is
straightforward to convert a standard cubic meter in kilograms, knowing the molar
mass (MM ) (2.2) [19]. The dependence on the molar mass links the definition of
this measurement unit to the composition of the considered gas. For this reason,
it is crucial to have a well-known mixture in every point where this conversion is
performed.
RT0 Sm3
vmol = = 23.645 · 10−3 (2.1)
p0 mol
where R = 8.314 J/mol K is the universal gas constant.
23.645 · 10−3
1 kg = Sm3 (2.2)
MM
Since the sensors integrated in the Supervisory Control And Data Acquisition (SCADA)
are used to measure flow-rates, the utilized unit is the standard cubic meter per hour
(Sm3 /h), that is then converted into kg/s.
In addition, in the SCADA, gauge pressures are acquired. In this work, only absolute
pressures are employed, so the data are converted, as explained in Section 3.3.1.
2.4 Gas model

In order to properly describe the natural gas, a state equation that considers the
non-ideal behaviour of the hydrocarbon mixture has to be taken into account. In
fact, in the transport network, it is affected by variations caused by the different
temperature (T ), pressure (p), and composition, represented by the vector x that
contains the molar fractions of the components. The typical pressures measured in
pipes are between 45 bar and 75 bar, and this leads to sensible variations from the
model of ideal mixture of ideal gases. For this reason, the compressibility factor
Z(T, p, x) is introduced. For a substance, it is defined by the equation:
ρid v
Z(T, p) = = (2.3)
ρ vid
where ρ represents the gas density, v the specific volume, and the subscript id the
parameters evaluated with the ideal gas assumption.
At high pressure, the intermolecular forces cause an increase of the density with
respect to the value predicted in the hypothesis of ideal gas in the same pressure
and temperature conditions. Besides, phenomena related to the non-idealities of the
mixture must be taken into account. The law that describes this aspect is:
pv = Z (T, p, x) R̂T (2.4)
where R̂ = R/MM .
Values in Table 2.1 are generated by the software REFPROP [25] using a mixture
2.4. GAS MODEL 13
T [K] p [bar] ρ [kg/m3 ] Z [-] cp [kJ/ (kg · K)]

288.15 1.0133 0.73959 0.99757 2.1248
288.15 45.000 36.718 0.89239 2.5213
288.15 75.000 66.158 0.82546 2.9089
Table 2.1: Calculations for CH4 (90%) and C2 H6 (10%) mix using GERG-2008
Source: [25]
of methane and ethane. It is evident how this behaviour moves away from the ideal
one.
It is not possible to use a constant value for the compressibility in the gas model.
So, Z has to be estimated using an additional equation. In literature, there are
two groups of state equations. The first one contains cubic relationships such as
Peng-Robinson, Redlich-Kwong (RK), and Redlich-Kwong-Soave (RKS). These are
not commonly used in the simulation and optimization field because, being third-
order equations, it is requested an additional computational effort to choose the
correct solution. Moreover, they are general purpose models that can be applied to
a variety of gases and situations, so they can be either inefficient or less accurate for
well-defined mixtures and conditions.
The other group is represented by empirical exponential equations. These relation-

ships, e.g. AGA and Papay, are widely used in researches in natural gas field for
their accuracy and numerical reliability [18, 11].
Another simplifying assumption is to consider the gas as isothermal [18, 10, 28, 38].
This is reasonable when the network is composed by long pipes, mainly underground.
In fact, modelling the heat exchange process between the gas and the ground would
require an expensive analysis of the soil and an additional computational workload
to evaluate the energy balance equation. Moreover, the thermal effects become
negligible in a short distance after the compression in a station. So, the natural
gas temperature T̄ is considered fixed in the network and equal to 15 ◦ C. For this
reason, equation 2.4 becomes
pv = Z (p, x) R̂T̄ (2.5)
This assumption is lifted when the gas enters a gas pumping station. In fact, the
behaviour of the compressor is strongly affected by the variability of the temper-
ature. In addition, at the inlet of every gas pumping unit there is a temperature
Figure 2.1: Plot containing the temperature of the natural gas at the inlet of the
machines of Malborghetto, Gallese, Tarsia, and Enna.
sensor for the gas flow. This information has been used in the two case-studies ana-
lysed in Chapter 4 to evaluate the value of Z at the actual temperature, considered
fixed during the simulation time interval. This assumption is justified by the data
provided by Snam, shown in Figure 2.1. It is possible to see how the values does
not substantially change during the considered interval. The behaviour of the Mal-
borghetto time series is motivated by the state of the machine: for the first 2000
second the unit is stopped and the gas is stagnant. After a short transient, the
temperature does not change anymore.
In other papers, the dependence of the compressibility factor from the mixture
composition x is neglected [22] and a value depending only from pressure is chosen,
or a constant mean value [12]. This is justified by the fact that the composition
of the mixture undergoes small variations during the transport in the network. In
particular, the molar fraction of methane, that is the most important component of
NG, is nearly constant.
In this work, the effects of the composition variations in the network are analysed
in a simplified way, in order to reduce the computational complexity during the
simulation. So, the compressibility factor has been described with a polynomial
that depends on the pressure [23]. It is obtained interpolating data from the soft-
ware ASPEN [16], simulating the gas behaviour in different conditions of pressure
and temperature. Three composition have been considered: the first, reported in
Table 2.2, is a mean value of the different compositions measured by Snam in the
different locations on the territory. The second one is the mixture that enters the
station of Enna, and it is mainly composed by the mixture imported from Algeria
and, secondary, from Libya. It is used for the first case-study. The third composition
is analysed near the inlet of the Malborghetto station, near the Austrian border,
and it is employed in the second case-study.
2.4. GAS MODEL 15
CH4 C2 H6 C3 H8 C4 H10 C5 H12

0.926 0.046 0.0078 0.002 0.00044
CO2 O2 N2 C6 H14 He
0.0077 2.00 · 10−5 0.00953 0.00012 0.00034
Table 2.2: Mean composition that has been considered in ASPEN simulations to
obtain the interpolating polynomial for Z
The correlation obtained is the 2.6, where pressures p are expressed in bar and
temperature T in ◦ C.
Z(T, p) = k0 + k1 · p + k2 · T + k3 · p · T + k4 · p2 + k5 · T 2 (2.6)
In Figure 2.2, the temperature is assumed constant T̄ = 15 ◦ C. The polynomial can

be simplified, keeping only the terms k0 and k1 , whose values are in Table 2.3
k0 k1
0.9858 -0.00174165
Table 2.3: Coefficients evaluated for the considered mixture at T̄ = 15 ◦ C.
The state equation 2.5, substituting the obtained linear relationship, becomes the
equation 2.7. p
pv = k0 + k1 5 R̂T̄ (2.7)
10
Figure 2.2: Linear interpolation of Z values obtained with ASPEN (T = 15 ◦ C) for

the mean composition.
Evaluated expression: Z = −0.0017 · p + 0.9859, R2 = 0.9994
2.5 Connectors
In Section 2.1, one of the described advantages of Modelica language is the En-
capsulation. In fact, it states that every model must connect to other components
through connectors. In the GasNetwork library, two types of connectors have been
employed.
The first group is represented by Causal connectors. These items have a predefined
direction of flow for the information. In fact, they are used as inputs or outputs. In
this work, only inputs from the Modelica Standard Library have been employed, e.g.
RealInput or IntegerInput. With these components, a model is able to receive
from another block a value, e.g. a real or integer number.
The second group contains the A-Causal connectors, commonly also called ports.
These objects represent a link between the models that does not have a predefined
direction. A port contains one effort variable and one flow variable. When two or
more ports are connected together using the connect construct, Modelica generates
a set of equations that states that:
All the effort variable must assume the same value;
The sum of all the flow variables, assumed positive when entering the com-
ponent, must be equal to zero.
In this work, three custom ports have been used.
PneumaticPort connector To connect the components that carry the natural

gas to each other, the PneumaticPort connectors are used (Fig. 2.3). Its effort
variable is the pressure and the flow one is the mass flow-rate. This model is then
extended to create two different icons of the same item. In this way, it is possible
to differentiate the inlet and the outlet ports.
connector PneumaticPort
flow T yp es .M as sF lo wR at e w " Mass flow - rate flowing into port " ;
T y p e s . A b s o l u t e P r e s s u r e p " Absolute pressure " ;
end PneumaticPort ;
Figure 2.3: Modelica code of the PneumaticPort connector.
MechPowerPort connector The connector MechPowerPort is used to mech-

anically connect compressors and gas turbines. In fact, its effort variable is the
2.6. PIPE MODELS 17
rotational speed expressed in rad/s and the flow one is the mechanical power (Fig.
2.4).
connector MechPowerPort
flow Types.Power P_mech " Mechanical Power at shaft " ;
M o d e l i c a . S I u n i t s . A n g u l a r V e l o c i t y rot_speed " Rotational speed of
the shaft in rad / s " ;
end MechPowerPort ;
Figure 2.4: Modelica code of the MechPowerPort connector.
PowerMeteringPort connector The PowerMeteringPort is a fictitious con-

nector that is used to automatically sum the overall thermal power consumption of
the gas turbines in the network (Fig. 2.5). In fact, in the gas turbine model (Sec.
2.8), this port is used to connect the machine to the system (Sec. 2.13.3). The flow
variable of this port is the thermal power and the effort one is a dummy entity that
is not employed.
Using the property of the connection between ports, in the PowerMeteringPort of

the system, all the values of the flow variables are summed. This means that, in the
port, it is possible to read the total power consumption of the entire network.
connector PowerM eteri ngPort

flow Types.Power P " Power entering the component " ;
Types.PerUnit dummy " Dummy effort variable " ;
end P owerMe tering Port ;
Figure 2.5: Modelica code of the PowerMeteringPort connector.
2.6 Pipe models

In order to describe the flow of the natural gas inside the network pipes, the following
assumptions are made:
Newtonian fluid: the NG respects Stokes’ hypoteses;
Compressible fluid: the density of the gas considerably varies in the network;
Single-phase fluid: the NG is not affected by phase transitions;
One-dimensional flow: the lenght of the pipes is much greater than its diamet-
ers;
Isothermal flow: the gas temperature is considered as fixed, as explained in

Section 2.4;
Turbulent flow: the velocity of the fluid in pipes and its density lead to a
more important inertial term with respect to the viscous contribution in the
Reynolds number.
The dynamic behaviour of a pipe in unsteady conditions is governed by Euler equa-

tions for compressible fluids [11].

∂ρA ∂w

 + =0
∂t ∂x

(2.8)
 ∂w + A ∂p + ρAg dz + cf ρ ω |u| u = 0


∂t ∂x dx 2
In 2.8, A and ω represent the area and the wet perimeter of the pipe, x is the
coordinate parallel to the pipe axis and z is the elevation, w is the mass flow-rate
and u the gas velocity, cf is the friction coefficient (defined in 2.9), and g is the
gravitational acceleration. The first equation represents the mass balance, while the
second one the momentum balance. Since the isothermal assumption has been made
(Sec. 2.4), the energy equation has been neglected.
The mass balance describes how the gas reacts to differences in the flow-rate chan-
ging its density, while momentum balance is the equation describing how phenomena
like friction, inertia and pipe inclination modify the gas pressure.
The friction has been considered in a simplified way too: indeed, conditions in the
pipes always guarantee very high values for the Reynolds number. This means that
the gas is always considered to be in turbulent flow conditions. As a result, the
first term in the Colebrook-White equation 2.9 becomes negligible, so that only the
second term is used in the model to evaluate the friction coefficient cf [18].

1 2.51 ε h ε i−2
√ = −4 log10 √ + → cf = −4 log10 (2.9)
cf 2 · Re · cf 3.71D 3.71D
Using this value, it is possible to obtain the distributed pressure loss:
cf L
∆pf rict = ρ ωu2 (2.10)
2 A
where L is the total pipe length.
In the second equation of 2.8, the kinetic term has been neglected. In fact, starting
2.6. PIPE MODELS 19
T [K] p [bar] c [m/s]

288.15 45.00 406.63
288.15 60.00 403.95
288.15 75.00 403.40
Table 2.4: Values of the sound velocity for the considered composition of natural
gas, evaluated with REFPROP.
from the complete momentum balance equation 2.11 and considering it in steady-
state condition (2.12), the term ∂ρAu2 /∂x can be rewritten as 2.13, where c is the
speed of sound in the fluid.
∂w ∂ρAu2 ∂p dz cf
+ +A + ρAg + ρ ω |u| u = 0 (2.11)
∂t ∂x ∂x dx 2
dρ̄Āū2 dp̄ dz cf
+ Ā + ρ̄Āg + ρ̄ ω̄ |ū| ū = 0 (2.12)
dx dx dx 2
∂ ρ̄Āū2 w̄2 w̄2 d ū2 dp

d
= =− ρ̄Ā = −A (2.13)
∂x dx ρ̄Ā ρ̄2 Ā2 dx c2 dx
dp̄
The final term in 2.13 can be compared to Ā dx from 2.12. The speed of sound at
the working conditions, evaluated using REFPROP [25], is close to 400 m/s (Tab.
2.4), while the usual gas velocities in the network are lower than 10 m/s. For this
reason, the assumption 2.14 is valid, and it leads to neglect the kinetic term also in
the dynamic case.
ū2 dp̄

dp̄
Ā 1 − 2 ' Ā (2.14)
c dx dx
2.6.1 Pipe transfer function

This paragraph aims to evaluate the time constants of the phenomena taking place
in pipes.
The mass balance equation can be linearized around the steady-state ∂ w̄/∂x = 0:
∂∆ρ ∂∆w
A + =0 (2.15)
∂t ∂x
Substituting ∆ρ with the value that can be extracted from (2.4) and integrating for
x from 0 to the total length L, equation (2.16) is obtained.
AL ∂∆p
= ∆win − ∆wout (2.16)
z T̄ R̂ ∂t
Considering the mass flow-rate deviation as an electrical current and the pressure
deviation as a voltage, it is possible to build the electrical equivalent of this model.
In particular, the fraction that multiplies ∂∆p/∂t is the capacity of this equivalent
component.
AL
C= (2.17)
z T̄ R̂
The same procedure can be applied to the momentum balance. The pipe is supposed
dz
to be horizontal to simplify the analysis, so the term ρAg dx vanishes. The steady-
state is
∂ p̄ cf ω w̄ |w̄|
A + =0 (2.18)
∂x 2A2 ρ̄
The linearization and the integration for x from 0 to L provides
L ∂∆w cf Lω cf Lω |w̄| w̄
+ (∆pout − ∆pin ) + 3
w̄∆w − ∆ρ = 0 (2.19)
A ∂t ρ̄A 2A3 ρ̄2
The term containing ∆ρ can be expressed as in 2.20: the numerator is equal to the
pressure loss in (2.10) and the denominator is the mean value of the pressure in the
pipe.
cf Lρ̄2 ū2 A2 ω cf ωū2 L ∆pf rict ∆pf rict ∆pf rict

∆p = ∆p = = p̄
∆p = ∆p (2.20)
2ρ̄2 A3 Z̄ R̂T̄ 2AZ̄ R̂T̄ ρ̄Z̄ R̂T̄ Z̄ R̂T̄ p̄
Z̄ R̂T̄
Finally, the value of ∆pf rict /p̄ can be considered small, because, in the real data
provided by Snam, the pressures are between 50 and 70 bar and the maximum value
of ∆pf rict is near to 15 bar. For this reason, the term depending on ∆ρ can be
neglected.
From the equation 2.19, two equivalent quantities can be evaluated: an inertance
and a resistance.
L cf Lω ∆pf rict
I= R= 3
w̄ = 2 (2.21)
A ρ̄A w̄
The overall equivalent circuit for a pipe is obtained connecting an infinite number of
2.6. PIPE MODELS 21
Figure 2.6: The equivalent RCI circuit of the infinitesimal building block for a pipe.
blocks similar to the one in Figure 2.6, each composed by the series of the inertance
( LI · dx) and the resistance ( R
L
· dx), with the capacity ( CL · dx) connected between
the end of the resistance and the ground.
This circuit is an example of the well-known Telegrapher’s equations problem; in

this case, the conductance to the ground is null, and I, R, and C are distributed
parameters.
In the literature, the exact solution of this problem is available [29]. It represents the
transmission line as a two-port linear network (Fig. 2.7), described by the matrix
ABCD (2.22).
" # " #" #
V1 A B V2
= (2.22)
I1 C D I2
The most common boundary condition, that a pipe can undergo in the considered
gas network, is to be connected to the outlet of a compressor station and to the
inlet of another one. This is equivalent to prescribe the inlet and outlet flow-rates
of the pipe. Considering varying the flow-rate of the first station, I1 also changes.
On the contrary, fixing the value of w in the second station keeps I2 = 0. Then, at
the outlet, the pressure (V2 ) is monitored. Solving the matrix product in 2.22, the
Figure 2.7: Two-port linear network described by the matrix ABCD.

Source: [29]
relationship between I1 and V2 is obtained.
1
I1 = C · V2 → V2 = · I1 (2.23)
C
The value for the coefficient C is reported in 2.24, and the parameters γ and Z0 are
in 2.25. In 2.26, the distributed parameters for the pipe are defined.
sinh (γL)
C= (2.24)
Z0
s
p (r + jω i)
γ= (r + jω i)(g + jω c) Z0 = (2.25)
(g + jω c)
R I C G
r= i= c= g= =0 (2.26)
L L L L
Finally, the transfer function is obtained substituting all the parameters into 2.23.
q
(r+jω i)
pout Z0 (jω c)
GRCI (jω) = = = p
win sinh (γL) sinh (r + jω i)(jω c) L
q q q q (2.27)
r i R
jω c
1 + jω r jω C
1 + jω RI
GRCI (jω) = √ q = √ q
i I
sinh jω rc 1 + jω r L sinh jω RC 1 + jω R
The presented RCI model is characterized by two states (p and w). Analysing the
transfer function, it is possible to understand at which time-scales the full RCI model
is necessary to correctly describe the behaviour of the pipe, and on which times a
simpler RC model is sufficent.
This simplification is profitable for two reasons:
The steady-state initialization for a big and complex system like the considered
network is not always a viable option. To circumvent this problem, the start
values for both the states are needed. Unfortunately, the data available for
the network contain only a few values of mass flow-rate. So, a sufficient ini-
tialization for the flow-rates is not available.
The computational cost increase caused by the increased number of states to

2.6. PIPE MODELS 23
L [m] ∆pf rict [Pa] ρ̄ [kg/m3 ] u [m/s] τRI [s]

17 114 37 796 43.3 2.6 25.49
20 948 7 986.8 61.1 0.8 64.50
66 724 153 173 51.9 3.0 33.57
96 205 254 554 54.1 3.7 37.83
148 325 675 096 54.2 4.9 29.18
189 382 544 357 46.9 4.2 34.26
Table 2.5: τRI time constants and parameters used for the evaluation for a set of
pipes extracted from the real network.
consider is notable, but it could not lead to a valuable increase of the quality
of the results, considering the purposes of the work.
The transfer function for the RC model can be directly obtained imposing I = 0 in
(2.27). q
R
jω C
GRC (jω) = √ (2.28)
sinh jω RC
q
The two transfer functions are equivalent when the term 1 + jω RI is equal to one.
This happens when ω 1/τRI , where τRI is the value defined by the equation 2.29.
I L w̄ ρ̄ ūL
τRI = = = (2.29)
R A 2 ∆pf rict 2 ∆pf rict
In the data provided for this work, the measurements are sampled with a period of
4 minutes. To be comparable, the time scales considered in the simulations must
be similar. In Table 2.5, the values of τRI for a set of pipes extracted from the real
network are calculated.
Since the time constants are lower than the reference value of 240 seconds, it is clear
how, for every case, the RC model is sufficient to describe all the relevant network
dynamics. In other words, the inertial term in the momentum balance equation can
be neglected.
2.6.2 Spatial discretization

As described in Section 2.6.1, pipes can be modelled combining a resistor, a capacitor
and an inductor. The integration process applied to obtain equations 2.16 and
2.19 approximates as uniformly distributed all the values of space-varying quantities
to their mean values. This is not acceptable while modelling a pipeline for gas
Figure 2.8: The RCR pipe building block
transport, because the spatial distribution of quantities like pressure and mass flow-
rate has crucial importance. The finite volume method has been used to solve
this problem: every pipe is divided into a finite number of elements with definite
length. In the previous section, the equivalent electric circuit for a pipe has been
represented as a series of R and C. Since the capacitance expression is obtained from
the mass balance equation, it is dependent from the pressure in the pipe. Assuming
that the fluid flows from the left-hand side to the right-hand side of Figure 2.8, it
encounters at first the resistance, undergoing the pressure loss, and then it enters
the capacitance. In this way, the density is evaluated with the pressure that the gas
has at the end of the pipe. Since the flow direction is not prescribed, components
can be placed with the gas flow entering the pipe from the outlet. In this case, the
fluid at first encounters the capacity: its density is evaluated using the pressure at
the beginning of the pipe. Since pipes can be long, leading to pressure losses of
several bar, the gas density and the mass storage in the capacity can substantially
change.
The solution developed for this project is to split the equivalent resistor into two
parts with resistance R/2N and to connect the capacitance between the middle
point and ground. This allows flipping the component, or inverting the gas flow,
preserving the value of the pressure drop on the entire pipe and improving the
density estimation in the volume. When more RCR blocks are connected to each
Figure 2.9: Example of RCR pipe, split in four finite volumes.

2.6. PIPE MODELS 25
Figure 2.10: The CRC pipe building block
other, the two resistances placed between two capacitances are connected in series.
So, their values can be summed (Fig. 2.9). In the following paragraphs, this model
will be called RCR pipe.
This model is preferable when the pipe is represented with a small number of finite
volumes. In fact, this model describes more accurately the pressure and density
profiles thanks to the higher number of equivalent resistances per finite volume.
This advantage is reduced when the number of elements is increased. In addition, if
the pipe is split into more blocks, the outermost resistors value becomes very small.
This can lead to numerical problems in some cases.
To solve this problem, it is possible to build the complementary model, called CRC
pipe, where a resistor is placed between two capacitors (Fig. 2.10). To avoid the
above-mentioned problems, in this phase of the work the CRC model has been
preferred.
The employment of the finite volume method can restore the spatial description
of the quantities, but it leads to an increase of the computational workload of the
solver. For this reason, the length of a finite volume must be pondered to balance
the spatial accuracy of the profile and the execution time of the simulation. To do
this, the analytical solution of the system is considered.
The equation that describes the dynamical behaviour of a CRC pipe coincides with
the well-known description of the thermal diffusion phenomenon, and it is reported
L2
in (2.30), where a = RC .
∂w ∂ 2w
=a 2 (2.30)
∂t ∂x
Its analytical solution is described in 2.31 [30].

Domain 0 ≤ x ≤ l
w = f (x) at t=0 (initial condition)
w = g1 (t) at x=0 (boundary condition)
w = g2 (t) at x=l (boundary condition)
∞
an2 π 2 t

2X nπx
w(x, t) = sin exp − 2 Mn (t) (2.31)
l n=1 l l
Z l
anπ t
Z 2 2
nπξ an π τ
Mn (t) = f (ξ) sin dξ + exp 2
[g1 (τ ) − (−1)n g2 (τ )] dτ
0 l l 0 l
As boundary condition, a flow-rate step is applied to the inlet (g1 (t) = 1), keeping the
outlet w constant (g2 (t) = 0). The initial profile considered is uniform (f (x) = 0).
The solution obtained is in 2.32.
∞
an2 π 2 t

2anπ X nπξ
w(x, t) = 2 sin 1 − exp − 2 (2.32)
l n=1 l l
In that expression, there is an exponential term, that depends on n. The coefficient

that multiplies the t in the exponent represents the time constants of the dynamical
phenomena of the RC pipe model. The dependence by n describes that the propaga-
tion of the flow-rate wave is composed by the superposition of infinite modes, each
with the time constant 2.33.
L2
τRC,n = 2 2 (2.33)
an π
A sinusoidal expression is multiplied to the exponential term, and it represents

the spatial modes of the propagation considered. The coefficient that multiplies ξ
represents the space constant of the mode (2.34).
L
ξRC,n = (2.34)
n
The term n in the series correspond to the number of finite volume of the model.
To have a suitable description of the dynamic behaviour, τRC,n should be smaller
than the sampling time of the measures. The correspondent ξRC,n represents the
maximum length for the finite volume. In the Table 2.6, the first five τRC,n are
evaluated for the same set of pipes considered in Table 2.5. Is it possible to see how,
choosing n = 5, the model is accurate for every pipe considered.
2.6. PIPE MODELS 27
L [m] τRC,1 [s] τRC,2 [s] τRC,3 [s] τRC,4 [s] τRC,5 [s]
17 114 6.27 1.57 0.70 0.39 0.25
20 948 2.48 0.62 0.28 0.16 0.10
66 724 74.46 18.62 8.27 4.65 2.98
96 205 139.72 34.93 15.52 8.73 5.59
148 325 439.83 109.96 48.87 27.49 17.59
189 382 614.59 153.65 68.29 38.41 24.58
Table 2.6: τRC,n time constants for the same set of pipes of Table 2.5.
2.6.3 Modelica models for the pipes
Figure 2.11: Icon for pipes in GasNetworks Modelica package
Pipes represent the largest group of components in network models. In the Modelica
package, they are described using models that extend one base partial class, that
contains the definition of interfaces and parameters, called BasePipe. This is used
to define every variable that describes physical parameters that are not different in
the two variants of the extended models. For example, the base model contains the
area and the wet perimeter of the pipe, the elevation of the inlet and the outlet, and
the roughness.
The extended models are:
PipeCRCDistributedConsumption
PipeRCRDistributedConsumption
A pipe model is characterized by
geometrical parameters: length, diameter, roughness, and elevation
number of finite volumes
initialization parameters
gas consumption distribution parameters
The equations are the simple translation into Modelica language of the mathematical
expressions described in the previous sections, repeated for every finite volume using
for loops.
// Mass balances : all the volumes have the same size

for i in 1 : N - 1 loop
der ( M [ i ] ) = w [ i ] - w [ i + 1 ] - wcons [ i ] ;
end for ;
M [ i ] = gas [ i + 1 ] .rho * V / ( N - 1) ;
end for ;
Figure 2.12: Modelica code for the mass balance equation in an RCR pipe.
// Momentum balances : all the volumes have the same size

u [ i ] = w [ i + 1 ] / ( gas [ i ] .rho * A ) ;
u2 [ i ] = M o d e l i c a . F l u i d . U t i l i t i e s . r e g S q u a r e ( u [ i ] , unom * delta ) ;
dp [ i ] = p [ i ] - p [ i + 1 ] ;
dp_frict [ i ] = cf / 2 * gas [ i ] .rho * omega * u2 [ i ] * l / A ;
dp_head [ i ] = gas [ i ] .rho * g * dz ;
dp [ i ] = dp_frict [ i ] + dp_head [ i ] ;
// Momentum balance , entering side , entire volume
end for ;
dp_head_tot = sum ( dp_head ) ;

dp_frict_tot = sum ( dp_frict ) ;
dptot = sum ( dp ) ;
Figure 2.13: Modelica code for the momentum balance equation in a CRC pipe.
In Figure 2.13, it can be seen the only peculiarity of this translation procedure: the
use of the regSquare function. It substitutes the u·|u| term providing a better numer-
ical consistency. It is equivalent to x^2 * sgn(x) when abs(x) >> unom * delta,
and it is equal to x * delta when abs(x) >> unom * delta, avoiding the prob-
lems that are risen by the zero derivative in x = 0. In fact, this would lead to a
singular Jacobian during the solution of non-linear equations involving the pressure
losses.
There are also equations that link pressure and flow-rate of the first and the last
elements to the quantities in the ports of the model.
2.6.3.1 Initialization
A critical aspect to consider when modelling components that contain differential

equations is the initialization. In fact, choosing wrong initial equations can lead to
simulation errors.
2.6. PIPE MODELS 29
Both models of pipes have five different initialization methods:
Steady-state: pressures and mass flow-rates in the pipe start from their equilib-
rium value. This method guarantees the absence of transients on the quantities
considered, but it is also difficult to evaluate, because the equilibrium must be
found at system level, and not only on the component. The code involved in
this case is:
initial equation
der ( M [ i ] ) = 0;
end for ;
Flow-rate profile: a guessed mass flow-rate profile is used to initialize the

states. It is obtained from wguess and considering the consumption placed on
the pipe;
Pressure profile: the states are initialized using a guessed pressure profile,
that is evaluated from the parameters pguessIn and wguess using an initial
algorithm;
Simple pressure profile: the pressure profile is evaluated using the value of
pguessIn and pguessOut. The values in the between are interpolated between
the values at the boundaries with the function linspace;
No initialization: the states initial values are not specified.
The guess values are parameters attributed to the model when it is instantiated.
2.6.3.2 Consumption distribution algorithm
Even if the main purpose of the network considered in this thesis is to transport
natural gas through the country, there are a lot of points in the middle of the pipes
where gas is extracted or injected. In several cases, these extractions represent
a consistent part of the total mass flow-rate, so moving them from their actual
position to the nearest node in the network can lead to important errors in the
friction evaluation.
For this reason, the script that optimizes the network before transforming it into
a Modelica model deals with this aspect. In fact, it stores into two designated
arrays the value of these mass flow-rates (consumptionArray) and their positions
(positionArray), using a linear coordinate on the pipe length.
Since the pipe is represented using a finite number of volumes, a procedure that pairs
every consumption with the most appropriate volume is required. For this reason,
the function BuildConsuptionMatrix, shown in Figure 2.14, has been developed. It
takes as input the number of volumes of the pipe Nvol, the number of consumption
placed on it Ncons, the total pipe length L and the positions array positions. Using
two nested for loops, the function saves, into a Nvol * Ncons matrix, a “1” in every
cell that matches a consumption with a finite volume.
function B u i l d C o n s u m p t i o n M a t r i x
input Integer Nvol ;
input Integer Ncons ;
input Types.Length L ;
input Types.Length positions [ Ncons ] ;
output Real matrix [ Nvol , Ncons ] ;
algorithm
matrix : = zeros ( Nvol , Ncons ) ;
for i in 1 : Nvol loop
for j in 1 : Ncons loop
if positions [ j ] >= L / Nvol * ( i - 1) and
positions [ j ] < L / Nvol * i then
matrix [i , j ] : = 1;
end if ;
end for ;
end for ;
end B u i l d C o n s u m p t i o n M a t r i x ;
Figure 2.14: Modelica code of the function that pairs the consumption positions
with the pipe finite volumes.
As a consequence, the output of this function is the matrix consMatrix that, mul-
tiplied with consumptionArray, gives the array wnom of mass flow-rate extracted
from every finite volume.
This procedure can also work when the consumptions in the network are considered
as time-varying. In fact, every cell in consumptionArray can contain the reference to
the correspondent block that generates the signal that describes the mass flow-rate
of the consumption.
2.7. COMPRESSOR MODEL 31
2.7 Compressor model

In the Italian transport network, centrifugal compressors are employed to move the
natural gas. This type of machines has a larger operational range with respect to
the axial compressors, and they are less prone to the risk of surge or choking when
the flow-rate is different from the design value. In addition, these appliances are
characterized by good efficiencies and high pressure ratios, and so, it is sufficient to
use only one stage.
To model these components, the polynomial relationships described in the manual of

the software SIRE2000 [6] are used. This software is currently used by Snam to run
static simulations of the pumping units placed in the network. The characteristic
curves of every component are described by two pairs of polynomials that describe
the generated adiabatic head (H) and the efficiency (η) as function of the number
of revolutions per unit of time (N ) and of the volumetric flow-rate of natural gas
Q = w/ρin . These quantities, to be used in the expressions, must be dimensionless
and they are marked with a star. The reference value for N is the nominal rotational
speed of the considered machine Nnom , for Q and H the values are hard-coded in
the SIRE2000 software and they are listed in Table 2.7.
N Q H
N∗ = Q∗ = H∗ = (2.35)
Nnom Qref Href
Href [m] Qref [m3 /s]

10 000 10 000/3 600
Table 2.7: Reference values for the compressors map.
The polynomials 2.36 and 2.37 describe the compressor head. The structure of the
equations is the same, but, depending on the operating conditions of the machine,
their coefficients are different. The equation 2.36 is valid when it is in the normal
range, 2.37 when it is near to the absolute choking. The manual refers to the latter
using the name “choking zone”.
H ∗ = k1n + k2n Q∗ + k3n N ∗ + k4n N ∗ Q∗ + k5n Q∗ 2 + k6n N ∗ 2 (2.36)
H ∗ = k1c + k2c Q∗ + k3c N ∗ + k4c N ∗ Q∗ + k5c Q∗ 2 + k6c N ∗ 2 (2.37)

The polynomials that describe the compressor efficiency have a different structure,
depending on the working region considered. The equation 2.38, that describes the
normal region, has five terms, while the 2.39, valid in the choking area, six. Both
expressions are then multiplied by a corrective coefficient Cη−corr , that considers the
drift of the machine from its maps.
" 2 #
Q∗ Q∗

η = k1n + k2n + k3n + k4n N ∗ + k5n N ∗ 2 · Cη−corr (2.38)
N∗ N∗
∗
∗ ∗2 Q ∗ ∗ ∗
η = k1c + k2c N + k3c N + k4c + k5c Q + k6c N Q · Cη−corr (2.39)
N∗
The application of these polynomials requires two variables between H, N , η and Q.

During the network simulation, the adiabatic head and the volumetric flow-rate are
usually known, and so, the solver can evaluate both η and N . The former is linked
to the pressure difference on the compressor ∆p by ∆his , the isentropic enthalpic
change, as showed in 2.40.
∆his = gH (2.40)
Knowing the flow-rate, the adiabatic head, and the efficiency, it is possible to eval-
uate the mechanical power that the compressor requires from the connected gas
turbine (2.41).
Q ρin ∆his
P = (2.41)
η
The enthalpy difference in 2.40 is related to the ∆p by the relationship that describes
the compression along a real gas isentropic transformation, which would require to
evaluate the deviations from the ideal behaviour [24]. To avoid the exact imple-
mentation of these calculations, a polynomial that interpolates the values of ∆his ,
obtained from the software ASPEN [16], is employed. These simulations, that use
the RKS state equation, have been performed for a range of pressure ratios. At first,
the natural gas at 15 ◦ C with the composition in Table 2.2 has been considered, ob-
taining the coefficients in Table 2.8. Then, to improve the quality of the results, a
set of polynomials has been developed, representing different mixtures at different
temperatures. The simulated points have been interpolated using the least squares
method. The fitting equation is reported in 2.42: it requires pressures expressed in
bar and it returns the value for ∆his in kJ/kg.
∆his = k0 + k1 pout + k2 pin + k3 pout pin + k4 p2out + k5 p2in (2.42)

Figure 2.15: Comparison between the values of ∆his obtained from the polynomial
and from ASPEN simulations.
k0 k1 k2 k3 k4 k5
6.0965 4.3899 -4.6004 -0.013894 -0.012502 0.028191
Table 2.8: Example of coefficients for the interpolating polynomial for ∆his
From the same ASPEN simulations, the outlet temperature data have been used
to develop the polynomial for the isentropic outlet temperature. The equation is
reported in 2.43 and, in Table 2.9, the values for the same mixture of natural gas
are shown. This expression requires the pressure in bar and gives as result the
temperature in ◦ C.
Tout,is = k0 + k1 pout + k2 pin + k3 pout pin + k4 p2out + k5 p2in (2.43)
In Figures 2.15 and 2.16, the results of the polynomials are compared with the
respective values from the ASPEN simulations.
k0 k1 k2 k3 k4 k5
20.5293 2.38631 -2.57524 -0.00269682 -0.00821598 0.0125114
Table 2.9: Example of coefficients for the interpolating polynomial for Tout,is
Figure 2.16: Comparison between the values of Tout,is obtained from the polynomial
and from ASPEN simulations.
Finally, using the value of η, the real temperature at the compressor outlet is estim-
ated with the 2.44. Then, the result is compared with the technical limit of 50 ◦ C,
imposed to avoid thermal stresses on the materials of the outlet pipes.
Tout,is − Tin
Tout = Tin + (2.44)
η
In Figure 2.17, an example of map of a centrifugal compressor is depicted: the

red line represents the surge limit, the green one the relative choking boundary (a
parabola that crosses the origin that divides the normal area from the choking one),
the blue curves represents lines where the efficiency is constant, the black ones the
trajectories with constant rotational speed. The map does not define the working
region with too high flow-rates because η, in those points, is poor due to the choking.
This behaviour is described by the black parabola.
To model this component in Modelica language, at first, the base class BaseCom-
pressor is created. In this item, the parameters, the variables, and the equations
that are in common with all the models of compressors are present. In fact, the
compressor that uses the polynomial relationships is too complex to be implemen-
ted in the MILP optimizer that will be implemented in the second phase of this
project. For this reason, a reduced model will be developed, using the base class
as starting point. The advantage of this programming approach is to guarantee the
external compatibilty of the future blocks with the actual architecture. In fact, the
interfaces are included in the base model, and so they will be intherited by the later
Figure 2.17: Example of map of a centrifugal compressor.

Source: Snam
extensions.
The model for the simulation is called CompressorPolynomials and it is represented

by the icon in Figure 2.18. The interfaces of this model are:
Inlet and outlet pneumatic ports (light blue circles): they are employed to
connect the compressor to the network and to allow the flow of the gas;
Mechanical port (yellow circle with red border): it is used to connect this
machine with the relative gas turbine. The angular velocity and the mechanical
power are exchanged through it, mimicking the physical connection given by
the shaft;
Integer input (orange triangle): it is used to determine if the appliance is active

and how many compressors working in parallel it is equivalent to. It is linked
to the variable Mp, and its purpose is further explained in Section 2.9.
The equations implemented in this model are quite complex, and the solver included
Figure 2.18: Icon for the CompressorPolynomials model.

in OpenModelica can encouter some difficulties. For this reason, the operator ho-
motopy and the custom function smoothStep have been used.
The homotopy operator [26] allows to transition from a simplified model to the full
version during the initialization procedure. This is very useful when the non-linear
equations have to be solved by iterative procedures and it can strengthen the solving
algorithm. In fact, the simplified model should be written in a way that leads the
solver to convergence even if the start values are imprecise.
This is obtained in the equation 2.45, summing up the two versions of the equations
with the variable lambda as weight. This value is then changed gradually by the
solver from zero to one during the iterations.
lambda · actual + (1 − lambda) · simplified (2.45)
In the considered case, the simplified equation (2.46) is the linearized version of
the parametric curves for the head, with the dimensionless variables.
H ∗ Href Q∗ Qref

= N∗ − −1 (2.46)
Hnom Qnom
The actual equation is the complete expression for adiabatic head. This contains
the function smoothStep in order to have a smooth transition between the normal
and the choking ranges.
The function smoothStep [31], using the function hyperbolic tangent, generates a
value that can be approximated to zero when the considered value is much lower
that the reference one, and it is similar to one in the opposite case. Multiplying
this modulating function with the expressions for H in the two operational areas,
the transition between the two is smoothed and the eventual discontinuities are
cancelled. The same consideration is applied to the equations for η.
In the considered fleet of compressors, several machines have only one operational
region. Therefore, they are modelled using only a series of coefficients for H and for
η.
2.8. GAS TURBINE MODEL 37
2.8 Gas turbine model

Every compressor placed in the stations is moved by a dedicated gas turbine. The
two machines are mechanically coupled since they are mounted on the same shaft.
For this reason, they rotate at the same speed and the turbine must satisfy the
compressor power request. The thermal input is provided by the extraction of a
small quantity of the gas that flows in the station. For this application, the most
appropriate machines are simple cycle aero-derivative gas turbines, thanks to their
small size and higher efficiency with respect to heavy-duty ones, especially when
power is lower than 70 MW [19].
The performances of these appliances are expressed as maximum available power

(Pmax ) and heat rate (HR) and they are reported in specific maps provided by the
manufacturers. The heat rate is the inverse of the gas cycle efficiency (2.47) and,
together with the produced power, it allows to estimate the required thermal input
(2.48) and the consequent mass flow-rate spilt from the network (2.49).
1
ηT G = (2.47)
HR
Pf uel = P · HR (2.48)
Pf uel
wGN = (2.49)
LHVGN
The operation maps are realized by the manufacturers that test the turbines in ISO
conditions, that are [17]:
Ambient temperature: 15 ◦ C;
Ambient pressure: 101 325 Pa;
Relative air humidity: 60 %;
Absence of pressure losses at inlet and outlet of the machine.
These diagrams are embedded in the software SIRE2000 and they are combined
with the corrective factors (fW T , fwp , fHR−T ) used to refer the real ambient condi-
tions to the ISO conditions. In addition, SIRE2000 applies two further coefficients
(Cpower−corr , CHR−corr ) to consider the performance worsening due to the components
aging and fouling.
Figure 2.19: Map of 23 MW gas turbine.

Source: Snam
The curves are described by equations that are written using adimensional quantit-
ies, marked with a star. Their definitions are:
N HR P Pmax Tamb
N∗ = ; HR∗ = ; P∗ = ∗
; Pmax = ; T∗ = (2.50)
Nnom HRref Pref Pref Tref
where N is the rotational speed, HR represents the heat rate, P is the power of the
gas turbine, Pmax is the maximum power available, and Tamb represents the ambient
temperature.
To be coherent with the Modelica code, where units from SI are used, the reference
values are converted as stated in Table 2.10. The nominal rotational speed Nnom
is expressed in RPM for simplicity, since it is natively supported by Modelica, and
derives from the machines datasheets.
In 2.51 and 2.52 the adimensional polynomial expressions, used to evaluate the heat
rate and the maximum power for the turbine in the actual working conditions, are
HRref [J/J] Pref [W] Tref [K]

10 000/3 600 107 288.15
Table 2.10: Reference values for gas turbine maps.

2.8. GAS TURBINE MODEL 39
reported.
∗
= k1 + k2 N ∗ + k3 N ∗2 · fW T · fW p · Cpower−corr

Pmax (2.51)
HR∗ = k1 +k2 N ∗ +k3 PISO

∗
+k4 N ∗ · PISO
∗
+k5 N ∗2 +k6 PISO
∗2

·fHR−T ·CHR−corr (2.52)
The corrective factors that consider the environment are expressed in 2.53, 2.54, and
2.55.
fW T = k1 + k2 T ∗ + k3 T ∗2 + k4 T ∗3 (2.53)
pamb
fW p = (2.54)
pamb−ref
fHR−T = k1 + k2 T ∗ + k3 T ∗2 + k4 T ∗3 (2.55)
The expression 2.54 needs, as parameter, the ambient pressure. It can be estimated
from the elevation using the barometric formula. This expression can be obtained
considering an infinitesimal volume of air, with height dz. The pressure difference
between the upper and the lower surfaces can be considered as completely hydro-
static 2.56. Then, using the ideal gas law to evaluate the air density 2.57, a first
order differential equation is obtained, and its general solution is 2.60. The boundary
condition 2.61 is then imposed, leading to the barometric formula 2.62.
dp = −ρg dz (2.56)
p MMaria
dp = − g dz (2.57)
RT
dp
Z Z
MMair
= − g dz (2.58)
p RT
MMair
ln(p) = − gz + C (2.59)
RT

MMair
p = K · exp − gz (2.60)
RT
p (z = 0) = 101325 Pa (2.61)
Figure 2.20: Icon for the GasTurbine model.

MMair gz
pamb = 101325 · exp − Pa (2.62)
RT
Finally, it is necessary to check if the gas turbine is working in an acceptable re-

gion. This means that the power produced must be lower than the maximum value,
referred to the same rotational speed and in the same conditions (ISO or ambient).
The gas turbine model is implemented in Modelica with the name GasTurbine.
As shown in Figure 2.20, the model has four interfaces to connect it to other items:
Mechanical port (yellow circle with red border): it is used to connect the gas
turbine with the relative compressor. The angular velocity and the mechanical
power are exchanged through it, mimicking the physical connection given by
the shaft. It is called powerPort.
Pneumatic port (light blue circle): it is employed to connect the turbine to

the network. Through this port, the amount of fuel consumed by the machine
is spilt from the network. Its name is fuelPort.
Integer input (orange triangle): it is used to determine if the appliance is active

and how many turbines working in parallel it is equivalent to. This value is
linked to the variable Mp and its purpose if further explained in Section 2.9.
Real input (blue triangle): it is utilized as input for the dimensionless rota-
tional speed of the machine. It is connected to the variable Nstar.
The variable Mp can assume only integer values and it describes the state of the gas
turbine. Mp = 0 represents a switched off machine, Mp = 1 a single active appliance,
Mp >= 2 multiple identical machines working in parallel.
The rotational speed input is not explicitly given when the gas turbine is embedded
in the GasPumpingUnit model (Sec. 2.9). In fact, in that case, the variable Nstar
2.9. GAS PUMPING UNIT 41
(a) Icon of the model (b) Internal diagram representing the con-
nections between the interfaces and the
components.
Figure 2.21: GasPumpingUnit graphical representations.
is evaluated backwards by the solver starting from the compressor flow-rate and
head, dependent from the pressure difference between inlet and outlet. The input
connector has been implemented to simplify the creation of simpler models, for
example for testing and validation purposes.
Inside the model, a hidden port of the class powerMeteringAdaptor is present and
it is employed to communicate the thermal input consumed by the turbine to the
System object. This is useful because, since there is only one System in the entire
network, it aggregates the total power consumption of every active machine.
2.9 Gas Pumping Unit

The pair composed by one (equivalent) gas turbine connected to one (equivalent)
compressor is described by a composite Modelica model called GasPumpingUnit. Its
external interfaces (Fig. 2.21a) are:
Two pneumatic ports, inlet and outlet, used to connect the unit to the gas
network, depicted with two light blue circles;
One integer input Mp (orange triangle);
One real input Nstar (blue triangle).
These ports are only used to route the flows and the signals from the external
environment to the connectors of the internal models, as represented in Figure 2.21b.
This container model is useful during the construction of the compression station
representations. In fact, during this process, the number of connections required

for each compressor-turbine couple is reduced from seven to four (six to three if
the connection for Nstar is not required). In addition, this guarantees that every
compressor is always connected to the relative turbine, and that both of them receive
the same Mp signal. This is crucial because, when Mp assumes the value zero, the
compressor and gas turbine pair is inactive. In their variables some placeholders
values are stored, to avoid numerical problems at the start-up phase. When Mp
is greater than zero, the machines are active, and its integer value represents the
number of active groups in parallel represented in the model. In fact, if two identical
compressors are coupled with two identical gas turbines, they constitute two twins
gas pumping units. To reduce the computational cost, it is very useful to aggregate
these machines into a single component, and then multiply the mass flow-rate, the
thermal power and fuel consumption by Mp.
2.10 Bypass Valve for Compressor Station

In the real compressor stations, the machines are connected to the network through
a complex system of valves and check-valves that routes the gas in the desired
direction. In particular, when there are no active units in the station, the gas flow
must have an alternative path to pass through the station. In the model, this system
is simplified introducing a bypass valve.
The component has three interfaces (Fig. 2.22): inlet and outlet pneumatic ports
and a Boolean input open (depicted by the purple triangle). This signal, when it is
false, stops the flow inside the valve, as described in the code in Figure 2.23.
The pressure loss is evaluated using a linear model instead of a quadratic one. In
fact, the latter could have lead to numerical problems without increasing the overall
precision at system level, since the pressure loss caused by an open bypass valve is
small.
Figure 2.22: Icon for the Bypass model.

2.11. COMPRESSOR STATION 43
w = if open then wNom / dpNom * dp else 0;

inlet.w = w ;
outlet.w = - inlet.w ;
dp = inlet.p - outlet.p ;
Figure 2.23: Modelica code of the Bypass model.
2.11 Compressor Station

The Compressor Station model is used to represent the plants where the gas that
flows in the network is compressed to compensate the pressure losses caused by the
motion in the pipes.
Every station in the network has its own configuration (numbers and types of ma-
chines, operational limits, etc.). For this reason in the GasNetwork package there
is not a CompressorStation model, but only the BaseCompressorStation. In this
object, there are only the components that are in common with every station. Then,
the Python code (Sec. 3.6.2) extends it, declaring the necessary items and connect-
ing them in the proper way.
The base class contains (Fig. 2.24):
The inlet and outlet pneumatic ports;
Three real input connectors, used to acquire the set-points for mass flow-rate,
minimum inlet pressure, and maximum outlet pressure;
The Bypass model;
A Joint model, used as flow-rate sensor for the ideal control;
The connections between these components.
Then, in every particular instance of compressor station, the following components

are added:
Many GasPumpingUnits as in the real plant;
One Mp input port per GasPumpingUnit;
The equations representing the ideal control of the station;
The equations that guarantee the equidistance criterion, used to distribute the
workload on the active units of the station;
(a) Icon of the model (b) Internal diagram representing the con-
nections between the interfaces and the
components.
Figure 2.24: BaseCompressorStation graphical representations.
The bypass valve control.
2.11.1 Ideal control for compressor station

The compressor stations are regulated by control system that monitors three values,
comparing them to three set-points:
The mass flow-rate that passes through the station wSP ;
The minimum pressure at the inlet port pin−SP ;
The maximum pressure at the outlet port pout−SP .
Since the goal of this work is not to accurately represent the behaviour and the
regulation system of the stations, a simplified version is employed in the model.
This logic uses two curvilinear dimensionless coordinates (r and s) that describe
the distance of the working point from the set-points (Fig. 2.25). The point where
r = 0 and s = 0 represents the situation where the monitored variables are on their
set-point values. When both r and s are greater than zero, the station is in its
normal working range: it is following and respecting its flow-rate set-point. When
one of the two variables is negative, the station has engaged one of the pressure
controls: w is reduced to keep the problematic quantity inside the boundaries.
In the actual Modelica code (Fig. 2.26) some additional parameters are added.
pInNom, pOutNom, and wNom are used to normalize the values used in the control
and to make r and s dimensionless for increased numerical robustness. In addition,
2.11. COMPRESSOR STATION 45
Figure 2.25: Graphical representation of the curvilinear coordinates used in the ideal
control.
w = if s >= 0 and r >= 0

then wSP elseif s < 0 and r > 0
then wSP + s * wNom else wSP + r * wNom ;
pout + ( w - wNom ) * impedanceCoeff * pOutNom / wNom =
if s < 0 and r > 0 then pOutMax
elseif s > 0 and r > 0 then pOutMax - s * pOutNom
else pOutMax - s * pOutNom ;
pin - ( w - wNom ) * impedanceCoeff * pInNom / wNom =
if s > 0 and r < 0 then pInMin
elseif r > 0 and s > 0 then pInMin + r * pInNom
else pInMin + r * pInNom ;
Figure 2.26: Modelica code for the ideal control equations.
the impedanceCoeff term is added putting a small impedance between the state
connected to the port and the internal state. In fact, a case that has to be considered
is the connection of a pipe to the inlet port of the station. Through the connect
declarations, this state is matched to the internal state of the compressor. When
the flow-rate control is active, the two states are independent; instead, when the
inlet pressure control is engaged, and in case that the impedance is not inserted,
the equation imposes an equivalence between the states, and they become a single
entity. This would lead to a problem of variable index that Modelica is not able to
manage. Placing that impedance, decouples the two states, avoiding this problem.
This control sets the parameters for the entire station, but it does not prescribe how
to split the load on the different active machines. This action is performed by the
equidistance criterion: it states that, when multiple compressors are active, they
must work at the same operational distance from the limit curve for surge, defined
in 2.63.
Q∗2 − Q∗2
surge
dop = ∗2 (2.63)
Qc − Q∗2surge
where Qsurge represents the volumetric flow-rate at the surge limit and Qc is the one
at absolute choke limit.
A number of equations equal to the number of active compressors minus one, stating
that their distances are the same, will thus be included in the station model.
This means that, if the compressors are different, the flow-rate is not equally divided.
Finally, in Figure 2.27, the code that manages the bypass valve is reported. It checks
the pressure difference between the valve inlet and outlet ports, and opens the valve
only if the inlet measure is higher than the outlet.
bypass.open = bypass.dp > 0;
Figure 2.27: Code that performs the control of the bypass valve.
2.12 Control valve

The transport network for natural gas does not have an immutable structure, but it
changes depending on the active import points and on the quantities injected and
extracted in the different areas.
To change the shape of the network, two types of valves are employed:
Gate valves are used to connect or isolate pieces of network. Their status can
be just either open or closed.
Control valves are operated to regulate either the pressure or the mass flow-
rate that passes through them. When they are completely closed, they act as
a closed gate valve.
In the model for the simulation and optimization, gate valves are not represented.
During the simplification process, these components are considered as open circuits
when closed, and as short circuits when open.
On the contrary, control valves are crucial to represent the real behaviour of the
network because they can regulate the state of the flow that passes through them.
2.12. CONTROL VALVE 47
For this reason, a branch of the network connected to the control valve outlet is
subjected to a flow considerably different to the one it would be affected if it was
directly connected to the rest of the network.
So, in this work, control valves are considered and modelled every time they connect
two branches of the network. On the contrary, they are neglected when they are
only used to reduce the pressure that exits from the network. In this case, they are
collapsed into the exit node.
There are two categories of control valves: the first one controls both the mass
flow-rate and the outlet pressure, the second one regulates only the outlet pressure.
2.12.1 Control valve with mass flow-rate and outlet pressure

set-points
This control valve is modelled with the code in Figure 2.28. Since this component is
short, the mass storage can be neglected, therefore the inlet and the outlet mass flow-
rates have always the same value. The ideal control is realized defining a curvilinear
coordinate s: when s > 0, the valve is controlling the flow-rate; when s < 0, it
regulates the outlet pressure. In addition, the control logic checks the pressure
difference on the valve: if it is greater than a saturation value, the valve can act
on the flow-rate; if not, w becomes a function of ∆p, corresponding to the valve
being fully opened and the control action being saturated. The different working
conditions are schematized in Table 2.11 and Table 2.12.
Flow-rate control Pressure control

s>0 s<0
∆p > ∆psat wSP wSP + s · wnom
wnom wnom
∆p < ∆psat · ∆p · ∆p + s · wnom
∆pmin ∆pmin
Table 2.11: Scheme of the possible different values for the w at the outlet.
Flow-rate control Pressure control

s>0 s<0
pmax − s · pnom pmax
Table 2.12: Scheme of the possible different values for the pout .
inlet.w = w ;
outlet.w = -w ;
outlet.p = pOut ;
dp = inlet.p - outlet.p ;
dpSat = dpMin * wSP / wNom ;
w = ( if dp > dpSat then wSP else wNom / dpMin * dp ) +

( if s > 0 then 0 else s * wNom ) ;
pOut + impedanceCoeff * w * pNom / wNom =
if s > 0 then pMax - s * pNom else pMax ;
assert ( s >= 0 or w >= 0 , " If pressure controller is engaged ,

flow - rate must be not negative " ) ;
Figure 2.28: Modelica code of the model for the control valve with set-points for
mass flow-rate and outlet pressure.
At first, the case where s > 0 is considered. The outlet pressure is not imposed
and its equation is used to evaluate the value of s. The mass flow-rate can be fixed
at either the value wSP , when the pressure difference is greater that the saturation
limit, or it can be proportional to ∆p. This condition emulates a simplified linear
characteristic curve for a valve, with the angular coefficient wnom/∆pmin . Then, the
case where s < 0 can be analysed. pout is fixed to its maximum value pmax . The
flow-rate is then used to evaluate s as dimensionless distance of the actual w from
wnom
wSP or ∆p min
· ∆p, depending on the value of ∆p.
The value of pOut is not directly impressed to the outlet port. In fact, doing so
could lead to a variable-index problem, not managed by Modelica. For this reason,
a fictitious impedance is placed in the equation and its value can be tuned acting
on impedanceCoeff.
2.12.2 Control valve with outlet pressure set-point

The purpose of this component is to set a maximum value for the pressure in the
following branch of network. In the code in Figure 2.29, the ideal control is imple-
mented using a curvilinear coordinate s. In this case, it links the variables ∆p and
pout . To prevent the issue with the solver caused by the variable-index problem,
even in this case, a small impedance is placed, whose value is tuned by the variable
impedanceCoeff.
When s > 0, the outlet pressure is limited by the control to the value pset and the
inlet pressure can be at any value grater than pset . In this case, s is obtained from
2.13. OTHER MODELS 49
inlet.p = pIn ;
outlet.p = pOut ;
Deltap = pIn - pOut ;
inlet.w = w ;
outlet.w = -w ;
Deltap - impedanceCoeff * w * DeltapNom / wNom =
if s > 0 then DeltapNom * ( 1 + s ) else DeltapNom ;
pOut + impedanceCoeff * w * pNom / wNom =
if s > 0 then pSet else pSet + s * pNom ;
Figure 2.29: Modelica code of the model for the control valve with outlet pressure
set-point.
the ∆p equation. On the contrary, when s < 0, the pressure is lower than the limit,
so the valve intervention is not required. A pressure loss ∆p, in the order of 0.1 bar,
is introduced to model the real behaviour of a fully-opened valve. The value of s is,
in this case, obtained from the pout equation.
2.13 Other models

2.13.1 Node
The Node (Fig. 2.30) is a fictional model: it has no physical meaning and it is used
only to simplify the writing of the network model performed by the Python script
(Sec. 3.6), since it is the Modelica counterpart of the graph’s node.
Figure 2.30: Icon for Node model.
It has only two equations (Fig. 2.31): the first one represents the value of mass
flow-rate that is spilled from the network in that point, while the second one is used
only during the initialization process. In fact, during the generation of the network
model, the pressure of the nodes is evaluated from the actual measures available.
Setting those start values for the node variable helps the solver.
2.13.2 Sources
The boundary conditions are used to set a point where the modelling process is
stopped and the rest of the environment is represented by a limited group of suitable
model Node
G a s N e t w o r k s . I n t e r f a c e s . P n e u m a t i c P o r t _ a port ;
Ty pe s. Ma ss Fl ow Ra te wcons ;
Types.AbsolutePressure p;
equation
port.w = wcons ;
port.p = p ;
end Node ;
Figure 2.31: Modelica code of the Node model.
Figure 2.32: Icons for the flow-rate source, the time-varying pressure source, and
the constant pressure source.
equations. In a problem that involves the flow of a fluid, the most common options
are to prescribe either the pressure or the flow-rate.
In the model of the gas network, these two options are used in different cases:
Prescribed pressure: it emulates the condition where the pressure at an import

point is guaranteed by the foreign vendor.
Prescribed mass flow-rate: the vendor does not guarantee a certain pressure
value, but it surely provides the gas quantity specified in the signed commercial
agreements.
These components are modelled using a pneumatic port, that connects the source
with the network, and a real input, that receives the signal that states the value of
the quantity to generate.
In particular, in the flow-rate source, the value of the input is assigned to the flow-
rate variable of the port, as described in Figure 2.33. At line 5, the minus sign
is used to respect the convention stating that a flow entering in a port must be
positive. In this case, the user expects a w flow-rate that exits from the port.
The same approach is followed in the model of the time-varying pressure source.
In addition, a constant pressure source is available in the GasNetworks package. In

this case, the value of the pressure is not received through an input port, but it is
2.13. OTHER MODELS 51
1 model FlowSource
2 I n t e r f a c e s . P n e u m a t i c P o r t _ b port ;
3 Modelica.Blocks.Interfaces.RealInput w;
4 equation
5 port.w = -w ;
6 end FlowSource ;
Figure 2.33: Modelica code of the FlowSource model.
directly declared using a parameter p.
2.13.3 System
The Modelica inner/outer construct is used to describe the environment in which
the components are placed. In a model M1, in the section where the variables are
declared, is possible to place also a model M2 marked as outer. Then, when M1 is
used in another model, for example a pipe is placed in the network, it connects to
the component of the network declared as inner and inherits its content.
In the GasNetworks package, the System model is used as container for the envir-
onment parameters. In fact it contains:
Gas temperature;
Number of finite volumes to split every pipe in;
Default initialization option chosen for the pipes in the network;
Boolean flag that activates the evaluation of the total pressure profile in the
pipes using Bernoulli’s principle.
Finally, the System contains the interface PowerMeteringPort (Sec. 2.5) that is
connected to every gas turbine in the network and it is used to evaluate the total
thermal power consumption.
2.13.4 System initializer

The initialization of the models that contain states is a critical issue. In Section
2.6.3.1, the different approaches to initialize a pipe are described. The steady-state
option is the most appealing one because it provides a starting point where the
component is at its equilibrium and so it avoids the relaxation swings that may
happen if the component is initialized in a non-optimal condition.
Sys. Init.
Figure 2.34: Icon for SystemInitializer model.
In some cases, this method leads to a singular initialization problem that arises at
system level. In fact, if none of components has an explicit initial value for the
pressure, the total mass of gas stored inside the pipes is unknown and cannot be
evaluated using the steady-state equations.
To solve this problem, a component as the one described in [4] is modelled in the
package. Its first equation sets that the mass flow-rate through its port is null, and
the second one imposes the desired initial value of the pressure on its port. This is
sufficient to remove the singularity in a system.
model System Initi alizer

parameter T y p e s . A b s o l u t e P r e s s u r e p_start ;
final parameter Ty pe s. Ma ss Fl ow Ra te w_b = 0;
G a s N e t w o r k s . I n t e r f a c e s . P n e u m a t i c P o r t _ a inletGasPort ;
equation
0 = inletGasPort.w + w_b ;
initial equation
inletGasPort.p = p_start ;
end S ystemI nitial izer ;
Figure 2.35: Modelica code for the SystemInitializer model.

Chapter 3
Modelling of the Italian Gas

Network
The Modelica package described in the previous chapter is employed to build the
model of the Italian network for gas transport.
3.1 Source files

To build the model of the Italian network, several different sources are required.
Network description XML file The first considered file contains the geometrical
and geographical data for every element. It is an XML file and its two main sections
describe the nodes and the edges of the network.
In Figure 3.1, all the parameters that can be extracted from the file are reported.
The Node class contains information about the points where two or more different
elements come together. The parameters are:
id and name are the univocal identifiers for the node. The first one is numer-
ical, the second one is alphanumeric;
x, y, and height contain the geographical coordinates of the node;
supply takes the value “1” when the node represents a source of gas for the
network.
53
54 CHAPTER 3. MODELLING OF THE ITALIAN GAS NETWORK
Figure 3.1: UML diagram describing the content of the XML file. In bold, the
parameters that are then imported into the graph are highlighted.
< NODE alias = " " height = " 70.0 " x = " 67 82999. 999926 411 " y = "
5942 999.99 993552 5 " name = " VA_601O " id = " 17875 " / >
Figure 3.2: Example of the definition of a node in the XML file.
< ELEMENT POMax = " 75 " pwm = " DIAMETER " Rout_unit = " 1 " dm = " base " pwp = " 1.0
" PIMin_unit = " barg " name = " POGGIO_RENAT " Rin_unit = " 1 " diameter = "
1378.40002 " subsystem = " 23698 " sizeFactor = " 35.0 " diameter_unit = "
mm " Rin = " 0.00000 " PIMin = " 37 " type = " compressor station " id = " 6 "
Rout = " 0.00000 " POMax_unit = " barg " >
< NODE symbol = " 1 " id = " 5777 " connection = " 0 " / >
< NODE id = " 5751 " connection = " 1 " / >
</ ELEMENT >
Figure 3.3: Example of the definition of an element in the XML file. Every element
contains two sub-items that define its terminal nodes.
3.1. SOURCE FILES 55
Figure 3.4: UML diagram describing the content of the “Mercato” CSV file.
" SERGNANO - CALUSCO "; 05/10/2020 15:28:00; 275670

Figure 3.5: Lines extracted from the “Mercato” CSV file.
The Element class contains a long list of objects, the most relevant are:
type: it is a string that shows which component is modelled with that element.
It can be pipe, compressor station, valve, control valve, non-return valve, joint,
or measuring station;
id and name identify the element univocally. The first one is an integer num-
ber, the second one is a string;
diameter, length, and roughness describe the geometrical features of the ele-
ment considered.
Flow-rates description CSV file The second file analysed is in Comma-Separated

Values (CSV) format. It is called Mercato, and contains the time series of the mass
flow-rate in every point, representing an entry or exit point for the network.
As described in Figures 3.4 and 3.5, every row of this file contains
An alphanumeric identifier for the element;
A string that contains the date and time information of the datum;
The value of the flow-rate expressed in Sm3 /h (Sec. 2.3).

I - TARS_CENT - ASPIRA__ - P |61.400| bar | |2 077 |7 5. 00 |7 4. 00 | ND |45.00| no

I - TARS_CENT - M3______ - STK | 00 0 0 00 0 0 00 0 0 01 0 0 || |2106|||||
Figure 3.6: Two lines extracted from an “NSI” file.
Figure 3.7: UML diagram describing the content of the “NSI” file.
Components status NSI file The next considered file is called NSI and contains
the logs generated by the SCADA system. These data are acquired with a time
interval of four minutes and then saved into a new file. Its structure, described in
Figure 3.7, is not strictly organized, as can be seen in Figure 3.6. In fact, in a row,
one can find:
A physical measurement, composed of its value and unit of measure;
A status, described by a binary number that can be interpreted using a key.
These elements have two common attributes: a numeric identifier an alphanumeric

code. In the latter, three features are saved in a “dash-separated value” format.
The first block represents the location of the machine, the second one the name of
its component and the third one the quantity measured. For example, in Figure 3.6,
the first line contains a pressure measurement at the inlet of the station in Tarsia
(CS). Differently, the second one is the status of the third machine in the station in
Tarsia. 2077 and 2106 are the numerical identifiers for these elements.
Aggregate file In the network model provided by Snam as the source for this
work, some complex configurations of valves, connected in parallel or in series, are
aggregated into a single and equivalent component. The statuses of these elements
are contained in the Aggregate CSV file. It is organized into three columns:
3.1. SOURCE FILES 57
" S_I - P A L M _ T E R M _ V 3 2 B _ V 2 2 B _ V 1 2 B ";03/03/2021 11:40:00;0

Figure 3.8: Example of rows extracted from the Aggregate file.
Figure 3.9: UML diagram describing the content of the “Aggregate” file.
The ID of the aggregated element;
The string that contains the date and time of the datum;
The saved status.
Figure 3.9 depicts the scheme of the file, while Figure 3.8 shows an example of some
lines extracted from it.
Associations REQ file The last file is used to match the nodes and the elements
extracted from the XML file with the data from the Mercato and NSI files. It is a
plain text file with a structure that can be called “space-separated values”, and its
extension is .REQ.
Figure 3.10: UML diagram describing the content of the REQ file.
CV VDC3 . SPO NSI #I - FAVA_TERM # REGOIONE # SP #1133340

NO TARVISIO_PRE . Q MERCATO # TARVISIO
CS TARSIA . TAMB NSI #I - TARS_CENT #I - TARS_CENT - CENTRALE - Tamb # CENTRALE_T
-12# Tamb #2087
Figure 3.11: Examples of rows from the REQ file.
As depicted in Figure 3.10, the file is composed of elements that contain three
attributes:
in the type, two characters are saved that associate the element to a category
of components.
– For nodes, they can be NS, representing points where the gas, usually,
enters the network, or NO, locations where the gas is extracted from the
network.
– For edges, they can be CS (compressor station), CV (control valve), MS

(metering station), or VA (valve).
the id is an alphanumeric string that identifies the element. It is composed

of two parts divided by a dot. The former is the element name, the latter
indicates which feature of the element is described;
the data attribute contains a composite string that can be split at every #
symbol.
The first part of the id string corresponds to the element name in the XML file, so
this represents the first link between the two files.
The number of tokens in the data attribute and their content can vary with the type
of the element, but the first one always links the element with either Mercato or NSI
files. When the element contains the link with an NSI record, the last token cor-
responds to the integer id of that record. Instead, when it describes the connection
with a Mercato element, the last token contains the string ID_ANAGRAFICA_MERCATO.
Finally, when it represents the link to Aggregate file, the last token is a string that
matches the OLD_ID column.
During the first part of the development of the script that manages these files, the
content of the XML and REQ files was stored into an XLSX file, organized in three
sheets: nodes and edges, that correspond to the NODE and ELEMENT part of the XML
file, and associations , that represents the REQ file. For this reason, in the code,
3.2. PYTHON LIBRARIES 59
there are several references to these names.
Using these relationships, it is possible to merge all the data coming from the con-
sidered files into a single container, well-organized and easier to access. The proced-
ure is described in section 3.3.2.
Machines datasheets Some files are not directly involved in the generation of the
gas network, but they are fundamental for the description of its functioning: they
are the datasheets of the machines, compressors and gas turbines. These files are
provided as PDF, generated by the program SIRE2000, and contain the numerical
coefficients for the polynomials described in Sections 2.7 and 2.8 and the machines
maps. This information is used during the instantiation of the compressor stations
models in the Modelica package. To automate this procedure, the numerical coeffi-
cients have been manually copied into two Excel spreadsheets. The first one contains
information about the compressors, the other one data about the gas turbines. They
are organized with one row per model of machine used in the network, and every
column contains one parameter that is used to model it.
3.2 Python libraries

Following the scheme in Figure 3.12, the libraries used in the different sections of
the program are described.
In the “Data import and cleaning” step, pandas and BeautifulSoup are employed.
Pandas [20] is a data analysis library, open-source since 2009. At first, it has been
chosen for the simple procedure that it offers to import CSV and Excel files. With
a one-line command, pandas is able to import the information contained in the
file and to organize it in a flexible and powerful data structure called DataFrame.
In a DataFrame is possible to apply operations to every value in a column in a
straightforward way, as represented in line 4 of Figure 3.13, or to keep only the rows
that satisfy a desired condition, as in line 5.
BeautifulSoup [32] is an HTML and XML parser library and its purpose is to make
the structure and the data contained in the file accessible to Python. Unfortunately,
no function in the library directly converts the XML file into a pandas DataFrame.
For this reason, a custom function called dataframeFromXMLItems had been de-
veloped. As reported in line 7 of Figure 3.14, the function returns a DataFrame and
takes as input:
Figure 3.12: Block diagram representing the procedure that elaborates the data
acquired, leading to the geographical map and the Modelica model.
1 import pandas as pd
2 nodi = pd . read_excel ( absolute_path_excel , engine = ' openpyxl '
, sheet_name = ' NODI ' )
3 mercato = pd . read_csv ( absolute_path_mercato , sep = ' ; ' ,
encoding = ' cp1252 ' , decimal = ' , ' )
4 mercato [ ' NR_VALORE ' ] = mercato [ ' NR_VALORE ' ] *
c o nv e r si o n Sm 3 h to K g s
5 mercato = mercato [ mercato [ ' D T _D A T A_ R I FE R I ME N T O ' ]==
selectTime ]
Figure 3.13: Examples of pandas instructions used in the program.

3.2. PYTHON LIBRARIES 61
1 from bs4 import BeautifulSoup as bs

2 with open ( xmlFilePath ) as fp :
3 xmlFile = bs ( fp , ' lxml ' )
4
5 items = xmlFile . nodes . find_all ( ' node ' )
6 attrList = [ ' height ' , ' x ' , ' y ' , ' name ' , ' id ' , ' supply ' ]
7 nodesDF = d a t a f r a m e F r o m X M L I t e m s ( items , attrList , ' nodes ' )
Figure 3.14: Examples of BeautifulSoup instructions used in the program.
The BeautifulSoup variable containing all the rows defined with tag node;
The list of attributes that must be imported from the XML file;
The category of the elements to import.
In the “Data union” phase, another function from pandas has been employed: merge.
Merge is a powerful function that, given two DataFrames, creates a new one that
contains the information coming from both inputs. The algorithm that performs
this operation is able to check if in two columns (one per DataFrame) there are rows
that have the same value: if so, it writes the information from both tables in the
same line in the output.
The pandas.merge function has an important parameter, that can be specified, that
changes the behaviour when the algorithm encounters a row in a table that has no
matches in the other one. To explain the possibilities, consider the tables in Figure
3.15 and the following command:
result = pd . merge { tableA , tableB , on = ' Name ' , how = ' inner ' }
The parameter how can be set to different methods:
inner: only the rows that are in both tables are written in the result;
outer: no rows are dropped;
left: it keeps all the rows of the first table, dropping the elements of the
second one that are not matched;
right: it executes the same procedure as left, swapping the first and the
second tables.
In the program, the left procedure has been used, as exampled in Figure 3.16, be-
Table A Table B
Name Value A Name Value B
Anna 15 Anna 8
Carl 4 Camille 12
John 33 John 21
Mary 27 Ted 18
Outer merge Inner merge

Name Value A Value B Name Value A Value B
Anna 15 8 Anna 15 8
Camille 12 John 33 21
Carl 4
John 33 21
Mary 27
Ted 18
Left merge Right merge

Name Value A Value B Name Value A Value B
Anna 15 8 Anna 15 8
Carl 4 Camille 12
John 33 21 John 33 21
Mary 27 Ted 18
Figure 3.15: Examples of the application of the merge algorithm with different
methods.
mergeNodes1 = pd . merge ( nodi , associazioniNodes , left_on = '

Attribute : name ' , right_on = ' ELEMENTO ' , how = ' left ' )
mergeNodes2 = pd . merge ( mergeNodes1 , mercato , left_on = '
MISURA ' , right_on = ' I D _ A N A G R A F I C A _ M E R C A T O ' , how = ' left ' )
dfNodes = pd . merge ( mergeNodes2 , nsi , left_on = ' MISURA ' ,
right_on = ' col5 ' , how = ' left ' )
Figure 3.16: Example of the use of pandas.merge in the program.

3.3. DATA ELABORATION 63
cause the first DataFrame, called nodi, contains the information about the network.
For this reason, every row in it must be preserved during the merging procedure.
In “Graph generation” and “Graph simplification” phases, the most widespread

package for the creation and manipulation of complex networks has been used:
NetworkX [8]. With this library, it is possible to represent the gas network as a
graph composed of nodes and edges. Edges must have two nodes as ends and can
join each other only through nodes. Both nodes and edges can be used to store
data inside their attributes. The main advantage offered by this solution is the
coherence of the stored information: in fact, since the nodes and the edges must
be created using the dedicated functions, the library automatically creates all the
needed relationships between the nodes and edges involved. In addition, it has some
useful tools to analyse the graph and to, for example, understand if two nodes are
connected to each other and to measure their distance.
The last package considered is used in the “Map export” phase. It is called Folium
[37], it is based on leaflet.js, and its purpose is to show Python data on maps that
can be opened in a web browser. Folium provides commands to add points on the
map with customizable markers, draw lines and shapes, add informative pop-ups on
these elements, and create heat-maps and other visualizations. In pop-ups, images,
plots and HTML formatted texts can be included, so a single webpage can represent
not only the shape of the network, but also the numerical data linked with it.
3.3 Data elaboration

The starting point for the building procedure is the information contained in the
files described in Section 3.1. To elaborate it, a Python program has been developed.
Python is an open-source, high-level, interpreted, and one of the most popular pro-
gramming languages. Thanks to this, it has an extensive catalogue of libraries that
covers a variety of fields of application, especially in the mathematical and scientific
area.
3.3.1 Data import and cleaning

The first set of operations that the program executes are related to the sources of
the information used to reconstruct the network and its behaviour.
At first, the script loads the XML file using the custom function based on Beautiful
Soup: importDataFromXML. It opens the file from the file system, and creates the
index of elements contained in the nodes section. Then, another custom sub-routine
parses every item copying its attributes into a DataFrame. Finally, the content of
every column that stores numerical data is converted from string into integer or float
type.
Later, importDataFromXML applies a similar procedure to the elements section.

The main difference in this part is that every element contains two node sub-items.
These objects embody the information about the connections of the element.
The following files, Mercato, NSI, and Aggregate, being in CSV format, are directly
loaded as DataFrames using the pandas function read csv.
The cleaning phase consists in preparing the data acquired for the subsequent op-
erations.
The proceeding starts with the removal of the nodes that are not actual nodes of
the network. In fact, the XML file contains, besides the actual components of the
networks, about 30 graphical elements. These objects are saved as joints, with the
respective nodes, and are used to contain the legend of the map in the program used
by Snam. So, in the script, these items are discarded and never represented.
The second block of operations is performed on mercato dataframe. Measures of

flow-rates are converted from Sm3 /h to kg/s using the conversion procedure ex-
plained in Section 2.3 and considering the average composition of the gas for the
considered area, the column containing strings that represent date and time are
transposed into timestamps, and the spaces in the names are swapped with under-
scores. Furthermore, to avoid performance degradation, the dataframe is reduced
keeping only the data relative to the time interval of the simulation. In fact, in some
cases, the CSV file contains records that cover an entire year. So, applying an initial
filter can speed up the program, avoiding useless operations.
Then, the pressure measurements contained in the NSI dataframe are considered.
In fact, the majority of them are expressed in barg. To be compatible with the
compressor’s equations, these values are converted into absolute bar.
Within the associazioni dataframe, all the rows with an empty cell in the SORGENTE
column are dropped, the MISURA column is created extracting the last token from
the data string, as explained in Section 3.1, and the remaining rows are split into
two dataframes. The first table contains the information related to the nodes,
and it is called associazioniNodes, the second one refers to edges, so its name
is associazioniEdges.
3.3. DATA ELABORATION 65
Figure 3.17: Diagram of the connections between items in different source files used
when merging data about nodes.
Finally, the aggregate dataframe is prepared by converting the strings in the

DT_DATA_ORA_RIFERIMENTO column in timestamps, and then, it is filtered keeping
only the rows with the desired time values.
3.3.2 Data merge

The merging phase is, from the point of view of the Python code, straightforward.
With a single line of code, the content of two different dataframes is copied into a
new table, where the information about an element is organized in a single row. The
operating principle of this algorithm is explained in Section 3.2.
The search of the common attributes to use as an index for the merging procedure is
more relevant. In Section 3.1, it was explained that the information about nodes is
obtained by analysing four files: XML, REQ, Mercato, and NSI. Starting from one
element extracted from the XML file, it is possible to follow the merge procedure,
as shown in Figure 3.17.
Consider an element of the XML file. It contains the geographical location of

the node, its identifiers and the tag that indicates that it is a supply node.
< NODE alias ="" height ="100.0" x = " 6 7 3 7 0 0 0 . 0 0 0 2 5 3 4 1 7 5 " y

= " 58 5 1 00 0 . 00 0 2 20 0 9 " name =" PANIGAGLIA " id ="4963" supply ="1" >
Look for the rows in the REQ file that has, as id, the same name string from
the XML element. There can be more than one result because one node could
be associated with more physical quantities. For example, in this case, there
are a flow-rate measurement, that can be found in the Mercato file, and a
pressure reading, linked to an NSI entry.
Figure 3.18: Diagram of the connections between items in different source files used
when merging data about edges.
1 NS PANIGAGLIA . Q MERCATO # PANIGAGLIA

2 NS PANIGAGLIA . PM NSI #I - SPEZI_SRG # CORTE # P # 1190022
Match the last token of the string in the last column of row 1 with the string
in the first column of the Mercato file. This procedure returns a series of rows
that complement the node with its measurements representing the flow-rate
that either enters or leaves the network.
" PANIGAGLIA " ;05/10/2020 13:36:00;372116

" PANIGAGLIA " ;05/10/2020 13:40:00;372517
" PANIGAGLIA " ;05/10/2020 13:44:00;373108
Finally, the procedure looks for the last token of the last string in row 2 in
the NSI file. The result is always unique, and, in this case, it represents a
pressure measurement.
I - SPEZI_SRG - CORTE___ - P |59.426| bar | | 1190022 |69.00|68.00| ND

|44.00| no
When considering the edges of the network, the method to aggregate the information
is very similar to the node’s one. The only difference is that there are no entries in
the Mercato file for edges, but the Aggregate table takes its place. In Figure 3.18,
this procedure is summarized, while in Figure 3.19 the Python code that performs
the actions is reported.
3.4. GRAPH HANDLING 67
# Merging dataframes for Nodi , Associazioni , Mercato and NSI

mergeNodes1 = pd . merge ( nodi , associazioniNodes , left_on = ' Attribute :
name ' , right_on = ' ELEMENTO ' , how = ' left ' )
mergeNodes2 = pd . merge ( mergeNodes1 , mercato , left_on = ' MISURA ' ,
right_on = ' I D _ A N A G R A F I C A _ M E R C A T O ' , how = ' left ' )
dfNodes = pd . merge ( mergeNodes2 , nsi , left_on = ' MISURA ' , right_on = '
col5 ' , how = ' left ' )
# Merging dataframes for Binari , Associazioni and NSI

mergeEdges1 = pd . merge ( binari , associazioniEdges , left_on = '
Attribute : name ' , right_on = ' ELEMENTO ' , how = ' left ' )
dfEdges = pd . merge ( mergeEdges1 , nsi , left_on = ' MISURA ' , right_on = '
col5 ' , how = ' left ' )
Figure 3.19: Python code that performs the merge operation for both nodes and
edges.
3.4 Graph handling
3.4.1 Graph generation

In Section 3.3.2, the data about nodes and edges have been aggregated and organized
in dataframes. This structure is powerful and flexible for data storage, but it is not
ideal to represent the connections and the relationships between the components.
In addition, the use of tables to store information does not guarantee the integrity
of the network. For example, if a node is deleted from its dataframe, all the edges
connected to it become “orphans”.
For this reason, NetworkX is used to build, elaborate and inspect the network. The
first step, described in this section, is the generation of the graph.
The first, fundamental, assumption is to consider the network topology and the gas
flow directions as fixed during the simulation or optimization time intervals. This
allows to generate the network graph using the information available in the instant
considered as the initial point.
The procedure starts with the declaration of the variable that contains the graph:
G = nx . Graph ()
This command creates a non-directed, simple graph. This type is used because the
gas flow inside the pipes is not known a-priori, and because in the Snam network
there are no multiple edges between the same pair of nodes. In fact, a simple graph
can not contain more than one edge between a considered pair of nodes. If the
program tries to add a new edge where another one is already present, NetworkX
does not raise any exception and it simply overwrites the old item with the new
one. Instead, in a non-directed MultiGraph, two nodes can be connected by more
than one edge. But, for this reason, the identifier of each edge requires an additional
information, called key, to distinguish between the multiple instances with the same
boundary nodes. This additional hurdle would have obstruct the development of
the simplification procedures, and so the MultiGraph is employed only at the end
of that phase.
For the asymmetrical components, e.g. compressor stations, the inlet node ID is
stored in an auxiliary variable to keep track of it.
Then, the program instantiates the nodes of the network. To do so, it iterates on
every row of the nodes dataframe, copying every attribute into the graph. Since
in the dataframe there could be more lines for the same node, containing different
attributes for it, the algorithm must check if an instance with the same id is already
present in the graph and merge the data.
In the code fragment in Figure 3.20, the try - except structure is used for this
purpose. The instruction at line 4 tries to load the node in the graph using the id
of the current dataframe line. If this node is not present, the instruction raises a
KeyError exception, that is caught by the statement at line 19. So the script reads
the attributes from the dataframe and, at line 38, creates the node. Instead, if the
node is found, the execution continues from line 5, and the routine adds the missing
information to the existing node.
In this step, a new attribute is introduced: tpMarket. In this variable of type list,
a tuple, containing the node identifier in the Mercato file and its category, is saved
and it will be exploited in the process that writes the Modelica model, explained in
Section 3.6.
When every row of the nodes dataframe has been copied into the graph, the pro-
gram can start to create the edges. Since there are several different types of edges,
the dataframe that contains their information, called dfEdges, is split into smaller
tables, each containing only one category.
The conceptual procedure is very similar to the one dealing with the nodes, and
it does not present relevant differences between the categories. In fact, after filter-
ing dfEdges preserving only the desired type of edges, the program performs the
1 for index , row in dfNodes . iterrows () :

2 try :
3 # Check if the node is already present
4 dummy = G . nodes [ row [ ' Attribute : id ' ]]
5 nodeId = row [ ' Attribute : id ' ]
6
7 # If so , add only the data that are missing
8 if G . nodes [ nodeId ][ ' P ' ] == ' None ' and row [ ' col3 ' ]== ' bar ' :
9 G . nodes [ nodeId ][ ' P ' ] = row [ ' col2 ' ]
10
11 if G . nodes [ nodeId ][ ' T ' ] == ' None ' and row [ ' col3 ' ]== ' ◦ C ' :
12 G . nodes [ nodeId ][ ' T ' ] = row [ ' col2 ' ]
13
14 if row [ ' SORGENTE ' ] == ' MERCATO ' :
15 if str ( row [ ' NR_VALORE ' ]) == ' nan ' :
16 row [ ' NR_VALORE ' ] = 0
17 G . nodes [ nodeId ][ ' market ' ] = row [ ' NR_VALORE ' ]
18 G . nodes [ nodeId ][ ' tpMarket ' ] = [( row [ ' MISURA ' ] , row [ '
Column1 ' ]) ]
19
20 except KeyError :
21 category = ' None '
22 pressure = ' None '
23 temperature = ' None '
24 misura = ' None '
25 tpMarket = []
26 if str ( row [ ' Column1 ' ]) != ' nan ' :
27 category = row [ ' Column1 ' ]
28
29 if str ( row [ ' NR_VALORE ' ]) == ' nan ' :
30 row [ ' NR_VALORE ' ] = 0
31
32 if str ( row [ ' MISURA ' ]) != ' nan ' :
33 misura = row [ ' MISURA ' ]
34
35 if row [ ' col3 ' ]== ' bar ' :
36 pressure = row [ ' col2 ' ]
37 elif row [ ' col3 ' ]== ' ◦ C ' :
38 temperature = row [ ' col2 ' ]
39
40 if row [ ' SORGENTE ' ] == ' MERCATO ' :
41 tpMarket = [( misura , category ) ]
42
43 G . add_node ( row [ ' Attribute : id ' ] , height = row [ ' Attribute :
height ' ] , x = row [ ' Attribute : x ' ] , y = row [ ' Attribute : y ' ] ,
name = row [ ' Attribute : name ' ] , supply = row [ ' Attribute : supply
' ] , market = row [ ' NR_VALORE ' ] , category = category , misura =
misura , P = pressure , T = temperature , tpMarket = tpMarket )
Figure 3.20: Excerpt of the program code that creates the nodes in the graph and
saves into them their attributes, read from the dataframe.
following steps:
Reads and stores into temporary variables all the attributes of the edge, using
the try - except construct when the information can be missing.
# Read ambient temperature and conversion from ◦ C to K

TAMB = pd . Series ( attrList [ attrList [ ' CARATT ' ]== ' TAMB ' ][ ' col2 ' ]) .
astype ( ' float ' ) . values [0] + 273.15
# Read outlet pressure set - point in bar
SPO = pd . Series ( attrList [ attrList [ ' CARATT ' ]== ' SPO ' ][ ' col2 ' ]) .
values [0]
# Read mass flow - rate set - point and conversion from Sm3 / h to kg
/s
SM = pd . Series ( attrList [ attrList [ ' CARATT ' ]== ' SM ' ][ ' col2 ' ]) . astype
( ' float ' ) . values [0] * c o n ve r s io n S m3 h t oK g s
Applies the correct unit conversion to the numerical attributes, checking the
unit of measure in the related column, if necessary. In the graph, length,
diameter and roughness are expressed in meters; flow-rate set-point and meas-
urement in kg/s; pressures set-points in bar; temperatures in Celsius degrees.
# Read the length and convert it from m to km , since it is the

most common case
length = pd . Series ( attrList [ attrList [ ' Attribute : subsystem ' ] ==
23698][ ' Attribute : length ' ]) . values [0]*1000
if pd . Series ( attrList [ attrList [ ' Attribute : subsystem ' ] ==
23698][ ' Attribute : length_unit ' ]) . values [0] == ' m ' :
# Manages the lines where the length is expressed in meters
length = length / 1000
Creates the edge in the graph, saving in it every attribute found.
G . add_edge ( node1 , node2 , edgeid = edgeid , name = name , category =

category , diameter = diameter , POMax = POMax , PIMin = PIMin ,
Tout = TO , Tamb = TAMB , SPI = SPI , SPO = SPO , MM = MM , SM = SM , MODE =
MODE , roughness = roughness , length = length )
At this point, the graph that represents the data acquired from the source files is
complete. All the nodes and edges are saved, and their attributes constitute the
gas network status at the considered time instant. The Figure 3.21 represents the
UML diagram of the graph. In the “Network Graph” class the methods used in the
Figure 3.21: UML diagram of the graph.
program are reported. They are the actions that the code performs on the graph.
The other classes are used to describe nodes and edges. Every edge is specified using
a particular class that contains its own attributes, but every edge class inherits some
common features from the superclass “Edge”.
3.4.2 Graph simplification

3.4.2.1 Description of the problem
In Figure 3.22, it is evident that the complexity of the graph is too high to achieve
an efficient simulation and optimization workflow.
This is due to several causes:

Figure 3.22: Detail of the map representing the original graph around the Enna
pumping station. Orange circles represent the NO nodes, green ones are NS nodes
and grey shapes are uncategorized ones. Yellow lines represent joints, blue ones are
valves, black segments are pipes, red ones stand for non-return valves, and purple
lines are compressor stations. A dashed blue line means that the valve is closed.
There is a large number of joints: they represent more than 2500 elements on
a total of 6800 edges, but they do not have a physical meaning. Joints are
fictitious components that are used to represent the graph on the map with a
clearer layout, and their meaning is to make a short-circuit connection between
their ends. For this reason, from the point of view of resources optimization,
the best approach is to collapse their nodes into one and delete them from the
network.
Long pipes are split into smaller chunks using valves, that, in normal condi-
tions, are fully opened and have a negligible pressure loss. Their purpose is
to guarantee safety in case of faults or to allow closing the section when some
kinds of maintenance work are carried out. Since the program should analyse
the network at the start of the day and consider it as immutable for the rest
of the 24-hour simulation, these valves can be simplified. The solution chosen
for this work is to consider each open valve as a short circuit, and to remove
them collapsing the ends into one. On the contrary, closed valves are removed
leaving the two nodes disconnected.
Pipes are often divided into shorter segments, even if there are no other ele-
ments in between. This happens for modelling reasons related to the tool used
in Snam to represent the network. For simplicity’s sake, these segments are
merged into a single tract, if their attributes like roughness and diameter are
the same. Since this procedure is carried out when the valves have been sim-
plified, the resulting pipes are interrupted only by nodes where multiple pipes
intersect or when they are connected to compressor stations or control valves.
Compressor stations and the elements that surround them are complex struc-
tures. To connect a pumping unit to the pipes, dozens of valves and check
valves are employed. The purpose of valves is to route the low-pressure gas to
the inlet of the compressor, and to direct the high-pressure flow towards the
outgoing pipes, while non-return ones are placed only for safety reasons. In this
work, the configuration of open and closed valves is considered as fixed, and,
neither during the simulation nor during the optimization, the flow direction
can be changed. Moreover, the pressure drop on a fully opened valve is negli-
gible. For these reasons, each station has been simplified directly connecting
its ports to the corresponding pipes.
Metering stations are elements that contain a group of sensors used to monitor
the quality of the gas flowing in the network, and there are about seventy of
them. They can not perform any action on the gas flowing in them, and
their pressure drop can be neglected. So, they are removed from the network,
collapsing their ends into a single node.
Control valves are elements that manage the properties of the gas flowing
through them. In the network, there are two different types of them, that can
be distinguished based on the quantity they control:
– Maximum pressure control: this valve monitors the gas pressure and
intervenes when it reaches a set-point. Its action is to prevent the outlet
pressure from going above its set-point (maximum value).
– Flow-rate control: this valve regulates the mass flow-rate of gas that
flows through it to match a given set-point. It also has a second set-point
that puts a maximum value for the outlet pressure, overriding the mass
flow-rate control when the maximum outlet pressure is reached.
These components are important, and so they must be preserved, when they
are used to regulate the gas flow in sections of the network that are simulated.
Instead, they can be neglected when they are employed to reduce the pressure
on an exit point, because, being their action visible only on the outlet, it does
not affect the simulated system. One example of the described situation is
represented in Fig. 3.23.
3.4.2.2 Adopted solutions
To create a simplified graph, considering all the assumptions made in Section 3.4.2.1,
a series of procedures has been developed in the code. The actions are executed in
the following order:
1. Remove the closed valves: these elements are recognized reading their MODE
attribute.
2. Remove the joints;
3. Remove the metering stations;
4. Remove the remaining valves;
5. Copy every edge that is not a pipe into a new multigraph;
6. Merge the series of small pipes into a single pipe, creating it in the new mul-
tigraph;
Figure 3.23: Examples of different placements of control valves. Arrow A points at

a control valve put inside the network, B indicates a control valve used to regulate
the pressure on an exit point, marked with C. In Section 3.5.2 the colour code for
the edges is explained.
7. Remove every remaining joint from the multigraph;
8. Remove the control valves that have a free end.
In steps 2, 3, 4, and 7, since the procedure is the same, a function called edge removal
is used, and it is described in the next paragraph.
edge removal function The edge removal function is used to remove every edge
of a specified category from the network graph and Figure 3.24 represents the steps
performed in the procedure. When an element is deleted, one of its nodes is cancelled
too. To avoid information loss, the attributes of the two ends are compared and
merged properly.
The function is called in the main program with a command as the following ex-
ample:
fn . edge_removal (G , ' joint ' , debug_flag )
where G is the variable containing the graph, ’joint’ is the category to remove,
and debug_flag is a boolean flag used to enable some optional printouts, useful to
check the algorithm behaviour.
Figure 3.24: Flowchart representing the sub-functions contained in edge removal

and their exit conditions.
1 for nd in network_graph . nodes :

2 if len ( network_graph . edges ( nd ) ) == 1 and network_graph . edges [
list ( network_graph . edges ( nd ) ) [0]][ ' category ' ] ==
edge_category :
3 node_to_keep = list ( network_graph . neighbors ( nd ) ) [0]
4 merge_adj_nodes ( network_graph , node_to_keep , nd )
5 removeNodes . append ( nd )
6 nod e_list _remov al ( removeNodes , network_graph )
Figure 3.25: Code that performs the removal of every edge that has one end not
connected to another element.
1 def i nt_edg es_rem oval ( network_graph , edge_category , debug_flag ) :

2 i n t _ e d g e _ s i m p l i f i c a t i o n ( network_graph , edge_category ,0 ,1 ,
debug_flag )
3 i n t _ e d g e _ s i m p l i f i c a t i o n ( network_graph , edge_category ,1 ,0 ,
debug_flag )
4 return
Figure 3.26: Code of the function int edges removal. Notice that the indexes in line
3 are swapped with respect to line 2.
The first sub-function of edge removal is called ext edges removal, and it is used to
take out every edge that has one end that is not connected to another element.
To achieve this goal, the algorithm iterates on every node and it does two checks:
that the target nodes are connected to only one edge, and that these edges’ category
is the desired one. If both conditions are satisfied, the other end of the edge is marked
as node_to_keep, and then the function that appropriately merges the attributes
is executed to allow deleting the first node without information loss. The actual
dismissal of the node is performed at the end of the loop to prevent the breaking of
the loop index. This step is repeated until the number of removed nodes is zero.
Then, the procedure executes the sub-function called int edges removal. Its beha-
viour is simple: it calls the function int edge simplification two times, switching the
indexes between the two executions. This is necessary because of some issues that
are analysed later in this paragraph.
In this routine, every edge of the considered category is analyzed following these
steps:
1. Consider the two end nodes of the edge: when the edge in between will be
removed, one of the two must disappear. It is called node_to_discard, the
other one is node_to_keep.
2. Check every edge that has node_to_discard as an end: if there is another

one with the same category of the edge to remove, skip this loop iteration and
start over with a new edge. If not, the algorithm continues.
3. The attributes of node_to_discard are merged into node_to_keep. The

primary feature that is preserved is the market value. In fact, the mass flow-
rates are summed, considering the sign convention too. The other performed
operations are explained in Section 3.4.2.2.
4. The edges that have node_to_discard as an end are copied changing it with
node_to_keep, managing the possible conflicts with the function safe add edge.
This procedure is described in Section 3.4.2.2.
5. node_to_discard is added to the list of nodes to remove from the graph.
At the end of the loop, a simple function takes the list of nodes to remove and
deletes them from the network graph. NetworkX automatically removes every edge
that has those nodes as an end.
The check performed in step 2 is fundamental to avoid errors when multiple ele-
ments of the same categories are placed one next to the other. This can be better
understood considering the following example, where the above-mentioned check is
skipped.
Taking as reference the scheme in Figure 3.27 and considering removing blue edges,
it can be seen that the (1,2) and (2,3) belong to the same category.
Figure 3.27: int edge simplification, initial condition.
When int edge simplification encounters the first edge, the procedure considers node 2
as node_to_discard. For this reason, both the edges (2,3) and (2,4) are copied with
1 as the new end, creating (1,3) and (1,4), as in Figure 3.28
Figure 3.28: int edge simplification, wrong procedure - step 1.
The next iteration of the for loop operates on the edge (2,3). The algorithm con-
siders 3 as node_to_discard, and so, it creates (1,2) as the copy of (1,3), and (2,5)
as the copy of (3,5).
At the end of the loop, when nodes 2 and 3 are removed, the only edge that is
preserved is (1,4). In Figure 3.30, it is evident how this procedure leads to a loss of
edges. In addition, one node is now isolated.
At this point, in the following example, the correct procedure is applied. In step 2,
the procedure “sees” that after the edge (1,2) there is another edge that belongs to
the same category. For this reason, (1,2) is skipped and the edge copy algorithm is
applied on edges that are connected in 3. This creates (2,5) as the copy of (3,5), as
depicted in Figure 3.31.
Figure 3.31: int edge simplification, correct procedure - step 1.
At this point, the loop is complete: node 3, together with the edge (3,5), is removed.
Then, a new cycle starts.
The procedure now considers the edge (1,2), and so 2 as node_to_discard. Edge
(2,5) is copied to (1,5), and (1,4) is created as the substitute of (2,4).
Finally, node 2 is deleted, and NetworkX takes care of removing (2,5) and (2,4). The
result is now correct: the two black edges have been preserved and the remaining
nodes are still connected to each other.
The other idea that must be clarified is the necessity to run int edge simplification
two times, exchanging the indexes. First of all, considering one edge, performing
this action is equivalent to take node_to_discard as node_to_keep and vice-versa.
The reason behind this choice is explained in the next example, considering the part
of the network in Figure 3.35, and supposing that the edges are defined as (1,2) and
(3,2).
Figure 3.35: Problematic configuration for int edge simplification function. The two
edges are defined as (1,2) and (3,2).
In the first execution of the function, the edge (1,2) has (2) as node_to_discard.
Since (3,2) is in the same category as (1,2), the first edge is skipped. Now the
function analyzes (3,2). In this case, node_to_discard is again (2), and this edge
is skipped too. The result is that this pair of edges is never removed.
Swapping the order of the elements allows solving this problem. In fact, in the
second execution, the edge (1,2) has (1) as node_to_discard. The edge can be
removed because it does not have on node (1) an adjacent edge belonging to the
same category. The same can be done on (3,2). This allows removing both the
target edges.
merge adj nodes function The merge adj nodes function is used when a node
must be removed from the graph, but the contained information has to be preserved.
It is executed using a command similar to the following one:
fn . merge_adj_nodes ( network_graph , node_to_keep , node_to_discard )
The performed actions are:

The name variable of node_to_discard is appended to the one of node_to_keep,

interposing @ as separator. The same procedure is applied to the attribute
misura.
The values of P, T, and height are copied from node_to_discard only if they
are not specified in node_to_keep.
The category of the node, together with supply, is important because it

defines the sign convention of the market values. These variables are managed
with a priority scheme:
– If at least one of the two nodes is NS, the merged item is NS too, and the
supply is 1. The market values are simply summed up if both nodes are
NS, instead, the NO node changes the sign of its market.
– If both nodes belong to the NO category, their market values are summed,
the merged category remains NO and the supply is nan.
– If one of the two nodes is neither NS nor NO, it is neglected and the other
item’s attributes are copied.
– If both nodes are neither NS nor NO, the merged node’s variables are
category = ’None’, market = 0, and supply = nan.
The list tpMarket of the result node contains both the tuples of the original
nodes. This method allows preserving the information about which consump-
tions are aggregated into a single node during the simplification process.
safe add edge function The safe add edge function secures the creation of a
new edge inside a graph. In fact, the add_edge function included in NetworkX does
not check if an edge is present in the graph when it is required to instantiate a
new one. If this happens, the old one would be completely overwritten by the new
one. The custom function safe add edge performs this check, and follows a priority
criterion to decide how to solve the conflict.
The algorithm takes as inputs the graph, node1 and node2 as ends of the new edge,
and a series of attributes that will be saved in the edge. It is designed as a series of
nested ifs:
The first step is to check if there is an edge connecting node1 and node2. If
this is not the case, the new edge is added copying all the attributes and the
procedure ends.
Figure 3.36: In step A the edge (28,30) is removed. So the edge (27,30) should
be moved to (27,28), as in step B, but there is an edge yet. Since they belong to
the same “high priority” category, they must be both preserved. So, in step C, the
dummy node (3) is created, and it is connected to (28) using the joint (3,28).
If the edge (node1, node2) exists, the function checks its category. If it belongs
either to pipe or compressor station categories, there are two sub-cases:
– If the new edge is neither a pipe nor a compressor station, it has a lower
priority and its creation is skipped. The function stops.
– If also the new edge is either a pipe or a compressor station, the inform-
ation contained in both the elements is important. Since the program is
working with a graph that does not allow multiple edges between a pair
of nodes, a dummy node is created to keep both elements. The new edge
is placed between node1 and the new node, and then it is connected with
node2 placing a joint. The joint will be then removed after the conver-
sion of the graph in a multigraph. Figure 3.36 shows an example of this
procedure.
If the existing edge is neither a pipe nor a compressor station, there are again
two sub-cases:
– If the new edge is a high priority element, the new edge overwrites the
existent one.
– If both the edges are neither pipes nor compressor stations, the creation
of the new element is skipped, preserving the actual one.
Valve removal procedure For the removal procedure of the valves, an additional
caution has been used. In fact, the elimination of the closed valves is the first
operation that the program performs during the simplification phase. This is done
because, in this way, the successive algorithms can follow every path they encounter,
without the risk that they are crossing a closed element. The removal process is
simple: a for loop checks every valve’s MODE, and deletes the edges whose value
Figure 3.37: Example of valve configuration that can be solved only using slowEle-
mentRemoval.
represents a closed state.
In addition, in the removal process for valves, another function is employed. In

fact, edge removal is not able to process some complex configurations of valves that
can be found in the network. In particular, this happens when these elements
form closed loops, like the example in Figure 3.37. To solve this problem, the
function slowElementRemoval has been developed. The algorithm considers only
the first edge found for every loop iteration; this allows neglecting the checks on the
neighbouring edges, explained in int edges simplification. After the first removal,
the loop is started over and the index of elements to remove is refreshed. This
approach is slower than the previous one, but it manages every configuration in
the network. The performance issue is mitigated by executing this function after
the edge_removal one: in this way, the set of elements that must be analysed is
reduced.
simplified network function At this point, only necessary elements are left in
the graph, but pipes are still split into hundreds of small segments. This condition
is not optimal for the simulation because the discretization in finite volumes should
be decided considering the procedure described in Section 2.6.2. The function sim-
plified network has been built to unify, where possible, the small parts. During this
process, the information about the position of the consumptions distributed along
the network must be preserved, in order to not modify the pressure profile in pipes.
This is provided by a list of consumptions and relative locations that is stored in
every edge. In addition, the procedure performs also some final modifications.
At this point, it is necessary to introduce a multigraph. In fact, in the original

map, the detailed description of every element and the presence of joints prevent
the creation of parallel edges between a pair of nodes. Instead, the goal of the
simplification process is to have the simplest representation of the components, as
Figure 3.38: Comparison between the original, detailed map and the simplified one.
shown in Figure 3.38.
The first operation that the function performs is to copy every edge that is not a pipe
into the new multigraph. It generates the list otherEdges, containing compressor
stations, control valves and a small number of joints created by the safe add edge
function. Then, a for loop browses it, instantiating the corresponding edges and
nodes, preserving all the information. These nodes are also saved in a list called
otherNodes. Then, it stores in the set startNodes every node of the original map
that has a number of connections different from two. After that, it merges them into
a single set and the contained nodes are used as the start point for the simplification
process. Finally, it creates the list seenNodes: it is used to trace the nodes that the
function has already elaborated.
At this point, a while loop iterates on the startNodes set, until it has remaining
elements (Figure 3.39). The primary operation is the execution of the function
Figure 3.39: Flowchart of the external loop around the function followPipe.
Figure 3.40: Flowchart of the internal structure of the function followPipe.
followPipes. It is a recursive algorithm that takes as inputs the graph G and the
multigraph M, the initial node node and the node of the previous call prvNode, the
startNodes set, the variable seenNodes, and the list path. The last item is used
to store the nodes encountered during the union procedure.
Once the path list is instantiated, the initial call is done, using as argument node
the first element of startNodes, and 0 for prvNode.
At first, the algorithm, whose flowchart is reported in Figure 3.40, checks if the
considered node is included in startNodes. If this is the case, there are two possib-
ilities: either prvNode is different from zero, and this means that the function has
reached a relevant point and it must be stopped, or prvNode is equal to zero. This
means that the function is at its first iteration. In this case, the list, called NeiList,
of nodes neighbours of node is generated. A while loop iterates the first element of
the list, performing three checks. Initially, it verifies if the first element neiList[0]
is in seenNodes: if this is true, the node has been considered before, and it is not
necessary to follow that path again. So, that node is removed from neiList and the
loop starts over. On the contrary, the edge between node and neiList[0] is con-
sidered: if it is a pipe, it must be followed, and so the function followPipe is called
again using neiList[0] as argument for node and node as argument for prvNode.
When the recursive call of the function is completed, node is removed from the list
and the loop can be restarted. Instead, if the edge is not a pipe, the procedure must
not follow that path. The node is removed from neiList and the loop is started
again. Finally, when neiList is empty, the loop can be stopped and the procedure
ends.
Going back to the first turning point in the algorithm, if the node considered is
not in startNodes, the procedure generates the list of neighbours. In this case,
this list contains necessarily two elements: prvNode and neiList[0]. The edge
(node, neiList[0]) undergoes two checks that can lead to three different out-
comes:
If the edge is not a pipe, the procedure is stopped without any other action.
If the edge is a pipe, but its attributes are different from the pipe between
(prvNode, node), the procedure is stopped and node is added to startNodes.
In fact, if there is a variation in pipes attributes in correspondence of a node,
that node must be considered as relevant and a starting point for comparePipes
function.
If both edges are pipes and their attributes are equal, the recursive proced-
ure must continue, calling followPipe on neiList[0], considering node as
prvNode.
The comparison between two pipes is performed using a function, called checkToler-
ances, that checks if their diameters and roughnesses are equal, with a 2% tolerance.
This approximation is introduced to avoid false negatives due to finite-precision rep-
resentations of the values that would have lead to create two separate pipe models
with a negligible difference in their geometrical features.
When the recursive search reaches an end point, the information about that node is
returned back to the main call. At this step, knowing the starting point, the ending
node, and the path, the function creates a new edge in the multigraph representing
the unified pipe. Its attributes are generated by the algorithm pipeMerginAttributes,
that, following the path, reconstructs:
The total length of the pipe;
The name string, created appending the segments name one to each other;
The array containing the consumption values that were placed on the nodes
merged into the unified pipe;
The array that places those consumptions in the correct point of the pipe,
using, as a reference, a linear coordinate that has its origin in node. This
information is stored in the attribute array_origin.
The array that contains the tpMarket IDs related to the above-mentioned
consumptions.
The simplified network function is then completed carrying out three simple ac-
tions. At first, the joints, created by safe add edge to avoid the overwriting of edges
between the same pair of nodes, are useless in a multigraph. So, they are simply
removed and the diverted edge is brought back to its original ends. Then, all the
nodes without connections, residual items from some of the simplification functions,
are deleted. Finally, there are compressor stations that are represented by two edges,
each containing only a part of the machines. Since the connections of each station
to the network coincides, it is more efficient to aggregate the two edges into a single
element.
At this point, the simplification of the network is complete. The multigraph can be
drawn on a map, as explained in Section 3.5, or can be exported as Modelica model,
as described in Section 3.6.
3.5 Map export

In every phase of the work, the representation of the network configuration, together
with the data describing the status of every component, is essential. For this reason,
the possibility to access this information in a expressive and convenient way is an
added value.
The best way to represent the considered network is to draw a map. In fact, showing
the approximate position of each component, and their mutual connections, is the
most natural method to understand the flow of the gas in the pipes.
In the following sections, the drawing process is explained.
3.5.1 Graphical representation of nodes

In the considered network, nodes are always represented with circles, and they are
organized into three categories:
NS nodes represent points where the gas, generally, enters the network. They
are depicted as light green circles.
3.5. MAP EXPORT 89
(a) The representation of the original net- (b) The representation of the simplified net-
work, as described in the source files work, as obtained after the manipulations
explained in Section 3.4.2
Figure 3.41
Figure 3.42: A real example of the representation of nodes in the map.

NO nodes are locations where the gas is extracted from the network, and they
are represented as orange circles.
The remaining nodes are not identified by an acronym. They represent points
where different edges interconnect to each other, and they are never related
to an entering or exiting mass flow-rate. They are rendered as slightly smaller
grey circles.
Figure 3.42 contains a representation of nodes from the three categories in the map.
In addition, using some options of the function generateNetworkMap, it is possible

to highlight some specific nodes. These elements are depicted using bigger, blue
circles.
These graphical objects contain textual and numerical information about the nodes.
In fact, when the mouse pointer hovers on the circle, a small label with the node
ID appears. Clicking on the shape, a pop-up that encloses every value of the node
attributes is shown. Pop-ups are explained in details in Section 3.5.3.
The code that manages the creation of the shapes of the nodes is reported in Figure
3.43.
3.5.2 Graphical representation of edges

In the network map, the edges are represented as segments of different colours, as
listed in Table 3.1. A graphical example is displayed in Figure 3.44
Edge category Colour

Compressor station Purple
Pipe Black
Valve Blue
Control valve Pink
Non-return valve Red
Joint Yellow
Metering station Dark green
Table 3.1: Colours of the segments representing the different edge categories.
When the map depicts a multigraph, also the line thickness can vary: in fact, it is
proportional to the number of edges in parallel between a pair of nodes.
When hovered, edges show a label containing their ID; when clicked, a pop-up with
complete information appears. This is explained in Section 3.5.3.
3.5. MAP EXPORT 91
1 for elem in network_graph . nodes :

2 nodePopup = genera teNod ePopup ( elem , network_graph . nodes [ elem ])
3 # If the node is in the highlight list -> blue and bigger
4 if elem in list :
5 radius = 13
6 weight = 3
7 color = ' blue '
8 else : # If not in list
9 if network_graph . nodes [ elem ][ ' category ' ] == ' NO ' :
10 radius = 9
11 weight = 3
12 color = ' orange '
13 elif network_graph . nodes [ elem ][ ' category ' ] == ' NS ' :
14 radius = 9
15 weight = 3
16 color = ' #00 dd00 '
17 else :
18 radius = 7
19 weight = 3
20 color = ' grey '
21 folium . CircleMarker (
22 c o n v e r t C o o r d i n a t e s F o r F o l i u m ( ' node ' , elem , network_graph ) ,
23 radius = radius , weight =2 , color = color , fill = True ,
24 popup = nodePopup , tooltip = f ' Node { elem } ' ) . add_to ( layer )
Figure 3.43: Code that manages the creation of nodes on the map.
In addition, a big, yellow circle is drawn around the compressor stations that have
at least one gas pumping unit active. This allows to immediately identify which
machines are effectively moving the fluid in the network. The code that draws these
shapes is reported in Figure 3.45.
Instead, a red circle appears in multigraphs when, between two nodes, there are
edges in parallel that belong to different categories.
3.5.3 Pop-ups
Every element on the map can show, when clicked, a pop-up that contains the values
of its attributes.
Since the amount of data that an item can enclose is big, the organization of the
pop-up content is substantial. In fact, badly organized values can be distracting
or misleading. For this reason, every pop-up is divided into two sections: the first
one, which contains the most relevant data, automatically appears when the item is
clicked, the second one is shown only when requested.
Figure 3.44: Part of the map containing every category of edges of the network.
The purple edge represents a compressor station, the black ones are pipes, the blue,
red, and pink elements depict valves, non-return valves, and control valves respect-
ively, yellow segments show joints, and dark green ones are metering stations.
Dashed lines stand for closed elements.
Mp = 0
for k in stationCoeff [ name ]. keys () :
Mp += fn . co u n tA c t iv e M a ch i n es ( stationCoeff [ name ][ k ][ ' Nomi macc . ' ] ,
stationNSI [ name ])
if Mp > 0:
folium . Circle ( c o n v e r t C o o r d i n a t e s F o r F o l i u m ( ' center ' , elem ,
network_graph ) , radius =10000 , color = ' yellow ' , fill_color = '
yellow ' , fill_opacity =0.2 , tooltip = f ' Active machines : { Mp } ' )
. add_to ( backLayer )
Figure 3.45: Code that generates the circles highlighting active machines on the
map.
3.5. MAP EXPORT 93
Nodes pop-up The most interesting attributes saved in nodes are height, category,
market and pressure. They are important because they are good indicators when
estimating the flow direction. For this reason, they are shown in the first part of
the pop-up.
The other variables that belong to every node are hidden in the second slice of the
page: x and y, temperature, supply, misura, name, tpMarket. In particular, the
last three items are put in this section not to overwhelm the user with very long
strings. In fact, during the merge adj nodes procedure, these variables are appended
to each other, growing in length.
Figure 3.46 shows examples of collapsed and expanded pop-ups referred to nodes.
Edges pop-up The shapes and the contents for edges pop-ups (Fig. 3.47) can be
different. So, the first step is to decide which parameters are important for every
category of edges, and for this reason, which ones must be always on the top of
the pop-up. Three of them are considered important for every category: category,
edgeid, and array_origin.
The other relevant attributes are:
for pipes: length, diameter, roughness, pguess;
for compressor stations: name, SPI, SPO, MM, SM;
for control valves: name, SPO, MM, SM, diameter;
for valves: name, MODE, diameter;
for non-return valves: name, diameter;
for metering stations: name, MM, diameter;
for joints: no additional attributes.
All these parameters are directly read from the graph. It contains only the data
that are also available in SIMONE, the software used by Snam, and so reported in
the associazioni file (Sec. 3.1).
Compressor stations have tens of attributes that are not stored in the graph, since
they are not listed in the associazioni file. The most important are the states
of the machines, but there is also information about the rotational speed, fuel con-
sumption, and ambient temperature. These data, extracted from the NSI file using
Figure 3.46: Example of a node pop-up collapsed and expanded.
Figure 3.47: Example of a pop-up referred to an edge. In this case, it describes the
attributes of a control valve.
3.5. MAP EXPORT 95
the function extractCompressorsInfo, which is explained in Section 3.6.2.5, are listed

on the second page of the pop-up (Fig. 3.48).
In addition, when drawing a multigraph, a single segment can represent multiple

edges, and so the pop-up must do the same. In fact, even if multiple edges can be
placed between the same pair of nodes, they can have different attributes. For this
reason, each edge has its dedicated page inside the pop-up, containing its values.
During the network simplification process, sections of pipes in series are merged
into a single, longer, pipe. In this process, many nodes disappear from the graph.
The consumptions that were associated with them are now stored in the variable
consArray inside the composite pipe. Their positions, expressed as linear coordinate
starting from array_origin, are memorized in posArray, while the list of their
tpMarket IDs is kept in tpMarketArray. To make these values more accessible,
they are organized in a table (Fig. 3.49).
The code that creates the maps (Fig. 3.50) is contained in the function called
generateNetworkMap. It takes as input the variables containing the graph, the
list of nodes to highlight, the name of the map as reported in the legend, and
the additional information about the compressor stations. This function returns a
Folium.FeatureGroup, a variable managed by the library Folium that contains a
layer of the map. This entity can not be directly opened in the web browser, but
it must be added to a Folium.Map object, and then exported. The advantage of
this approach is that the Map entity can contain several layers that can be later
exported into a single HTML file. In this way, the user can choose the layer to
analyse in the web browser without neither generate a new file nor reloading the
page. This possibility is exploited in the main program: both the original map
and the simplified one are added as a layer to the Map object. In Figure 3.41, it is
possible to see the two different layers and, in the top right corner, the box that
allows switching between the two.
A quicker alternative is offered by the function exportGraphToMap: using this com-

mand, an HTML file is produced, containing a single layer representing the selected
graph.
Figure 3.48: Example of a compressor station pop-up: the basic information page
on the left, the advanced data on the right.
Figure 3.49: Example of a pop-up containing the information about a composite

pipe.
3.6. MODELICA MODEL GENERATION 97
1 map . add_child ( mp . ge ne ra te Net wo rk Ma p (G , name = ' Original

Network ' , valveModeClose = valveModeClose ) )
2 map . add_child ( mp . ge ne ra te Net wo rk Ma p (M , name = ' Simplified
network ' , stationCoeff = stationCoeff , stationNSI =
stationNSI ) )
3 map . save ( ' combinedMap . html ' )
4 mp . exportGraphToMap (M , ' map . html ' )
Figure 3.50: Examples of the code that generates the map. In line 1, the layer
containing the original map is rendered. The third argument specifies which valves
are closed. In this way, they will be drawn with a dashed line. In line 3, the third
and fourth arguments are used to specify additional information about the stations.
In line 4, the combination of the two layers is exported to the file combinedMap.html.
In line 5, the function that directly creates the HTML file from the graph variable
M is used.
3.6 Modelica model generation

The graph of the network obtained at the end of Section 3.4.2 has been simplified
with the purpose of facilitating the simulation. At this point, the program can gener-
ate the file containing the Modelica model. To do so, Python writes a plain text file,
saved with extension .mo, that contains the instructions that, later, OpenModelica
will execute.
This file is a key in the entire process: all the data, collected and elaborated in
the previous steps, are used in it. For this reason, the function writeModelicaModel
performs a lot of operations and employs several sub-functions. Its working principle
is explained in Section 3.6.1.
3.6.1 writeModelicaModel function

The function writeModelicaModel generates the .mo file that contains the Modelica
model for the network. It performs the following operations:
Extracts the information about the topology and the components from the
multigraph that represents the simplified network;
Imports the data to assign to the items from mercato and NSI files;
Instantiates the model components, defining their parameters and their start
values, used during the initialization of the simulation;
Declares the sources for the signals describing consumptions, set-points for
machines and control valves, and boundary conditions;
Connects every component, following the graph;
Sets the simulation parameters, like the start and stop time and the number
of intervals for the temporal discretization.
To do so, it takes several parameters as input:
The network graph;
The destination path for the .mo file;
Three variables that contain the compressor station information;
The dataframe containing the control valves statuses and set-points;
The dataframe containing the data about the mass flow-rates entering or ex-
iting the network;
The path to the folder that contains the series of NSI files;
The list of the boundary conditions, expressed as pressure or flow-rate sources

connected to specified nodes.;
The number of finite volumes that the pipes must be split into;
The start and stop times for the simulation;
The steady-state flag, that activates the steady-state initialization of the sys-
tem;
The constantInput flag: if it is set to True, every signal in the simulation is

kept at its initial value;
The relaxationTransient flag, that activates an initial transient time before

the actual simulation interval.
The constantInput option is implemented creating different lines of code in the

model. In fact, when a constant input is requested, the signal source is set to be a
Modelica Standard Library Constant or IntegerConstant block (Sec. 3.6.2.1), and
its value is directly read from the function variables. Instead, when a time-varying
input is required, the signal is generated by a Modelica Standard Library TimeTable
or IntegerTable block. Its time series is extracted from the NSI files by the func-
tion importTimeSeriesNSI (Sec. 3.6.2.2), and then a string containing the values in
a format compatible with Modelica is prepared by convertTimeSeriesForModelica.

This function manages also the option relaxationTransient, and its behaviour is
explained in Section 3.6.2.4.
This function can be also used to generate a model further simplified, in fact it is
possible to integrate a function that merges a group of parallel pipes, placed between
the same pair of nodes, into a single equivalent pipe. This function is useful for the
generation of the model for the MILP optimizer.
3.6.2 generateModelicaStationModels function

The function writeModelicaModel instantiates the compressor stations using only a
bunch of lines of code. This is possible because, in the first part of that procedure,
the ancillary function generateModelicaStationModels is called. This algorithm cre-
ates in the .mo files a series of models that describes the internal structure of the
station, using the components described in Chapter 2.
Every station has the same external interfaces: in fact, they always have inlet and
outlet ports for the gas flow, three inputs for the flow-rate and pressures set-points,
and a variable number of integer inputs for the signal that activates the different
gas pumping units.
Inside, as described in Section 2.11, the following components are instantiated:
Inlet and an outlet port for the gas;
A variable number of gas pumping units, connected with the ports;
Input connectors for the set-points;
Equations for the ideal control, conforming with the number of units present
in the considered station;
Integer input ports, that are used by the ideal control to determine how many
machines are required at every time instant;
A bypass valve, that allows the gas to flow unimpeded through the station.
The considered function can procedurally generate these models, that extend the
base class BaseCompressorStation, obtaining the parameters and writing them
in a format compatible with Modelica (that preferably uses SI units) using the
procedures describes in the next sections.
table = [0 , 1; 5 , 2; 7 , 1]
Figure 3.51: Comparison between the output obtained using a TimeTable (red line)
and an IntegerTable (blue line) with the same parameter table.
3.6.2.1 Modelica signal blocks
When creating the network Modelica model, four types of signal generators are used.
These items can be constant or time-varying, and with real or integer output.
The two constant sources have one parameter, called k, that represents the value
that the block output y assumes. Their names are Constant and IntegerConstant.
The time-varying integer signal source takes as parameter a table, called table, and
it is set to keep the output y at the last assigned value, until a new one comes up.
Its name is IntegerTable.
The time-varying real signal source has table as parameter, but the values are here
linearly interpolated when output to y. Its name is TimeTable.
A graphical comparison between the two variable signals is shown in Figure 3.51.
3.6.2.2 importTimeSeriesNSI function
The function importTimeSeriesNSI is used to extract the time trend of a measure-

ment. It analyses a dataframe that contains multiple NSI files (Sec. 3.1), filter-
ing only the rows with the desired NSI ID. Since these files are generated every 4
minutes, the signal is sampled with that period. Then, the function returns a list
that contains, in turn, another list formed by:
Temporal information, that can be expressed as timestamp or time difference,

in seconds, from the initial time;
The value contained in the row.
3.6.2.3 importTimeSeriesMarket function
The importTimeSeriesMarket function is used to extrapolate from the Mercato data-

frame the temporal evolution of mass flow-rates that are extracted or injected in
correspondence of a node.
Since during the network simplification process several nodes can be aggregated into
a single element, the input of this function is a list of the tpMarket tuples.
In the example of Figure 3.52, two different sets of measurements must be con-
sidered together, because the node in location_1 has been merged with the one in
location_2. Considering that the sign convention for NS nodes is opposite to the
one for NO elements, the values obtained for location_1 are multiplied by -1 before
the sum.
In addition, for the most part, measurements are taken with a time interval of 4
minutes, the same as NSI series, but some of them are sampled every hour. To have
a coherent representation that simplifies the sum of different time series, the data
with a longer acquisition period are repeated 15 times each, mimicking a sampling
interval of 4 minutes. Finally, the information is returned as a list of lists, each
containing
Temporal information, that can be expressed as timestamp or time difference,

in seconds, from the initial time;
The value of flow-rate evaluated in the procedure.
tpMarket = [( ' location_1 ' , ' NS ' ) , ( ' location_2 ' , ' NO ' ) ]
Figure 3.52: Example of tpMarket variable taken as input by importTimeSeries-

Market.
3.6.2.4 Conversion functions
During the creation of the Modelica model, a big quantity of data has to be trans-
ferred from the Python variables to the Modelica components inputs and parameters.
To facilitate this procedure, two functions have been developed: convertArrayFor-
Modelica and convertTimeSeriesForModelica.
convertArrayForModelica function This function takes as input a Python

list, containing some numerical values, and returns a string filled the same data in
an array compatible with Modelica (Fig. 3.53). The principal modification that
is performed is changing the parentheses that contain the values: in fact, Python’s
lists are delimited by square brackets, while Modelica arrays by braces. The function
also converts pressures from bar to Pa, since Modelica preferably uses SI units.
[4 , 1.12 , 5 , 13.3 , 7] { 4 , 1 .12 , 5 , 13 .3 , 7 }
(a) Python list (b) Modelica array
Figure 3.53: Example of the conversion performed by the function convertArrayFor-

Modelica
convertTimeSeriesForModelica function This procedure takes is used to

convert a time series, usually generated by importTimeSeriesNSI or importTimeSer-
iesMarket, into a table compatible with the Modelica signal generators. In addition
to the simple change of format, the function is able to perform a multiplication
between the data and the parameter conversionParam, and it can prepare the time
series for the relaxationTransient. This transient consist of a time interval that
precedes the actual simulation time. In this period, the system state, initialized with
the values obtained from the graph, relaxes and reaches a more stable condition.
This step is necessary because the available information about the initial condition
is not good enough to initialize the Modelica model in a fully consistent way. To
achieve this, every signal must be kept constant. The TimeTable block interpolates
the values between the specified points, and this can lead to an unwanted drift when
the simulation exits from the table boundary, as in this case. To avoid this problem,
the function adds a copy of the first value delayed by a small-time interval. Doing so,
the block sees two contiguous entries in the table and, when they are interpolated,
they lead to a constant value. An example of this behaviour is depicted in Figure
3.54.
3.6.2.5 extractCompressorsInfo function
In the procedure that merges the different files into a single dataframe, only a small
group of attributes about the compressor stations is preserved. In fact, the REQ file
contains the associations only for the inlet and outlet pressures and their set-points,
the value and the set-point for the mass flow-rate.
Figure 3.54: Example of the behaviour of the convertTimeSeriesForModelica func-

tion with active relaxationTransient option.
These parameters are not enough to properly control the compressor stations. For
this reason, the function extractCompressorsInfo analyses the NSI dataframe, match-
ing the rows that refer to a station comparing the first part of the name with the
same string from a known datum. Then, it organizes a dictionary for every station
that uses as key the name of the attribute, and as value a tuple that contains:
The actual value of the attribute;
Its unit of measure;
Its misura ID.
Finally, every dictionary is saved into a container dictionary, that uses as keys the
names of the compressor stations.
3.6.2.6 importCompressionStationCoefficients function
This procedure is used when instantiating the model of a compressor station. It

acquires the coefficients for the equations that model the compressors and the gas
turbines. It reads two different Excel files, importing the attributes for every machine
type present in the network. Then reads a third file that associate every station with
the correct models of units. Then saves the matched data into a dictionary organized
by station and gas pumping unit.
3.6.2.7 countActiveMachines function
The countActiveMachines function can scan a dictionary generated by extractCom-

pressorsInfo, determining how many active gas pumping units it contains. It is used
both in the map generation process, to highlight the active stations, and while cre-
ating the signal source that determines how many machines are currently running
in the station.
The model resulting at the end of the described procedure is then validated in two
case-studies reported in Chapter 4.
Chapter 4
Model validation
The validation of the model obtained from the complete toolchain has been per-
formed on the longest, complete, and consistent dataset provided by Snam at the
time of this writing.
The dataset describes the network behaviour for six hours, starting from the 1st
of March 2021 at 20:00 to the 2nd of March at 01:56. In the reported period, five
stations are active, and three of them are located in the central-southern part of
Italy. This configuration is engaging because it allows to examine the pressure
losses in the pipes and how they are recovered by the compressing units. This work
is performed in the first case-study in Section 4.1. These data contain also the
start-up transient of the Malborghetto station, an interesting phenomenon that is
analysed in the case-study in Section 4.2.
105
106 CHAPTER 4. MODEL VALIDATION
4.1 Case-study 1: Central-Southern Italy

The first case-study is focused on the part of the network placed in the central-
southern part of Italy. This area, depicted in Figure 4.1, contains:
Import points of Mazara del Vallo (TP), Gela (CL), and Melendugno (LE);
Storage site of Fiume Treste - San Salvo (CH);
Compressor stations of Enna, Messina, Tarsia (CS), Montesano (SA), Be-

nevento (BN), and Gallese (VT);
More than 3660 km of pipes.
Figure 4.1: The map of the considered portion of network. The names of import
points, storage site, and compressor stations are reported.
4.1.1 Description of the scenario

In the considered time interval, Enna, Tarsia, and Gallese are the active stations,
while the remaining ones are bypassed. The first station required a manual interven-
tion to disable the ideal control described in Section 2.11.1, in fact the data shows
that the only active pumping unit is saturated at its maximum speed. The other two
4.1. CASE-STUDY 1: CENTRAL-SOUTHERN ITALY 107
Figure 4.2: Evolution of the mass flow-rates at the injection points in the network.
stations are correctly working controlling the inlet pressure. Their working mode
has been automatically selected by the simulator, without any external contribution.
In addition, also the opening of the bypass valve in the two inactive stations has
been automatically managed by the simulator.
All the import nodes are active, and the storage site is injecting its flow into the
network. The values of these mass flow-rates are reported in Figure 4.2.
Since most part of the natural gas of the network is imported in Mazara del Vallo
and Gela terminals, and since these two flows mixes at the inlet of the Enna station,
a good hypothesis for the gas composition estimation is to use the measurement
performed by the gas chromatograph installed in that station.
The boundary conditions are the following:
Mazara del Vallo: pressure source, followed by a control valve that regulates
the flow-rate and that monitors the maximum outlet pressure value;
Gela, Melendugno, and Fiume Treste - San Salvo: flow-rate source with time-
varying prescribed values read by the SCADA system;
Gallese: time-varying pressure source at the outlet of the station, with pre-
scribed values acquired from the SCADA log.
Enna station Simulation SCADA Difference Relative error

pout [bar] 69.36 69.15 0.21 0.30%
pin [bar] 56.95 57.26 -0.31 -0.54%
w [kg/s] 634.60 637.35 -2.75 -0.43%
Table 4.1: Numerical comparison between the simulation and the SCADA measures
for the station of Enna.
Figure 4.3: Comparison between the SCADA measurements and the simulated val-
ues for the inlet pressure of the Enna station.
4.1.2 Results
The first analysed quantity is the pressure at the inlet of the Enna station. As
shown in Figure 4.3, the initial values are different by 0.4 bar. This mismatch is
caused by lack of precision of the information used for the initialization. In fact, the
trend is well reproduced by the simulation, and the final ∆p between the SCADA
and the simulated values is reduced to 0.3 bar. This happens because, during the
simulation, the initialization error is balanced out.
The numerical data about the Enna station are reported in Table 4.1. There, the
values at the end of the simulation (01:56) are considered.
The second considered group of quantities contains the pressures of the outlet of
Mazara terminal, of Gela terminal, and the Enna inlet is used as reference. Since
both terminals pressures are not controlled, their values is determined by the pres-
sure losses on the pipes that connect them to the station. The real trend of the
quantities is well represented, but the simulation of the Gela pressure is underes-
timated by 0.7 bar. This difference is caused by the wrong temperature considered
for the gas. In fact, for the modelling of the gas in the pipes, the temperature is
assumed constant at the values of 15 ◦ C. The actual value is near to 8 ◦ C.
This mismatch is not propagated inside the compressor station because, as explained
in Section 2.7, the gas at the inlet is considered with its real temperature, read by
the sensor.
The numerical data are reported in Table 4.2 and they are plotted in Figure 4.4.
Pressure [bar] Simulation SCADA Difference Relative error

Mazara del Vallo 66.16 66.29 -0.13 -0.20%
Gela 60.37 61.13 -0.76 -1.24%
Enna (inlet) 56.95 57.26 -0.31 -0.54%
for pressures at Mazara del Vallo terminal, Gela terminal, and Enna station inlet.
ues for the pressures at Mazara terminal outlet, Gela terminal, and Enna inlet.
Pressure [bar] Simulation SCADA Difference Relative error

Tarsia (outlet) 63.94 64.58 -0.64 -0.99%
Melizzano (inlet) 57.31 57.42 -0.11 -0.19%
Gallese (inlet) 49.04 48.99 0.05 0.10%
for pressures at Tarsia station outlet, Melizzano station inlet, and Gallese station
inlet.
ues for the pressures at Tarsia station outlet, Melizzano station inlet, and Gallese
station inlet.
The third group of considered measurements is represented in Figure 4.5. It contains

the pressures at Tarsia outlet, Melizzano inlet and Gallese inlet. The figure shows
that the trend is followed by the three measurements, but the one in Tarsia is slightly
underestimated. The relative error, reported in Table 4.3, is 1%.
Then, the mass flow-rate in the bypass of the station of Melizzano is considered
(Tab. 4.4). This is the only available measurement of a non controlled flow-rate
that is neither exiting nor entering the network. In fact, in active stations, the flow-
rate is influenced by the control strategy, and, in stations like Messina, the bypass
flow is not monitored. Figure 4.6 shows how, after a first phase where the measured
data are affected by uncertainty, the simulation and the SCADA are very similar.
The last group of measurements that has been analysed concerns the active stations.
Mass flow-rate [kg/s] Simulation SCADA Difference Relative error

Melizzano 312.80 310.83 1.97 0.63%
for the mass flow-rate in Melizzano station.
ues for the mass flow-rate in Melizzano station.
In Figure 4.7, the trend of the dimensionless rotational speed of the stations of Enna,
Tarsia, and Gallese are depicted. In Table 4.7, their final values are reported. The
Enna simulation value is fixed, as described in Section 4.1.1. In Figure 4.8, the fuel
mass flow-rates are represented. Finally, in Table 4.6, the numerical values at the
last time interval of the simulation are reported.
The time requested for the execution of the entire toolchain on a laptop (Intel i7
4700HQ 2.40 GHz, RAM 16 GB, SSD for programs and data storage) is near to 5
minutes. In Table 4.5, the times and the simulation parameters are reported.
Phase Time [s] Parameter Value

Python script 165 Start time -3600 s
Model translation 92 Stop time 21360 s
Model compilation 48 Interval 30 s
Model simulation 23 Number of equations 16432
Number of variables 16432
Total 328
Table 4.5: Toolchain execution times and simulation parameters in the Case-Study 1.
ues for the dimensionless rotational speed of the active units in Enna, Tarsia, and
Gallese stations.
Dimensionless
Simulation SCADA Difference Relative error
Rot. Speed [-]
Enna 1.050 1.051 -0.001 -0.1%
Tarsia 0.840 0.863 -0.023 -2.7%
Gallese 0.864 0.824 0.040 4.8%
for the dimensionless rotational speed of the active units in Enna, Tarsia, and Gallese
stations.
ues for the fuel mass flow-rate of the active units in Enna, Tarsia, and Gallese
stations.
Fuel Mass
Simulation SCADA Difference Relative error
Flow-Rate [kg/s]
Enna 1.282 1.198 0.084 7.1%
Tarsia 1.082 1.053 0.029 2.8%
Gallese 1.157 1.002 0.155 15.5%
for the fuel mass flow-rate of the active units in Enna, Tarsia, and Gallese stations.
4.2 Case-study 2: Malborghetto station start tran-

sient
The second case study considers the region of the network that starts from the
Austrian border, contains the station of Malborghetto and Istrana, and arrives near
the city of Vicenza. It contains also two import points: Tarvisio and Gorizia. The
total length of the pipes is near to 760 km. The map of this part of network is
depicted in Figure 4.9.
4.2.1 Description of the scenario

This part of network has been considered because, during the time interval recorded
in the SCADA files, one gas pumping unit in the station of Malborghetto is switched
on. This event is significant because it gives the chance to verify if the model of the
network reacts with the same time constants of the real system when it is subjected
to a sudden variation.
To highlight the effect of the start-up of the machine, avoiding the influence of
any other source, the pressure sources follow the profile of the pressure measured
in the SCADA file. Then, the simulation is interrupted when the pressure wave
reaches the final points (Node 5394 and Node 5404 ). This is recognized because the
pressure profile, that was nearly constant, starts to rise. In Figure 4.10, the green
line indicates the start-up of the machine (3120 seconds), while the time instant in
which the pressure wave reaches the final nodes (5040 seconds) is highlighted by the
red line.
The following boundary conditions are applied:
Tarvisio, Node 5394, and Node 5404 : pressure source;
Gorizia: it is modelled like a common node because the import point is inactive
in the considered time span.
In this simulation, the reference gas composition is the one obtained from the gas
chromatograph installed in the station of Malborghetto.
4.2.2 Results
The first phenomenon considered is the propagation of the pressure wave generated
by the activation of the compressor station. To see it, the measurements and the
4.2. CASE-STUDY 2: MALBORGHETTO STATION START TRANSIENT 115
Figure 4.9: The map of the considered portion of network. The names of import
points, compressor stations, and relevant nodes are reported.
Figure 4.10: Pressure profile of the nodes 5394 and 5404.

The green line, placed at 3120 seconds, indicates the time instant in which the
Malborghetto compressor station is started.
The red line, placed at 5040 seconds, indicates the time instant in which the pressure
wave reaches the nodes.
simulated value are plotted in Figure 4.11, showing the pressure at the outlet of
Malborghetto station, at the node 21349 and at the Istrana station, that is inactive.
The simulation correctly reproduces the pressure trends. In particular, the pressure
rapidly grows near the compressor station, while the farther nodes are affected by a
very small variation because the effect is delayed and smoothed by the distance.
The second interesting point to analyse is the behaviour of the station during its
start-up. In Figure 4.12, three different values of mass flow-rate are shown.
Before the activation of the gas pumping unit, the entire flow of natural gas is
carried by the bypass valve (yellow line). At 3120 seconds, the start-up transient of
the machine begins: the blue line, that represents the flow inside the compressor,
rises. At first it increases to 330 kg/s in 200 seconds, then it grows gradually for
more than 1000 seconds, until it reaches the desired value. These values are read
from the SCADA and they are set as flow-rate set-point for the station. When the
compressor is started, the quantity of gas that passes through the bypass reduces
and, when the outlet pressure becomes greater than the inlet one, the bypass valve
is closed, as stated by the description in Section 2.11.1. The orange line represents
the mass flow-rate that passes through the the station. The spike at 3120 seconds
is an artefact caused by the finite-volume approximation, whose duration is much
shorter than one second, thus irrelevant in terms of the overall mass of pumped gas.
In this case, the execution time of the entire toolchain on a laptop (Intel i7 4700HQ
2.40 GHz, RAM 16 GB, SSD for programs and data storage) is near to 3 minutes.
In Table 4.8 are reported in detail the simulation settings and the execution times.
Phase Time [s] Parameter Value

Python script 125 Start time -3600 s
Model translation 25 Stop time 5040 s
Model compilation 16 Interval 0.5 s
Model simulation 22 Number of equations 4358
Number of variables 4358
Total 188
Table 4.8: Toolchain execution times and simulation parameters in the Case-Study 2.
4.2. CASE-STUDY 2: MALBORGHETTO STATION START TRANSIENT 117
Figure 4.11: Comparison between the SCADA measurements and the simulated
values for the pressure at Malborghetto station outlet, at Node 21349, and at Istrana
station.
Figure 4.12: Behaviour of the simulated mass flow-rates that flows in Malborghetto
compressor station, in its machine and in its bypass valve.
Chapter 5
Conclusions
In this work, the first steps of the development of a complete suite for the simulation
and optimization of the Italian gas network have been done.
At first, as described in Chapter 2, the GasNetworks package, containing the Mod-

elica models required to describe the network components, has been developed.
During this phase, several simplifying assumptions have been made, but the modu-
larity guaranteed by the choice of Modelica language will be exploited to easily lift
some of them.
In second place, a Python script has been constructed. Its purpose is to aggreg-
ate the information provided by Snam in several files with different formats, then
to elaborate them to obtain a Modelica model of the entire Italian network. To
achieve this result, it creates a graph that contains the parameters of every physical
component, then it elaborates and simplifies it, with a customizable level of detail.
Finally, it uses the components in the GasNetworks package to build the model of
the network.
Finally, the models resulting from the previous procedure have been simulated in
OpenModelica, and the obtained values have been compared with the logs of the
SCADA system.
The results of this thesis work match the initial objectives, even if on limited regions
of the network: the Modelica package is able to represent a large variety of cases and
configurations with a good precision, the Python script can autonomously obtain
working models from the provided data, and the obtained values match the real
data with a good precision. In addition, the assumptions made allow to execute the
routine in reasonable times, without the requirements of high-performance machines.
119
120 CHAPTER 5. CONCLUSIONS
5.1 Future work

The results obtained in this work are far from being considered conclusive: to be able
to reach a functioning model in the time span available for the thesis development,
several approximations have been made. As explained in the previous chapters,
every simplification is justified and its impact on the final result is taken into account.
The future developments of the project, integrated into the collaboration between
Politecnico di Milano and Snam, will consider:
Improvements in the network simplification process. The two models

considered in the case studies have been obtained in a fully automated pro-
cedure, but, in the Northen Italy regions, there is still some uncertainty in the
network configuration, probably due to incomplete or incorrect data export.
In addition, the algorithm that elaborates complex and extended system, like
the considered network, must take into account a variety of cases and excep-
tions that requires to analyse a lot of real data. For example, the impact of the
consumption seasonality on the network configuration has not been considered
due to the lack of data. Furthermore, the natural gas supply from different
locations can not be considered as constant: the opening of the Trans-Adriatic
Pipeline in the last months of 2020 is the perfect example of this subject.
Variable gas composition in the network. In this stage of the develop-

ment, the gas composition of the network is defined in a single place, and it is
used in every pipe. To better represent the pressure losses caused by the gas
motion in the network, a more accurate approach should be followed. This
could also improve the conversion of the consumption distributed on the net-
work from Sm3 /h to kg/s. In fact, the first measurement unit is influenced by
the gas composition. The partial mass balance could be implemented, and the
transport and mixing of each component could be modelled. This leads to two
main problems: the number of state variables is multiplied by a factor 2-4,
depending on how many gas constituents are considered, and so the computa-
tional cost. In addition, the geographical distribution of the composition must
be estimated at the simulation start time. This requires complex estimation
algorithms and lots of data. This approach is also not compatible with MILP
optimization, because the problem becomes non-linear. The viable alternative
could consider a static map of the composition, updated every day, that states
the mixture in every point of the network, considered as constant.
Improve the temperature modelling in pipes. During the development

5.1. FUTURE WORK 121
work, it was discovered how big is the temperature effect on the behaviour of
the compressors, and, for this reason, it has been taken into account. In the
future work, a simplified modelling of the temperature in pipes could improve
the performances of the pipes. The heat exchange phenomena between the
gas in pipes and the soil can not be accurately modelled: this approach would
require an explicit model of the interaction between every type of pipe with the
terrain it is buried into, and these data are not collected by Snam. The viable
idea could be considering the temperature measured from the sensors present
in the import terminals and at the inlet of the stations to simply estimate the
temperature of the gas.
Use a more accurate model for the natural gas. In the current version
of the Modelica package, the gas model is described using a set of simplified
equations that have been tuned on the same amount of predefined gas mix-
tures. An accurate model of gas could manage all the variations supposed in
the previous points, but at an high computational cost. The improvement on
the results should be compared with the loss in terms of speed of execution of
the entire toolchain.
Advanced representation of compressor stations. In the current version

of the CompressorStation model, the pressure losses caused by the piping that
connects the gas pumping units to the network are not estimated. This ap-
proach could be refined analysing numerous working conditions and estimating
some coefficients to represent these losses as function of flow-rate and number
of active units in the station. In addition, the ideal control that is currently
used could be substituted by an explicit controller that can follow more ac-
curately the control strategy employed by Snam. In particular, it is essential
to add the operating mode whereby the rotational speed hits the maximum
value, so that none of the three set-points (flow-rate, inlet pressure, outlet
pressure) is actually followed. Finally, some compressors and gas turbines are
modelled using the maps provided by the constructor. These representations
underestimate the performances of the real appliances and lead to mismatches
between the simulations and the SCADA logs.
Storage sites. Currently, the compressors placed in storage sites are rep-
resented in a very simplified way: their consumptions are considered to be
proportional to the flow-rate. To improve the quality of the optimization res-
ults, they must be modelled with accurate equations, similar to the ones used
in the compressor stations along the network.
122 CHAPTER 5. CONCLUSIONS
Pipes spatial discretization. The pipe models used in the case studies
are split in the same number of finite volume, not considering their length.
This number has been chosen to have a reliable result when considering the
longest pipes in the network. Modelling shorter pipes with a reduced number
of elements, proportionally to their length, can improve the efficiency of the
simulation.
Bibliography
[1] ARERA. Criteri di regolazione tariffaria per il servizio di trasporto e misura

del gas naturale per il quinto periodo di regolazione (2020-2023). Deliberazione
28 marzo 2019 114/2019/R/GAS, allegato A. 2019. url: https : / / www .
arera.it/allegati/docs/19/114-19alla.pdf.
[2] Modelica Association. Language Specification. Version 3.4. 10th Apr. 2017.
Chap. 13, p. 177. url: https://www.modelica.org/documents/ModelicaSp
ec34.pdf.
[3] Modelica Association. Modelica Language. 2021. url: https://www.modelic
a.org/modelicalanguage.
[4] Francesco Casella. “On the Formulation of Steady-State Initialization Prob-
lems in Object-Oriented Models of Closed Thermo-Hydraulic Systems”. In:
Nov. 2012, pp. 215–222. isbn: 9789175198262. doi: 10.3384/ecp12076215.
[5] Francesco Casella and Alberto Leva. “Object-Oriented Modelling Simulation
of Power Plants with Modelica”. In: Proceedings of the 44th IEEE Conference
on Decision and Control. 2005, pp. 7597–7602. doi: 10 . 1109 / CDC . 2005 .
1583388.
[6] A. Cerè. Manuale Tecnico Modellistico. Snam, 2007.
[7] European Commission. COMMISSION REGULATION (EU) No 312/2014 of
26 March 2014 establishing a Network Code on Gas Balancing of Transmission
Networks. 26th Mar. 2014. url: https://eur- lex.europa.eu/eli/reg/
2014/312/oj.
[8] NetworkX Developers. Software for Complex Networks. Version 2.5. 22nd Aug.
2020. url: https://networkx.org/documentation/stable/.
[9] Celestina Dominelli. Sistema gas più flessibile, Snam lancia le centrali ibride.
Il Sole 24 Ore. 5th Feb. 2019. url: https://www.ilsole24ore.com/art/
sistema-gas-piu-flessibile-snam-lancia-centrali-ibride-AFa9iFG.
[10] Björn Geißler et al. “Mixed integer linear models for the optimization of dy-
namical transport networks”. In: Mathematical Methods of Operations Re-
123
124 BIBLIOGRAPHY
search 73.3 (June 2011), pp. 339–362. doi: 10.1007/s00186- 011- 0354- 5.
url: https://link.springer.com/article/10.1007/s00186-011-0354-5.
[11] Giulio Guandalini, Paolo Colbertaldo and Stefano Campanari. “Dynamic mod-
eling of natural gas quality within transport pipelines in presence of hydro-
gen injections”. In: Applied Energy 185 (2017). Clean, Efficient and Afford-
able Energy for a Sustainable Future, pp. 1712–1723. issn: 0306-2619. doi:
https : / / doi . org / 10 . 1016 / j . apenergy . 2016 . 03 . 006. url: http :
//www.sciencedirect.com/science/article/pii/S0306261916303178.
[12] Martin Gugat et al. “MIP-based instantaneous control of mixed-integer PDE-
constrained gas transport problems”. In: Computational Optimization and Ap-
plications 70 (May 2018), pp. 1–28. doi: 10.1007/s10589-017-9970-1.
[13] IEA. Gas 2020. 2020. url: https://www.iea.org/reports/gas-2020.
[14] IEA. Gas demand by region and scenario, 2018-2040. 2019. url: https://
www.iea.org/data- and- statistics/charts/gas- demand- by- region-
and-scenario-2018-2040.
[15] IGU, BNF and Snam. Global Gas Report 2020. 2020. url: https : / / www .
snam . it / export / sites / snam - rp / repository / file / gas _ naturale /
global-gas-report/global_gas_report_2020.pdf.
[16] Aspen Technology Inc. Aspen. 2021. url: https://www.aspentech.com/en.
[17] ISO 2314:2009. Gas turbines — Acceptance tests. Standard. Geneva, CH: In-
ternational Organization for Standardization, Dec. 2009. url: https://www.
iso.org/standard/42989.html.
[18] Kai Liu, Lorenz T. Biegler, Bingjian Zhang and Qinglin Chen. “Dynamic
optimization of natural gas pipeline networks with demand and composition
uncertainty”. In: Chemical Engineering Science 215 (2020), p. 115449. issn:
0009-2509. doi: https://doi.org/10.1016/j.ces.2019.115449. url: http
s://www.sciencedirect.com/science/article/pii/S000925091930939X.
[19] Giovanni Lozza. Turbine a gas e cicli combinati, terza edizione. Società editrice
Esculapio, 2016.
[20] Wes McKinney and the Pandas Development Team. pandas documentation.
Version 1.2.2. 9th Feb. 2021. url: https : / / pandas . pydata . org / docs /
index.html.
[21] McKinsey. The future of liquefied natural gas: Opportunities for growth. 2020.
url: https://www.mckinsey.com/industries/oil-and-gas/our-insight
s/the-future-of-liquefied-natural-gas-opportunities-for-growth.
[22] Sidhant Misra, Marc Vuffray and Anatoly Zlotnik. “Monotonicity Properties
of Physical Network Flows and Application to Robust Optimal Allocation”.
BIBLIOGRAPHY 125
In: Proceedings of the IEEE PP (Aug. 2020), pp. 1–22. doi: 10.1109/JPROC.
2020.3014069.
[23] Mohamad Mohamadi-Baghmolaei, Reza Azin, Zeinab Zarei and Shahriar Os-
fouri. “Presenting decision tree for best mixing rules and Z-factor correla-
tions and introducing novel correlation for binary mixtures”. In: Petroleum
2.3 (2016), pp. 289–295. issn: 2405-6561. doi: https://doi.org/10.1016/
j.petlm.2016.05.003. url: https://www.sciencedirect.com/science/
article/pii/S240565611630027X.
[24] Pim Nederstigt. “Real Gas Thermodynamics: and the isentropic behavior of
substances”. Master Thesis. Delft University of Technology, 2017. url: http:
//resolver.tudelft.nl/uuid:ee16f7e5-4251-4629-9192-8f4a2e3d599b.
[25] NIST. NIST Reference Fluid Thermodynamic and Transport Properties Data-
base (REFPROP). Version 10. 16th Nov. 2020. url: https://www.nist.gov/
srd/refprop.
[26] OpenModelica. .ModelicaReference.Operators.’homotopy()’. 2021. url: https
://build.openmodelica.org/Documentation/ModelicaReference.Operat
ors.%27homotopy()%27.html.
[27] OpenModelica. OpenModelica - Introduction. 2021. url: https://www.openm
odelica.org/.
[28] Andrzej J. Osiadacz and Maciej Chaczykowski. “Comparison of isothermal and
non-isothermal pipeline gas flow models”. In: Chemical Engineering Journal
81.1 (2001), pp. 41–51. issn: 1385-8947. doi: https://doi.org/10.1016/
S1385-8947(00)00194-7. url: https://www.sciencedirect.com/science
/article/pii/S1385894700001947.
[29] Pedro L. D. Peres, Ivanil S. Bonatti and Amauri Lopes. “Transmission Line
Modeling: A Circuit Theory Approach”. In: SIAM Review 40.2 (1998), pp. 347–
352. issn: 00361445. url: http://www.jstor.org/stable/2653345.
[30] Andrei Polyanin. Handbook of Linear Partial Differential Equations for En-
gineers and Scientists. Jan. 2002. isbn: 9781584882992. doi: 10.1201/97814
20035322.
[31] Jörg Rädler. Smooth transition between functions with tanh(). 2010. url: h
ttps : / / www . j - raedler . de / 2010 / 10 / smooth - transition - between -
functions-with-tanh/.
[32] Leonard Richardson. Beautiful Soup Documentation. 2020. url: https : / /
www.crummy.com/software/BeautifulSoup/bs4/doc/.
[33] Snam. Network Code. Rev. LXIV. 29th Dec. 2020. Chap. 2 - Network De-
scription and Management. url: https : / / www . snam . it / export / sites /
126 BIBLIOGRAPHY
snam-rp/repository-srg/file/en/business-services/network-code-
tariffs/Network_Code/Codice_di_Rete/64.CdR_LXIV/Rev.64_PDF/02_
network_description_management_Rev_LVIII_ENG.pdf.
[34] Snam. Snam e idrogeno. 2021. url: https://www.snam.it/it/transizione_
energetica/idrogeno/snam_e_idrogeno.
[35] Snam. Snam’s infrastructures. 15th Mar. 2019. url: https://www.snam.it/
en/about-us/snam-infrastructures/index.html.
[36] Snam. Sustainability report 2019. 2020. url: https://www.snam.it/export/
sites/snam- rp/repository/ENG_file/investor_relations/reports/
annual_reports/2019/2019_sustainability_report.pdf.
[37] Rob Story. Folium - Python data, leaflet.js maps. 2021. url: https://pytho
n-visualization.github.io/folium/index.html.
[38] Yue Qiu, Sara Grundel, Martin Stoll and Peter Benner. “Efficient numerical
methods for gas network modeling and simulation”. In: Networks & Hetero-
geneous Media 15.1556-1801-2020-4-653 (2020), p. 653. issn: 1556-1801. doi:
10 . 3934 / nhm . 2020018. url: http : / / aimsciences . org / /article / id /
6b094051-0354-42c7-9327-7d411e69895b.

02 Final

Caricato da

Copyright:

Formati disponibili

02 Final

Caricato da

Informazioni sul documento

Titolo originale

Copyright

Formati disponibili

Condividi questo documento

Condividi o incorpora il documento

Opzioni di condivisione

Hai trovato utile questo documento?

Questo contenuto è inappropriato?

Copyright:

Formati disponibili

02 Final

Caricato da

Copyright:

Formati disponibili

SCHOOL OF INDUSTRIAL AND INFORMATION ENGINEERING

Master of Science in Automation and Control Engineering

Object-Oriented Modelling and Simulation

Thesis advisor: M.Sc. Thesis of:

Academic Year 2019-2020

In questo lavoro è stato sviluppato un pacchetto Modelica che contiene i modelli

Il risultato è un pacchetto di strumenti, basati solo su software open-source, che è in

Il primo ringraziamento va ai proff. Casella e Martelli per avermi guidato in questo

A Snam va la mia gratitudine per avermi dato l’opportunità di contribuire a questo

2 Modelling of the Gas Network Components 9

3 Modelling of the Italian Gas Network 53

4 Model validation 105

1.1 Natural gas price evolution in the 2014-2020 time interval. . . . . . . 1

2.17 Example of map of a centrifugal compressor. . . . . . . . . . . . . . . 35

3.1 UML diagram describing the content of the XML file. . . . . . . . . . 54

3.14 Examples of BeautifulSoup instructions used in the program. . . . . . 61

3.42 A real example of the representation of nodes in the map. . . . . . . . 89

4.1 Map of the Central-Southern Italy network. . . . . . . . . . . . . . . 106

4.8 Comparison between the SCADA measurements and the simulated

3.1 Colours of the segments representing the different edge categories. . . 90

4.5 Toolchain execution times and simulation parameters in the Case-

1.1 Description of a natural gas transport net-

In the context of European regulations, companies that manage the NG transport

Figure 1.3: Structure of a transport network for natural gas.

The current environmental policies in Europe promote the reduction of greenhouse

1.2 Thesis goals

 Development of a script, programmed in Python, that reconstructs the network

 Simulation and validation of the model.

The Regional Network is composed of the remaining part of pipelines managed by

 2251 gate valves;

Figure 1.5: Snam network map, early 2020.

1.3 Simulation program

possible to overlook a precise comprehension and representation of the phenomena

1.4 Network data analysis and simplification

Modelling of the Gas Network

2.1 Modelica language

The main advantages of this language are: [5]

 A-causal, declarative modelling: the model of every component is composed

 Code transparency: equations can be written as they appear in textbooks, the

 Encapsulation: every model has rigorously defined connectors that allows to

 Inheritance: models can be defined with a hierarchical approach. One basic

 Multi-physics modelling: being Modelica an object-oriented language, it is

2.2 GasNetworks package

The hierarchical structure allowed by Modelica packages eases the organization of

The GasNetworks package contains nine sub-packages:

 Interfaces: the connectors used to describe the ports of the components;

 Components: models of items, used to create other models;

 Utilities: contains the custom functions created for this package.

In this thesis, the content of Interfaces, Sources, BaseClasses, Components, Media,

2.3 Methodological introduction

2.4 Gas model

pv = Z (T, p, x) R̂T (2.4)

T [K] p [bar] ρ [kg/m3 ] Z [-] cp [kJ/ (kg · K)]

The other group is represented by empirical exponential equations. These relation-

pv = Z (p, x) R̂T̄ (2.5)

CH4 C2 H6 C3 H8 C4 H10 C5 H12

In Figure 2.2, the temperature is assumed constant T̄ = 15 ◦ C. The polynomial can

Table 2.3: Coefficients evaluated for the considered mixture at T̄ = 15 ◦ C.

Figure 2.2: Linear interpolation of Z values obtained with ASPEN (T = 15 ◦ C) for

Development of a script, programmed in Python, that reconstructs the network

Simulation and validation of the model.

2251 gate valves;

A-causal, declarative modelling: the model of every component is composed

Code transparency: equations can be written as they appear in textbooks, the

Encapsulation: every model has rigorously defined connectors that allows to

Inheritance: models can be defined with a hierarchical approach. One basic

Multi-physics modelling: being Modelica an object-oriented language, it is

Interfaces: the connectors used to describe the ports of the components;

Components: models of items, used to create other models;

Utilities: contains the custom functions created for this package.

All the effort variable must assume the same value;

Newtonian fluid: the NG respects Stokes’ hypoteses;

Single-phase fluid: the NG is not affected by phase transitions;

Isothermal flow: the gas temperature is considered as fixed, as explained in

The computational cost increase caused by the increased number of states to

geometrical parameters: length, diameter, roughness, and elevation

number of finite volumes

gas consumption distribution parameters

Flow-rate profile: a guessed mass flow-rate profile is used to initialize the

No initialization: the states initial values are not specified.

Integer input (orange triangle): it is used to determine if the appliance is active

Ambient pressure: 101 325 Pa;