Dynamic simulation and control of the heat supply system of a civic building with thermal energy storage units

The integration of low-carbon energy solutions into heating and cooling systems is a growing topic of interest towards the decarbonisation of the energy sector. The complexity of interactions between thermal-hydraulic components in these systems requires moving from a steady-state analysis to a dynamic approach. Commercial dynamic process simulators provide a good platform to model and simulate the physics phenomena of these components. However, co-simulation between software engines is often necessary to incorporate process control systems and regulate key system variables. This paper presents the co-simulation of a real heating application from a civil building for health services integrating thermal energy storage (TES) units. A detailed dynamic model of the heating system under study is developed in Apros, a commercial dynamic process simulator. The process control system responsible for regulating the hydraulic pumps, valves, and heat sources is implemented in MATLAB/Simulink. The TES system is operated by the process control system based on a daily schedule, achieving a reduced utilisation of the heat generation units during peak-demand hours.


INTRODUCTION
Heating and cooling systems deliver thermal energy to a facility or a set of facilities.These are normally hydronic systems using liquid water as a heat transfer fluid (HTF) to transport and store thermal energy.Common applications include heating, ventilation and air conditioning (HVAC) and domestic hot water (DHW) systems in buildings [1].
In recent years, thermal systems have attracted interest from the international community for the decarbonisation of the energy sector.Thermal energy demand accounts for around half of the global energy consumption and over 40% of the CO 2 emissions [2].The scalability of thermal systems from smallscale facilities (e.g. a single dwelling or set of buildings) to largescale systems (e.g.city-wide networks) and their capabilities to incorporate diverse thermal energy generation technologies are crucial towards a transition to efficient, low-carbon renewable energy systems [3].
There are relevant examples of thermal systems considering different generation technologies reported in the literature.For instance, the solar-assisted HVAC system presented in [4] This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.© 2022 The Authors.IET Generation, Transmission & Distribution published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology is used to meet a significant part of the heating and cooling demands of a university building in New Mexico.A combined heat and power (CHP) generation system coupled with an absorption system to chill food cabinets in a supermarket proved to be highly efficient [5].Different solar-thermal and heat pump-based systems for net-zero energy homes were studied and compared in [6] to meet the total production of DHW with renewable energy.At a city-wide level, intelligent electric heaters may be used to integrate wind power into several thermal systems [7].District thermal systems also offer a solution to this end [8], additionally supporting the recycling of large volumes of thermal energy from a variety of sources [9].
The efficiency of thermal systems may be improved by incorporating thermal energy storage (TES) units.The effective management of thermal systems with TES devices enables the use of load shaping techniques to shift thermal energy demand in time to reduce operating costs.For instance, the implementation of a TES unit is proposed in [10] to capture the surplus solar photovoltaic (PV) generation in a factory outside working hours to improve energy utilisation.Dedicated energy storage devices in thermal networks can provide aggregated environmental and economic benefits when coupled with other energy networks [11].In [12], for example, a TES system is used to decouple the electrical and heat demand in an integrated energy system (IES) that includes a CHP generation plant; thus, allowing a flexible operation of the CHP plant to minimise generation costs.A similar work is presented in [13].The adoption of TES technologies, however, must be accompanied by appropriate control techniques and operation strategies depending on the storage technology and thermal energy sources available.Further details are provided in [14], where a review of the different control techniques for TES systems integrated into buildings is provided.
Due to their complex dynamics, the design and implementation of heating and cooling systems need to be done carefully.Computer simulations based on mathematical models provide a safe environment to perform tests to study the behaviour of specific technologies and different system architectures [15].The design, test, and verification of automatic control techniques also require mathematical models that suitably reproduce the thermal-hydraulic dynamics of the system [16,17].Similarly, the effective formulation of operational optimisation algorithms requires models that accurately characterise the dynamics of the system and that facilitate the quantification of its load shaping capacity [11].For example, an optimisation technique is presented in [18] to consider the time delays of temperature transport in a district heating network.The dynamic characteristics of a similar network are studied in [19].In [20], a full dynamic model is used to characterise the operational flexibility of buildings provided with TES units.Simplified models of thermal systems considering steady-state hydraulic and thermal energy balance are available in the literature [21,22].However, these models only provide information on the final state of the system under study for a specific operating point, disregarding its transient behaviour.Conversely, dynamic models based on differential equations accurately reproduce the dynamics of a thermal system but at a higher computational cost.
Several commercial software packages for modelling of dynamic systems, also called dynamic process simulators, provide accurate predefined one-dimensional dynamic representations for the main elements of thermal-hydraulic systems [23,24].Models available in these simulators, however, are generally focused on the needs of specific industries.For instance, the multi-domain software ESP-r was originally intended for the modelling and simulation of HVAC systems, thermal dynamics of buildings and inter-zone airflows [25].Thus, the predefined models in a dynamic process simulator may limit the software capabilities to represent a thermal system with the desired fidelity for holistic system-level studies [26,27].In addition, most of these software packages lack a practical framework to design and incorporate advanced process controllers into the simulation environment and, without these, understanding the limitations and potential of a control system may not be possible.A control engineer developing a new control algorithm may use specific software engines to this end, such as MATLAB [28,29]-widely used to model, simulate and analyse control systems.
To circumvent the limitations of existing dynamic process simulators, co-simulations are commonly conducted between a software engine with different modelling domains and suitable software for implementing process controls.Relevant examples can be found in the literature.In [26], a co-simulation between ESP-r and TRNSYS is carried out to model the heating system of a building powered with solar energy, with the building's geometry and the thermal characteristics of the building envelope being considered.A co-simulation architecture is proposed in [27] for simulating the thermochemical and electrical domains of a nuclear power plant through a master program in MAT-LAB.A co-simulation between nxtStudio and Apros (a commercial dynamic process simulator) is adopted in [28] to incorporate an event-based control system to the model of a district heating network.A similar approach is used in [29], where a controller developed in NAPCON is incorporated to the model of a refinery furnace in Apros.
Although the co-simulation methodologies presented in [26][27][28][29] are effective to simulate or incorporate process control systems into thermal networks, fictional case studies or simplified models of the thermal systems have been adopted in the references.For instance, part of a real district heating network was modelled in [28] considering an ideal heat source with constant pressures and temperatures in the supply and return lines, but not accounting for the pressure and supply temperature control systems in the network.Due to confidentiality restrictions, models were not reported in [27].In all cases, interactions with heat dissipating stations were considered either as ideal or modelled based on parameters such as the outdoor temperature.The use of historical data to assess the performance of the thermal networks under real operating conditions was not considered.
MATLAB/Simulink has been used to conduct cosimulations alongside Apros.In [30], the modelling and control of a counter-flow heat exchanger is performed in MATLAB and then corroborated through a co-simulation with a similar component available in Apros.Similarly, the design and verification of a state-of-charge estimator of a TES unit was reported in [31], where a high-fidelity water tank component available in Apros is adopted as the TES system but the estimator and feedback control system are implemented in MATLAB/Simulink.Although the co-simulation approach adopted in [30]a n d [31] is useful for control system design and to analyse the performance of individual thermal-hydraulic components via the functionalities available in MATLAB/Simulink, real data is not employed to verify the mathematical models and, more importantly, the focus is on specific components rather than on the operation of the heating network.
This paper presents a co-simulation framework between Apros and MATLAB/Simulink to incorporate a process control system into the detailed dynamic simulation of a real complex heating network from a public facility dedicated to health services.The thermal network under investigation consists of multiple thermal-hydraulic components and actuators.The dynamic model of the network was developed in Apros using available library models.The process control system includes several loops and logic operational orders and was implemented in MATLAB/Simulink.Heat supply and demand profiles based on historical data were included for one week of simulation time.These heat flows have their own dynamics, which in turn affect the operation of the thermal network.
As a case study, the co-simulation framework was applied to analyse the effectiveness of incorporating TES units operated by the process control system into the facility's heating system to provide it with load shifting capabilities.It is shown that a timeschedule control approach based on the heat load profiles of the building allows the TES units to shift the heat demand from low-demand to high-demand periods.Consequently, a reduced utilisation of the heat generation units is achieved during peakdemand hours.

MODELLING OF THERMAL-HYDRAULIC SYSTEMS
The dynamic models of the components in a thermal-hydraulic system are based on differential equations describing the HTF transport phenomenon (i.e.hydraulics) and the propagation of temperature (i.e.thermal).Figure 1 shows the most common elements in a heating or a cooling system: pipelines, valves, hydraulic pumps, heat exchangers, hydraulic separators and tanks.Dynamic process simulators provide one-dimensional models for each of these elements.

Modelling approach for hydraulic components
The dynamic model of a hydraulic component is formulated in a general form by the Reynolds transport theorem expressed for linear momentum and assuming one-dimensional flow and homogeneous-incompressible fluid conditions [32]: where  In (2), Δp f is given by where D h [m] is the hydraulic diameter of the flow path within the component and L [m] its internal length.The dimensionless friction factor f approximates the friction effect produced by the flow velocity of the fluid and the component's geometrical and physical characteristics.According to the flow regime (dictated by the value of the Reynolds number), f is approximated by interpolating the value between the maximum friction factor for laminar flow, and the Blasius and Colebrook equations for turbulent flow [33].For all hydraulic components except for hydraulic valves, Δp L is expressed as where K L is the loss coefficient-a dimensionless number inherent to the hydraulic component.Although K L depends on the flow regime, it is approximated as a constant value in practical applications [32].For hydraulic valves, Δp L is expressed by where the flow coefficient ] is dependent on the opening area of the valve A op [m 2 ] and its curve is experimentally characterised.Both linear and quadratic approximations of K v are used in dynamic process simulators.
Pressure boost Δp B is only included in models of hydraulic pumps and is given by where H x [m] is the pressure head produced by the mechanical input of the pump, which depends on the pump's rotational speed  x [rad/s] and volumetric flow q x [m 3 /s].The relation between these three variables is often characterised by generating two-dimensional look-up tables from experimental data at each operating condition of the pump [17].Alternatively, dynamic process simulators provide quadratic approximations based on the maximum and nominal values of the pumps' head (H max [m] and H nom [m]), rotational speed  nom [rad/s] and volumetric flow q nom [m 3 /s] as The models of heat exchangers and hydraulic separators are conceptually different from the previously mentioned as they are the interfaces connecting two or more hydraulic circuits.Heat exchangers consist of a primary and a secondary hydraulic channel, and each channel is modelled as an individual hydraulic component [32].Hydraulic separators, on the other hand, are modelled as single hydraulic components with an equal number of inlet and outlet ports.

Modelling approach for thermal components
The dynamic model of a one-dimensional thermal component with a volume V [m 3 ] and containing an incompressible singlephase HTF is given by [34] where where is the surface of the contact layer through which heat is transferred to and from the surroundings by convection or conduction, h [W/(m 2• C)] and k [W/(m • C)] are the convection and conduction heat transfer coefficients, and Δx s [m] is the width of the contact layer through which heat is transferred to and from the surroundings by conduction.Both h and k depend on the medium's temperature, and h is also dependent on the flow regime [34].
As when describing their hydraulic dynamics, heat exchangers are modelled as two separate thermal components exchanging thermal energy mainly by convection [30,35].Similarly, hydraulic separators are modelled as single thermal components with an equal number of inlet and outlet ports.

Modelling of thermal energy sources and demand
The dynamics of thermal generation plants are complex and may be modelled in dynamic process simulators by integrating several dynamic components [36][37][38].However, these plants are often modelled as simple heat or temperature sources interacting with the thermal-hydraulic network [17,28,39,40].
The thermal energy demand of a heating or cooling network consists of two components: the space heating or space cooling demand and the DHW demand.The space heating or space cooling demand can be approximated by modelling a building's thermal dynamics as a set of thermal capacitances and resistances interacting with the ambient temperature [41].DWH demand, on the other hand, cannot be easily modelled physically.Instead, DHW consumption profiles are often generated based on technical standards, statistical or stochastic models, or directly implemented from empirical data [42].In this paper, thermal energy demand is assumed as a single external heat flow consuming thermal energy from the network according to historical heat consumption data (as done in [39,40]).

Implementation of dynamic models
The heating system under study (see Section 4 for further details) was implemented and simulated in Apros, a dynamic process simulator based on the 'simantics constraint language' developed by VTT and Fortum.Apros is commonly used for simulation of power plants and thermal-hydraulic systems [27-29, 33, 36].The thermal-hydraulic library provides predefined one-dimensional models of the components discussed in Section 2.
The hydraulic behaviour of pipelines and water tanks integrated into the thermal-hydraulic model is described by Equations (2), (3)a n d( 4).Pipeline elements were dimensioned according to the length and diameter specifications of each section of the heating system.Water tanks were dimensioned based on TES units commercially available.Due to the lack of FIGURE 2 Implementation of heat generation and heat consumption units information on the pressure loss for each system section, a constant loss coefficient K L = 0.1 was adopted for all components.The thermal behaviour of pipelines and water tanks described by (8) was limited to the effect of mass exchange with the subsequently connected components.Discretised long pipelines and water tanks components were used with an average of 0.1 m 3 per calculation node.This enabled to split the components into several water volumes with a local temperature calculation [23,43].
The hydraulic behaviour of valves as described by ( 2), ( 3) and ( 5) was implemented considering a linear approximation of the flow coefficient K v at the maximum valve opening provided by the manufacturer.For pump components described by ( 2), ( 3), ( 4)and( 6), the hydraulic dynamics were characterised using experimental curves found in manufacturer's datasheets relating the pump's head (H x ) and volumetric flow (q x ) at its nominal speed ( nom ) and interpolated with (7) for different values of the pump's speed ( x ).The thermal behaviour of valves and pumps was limited to the effect of mass exchange with the subsequent connected components.
A relevant aspect to consider when implementing thermalhydraulic systems is heat loss.All the components in a heating system are covered with an insulating envelope to minimise heat exchange with the environment.However, significant heat losses are present in long pipeline sections, generating a drop in the outlet temperature with respect to the inlet temperature.The heat loss of long pipelines was modelled using (10)a sa conduction heat transfer between a pipeline's internal and external walls.According to the specifications of the heating system, a 40 mm thick insulating layer with a constant conduction heat transfer coefficient value k = 0.025 W/(m • C) was considered.Due to lack of information on the environmental temperature, a constant outdoor temperature of 10 • Cw a s considered.
Heat generation and consumption units responsible for heat supply and demand are found in secondary thermal-hydraulic loops.As shown in Figure 1, these are coupled to the main heating network through hydraulic separators or heat exchangers.For simplicity, these thermal-hydraulic couplers and secondary hydraulic loops were omitted in this work.Instead, heat supply and demand were implemented as single external heat flows generating or consuming heat from a pipeline element, as illustrated in Figure 2. Thus, the thermal model of these pipeline components is given by where Qext [W] is the external heat flow with positive values for heat supply and negative values for heat consumption.

Implementation of the process control system and external heat flows
The operation of the heating system model is regulated by a process control system consisting of several logical operational orders and feedback loops.This control system manipulates the rotational speed of pumps, the opening area of valves and the heat flow supplied by the heat generation units to regulate pressure, mass flow rate and temperature in different sections of the system.The process control system was implemented and simulated in MATLAB/Simulink, linking it to the system model in Apros through the open platform communications (OPC) client/server interface (see Figure 3).The Simulink OPC toolbox was used for reading and writing the controlled and manipulated variables in real-time from and to the model in Apros.
The external heat flows that were not manipulated by the process control system were adjusted using historical data from the real heating system (see Section 4).Operational data of the system was collected for one week and supplied to the system model through MATLAB/Simulink.These external heat flows considered the demand from consumption units and the heat supplied by uncontrolled generation units.

System description
The system under study is the heating distribution network from the Queen Elizabeth Hospital King's Lynn, a civic building dedicated to health services in the Norfolk County in England.The facility is a two-storey building within an approximate area of 25,000 m 2 , hosting a staff of around 2,800 people plus daily visitors.The facility is segmented into zones that consume heat from a primary heating network, namely heating zones (HZs), and distribute it locally through secondary pipeline networks.Each of these zones is considered as a single heat consumption unit in the primary heating system.The architecture of the primary heating system consists of 1200 m of pipelines integrating six HZs and five heat generation units through the main and the boiler pipeline loops.The main loop connects the six HZs to the boiler loop with supply and return pipelines, as shown in Figure 4.A motorised control valve located at the inlet port of each HZ adjusts the water flow fed to the zone-and thus the heat flow.The boiler loop integrates the waste heat production of two electricity-driven CHP units and four gas boilers located in HZs 1, 2, 4 and 5 through a single pipeline loop, as illustrated in Figure 5.Each gas boiler has a capacity for supplying up to 2 MW of heat to the network.A pair of hydraulic pumps is located at the end of the boiler loop to circulate water through both hydraulic loops.
The TES system designed for the case study consists of four hot water tanks (HWTs) connected in series, as shown in Figure 6, with a total capacity of 100,000 L. A set of motorised control valves alternates between the operations of charging and discharging of the TES system.A hydraulic pump circulates water through the HWTs for charging.For discharging, water is circulated by the pumps in the main loop.The bottom of HWT FIGURE 6 TES system diagram illustrating temperature and mass flow rate in each pipeline section 4 is connected to the return pipeline of the main loop, as shown in Figure 4, close to the outlet port of HZ 1.The top of HWT 1 is connected to the boiler loop (see Figure 5) after the gas boiler in HZ 1 for charging.For discharging, the TES system is connected to the boiler loop before the gas boiler in HZ 1.This arrangement is intended to avoid disturbing the supply temperature in the system when HWTs are discharging water at a lower temperature.

Heating system control and dispatch scheme
The process control of the heating system under study consists of four control elements: the marginal differential pressure control, the supply temperature control, the return temperature control, and the TES system control.Details on the design and implementation of the different automatic control systems for each one of these control elements can be found in the Appendix.
The marginal differential pressure control maintains a constant pressure differential at the periphery of the heating network.Thereby, enough pressure boost is generated to circulate water through the whole pipe network.A feedback control loop with a proportional-integral (PI) controller was implemented on the circulation pumps in the boiler loop, adjusting the pumps' velocity to maintain a constant pressure differential p marg = 0.5 bars between the inlet and outlet ports of HZ 4.
The supply temperature control is based on a decentralised dispatch scheme of the different heat generation units connected to the boiler loop.It is used to maintain a constant temperature T sup = 80 • C in the supply pipeline of the main loop.The CHP units constantly supply the network with a variableuncontrolled heat flow, raising the water temperature in the boiler loop from a T ret [ • C] to a T CHP [ • C] (see Figure 5).However, the gas boilers are fully controlled and activate when the temperature at the inlet port of a boiler drops below its setpoint.This only happens when the heat demand is greater than the heat supplied by the CHP units.
To properly regulate T sup , a feedback control loop with a PI controller was implemented for each gas boiler.An additional feedforward loop with information on each boiler's mass flow rate and inlet temperature was included to improve the tracking accuracy.The temperature setpoints used for the gas boilers in

TES system operation and control
To control the operation of the TES system in an effective way, minimum and maximum operating temperatures need to be defined.These temperature constraints help to determine the minimum and maximum state-of-charge (SoC) of the water tanks.To accurately calculate the SoC of TES system, several temperature sensors in the storage vessels are normally required [44].However, for simplicity, the operation of the TES system investigated in this paper is limited to detecting a minimum temperature at the top of HWT 1 and a maximum temperature at the bottom of HWT 4 to determine if the TES units are either fully discharged or fully charged.The minimum and maximum operating temperatures established are T TES ,min = 75 • Ca n d T TES ,max = 83 • C. Thus, the storage capacity of the TES system is around 800 kWh.

Charging operation
Charging is carried out based on a pre-programmed daily schedule.An average heat demand value is calculated from the daily heat demand profile of the HZs to identify the high and lowheat demand periods (see Figure 7).When heat demand is expected to be lower than the daily heat demand average for more than a half-hour, a charging order is commanded by the process control system.The HWTs are charged at a constant mass flow rate ṁcharg = 20 kg/s-ensured by implementing a feedback control loop with a PI controller on the hydraulic pump of the TES system.
Charging is also performed at a constant temperature.To fully charge the TES units, the supply temperature setpoints of the gas boilers in HZs 1 and 2 are temporarily raised to T GB1 = T GB2 = 83 • C. The charging time of the TES system is limited either by the daily charging schedule (i.e. the time during the charging order is active) or by the maximum capacity of the TES system (i.e. when the bottom of HWT 4 reaches the maximum operating temperature).

Discharging operation
The discharging operation aims to reduce the heat supply of the heating system during the daily peak-demand hours.To this end, the discharging of HWTs is triggered based on the heat peak-demand calculated from the daily heat demand profile.Whenever the total heat demand of the heat generation units (i.e. the heat supplied by the CHP units and gas boilers) exceeds 85% of the heat peak-demand value, a discharging order is commanded by the process control system.Table 1 shows the daily heat peak-demand of HZs.HWTs are discharged at a constant heat supply of 600 kW.Thus, the fully charged TES system provides approximately 1 h and 15 min of continuous heat supply.To this end, a feedback control loop with a PI controller is implemented on the discharging motorised valves, varying the opening of the valves to regulate the discharging mass flow rate ṁdisch .The setpoint of the control loop is calculated based on discharging inlet and outlet port temperatures T disch,in and T disch,out as

RESULTS AND DISCUSSION
The heating system was simulated for one week with and without the TES system in place.Figure 8 shows the behaviour of the marginal differential pressure Δp marg and the supply and return temperatures of the boiler loop, T sup and T ret .Δp marg and T sup were kept constant by the marginal differential pressure and the supply temperature controls.However, T sup reduced several times over the weekend (i.e.Friday, Saturday and Sunday)coincidently with the highest heat peak-demands of the week (see Figure 7 and Table 1).The decrease in T sup is reflected in T ret , which excessively drops at the same time.This transient behaviour is also observed in the return temperatures of the HZs, shown in Figure 9, where the return temperature control in HZ 2 failed to maintain its temperature setpoint.To prevent these decrements in the system's temperatures, the mass flow crossing through HZ 2 can be increased by modifying the setpoint of the marginal differential pressure control.The number of reductions in T sup is decreased for the heating system when the TES system is in operation (see Figure 8c).This shows that the TES system supports the heat generation units during heat peak-demand periods.Figure 10 shows the heat requested to the gas boilers during the week.A shift in the heat demand pattern of the gas boilers in HZs 1 and 2 is observed due to the presence of the TES system, increasing the heat demand during low-demand periods but decreasing it for peak-demand periods.This is especially noticed in the utilisation of the gas boiler in HZ 4, with the number of turns-on being reduced during heat peak-demand periods.
Table 2 shows the daily heat peak-demand observed by the heat generation units and the reduction of this demand resulting from incorporating the TES system.A significant reduction in the daily heat peak-demand is achieved on Monday, Tuesday, Wednesday, and Sunday.However, this is not significantly reduced from Thursday to Saturday due to the number of the daily peak-demand periods.Although a reduction of the total heat demand in the system is observed on Thursday and Saturday (see Figure 11a), this is not observed on Friday because the TES system remained fully discharged during the peak-demand periods (see Figure 11b).
The case study was simulated considering one-week heat supply and demand profiles due to the limited availability of histor-ical data.Weekly heat supply and demand patterns often capture the impact of periodic social activity in the thermal systems of buildings [45].However, these patterns may change seasonally due to variations in the outdoor temperature.Should data availability not be a constraint, further simulations with larger time scales could be conducted to understand the impact of the seasonal variability in heat supply and demand patterns and, as a result, on the system's performance.This may lead to observe a reduced utilisation of the gas boilers and hydraulic pumps during summer.Conversely, reduced supply and return temperatures may be observed in the thermal network during winter if the substation controls fail to maintain their setpoints, as occurred in the simulation results presented.The TES system may not be significantly affected by seasonal variation in the heat supply and demand patterns since its operation is based on the periodicity of the heat demand.However, disturbances to the periodic activity of the facility, such as holidays, may influence the operation.

Scalability
The co-simulation approach adopted in this paper to model, implement and simulate real thermal networks is scalable.It enables conducting simulation studies to analyse the dynamic performance of larger and more complex network architectures.Also, the synergies between thermal networks with other energy systems may be modelled and analysed.However, this may require considering the time scales of different energy vectors within an IES system.For example, electrical systems present dynamics in a time frame of milliseconds to seconds [46], while for thermal systems they range from minutes to hours [34].Thus, a smaller time step will be necessary to accurately solve the electrical part of the model of a thermal-electrical energy network if this is of interest.Co-simulations between different software engines where electrical and thermo-mechanical dynamics of an IES are considered have been reported in [27].The system is solved using a two-time scale method.A similar approach is reported in [47] to solve the model of a hybrid natural gas and electricity system in MATLAB/Simulink.This example is representative as the dynamics of the gas system are significantly slower than the electrical ones.However, integrating the electricity network to the thermal network of the system under study falls out of the scope of this work.

Limitations
An important limitation when modelling large and complex thermal systems is the amount of technical information required to parametrise all their dynamic components (e.g.size, heat transfer coefficients, pressure loss coefficients).Although datasheets from manufacturers often provide relevant component parameters, this information is usually ideal-and for old installations datasheets may no longer be available.Thus, obtaining this data represents a challenging task in the modelling process.
Although complex systems with many thermal-hydraulic components can be modelled in dynamic process simulators such as Apros, longer simulation times may be observed due to the increased number of differential equations to be solved.When conducting simulations for an extended time scale (e.g.months or a full year), this problem may be circumvented by reducing the accuracy of the dynamic calculations.Alternatively, high-performance computers may be adopted to maintain accuracy.
As highlighted in Section 5, the simulation time scales may be restricted by the availability of historical data.To investigate the effect of seasonal variability, a larger dataset would be required.If this is not possible, representative datasets available in the literature could be adopted or new ones artificially created-at the expense of decreasing the practical value of the simulation results.

CONCLUSIONS
The design, test and verification of heating and cooling applications require studying the dynamic interactions between their components in a safe environment.Dynamic process simulators serve as a platform for modelling and simulating the thermal-hydraulic behaviour of each element, reducing time and capital investments for developing systems.An adequate co-simulation framework between software engines is often necessary to include advanced process control systems and regulate relevant variables of the dynamic model.In this line, a co-simulation framework between MATLAB/Simulink and Apros was adopted to implement the process control system and external heat flows in thermal-hydraulic systems.
The co-simulation approach was demonstrated by analysing the implementation of a TES system in a conventional heating network from a public facility.A dynamic model of the heating system was developed in Apros using predefined library models of pipelines, water tanks, valves and pumps.Conversely, the process control system, responsible for regulating the operation of pumps, valves and heat generation units, was implemented in MATLAB/ Simulink.Communication with Apros was establishing using an OPC server.
The heating system was simulated for a week, from Monday to Sunday.The operational control of the TES system, based on a daily schedule, proved to help reduce the utilisation of the heat generation units during daily heat peak-demand periods.Through the presented case study, it has been shown that the co-simulation framework is useful to test and verify the dynamic performance of a heating system for different operating conditions and configurations.For instance, different charging schedules, charging mass flow rate and discharging heat supply values of the TES system may be examined.Similarly, the effect of the feedback control loops on the overall system performance can be studied and modified if required.
The simulation results unveiled the possibility for enhancing the performance of the heating system under study.As it has been observed, the dynamics of HZ 3 affect the overall performance of the system-reducing the return and consequently the supply temperature in the thermal network.Although the implementation of the TES units helps to reduce the number of times temperatures drop in the system, extra measures can be taken from the process control system to prevent this behaviour.

APPENDIX Control of the marginal differential pressure
The control of the marginal differential pressure described in Section 4.2 considers a single feedback control loop on the circulation pumps of the heating system, as shown in Figure A1.A PI control structure is adopted.The PI controller, with proportional and integral gains k P,m and k I ,m , modifies the velocity of the circulation pumps  The non-linear hydraulic dynamics of the circulation pumps and of the heating system make the analytical design of the PI control system a challenging task which falls beyond the scope of this paper.Therefore, k P,m and k I ,m were obtained heuristically, with a response settling time between two and three minutes being considered.

Control of the supply temperature
The supply temperature control of the heating system consists of four distributed control loops implemented on each one of the gas boilers.The gas boilers are modelled as a single heat flow QGB [W] connected to the heating system through a hydraulic separator, as illustrated in Figure A2 where C p,H [J/(kg • C)] and  H [kg/m 3 ] are the specific heat coefficient and the density of the water at T H . From a control perspective, T H and QGB are the control and manipulated variables of the dynamic process described above.However, since ṁH and T H ,ret are independent system variables, the first term in (13) (i.e.ṁH C p,H (T H ,ret − T H )) is seen as a system disturbance.Thus, the effects of the manipulated variable and the disturbance on the control variable are modelled separately as where g H (s) is the transfer function of the dynamic process, g d (s) is the transfer function of the disturbance and ΔT H is the To effectively control T H , a feedback control loop with a PI controller was implemented.An additional feedforward loop was included to reduce the effect of the disturbances in the system, as shown in Figure A3.Here, the feedforward controller k ff is given by The proportional and integral gains of the PI controller implemented on each gas boiler, k P,GB and k I ,GB , were defined using classic control techniques [48] for a settling time of ten minutes.
It should be noted that the approach adopted in this paper to model and control the supply temperature in each gas boiler considers an ideal heat generation (i.e.QGB ).Thus, the effects related to the boilers' dynamics and controls are neglected.

Control of the return temperature
The return temperature control consists of six distributed control loops implemented on the HZs of the hospital's main building.These control loops prevent an excessive drop in the return temperature.The heat demand in each HZ is modelled as a single heat flow Qd  where C p,H [J/(kg • C)] and  H [kg/m 3 ] are the specific heat coefficient and the density of the water at T H . Similar to the gas boiler control scheme, the control variable of the dynamic process is T H .However, for this control scheme, the manipulated variable is ṁH and the disturbance is given by Qd .Thus, the effects of the manipulated variable and the disturbance on the control variable are modelled separately as g * H (s) is the transfer function of the dynamic process, ΔT H is the differential temperature across the gas boiler's hydraulic separator (Δ T H = T H ,sup − T H ), and g * d (s)i st h e transfer function of the disturbance.
Although ṁH could be firstly seen as the manipulated variable, it depends on the motorised control valve opening area O v [m 2 ] and pressure conditions (see Section 2.1).Due to the highly nonlinear behaviour of the motorised control valve and its effect on T H , a conventional PI control approach was used [49].A feedback control loop with a PI controller was implemented, as shown in Figure A5, manipulating O v to directly control T H .The proportional and integral values of the PI controller, k P,HZ FIGURE A6 Control schemes for charging and discharging the TES system and k I ,HZ , were obtained heuristically, where a response settling time of approximately ten minutes was considered.

Control of the TES system
The control of charging and discharging operations of the TES system are performed by two different feedback control loops, shown in Figure A6.
In the charging operation, hot water is circulated into the HWTs of the TES system by a local hydraulic pump.A PI controller with proportional and integral gains k P,ch and k I ,ch manipulates the rotational speed of the TES system's pump  x,TES [rad/s]; thus, providing sufficient pressure boost ΔP B,TES [bar] to generate a charging mass flow rate ṁcharg [kg/s] in the HWTs.ṁcharg was established as the control variable and a feedback loop was implemented to maintain a mass flow rate setpoint ṁre f ,ch = 20 kg/s.
In the discharging operation, a mass flow rate ṁdisch [kg/s] is generated within the HWTs by opening a pair of motorised control valves.An opening area of the valves O v,TES [m 2 ] produces a pressure drop ΔP v,TES [bar] across the TES system due to operating pressures in the heating network.Thus, a PI controller with proportional and integral gains k P,disch and k I ,disch manipulates O v,TES to maintain a mass flow rate setpoint ṁre f ,disch [kg/s].In turn, ṁre f ,disch is determined from the differential temperature in the TES system ΔT TES [ • C], given by (12), to maintain a constant heat flow QTES ,disch = 600 kW.
The proportional and integral gains of the charging and discharging PI controllers (i.e.k P,ch , k I ,ch , k P,disch and k I ,disch )were determined heuristically, considering a response settling time between two and three minutes.

FIGURE 1
FIGURE 1 Common thermal-hydraulic elements in heating or cooling systems [kg/m 3 ] is the fluid's density, d⃗ v∕dt [m/s 2 ] the rate of change of the fluid's linear velocity, and ∑ ⃗ F [N] the summation of all the external forces acting on the volume V [m 3 ]of the hydraulic component.∑ ⃗ F is expressed by the product of the flow area A f [m 2 ] and the pressure differential Δp C [Pa] across the component as ∑ ⃗ F = A f Δp C .For each hydraulic component, Δp C is given by ) where the first two terms on the right-hand side (i.e.1∕2(⃗ v 2 0 −⃗ v 2 ) +g(z 0 − z )) represent the pressure loss due to the hydrostatic and dynamic pressure differentials across the component, g [m/s 2 ] is the gravitational constant, z 0 and z [m] are the heights of the inlet and outlet ports of the component, and ⃗ v 0 [m/s] the velocity of the fluid in the inlet port of the component.Terms Δp f , Δp L and Δp B [Pa] are, respectively, the pressure loss due to the fluid's friction against the component's internal walls, the pressure loss due to sudden changes in the flow path within the component (i.e. in the flow area, fittings and turns) and the pressure boost generated by the component.
ḢC [W] is the rate of change in the enthalpy of the fluid contained in the component, ṪC [ • C/s] is the rate of change in the temperature of the fluid, and  [kg/m 3 ]andC p,C [J/(kg • C)] are the density and specific heat value of the fluid at a temperature T C [ • C].The first term on the right-hand side of (8) (i.e.ṁ(C p,in T in − C p,C T C )) is the rate of change in the component's enthalpy owing to the mass exchange with other thermal components.The rate of mass flow exchange is given by ṁ [kg/s], where T in [ • C] is the temperature of the fluid entering the component and C p,in [J/(kg • C)] the specific heat value of the fluid at T in .∑Qconv and ∑ Qcond [W] are the summation of the different heat transfer events between the component and its surroundings (i.e. with the environment or other thermal components) by convection and conduction.These thermal dynamics are individually given by

FIGURE 3
FIGURE 3 Real-time data exchange between MATLAB/Simulink and Apros through the OPC server

FIGURE 4
FIGURE 4 Main loop diagram with specification for the temperature level in each pipeline section

FIGURE 5
FIGURE 5 Boiler loop diagram illustrating temperature and mass flow rate in each pipeline section

FIGURE 7
FIGURE 7Weekly heat demand profile, daily heat demand average and TES system charging schedule

FIGURE 8
FIGURE 8 One-week simulation performance of: (a) Marginal differential pressure Δp marg ; (b) Return temperature in the boiler loop T ret ; (c) Supply temperature in the boiler loop T sup

FIGURE 9 FIGURE 10
FIGURE 9 Return temperatures in HZs 1 to 6 obtained from the one-week simulation of the heating system with the TES system in place: (a) T HZ 1,ret , T HZ 2,ret , T HZ 4,ret , T HZ 5,ret and T HZ 6,ret ;(b)T HZ 3,ret

FIGURE 11 (
FIGURE 11 (a) Total heat demanded to the heat generation units; (b) TES system inlet and outlet temperatures FIGURE A1Marginal differential pressure control scheme . In turn, the hydraulic separator is represented as a single pipeline with a volume V H [m 3 ].The heat supplied by a gas boiler raises the temperature of the water flowing through the hydraulic separator with a mass flow ṁH [kg/s] from an inlet temperature T H ,ret [ • C] to an outlet temperature T H [ • C].Thus, the rate of change in the supply temperature of the hydraulic separator ṪH [ • C/s] is modelled by [W] consuming heat from a hydraulic separator of volume V H [m 3 ], as illustrated in Figure A4.I nt u r n , the hydraulic separator is fed by a motorised control valve with a water mass flow rate ṁH [kg/s].The heat consumed by the HZ decreases the temperature of the water flowing through the hydraulic separator from an inlet temperature T H ,sup [ • C] to an outlet temperature T H [ • C].Thus, the rate of change in the

FIGURE A5
FIGURE A5Control scheme of the return temperature of an HZ

TABLE 1
Daily heat peak-demand of HZs

TABLE 2
Daily heat peak-demand observed by the heat generation units and reduction of the peak-demand obtained with TES system Daily heat peak-