Experiences from City-Scale Simulation of Thermal Grids

: Dynamic simulation of district heating and cooling networks has an increased importance in the transition towards renewable energy sources and lower temperature district heating grids, as both temporal and spatial behavior need to be considered. Even though much research and development has been performed in the ﬁeld, there are several pitfalls and challenges towards dynamic district heating and cooling simulation for everyday use. This article presents the experiences from developing and working with a city-scale simulator of a district heating grid located in Luleå, Sweden. The grid model in the case study is a physics based white-box model, while consumer models are either data-driven black-box or gray-box models. The control system and operator models replicate the manual and automatic operation of the combined heat and power plant. Using the functional mock-up interface standard, a co-simulation environment integrates all the models. Further, the validation of the simulator is discussed. Lessons learned from the project are presented along with future research directions, corresponding to identiﬁed gaps and challenges.


Introduction
Sustainable energy systems rely on a wide range of energy sources such as biomass, wind, solar energy, combustion of waste fuel, and recovered energy from industrial plants and data centers [1,2]. An integral part of the utilization of renewable energy sources is to use the available energy as efficiently as possible. Several reports [3,4] point out that an increased district heating and cooling capacity in the European Union is key to reduced energy consumption and reduced green house gas emissions, where the goal is a reduction of the energy related CO 2 emissions by 40% in 2030 compared to the 1990 levels [5]. Currently, half of the total energy consumption in the European Union is for the purpose of heating and cooling [6], and it has been concluded that there is enough excess heat in the European Union to cover the heat demands of all buildings from the service sector to households. Previous research [7,8] has shown that from an economic perspective, the competitiveness of district heating as a means of decarbonization is sensitive to the density of the cities. District heating is mainly favorable in high density urban districts, whereas other technologies such as electrical heat pumps can be favorable in rural areas or less dense cities.
With current policies, the total energy supplied by district heating is estimated to increase by around 50% by 2050. Within the European Union, in the Nordic countries, Baltic countries, and Poland, district heating serves more than 50% of the households and around 12% European Union wide. To meet the demands of an increased utilization of renewable heat sources, the Fourth-Generation District Heating (4GDH) [9] networks include lower supply temperatures with the ability to recycle heat from low temperature sources and forming a smart energy system with the integration of, e.g., electricity and gas networks. This also has implications for how district heating networks are controlled, where it was shown, e.g., in [10] that large energy savings are possible using more advanced control strategies and efficient utilization of thermal storage units and building inertia.
To explore the possibilities and implications of novel types of district heating and cooling networks, the impact of renewable energy sources, and the corresponding control strategies, dynamic simulation of the thermal grid is considered an important aid in the process. Early results in this area are from 1990s and early 2000s, with examples such as [11,12]. More recently, several European Union projects have employed district heating and cooling simulation for their use cases, such as the Optimisation of District Heating and Cooling Systems project (OPTi) [13], which is used as a case study in this paper, and the New Generation of Intelligent Efficient District Cooling Systems project (INDIGO) [14], where a simulator is used for district cooling networks. Extensive research and development are also performed within the International Building Performance Simulation Association (IBPSA) [15] aiming to create a complete open source Modelica Framework for building and community energy system design and operation. Furthermore, several articles such as [16][17][18] deal with the concept of dynamic district heating models and simulation.
Common use cases for the simulation of district heating and cooling networks are operational optimization with regard to short-term production planning such as [18], considering the production rates and scheduling of production units, including optimization of thermal storage units. In [19], a dynamic simulator was used to explore the impact of supply temperature on a combined heat and power (CHP) plant. The integration and operation, including optimization of operation, of a seasonal storage were examined in [20]. A slightly different use case is [21], where a simulator was used for the evaluation of the safety risks with regard to safety in the case of network failures. Other articles deal with the demand response for district heating using district heating simulation tools [22]. In the previously mentioned INDIGO project, the simulator was used for the design of district cooling networks.
In this paper, the simulator developed within the Horizon 2020 OPTi project is presented as a case study, including a validation example and a discussion of the challenges and limitations encountered during the project. On the basis of the case study and references, research gaps and proposed relevant research directions are identified and analyzed. The concept of a digital twin is then presented, and it is shown that the case study fulfills the properties of a digital twin and how these properties relate to the challenges presented in the paper. Finally, conclusions are given.
The paper's main contribution is to present results and lessons learned from a case study on dynamic city-scale district heating modeling and simulation. This includes a thorough discussion of modeling and simulation paradigms, automatic model generation, fidelity, computational performance, validation, and model calibration. In addition, this work complements the existing reviews of digital twins by identifying the requirements and specifics of digital twins for district heating and cooling systems.

Case Study: Dynamic City-Scale District Heating Simulation
As a practical example used in the article to illustrate the challenges of district heating and cooling simulation, the simulator developed within the OPTi project [13] is used. The project started in 2015 as a European Union Horizon 2020 project and was coordinated by Luleå University of Technology, Sweden. The project objective was to analyze and rethink the way district heating and cooling systems are architected and controlled. Within the project, OPTi-Sim, a dynamic city-scale district heating simulator, was developed to explore use cases like: • Validation of optimization based control methods; • Design of automated demand response schemes; • Exploiting the potential of passive-heat storage.
From the width of the scope, the size of the district heating grid, and the lack of measurements available, it was determined that a first principle model building on known physical relations was needed. Due to many different stakeholders in the project, this was also seen as an aid to facilitate mutual understanding. This choice had large implications with regard to the choice of methods and tools, accuracy, computational performance, and validation, where the specific challenges encountered are presented in the following sections.
To put things in perspective, the International Renewable Energy Agency (IRENA) categorizes power grid models using five different hierarchical levels [23], depending on the time scale, detail, and width of scope. The levels were originally aimed at electrical power grids, but can be easily adapted to the district heating and cooling context. A summary with some minor adaptions is found in Table 1. Table 1. District heating and cooling grid modeling levels, based on IRENA power grid modeling levels.

Hierarchical Level Time Resolution
Planning and design models Seasons-years Economic models Hours-seasons Static grid models Single point Dynamic thermal grid models Minutes-hours Dynamic hydraulic and electrical grid models Milliseconds-minutes Using the defined levels, the case study encompasses both the hydraulic and thermal dynamics and thus combines the two lowest levels. The experiences drawn from the case study relate mainly to this particular level of modeling, whereas other levels are touched on more briefly. The simulator was validated by replicating real-life scenarios, where the main objective was exploring alternative production and demand control strategies. Modelica was used as the modeling language, with industrial process components such as pumps and valves provided by commercial simulation libraries.
Three main data sources were used for the models: 1. Models based on Geographic Information System (GIS) data; 2.
Models based on piping and instrumentation (P&I) diagrams and design data; 3.
Models based on historical measurement data (data-driven models).
The piping of the thermal grid was automatically generated from GIS data, whereas production units and pumping stations were manually modeled and configured from P&I diagrams and data. Data-driven and gray box models for the heat load of consumers were generated from historical data. The GIS data of the DHC grid of Luleå span over >9000 consumers and >44,000 single pipes.
To enhance the computational performance, a topological optimization, known as aggregation, of the grid was performed based on the so-called German method described in [24]. The complete topological optimization procedure for the case study consisted of a pre-processing by pruning of erroneous data such as empty nodes or singular pipes, combined with several aggregation steps of merging nearby nodes and consumers, merging serial pipes, and merging of branches. Notably, this part of the case study consumed a large portion of the development resources.
To control the preprocessing and aggregation process, a graphical user interface (GUI) in MATLAB was developed, as seen in Figure 1. In the GUI, the thermal grid is represented as a network plot where the effect of the preprocessing is visualized. With the topological optimization steps and preprocessing described above, the full Luleå grid consisting of 44,752 single pipes and 9533 consumers was aggregated down to 3149 single pipes and 494 customers before simulation. The total number of simulated components was reduced by around 90%. After optimization, the full city-scale grid could be simulated on a Core i7 Gen 5 laptop computer at around 2-3 times real-time speed.
Examples of the automatically generated grid can be seen in Figure 2, and an example of a production unit can be seen in Figure 3.
The buildings were clustered into seven building categories by the size and type of building, validated by measurements for each building category, for simulation of the consumer heat load. Due to the large amount of buildings and lack of measurements, it was deemed unfeasible to model each building individually.
The manual and automatic control of, e.g., the power plant, pump stations, and feed forward supply temperature was replicated from available diagrams and measurements. As an example, the pump control for the pumping station consisted mainly of two parts, one for safety and the other for regulation. The safety part ensures that the pump's operational pressures do not violate the limits. The control system ensures that the pressure before the pump will not go below a certain limit and the pressure after the pump will not be higher than a threshold limit. This is achieved by the upper two PID controllers shown in Figure 4. Furthermore, the safety part will switch off the pump if the speed reaches the lowest allowed speed. The tracking of the optimal differential pressure of a selected critical point in the grid is performed by the regulation part.  The model is automatically generated from GIS data, including district heating pipes, nodes, and consumers, where the size of the consumer reflects the annual heat load.  To highlight some of the challenges encountered during the validation process, a period of 10 days in January exhibiting an outdoor temperature drop from −5 • C to −30 • C and a rebound to −15 • C was used. Large temperature changes during winter, when the space heating demand is high, can be seen as the worst case scenario with regard to model validation, exhibiting the dynamic behavior of space heating demand and possibly exciting unmodeled behavior at the production side. Notably, due to the limited computational performance, simulating more than 10 calendar days was impractical.
The validation process was based on measurements from the Luleå district heating network. The validation example is focused on the required generated power, collected from time series from the plant data acquisition system. The total power equals the power supplied to the grid from the main combined heating and power plant, plus auxiliary energy units that are sometimes needed during the winter months. The inputs to the simulation were the outdoor temperature and time, where the latter was used in the consumer simulation. In Figure 5, the outdoor temperature, the actual total power, and the simulated total power are shown. As is apparent from the figure, for certain periods of time, there is a good correspondence between measured and simulated power, whereas for other time periods, the deviation can be quite large, raising several important issues, such as:

1.
What is the needed fidelity for the use case? 2.
What is a good metric to evaluate the accuracy? 3.
How can we account for human behavior such as manual control of production units? 4.
How can we ensure that the model is up-to-date with the real process?
These questions do not have simple answers and largely depend on the use case, which means that the needed fidelity and the used metrics have to be defined in relation to the use case. This is directly related to the model validation and model update. In Section 5, the validation and update aspects are discussed, as well as which methods can be employed. For this particular case, it is hypothesized that the deviation stems mainly from the building models being too coarsely modeled, as not all buildings have measurements available for data-driven modeling or gray-box modeling or are only part of an aggregated building set. Thus, identifying the root cause of such a deviation is non-trivial, especially since there are relatively few measurements available in the grid.
The following sections conclude about the experiences from the modeling and simulation of the Luleå grid, including a discussion and proposed research directions with regard to the challenges encountered during the conducted work.

Modeling and Simulation
For any simulation project, there are several important decisions regarding modeling and simulation that need to be made. Here, the main choices are presented along with experiences, appropriate references, and proposed research directions.

Co-Simulation
The general approach to the modeling and simulation of large-scale systems can be distinguished between monolithic and co-simulation approaches [25]. In a monolithic simulation, the entire system is modeled and simulated in a single tool, whereas in the co-simulation approach's separate tools for the respective subsystems are coupled together. The advancements in co-simulation, mostly using the functional mockup interface (FMI) standard [26], have greatly increased the capability of interconnecting models from different tools.
For the case study, the grid model was compiled as a so-called functional mockup unit (FMU) according to the FMI standard and was used seamlessly in conjunction with the representation of the control system and of the operator behavior simulation in MATLAB. The approach was considered essential due to many stakeholders with different backgrounds and areas of expertise. The integration of components was mostly straightforward, but some care needed to be taken-especially with regard to physically strongly coupled components. Here, such cases were avoided by compiling strongly physically coupled components, such as the district heating grid, into a single FMU. Otherwise, methods such as transmission line modeling [27] (TLM) might need to be used.
The cost of multiple software licenses provided a large hurdle with regard to prolonged use of the simulator outside of academia. Novel business models might be needed to deal with this issue when co-simulation finds a wider adoption.

Physical or Data-Driven Models
Recently, data-driven machine learning (ML) methods such as artificial neural networks (ANNs) have gained much popularity and attention. With an application to district heating, most machine learning models have been used for forecasting of heat demands, e.g., [28]. For the simulation and prediction of complex large-scale systems such a district heating and cooling network, there are several limitations of the most popular machine learning approaches that need to be addressed to make models fit aspects such as the interpretability of results, how to train a model for specific scenarios or designs not yet realized, and how persistent excitation properties can be guaranteed for proper identification of the system. A lack of measurements in the grid, as in the case study, also limits the feasibility of a data-driven approach.
A middle ground is to use so-called gray-box models, where a data-driven model can be used to complement a physics based model. Yet another interesting option is to use the physics based model as a reference to train a data-driven model, for cases where the data-driven model has desirable properties such as the propagation of uncertainties or fast execution time [29]. Needless to say, this requires that the physics based model is in place and correct in the first place.
For the case study, a physics based grid model was used. There were few measurements available in the grid, thus data scarcity, while the physics of water flow is well studied. The model was also used as an aid to augment the understanding of the physical process. On the consumer side, a mix of physics based, gray-box, and data-driven models was used.
To summarize, the data-driven approach provides interesting methods that are often used, e.g., for heat load prediction, but to the authors' knowledge, there is no feasible way to model a whole district heating or cooling grid without explicitly using the known physics.

Acausal and Causal Modeling
Another important choice of modeling paradigms is between acausal and causal modeling. For causal modeling, represented by tools such as e.g., Simulink [30], there is always a direct causality as the models are functions with predefined outputs and inputs. For acausal modeling, however, there are also relations, such as the relation between the pressure drop and flow rate in a pipe. Equation based models are expressed in a way that is relatively easily interpretable for domain experts and engineers. Acausal and causal modeling for district heating simulation were thoroughly described and compared in [25].
For the case study, using acausal first principles methods was considered crucial by the participants, providing a common ground for discussions and the interoperability of tools. By using acausal equation based modeling, the reusability, adaptability, and extendability of models is increased compared to the causal modeling of dynamical systems.
It is the authors' belief that equation based acausal modeling will continue to play an important role in the world of district heating and cooling simulation. There are however circumstances when a causal model is a better choice, such as for the control system, being causal in its nature, or when high quality causal models are readily available for some parts of a system. Notably, combining acausal and causal models from different modeling tools is enabled by a co-simulation approach.

District Heating Simulation Tools
The Modelica model language is widely used within both industry and the scientific community for modeling of district heating and cooling networks. Other non-domain specific tools such as MATLAB/Simulink [30] can also be used for district heating simulations, but are limited to causal modeling. Using Julia [31] for simulation has recently gained popularity due to its performance focused approach, where Modia [32] provides a Modelica-like domain specific extension of Julia for the modeling and simulation of physical systems, although in an early phase of development. The tool IDA Indoor Climate and Energy (IDA ICE) [33] has also been used for district heating simulation with better scalability performance than Modelica and Simulink for the presented use case in [34], but is lacking support for FMI. A comparison of the tools can be found in [25].
There are also several domain specific tools for district heating simulation, where a comprehensive overview of available simulation tools within district heating, and their features and limitations, can be found in [35,36]. Tools such as Termis [37], PSS SINCAL [38], TRNSYS [39], and Netsim [40] are widely used for the simulation of DHC networks. Most available domain specific tools are at the time of writing either (i) static, i.e., do not provide a representation of the dynamic behavior of the DHC system, (ii) specialized for a limited amount of use cases, or (iii) do not allow for co-simulation/interfacing with other software.
For the case study, the choice of the the Modelica model language was motivated by cosimulation capabilities and well established and validated model libraries. While Modelica has been the go-to language for modeling of district heating and cooling networks for some time, the engineer and researcher of tomorrow will have a wider array of interesting choices.
It is often mentioned that fully dynamic simulations are needed for the simulation of 4GDH systems, with a larger proportion of renewable and highly fluctuating energy sources [35]. However, this and other accuracy and fidelity related questions should not be assumed, but chosen in accordance to the project requirements; static models might well be appropriate in many cases. Comparisons between static and dynamic models for different scenarios and time and spatial resolutions would provide a valuable contribution to the scientific community.

Simulation Models
Components such as pipes, pumps, and valves are well represented in commercial and open-source industrial simulation libraries. While the equations describing these components are well known and verified, it should be noted that the oversimplification of, e.g., valve or pump characteristics can have a profound effect on the accuracy of the simulation. For the case study, the accuracy of the models was not considered as a limiting factor, instead the computational performance of the simulation was experienced as the major bottleneck.
Research with regard to individual simulation models for use in district heating simulation has been mainly focused on efficient and accurate pipe models, with examples in [41,42]. The pipe is a crucial model in fully dynamic district heating grid simulations, where the long pipe lengths require a spatial distribution, and as a consequence, model improvements resulting in increased computational performance or accuracy are highly important for the overall simulation. A comparison of different pipe models and their respective accuracies compared to actual measurements can be found in [43].
For the case study, a commercial model library, providing validated models, was used. While this provides convenience, there is also a trade-off between how general or specific the model is with regard to the use case. The scope of simulating a district heating and cooling system is relatively narrow compared to that of a general purpose process industrial model library, and using more specialized models might yield better performance at the cost of increased development time.
It is the authors' experience that accuracy demands need to be carefully investigated with regard to the use cases, with model libraries chosen as a result of the imposed demands.

Simulation of Heat Consumption and Production
The simulation of heat load at the consumption and production side, including the simulation of produced heat from fluctuating sources that can depend on a multitude of different factors such as wind speed or server utilization, are important for realistic simulation results. There are many published methods for heat load prediction that can be employed for simulation such as [44,45]. Since the patterns are subject to social behavior, individual buildings can show erratic patterns that are hard to predict. Furthermore, specific consumers such as industrial plants are not always possible to predict from known data.
In the case study, consumer models both utilizing first principles and data-driven methods were used, rendering seven separate categories of consumers. There was however a significant scaling problem: while it might be possible to accurately predict the heat load of a single building, it is challenging to scale up to more than 9000 consumers. Based on this issue, the consumers were categorized into seven different types depending on the size and type of building.
A less researched area is how to accurately model production units affected by unknown control and safety structures and manual operation. How, e.g., operators choose to start peak load boilers can have a profound effect on the grid dynamics, but information about these procedures might not be available. While not as straightforward as heat load prediction, research in this direction might provide valuable insights.

Automatic Model Generation
Historically, most modeling efforts have involved manual efforts of replicating the physical reality from data and diagrams. This quickly becomes inefficient and error prone. The acquisition of relevant data can be a tedious process in itself and often requires a significant effort of preprocessing or converting the data in order to make them suitable for parameterizing the simulation models.
There exist several different formats for the data depending on the domain. Automatic model generation from GIS, as in the case study, and BIM data are covered in several articles [46][47][48], with an extension to also generate models from CAD diagrams in [16]. A promising approach, e.g., used for the forecasting of heating and cooling demands [49], is using CityGML: 3D city models. Native CityGML lacks the representation of energy relevant building parameters [49]. However, CityGML supports specific extensions, which are called Application Domain Extension (ADE) [50], and have been applied in various studies to calculate the energy demand at different scales [51].
Despite being an active research area and that the feasibility of automatic modeling has been shown, the lack of common standards make most solutions specific for the use case. Promising results within the standardization of the model specification using a metamodel concept such as SysML [52] and AutomationML [53] are proposed. In the general case, the choice of the data format is not a choice for the modeling engineer or researcher, but needs to be dealt with already at the planning and management stage.
Regardless of the data source, it can usually be assumed that not all relevant design data can be acquired. This raises an important question: How can the validity of a digital twin be ensured without complete underlying data? This process is called data enrichment and was discussed for building models, e.g., in [54,55]. It should be noted that the methods are dependent on standardized data, might require project specific adaptions, and there are no model-agnostic answers to the question.
The case study used GIS data for the automatic generation of the grid and consumers, whereas production units were modeled manually. It was experienced from the case study that the data acquisition and preprocessing of data consumed a large amount of time, and the results can seldom be reused. This is considered crucial for future usability and effectiveness with further research and industry adoption of standards, but also sharing of code and algorithms within the scientific community.

Topological Optimization
Before the thermal grid model can be generated, the topology is usually optimized for the proposed use case, with regard to the resolution of the grid, as described previously for the case study. This process, called aggregation, of pipes or consumers in the district heating grid has a profound impact on performance [24] and with an increased level also on the accuracy. The most popular methods for aggregation are the so-called Danish and German methods, respectively. A comparison of both methods and a discussion of the potentials and limitations of the aggregation methods on 4GDH systems can be found in [56]. The methods are focused mainly on static simulations, and in the Danish case, there is no support for multiple production units. It is also not straightforward to derive any objective properties with regard to the spatial or temporal resolution from the aggregation process for the presented methods.
For the case study, the algorithm was loosely based on the German method, where a significant amount of time was consumed for tweaking the algorithm for the use cases and finding a feasible trade-off between performance and accuracy.
It is suggested by the authors that aggregation methods be further researched with respect to the validity for 4GDH networks with fully dynamic simulation and multiple production units. Effective numerical methods can to some extent mitigate the need for aggregation methods from a computational performance standpoint.

Fidelity
The fidelity, such as the time scales of interest and the spatial resolution, needed for the use case has large implications on the requirements on the model and simulation related choices. Individual components are usually not specified with regard to fidelity, but rather by what features or simplifications that should be included. Relating these parameters and settings to, e.g., the time scales of interest can be challenging. With appropriate planning and requirement specification, these problems can be mitigated for many projects. However, this requires extra work and possibly different simulators for different parts of the project.
For the case study, it was apparent that the different stakeholders of the projects had different and sometimes contradictory requirements. On the substation level for example, dynamics with a time resolution of seconds can be of interest, while the spatial resolution of the grids is relatively unimportant. However, the degradation in computational performance caused by the stiffness of modeling fast dynamics can render the simulator useless for simulation scenarios ranging over long time spans.
Ideally, the same simulator could be used for different phases and stakeholders of projects. This would require that the model be automatically adapted to, e.g., the required spatial and time step resolution. For linear systems, different model reduction techniques are well developed [57], where, e.g., fast dynamics can automatically be reduced to static relationships. Model reduction is far less researched for differential algebraic systems, where the survey in [58] can serve as an introduction. To the authors' knowledge, there are no implementations for commonly used modeling tools, and it is perceived that more research in this direction would be a valuable contribution.

Computational Performance
The largest obstacle with regard to usability in the case study was the computational performance. With scenarios possibly spanning over weeks, months, or even years, several orders of magnitude of increased simulation performance would have been desirable. Debugging models with regard to computational performance is often a tedious task that requires extensive knowledge from the developer.
As an example to highlight the volatility of computational performance, a grid model with 84 consumers was generated, running a simulation of 10,000 s within Dymola using the variable step size solver LSODAR. The inputs to the model were the heat load and return temperature of the consumers, varying slowly over time. The consumer model used a simplified approach where the mass flow rate was directly calculated from the inputs.
The first example used continuous input signals generated within the simulation model. The second examples used the same signals as in the first example, but sampled with a sampling time of 1 s, to emulate how input signals are often provided to the model. The third example was similar to the second example, except that for one of the consumers, the power input signal was filtered with a first-order filter with the time constant τ = 3 s, emulating how components with faster dynamics could be introduced to the system.
The results of the simulations are provided in Table 2. For the modeling expert or someone well versed in numerical computing, this should not come as a surprise: the second example limited the maximum step size for the variable step size solver; the third example made the system stiff, causing a severe degradation in computational performance. For the engineer or new practitioner, however, this might come as unpleasant news when small changes (from an engineering perspective) cause a large degradation in performance.

Numerical Methods
The nature of the problem is that a district heating or cooling grid poses several numerical challenges. Fluid networks are non-linear, and the problem is stiff due to the difference in the time constants between the pressure and temperature dynamics, as well as the faster dynamics that might appear at the consumer or production level. The model is hard to parallelize due to the strongly coupled nature of the system. Usually, there are also components that are discrete such as the simulation of the distributed control system (DCS).
The choice of solvers is usually limited by the choice of the tool used for modeling. For the Modelica language, using sparse solvers [59] or more efficient compilers as in [60] has been proposed. Using differential algebraic equation (DAE) system solvers rather than ODE solvers has also shown large improvements for power grid models [61]. Yet another promising approach is using quantized state system solvers (QSS) [62]. The stiffness of the problem can be addressed by using multirate solvers, where a summary of the state-ofthe-art solvers and tools, including multirate techniques, for the Modelica language was published in [63]. However, as the article pointed out, while many of these concepts are promising even with current state-of-the-art tools and solvers, the performance is very limited as the model size grows to even moderately large.
Despite promising concepts, from an engineering perspective, an integrated approach is needed where all major problems causing numerical problems are addressed to a point where accuracy, robustness, and simulation speed are good enough for the employed tools and use cases. Using a co-simulation approach provides the possibility for using state-of-the-art solvers and methods from different tools-provided that they support co-simulation and that the co-simulation itself does not introduce other issues. Following up on promising concepts, it seems that the step from concept to deployment is large, likely due to the diversity and complex nature of hybrid DAE systems.
For the case study, a classic fixed step explicit fourth-order Runge-Kutta method as implemented in Dymola was eventually used. Despite many theoretical advantages of variable step size methods, in practice, it turned out that the overhead, and in some cases unpredictable behavior, of these methods was too large to provide an acceptable computational performance.
From the case study, as well as other projects conducted by the authors and the body of scientific work in this area, it was perceived that most work from previous projects could not be readily reused, so that the wheel must be reinvented again and again. A reference model that could be used for benchmarking numerical methods, aggregation methods, and so forth could be one way forward towards better knowledge transfer between projects, and in the longer run towards better performance.

Validation
A crucial step in modeling and model calibration is the validation of the resulting model. As stated in the literature, like, e.g., [64], model validation is the process of assessing the eligibility of a model and to what degree it represents the true model for a specific use. The validation data may not have been used in the modeling or the calibration and is in that sense new data for the model.
The basic setting for model validation is depicted in Figure 6, where the estimate of the output model can be used to assess the validity of the model. Using a metric on the deviation between the model or deriving a model error model [65] not only provides an understanding of the validity, but also a means to characterize the root cause for deviation and to compare different models with each other. Reflecting back to the case study, validation and model calibration were shown to be time-consuming tasks. It was often experienced that there was a discrepancy between what the engineers considered important and the accuracy metrics. The accuracy of a simulation is often quantified by metrics such as the root mean squared error (RMSE) and mean absolute percentage error (MAPE), measuring the deviation of the simulated value compared to the measured value at fixed time instances, namely y −ŷ. Reaching a sufficiently high accuracy for a complex system can in many cases be a difficult task. There exists a multitude of exogenous inputs d and design parameters that affect the dynamic behavior of a district heating network, which can only in part be acquired from data or modeled, denoted asd. Examples are detailed pipe geometries, soil temperatures around the piping, consumer and operator behavior, and the automatic control of production units.
An important question raised in [66] was whether a higher accuracy results in a more useful simulator. Higher accuracy also goes hand-in-hand with higher complexity and usually comes with a loss of computational performance or with additional parameters that require calibration, making the simulator less useful from an engineering perspective. The required accuracy is determined by the purpose or the use case, and from a scientific perspective, the law of parsimony should be applied, meaning the simplest model for the intended use is also the right one, provided it gives the same answers as a model of higher accuracy and complexity.
Accuracy metrics such as RMSE do not take into account dynamic behavior, so that moderately delayed fast dynamics might be considered by the metric as a large deviation. There exist approaches to mitigate these shortcomings such as nonlinear alignment techniques (also referred to as dynamic time warping) [67], although not widely adopted by the scientific community. Furthermore, discrete events for cases like "when A happens, B should happen within a specified time frame" also cannot be captured by just measuring the deviation between the measurement and the simulation. Using other distance measures, such as the Levenshtein distance [68], which can be used to evaluate the order by which events occur, is one option. In the end, the lack of widespread adoption is a problem in itself with using any metric beyond the most common ones. Even when there are methods that are valid for individual components, the question of how to evaluate a complex model such as a district heating grid is largely unanswered.
Another validation related issue is the assessment of the need to update the model. The decision on the update can be based on the outcome from a simulation run in relation to the system it replicates, as shown in Figure 6. There, the observations y andŷ can be used to establish a metric for the validity of the model, like, e.g., the root mean squares of y −ŷ.
While this approach seems straightforward, there are a number of factors that drive the deviation y −ŷ: • Unknown exogenous effects: A minority of the factors that affect the replicated system can be measured or observed, meaning that there is a deviation between d andd, rendering a deviation between y andŷ. While there is often a qualitative understanding of these effects, a quantitative description is usually lacking. • Modeling errors: Any model is only valid to a certain degree of accuracy and is therefore not perfectly replicating the real-life system, again rendering a deviation between y andŷ. • Wear and tear: Any component operating in a real-life situation will be affected by wear and tear over its life-cycle, leading to maintenance action and eventual component replacement. Such degradation phenomena render similar effects like modeling errors, but occur gradually. • Undocumented realizations: In line with the modeling errors are undocumented realizations, which means that the real-life system deviates in realization from its documented specification. Such effects are already present when an initial validation of a model is performed, but will nevertheless affect the monitoring performance. • Unreported system changes: A system update is clearly necessary whenever a change in the real-life system occurs, like, e.g., replacement of components. The consequence is similar to the modeling errors.
Referring back to Figure 5, it can be seen that the overall behavior in the long term was well replicated, but in several instances, there were quite large deviations. Applying a typical measure for the fit like the normalized root mean squared error indicated a poor fit value, less than 90%, as expected. Nevertheless, from an engineering perspective on designing a control strategy for the supply of heat to the grid, the accuracy was deemed sufficient by practitioners in the field. Some of the reasons for the deviations were model errors in the building models, assuming that individual buildings for a certain category had the same heat load pattern and assuming a disturbance-free scenario. It is well known that individual buildings' heat load pattern exhibits a stochastic behavior in the short term.
Consequently, a monitoring procedure needs to be resilient to these factors and render an appropriate decision on the update needs. An additional challenge in most actual district heating networks is the fact that there are few measurements available from the grid, often restricted to a handful of temperature and pressure sensors. While there is a tendency in the industry in general to include more sensors, in line with the vision of Industry 4.0 and the smart sensors concept [69], the current lack of measurements makes model calibration, validation, and update a difficult task.
Further complicating the calibration and update is the composition of the real-life system that is replicated. As seen in the principle sketch in Figure 7, the real-life system could be composed of physical system components that are well understood and for which models are available (white boxes) and physical components where partial or no knowledge of the behavior is available (gray and black boxes). Naturally, all those will interact and also be connected to software components like a control system. Already when performing model calibration, the interconnection between components plays a role in the quality of the calibration if it is primarily based on measurement data. Similarly, the interconnections need to be considered in the model update. Thus, to be able to update the models with sufficient accuracy, the measurements need to be monitored, and informative data need to be used to update the model. Data from faulty sensors, equipment under maintenance, and un-modeled modes of operation need to be discarded-a non-trivial data analysis task. Furthermore, there are few methods for updating or calibrating DAE models from known data. It can be assumed that without model updates, the model accuracy will deteriorate with time due to, e.g., equipment wear or changes or new operating procedures. Promising work in this area, using the aforementioned AutomationML and FMI, was conducted in [70].
Methods for model calibration are well established, but are often focused on methods for a limited set of specified models (e.g., system identification for autoregressiveexogenous (ARX) and state space models) and are not suitable for calibration of models on DAE form. Some promising steps were undertaken by, e.g., [71,72]  Essentially, the model update requires a number of steps where there are several methodologies available that could be used in combination with each other. Each of these methods depends user selections that affect the performance of the methods. As soon as those methods need to be used in a tool chain, their effects on each other need to be understood and quantified. In this respect and as far as the authors are aware, there are no studies that have investigated the aggregation of established methods into tool chains automating the model update.
Clearly, there is a need for research in this direction and to propose automated methodologies and tool chains that reliably monitor and update models.

The Digital Twin
The concept of a digital twin is attributed to Michael Grieves in 2003 [73] and was popularized by the inclusion in NASA's Modeling, Simulation, Information Technology & Processing Roadmap [74]. There, the following definition is provided "A Digital Twin is an integrated multiphysics, multiscale simulation of a vehicle or system that uses the best available physical models, sensor updates, fleet history, etc., to mirror the life of its corresponding flying twin. The Digital Twin is ultra-realistic and may consider one or more important and interdependent vehicle systems, including propulsion/energy storage, avionics, life support, vehicle structure, thermal management/TPS, etc." While this is a useful vision, a more applicable approach for an industrial process was described in [75], comprising the following properties of a digital twin 1.
The digital twin is the linked collection of the relevant digital artifacts including engineering data, operation data, and behavior descriptions via several simulation models. The simulation models making up the digital twin are specific for their intended use and apply the suitable fidelity for the problem to be solved.

2.
The digital twin evolves along with the real system along the whole life-cycle and integrates the currently available knowledge about it.

3.
The digital twin is not only used to describe the behavior, but also to derive solutions relevant for the real system, i.e., it provides functionalities for assisting systems to optimize operation and service. Thus, the digital twin extends the concept of model based systems engineering (MBSE) from engineering and manufacturing to the operation and service phases.
While not strictly required to fulfill this definition, enabling features to achieve the digital twin are that the model can be automatically generated from the data and that the model can be automatically updated from the data. The evolution of the digital twin can be seen as a way of addressing the weaknesses of many simulators, where the simulator provides a useful tool at certain stages, but often remains underutilized for early and late stages of the life-cycle. In [76], it was emphasized how data-driven methods used in conjunction with first principles model techniques can enhance models to adapt to changing dynamics and environmental characteristics. It is also concluded that error quantification and plug-and-play functionality in a modular approach are open research topics that are fundamental for the future of complex systems' models.
Based on the definition above, we can distinguish three levels of the digital twin: (i) the digital twin that is a snapshot of the real process at a certain stage, (ii) the digital twin that can be regenerated from underlying data when the underlying process changes, and (iii) the digital twin that is continuously updated and regenerated along with the process as it evolves. The main components of the digital twin, including all previously defined levels, are: 1.
Models of the components, including control and user behavior, 3.
Numerical methods for simulation, and 4.
Model update functionality.
Concerning the current case study, the simulator fulfilled the properties of a digital twin having the focus on the analysis and optimization of both control and operation. While the simulator can be understood as one model, there were multiple simulation models for production, distribution, and consumer-side present and fully integrated to reflect the complete system behavior. The models were derived using automated model generation, and the models reflected all components including the control system, the user behavior, and the behavior of the human operators. The implementation of the simulator using FMUs enabled the integration of models from different modeling paradigms and the simulation of different operating scenarios, optimization schemes, and climatic conditions. Such a digital twin has the purpose of the so-called what-if analysis with a high level of flexibility in terms of the analyzed scenarios.
A lacking functionality is the integrated model update and the ability of the digital twin to track the real-life system properly in parallel during real-time operation. Further, the fidelity of the models varies in relation to the scenarios that can be analyzed. To give an example, the analysis of the energy used for heating on the overall system level is well replicated by the digital twin, while the energy used for heating at a specific individual building might not.
In order to reach a full city-scale digital twin with all the above-mentioned properties, several of the identified challenges in the paper still need to be addressed.

Conclusions
In this article, experiences from the modeling and simulation of district heating and cooling grids are presented, including a case study with practical examples. Research gaps and proposed relevant research directions are presented, drawing from the experiences.
It is shown how the concept of a digital twin relates to simulation and modeling in general and that it in its essence is a way to deal with the challenges and limitations identified in the paper.
It is perceived by the authors that while each individual component will need research and development for a simulator suitable for daily use, or a complete digital twin, the integration of the concepts is an equally important issue. The adoption of standards within the industry, by tool vendors, and by academia is considered crucial.
While some of the presented research and suggested research directions are strictly from a district heating and cooling perspective, the main body is applicable to process industrial modeling and simulation in general. It is the authors hope that the growing interest in the digital twin will result in moving the industry forward towards both digital and environmental sustainability.