# Grey-Box Modelling of District Heating Networks Using Modified LPV Models

Olamilekan Tijani, Sylvain Serra, Patrick Lanusse, Rachid Malti, Hugo Viot, Jean-Michel Reneaume

► To cite this version:

Olamilekan Tijani, Sylvain Serra, Patrick Lanusse, Rachid Malti, Hugo Viot, et al. Grey-Box Modelling of District Heating Networks Using Modified LPV Models. Energies, 2025, 18 (7), pp.1626. 10.3390/en18071626. hal-05069780

HAL Id: hal-05069780

https://univ-pau.hal.science/hal-05069780v1

Submitted on 15 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Grey-Box Modelling of District Heating Networks Using Modified LPV Models

Olamilekan E. Tijani <fnref n="1" />,*, Sylvain Serra <fnref n="1" />,*, Patrick Lanusse <fnref n="2" />, Rachid Malti <fnref n="2" />, Hugo Viot <fnref n="3" />
and Jean-Michel Reneaume <fnref n="1" />

<figure><img src="image_3.png" /></figure>

check for
updates

Academic Editors: Jaroslaw Krzywanski, Marcin Sosnowski, Karolina Grabowska, Dorian Skrobek and Ghulam Moeen Uddin

Received: 21 February 2025

Revised: 14 March 2025

Accepted: 18 March 2025

Published: 24 March 2025

**Citation:** Tijani, O.E.; Serra, S.; Lanusse, P.; Malti, R.; Viot, H.; Reneaume, J.-M. Grey-Box Modelling of District Heating Networks Using Modified LPV Models. *Energies* **2025**, *18*, 1626. https://doi.org/10.3390/en18071626

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The International Energy Agency (IEA) 2023 report highlights that global energy losses have persisted over the years, with 32% of the energy supply lost in 2022 alone. To mitigate this, this research adopts optimisation to enhance the efficiency of district heating networks (DHNs), a key global energy supply technology. Given the dynamic nature of DHNs and the challenges in predicting disturbances, a dynamic real-time optimisation (DRTO) approach is proposed. However, this research does not implement DRTO; instead, it develops a fast grey-box linear parameter varying (LPV) model for future integration into the DRTO algorithm. A high-fidelity physical model replicating theoretical time delays in pipes serves as a reference for model validation. For a single pipe, the grey-box model achieved a 91.5% fit with an $R^2$ value of 0.993 and operated 5 times faster than the reference model. At the DHN scale, it captured 98.64% of the reference model's dynamics, corresponding to an $R^2$ value of 0.9997, while operating 52 times faster. Low-fidelity physical models (LFPMs) were also developed and validated, proving to be more precise and faster than the grey-box models. This research recommends performing dynamic optimisation with both models to determine which better identifies local minima.

**Keywords:** district heating networks; grey box; linear parameter varying; modelling and simulation

## 1. Introduction

### 1.1. Background of the Study

The International Energy Agency (IEA) 2023 report indicates that, since 1990, more energy has consistently been supplied to compensate for the energy losses [1]. As a result, an average of approximately 31% of the global energy supplied has been lost over the years [1]. Thus, the IEA 2023 report highlights energy efficiency enhancement as one of the key strategies for reducing global energy waste [1]. Among the various sectors using the global energy supply, approximately 50% is allocated to thermal energy generation (TEG) [1]. Therefore, improving the efficiency of the thermal energy generation sector can significantly enhance the global energy supply efficiency and reduce global energy losses. Research has shown that optimisation effectively improves efficiency [2]. Consequently, this research focuses on developing a methodology for optimising the district heating networks (DHNs) of the global thermal energy generation sector to improve the global energy supply efficiency and contribute to global energy and environmental sustainability. A district heating network (DHN) generates energy in a central station and distributes it to

<footnote n="1">¹ LaTEP, Université de Pau et des Pays de l'Adour, 64075 Pau, France; jean-michel.reneaume@univ-pau.fr</footnote>

<footnote n="2">² Université de Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, 33400 Talence, France; patrick.lanusse@bordeaux-inp.fr (P.L.); rachid.malti@u-bordeaux.fr (R.M.)</footnote>

<footnote n="3">³ Nobatek, 67 rue de Mirambeau, 64600 Anglet, France; hviot@nobatek.com</footnote>

<footnote>* Correspondence: tijani.oe@etud.univ-pau.fr (O.E.T.); sylvain.serra@univ-pau.fr (S.S.)</footnote>
consumers according to their dynamically changing power demands [3]. So, the operation of a DHN is dynamic, demanding the implementation of dynamic optimisation. Also, this network operation is influenced by unpredictable meteorological conditions [4], which are best forecasted in the short term (real time). Consequently, this paper is the first phase of the research, which aims to perform the dynamic real-time optimisation of a DHN to enhance the network efficiency.

A dynamic real-time optimisation (DRTO) algorithm performs repetitive sequential optimisations based on updated system disturbances and current state conditions [5,6]. This approach requires low-computation models to allow frequent optimisations as disturbances are predicted [7,8]. Models are broadly categorised into three types: physical (white box), black box, and grey box [9]. The physical models are high-order detailed models with high precision but are computationally expensive [9]. Thus, they may be unsuitable for real-time operation [10], particularly for complex physical systems with intricate models. Black-box models, which are based on historical data, suffer from low extrapolation capability [11]. This extrapolation capacity is indispensable for optimisation algorithms, so this research does not consider black-box models. In contrast, the grey-box models have good interpolation capability because they blend physical model structures with data, offering a balance of speed and reliability [12]. Moreover, grey-box models make it possible to estimate hidden states, which can be employed for other predictions [13]. Hence, this paper proposes a grey-box model for DHNs. While no optimisation is implemented here, the proposed grey-box model can serve as a foundation for facilitating real-time DHN optimisation. Nevertheless, the grey-box model is validated with the simulation result of a nonlinear physical DHN model to ensure its compatibility with the network operation. Thus, the principal objectives of this research are the physical and grey-box modelling of a DHN as well as the validation of the grey-box model.

## 1.2. *Literature Review on Grey-Box Models*

Several grey-box models can be formulated from the physics of a system. The research study by [10] proposed a grey-box model to estimate the hydraulic parameters of a tree-shaped DHN, aiming to ameliorate pumping efficiency during network control. In conventional DHN modelling, the hydraulic model relies on the pressure losses within the network pipes [14]. However, the research of [10] developed a grey-box hydraulic model that depends solely on the pressure measurements at two consumer substations in the network [10]. The hydraulic model is then validated through a case study simulation of a DHN, achieving an estimation measurement error of less than 1%.

A common grey-box model for thermal systems adopts the resistance-capacitance (RC) circuit analogy from electrical systems to define its structure [15]. This modelling principle is called the RC modelling [15]. The model parameters are the resistances and capacitances of the circuit, which can be identified through system identification techniques [15]. The research study by [16] developed a nonlinear grey-box model to capture the heat dynamics of a building, enabling predictions of indoor air temperature, return water temperature, and heat load. In this setup, the building receives heat from a district heating network and distributes it to its rooms via radiators. The room modelling process is divided into stages. First, a state-space model structure is formulated from an RC circuit representation of a room, expressing indoor temperature as a function of heat load. Second, an RC circuit representation is developed for the radiators in the room. These two models are then combined to derive a generalised model structure for predicting the building heat dynamics. Given the potential uncertainty in predicting building heat dynamics [16], process and measurement noises with Gaussian densities are incorporated into the deterministic model (the generalised model structure), transforming it into a stochastic nonlinear state
space model. Consequently, maximum likelihood, combined with a continuous discrete extended Kalman filter (CDEKF), is employed to identify the model parameters. The research concluded that the nonlinear grey-box model accurately predicted the indoor air temperature, return temperature, and heat load. The stochastic differential equation (SDE) used in this research was similarly developed by [13] to model the heat dynamics of a solar thermal collector. Maximum likelihood was also employed to estimate the grey-box model parameters. The research concluded that the grey-box model agreed with the EN 12975 Quasi Dynamic Test Method (CEN, European Committee for Standardisation, 2006) modelling of the same collector. This validated grey-box model was then employed in the research work of [17] to perform a nonlinear predictive control of a solar thermal system in a district heating system.

The RC modelling principle was also employed in the research work of [18] to develop a grey-box model. This grey-box model, represented by stochastic differential equations, was employed for real-time forecasting and prediction of the transformer temperatures in a power distribution network. This grey-box model accurately predicted the transformer temperatures with a mean absolute error ranging from 0.4 to 0.6 °C.

However, the main challenge in grey-box model development lies in identifying the best model structure that effectively captures the system dynamics at a low computation cost. This challenge can be addressed by explicitly capturing the nonlinear dynamics using linear parameter varying (LPV) modelling techniques [19].

A linear parameter varying (LPV) model is a linear model whose behaviour depends on time-varying parameters called the scheduling parameters or variables [20]. The research study by [19] developed a continuous-time LPV model to predict the outlet temperature of a cross-flow heat exchanger. In this model, the outlet temperature is influenced by the input (inlet temperature) and the scheduling variable (mass flow rate). The LPV model structure is formulated such that the outlet temperature depends on the inlet temperature and the model parameters, while the model parameters are functions of the scheduling variable (mass flow rate). To develop and validate this LPV model, three principal datasets were generated from the simulations of the heat exchanger physical model. The first dataset was used to train the LPV model, identifying both the model parameters and their dependency on the mass flow rate. Initial validation with the second dataset, which had constant mass flow rates, demonstrated over 90% accuracy. However, during the validation with the third dataset, which consisted of continuous mass flow rate profiles, the LPV model failed to capture the system dynamics accurately. To address this limitation, a quasi-LPV model was introduced by modifying the model structure to incorporate additional constant model parameters. This refined model expressed the heat exchanger outlet temperature as a function of the inlet temperature, mass flow rate, and both varying and constant model parameters, and it was validated against the third dataset, achieving 91.56% accuracy. Similar research was conducted by [20] to identify an LPV model that captures the dynamics of a counter-flow water-to-oil heat exchanger. The LPV model is termed a Multiple Input Single Output (MISO) quasi-LPV model [20], and it effectively models the dynamics of the oil heat exchanger under the simultaneous excitation of the inlet temperatures and mass flow rates of both water and oil.

Furthermore, LPV models have been integrated into optimisation and control algo-
rithms to leverage their speed and reliability [21]. The research work of [21] developed
a grey-box linear parameter varying (LPV) model to capture the dynamics of a gas-lift
system. The model was incorporated into a nonlinear model predictive control (NMPC) al-
gorithm and compared with a physical model-based NMPC algorithm. The results showed
that the LPV-based NMPC was three times faster than its physical model counterpart. In
addition, the research study by [22] developed an LPV model that can be integrated into
an optimisation algorithm for the optimal design of wastewater treatment plants. In this LPV modelling approach, the scheduling variables include the water and air flow rates as well as ammonium concentration. To address the local minima convergence challenge in the optimisation algorithm, a new initialisation strategy (reinitialised partial moment (RPM) filter)) was employed. The LPV model parameters were then identified through the optimisation process and subsequently validated, achieving a 65.45% fit.

## 1.3. *Research Contributions and Novelty*

Different classes of grey-box models have been developed in previous research, including the RC and the LPV modelling approaches. Some studies have incorporated process and measurement noises into deterministic models to account for system uncertainty. The research work of [23] applies the Calinski–Harabasz index to address the uncertainty in renewable energy systems, particularly their intermittency challenges. Nevertheless, uncertainty is beyond the scope of this paper. Instead, the proposed LPV grey-box model is directly validated against the simulation results of a DHN physical model.

Furthermore, while some studies have developed grey-box models for the hydraulic modelling of DHNs, to the best of our knowledge, no published research has proposed a methodology for the thermal modelling of DHNs using a modified linear parameter varying (LPV) model. Therefore, this study introduces a novel methodology for LPV-based thermal modelling of DHNs.

In contrast to the RC modelling techniques, where sizing parameters may be lost during model formulation, the proposed LPV model preserves the network sizing parameters to facilitate future DHN design optimisation. In addition, the LPV model is developed under stringent conditions to ensure robustness and computational efficiency, particularly in handling discontinuous input variables. Such discontinuities may arise from emergency system shutdowns.

## 1.4. *Research Objectives and Scope*

The objectives of this paper are to develop a physical model of a district heating network, create a grey-box model using LPV techniques, and validate the grey-box model against simulation results for both a single DHN pipe and the entire network. The scope of the paper is outlined as follows: Section 2 discusses the DHN process description, physical modelling, and grey-box model development. Section 3 presents the validation of the grey-box model and provides a discussion of the results. Finally, Section 4 concludes with insights on incorporating the grey-box model into future dynamic real-time optimisation of DHNs.

# 2. System Description and Modelling

A DHN consists of two main sections: energy generation (production) and distribution. Thermal energy in the form of hot water is generated at a central heat generation station and distributed through insulated pipes to multiple consumers based on their dynamically changing power demands, as depicted in Figure 1.

The production section begins with the main production pipe *pmp*, which directly connects to the generation station and transports hot water to the production splitter. The production splitter then distributes the hot water to branch inlets according to their designated mass flow rates.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A district heating network with a single generation station, multiple branches, and no thermal storage.</figcaption></figure>

The distribution section consists of *nb* branches, each containing *ncb* consumer components. Each consumer component has a main pipe that receives hot water and a splitter that diverts a portion of it through an extended-in pipe (also called a lateral-in pipe) *epin* to its substation to satisfy the consumer's demand. The remaining hot water continues to the next consumer component via the main pipe. After fulfilling the consumer's demand, the water cools and exits through an extended-out pipe (lateral-out) *epout* into the return pipe *rp*. This return pipe channels the cooled water to the consumer mixer, where it mixes with the branch return flow. The mixed return flows from all the branches are mixed in the production mixer before being returned to the heat-generation station via the production return pipe *prop*.

The subsequent section presents the physical and grey-box modelling of both the generation and distribution sections of a DHN.

## 2.1. *Physical Modelling of the Generation Section*

As described, the generation section comprises the production splitter, production mixer, the heat-generation station with power $P_{gen}(t)$, the production main pipe $p_{mp}$, and the production return pipe $p_{rp}$, as depicted in Figure 2. The modelling of the production pipes is addressed in Section 2.2. The physical model is developed by applying the principle of mass and energy balance to the components of the DHN.
<figure><img src="image_5.png" /><figcaption>Figure 2. The generation section of a DHN.</figcaption></figure>

Production splitter model:

$$ \text{Mass balance} \qquad \dot{m}_{pmp}(t) = \sum_{b}^{nb} \dot{m}_{b,c1,mp}(t) \qquad (1) $$

$$ \text{Temperature equality} \qquad T_{pmp,f}(t) = T_{b,c1,mp,0}(t) \quad \forall b \qquad (2) $$

Production mixer model:

$$ \text{Mass balance} \qquad \dot{m}_{prp}(t) = \sum_{b}^{nb} \dot{m}_{b,c1,rp}(t) \qquad (3) $$

$$ \text{Energy balance} \qquad T_{prp,0}(t) = \frac{\sum_{b}^{nb} \dot{m}_{b,c1,rp}(t) T_{b,c1,rp,f}(t)}{\dot{m}_{prp}(t)} \qquad (4) $$

Equation (4) assumes that the specific heat capacity of liquid water varies slowly with temperature [24]. Consequently, it is treated as constant in the modelling process. Additionally, it is considered constant for all the mixer models in this paper.

Combined heat-generation station:

$$ \text{Energy balance} \qquad P_{gen}(t) = \dot{m}_{pmp}(t)c_p(T_{pmp,0}(t) - T_{prp,f}(t)) \qquad (5) $$

## 2.2. Physical Modelling of the Distribution Section

The distribution section consists of nb branches, and each branch b has a number of consumer components, nc<sub>b</sub>. This paper models a consumer component c in a branch b for easy generalisation. As described, a consumer component c has a main pipe mp, an extended-in pipe ep<sub>in</sub>, an extended-out pipe ep<sub>out</sub>, a return pipe rp, a splitter, and a mixer, which are connected as shown in Figure 3.
<figure><img src="image_6.png" /><figcaption>Figure 3. The modelling of a consumer component c in branch b of a DHN distribution section.</figcaption></figure>

Distribution splitter model:

$$ \text{Mass balance} \qquad \dot{m}_{b,c,mp}(t) = \dot{m}_{b,c+1,mp}(t) + \dot{m}_{b,c,ep_{in}}(t) \qquad (6) $$

$$ \text{Temperature equality} \qquad T_{b,c,mp,f}(t) = T_{b,c+1,mp,0}(t) = T_{b,c,ep_{in},0}(t) \qquad (7) $$

Distribution mixer model:

$$ \text{Mass balance} \qquad \dot{m}_{b,c,rp}(t) = \dot{m}_{b,c,ep_{out}}(t) + \dot{m}_{b,c+1,rp}(t) \qquad (8) $$

$$ \text{Energy balance} \qquad T_{b,c,rp,0}(t) = \frac{\dot{m}_{b,c,ep_{out}}(t)T_{b,c,ep_{out},f}(t) + \dot{m}_{b,c+1,rp}(t)T_{b,c+1,rp,f}(t)}{\dot{m}_{b,c,rp}(t)} \qquad (9) $$

Substation (heat exchanger) model:

$$ \text{Mass balance} \qquad \dot{m}_{b,c,ep_{out}}(t) = \dot{m}_{b,c,ep_{in}}(t) \qquad (10) $$

$$ \text{Energy balance} \qquad P_{\text{demand},c}(t) = \dot{m}_{b,c,ep_{in}}(t)C_p(T_{b,c,ep_{in},f}(t) - T_{b,c,ep_{out},0}(t)) \qquad (11) $$

It can be observed that the pipes in both the distribution and generation sections have not been modelled. However, their modelling approach remains the same, irrespective of their location in the network configuration. A pipe pk, as shown in Figure 4, is modelled based on the assumptions stated in [14,25,26], which are outlined below:

* Heat transfer occurs solely in the radial direction.

* Conduction heat transfer is considered through the pipe, insulation, casing, and soil.

* Material properties are constant and independent of temperature.

* Thermal interactions between the supply and return pipes are not accounted for.

* The thermal inertia of the pipes, casing, and insulation is neglected.

* Bent pipes are treated as straight pipes of the same length.

where $\forall pk = pmp, prp, ep_{in}, ep_{out}, mp, rp$.
<figure><img src="image_7.png" /><figcaption>Figure 4. The cross-sectional view of a pipe carrying hot fluid.</figcaption></figure>

$$ \frac{\partial T(x,t)}{\partial t} + \frac{\dot{m}(t)}{\rho A} \frac{\partial T(x,t)}{\partial x} = B(t)(T_{\text{soil}}(t) - T(x,t)) \quad B(t) = \frac{\pi D_i U(t)}{\rho c_p A} \quad (12) $$

The pipe model given in (12) is a two-dimensional (2D) partial differential equation (PDE) in space and time, which is best solved by discretising one or both independent domains. In this research, the spatial temperature derivatives are approximated using second-order finite differences to balance model complexity and computation time. Specifically, the second-order central difference method is applied to approximate all the spatial temperature derivatives except the last point, as shown in (13) and (14), while the second-order backward difference method is used for the last spatial temperature derivative, as shown in (14). This spatial discretisation transforms the physical model from a PDE into a set of ordinary differential equations (ODEs).

$$ \frac{dT_i(t)}{dt} = -\left(\frac{\dot{m}(t)}{2\rho A \Delta x}\right)(T_{i+1}(t) - T_{i-1}(t)) + B(t)(T_{\text{soil}}(t) - T_i(t)) \quad \forall i = 2 \cdots n-1 \quad (13) $$

$$ \frac{dT_n(t)}{dt} = -\left(\frac{\dot{m}(t)}{2\rho A \Delta x}\right)(T_{n-2}(t) - 4T_{n-1}(t) + 3T_n(t)) + B(t)(T_{\text{soil}}(t) - T_n(t)) \quad (14) $$

Having developed the physical model of a DHN, the grey-box linear parameter varying (LPV) modelling of the network generation and distribution sections will be introduced.

## 2.3. Grey-Box Modelling of the Generation and Distribution Section

The model equations for the mixers, splitters, and substations (heat exchangers) in the generation and distribution sections are purely algebraic and, therefore, remain unchanged during the grey-box modelling. However, the primary computation challenge in the DHN physical modelling lies in the pipe models. Thus, a grey-box model is developed specifically for the pipes in the network. Since the modelling approach for all pipes is the same, a grey-box LPV model is formulated for a single pipe $pk$.

### 2.3.1. LPV Modelling of a Pipe

An LPV model is a linear model whose parameters depend on scheduling variables, which may vary over time. The principal objective of the LPV modelling technique is to
explicitly capture the system dynamics while maintaining computational efficiency. The general formulation of a single input (u), single output (y) (SISO) linear model can be mathematically shown in the time domain in (15).

$$ 
\begin{aligned}
& a_{\alpha} \frac{d^{\alpha} y(t)}{dt^{\alpha}} + a_{\alpha-1} \frac{d^{\alpha-1} y(t)}{dt^{\alpha-1}} + \cdots + a_1 \frac{dy(t)}{dt} + a_0 y(t) \\
& = b_{\beta} \frac{d^{\beta} u(t - T_d)}{dt^{\beta}} + b_{\beta-1} \frac{d^{\beta-1} u(t - T_d)}{dt^{\beta-1}} + \cdots + b_1 \frac{du(t - T_d)}{dt} + b_0 u(t - T_d)
\end{aligned}
\quad (15) $$

where $a = [a_{\alpha}, a_{\alpha-1}, \cdots, a_0]$, $b = [b_{\beta}, b_{\beta-1}, \cdots, b_{\beta}]$, and delay $T_d$ are the model parameters dependent on the scheduling variables of an LPV model. This high-order linear model can be transformed into a state-space representation or expressed in the Laplace domain, as shown in (16), to facilitate its resolution.

$$ \frac{y(s)}{u(s)} = \frac{b_{\beta}s^{\beta} + b_{\beta-1}s^{\beta-1} + \cdots + b_1s + b_0}{a_{\alpha}s^{\alpha} + a_{\alpha-1}s^{\alpha-1} + \cdots + a_1s + a_0} e^{-\tau s} = G(s)e^{-T_d s} \quad (16) $$

where $\alpha$ and $\beta$ are the number of poles and zeros, respectively.

The principal assumption of this LPV model is that the model parameters $a$, $b$, and $T_d$ must be constant with time. If this assumption is violated, a slight modification must be introduced to the LPV model to ensure its feasibility and predictive capability for continuous dynamic systems. LPV models incorporating such modifications are referred to as modified or quasi-LPV models. The steps for formulating this LPV model are now outlined.

### 2.3.2. Steps for Formulating LPV Models

The following are the steps for formulating an LPV model:

#### Step 1: Scheduling parameter or variable identification

The scheduling parameters, also known as the scheduling variables, are part of the input variables of a nonlinear physical model. When these variables remain constant, the nonlinear model becomes linear. During the development of an LPV model, these scheduling variables are no longer part of the input variables.

#### Step 2: Linear model development

Develop a linear model that shows the relationship between an output variable and the input variables at given values of the scheduling variables. This linear model development can be achieved through global and local approaches [19]. In the global approach, the dataset for the model development is obtained by exciting all the input variables simultaneously, so a Multiple Input Single Output (MISO) development is identified. In the local approach, Single Input Single Output (SISO) models are independently built and combined to form a MISO model by leveraging the superposition advantage of linear models.

#### Step 3: Generating data and identifying the relationship between model parameters and the scheduling variables

Generate a dataset that records the linear model parameters at given values of the scheduling variables. Then, machine learning, curve fitting, or system identification techniques are employed to find the relationship between the model parameters and the scheduling variables.

#### Step 4: Finalise the LPV model

Express the output variable as a function of the input variables and model parameters. Then, establish the relationship between the model parameters and the scheduling variables.
The outlined steps for formulating an LPV model are followed in this paper to develop the LPV model of a pipe.

### Step 1: Scheduling variable identification

Examining the nonlinear physical pipe model in (13) and (14), it can be observed that the input variables are the pipe inlet temperature $T_{inlet}(t)$, soil temperature $T_{soil}(t)$, and mass flow rate $\dot{m}(t)$, and the output variable is the pipe outlet temperature $T_n(t) = T_{outlet}(t)$. If the mass flow rate remains constant over time, the nonlinear model becomes linear. Therefore, the mass flow rate is the scheduling variable and is excluded from the set of input variables. Consequently, the input variables to be considered during the LPV model development of a pipe are the pipe inlet temperature and soil temperature.

### Step 2: Linear Model Development

Here, a linear model that shows the relationship between the output variable (pipe outlet temperature) and the input variables (pipe inlet temperature and soil temperature) is developed. The local approach was considered for the LPV modelling of a heat exchanger in [20] because the simultaneous long-range excitation of all the input variables of the heat exchanger is not possible in practice. Similarly, in a real-life district heating system, the simultaneous long-range excitation of all the input variables of a pipe is not feasible, so the local approach is also adopted for the linear model development in this paper.

As discussed earlier, the local approach involves developing an independent SISO model for the relationship between the output variable and each input variable. Therefore, two SISO models are developed: one showing the relationship between pipe outlet temperature $T_{outlet}(t)$ and pipe inlet temperature $T_{inlet}(t)$, represented by $T_{outlet}(t) = f(T_{inlet}(t), \theta_{inlet})$, and the other showing the relationship between the pipe outlet temperature $T_{outlet}(t)$ and soil temperature $T_{soil}(t)$, represented by $T_{outlet}(t) = f(T_{soil}(t), \theta_{soil})$. Finally, these two SISO models are combined to obtain a MISO model that shows the relationships between the pipe outlet temperature, inlet temperature, and soil temperature, as given in (17):

$$T_{outlet}(t) = f(T_{inlet}(t), T_{soil}(t), \theta) \quad (17)$$

where $\theta_{inlet}$ and $\theta_{soil}$ are sets of model parameters that depend on the scheduling variable (mass flow rate), and $\theta$ is a set of $\theta_{inlet}$ and $\theta_{soil}$.

Let us now develop the SISO models, starting with $T_{outlet}(t) = f(T_{inlet}(t), \theta_{inlet})$.

The SISO linear model that shows the relationship between $T_{outlet}(t)$ and $T_{inlet}(t)$.

The general SISO linear model given in (15) is employed to develop the linear relationship between $T_{outlet}(t)$ and pipe inlet temperature $T_{inlet}(t)$.

Relative to the physical model, pipe discretisation can also be performed during the LPV model development to find the best SISO linear model for each discretised pipe segment. Thus, the linear model for each discretised pipe model can be combined to produce a single-pipe model. In addition, the physics of fluid flow in a pipe states that a fluid element into a pipe segment will take some time before it arrives at the pipe segment outlet. This time of travel is known as the delay. This implies that a time delay not equal to zero ($T_d \neq 0$) is possible in (15). Therefore, linear models with $T_d \neq 0$ and $T_d = 0$ are tested. To minimise model complexity, the number of zeros is assumed to be zero ($\beta = 0$), and the number of poles (also known as the order of the linear model) is assumed to be less than two ($\alpha < 2$). Consequently, three case study models are possible to develop the relationship between $T_{outlet}(t)$ and $T_{inlet}(t)$. These models are shown in the Laplace domain in Figure 5, Figure 6, and Figure 7, respectively. The first-order model without delay in Figure 5 is regarded as the linearised model (model A) in this research, the zero-order
model with delay in Figure 6 is regarded as the globalised linear model (model B), and the
first-order model with delay in Figure 7 is called the distributed linear model (model C).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** LPV model for outlet temperature as a function of inlet temperature for $\beta = 0, \alpha = 1$, $T_d = 0$—Model A.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** LPV model for outlet temperature as a function of inlet temperature for $\beta = 0, \alpha = 0$, $T_d \neq 0$—Model B.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** LPV model for outlet temperature as a function of inlet temperature for $\beta = 0, \alpha = 1$, $T_d \neq 0$—Model C.</figcaption></figure>

It should be noted that the numbers of discretised pipe segments in the physical and
grey-box models are not equal: $n \neq N$.

At this stage, it is possible to identify the model parameters $k_{inlet}$, $\tau_{inlet}$, $T_d$, $\tau_{in}$, and $T_{d,in}$ since the model structures have been formulated. These model parameters are members of the set $\theta_{inlet}$. According to the LPV modelling steps, these model parameters should be obtained at given values of the scheduling variables (mass flow rate) bounded by $\dot{m}_{min} < \dot{m} < \dot{m}_{max}$ by the exciting input variable ($T_{inlet}$) for all the three possible model structures. However, this could be somewhat cumbersome. As a result, this research decided to obtain the values of some model parameters from the physical model analytically. Firstly, the spatial temperature derivative in the physical pipe model in (12) is approximated with the first-order backward difference formula, as shown in (18), with the purpose of developing a reduced order model. Laplace transformation is then imposed on the physical model at constant soil temperature to yield (19). Finally, (19) is rearranged to yield (20), which is then compared with model A (first-order LPV model without delay) in Figure 5 to identify some model parameters ($k_{inlet}$ and $\tau_{inlet}$). The time delay $T_d$ can be computed
using the mass continuity equation. Therefore, three out of the five model parameters ($k_{inlet}$, $\tau_{inlet}$, $T_d$) have been successfully identified, as shown in (21), (22), and (23).

$$ \frac{dT_i(t)}{dt} = -\left(\frac{\dot{m}(t)}{\rho A \Delta x} + B(t)\right) T_i(t) + B(t) T_{soil}(t) + \frac{\dot{m}(t)}{\rho A \Delta x} T_{i-1}(t) \quad (18) $$

$$ sT_i(s) - T_i(0) = -\left(\frac{\dot{m}}{\rho A \Delta x} + B\right) T_i(s) + \frac{\dot{m}}{\rho A \Delta x} T_{i-1}(s) \quad (19) $$

$$ T_i(s) = \frac{\frac{\dot{m}}{\rho A \Delta x} / \left( \frac{\dot{m}}{\rho A \Delta x} + B \right)}{\frac{1}{\left( \frac{\dot{m}}{\rho A \Delta x} + B \right)^s} + 1} T_{i-1}(s) \quad (20) $$

Therefore, the model parameters are

$$ k_{inlet} = \frac{\frac{\dot{m}}{\rho A \Delta x}}{\frac{\dot{m}}{\rho A \Delta x} + B} \quad (21) $$

$$ T_d = \frac{\rho A \Delta x}{\dot{m}(t)} \quad (22) $$

$$ \tau_{inlet} = \frac{1}{\frac{\dot{m}}{\rho A \Delta x} + B} \quad (23) $$

$$ N = \frac{L}{\Delta x} \quad (24) $$

This model identification technique implies that model A is similar to the first-order physical pipe model, which implicitly accounts for the delay. While model B is a delayed model without any time-constant dynamics, model C incorporates time-constant ($\tau_{in}$) dynamics and explicitly accounts for the delay with $T_{d,in}$. So, $T_{d,in} < T_d$, because $T_d$ represents the complete delay in a pipe segment, part of which has already been provided by $\tau_{in}$. Therefore, $T_{d,in}$ and $\tau_{in}$ can be computed using $T_{d,in} = F_1 T_d$ and $\tau_{in} = F_2 T_d$, where $F_1$ and $F_2$ are constants to be determined. Thus, the identification of the $T_{d,in}$ and $\tau_{in}$ for model C depends on the values of $F_1$ and $F_2$.

$F_1$ and $F_2$ are determined using system identification techniques. First, a dataset is generated by simulating the nonlinear physical pipe model with given profiles of the input and scheduling variables. A pseudo-step change, given in (30), is applied to the pipe inlet temperature $T_{inlet}$ at a specified step-change time, while a step change in mass flow rate is imposed at a later time, and a constant soil temperature is considered. The pipe outlet temperature profile obtained from the simulation of the physical model is regarded as the actual pipe outlet temperature. Model C is then implemented in Simulink with initial guess values for $F_1$ and $F_2$. A grid search optimisation algorithm is used to obtain the optimal values of $F_1$ and $F_2$ that best minimise the objective function. The objective function is defined as the sum of the relative errors between the actual pipe outlet temperature and predicted pipe outlet temperature from the Simulink model. The optimal values of $F_1$ and $F_2$ obtained from the optimisation process are then used to identify $T_{d,in}$ and $\tau_{in}$. As a result, the set of model parameters $\theta_{inlet}$, which defines the relationship between the pipe outlet temperature $T_{outlet}(t)$ and pipe inlet temperature $T_{inlet}(t)$, has been successfully identified.

Figures 5–7 can be then transformed into mathematical formulas to derive the relationship between the outlet temperature and the inlet temperature for models A, B, and C, as shown in (25), (26), and (27), respectively.
$$T_{\text{outlet}}^{\text{inlet}}(s) = \left[ \frac{k_{\text{inlet}}}{\tau_{\text{inlet}s} + 1} \right]^N T_{\text{inlet}}(s) \qquad (25)$$

$$T_{\text{outlet}}^{\text{inlet}}(s) = \left[ k_{\text{inlet}}e^{-T_d} \right]^N T_{\text{inlet}}(s) \qquad (26)$$

$$T_{\text{outlet}}^{\text{inlet}}(s) = \left[ \frac{k_{\text{inlet}}}{\tau_{\text{inlet}}s + 1} e^{-T_{d,in}} \right]^N T_{\text{inlet}}(s) \qquad (27)$$

where $T_{\text{outlet}}^{\text{inlet}}$ denotes the influence of the inlet temperature on the pipe outlet temperature.

Having completed the model structure development and parameter identification, the best model structure among models A, B, and C can be selected through model evaluation using the coefficient of determination ($R^2$) and percentage fit formulas given in (28) and (29), respectively, as the evaluation metrics. The simulation result of each grey-box model (A, B, and C) is validated against that of the physical model of a single pipe. Consider a first case study where a non-adiabatic pipe of 0.31 m diameter, 200 m length, having thermophysical properties listed in Appendix A.2, transmits hot water at the inlet temperature described in (30). The inlet temperature undergoes a pseudo-step change at 60 s from 90 °C to 95 °C. The pipe operates with a constant mass flow rate 80 kg/s and is influenced by a constant soil temperature of 0 °C. A 0.5 m length pipe segment is used for the physical model spatial discretisation, ensuring that the delay from the simulation result matches the physical delay calculated using the mass continuity equation. The simulation runs for 2400 s, and the physical model is simulated in MATLAB with ode15s, a variable step-variable order (VSVO) solver employing the numerical differentiation formula (NDF). The absolute and relative tolerances are set to $1 \times 10^{-10}$ and $1 \times 10^{-9}$, respectively. The N values of the grey-box models, sufficient to capture the pipe dynamics without numerical problems, are 100, 180 (calculated using a Courant number of 1 for explicit models), and 28 for models A, B, and C, respectively. In this example, the general hypothesis of an LPV model is not violated because the scheduling variable (mass flow rate) is constant with time. Therefore, models A, B, and C are used as linear parameter-varying models without any modification. These models can be directly simulated in MATLAB using the 'lsim' function by declaring them as transfer function models. After model resolution, the pipe outlet temperatures obtained from simulating the physical model and the three-case study LPV models are illustrated in Figure 8.

$$R^2_{pk} = 1 - \frac{\sum_{ti}^{nti} (T_{pk}(t_{ti}) - \hat{T}_{pk}(t_{ti}))^2}{\sum_{ti}^{nti} (T_{pk}(t_{ti}) - \bar{T}_{pk})^2} \qquad (28)$$

$$\% fit_{pk} = 1 - \sqrt{\frac{\sum_{ti}^{nti} (T_{pk}(t_{ti}) - \hat{T}_{pk}(t_{ti}))^2}{\sum_{ti}^{nti} (T_{pk}(t_{ti}) - \bar{T}_{pk})^2}} \qquad (29)$$

$$T_{\text{inlet}}(t) = T_{\text{inlet\_start}} + \frac{T_{\text{inlet\_step}}}{1 + \exp\left(-\frac{(t - t_{\text{step\_change}})}{3.6}\right)} \qquad (30)$$

where $T_{\text{inlet\_start}}$ is the initial value of the inlet temperature, $T_{\text{inlet\_step}}$ is the step change of the inlet temperature, and $t_{\text{step\_change}}$ is the time at which the step change is imposed on the inlet temperature.
<figure><img src="image_11.png" /><figcaption>Figure 8. Outlet pipe temperatures for inlet temperature variation.</figcaption></figure>

It is important to note that the percentage fits and R² values are calculated over the region where significant changes occur in the input variables and where these changes have a measurable impact on the output variable. Therefore, it can be observed that models A, B, and C fit 84.4%, 99.51%, and 90.63% of the physical model results, respectively, with R² values of 0.976, 1, and 0.991. Since the percentage fit and R² of model A are considerably lower than those of models B and C, it can be eliminated from consideration. The best model between models B and C can be selected by exciting the mass flow rate from 40 kg/s abruptly to 80 kg/s at 1200 s, along with the pseudo-step changed inlet temperature in the first case study at a constant soil temperature of 0 °C. In this case study, the LPV modelling assumption has been violated because the variation in the scheduling variable (mass flow rate) causes the model parameters to change with time. Therefore, models B and C are employed as modified LPV models. The thermal conductivities of the insulation and casing in Appendix A.2 were changed to that of the pipe material to increase the heat losses and assess the model’s efficiency under stringent conditions. Since varying the model parameters of a Laplace model is more effortless in Simulink, the simulation was transferred to Simulink. Thus, the Simulink representation of model B can be found in Figure 9. In addition, to allow the possibility of changing the number of discretisation N for any pipe length, model C is converted to a delayed-state space model ((31) to (36)), with the Simulink representation in Figure 10.

$$ \tau_{in} T'_{gb}(t) = AT_{gb}(t) + BT_{inlet}(t) + aT_{gb}(t - T_{d,in}) + bT_{inlet}(t - T_{d,in}) \quad (31) $$

$$ T_{outlet}^{inlet}(t) = CT_{gb}(t) \quad (32) $$

where $T_{gb} \in R^{N \times 1}$ is a vector of the discretised temperature that shows the influence of the inlet temperature on the pipe outlet temperature; A, a, B, b, C are the modified LPV model parameters, $B = 0$, $B \in R^{N \times 1}$.

$$ A = -I, \quad A \in R^{N \times N} \quad (33) $$
$$a = k_{\text{inlet}} \begin{bmatrix} 0 & 0 & 0 & \dots & 0 \\ 1 & 0 & 0 & \dots & 0 \\ 0 & 1 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix}, a \in \mathbb{R}^{N \times N} \qquad (34)$$

$$b = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix}, b \in \mathbb{R}^{N \times 1} \qquad (35)$$

$$C = \begin{bmatrix} 0 & \dots & 0 & 1 \end{bmatrix}, C \in \mathbb{R}^{1 \times N} \qquad (36)$$

<figure><img src="image_12.png" /><figcaption>Figure 9. Simulink representation of model B.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10. Simulink representation of model C.</figcaption></figure>

The pipe outlet temperature for this second case study is illustrated in Figure 11 for models B and C. The R² values and percentage fit of models B and C at the region of mass flow rate variation only are (−1.053, −43.27%) and (0.996, 93.63%), respectively.
This implies that model B could not capture the dynamics induced by the mass flow rate variation. Therefore, model C is the best LPV model, as it effectively demonstrates the relationship between the pipe inlet and outlet temperatures. This research proceeds to develop the LPV model that shows the relationship between the pipe outlet temperature and the soil temperature. The transfer function model in (27) and the state space models in ((31) to (36)) represent the relationship between the inlet temperature and the pipe outlet temperature.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Outlet pipe temperatures for inlet temperature variation and mass flow rate excitation.</figcaption></figure>

The SISO linear model that shows the relationship between $T_{\text{outlet}}(t)$ and $T_{\text{soil}}(t)$.

Relying on the physical knowledge of a pipe operation, the soil temperature influence on the pipe is distributed over the pipe length. Therefore, the LPV model structure is formulated as a distributed delay model, as shown in Figure 12.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Distributed delay for the soil temperature influence on the pipe outlet temperature.</figcaption></figure>
Figure 12 shows that this SISO (soil temperature influence) model depends on both the inlet and soil temperature influences. Therefore, these influences are captured separately.

Although model C has been identified as the best LPV model for describing the relationship between pipe outlet and inlet temperature, it is also important to assess the feasibility of using model B to capture the dynamics of the inlet temperature influence in this SISO model. The reason is that model B is less complex than model C, and it could leverage the distribution delay effect to account for the dynamics of the inlet temperature influence. Consequently, two mathematical models are developed for the SISO model that demonstrate the relationship between the pipe outlet and soil temperature. Thus, these models are referred to as soil influence model B (37) and soil influence model C (39).

For soil influence model B, Figure 12 can be translated to a mathematical equation in (37):

$$ \frac{T_{\text{soil}}^{\text{soil}}(s)}{T_{\text{soil}}(s)} = \frac{k_{\text{soil}}}{1 + \tau_{\text{soil}}s} \sum_{i=0}^{N-1} k_{\text{inlet}}^i e^{-iT_{\text{soil}}s} \quad (37) $$

$\sum_{i=0}^{N-1} (k_{\text{inlet}})^i e^{-iT_{\text{soil}}s}$ corresponds to the sum of a geometric series with a common ratio of $k_{\text{inlet}}e^{-T_{\text{soil}}s}$, and first term equals 1. This term can be replaced with the general formula for summing a geometric series to give (38):

$$ \frac{T_{\text{soil}}^{\text{soil}}(s)}{T_{\text{soil}}(s)} = \frac{k_{\text{soil}}}{1 + \tau_{\text{soil}}s} \frac{1 - (k_{\text{inlet}})^N e^{-NT_{\text{soil}}s}}{1 - k_{\text{inlet}}e^{-T_{\text{soil}}s}} \quad (38) $$

Also, the soil influence model C can be formulated from Figure 12, as expressed in (39),

$$ \frac{T_{\text{soil}}^{\text{soil}}(s)}{T_{\text{soil}}(s)} = \frac{k_{\text{soil}}}{1 + \tau_{\text{soil}}s} \sum_{i=0}^{N-1} \left( \frac{k_{\text{inlet}}}{1 + \tau_{\text{in}}s} \right)^i e^{-iT_{\text{soil}}s} \quad (39) $$

and its transformation yields (40):

$$ \frac{T_{\text{soil}}^{\text{soil}}(s)}{T_{\text{soil}}(s)} = \frac{k_{\text{soil}}}{1 + \tau_{\text{soil}}s} \frac{1 - \left( \frac{k_{\text{inlet}}}{1 + \tau_{\text{in}}s} \right)^N e^{-NT_{\text{soil}}s}}{1 - \frac{k_{\text{inlet}}}{1 + \tau_{\text{in}}s} e^{-T_{\text{soil}}s}} \quad (40) $$

Having developed the model structures, it is imperative to find the relationship between model parameters $k_{\text{soil}}$ and $\tau_{\text{soil}}$ and the scheduling parameter (mass flow rate). This relationship can be obtained by applying Laplace transformation to (18) at a constant pipe inlet temperature to yield (41):

$$ sT_i(s) - T_i(0) = -\left(\frac{\dot{m}}{\rho A \Delta x} + B\right) T_i(s) + BT_{\text{soil}}(s) \quad (41) $$

$$ T_i(s) = \frac{B/(\frac{\dot{m}}{\rho A \Delta x} + B)}{\frac{1}{(\frac{\dot{m}}{\rho A \Delta x} + B)^s} + 1} T_{\text{soil}}(s) \quad (42) $$

$$ \text{Therefore, } k_{\text{soil}} = \frac{B}{\frac{\dot{m}}{\rho A \Delta x} + B} \text{ and } \tau_{\text{soil}} = 1/(\frac{\dot{m}}{\rho A \Delta x} + B). $$

### Step 3: Generating the dataset and identifying the relationship between model parameters and scheduling variables

During the development of the SISO linear models, two classes of model parameters were identified:
* Parameters that vary with the scheduling variables (mass flow rate), which include $k_{inlet}$, $k_{soil}$, $\tau_{inlet}$, $\tau_{soil}$, $T_d$, $\tau_{in}$, and $T_{d,in}$.

* Parameters that remain constant with the scheduling variable, which are $F_1$ and $F_2$.

Since the relationships between the model parameters and the mass flow rate have been directly identified from the physical knowledge of the system, implementing this step is unnecessary.

### Step 4: LPV Model Finalisation

Here, the developed LPV SISO models can be combined into a single LPV MISO model using the superposition principle:

$$T_{outlet}(s) = T_{outlet}^{inlet}(s) + T_{outlet}^{soil}(s) \quad (43)$$

As the optimal soil influence model is yet to be determined, two pipe models are compared: model CC (44), which uses model C for both inlet and soil temperature influences, and model CB (45), which uses model C for the inlet temperature influence and model B for the soil temperature influence.

$$T_{outlet}(s) = \left[ \frac{k_{inlet}(\dot{m})}{1 + \tau_{in}(\dot{m})s} e^{-T_{d,in}(\dot{m})s} \right]^N T_{inlet}(s) + \frac{k_{soil}(\dot{m})}{1 + \tau_{soil}(\dot{m})s} \frac{1 - \left( \frac{k_{inlet}(\dot{m})}{1 + \tau_{in}(\dot{m})s} \right)^N e^{-NT_{d,in}(\dot{m})s}}}{1 - \frac{k_{inlet}(\dot{m})}{1 + \tau_{in}(\dot{m})s} e^{-T_{d,in}(\dot{m})s}} T_{soil}(s) \quad (44)$$

$$T_{outlet}(s) = \left[ \frac{k_{inlet}(\dot{m})}{1 + \tau_{in}(\dot{m})s} e^{-T_{d,in}(\dot{m})s} \right]^N T_{inlet}(s) + \frac{k_{soil}(\dot{m})}{1 + \tau_{soil}(\dot{m})s} \frac{1 - (k_{inlet}(\dot{m}))^N e^{-NT_d(\dot{m})s}}{1 - k_{inlet}(\dot{m})e^{-T_d(\dot{m})s}} T_{soil}(s) \quad (45)$$

A third case study is conducted to determine the optimal single-pipe model between models CB and CC. The pipe inlet temperature undergoes a pseudo-step change from 90 °C to 95 °C at 60 s using the sigmoid function in (30), and the soil temperature and mass flow rate undergo a step change from -5 °C to 35 °C at 600 s and 40 kg/s to 90 kg/s at 1200 s, respectively. Also, the thermal conductivity of the pipe material is assumed to be the same as that of the insulation and casing. The simulation spans from 0 to 2400 s, using $N = 28$ for the grey-box models and $\Delta x = 0.5$ m for the physical model. The grey-box model is solved using a fixed-step ode3 solver, based on the Bogacki-Shampine algorithm, using a fixed time step of 1 s. The computation times for the physical model and models CB and CC are 1.86 s, 0.38 s, and 0.39 s, respectively, conducted on a Windows 10 DELL PC with 32 GB RAM, a 12th Gen Intel(R) Core (TM) i7-12700H processor, and a clock speed of 2.30 GHz. After simulation, the pipe outlet temperatures obtained from the physical model and models CB and CC are shown in Figure 13, with both BC and CC achieving an overall average $R^2$ value of 0.993 and a percentage fit of 91.5% against the physical model.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Pipe outlet temperature for input variable variations.</figcaption></figure>

Since models CB and CC have similar computation times and percentage fits against the physical model, either can serve as the grey-box single-pipe model. However, the less complex model, in terms of the model order, is recommended for optimisation tasks. Thus, model CB is the preferred grey-box pipe model for predicting or estimating a pipe outlet temperature as a function of the inlet temperature, soil temperature, mass flow rates, and thermophysical properties. A large step change was imposed on the soil temperature and mass flow rate, despite such a change being unrealistic in a real DHN, to explicitly select the best model under stringent conditions. However, model CB will always be preferred over model CC due to its simplicity and robustness.

Since the grey-box pipe model has been validated, it can be used to simulate a DHN to test its feasibility with the network operation.

**3. Model Validation Using DHN Simulation**

To simulate a DHN, it is imperative to perform a degree of freedom (DOF) analysis to determine the number of variables that must be specified for the simulation. A presumed existing DHN is simulated in this research. Since the network is assumed to be pre-existing, the pipe dimensions, thermophysical properties of the pipe and water, and the pipe insulation properties are considered known parameters. The distribution pipes are represented as $k = \{mp, ep_{in}, ep_{out}, rp\}$, meaning that $nk = 4$ for all DHNs, while the production pipes are represented as $p = \{pmp, prp\}$.

Thus, the degree of freedom can be calculated using the information in Table 1 as
$$DOF = N_{variables} - N_{equations} = 2nkc + 2 - nc - 5nb - \sum_{b=1}^{nb} 5(nc_b - 1)$$
for an existing DHN with a central heat-generation station.
<table><caption>**Table 1.** Degree of freedom (DOF) analysis of an existing DHN with a central generation station.</caption><thead><tr><th>Variables</th><th>No. of Variables</th><th>Equations</th><th>No. of Equations</th></tr></thead><tbody><tr><td>T<sub>b,c,k,i</sub></td><td>∑<sub>b=1</sub><sup>nb</sup> ∑<sub>c=1</sub><sup>nc<sub>b</sub></sup> ∑<sub>k=1</sub><sup>nk</sup> nd<sub>b,c,k</sub></td><td rowspan="3">Discretised pipe equations</td><td>∑<sub>b=1</sub><sup>nb</sup> ∑<sub>c=1</sub><sup>nc<sub>b</sub></sup> ∑<sub>k=1</sub><sup>nk</sup> nd<sub>b,c,k</sub> - 1</td></tr><tr><td>T<sub>p,i</sub></td><td>nd<sub>pmp</sub> + nd<sub>prp</sub></td><td>nd<sub>pmp</sub> + nd<sub>prp</sub> - 2</td></tr><tr><td>T<sub>soil</sub></td><td>1</td><td></td></tr><tr><td>P<sub>demand,c</sub></td><td>nc</td><td>Substation</td><td>2nc</td></tr><tr><td>ṁ<sub>b,c,k</sub></td><td>nk × nc</td><td>Distribution splitters</td><td>∑<sub>b=1</sub><sup>n<sub>b</sub></sup> 3(nc<sub>b</sub> - 1)</td></tr><tr><td rowspan="3">ṁ<sub>p</sub></td><td rowspan="3">2</td><td>Splitter continuity</td><td>2nb</td></tr><tr><td>Distribution mixers</td><td>∑<sub>b=1</sub><sup>n<sub>b</sub></sup> 2(nc<sub>b</sub> - 1)</td></tr><tr><td>Mixer continuity</td><td>2nb</td></tr><tr><td rowspan="2">P<sub>gen</sub></td><td rowspan="2">1</td><td>Production splitter</td><td>nb + 1</td></tr><tr><td>Production mixer</td><td>2</td></tr><tr><td>N<sub>variables</sub></td><td>{∑<sub>b=1</sub><sup>nb</sup> ∑<sub>c=1</sub><sup>nc<sub>b</sub></sup> ∑<sub>k=1</sub><sup>nk</sup> nd<sub>b,c,k</sub>} +<br>nd<sub>pmp</sub> + nd<sub>prp</sub> +<br>(nk + 1)nc + 4</td><td>N<sub>equations</sub></td><td>{∑<sub>b=1</sub><sup>nb</sup> ∑<sub>c=1</sub><sup>nc<sub>b</sub></sup> ∑<sub>k=1</sub><sup>nk</sup> nd<sub>b,c,k</sub> - 1 + ∑<sub>b=1</sub><sup>nb</sup> 5(nc<sub>b</sub> - 1)} +<br>nd<sub>pmp</sub> + nd<sub>prp</sub> + 2nc + 5nb + 2</td></tr></tbody></table>

Let us perform a case study simulation of an existing DHN with two branches $nb = 2$, five consumers in each branch $nc_1 = 5$, $nc_2 = 5$. So, $nc = 10$. The pipe diameters and lengths in the network are given in Table 2.

<table><caption>**Table 2.** Pipe diameters and lengths for the case study DHN simulation.</caption><thead><tr><th>Consumer</th><th>Φ<sub>mp</sub> (in)</th><th>L<sub>mp</sub> (m)</th><th>Φ<sub>ep</sub> (in)</th><th>L<sub>ep</sub> (m)</th></tr></thead><tbody><tr><td>C1</td><td>8.00</td><td>279.93</td><td>5.00</td><td>95.21</td></tr><tr><td>C2</td><td>6.00</td><td>720.07</td><td>4.00</td><td>150.00</td></tr><tr><td>C3</td><td>5.00</td><td>176.46</td><td>2.00</td><td>40.00</td></tr><tr><td>C4</td><td>5.00</td><td>124.69</td><td>3.50</td><td>30.00</td></tr><tr><td>C5</td><td>1.50</td><td>397.32</td><td>1.50</td><td>100.00</td></tr><tr><td>C6</td><td>10.00</td><td>190.01</td><td>5.00</td><td>25.00</td></tr><tr><td>C7</td><td>8.00</td><td>97.04</td><td>3.50</td><td>10.00</td></tr><tr><td>C8</td><td>6.00</td><td>268.79</td><td>3.00</td><td>90.20</td></tr><tr><td>C9</td><td>5.00</td><td>1280.60</td><td>3.00</td><td>200.00</td></tr><tr><td>C10</td><td>5.00</td><td>198.28</td><td>4.00</td><td>50.00</td></tr><tr><td>Production</td><td>12.00</td><td>200.00</td><td></td><td></td></tr></tbody></table>

The DOF of the case study DHN simulation is $DOF = 22$. This implies that twenty-two variables must be specified to simulate the DHN. The profiles of these variables must be defined as functions of time. This research uses the supply mass flow rates $ṁ<sub>c,k</sub>$ (10 variables) of consumer component c, the power demands $P<sub>demand,c</sub>$ (10 variables) of consumer component c, the soil temperature $T<sub>soil</sub>$ (1 variable), and the plant generation
temperature $T_{pmp,0}$ (1 variable). These declared variables are the simulation input variables, and their time profiles are shown in Figures 14–16. The power demand profile and soil temperature were sourced from [26,27], respectively. The mass flow rate, adapted from [27], was modified to stay within the pipes’ maximum allowable limits to prevent erosion. A second-degree polynomial function was arbitrarily chosen for plant generation temperature to test the grey-box model under varying plant temperature conditions.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Supply mass flow rates of the consumers.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Plant generation and soil temperatures.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Power demands of the consumers.</figcaption></figure>
The simulation environment consists of MATLAB for the physical model and Simulink for the grey-box model, covering a duration of 1 to 24 h. The pipes in the DHN are discretised using an initial discretised segment length of $\Delta x_{initial} = 0.5$ m and $\Delta x_{initial} = 8$ m for the physical and grey-box models, respectively. The number of discretised segments of each pipe is initially calculated with the $\Delta x_{initial}$ using the piecewise functions in (46) and (47) for the physical and grey-box model simulations, respectively. Different piecewise functions are employed due to the distinct discretisation techniques used in the physical and grey-box models. In the grey-box model, the pipe inlet temperatures are treated directly as input variables, necessitating at least one discretised point to estimate the pipe outlet temperatures. Conversely, a second-order discretisation is applied to the spatial domain of the physical model, requiring a minimum of three discretisation points (including the inlet temperature) to estimate the pipe outlet temperatures. To ensure that a pipe is evenly discretised along its length, the discretised segment length of a pipe $pk$ to be used in the simulation is recalculated using $\Delta x_{pk} = \frac{L_{pk}}{n_{pk}-1}$ and $\Delta x_{pk} = \frac{L_{pk}}{N_{pk}}$ for the physical and grey-box model simulations, respectively. Consequently, $\Delta x_{pk}$ varies with the pipe lengths. However, the initially assigned segment lengths $\Delta x_{initial}$ are still referenced when describing the DHN simulation setup.

$$n_{pk} = \begin{cases} 3 & \text{if } \text{round}\left(\frac{L_{pk}}{\Delta x_{initial}}\right) + 1 < 3 \\ \text{round}\left(\frac{L_{pk}}{\Delta x_{initial}}\right) + 1 & \text{(includes pipe inlet temperature), otherwise} \end{cases} \quad (46)$$

$$N_{pk} = \begin{cases} 1 & \text{if } \text{round}\left(\frac{L_{pk}}{\Delta x_{initial}}\right) = 0 \\ \text{round}\left(\frac{L_{pk}}{\Delta x_{initial}}\right) & \text{otherwise} \end{cases} \quad (47)$$

The physical model uses the ode15s solver in MATLAB, a variable-step, variable-order algorithm, with a maximum time step of 10 s and a relative and absolute tolerance of $1 \times 10^{-4}$ and $1 \times 10^{-6}$, respectively. Higher maximum time steps led to unrealistic simulation results, so the maximum time step was constrained to 10 s. Also, the Jacobian sparsity was supplied to the ode15s solver to hasten the simulation. The grey-box model in Simulink uses the fixed-step ode3 solver, based on the Bogacki–Shampine algorithm, with a fixed time step of 1 s.

The network temperatures after simulation are shown in Figures 17–24. The physical model simulation results are represented with continuous lines, while the grey-box model simulation results are represented with asterisks. The DHN grey-box model achieves minimum and average $R^2_{DHN}$ values of 0.9983 and 0.9997 using (48) and (49), with corresponding minimum and average percentage fits of 95.92% and 98.64% using (50) and (51), respectively, when compared with the physical model. These results demonstrate the robustness of the grey-box LPV pipe model in capturing pipe dynamics without numerical issues.

$$\text{minimum } R^2_{DHN} = \min\{R^2_1, R^2_2, \dots, R^2_{npk}\} \quad (48)$$

$$\text{average } R^2_{DHN} = \frac{\sum_{pk}^{npk} R^2_{pk}}{npk} \quad (49)$$

$$\text{minimum \% fit}_{DHN} = \min\{\% \text{ fit}_1, \% \text{ fit}_2, \dots, \% \text{ fit}_{npk}\} \quad (50)$$

$$\text{average \% fit}_{DHN} = \frac{\sum_{pk}^{npk} \% \text{ fit}_{pk}}{npk} \quad (51)$$
<figure><img src="image_20.png" /><figcaption>Figure 17. Main pipe outlet temperatures for consumers 1–5 and production.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 18. Main pipe outlet temperatures for consumers 6–10.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 19. Extended-in pipe outlet temperatures for consumers 1–5.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 20. Extended-in pipe outlet temperatures for consumers 6–10.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 21. Extended-out pipe outlet temperatures for consumers 1–5.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 22. Extended-out pipe outlet temperatures for consumers 6–10.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 23. Return pipe outlet temperatures for consumers 1–5.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 24. Return pipe outlet temperatures for consumers 6–10 and production.</figcaption></figure>

Simulating the DHN with the physical model took approximately 165.9 min for a pipe segment length $\Delta x$ of 0.5 m with a maximum Courant number of 30.5 compared with 3.2 min with the grey-box model, which had a maximum Courant number of 0.19 and a maximum integration time step of 1 s, making the grey-box model about 52 times faster.

The DHN model validation results may seem counterintuitive, as the average percentage fit is higher for the DHN (with multiple pipes) case than for the single-pipe model. However, this discrepancy arises from the differing validation conditions. The single-pipe model was tested under strict conditions, including step changes in mass flow rate (50 kg/s) and soil temperature (40 °C), along with a pseudo-step change of 5 °C in inlet temperature. In contrast, the DHN validation involved smoothly varying profiles for input variables such as supply mass flow rates, power demands, plant temperature, and soil temperature.

Up to this point in the paper, the accuracy of the model simulations has been assessed based on the physical time delay, which is the time it takes for a fluid particle to travel from the pipe inlet to the outlet. This rationale led to the simulation of the physical and grey-box models using discretised pipe segment lengths $\Delta x_{intial}$ of 0.5 m and 7.14 m, respectively. However, a key question arises: could larger discretised pipe segment lengths achieve nearly the same accuracy as the high-fidelity physical model (simulated with $\Delta x_{intial}$ = 0.5 m)? To investigate this, the physical and grey-box models are simulated at several $\Delta x_{intial}$ values and compared with the high-fidelity physical model, as shown in
Tables 3 and 4, using the $R^2_{DHN}$ in (48) and (49) and the percentage fit formula in (50) and (51).

<table><caption>**Table 3.** Computation time and fits for physical model simulated at varying $\Delta x$ values.</caption><thead><tr><th>Δx<sub>initial</sub> (m)</th><th>1</th><th>2</th><th>5</th><th>8</th><th>10</th><th>15</th></tr></thead><tbody><tr><td>Computation time (s)</td><td>24.03 × 60</td><td>3.72 × 60</td><td>29.92</td><td>8.36</td><td>4.31</td><td>2.26</td></tr><tr><td>Minimum % fit<sub>DHN</sub> (%)</td><td>99.97</td><td>99.91</td><td>99.83</td><td>99.7</td><td>99.66</td><td>99.50</td></tr><tr><td>Minimum R<sup>2</sup><sub>DHN</sub></td><td>1.0000</td><td>1.0000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td>Average % fit<sub>DHN</sub> (%)</td><td>99.99</td><td>99.98</td><td>99.95</td><td>99.92</td><td>99.90</td><td>99.86</td></tr><tr><td>Average R<sup>2</sup><sub>DHN</sub></td><td>1.0000</td><td>1.0000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td>Maximum integration time step (s)</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>Maximum Co</td><td>15.25</td><td>7.71</td><td>3.05</td><td>2.46</td><td>2.46</td><td>2.46</td></tr><tr><td>Δx<sub>initial</sub> (m)</td><td>20</td><td>50</td><td>100</td><td>150</td><td>200</td><td>500</td></tr><tr><td>Computation time (s)</td><td>1.82</td><td>0.99</td><td>0.83</td><td>0.81</td><td>0.79</td><td>0.78</td></tr><tr><td>Minimum % fit<sub>DHN</sub> (%)</td><td>99.34</td><td>98.08</td><td>95.34</td><td>92.32</td><td>88.79</td><td>82.38</td></tr><tr><td>Minimum R<sup>2</sup><sub>DHN</sub></td><td>1.000</td><td>1.000</td><td>0.998</td><td>0.994</td><td>0.987</td><td>0.967</td></tr><tr><td>Average % fit<sub>DHN</sub> (%)</td><td>99.82</td><td>99.56</td><td>99.10</td><td>98.75</td><td>98.39</td><td>97.71</td></tr><tr><td>Average R<sup>2</sup><sub>DHN</sub></td><td>1.000</td><td>1.000</td><td>1.000</td><td>0.999</td><td>0.999</td><td>0.998</td></tr><tr><td>Maximum integration time step (s)</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>Maximum Co</td><td>2.46</td><td>2.46</td><td>2.46</td><td>2.46</td><td>2.46</td><td>2.46</td></tr></tbody></table>

<table><caption>**Table 4.** Computation time and fits for grey-box model simulated at varying $\Delta x$ values.</caption><thead><tr><th>Δx<sub>initial</sub> (m)</th><th>8</th><th>10</th><th>15</th><th>20</th><th>50</th><th>100</th><th>150</th><th>200</th><th>500</th></tr></thead><tbody><tr><td>Computation time (s)</td><td>43.13</td><td>26.73</td><td>11.23</td><td>9.38</td><td>6.37</td><td>2.89</td><td>2.71</td><td>2.67</td><td>2.34</td></tr><tr><td>Minimum % fit<sub>DHN</sub> (%)</td><td>95.23</td><td>94.98</td><td>95.00</td><td>95.18</td><td>92.13</td><td>85.97</td><td>85.59</td><td>79.48</td><td>40.23</td></tr><tr><td>Minimum R<sup>2</sup><sub>DHN</sub></td><td>0.998</td><td>0.998</td><td>0.998</td><td>0.998</td><td>0.994</td><td>0.980</td><td>0.979</td><td>0.958</td><td>0.643</td></tr><tr><td>Average % fit<sub>DHN</sub> (%)</td><td>98.54</td><td>98.54</td><td>98.43</td><td>98.43</td><td>97.56</td><td>95.79</td><td>95.04</td><td>92.32</td><td>75.87</td></tr><tr><td>Average R<sup>2</sup><sub>DHN</sub></td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>0.999</td><td>0.997</td><td>0.996</td><td>0.991</td><td>0.913</td></tr><tr><td>Maximum integration time step</td><td>58.18</td><td>66.93</td><td>90.76</td><td>142.78</td><td>343.38</td><td>523.66</td><td>731.72</td><td>907.81</td><td>1004.31</td></tr><tr><td>Maximum Co</td><td>6.04</td><td>8.15</td><td>6.17</td><td>7.09</td><td>12.53</td><td>32.28</td><td>59.18</td><td>51.94</td><td>77.43</td></tr></tbody></table>

The simulation solver of the grey-box model in Simulink is changed to ode15s, similar to the physical model, with a relative tolerance of $1 \times 10^{-3}$ and with the Jacobian method being set to sparse analytical.

For all simulations in Tables 3 and 4, the maximum velocity remains constant at 1.53 m/s due to the identical mass flow rate profiles. The Courant numbers are calculated using (52),

$$Co = u_f \Delta t / \Delta x_{pk} \qquad (52)$$

where $u_f$ is the fluid velocity and $\Delta t$ is the variable integration time step.

The maximum Courant numbers reported in Tables 3 and 4 are determined by taking the maximum of all Courant numbers throughout the simulation. These values may not match the product of the maximum velocity and maximum time step divided by the pipe segment length. The Courant numbers exceed 1 due to the implicit integration methods
employed by the numerical differentiation formula of the ode15s solver. Additionally, the minimum Courant number for all simulations is zero, as the integration time step at the initial condition is always zero.

From Tables 3 and 4, it can be inferred that the simulation results of the low-fidelity physical model (obtained by simulating the physical model at an initial pipe segment length greater than $\Delta x_{initial} = 0.5$ m) consistently outperform the grey-box models in terms of fit to the high-fidelity physical model. This is reflected in the relatively high % fit<sub>DHN</sub> and $R^2_{DHN}$ values regardless of the initial pipe segment length $\Delta x_{initial}$. Moreover, the low-fidelity physical model simulations are consistently faster, demonstrating higher precision compared with the grey-box model.

The low-fidelity physical model with $\Delta x_{initial} = 50$ m is recommended for future optimisation, as it achieves an average fit of 99.56% to the high-fidelity model with an $R^2_{DHN}$ value of 1. Additionally, the reduction in computation time for simulations with $\Delta x_{initial} > 50$ m is relatively small when compared with the $\Delta x_{initial} = 50$ m. By the same criteria, the grey-box model with $\Delta x_{initial} = 50$ m can also be used for future optimisation.

Although the low-fidelity physical model proves to be faster and more precise than the grey-box model, it is still essential to perform optimisation with both models to determine the best one. The reason for this is that the nature of an optimisation model influences the optimisation result, particularly with regard to local minima. In this case, the differing characteristics of the physical and grey-box models may present an opportunity to explore the advantages offered by each.

Finally, the simulation time in this paper was measured using MATLAB's tic-toc function, which may differ from the actual CPU time. However, to minimise external influences, all non-MATLAB and Simulink activities on the laptop were closed before running the simulations. Additionally, each simulation case was executed three times to ensure accurate computation time measurement.

## 4. Conclusions

This study addressed the challenge of optimising district heating networks (DHNs) by developing a computationally efficient grey-box model suitable for dynamic real-time optimisation (DRTO). Unlike traditional physical models, which are accurate but computationally expensive, and black-box models, which lack extrapolation ability, the proposed grey-box linear parameter varying (LPV) pipe model balances accuracy and efficiency. The LPV pipe model was modified to allow DHN simulations with continuous mass flow rate profiles, expressing pipe outlet temperature as a function of inlet and soil temperatures, varying model parameters (dependent on mass flow rate), and constant parameters (determined through grid-search optimisation). Validation against a high-fidelity physical model (small discretisation) for a single pipe showed a 91.5% fit with an $R^2$ of 0.993 while being 5 times faster. Further validation in a 10-consumer DHN demonstrated an average fit of 98.64% and an $R^2$ of 0.9997, with the grey-box model running approximately 52 times faster than the physical model. Low-fidelity physical models (LFPMs), developed by increasing pipe segment lengths, remained more accurate and computationally efficient than the grey-box model when validated against the high-fidelity reference. Given the structural differences between grey-box and physical models, optimisation must be performed with both to determine which identifies better local minima. This research provides a foundation for future studies on DRTO implementation, exploring trade-offs between model accuracy, computational cost, and optimisation performance.

**Author Contributions:** Conceptualisation, O.E.T., S.S., P.L., R.M., H.V. and J.-M.R.; data curation, O.E.T., P.L. and R.M.; formal analysis, O.E.T., S.S., P.L. and R.M.; investigation, O.E.T., P.L., R.M. and S.S.; methodology, O.E.T., S.S., P.L. and R.M.; software, O.E.T., P.L. and R.M.; validation, O.E.T., S.S.
P.L., R.M., H.V. and J.-M.R.; writing—review and editing, O.E.T., S.S., P.L., R.M., H.V. and J.-M.R.; project administration, S.S., P.L., R.M., H.V. and J.-M.R.; funding acquisition, S.S., P.L., R.M., H.V. and J.-M.R.; supervision, S.S., P.L., R.M. and J.-M.R.; resources, S.S., P.L., R.M., H.V. and J.-M.R.; writing—original draft, O.E.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Regional Council of Nouvelle Aquitaine, France, within the framework of the RESEAUDATA project.

**Data Availability Statement:** The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding authors.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Nomenclature

### Abbreviations

<table><tr><td>DHN</td><td>District heating network</td></tr><tr><td>DRTO</td><td>Dynamic real-time optimisation</td></tr><tr><td>IEA</td><td>International Energy Agency</td></tr><tr><td>LPV</td><td>Linear parameter varying</td></tr><tr><td>LFPM</td><td>Low-fidelity physical model</td></tr><tr><td>MISO</td><td>Multiple Input Single Output</td></tr><tr><td>MPC</td><td>Model predictive control</td></tr><tr><td>NMPC</td><td>Nonlinear model predictive control</td></tr><tr><td>RPM</td><td>Reinitialised partial moment</td></tr><tr><td>SISO</td><td>Single Input Single Output</td></tr><tr><td>SDE</td><td>Stochastic differentiation equation</td></tr><tr><td>TEG</td><td>Thermal energy generation</td></tr></table>

### Indices

<table><tr><td>pmp</td><td>Production main pipe</td></tr><tr><td>prp</td><td>Production return pipe</td></tr><tr><td>mp</td><td>Main pipe</td></tr><tr><td>ep<sub>in</sub></td><td>Extended-in pipe</td></tr><tr><td>ep<sub>out</sub></td><td>Extended-out pipe</td></tr><tr><td>rp</td><td>Return pipe</td></tr><tr><td>b</td><td>Network branch</td></tr><tr><td>nb</td><td>Number of branches in the network</td></tr><tr><td>nc<sub>b</sub></td><td>Number of consumers in branch b</td></tr><tr><td>nc</td><td>Total number of consumers in the network</td></tr><tr><td>k</td><td>Distribution pipe (pipes in the distribution section), k = mp, ep<sub>in</sub>, ep<sub>out</sub>, rp</td></tr><tr><td>p</td><td>Production pipe p = pmp, prp</td></tr><tr><td>npk</td><td>Total number of pipes in the DHN</td></tr><tr><td>pk</td><td>Production or distribution pipe, pk = pmp, prp, mp, ep<sub>in</sub>, ep<sub>out</sub>, rp, ∀pk = 1,2...npk</td></tr><tr><td>c</td><td>A consumer component (distribution pipes, splitter, mixer, and substation)</td></tr><tr><td>nd<sub>b,c,k</sub></td><td>Number of discretised segments of distribution pipe k in consumer component c in branch b</td></tr><tr><td>nd<sub>p</sub></td><td>Number of discretised segments of production pipe p</td></tr><tr><td>i</td><td>Discretisation index ∀i = 1,2...f, where f = nd<sub>b,c,k</sub> for distribution pipes and f = nd<sub>p</sub> for production pipes</td></tr><tr><td>ti</td><td>Time instance ∀ti = 1,2...nti, nti-total number of time instances</td></tr><tr><td>pm</td><td>Pipe and its protective materials, pm = pipe, casing, and insulation</td></tr></table>

### Symbols

<table><tr><td>T<sub>b,c,k,i</sub></td><td>Discretised temperatures of distribution pipe k in consumer component c in branch b (°C)</td></tr><tr><td>T<sub>p,i</sub></td><td>Discretised temperatures of production pipe p (°C)</td></tr></table>
<table>
  <tr>
    <td>T<sub>soil</sub></td>
    <td>Soil temperature (°C)</td>
  </tr>
  <tr>
    <td>T<sub>inlet</sub></td>
    <td>Temperature of the fluid at the pipe inlet (pipe inlet temperature) (°C)</td>
  </tr>
  <tr>
    <td>T<sub>outlet</sub></td>
    <td>Temperature of the fluid at the pipe outlet (pipe outlet temperature) (°C)</td>
  </tr>
  <tr>
    <td>P<sub>demand,c</sub></td>
    <td>Power demand of the substation in consumer component c (W)</td>
  </tr>
  <tr>
    <td>ṁ<sub>b,c,k</sub></td>
    <td>Mass flow rates in distribution pipe k in consumer component c in branch b (kg/s)</td>
  </tr>
  <tr>
    <td>ṁ<sub>p</sub></td>
    <td>Mass flow rate in production pipe p (kg/s)</td>
  </tr>
  <tr>
    <td>P<sub>gen</sub></td>
    <td>Power generated in the combined generation station (W)</td>
  </tr>
  <tr>
    <td>N<sub>variables</sub></td>
    <td>Total number of variables in a presumed existing DHN model</td>
  </tr>
  <tr>
    <td>N<sub>equations</sub></td>
    <td>Total number of equations formulated from a presumed existing DHN model</td>
  </tr>
  <tr>
    <td>λ<sub>pm</sub></td>
    <td>Thermal conductivity of a pipe and its protective materials (W/(mK))</td>
  </tr>
  <tr>
    <td>λ<sub>s</sub></td>
    <td>Thermal conductivity of soil (W/(mK))</td>
  </tr>
  <tr>
    <td>λ<sub>f</sub></td>
    <td>Thermal conductivity of the heat transfer fluid (water) (W/(mK))</td>
  </tr>
  <tr>
    <td>Re</td>
    <td>Reynolds number</td>
  </tr>
  <tr>
    <td>Pr</td>
    <td>Prandtl number</td>
  </tr>
  <tr>
    <td>u<sub>f</sub></td>
    <td>Fluid velocity (m/s)</td>
  </tr>
  <tr>
    <td>d<sub>b</sub></td>
    <td>Depth of buried pipes in the ground (m)</td>
  </tr>
  <tr>
    <td>e<sub>pm</sub></td>
    <td>Thickness of a pipe and its protective materials (m)</td>
  </tr>
  <tr>
    <td>L</td>
    <td>Pipe length (m)</td>
  </tr>
  <tr>
    <td>L<sub>pk</sub></td>
    <td>Length of a pipe pk (m)</td>
  </tr>
  <tr>
    <td>D<sub>0</sub></td>
    <td>Pipe internal diameter (m)</td>
  </tr>
  <tr>
    <td>D<sub>pm</sub></td>
    <td>External diameter of a pipe and its protective materials (m)</td>
  </tr>
  <tr>
    <td>S</td>
    <td>Conduction shape factor (m)</td>
  </tr>
  <tr>
    <td>ρ</td>
    <td>Fluid density (kg/m³)</td>
  </tr>
  <tr>
    <td>C<sub>p</sub></td>
    <td>Specific heat capacity of fluid (J/(kgK))</td>
  </tr>
  <tr>
    <td>μ</td>
    <td>Dynamic viscosity of fluid (kg/(ms))</td>
  </tr>
  <tr>
    <td>y</td>
    <td>Output of a general SISO system</td>
  </tr>
  <tr>
    <td>u</td>
    <td>Input of a general SISO system</td>
  </tr>
  <tr>
    <td>β</td>
    <td>Number of zeros of a general SISO system</td>
  </tr>
  <tr>
    <td>α</td>
    <td>Number of poles of a general SISO system</td>
  </tr>
  <tr>
    <td>y<sub>i</sub></td>
    <td>Real values of y at i sampling</td>
  </tr>
  <tr>
    <td>ŷ<sub>i</sub></td>
    <td>Predicted values of y at i sampling</td>
  </tr>
  <tr>
    <td>θ<sub>inlet</sub></td>
    <td>Set of model parameters for the relationship between T<sub>outlet</sub> and T<sub>inlet</sub></td>
  </tr>
  <tr>
    <td>θ<sub>soil</sub></td>
    <td>Set of model parameters for the relationship between T<sub>outlet</sub> and T<sub>soil</sub></td>
  </tr>
  <tr>
    <td>T<sub>pk</sub>(t<sub>ti</sub>)</td>
    <td>Pipe pk outlet temperature using the high-fidelity physical model at time instance ti (°C)</td>
  </tr>
  <tr>
    <td>T̂<sub>pk</sub>(t<sub>ti</sub>)</td>
    <td>Pipe pk outlet temperature using the grey-box or low-fidelity physical models at time instance ti (°C)</td>
  </tr>
  <tr>
    <td>T̄<sub>pk</sub></td>
    <td>Mean of the pipe pk outlet temperature using the high-fidelity physical model (°C)</td>
  </tr>
  <tr>
    <td>t<sub>ti</sub></td>
    <td>Time value at time instance ti (s)</td>
  </tr>
  <tr>
    <td>R<sup>2</sup><sub>pk</sub></td>
    <td>R-squared or coefficient of determination for comparing the outlet temperatures of pipe pk using grey-box or low-fidelity physical models to the high-fidelity physical model</td>
  </tr>
  <tr>
    <td>% fit<sub>pk</sub></td>
    <td>Percentage fit of pipe pk outlet temperature from the grey-box or low-fidelity physical model against the high-fidelity physical model simulation result</td>
  </tr>
  <tr>
    <td>R<sup>2</sup><sub>DHN</sub></td>
    <td>R-squared for comparing the outlet temperature of all the pipes in DHN using grey-box or low-fidelity physical models to the high-fidelity physical model</td>
  </tr>
  <tr>
    <td>% fit<sub>DHN</sub></td>
    <td>Percentage fit of all the pipe outlet temperatures in the DHN from the grey-box or low-fidelity physical models against the high-fidelity physical model simulation result.</td>
  </tr>
  <tr>
    <td>θ</td>
    <td>Model parameters to be identified</td>
  </tr>
  <tr>
    <td>T<sub>inlet_start</sub></td>
    <td>Initial value of the inlet temperature (°C)</td>
  </tr>
  <tr>
    <td>T<sub>inlet_step</sub></td>
    <td>Step increase on the inlet temperature (°C)</td>
  </tr>
  <tr>
    <td>t<sub>step_change</sub></td>
    <td>Time at which step change is imposed on the inlet temperature (s)</td>
  </tr>
  <tr>
    <td>k<sub>inlet</sub></td>
    <td>Gain of the pipe outlet temperature resulting from inlet temperature variation</td>
  </tr>
</table>
<table><tr><td>k<sub>soil</sub></td><td>Gain of the pipe outlet temperature resulting from soil temperature variation</td></tr><tr><td>τ<sub>inlet</sub></td><td>Time constant based on inlet temperature variation (s)</td></tr><tr><td>τ<sub>in</sub></td><td>Modified time constant based on inlet temperature variation (s)</td></tr><tr><td>τ<sub>soil</sub></td><td>Time constant based on soil temperature variation (s)</td></tr><tr><td>T<sub>d,inlet</sub></td><td>Time taken for the fluid to move through the pipe length (s)</td></tr><tr><td>T<sub>d</sub></td><td>Time delay of a general SISO system or the time taken for the fluid to move through a given pipe segment length (s)</td></tr><tr><td>T<sub>d,in</sub></td><td>Modified time taken for the fluid to move through a given pipe segment length (s)</td></tr><tr><td>n</td><td>Number of pipe discretisation segments during physical model development</td></tr><tr><td>n<sub>pk</sub></td><td>Number of pipe pk discretisation segments during physical model development, n<sub>pk</sub> = f</td></tr><tr><td>N</td><td>Number of pipe discretisation segments during grey-box model development</td></tr><tr><td>N<sub>pk</sub></td><td>Number of pipe pk discretisation segments during grey-box model development, N<sub>pk</sub> = f</td></tr><tr><td>T<sub>inlet</sub><br>T<sub>outlet</sub></td><td>Influence of inlet temperature on the pipe outlet temperature (°C)</td></tr><tr><td>T<sub>soil</sub><br>T<sub>outlet</sub></td><td>Influence of soil temperature on the pipe outlet temperature (°C)</td></tr><tr><td>I</td><td>Identity matrix</td></tr><tr><td>T<sub>gb</sub></td><td>Vector of the discretised temperature that shows the influence of the inlet temperature on the pipe outlet temperature (°C)</td></tr><tr><td>Δx<sub>initial</sub></td><td>Initially assigned pipe segment lengths for a DHN simulation (m)</td></tr><tr><td>Δx</td><td>Length of a discretised pipe segment (m)</td></tr><tr><td>Δx<sub>pk</sub></td><td>Discretised segment length of a pipe pk (m)</td></tr><tr><td>Co</td><td>Courant number</td></tr></table>

# Appendix A

## Appendix A.1. Thermophysical Properties of the Heat Transfer Fluid (Water)

The research study of [4] stated that the current district heating networks (DHNs), which are the 3G (third generation), operate at a supply temperature of less than 100 °C and a return temperature of less than 45 °C. Therefore, the thermophysical properties of water and the pipe materials are estimated at an average temperature of 72.5 °C. The density of water is ρ = 976.18 kg/m³, the specific heat capacity of water is C<sub>p</sub> = 4191.5 J/(kgK), the dynamic viscosity of water is μ = 0.000391 kg/(m·s) [28], and the thermal conductivity of water is k<sub>f</sub> = 0.661W/(m·K) [29].

## Appendix A.2. Thermophysical Properties of the Pipe and Its Protectives

The research work of [26] conducted linear regressions on industrial data to obtain the relationship between pipe, insulation, and casing thickness for discrete pipe diameters in the range $20 \le D_o \le 500$ mm.

Pipe thickness:

$$e_{pipe} = 0.009D_o + 0.0024 \quad (A1)$$

Insulation thickness:

$$e_{ins} = 0.1305D_o + 0.022 \quad (A2)$$

Casing thickness:

$$e_{cas} = 0.0091D_o + 0.0025 \quad (A3)$$

Also, from the research study of [26], the thermal conductivities of pipe pk in DHN and its protectives are shown in Table A1.
<table><caption>**Table A1.** Thermal conductivities of pipe material, insulation, casing, and soil.</caption><thead><tr><th></th><th>Pipe</th><th>Insulation</th><th>Casing</th><th>Soil</th></tr></thead><tbody><tr><th>Material</th><td>Stainless steel</td><td>Polyurethane</td><td>Polyethylene</td><td></td></tr><tr><th>λ (W/(m.K))</th><td>15.1</td><td>0.025</td><td>0.47</td><td>1</td></tr></tbody></table>

Overall heat transfer coefficient $U_o$:

$$ \frac{1}{U_o(t)} = \frac{1}{h_o(t)} + \frac{D_o \ln \frac{D_{pipe}}{D_o}}{2\lambda_{pipe}} + \frac{D_o \ln \frac{D_{ins}}{D_{pipe}}}{2\lambda_{ins}} + \frac{D_o \ln \frac{D_{cas}}{D_{ins}}}{2\lambda_{cas}} + \frac{\pi D_o L}{S\lambda_s} \quad (A4) $$

where $\lambda_{pipe}$ is the thermal conductivity of steel pipe; $\lambda_{ins}$ is the thermal conductivity of insulation; $\lambda_s$ is the thermal conductivity of soil; $d_b$, the depth of buried pipes in the ground, is 1 m (assumed to be uniform for all the pipes in the network); the external diameter of pipe is $D_{pipe} = D_o + 2e_{pipe}$; the external diameter of the insulation is $D_{ins} = D_{pipe} + 2e_{ins}$; and the external diameter of the casing is $D_{cas} = D_{ins} + 2e_{cas}$.

The conduction shape factor S can be derived from [30] for $L > D_o$.

$$ S = \frac{2\pi L}{\cosh^{-1}\left(\frac{2d_b}{D_o}\right)} \quad (A5) $$

The convective heat transfer coefficient $h_o(t) = \frac{k_f N_u(t)}{D_o}$.

The average Nusselt number $N_u$ can be derived from correlations based on the Reynolds's and Prandtl numbers.

The Gnielinski correlation for $3000 \le R_e(t) \le 5 \times 10^6$ and $0.5 \le P_r \le 2000$ [26] is

$$ N_u(t) = \frac{\left(\frac{f(t)}{8}\right)\left(R_e(t) - 1000\right)P_r}{1 + 12.7\left(\frac{f(t)}{8}\right)^{0.5}\left(P_r^{\frac{2}{3}} - 1\right)} \quad (A6) $$

Dittus-Boelter correlation for $R_e(t) \le 1 \times 10^4$ and $0.6 \le P_r \le 160$ [14] is

$$ N_u(t) = 0.023 R_e(t)^{\frac{4}{5}} P_r^n \quad (A7) $$

where $n = 0.3$ because the pipe outlet temperature is less than the pipe inlet temperature [30], the Reynolds number $Re = \frac{\rho u_f(t)D_o}{\mu}$, Prandtl number $Pr = \frac{C_p\mu}{\lambda_f}$, and fluid velocity $u_f(t) = \frac{\dot{m}(t)}{\rho A}$.

## References

1. IEA Report. Heat—Renewables 2023—Analysis—IEA. Available online: https://www.iea.org/reports/renewables-2023/heat (accessed on 12 November 2024).

2. Pezzini, P.; Gomis-Bellmunt, O.; Sudrià-Andreu, A. Optimization techniques to improve energy efficiency in power systems. *Renew. Sustain. Energy Rev.* **2011**, *15*, 2028–2041. [CrossRef]

3. Rafati, A.; Shaker, H.R. Predictive maintenance of district heating networks: A comprehensive review of methods and challenges. *Therm. Sci. Eng. Prog.* **2024**, *53*, 102722. [CrossRef]

4. Jiang, M.; Rindt, C.; Smeulders, D.M.J. Optimal Planning of Future District Heating Systems—A Review. *Energies* **2022**, *15*, 7160. [CrossRef]

5. Untrau, A.; Sochard, S.; Marias, F.; Reneaume, J.-M.; Le Roux, G.A.; Serra, S. Analysis and future perspectives for the application of Dynamic Real-Time Optimization to solar thermal plants: A review. *Sol. Energy* **2022**, *241*, 275–291. [CrossRef]

6. Ellis, M.; Durand, H.; Christofides, P.D. A tutorial review of economic model predictive control methods. *J. Process. Control.* **2014**, *24*, 1156–1178. [CrossRef]
7. Hotvedt, M.; Hauger, S.O.; Gjertsen, F.; Imsland, L. Dynamic real-time optimisation of a CO₂ capture facility. *IFAC-PapersOnLine* **2019**, *52*, 856–861. [CrossRef]

8. Untrau, A.; Sochard, S.; Marias, F.; Reneaume, J.-M.; Le Roux, G.A.; Serra, S. Dynamic Real-Time Optimization of a solar thermal plant during daytime. *Comput. Chem. Eng.* **2023**, *172*, 108184. [CrossRef]

9. Serasinghe, R.; Long, N.; Clark, J.D. Parameter identification methods for low-order gray box building energy models: A critical review. *Energy Build.* **2024**, *311*, 114123. [CrossRef]

10. Agner, F.; Kergus, P.; Pates, R.; Rantzer, A. Hydraulic Parameter Estimation in District Heating Networks. *IFAC-PapersOnLine* **2023**, *56*, 5438–5443. [CrossRef]

11. Wang, J.; Kong, F.; Pan, B.; Zheng, J.; Xue, P.; Sun, C.; Qi, C. Low-order gray-box modeling of heating buildings and the progressive dimension reduction identification of uncertain model parameters. *Energy* **2024**, *294*, 130812. [CrossRef]

12. Chen, Y.; Guo, M.; Chen, Z.; Chen, Z.; Ji, Y. Physical energy and data-driven models in building energy prediction: A review. *Energy Rep.* **2022**, *8*, 2656–2671. [CrossRef]

13. Bacher, P.; Madsen, H.; Perers, B. Models of the Heat Dynamics of Solar Collectors for Performance Testing. In Proceedings of the ISES Solar World Congress 2011, Kassel, Germany, 28 August–2 September 2011. Available online: www.imm.dtu.dk (accessed on 12 November 2024).

14. Nova-Rincon, A.; Sochard, S.; Serra, S.; Reneaume, J.-M. Dynamic simulation and optimal operation of district cooling networks via 2D orthogonal collocation. *Energy Convers. Manag.* **2020**, *207*, 112505. [CrossRef]

15. Li, Y.; O'Neill, Z.; Zhang, L.; Chen, J.; Im, P.; DeGraw, J. Grey-box modeling and application for building energy simulations—A critical review. *Renew. Sustain. Energy Rev.* **2021**, *146*, 111174. [CrossRef]

16. Thilker, C.A.; Bacher, P.; Bergsteinsson, H.G.; Junker, R.G.; Cali, D.; Madsen, H. Non-linear grey-box modelling for heat dynamics of buildings. *Energy Build.* **2021**, *252*, 111457. [CrossRef]

17. Svensen, J.L.; Bergsteinsson, H.G.; Madsen, H. Modeling and nonlinear predictive control of solar thermal systems in district heating. *Heliyon* **2024**, *10*, e31027. [CrossRef]

18. Blomgren, E.; D'ettorre, F.; Samuelsson, O.; Banaei, M.; Ebrahimy, R.; Rasmussen, M.; Nielsen, N.; Larsen, A.; Madsen, H. Grey-box modeling for hot-spot temperature prediction of oil-immersed transformers in power distribution networks. *Sustain. Energy Grids Networks* **2023**, *34*, 101048. [CrossRef]

19. Mercere, G.; Palsson, H.; Poinot, T. Continuous-time linear parameter-varying identification of a cross flow heat exchanger: A local approach. *IEEE Trans. Control. Syst. Technol.* **2011**, *19*, 64–76. [CrossRef]

20. Chouaba, S.E.; Chamroo, A.; Ouvrard, R.; Poinot, T. A counter flow water to oil heat exchanger: MISO quasi linear parameter varying modeling and identification. *Simul. Model. Pract. Theory* **2012**, *23*, 87–98. [CrossRef]

21. Naspolini, A.; Morato, M.M.; Flesch, R.C.C.; Normey-Rico, J.E. Enhanced gas-lift system operation using LPV nonlinear model predictive control. *Geoenergy Sci. Eng.* **2024**, *239*, 212969. [CrossRef]

22. Boutourda, F.; Ouvrard, R.; Poinot, T.; Mehdi, D.; Mesquine, F.; De Tredern, E.; Jauzein, V. A global approach to estimate continuous-time LPV models for wastewater nitrification. *IFAC-PapersOnLine* **2024**, *58*, 61–66. [CrossRef]

23. Shayan, M.E.; Petrollese, M.; Rouhani, S.H.; Mobayen, S.; Zhilenkov, A.; Su, C.L. An innovative two-stage machine learning-based adaptive robust unit commitment strategy for addressing uncertainty in renewable energy systems. *Int. J. Electr. Power Energy Syst.* **2024**, *160*, 110087. [CrossRef]

24. The Engineering Toolbox. Water—Specific Heat vs. Temperature. Available online: https://www.engineeringtoolbox.com/specific-heat-capacity-water-d_660.html (accessed on 13 January 2025).

25. Duquette, J.; Rowe, A.; Wild, P. Thermal performance of a steady state physical pipe model for simulating district heating grids with variable flow. *Appl. Energy* **2016**, *178*, 383–393. [CrossRef]

26. Elhafaia, M.H.; Nova-Rincon, A.; Sochard, S.; Serra, S.; Reneaume, J.-M. Optimization of the pipe diameters and the dynamic operation of a district heating network using a robust initialization strategy. *Appl. Therm. Eng.* **2025**, *258*, 124533. [CrossRef]

27. Tijani, O.E.; Serra, S.; Sochard, S.; Viot, H.; Henon, A.; Malti, R.; Lanusse, P.; Reneaume, J.M. Transient Modelling and Simulation for Optimal Future Management of a District Heating Network. *31ème Congrès Français Therm* **2023**. [CrossRef]

28. ThermExcel. Physical Characteristics of Water. 2003. Available online: https://www.thermexcel.com/english/tables/eau_atm.htm (accessed on 14 November 2024).

29. The Engineering Toolbox. Water—Thermal Conductivity vs. Temperature. Available online: https://www.engineeringtoolbox.com/water-liquid-gas-thermal-conductivity-temperature-pressure-d_2012.html (accessed on 14 November 2024).

30. Bergman, T.L.; Lavine, A.S.; Incropera, F.P.; Dewitt, D.P. *Fundamentals of Heat and Mass Transfer*, 7th ed.; John Wiley & Sons: Hoboken, NJ, USA, 2011.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.