Modeling Supermarket Refrigeration Systems for Demand-Side Management

Modeling of supermarket refrigeration systems for supervisory control in the smart grid is presented in this paper. A modular modeling approach is proposed in which each module is modeled and identified separately. The focus of the work is on estimating the power consumption of the system while estimating the cold reservoir temperatures as well. The models developed for each module as well as for the overall integrated system are validated by real data collected from a supermarket in Denmark. The results show that the model is able to estimate the actual electrical power consumption with a high fidelity. Moreover a simulation benchmark is introduced based on the produced model for demand-side management in smart grid. Finally, a potential application of the proposed benchmark in direct control of the power/energy consumption is presented by a simple simulation example.


Introduction
The increase of renewable energy sources and distributed generation, like wind farms/PVs, creates new challenges for the power system. Power generations with intermittent and unpredictable characteristics create a significant balancing strain and are therefore difficult to integrate in large capacity in the electrical grid. To compensate, the power system will be required to manage actively not only the generation, but also the consumption.
The thermal mass existing in foodstuffs in supermarket refrigeration systems (SRS) makes it possible to store some amount of electrical energy as coldness by lowering the temperatures of the display cases Figure 1. A typical control system structure for connecting supermarket refrigeration systems to the smart grid.
As mentioned above, numerous studies in the literature have been dedicated to modeling refrigeration systems with various levels of details and emphasis on different parts of the system behavior. In the sequel, we shall describe a few, which have inspired the present work in particular. A mathematical model for industrial refrigeration plants has been proposed by [7] for simulation of food refrigeration processes. Only the air temperature of the cold rooms is estimated by the model and there is discrepancy between predicted and measured air temperatures. The development of a Modelica library for CO 2 refrigeration systems have been presented by Pfafferott and Schmitz [8]. The heat exchanger model has been explained, and only the steady-state conditions have been validated by experimental data. Subsequently, they extended their model for transient simulation and validation provided in [9]. The pressures are well estimated, but the air temperature of cold rooms is not reported, and estimation of the mass flow rate is not very accurate. The study of a thermodynamic model supported by experimental analysis has been accomplished in [4]. The method has been proposed to improve the coefficient of performance (COP) by matching the compressor capacity to the load, replacing the classical thermostatic control. A similar work has been performed by [10] to provide an excellent transient characteristics by using a decoupling model. Although the latter two methods can be applied to a single vapor compression cycle (VCR), they are not applicable to the multi-evaporator systems in supermarkets. Reference [5] presents a dynamic mathematical model for coupling the refrigeration system and phase change materials. But this model can only predict the refrigerant states and dynamic COP. A more generalized model has been introduced by [11] based on a two-phase fluid network. The proposed model can simulate different kinds of complex refrigeration systems in different operating modes and conditions, but it can only estimate the steady-state values in different operating modes. Applications of the developed model to estimate the different operating points of the multi-unit inverter air conditioner, heat pump with domestic hot water and multi-unit heat pump dehumidifier have been demonstrated in [12]. It should be noted that the steady-state models without considering cold room dynamics cannot be employed for demand-side managements, in which playing with the cold room temperatures is required for energy solutions. Another such steady-state model is proposed in [6].
The model presented in [13] provides an analytical method for COP optimization. This model does not simulate the system operation and has applicability to the design of the refrigerating machines in terms of selection, design and optimization of the main parameters. With focus on local controls, Schurt et al. proposed a model-driven linear controller as well as its assessment for a single vapor compression cycle in [14,15]. A dynamic model of a transcritical CO 2 SRS has been developed by [16] using the equation-based modeling language Modelica, which can be simulated e.g., in Dymola. Only gas cooler equations and model validation are provided in that paper. As a different perspective, a numerical model has been recently developed by [17] for the evaluation of the use of thermoeconomic diagnosis in transcritical refrigeration.
None of the integrated models in the above cited works can predict the power consumption of compressors that is necessary for demand-side management of supermarket refrigeration systems. There is however a model developed by [18] for simulating the energy consumption of supermarket refrigeration systems. In this work, however, the cold room dynamics have not been modeled and the map-based routing proposed for estimating the compressor power consumption is not very accurate. Considering power consumption in refrigeration systems, Reference [19] has proposed a modeling for optimization purposes. There is however no identification and model validation reported by that work, and the first order model presented for the cold room regarding its air temperature results in storing the energy in the air instead of foodstuffs.
By analyzing the thermodynamics of the system, we develop a model that simulates the electrical behavior. This enables us to control the electrical power consumption (needed for demand-side management) by controlling the cooling capacity (utilizing the existing thermal mass). The proposed model possesses most advantages of the preceding researches like (i) simplifying the behavior of supermarket refrigeration systems as a simulation benchmark for supervisory control; (ii) being very accurate despite its simplicity; (iii) having modularity and reconfigurability; and (iv) being validated by real data. Our modeling is the first principle approach using the physical insight to define the model structure.
The modeling exercise is performed by a modular approach in which the system is separated into different subsystems (modules), each modeled and validated separately. This modularity leaves open the possibility of modeling refrigeration systems with different configurations that are adopted in various supermarkets. Finally, an illustrative example is provided to show how the developed model can be utilized as a simulation benchmark for implementing a simple demand-side management scenario.

System Description
In this work, the case study is a large-scale CO 2 refrigeration system including several display cases and freezing rooms as well as two compressor banks configured in a booster layout as shown in Figure 2(a). The corresponding pressure-enthalpy (P-H) diagram describing the thermodynamic cycle of the system is provided in Figure 2(b). Due to low ambient temperature, the high pressure side operates below the critical value, which keeps the operating condition in the subcritical cycle, as shown in the figure.  Starting from the receiver (REC), two-phase refrigerant (mix of liquid and vapor) at point "8" is split out into saturated liquid ("1") and saturated gas ("1b"). The latter is bypassed by BPV, and the former flows into expansion valves where the refrigerant pressure drops to the medium ("2") and low ("2 ") pressures. The expansion valves EV MT and EV LT are driven by hysteresis ON/OFF controllers to regulate the temperature of the fridge display cases and freezing rooms, respectively. Flowing through medium and low temperature evaporators (EVAP MT and EVAP LT), the refrigerant absorbs heat from the cold reservoir while a superheat controller is also operating on the valves. That is to make sure the refrigerant leaving the evaporators toward compressors is completely vaporized (only in gas phase). Both pressure and enthalpy of the refrigerant are increased by the low stage compressor rack (COMP LO) from "3 " to "4 ". All mass flows at the outlet of COMP LO, EVAP MT and BPV are collected by the suction manifold at point "5" where the pressure and enthalpy are increased again to the highest point, "6", by the high stage compressors (COMP HI). Afterward, the gas phase refrigerant enters the condenser to deliver the absorbed heat from cold reservoirs to the surrounding where its enthalpy decreased significantly from "6" to "7" accompanied by a small pressure drop. At the outlet of the high pressure control valve (CV HP), the pressure drops to an intermediate level and the refrigerant, which is now in two-phase, flows into the receiver and the cycle is completed.

Modeling
There are three major subsystems to be modeled, including display cases, suction manifold together with high stage compressors, and condenser. The modeling of low temperature section is similar to the medium temperature and will be addressed at the end of the modeling section. The operation of the high pressure valve and the receiver are not modeled since the intermediate pressure at the outlet of the receiver is assumed constant.

Display Cases
In display cases, heat is transferred from foodstuffs to evaporator,Q f oods/dc , and then from evaporator to circulated refrigerant,Q e , also known as the cooling capacity. There is however heat load from the environment,Q load , formulated as a variable disturbance. Here, we consider the measured air temperature at the evaporator outlet as the display case temperature, T dc . Assuming a lumped temperature model, the following dynamical equations are derived based on energy balances for the foregoing heat transfers.
where M Cp denotes the corresponding mass multiplied by the heat capacity. The energy flows arė where U A is the overall heat transfer coefficient, h oe and h ie are enthalpies at the outlet and inlet of the evaporators and are nonlinear functions of the evaporation temperature, and T indoor is the supermarket indoor temperature. The heat transfer coefficient between the refrigerant and the display case temperature, U A e , is described as a linear function of the mass of the liquefied refrigerant in the evaporator [21], where k m is a constant parameter. The refrigerant mass, 0 ≤ M r ≤ M r,max , is subject to the following dynamic [20], whereṁ r,in andṁ r,out are the mass flow rate of refrigerant into and out of the evaporator, respectively. The entering mass flow is determined by the opening degree of the expansion valve and is described by the following equation:ṁ where OD is the opening degree of the valve with a value between 0 (closed) to 1 (fully opened), P rec and P suc are receiver and suction manifold pressures, ρ suc is the density of the circulating refrigerant, and KvA denotes a constant characterizing the valve. The leaving mass flow is given bẏ where ∆h lg is the specific latent heat of the refrigerant in the evaporator, which is a nonlinear function of the suction pressure. When the mass of refrigerant in the evaporator reaches its maximum value (M r,max ), the entering mass flow is equal to the leaving one.

Suction Manifold
The suction manifold is modeled by a dynamical equation by using the suction pressure as the state variable and employing the mass balance as [21], where the compressor bank is treated as a big virtual compressor,ṁ dc is the total mass flow of the display cases,ṁ dist is the disturbance mass flow including the mass flow from the freezing rooms and the bypass valve, and V suc is the volume of the suction manifold.V comp is the volume flow out of the suction manifold,V where f comp is the virtual compressor frequency (total capacity) of the high stage compressor rack in percentage, V d denotes the displacement volume, and η vol is clearance volumetric efficiency approximated by with constant clearance ratio c, and constant adiabatic exponent γ [1]. P c is the compressor outlet pressure.
Since the main purpose of this modeling is to produce a suitable model for demand-side management, we need to estimate the power consumption of the compressor bank. The following equation estimates the electrical power consumption,Ẇ whereṁ ref is the total mass flow into suction manifold, and h o,comp and h i,comp are the enthalpies at the outlet and inlet of the compressor bank. These enthalpies are nonlinear functions of the refrigerant pressure and temperature at the calculation point. The constant η me indicates overall mechanical/electrical efficiency considering mechanical friction losses and electrical losses [1]. The enthalpy of refrigerant at the manifold inlet is bigger than that of the evaporator outlet (h i,comp > h oe ) due to disturbance mass flows. The outlet enthalpy is computed by in which h o,is is the outlet enthalpy when the compression process is isentropic, and η is is the related isentropic efficiency given by [6] (neglecting the higher order terms), where c i are constant coefficients.

Condenser
Most of the models developed for condensers need physical details like fin and tube dimensions [2,22,23] and thus are not directly applicable here since our modeling approach is mainly based on general knowledge about the system. So, neglecting the condenser dynamics, the steady-state multi-zone moving boundary model developed in [6] is utilized here with further considerations.
The condenser is supposed to operate in three zones (superheated, two-phase, and subcooled). A pressure drop is assumed to take place across the first zone (superheated) given by [6], where A c is the cross-sectional area of the condenser, and P cnd and ρ cnd are the pressure and density at the outlet of the superheated zone. The first term at the right hand side of Equation (16) indicates acceleration pressure drop and the last term stands for the frictional pressure drop (∆P f ) assumed constant. The Rate of the heat rejection is described by Equation (17) for superheated (first) and subcooled (third) zones, and the following is for the two-phase (second) zone, where U A c is the overall heat transfer coefficient of the corresponding condenser zone, T i and T o are the refrigerant temperature at the inlet and outlet of each zone, and T outdoor is the outdoor temperature. Note that the inlet and outlet temperatures of the two-phase zone are the same when pressure does not change across it. The heat transferred by refrigerant flow across the kth zone is provided by the following energy balance equation:Q in which h i and h o are enthalpies at the inlet and outlet of the kth zone. Accordingly, the total rate of the heat rejected by the condenser is:Q

Parameter Estimation
The system under study is a large-scale CO 2 refrigeration system operating in a supermarket in Denmark. The system consists of seven display cases, four freezing rooms and two stages of the compressor banks configured in the booster setup explained in Section 2. The system is not a test setup and needs to keep its normal operation. Thus, we identify the model from the data collected by every minute logging the measurements. An off-line identification is performed to estimate the constant parameters and coefficients introduced before. The modeling error, computed by dividing the maximum absolute error over the maximum amplitude of variation of the measured signal, is provided on each plot. We use two sets of data: -Training set: This contains the measured input/output data required for the estimations. The data are selected from an interval during the day time when no defrost cycle takes place.
-Validation set: This contains the measured input/output data required for validating the model after estimation. The data are selected in the same hours of interval as the previous one but from a different working day.
After estimating the model using the training data, we plot the simulated response of the estimated model using the validation data for comparison. Thus far we have used the first principle modeling approach to describe the system with mathematical equations. For system identification, we make a nonlinear model using the foregoing equations with unknown parameters. An iterative prediction-error minimization (PEM) method, implemented in System Identification Toolbox of MATLAB, is employed to estimate the model parameters [24]. A modular parameter estimation approach is introduced in which the parameters of each subsystem are identified by providing related input-output pairs from measurement data. An example of estimating nonlinear models using PEM is given by [25].
(1) Display cases estimations: In this subsystem, the model should be able to estimate mass flow and display case temperatures. The mass flows are controlled by thermostatic actions as well as superheat control. These two result in a specific opening degree for the each expansion valve. So the input and output vectors used for estimation are and Y dc = ṁ dc T dc,1 · · · T dc,7 T Estimated parameters for the total seven display cases are collected in Table 1 assuming constant superheat T sh = 5 [ • C], and constant receiver pressure P rec = 38 × 10 5 [pascal]. A frame of five hours' data sampled every minute is used. The inputs are presented in Figure 3. The suction pressure, P suc , is regulated around 26 × 10 5 [pascal] by the high stage compressors, and the indoor temperature, T indoor is regulated around 26 [ • C] by an air-conditioning system. The ON/OFF behavior of the OD signals is caused by the simple hysteresis control actions. During the ON period of action, the superheat controllers operate on the valves to ensure the completed vaporization of the refrigerant at the evaporator outlets. Incorporating the superheat control into the display case models needs a more complicated (e.g., a moving boundary) model for the evaporators. This incorporation seems unnecessary for our control purposes in the supervisory level. Thus for simplicity, we assume a constant superheat temperature and identify the plant directly instead of the closed loop system.  Estimation results for the display case temperatures and the overall mass flow from the expansion valves are illustrated in Figure 4 and Figure 5, respectively. For a fair comparison, all temperature plots have the same scale. Also, the value of the modeling error provided on each plot shows the best and worst fit for the 5th and 7th display case, respectively. The high fluctuation in the overall mass flow is the result of ON/OFF (thermostatic) actions of the expansion valves.
Despite the use of a simple model and relatively large number of estimated parameters (42 parameters), the display case models show satisfactory results. Figure 4. Estimation of display case temperatures. The 5th display case shows the best fit and the worst fit is related to 7th display case, which is still a good estimation. (a) Estimation using training data; (b) Estimation using validation data.  (2) Suction manifold estimations: Although one of the main goals of modeling in this section is to estimate power consumption of the compressor bank, we also need to estimate P suc for this purpose as stated in the previous section.
Suction pressure is computed from Equation (10) by entering the following inputs and output into the identification process and assuming V suc = 2 and knowing V d = (6.5 × 70/50 + 12.0)/3600 from compressor label.
This results in estimating parameters required for volumetric efficiency in Equation (12) as c = 0.65 and γ = 0.47.  It is worth mentioning that the compressor speed is regulated by a local PI controller having a transient faster than the one minute used for data log. So the closed loop response of the compressor speed is considered in estimating the suction pressure and the compressor power consumption. The following inputs and output are chosen to estimate the power consumption of the compressor bank.
The parameters needed for calculating isentropic efficiency are estimated as c 0 = 1, c 1 = −0.52 and c 2 = 0.01, and the mechanical/electrical efficiency is also obtained as η me = 0.65. Figure 6 shows estimation results for the suction pressure and power consumption. Even though the simple first order model introduced for suction manifold cannot generate the high frequency parts of the pressure signal, it can fairly estimate a low pass filtered version of the suction pressure. The bottom plots show a very good estimation of the power consumption in so far as the validation result shows even a better estimation.
(3) Condenser estimations: In the condenser model, the corresponding parameters should be estimated such that the heat transfer generated by the steady-state model has to be equal to the heat transfer delivered by the refrigerant mass flow. The speed of the condenser fan has been fixed to the maximum value, so the developed static model can fairly describe its behavior without incorporating the fan speed. The input vector used for identification is where T i,cnd is the refrigerant temperature at the condenser inlet, which is also the inlet of the first zone. Enthalpies h o,1 = h i,2 and h o,2 = h i,3 are the enthalpies of saturated vapor and saturated liquid at the pressure P cnd , respectively. The output temperature is also calculated by assuming 2 • C constant subcooling. The desired output for estimation is Estimation results are shown in Figure 7, and the estimated parameters are provided in Table 2. The associate result for the pressure drop can justify the assumption that it mainly takes place in the first (superheated) condenser zone. Despite considering a simple steady-state model for condenser and also not using any of its physical details, the estimation is still acceptable and valid.

System Integration and Simulation Benchmark
Thus far, three separate models have been developed for different subsystems using the corresponding measured input and output pairs. In this section, we integrate all subsystems to build a complete model for supermarket refrigeration systems ready for use as a simulation benchmark.
The inputs used in the model are the opening degree of expansion valves (OD i ) and the running capacity of the compressor bank (f comp ).
The disturbance vector is In order to simulate the control strategy depicted in Figure 1, the SRS model should be able to estimate display case temperatures and compressor power consumptions with satisfactory degrees of accuracy. Figures 8 and 9 show the results of running the model by Equations (27) and (28) using the data set used for both identification and validation in the previous section. As can be seen from these results, estimation errors do not increase significantly and the results are still convincing and can satisfy our expectation to have a model as a simulation benchmark. Figure 8. Estimation of display case temperatures after system integration. The estimation error increases slightly due to modeling error associated with each subsystem. (a) Estimation using training data; (b) Estimation using validation data. The modularity of the model makes it possible to simply add the freezing rooms and low stage compressors to the model. There is no need to do changes in the model since the mass flow of this section is already included in the disturbance mass flow in Equation (10). The estimated parameters for the total four freezing rooms are collected in Table 3.  Figures 10 and 11 show the estimation results for the freezing room temperatures and the corresponding power consumption of the low stage compressor bank. The pulse shape of power consumption in the low stage compressor rack is due to the fact that it contains several ON/OFF controlled compressors, in contrast with the high stage rack that contains a continuous frequency control compressor as well as other ON/OFF control types.
So far we have applied the expansion valve control signals (OD) and the compressor capacity (f comp ) to the model from the values obtained from data to have a fair comparison with measured outputs. In order to complete the simulation benchmark, the model should contain these local controllers to be ready for demand-side management as will be explained in the next section. As realized from the measurements, condenser fans are working in their full speed and thus the high-side pressure is not regulated to a specific value in this system. The expansion valve control is a hysteresis controller operating between hysteresis bounds defined based on the temperature limits. We can estimate the compressor capacity such that the suction pressure and low-side pressure are regulated to predefined set-points. For this purpose, the compressor capacity should provide a volume flow in Equation (11) that keeps Equation (10) equal to zero at the desired pressure set-point. The result illustrated in Figure 12 justifies that the proposed estimation method is nearly the same as what is used in the real system.
The disturbances can be modeled by sinusoidal signals with the amplitudes and frequencies equal to the average values for real disturbances. The simulation benchmark is now largely ready to be employed by a supervisory control in order to develop different demand-side management scenarios.

Demand-Side Management
There are two major approaches for demand-side management in smart grid: direct control and indirect control. In the case of direct control, consumers in a smart grid should follow the power reference sent by the grid in accordance with their storage and consumption limitations [27]. On the other hand, in the indirect approach, consumers receive real-time price signal from the grid and manage their consumption such that the operating cost in a specific time interval is minimized [28].
Discussing the flexibilities, and the corresponding storage and consumption limits, and other details in this context is however out of the scope of this paper. In the following, we will show by a simple example how the produced model can be utilized in implementing direct control with the structure explained in Figure 1.
The supervisory controller (see Figure 13) includes a PI controller that regulates the power consumption (sum of low stage and high stage compressor powers) to the reference level received from the grid. The output of the PI controller is downscaled by the gains G i to provide the required differences ∆T i for the set-points for each cold room. This difference values are added to predefined set-points T i and the results are applied as new set-points to the refrigeration system model. The suction and low-side pressures are regulated to constant levels. Figure 14 shows the power consumption in a normal operation when no supervisory control affects the system on the top, and the related direct control at the bottom. The system is simulated for one day. The 60-minute moving average of the power is shown in the plot. A sinusoidal shape of change in the average consumption in normal operation is because of the sinusoidal change assumed for the outdoor temperature.
The objective of the control is to have the average power consumption of refrigeration system follow the power reference while respecting the temperature limits in cold reservoirs. The corresponding cold reservoir temperatures are also depicted in Figure 15. The regulator gains G i are defined such that each display case or freezing room temperature remains in the permissible range required for no food damage. Figure 13. A simple direct control structure. T i are temperature set-points for normal operation of the system, and ∆T i are temperature changes applied for consumption regulation. It is worthwhile to mention that this is only a simple example to show how the produced model can be employed for demand-side management in smart grid. The design of advanced controllers with predictive ability and conducting COP optimization is left for the future works. In order to respect the limitations on food temperatures, a predictive controller can be supported by a state estimator for estimating food temperatures.

Conclusions
A supermarket refrigeration system suitable for supervisory control in smart grid is modeled. The system was divided into three subsystems, each modeled and validated separately. The proposed modular modeling approach leaves open the possibility of modeling refrigeration systems with different configurations and operating conditions. A first principle approach was taken and supported by the PEM method for parameter estimation. The models give satisfactory results in which the electrical power consumption is estimated very accurately by the accurate estimation of the thermal states. These subsystem models were finally integrated to make a booster configuration and the corresponding results confirmed the effectiveness of the proposed modeling approach. It was also explained how a simulation benchmark including the integrated model and local controllers can be configured for the purpose of supervisory control implementations. Using the produced model, appropriate control algorithms can be developed to govern the electrical power consumption (required for demand-side management) by controlling the cooling capacity. At the end, a simple simulation example was provided to demonstrate the utilization of the developed benchmark in demand-side management in smart grid within a direct electrical power control scenario.