Using Multi-Physics Simulation to Estimate Energy Flexibility for Local Demand Response Strategies in a Microgrid

This work discusses the development of a multi-physics simulated model, in the frame of the decarbonization and energy efficiency objectives of the European Commission. Its central feature is the interconnection, through a micro-grid, of a distributed PV installation and of several electric dispatchable loads, thus powering a Collective Self-Consumption network. The simulator presented within this document aims to serve as a technological enabler for the design and testing of On-Site DR strategies, which actuate directly on the connection status of the loads, before their deployment on the target, real-world systems. The simulator supports the design and validation of such strategies by generating realistic simulated data of certain loads that present monitoring difficulties, taking into account online, real external weather conditions. All the elements described and modeled in the current work belong to a real-world installation, which is a university campus —ESTIA, Bidart, France— composed by several buildings with DER.


Introduction
The non-dispatchable nature of renewable sources usually leads to remarkable differences between generation and demand profiles in microgrids and general power grids.These differences must be solved somehow, as a uninterruptible balance between generation and demand must be ensured at all times, mainly by the EMS ruling these grids.Historically, the most common solution to this issue is to use storage systems in order to shift the time of consumption of the energy generated in the grid, giving the system a limited, yet dispatchable, energy source (Jurasz et al., 2020), (Angenendt et al., 2019), (Marańda, 2019).
The renewable generation cannot be largely modified, since these power generators are mostly dependant of weather conditions, which present a stochastic behaviour.The use of short-term energy buffering in storage systems is studied by Marańda et al. (Marańda, 2019), for its application into different scenarios.The figure 1 depicts the typical generation and demand profiles -P P (t) and P L (t) respectively-of residential PV installations.In the case depicted in (b) chart, the storage is considered to cause the buffering of the energy surplus for later use.Following the notation proposed by Marańda et al., E −  B is the amount of energy absorbed by the storage, which is consumed by the system when needed.Contrarily, the supply of previously stored energy into the system is referred as E + B . [a] [b] It is previewed that the deployment of smart systems and the Internet of Things (IoT) will largely contribute to the success of DR programs, since this technology is previewed to allow the remote control of electrical loads (Pop et al., 2018).Duman et al. (Duman et al.) and Romanchenko et al. (Romanchenko et al., 2021) analyze the impact of DR actions achieved through deviations of the set-point of the so-called smart thermostats.Li et al. present an optimal management through DR actions of multi-stakeholder systems in low-carbon communities, considering the carbon emission restrictions through the carbon tax (Li and Yu, 2020).
However, all the mentioned techniques are designed to influence the demand on a macro-scale, or aggregating several end-users.Micro-scale is not the most common aim of the existing DR programs.To overcome the engagement issues of residential DR, some EMSs are designed to shape the performance of the loads, with previous agreement of the user, actuating directly over them to follow more precisely the actions given by the management strategies.These strategies usually generate a schedule for each day and for each load of the system, synchronizing them to reach certain energy and money saving objectives (Vahedipour-Dahraie et al., 2020).The EMS used for these cases computes certain parameters of the loads like ToU (Time of Use), power, etc.Also some user defined parameters, in order to respect the comfort of the users and evaluate the available energy resources to manage the demand levels, scheduling the loads' activation (Mohsenian-Rad et al.).
Nolan et al. point the main barriers and challenges for the deployment of DR programs in their study, giving a notorious importance to the modeling of the physical characteristics of such strategies (Nolan and O'Malley).Another main barrier in this study is the lack of data regarding the disaggregated consumption of some dispatchable loads.This is mainly due to the fact that the energy consumption is usually measured at the installations point of interconnection (POI) with the general grid, while each loads consumption is usually left unmonitored.This is a crucial information for EMS algorithms that need to evaluate at any time the amount of power available to whether increase or decrease the demand levels in order to produce the required DR actions.This is explored by Azari et al. focusing onto the data uncertainty of several DR programs (Azari et al.).Through the use of the simulator presented within this document, we aim to produce the missing data heuristically, thus overcoming the lack of data mentioned before, which is considered a barrier to deploy some DR programs successfully.
In summary, the state of the art in the available literature mainly refers to macro-scale DR strategies evaluated on off-line, fixed weather conditions.A clear gap in this state of the art refers therefore to generic simulation environments, taking into account online weather conditions, aimed at the benchmarking and optimization of generic algorithms for DR, focusing in particular on micro-scale, direct load management.The present contribution has the ambition to address this gap in the literature, by quantitatively validating the performance of the approach presented here, demonstrating the value of micro-scale, online-weather, simulation-based benchmarking and tuning environments for direct load management strategies in a specific real-world case.Furthermore, since the ther-mal devices are external conditions dependant, the present simulation environment allows to evaluate dynamically the energy flexibility that can be achieved through deviations of the temperature set-point.

Description of the Scenario
This study analyses the case of a certain university campus and its electric system interconnecting energy generators and loads.The École d'ingénieurs ESTIA is located in a technology campus called Technopole Izarbel (see figure 2), located on the outskirts of Bidart, near Biarritz and Bayonne, in South-West France.Considering the structure of a microgrid, we can conceive the campus electric system as a microgrid, since it has its own generation and consumption systems, both connected by a distribution network.There are several buildings (called ESTIA1, 2 and 3) that need to be supplied with electric power.These three buildings are foreseen to have their own generation too; so, at some point, the power generated by their generators could be shared between each other, when the EMS of any of the buildings measures a generation surplus, forming a CSC network in the campus.Following the economic constraint mentioned before, the energy is preferentially used by another building of the CSC network instead of injecting it to the general grid.However, we must note that the energy transactions between buildings is also taxed in the French system, since these transactions use the general grid's distribution network to share the energy between buildings.Thus, the energy should be used preferentially by the building where the energy has been generated.However, these taxes belong to a special category -conceived for CSC networks-with reduced taxes, avoiding to constrain the profitability of these CSC networks by its users.

Energy production
The described campus production relies on a local renewable energy source, as mentioned before, which is a PV installation.This installation is composed, at the moment of this study, by 32 PV modules with a maximum rated power of 175 W per module.We could not be provided with further information regarding the characteristics of the certain model of the mentioned modules, inasmuch as this model is not available in the market anymore.Thus, we have used a commercially available PV model that is due to be installed in the future on the campus and whose technical characteristics are provided by the manufacturer.This model is TARKA 120 VSMS 300 which presents a maximum power -on NOCT conditions -of 244.9 W.

Energy demand
ESTIA does not have a very different energy demand profile from an office building, which is composed of typical electrical devices like space heating, lighting etc.The buildings' installation is not equipped with smart meters for each load type, just a general meter at the grid POI and another one dedicated to PV generation.Therefore, it is impossible to find, by measuring methods, which fraction of the total energy demand is requested from each load group (lighting, air conditioning, etc.).It must be noted that the space heating and cooling of this campus are powered by electric air-source heat pumps and not by a natural gas-based system.Thus, as mentioned, both heating and cooling systems rely on electric energy to produce the required heat transfer rates to maintain the indoor air between certain comfort ranges.
The datasets used in this project were obtained from ESTIA1 building's general smart meter, regarding the general demand of the building during two separated periods, with a step time of 1 hour.More specifically, the winter dataset comprises data from January and February 2020, while the summer dataset regards July, August and September 2020.As seen in figure 3, we have computed the average value and standard deviation of each hour, making the distinction between summer and winter periods.Also, we must note that due to the limited available data we will only focus on the modeling of ESTIA1 building.

Dispatchable loads
From a DR perspective, we must identify the loads that could be disconnected from the supply at certain moments, in order to shape the behavior of the building's general demand without disrupting the users' comfort.The so-called dispatchable loads are those that, when required, can totally or partially reduce or increase their energy demand.
In the case of ESTIA1, two loads have been identified as the candidates to perform the DR actions, due to its large power demand and controllability.These loads are: HVAC, water heating system and an electric car battery (BT) charger.

Load Type
Power The ToU of the BT charger is variable between different car models, depending on their BT type and capacity.Also, the charging process of the BTs does not usually have a constant power profile, but again, a constant profile has been used for simplification matters.Instead, the BTs are usually charged with a constant current profile, followed by a constant voltage charge once the BT voltage has reached the set voltage.So, we have assumed an indicative ToU of 3 h and a constant charging power of 11 kW, this one provided by ESTIA.
As we can see in table 1, the ToU of the HVAC system is denoted as variable.This is because the behaviour of these systems depends on the thermal losses of the building, which also depend on the external temperature, a stochastic variable.The control of the heating systems is typically a thermostatic ON/OFF control system with hysteresis thresholding, whose actuation will also vary along with the thermal losses of each day.Furthermore, their thermal contribution varies along with external temperature too.Besides, the HVAC system relies on an air-source heat pump to whether heat or cool the space.Due to this heat recovery from the exterior, its electric consumption is external air temperature dependant, with a maximum consumption at 15 kW.

Modeling of the energy systems
In the present project, we have developed a simulated model of the energy system described in the previous section.The model has been generated using Open-Modelica, an open source simulated environment, mainly powered by the Modelica modeling language, also using the libraries Buildings (Wetter, 2009), PhotoVoltaics (Brkic et al., 2019), and OpenModelica native PVSystems (OpenModelica, 2021).The model, once generated, is subjected to different external conditions (solar irradiation, external air temperature and wind speed), in order to design and test the DR strategies.

PV generation
The PV generation is estimated by its own model separated from the main model, which includes the major part of the energy system of the campus, with the aim of lightening the execution of the main model.This simulated model regarding the PV installation is shown in figure 4. In this model, we are simulating the production of a single module connected to the 400 V side microgrid of the campus via an ideal DC/DC converter, as seen in the diagram.
where  Once T c is determined, the PV module block is fed with this data -in Kelvin-and the irradiance from the me-teoData dataset.The irradiance is limited to a defined range -[0, 1500]-previous to the connection with the PV module, in order to avoid measurement errors, which, in the case of negative values, could cause an error in the electric section.Then, the PV module block injects an electric current into the circuit depending not only on the mentioned inputs but on the voltage on its terminals too.Thus, we are using a mpTracker, along with a power sensor, aiming to find the Maximum Power Point (MPP) at any time, which changes along with the weather conditions.The MPP is a certain voltage value for each weather state, therefore the controller aims to track this value through the whole voltage range.Even though there are many ways to produce this tracking, the tracker used in this model is based on the widely used Perturb and Observe technique (Putri et al., 2015).

Microgrid structure
The electric section is the one in charge of simulating the main loads of the system and, additionally, the energy supply of these loads.This energy supply comes from two different and compatible sources.On the one hand, the microgrid has a distributed energy source with the PV modules and, on the other hand, the system's energy supply also relies on a POI with the general electric grid.We assume that the energy can flow in both directions of this POI, whether demanding energy from the grid or injecting the surplus of the PV generation, when necessary.Both sources and the whole electric section are displayed in detail on the figure 5.The actual electric system uses 50 Hz alternated current (AC), and the power converters required for an AC system normally work at a rate of 22 kHZ, or higher.Also, the iteration frequency of the simulation is recommended to be set 10 times bigger than the highest frequency of the system, which would suppose to work with a time step of T SIM = 4.54 • 10 −6 s.This fact results in an unpractical high computational cost when simulating the system in the order of days, weeks or months, as desired for this study.
Since the analysis of the dynamics of the system is not the aim of this study, the simulated microgrid has been set as a direct current (DC) circuit, allowing to set a time step as high as 15 s.Thus, the results obtained by the simulator regard the steady state of the system, which is fully valid at the context of this study.
We have selected a voltage of 400 V, being this value the typical RMS voltage value of the European threephase grid connection, and assuming that ESTIA is connected to the general grid using this connection type.Based on the fixed voltage of the microgrid, each load has been modeled using a fixed resistor (equivalent to the nominal power consumption of the loads) connected to the 75 DC circuit and it is activated or deactivated using an opening or closing switch, depending on the case.The value of each equivalent resistor has been obtained combining the Ohm's Law and the electric power calculation, getting as result equation 2: where , and P is electric power [W].
We are representing the consumption of all the loads taken into account in this study, which are: the dispatchable loads shown in table 1 -with two HVAC systems due to its use in two different spaces-, and a fixed resistor representing resting fixed consumption, whose calculation is defined next.
As it might be seen in figure 3, the demand average curve never crosses a certain lower threshold.This can be clearly seen at nights, where the demand is closer to this threshold.This lower threshold in the demand can be due to certain parasitic loads that are not switched off during nights or weekends, when the building is empty.We included a fixed resistor to simulate this effect, equivalent to this fixed consumption, computing the average value of the time where the building is unused -from 19 to 6 on weekdays and the 24 hours of weekends and holidays-.We use the same expression as for the rest of the loads to calculate the equivalent resistor (equation 2).

Physical environment of the loads
As mentioned previously, we have considered three load types.Between them, the two heaters consume electric energy to convert it into thermal energy, with the difference that the HVAC relies on a heat recovery system to leverage the thermal energy contained in the external air.Contrarily, the BT charger consumes electric energy from the circuit to subsequently inject it into the car's BT, after converting and adapting the electric supply to certain conditions.We have considered the charger as a constant power consumption that can be switched on and off when needed, so, it can be modeled with a fixed resistor controlled by a switch to model the dispatchable consumption.
In order to model the control strategy of the heating systems, we have used the so-called hysteresis thresholding.This control technique switches the controlled thermal devices on and off, aiming to keep the indoor temperature inside a defined range that can be set according to the user's preferences.This type of control avoids a continuous intermittent switching on and off by setting a temperature range sufficiently wide that allows the temperature to increase or decrease before switching the system again.Thus, the temperature of the plant will cycle between the hysteresis upper and lower limits, ideally never surpassing them.This, in addition to avoiding a continuous and relatively fast switching, saves energy while it ensures the users comfort requirements.

HVAC System's environment
The HVAC system actuation depend on the external weather conditions, both in terms of thermal contribution and electric consumption.Therefore, we included the influence of the external air into the model.For the internal section, we can find room1 and room2 block simulating the internal air closed volume of two separated spaces of the building, along with a heatCapacitor for each of them, simulating the thermal capacitance of any object in contact with the internal air.Each room block has been set with a rough approximation of the volume of the spaces.This approximation has been made using satellite images to measure the side lengths of the building (65x32 m).Besides, assuming a 3 m floor height, we got a volume of 6 560 m 3 for each room.There are three thermal sources associated to each room.The upper heat source -AC ROOM x -represents the main source, which is the HVAC contribution of each room, independent between them.The thermal contribution is interpolated through the data presented in the appendix tables 4 and 5, obtained from (Priarone et al., 2020), which represent the performance of a certain airsource heat pump.Thus, for each external temperature we get a different thermal contribution, along with a electric consumption that is represented in the electric circuit of the model, related to the Coefficient Of Performance (COP) of the HVAC system.The second thermal contribution is related to the irradiance coming from the sun, which is scaled for each room, aiming to simulate different orientations of the rooms.The third contribution comes from the thermal losses with the exterior.Here we compute the external temperature data in K, connecting it with the thermal circuit through a thermalConductor block that represents the thermal insulator of the building's walls.The value of the convection constant G has been assumed taking a fixed amount of heat power losses (1 000 W) at a certain temperature difference with the exterior (3 K).

Environment Programming Interface
The tool presented within this document aims to serve as a intermediate step for any DR strategy previous to its deployment into the actual installations.These strategies are more likely to be developed in a programming language different from Modelica, such as Python or similar.Therefore, we identified the need for a programming infrastructure interconnecting the simulator with the mentioned strategies, which in the context of this project is referred as the Environment Programming Interface, illustrated in the figure 7.As seen in the flowchart of the interface, the starting point is the online weather data obtained from local weather station, which could be historical or forecasted data, depending on the use of the simulator.The download of this data is made by the script Weather Data Downloader, which stores it into a text file to be read by the simulator through a combiTimeTable block.Once this data is stored into the text file, the script Simulation Manager controls the execution of the simulation models through the API conceived for such use called OMPython (Ganeson et al., 2012).Through this API we are able to send the parameters of the simulation to OpenModelica, execute the simulation, and store the results into a user defined CSV file.This file is subsequently read by the Simulation Manager, which extracts the estimated PV current to store it into a text file.Finally, the DR strategies contained in the EMS can be tested controlling the dispatchable loads contained in the MainModel.

Results and Validation
The simulator presented within this document aim to replicate the performance of the energy system analyzed in previous sections, in order to produce the missing data heuristically.Since we were provided with data regarding this installation, we can proceed to evaluate the correctness of this modeling, interpreting the results of the simulations while comparing them with real data of the energy installation.

PV production
In the case of the PV production we are conscious about the differences between the estimation through simulation and the reality of the installation, since we are neglecting the losses due to: non-optimal orientation of the modules on the one hand, and electric inefficiencies on the other hand.Thus, a comparison with the real measured data is crucial to adapt our model more precisely to reality.
The data we are using to validate this model comprises PV production from the period between the 6 th April and 24 th May 2021.Having this dataset, we are executing a simulation of the PV model using weather data relative to the same period, to subsequently save the production data into a CSV file.
The principal aim of this process is to tune a direct relation between data obtained through simulation and real measured data.This is needed to complete the PV simulator since, as mentioned before, the simulator estimates the maximum theoretical production for given weather conditions, neglecting certain power losses.Thus, in order to produce a direct relation, the methodology applied was to tune a polynomial curve that relates the production obtained through simulation and the real data.
Since the simulations are close to be continuous processes, the data obtained from them are close to be continuous too.Thus, each line of the figure 8 represents the trajectory of the PV production power through a whole day.However, as it can be observed in the graph, the related data present high levels of noise.As is evident from figure 8a, sudden changes occur in the production, since there are many lines with horizontal and vertical sectors.This is due to changes in one variable that are not reflected in the other one, which produce the sudden changes in the trajectories.These changes in the production data, whether measured or simulated, must be due to differences in weather conditions too.The only weather variable able to present such high changes is the irradiation, which is highly affected by clouds.Furthermore, the used weather data does not belong to the exact location of the modules -it is measured by a station less than 10 km away-, thus, the clouds would not cover the sunlight the same way in both locations, causing the differences we are observing here.Finally, we tried to fit a polynomial curve to this noisy dataset, getting as result the orange curve of the figure 8a.This approximation can be visually discarded, since it does not follow the pattern of the major part of the trajectories.Considering the clouds a highly chaotic system, instead of taking its influence into account to tune the polynomial curve, we selected several optimal days where the production curves are not distorted by the clouds.These optimal production curves are shown in the figure 8b, where we can see that they follow the pattern of the unfiltered dataset mentioned before.Another fact to mention is that every optimal day used in this graph presented the same two distortion points -located around 3000 and 5000 W respectively in the horizontal axis-, where there is a sudden change on each of them too.Nevertheless, due to its repetition on different days' patterns, we concluded that they are the result of other objects' or buildings' shadows into the PV modules or the irradiance sensors.Splitting the dataset into morning (first half of each day) and afternoon, and using a 4 th degree polynomial fitting, we get the curves presented in the figure 8b, which follow precisely the mentioned pattern.These curves form a function defined in two parts, which are related to morning and afternoon periods and presented in the next expression: where x is the estimated theoretical maximum production data, h is hour, and the constant parameters (A P , G P ) are defined in table 2.
It must be noted that the obtained function is tuned for the period with available data, which is April and May 2021.For a more precise function we should use a dataset comprising the production of the whole year, since the sunlight pattern changes along with the period of the year.However, the period used to tune this function is one of the most productive periods of the year in this certain climate zone, due to high irradiation and cool temperatures.Therefore, it is highly probable that any estimated value at any other time of the year belongs to the range of this function.

Dispatchable loads demand
One of the main objectives of the simulator was to reproduce the performance of certain electric loads.In this section we aim to validate this simulated performance through a comparison with real data.Anyway, as mentioned in previous sections, we cannot be provided with actual data of the modeled loads.Instead, we are comparing the data obtained through simulation with the general demand data.Having in mind that the HVAC systems usually require the major fraction of the energy consumed in buildings -around 67 % of the total energy demand of buildings in France (Enerdata, 2021)-, we can assume that its particular consumption is highly correlated with the general demand measurements.In this case, we are executing a simulation of the general model of the installation, putting the focus onto HVAC demand.In this simulation we are using weather data of the same days as the general demand data of the actual installation, which is previewed to present a similar profile to the simulated data.Figure 9 displays this comparison, showing the data relative to 10 certain days between those with available data.These days belong to the period where the consumption is the highest, which is the winter period, as seen in figure 3. Furthermore, this period comprises several weekdays and a whole weekend, in order to simulate different scenarios.Also, we are adding a threshold to the HVAC consumption data, which is the average value of the consumption when the building is unused, as explained in section 3.2.As shown in figure 9, the general demand profile is highly influenced by the particular demand of the HVAC, as expected.We have computed the Pearson correlation coefficient between these two variables, giving as result ρ Sim,Real = 0.852, which represents a very strong correlation.

Future work
The simulation environment presented within this document aims to replicate the behavior of several loads of a campus.As mentioned, these loads are external conditions dependant as well as several constructive parameters such as the thermal envelope resistance, and the building thermal capacitance.These parameters have been estimated for this implementation of the simulator.However, future work will deal with this issue, estimating these parameters through several on-site experiments which will involve Design of Experiments procedures along with analysis of the resulting data.This calibration process will be based on recent studies on this issue (Wüllhorst et al., 2022).
The results of the present study will be embedded in a OpenADR structure, estimating dynamically the energy flexibility offer for the local consumption controller, or Virtual End Note (VEN) using OpenADR terminology, with the constraint of respecting users comfort in such flexibility offer.The implementation of this structure is an ongoing study which is linked to the present work.

Figure 1 .
Figure 1.Typical generation and demand profiles of residential applications with (b) and without (a) energy storage ((Marańda, 2019))

Figure 3 .
Figure 3. ESTIA's mean power demand profiles with shaded standard deviations' range.
The PV module block computes two inputs to calculate the electrical power produced by them, which are solar irradiance and cell temperature.Solar irradiance is provided by local weather station, while cell temperature is determined computing both external air temperature and wind speed, also provided by local weather station.The cell temperature calculation follows the model presented by Duffie and Beckman in(John A. Duffie, 2013), which is the most accurate model of the analysis made byYang et al. in (Yang et al., 2019) and shown in the equation 1: wind is wind speed [ms], T c NOCT and T a NOCT are cell and ambient temperatures respectively at NOCT conditions [ºC], η c is the conversion efficiency of the PV module and τ α is the product of transmittance-absorbance.The value of these parameters, for the PV module used in the simulator, is shown in table 3 of the Appendix.

Figure 4 .
Figure 4. Block diagram of the PV generation model

Figure 5 .
Figure 5. Block diagram of the electric section of the model

Figure 7 .
Figure 7. General flowchart of the Environment Programming Interface

Figure 8 .
Figure 8. Correlation between estimated and measured PV production.a: using the whole dataset unfiltered, b: selecting two days with optimal production profiles

Figure 9 .
Figure 9.Estimated HVAC demand vs measured general demand comparison

Table 2 .
Parameters of the equation 3