Energy Management of Refrigeration Systems with Thermal Energy Storage Based on Non-Linear Model Predictive Control

: This work addresses the energy management of a combined system consisting of a refrigeration cycle and a thermal energy storage tank based on phase change materials. The storage tank is used as a cold-energy buffer, thus decoupling cooling demand and production, which leads to cost reduction and satisfaction of peak demand that would be infeasible for the original cycle. A layered scheduling and control strategy is proposed, where a non-linear predictive scheduler computes the references of the main powers involved (storage tank charging/discharging powers and direct cooling production), while a low-level controller ensures that the requested powers are actually achieved. A simpliﬁed model retaining the dominant dynamics is proposed as the prediction model for the scheduler. Economic, efﬁciency, and feasibility criteria are considered, seeking operating cost reduction while ensuring demand satisfaction. The performance of the proposed strategy for the system with energy storage is compared in simulation with that of a cycle without energy storage, where the former is shown to satisfy challenging demands while reducing the operating cost by up to 28%. The proposed approach also shows suitable robustness when signiﬁcant uncertainty in the prediction model is considered.


Introduction
Vapour-compression refrigeration systems represent the most used technique worldwide to provide cooling power, including freezing, cooling, and air conditioning applications. The amount of energy involved in such processes is huge, both for industrial and commercial/residential purposes [1]. It is reported that up to 30% of total energy needs around the world are due to Heating, Ventilating, and Air Conditioning (HVAC) systems, among which refrigeration processes represent a significant share [2].
Energy efficiency of current mechanical compression refrigeration systems has been highly improved over time, while their environmental impact has been reduced due to diverse factors. Firstly, more effective design of heat exchangers, variable-speed compressors, and electronic expansion valves (EEV) have improved heat transfer efficiency and controllability. Secondly, more efficient control strategies have been applied to satisfy cooling demand while rejecting disturbances and achieving the highest energy efficiency possible with a direct environmental impact. For instance, Jain and Alleyne developed a with the refrigeration cycle. The prescribed delivery route profile, along with traffic and environmental information, are incorporated into the algorithm to predict the cooling demand. Furthermore, Schalbart et al. present an MPC strategy to the management of an ice-cream warehouse refrigeration system coupled to a TES tank [23]. The controller is based on a steady-state refrigeration cycle model and energy balances on the TES tank and the warehouse.
In this work, a novel hybrid configuration, that comprises a PCM-based TES tank especially designed to complement an existing refrigeration system, is analysed from the point of view of energy management. The layout, shown in Figure 1, is slightly different to that usually applied in the aforementioned works, due to the features of the original refrigeration facility. As depicted in Figure 1, the fluid to be cooled (hereinafter referred to as secondary fluid) is pumped from a tank, which represents the refrigerated chamber, both to the evaporator and to the TES tank, and it is then recirculated to the chamber. An electric resistance is used at this tank to produce heat and simulate the cooling demand, which must be satisfied by both secondary fluid streams. The vapour-compression cycle transfers cold energy to the evaporator by cooling the secondary fluid, but the refrigerant also circulates through the TES tank while charging it, being the latter the cold HTF in this case. The secondary fluid also circulates through the TES tank while discharging it, and thus the warm HTF does not match the cold HTF, as a very relevant difference with respect to the packed bed technology. The layout of the TES tank comprises a number of PCM cylinders and two bundles of pipes, one corresponding to the refrigerant and the other to the secondary fluid, all of them bathed in the so-called intermediate fluid. This TES setup is very similar to that described in a recent work by Bejarano et al. [24], the PCM encapsulation being the only relevant difference. To the authors' knowledge, this setup is novel and the cold-energy management and economic potential offered to industrial refrigeration facilities is first analysed in this work, which represents one of the main contributions of the article. This configuration involves several operating modes, according to the manipulation of the expansion valves and pumps, that enable/disable the different fluid streams and make the system hybrid, in addition to its inherent non-linear features.
As previously mentioned, cold-energy storage allows decoupling of cooling demand and production. Moreover, given the parallel arrangement between the evaporator and the TES tank, the cooling demand at the refrigerated chamber may be satisfied using both secondary fluid streams. It implies that the TES-backed-up refrigeration system is able to satisfy peak cooling demand that might be infeasible for the original standard refrigeration cycle, due to the double contribution provided by both secondary fluid streams. The economic and feasibility potential offered by the addition of the TES tank to the original refrigeration cycle is explored by proposing a scheduling and control strategy based on non-linear MPC to satisfy the cooling demand while minimizing the daily operating cost.
The main contribution of this work regarding the scheduling is the application of nonlinear model predictive control (NMPC) to this hybrid system, where the computational cost of the prediction model is also considered. A non-linear model describing the dominant dynamics of the TES tank is proposed to be used as the prediction model. The latter is developed from the frequency analysis performed in the work by Rodríguez et al. [25], where it is stated that two very different time scales arise in the combined system: one related to the fastest dynamics of the refrigeration cycle, and another, slower one, related to heat transfer within the TES tank. In that work, a detailed model of the TES-backed-up refrigeration cycle was presented, focusing on the faster dynamics caused by the refrigerant circulation, which must be integrated using a small sampling time, in the order of a few seconds. However, such complexity is not affordable or even necessary when addressing the scheduling stage, according to the desired sampling time of several minutes, given the slower dynamics related to heat transfer within the TES tank. Therefore, in this work a simplified non-linear model is proposed, focusing on the dominant dynamics related to the TES tank, which turns out to be far from trivial but more suitable to model-based predictive strategies, especially concerning computational load. This is one of the main contributions of the work, since the reduced computational cost of the prediction model allows direct application of a non-linear model predictive control strategy while ensuring reasonable computation times. The performance of the proposed scheduler for the TES-backed-up system is compared in simulation with the refrigeration cycle without energy storage, while an analysis on the operating cost and constraint meeting is also performed, even when considering significant parametric uncertainty in the prediction model. Therefore, the main contributions of the work are summarised below: • A novel setup of a TES-backed-up refrigeration system is presented, where the TES tank is inserted within the refrigeration cycle, causing the heat transfer fluid to differ between charging and discharging processes. This is an essential and relevant difference with respect to the dominant packed bed technology [13,14]. • A layered NMPC-based scheduling and control strategy is applied to the TES-backedup refrigeration system, which explores the economic and feasibility potential of such a setup, where the control of the compressor and expansion valves is directly affected by the embedding of the TES tank within the refrigeration cycle, in contrast with the packed bed technology, where the refrigerant circulation is not affected by the state of the storage tank [13]. Thus, the existing control strategies previously detailed [17,18,[21][22][23] cannot be applied to the configuration under study. • A simplified and computationally efficient non-linear prediction model describing the dominant dynamics of the system is presented, which allows the application of the previously mentioned non-linear predictive strategy with reasonable computational load, when compared with the accurate but computationally expensive existing dynamic model [25].
The work is organised as follows. Section 2 describes the main characteristics of the designed TES tank and its embedding in the existing facility. Section 3 addresses the modelling of the combined system, focusing on the dominant dynamics related to heat transfer within the TES tank. The overall control strategy is presented in Section 4, focusing on the NMPC-based scheduling algorithm. Section 5 describes a case study for a demanding load profile, where the need for operating mode scheduling is justified according to the power constraints. The economic performance and demand satisfaction of the scheduling strategy are compared to that of the original cycle without energy storage in simulation. Finally, Section 6 summarises the main conclusions and some future work is proposed.

Terminology
All abbreviations, symbols, and subscripts/superscripts used in the work are detailed in Table 1. The CoolProp tool [26] has been used to compute the thermodynamic properties of all involved fluids.

TES Tank Embedding
The designed TES tank is intended to complement a two-compression-stage, twoload-demand experimental refrigeration facility located at the Department of Systems Engineering and Automatic Control at the University of Seville (Spain) [25,27], which works with R404a as refrigerant and a 60% in volume propylene glycol aqueous solution as the secondary fluid. The models and parameters of the refrigeration cycle elements are described in detail in [28]. Figure 2 shows only a section of this facility, including a one-stage, one-load-demand refrigeration cycle where the TES tank is arranged in parallel with the evaporator. Two bundles of pipes run through the TES tank, corresponding to the refrigerant and the secondary fluid, respectively. An additional expansion valve is added to drive the refrigerant from the condenser (labelled as TES EEV in Figure 2), whereas an additional pump is needed to impulse the secondary fluid from the refrigerated chamber (labelled as TES Pump). Notice that the elements defined in [28] as the main compressor, the air condenser, the evaporator related to the refrigerated chamber at −20 • C, and the corresponding expansion valve EEV2 are considered in this configuration, whereas the TES EEV is assumed to be identical to EEV2. Further technical details about the embedding of the TES tank are described in [24].

TES Tank Configuration
The TES tank setup is shown in Figure 3, where the relevant inputs, outputs, and states are highlighted. A number of steel cylinders enclose the PCM, while two bundles of pipes distribute the refrigerant and secondary fluid streams in a counter-current configuration, in such a way that homogeneous heat transfer is achieved. The intermediate fluid (featuring high thermal conductivity and low heat capacity) bathes all cylinders and pipes to homogenise heat transfer.

TES Tank
PCMṁ T ES,seċ m T ES,sec T T ES,sec,in The physical and thermodynamic parameters of the TES tank are included in Table 2, while the thermodynamic properties of the PCM are shown in Table 3. The intermediate fluid is a 60% in volume ethylene glycol aqueous solution featuring high thermal conductivity.  Regarding the inputs, the description of the inlet flows includes extensive (ṁ TES anḋ m TES,sec ) and intensive variables (the {P TES,in -h TES,in } pair for the refrigerant and T TES,sec,in for the secondary fluid). The latter describe the thermodynamic state of the corresponding fluid, whereas the surroundings temperature T surr represents a measurable disturbance. Concerning the outputs, the thermodynamic state of the outlet flows is determined by T TES,sec,out and the {P TES,out − h TES,out } pair. The variablesQ TES andQ TES,sec refer to the cooling powers transferred between the refrigerant/secondary fluid and the intermediate fluid. The TES tank state vector x TES is made up of the temperature of the intermediate fluid T int and an enthalpy description of the cold energy stored inside every PCM cylinder. Notice that the cold-energy storage ratio γ TES (also called as charge ratio) is not actually a state variable, since it can be computed from the cold-energy distribution of the PCM cylinders.

Simplified Dynamic Modelling
In this section a simplified dynamic model describing the dominant plant dynamics as straightforwardly as possible is intended to be developed, so that it can be used as the prediction model in the cooling power scheduling strategy. It is to be noted that the simplified model is used only for the sake of scheduling, while the plant is simulated using the detailed and more realistic model described in [25].
In that work, it was shown that two different time scales arise in the combined system, one related to the faster intrinsic dynamics of the refrigeration cycle, and a slower one related to the states of the TES tank. Indeed, since it is the intermediate fluid that transfers heat to the refrigerant pipes, the secondary fluid pipes, and the PCM cylinders, its low heat capacity is the main cause of those slower dynamics. A detailed model of the combined system, focused on the fastest dynamics caused by the refrigerant circulation, was proposed in [25]. That model also describes the slower dynamics related to heat transfer within the TES tank, but the sampling time needed to be small enough to properly describe the refrigerant dynamics.
Since the scheduling strategy is intended to be performed using a much greater sampling time that fits the dynamics of the cooling demand profile and the TES tank, such a detailed model becomes unsuitable and computationally unaffordable for the scheduling strategy. Therefore, in this work we are interested in developing a simplified version of the previous detailed model, focusing on the thermal behaviour of the intermediate fluid and the energy distribution of the PCM cylinders, that are indeed the slow TES states of the detailed model.
Concerning the slower time scale, since the refrigeration cycle intrinsic dynamics become negligible, the whole TES-backed-up vapour-compression cycle can be statically modelled. Then, the intermediate fluid and its heat transfer to the PCM cylinders represent the dominant dynamics. As a consequence, in the simplified model, the TES tank chargingQ TES and dischargingQ TES,sec cooling powers are considered as virtual manipulated inputs, whereas the refrigeration cycle is implicit and supposed to somehow provide these cooling powers, within some feasible ranges. The simplified model is schematically described in Figure 4, where the surroundings temperature T surr is considered as a measurable disturbance.

Virtual manipulated variables
Simplified TES tank dynamic model T surṙ Then, the simplified model just describes the non-linear dynamics of the intermediate fluid and its heat transfer to the PCM cylinders. The highly time-efficient dynamic modelling approach proposed in [29] is applied, due to the high speed requirements imposed by the predictive scheduler, since the model is expected to run countless times over a certain horizon within the optimization procedure.
In the modelling approach presented in [29], the PCM elements (in this case, cylinders) are divided into a certain number of layers n lay , each one of them featuring the same mass m lay . The number of considered layers represents a trade-off between accuracy and computational load. Some assumptions have been made to develop the time-efficient discrete model. On the one hand, the intrinsic dynamics of the layers in the sensible zone are assumed to be negligible, thus their temperatures are computed according to a steady-state thermal conduction model. On the other hand, cold energy is considered to be transferred at constant rate, provided that the set of layers in latent and sensible zone does not vary.
As mentioned above, the TES tank state vector x TES is made up of the temperature of the intermediate fluid T int and an enthalpy description of the cold energy stored inside every PCM cylinder, which may be defined as indicated in Equation (1): where the specific enthalpies h pcm,j ∀j = 1, ..., n lay of all the cylindrical layers are included in the state vector. The description of the recursive algorithm implementing the simplified dynamic model of the TES tank, given the TES tank chargingQ TES and dischargingQ TES,sec cooling powers as virtual manipulated inputs and the surroundings temperature T surr as a measurable disturbance, is given in a step-by-step sketch, similarly to the original timeefficient discrete model described in [29]: 1. From a given state vector x TES (t), an inward scanning sequence of the PCM layers is performed, looking for the first layer j 0 in latent zone, as indicated in Equation (2): 2.
The thermal resistance R cond pcm,j 0 , corresponding to the spherical shell including all PCM layers exterior to layer j 0 : j ∈ {j 0 +1, . . . , n lay }, is computed as if it was a single, clustered thermal-conduction layer. 3.
The cooling power transferred between the intermediate fluid and j 0 layer is computed as shown in Equation (3): where R cond,wall pcm refers to the thermal conduction resistance due to the PCM cylinder coating, and R conv,ext pcm represents the thermal resistance caused by external natural convection with the intermediate fluid.

4.
Layer j 0 is initially considered to remain in latent zone during the complete time period ∆t, in such a way that cold energy is assumed to be transferred at the constant rateQ pcm,j 0 (t), as shown in Equation (4):

5.
The specific enthalpy of layer j 0 is updated accordingly at the end of the considered time period ∆t, as indicated in Equation (5): 6.
Having reached this point, two possibilities may arise: • Layer j 0 remains in latent zone: h lat− pcm ≤ h pcm,j 0 (t +∆t) ≤ h lat+ pcm . Therefore, the assumption made in Step 4 is confirmed and the specific enthalpy of the layers interior to j 0 is not affected during the time period under study, as described in Equation (6): However, layers exterior to j 0 are shown to be already in the sensible zone.
To compute their specific enthalpy, their temperature is first estimated as indicated in Equation (7), applying a steady-state thermal conduction model, where R cond pcm,j for every layer j is computed in a similar way to R cond pcm,j 0 in Step 2: Then, the specific enthalpy of the layers exterior to j 0 is computed at the end of the time period under study, depending on whether the layer is superheated or subcooled, as shown in Equation (8): The cooling power that the whole PCM cylinder transfers to the intermediate fluid,Q pcm (t), matches the cooling power transferred between the intermediate fluid and j 0 layer,Q pcm,j 0 (t). Eventually, the temperature of the intermediate fluid is updated using the discrete approximation given in Equation (9): where the following sign criteria is considered: • Layer j 0 leaves latent zone. In other words, the latent energy of layer j 0 run out some time before the period under study finished, ∆t j 0 ≤ ∆t, and then the layer j 0 got into sensible zone. The precise instant is estimated, considering a constant rate,Q pcm,j 0 (t), and the remaining latent energy of layer j 0 at instant t, as shown in Equation (10) for a charging process or Equation (11) for a discharging process: Then, a shorter time interval, ∆t j 0 −1 , can be computed as shown in Equation (12): This time interval ∆t j 0 −1 is used for successive estimation, in which the next inner layer, j 0 −1, is considered as the new first layer in the latent zone, and the algorithm restarts from Step 2, but using the corresponding time period ∆t j 0 −1 , instead of ∆t.
Then, the cooling power transferred between each PCM cylinder and the intermediate fluid,Q pcm eventually has several successive contributions along the original time period ∆t, as shown in Equation (13): j 0 −i being the most inner layer found in the latent state during the time period under study ∆t. Obviously, the addition of the partial time intervals matches the complete interval, as shown in Equation (14): For each one of the partial time intervals, T int is updated using Equation (9), giving rise to the succession shown in Equation (15): The proposed model can be expressed as shown in Equation (16), where the function g stands for the complete discrete model, also shown in Figure 4: As stated before, since the dynamic model presented in [25] involves a small sampling time to properly describe the fastest system dynamics, it does not turn out to be suitable or even computationally affordable to be used as the prediction model within a predictive scheduling strategy. This is where the simplified dynamic model proposed in this section, which retains the dominant dynamics while allowing a much greater sampling time, comes handy. Therefore, the computational complexity of the proposed model is not compared with the one in [25] due to impracticability, but the benefits provided in terms of computational cost are evident.

Operating Modes
Diverse operating modes can be defined, according to all possible combinations of the three main cooling powers involved in the problem: the cooling power provided to the secondary fluid at the evaporatorQ e,sec , the TES charging cooling powerQ TES , and the TES discharging cooling powerQ TES,sec [25]. Up to eight operating modes have been defined and discussed in the aforementioned work, though the most useful operating modes might be modes 1 to 4, since modes 5 and 8 require no cooling demand, whereas modes 6 and 7 involve simultaneous TES tank charging and discharging, which may not make sense from the operating point of view. The most interesting operating modes are graphically described in Figure 5.  Mode 1 may be scheduled when the cooling demand can be satisfied by using only the refrigeration cycle and it is interesting to simultaneously charge the TES tank. Mode 2 might be scheduled in the same situation, but when the TES tank is fully charged or it is not economically interesting to continue charging it. Mode 3 is useful when the cooling demand is greater than the maximum power achievable by using only the refrigeration cycle. Mode 4 requires that the demand can be satisfied by only discharging the TES tank, thus the refrigeration cycle is stopped. Figure 6 represents the proposed scheduling and control strategy. The main objective is to track a certain reference of the chamber temperature T re f chamber , which, by means of the so-called outer controller, may result in a reference of the overall cooling power that must be provided to the secondary fluidQ re f sec , since an electric resistance transfers a certain heaṫ Q R to the secondary fluid to simulate the thermal load. This controller is not addressed in this work due to its simplicity, since it tackles a single reference T re f chamber and a single control actionQ re f sec . A simple PID control could be applied without difficulty, thus the development of such a control law has not been considered of interest for this work.  TES,sec , according to efficiency, economic, and feasibility criteria. The cooling power controller is then responsible for getting the TES-backed-up refrigeration system to actually provide the three required cooling powers by manipulating the four manipulated variables available, namely, the compressor speed N, both expansion valve openings A v and A v,TES , and the TES pump. In Figure 6, the secondary mass flow that circulates through the TES tankṁ TES,sec is considered as a virtual manipulated variable, while the actual one is the TES pump speed/power. Anyway, for a certain required value ofṁ TES,sec , a simple regulator is intended to drive the TES pump to actually provide the desired secondary mass flow.

Overview
Therefore, four manipulated variables are available in the cooling power control, whereas only three references must be tracked. These three set points define the desired operating mode, according to Figure 5, in such a way that if one is zero, it involves one of the extensively manipulated variables (A v , A v,TES , orṁ TES,sec ) also being zero. However, when the refrigerant circulates through the compressor (all operating modes considered in Figure 5 except mode 4), the degree of superheating at the compressor intake T SH must be positive, in order to avoid mechanical breakdowns. Since the fourth manipulated variable is thus responsible for getting the degree of superheating T SH over a certain security limit, in this case the control problem is fully actuated.
Concerning the scheduler, given a certain cooling demand forecast, a NMPC-based strategy is applied. In view of the demand forecast and the expected energy prices throughout the day, an operating mode scheduling is proposed, and then the cooling power references are scheduled by solving a non-linear optimization problem. The simplified TES tank dynamic model proposed in Section 3.1 is used as the prediction model, since the scheduling strategy is focused on demand satisfaction and the TES tank cold-energy management, not on the actual production of the required powers, to be addressed by the cooling power controller. Nevertheless, some constraints on the scheduled powers (that turn out to be the decision variables of the optimization problem) must be imposed in order to ensure that the cooling power control problem is feasible and the references are indeed achievable. Additional constraints aimed to keep the PCM cylinders within latent zone are imposed, while the objective function of the optimization problem is based on economic and efficiency criteria.
One important issue is related to time scales. The scheduling strategy is intended to be solved with a sampling time suitable to fit the demand profile, whereas the cooling power control is intended to be performed at the baseline sampling rate, fitting the fastest dynamics related to the refrigeration cycle.

Cooling Power Control
As stated in Section 4.1, the cooling power control is a multivariable problem with four inputs (N, A v , A v,TES , andṁ TES,sec ) and four outputs to be controlled (Q e,sec ,Q TES , Q TES,sec , and T SH ). The decentralised strategy presented in [25] is applied, thus only a few details will be remarked below.
Firstly, a cascade strategy is applied to manipulate the expansion valves, in such a way that the refrigerant mass flowsṁ e andṁ TES are used as virtual manipulated variables. The same strategy would be applied if considering the TES pump power/speed as the actual manipulated variable, in such a way that the secondary mass flow that circulates through the TES tankṁ TES,sec would be considered as a virtual manipulated variable. Secondly, a supervisory strategy is applied for the degree of superheating T SH , when applicable (modes 1-3), to ensure that it never decreases under a certain security limit. Regarding energy efficiency issues, the set point T re f SH is computed following the strategy applied to the original refrigeration facility described in [30], where the compressor speed N is forced to be as low as possible while holding T SH over a certain security limit. Finally, a decoupling strategy is applied to the control of the most coupled cooling powers, i.e.,Q e,sec andQ TES in mode 1.

Scheduling Strategy
As stated in Section 4.1, the scheduling problem is posed as a non-linear optimization, whose main features are detailed below.  (17): • Constraints: -Cooling demand: An overall power forecast is known,Q re f sec (t − 1 + k) ∀k ∈ [1, PH]. Therefore, the constraint described in Equation (18) must be imposed along the prediction horizon: Power feasibility: The maximum and minimum values of all cooling powers must be imposed along the prediction horizon, according to the selected operating mode, resulting in the constraints shown in Equation (19): Notice that, when a certain cooling power is forced to be zero at a given operating mode, the minimum and maximum values for the corresponding cooling power are also zero in Equation (19). The maximum and minimum values of the cooling powers are also shown to depend in some cases on the TES tank state vector [25]; thus, they may vary along the prediction horizon due to the evolution of x TES . -Latent zone: The PCM cylinders are not desired to reach subcooled solid or superheated liquid state at any instant, thus the charge ratio constraints indicated in Equation (20) are also imposed: As mentioned in Sections 2.3 and 3.1, the overall charge ratio can be computed from the specific enthalpies of the cylindrical layers included in the TES tank vector x TES defined in Equation (1). Specifically, Equation (21) allows computation of γ TES from some of the state variables included in x TES : • Objective function: The function J to be minimized in the optimization problem includes a term related to economic cost of cooling power production, as well as an additional term related to the charge ratio, that allows promotion of the TES tank charging/discharging according the operating mode scheduling, to be proposed once the cooling demand profile and the energy price variation have been analysed. Both terms included in the objective function J are detailed below:

Case Study
This section is devoted to the analysis of a case study, focusing on the scheduling problem. A challenging load profile is proposed, where the need for operating mode scheduling is highlighted. The presented simulation results compare the performance of the proposed scheduling strategy to a that of the same refrigeration cycle without TES, regarding economic cost and constraint meeting. A sensitivity analysis is also performed, where the performance of the proposed strategy is assessed when significant parametric uncertainty in the prediction model is considered.

Demand Profile
Consider the cooling demand profileQ re f sec represented in Figure 7. The shape could be similar to a typical supermarket load profile, but the specific power values have been tailored to the system under study described in Section 2. Please notice that the time window has been set to 12 h instead of considering a complete day due to the maximum charging/discharging periods for which the TES tank has been designed. As stated in [24], the latter has been devised to ensure full charge/discharge in 3 or 4 h periods at full charging/discharging power, that were regarded as most desirable due to the research aim of the refrigeration facility. Therefore, longer charging/discharging periods cannot be fulfilled by this TES tank, which would be necessary to address actual 24-h time periods. However, notice that it is only a scaling factor, since, provided that the TES tank volume was high enough, 24-h time periods could be addressed following the same control strategy. Table 4 gathers the maximum and minimum values of the different cooling powers for the system under study [25], regarding operating modes 1 to 4, which have been justified in Section 3.2 to be the simplest and most likely to be scheduled. Notice that a minimum admissible value T min SH = 2°C has been imposed on the degree of superheating when obtaining the cooling power limits in mode 1.

Cooling Power Ranges
As discussed in [25], most cooling power limits depend on γ TES , since the cylindrical shell in sensible zone grows as the TES is charged/discharged. This growth involves higher thermal resistance, that modifies the achievable cooling powers as γ TES evolves. Power ranges when the freezing/melting boundary is located approximately at the cylinder edge, halfway, and at the centre are detailed in Table 4.
It must be remarked that, in the case of mode 2, the cooling demand is provided exclusively at the evaporator, and thus the TES tank is not involved and that is why the cooling power range does not depend on the freezing/melting boundary location. However, in mode 4,Q TES,sec depends on the charge ratio, which decreases as the TES is discharged. Operating mode 3 is merely a combination of the previously described modes, where two independent secondary fluid circuits are used, and thus the cooling power ranges are the same as those presented for modes 2 and 4. Nevertheless, when operating in mode 1,Q e,sec andQ TES are strongly correlated, because the refrigerant circulates through the evaporator and the TES tank, then both flows merge at the compressor intake (point A in Figure 1). Therefore, the ranges indicated in Table 4 are not completely rigorous, since they include some unachievable points. As an illustrative example, Figure 8 shows the steady-state map betweenQ e,sec andQ TES when the freezing boundary locates at the PCM cylinder centre, where the continuous lines represent the overall limits and the small crosses indicate a number of steady-state achievable points.   Similar plots have been obtained for different locations of the freezing boundary, which are not included here for the sake of brevity. Notice that the cooling demandQ sec in operating mode 1 must be completely satisfied byQ e,sec , sinceQ TES,sec is zero. Therefore, the maximum charging powerQ TES is shown to be constrained by the cooling demand.

Energy Price
Realistic energy price evolution throughout the day is applied to the optimization problem. Figure 9 depicts the electricity prices corresponding to a given day (29 May 2019) in the Spanish spot market [31].
These prices, once scaled, are considered as the cost of the cooling power production at the evaporator, w e,sec in Equation (22), as well as the cost of the TES charging power, w TES . However, w TES,sec is assumed to be zero, given that no extra cost is involved when using the previously produced cold energy, beyond the electric power devoted to pumping the secondary fluid, which is not considered in the objective function detailed in Equation (22).

Operating Mode Scheduling
If the demand profile shown in Figure 7 is analysed together with the cooling power constraints shown in Table 4 and the energy prices shown in Figure 9, it can be observed that, approximately, during the first 5 h, the cooling demandQ sec can be satisfied by using onlyQ e,sec , given the power ranges for modes 1 and 2 presented in Table 4. It may be interesting to charge the TES tank during this period, taking advantage of the low energy price, but in any case the maximum charging power is constrained by the cooling demand, which must be satisfied in mode 1 by the refrigeration cycle. However, from t = 5 h, the cooling demandQ sec is too high to be satisfied only byQ e,sec , being necessary to turn tȯ Q TES,sec to complement the cooling power produced at the evaporator. Then, during the peak period until t = 8 h, the TES tank is intended to be discharged, avoiding as much as possible the direct cooling power production due to the high price corresponding to this period. That is when the high charge ratio achieved previously comes handy. Finally, at the end of the day, the cooling demand is again low enough to allow modes 1 and 2 to be scheduled, thus the TES tank might be also charged during this off-peak period.
In the light of the previous analysis, a certain operating mode scheduling must be proposed to make use of the TES tank as a cold-energy storage unit, as efficiently as possible, according to economic criteria. This operating mode scheduling affects the optimization problem posed in Section 4.3 through the feasibility constraints described in (19), as schematically shown in Figure 10.   Figure 11 shows a proposed operating mode scheduling for the specific demand profile introduced in Figure 7, according to the power constraints analysed in Section 5.2 and the energy price profile shown in Figure 9.   As observed in Figure 11, the TES tank is scheduled to be charged during the off-peak periods at the beginning and at the end of the day, when the cost of cooling power production is lower. Given that the demand is low but not zero during these periods, operating mode 1 is scheduled, in such a way that the cooling demand is satisfied exclusively at the evaporator. Furthermore, the cold-energy stored in the TES tank is scheduled to be released during the peak period, complementing the cooling power provided to the secondary fluid at the evaporator to meet the cooling demand (operating mode 3). Two intermediate periods when the TES tank is not charged nor discharged, but in stand-by (operating mode 2), have been inserted between the charging and discharging periods, given that the cooling demand can be satisfied at the evaporator and it is certainly more interesting from an economic point of view to ensure the cooling demand satisfaction during all the peak period, when the contribution of the TES tank is imperative. Moreover, once the peak period is finished, it might be more interesting to concentrate the TES tank charging during the central hours of the off-peak period, when the economic cost of cooling power production is lowest.

Simulation Results
Some simulation results of the proposed scheduling strategy are presented in this subsection, where they are also compared to the performance of the original refrigeration system without energy storage, regarding both economic cost and constraint meeting.
Given the cooling demand shown in Figure 7 and the energy price profile shown in Figure 9, the references of the cooling powers are computed every hour. However, the dynamics related to heat transfer within the TES tank require an internal sampling time of the prediction model of 10 min. The latter could be reduced to increase accuracy, at the expense of greater computational load of the predictive scheduling strategy. Safety limits γ min TES = 0.05 and γ max TES = 0.95 have been considered, while a prediction horizon PH = 4 is set. At the initial state of the TES tank, the intermediate fluid is assumed to be in thermal equilibrium with the PCM cylinders (T int (t = 0) = T lat pcm ), whereas the initial enthalpy distribution of the latter is such that γ TES (t = 0) = 0.35. Figure 12 shows the performance of both the proposed scheduling strategy for the TES-backed-up system and the refrigeration cycle without TES, concerning the satisfaction of the cooling demand. Furthermore, Figure 13 shows the three references of the relevant cooling powers, whereas Figure 14 represents the charge ratio of the TES tank in the proposed scheduling strategy. Moreover, Figure 15 shows the temperature of the intermediate fluid, as well as the distribution of the temperature field of the PCM cylinders in the cooling, charging, and discharging processes. Figure 12 shows that the refrigeration cycle without energy storage is not able to satisfy the cooling demand peak from t = 5 h to t = 8 h, thus incurring in some cooling power deficit, which is remarked in Figure 12. However, the TES-backed-up refrigeration cycle operated according to the proposed scheduling strategy is shown to provide the required cooling power to the secondary fluid by combining the power produced at the evaporator and the TES tank contribution during the peak period. To do this, the TES tank is charged during the off-peak periods, as shown in Figure 13b. The charge ratio is shown to be kept within the safety limits, as depicted in Figure 14, whereas a great deal of the energy capacity of the TES tank is used along the day. Notice in Figure 15a that, during the charging processes, the temperature of the intermediate fluid is below T lat pcm , being over T lat pcm in the discharging process; in the stand-by processes, the thermal inertia of the intermediate fluid causes that, while approaching the thermal equilibrium with the PCM cylinders at T lat pcm , the latter continues charging/discharging while the residual cold energy is transferred. The distribution of temperatures within the PCM cylinders represented in Figure 15b show how the layers are quitting the latent zone from the outermost layer to the innermost, both in charging and discharging processes, as their latent energy depletes.
The performance of the low-level cooling power controller is not analysed in this work for the sake of brevity, since it is analysed in depth in [25], thus no plots of the actual control actions are included in this work. In any case, some important issues concerning this control layer will be remarked upon. The settling time of the cooling power closed loop turns out to be small enough to assume, given the scheduling time, that the required cooling powers are almost immediately provided, as long as the computed references are feasible, which is ensured by the power constraints met by the scheduler. The hypothesis about the separation between the time scales considered in the modelling stage is confirmed in the aforementioned work. Concerning the charging and discharging cooling powers, constant cold energy supply/release is achieved by increasing the mass flow of the refrigerant/secondary fluid as the thermal resistance caused by the cylindrical shell in sensible zone grows. Further details about the performance of the cooling power controller can be found in [25].
Regarding economic cost, the overall energy costs along the day of both the proposed scheduling strategy for the TES-backed-up cycle and the original system are gathered in Table 5, as well as the achieved reduction. Notice that, in the case of the refrigeration cycle without storage, it is assumed that the energy deficit shown in Figure 12 must be bought at the market price detailed in Figure 9.     Obviously, the achieved reduction on the operating cost would come at a certain price in terms of CAPEX (capital expenditure). The economic analysis performed in this work is not intended to be comprehensive, but to throw some light on the high economic potential of such energy storage systems, which not only can ensure the satisfaction of demanding cooling load profiles, but they can also lead to significant cost reduction when properly operated. However, the simulation results highlight the need for suitable predictive scheduling strategies to manage the stored cold energy within the latent zone and fit the cooling power limits of a given facility. Indeed, the latter have been shown to vary significantly when the energy storage system is added to the original refrigeration cycle, both in charging and discharging processes.

Sensitivity Analysis
Some simulations of the proposed scheduling strategy considering parametric uncertainty are presented in this subsection. Given the design of the TES tank described in Section 2, minor uncertainty is expected to be devoted to geometry, constructive materials, working fluids, and their thermodynamic properties. Nevertheless, a simple heat transfer model is used to compute the cooling power transferred between the PCM cylinders and the intermediate fluidQ pcm , where a probably uncertain heat transfer coefficient computed by classical correlations for natural convection [32] has been considered to compute the thermal resistance R conv,ext pcm in Equation (3). Moreover, the assumption of homogeneous heat transfer between the PCM cylinders and the intermediate fluid could lead to some degree of uncertainty regarding the cooling power calculation and then the evolution of the temperature of the intermediate fluid T int . In order to assess the robustness of the proposed scheduling strategy against this kind of uncertainty, the latent temperature T lat pcm considered by the prediction model will be artificially modified. The latter is usually precisely known by design, but if an uncertain value was considered in the prediction model, it would strongly affect the calculation of both TES charging and discharging cooling powers, which in turn determine how the intermediate fluid evolves. Therefore, the consideration of some uncertainty in the latent temperature T lat pcm is just a simplified way of inducing uncertainty in all cooling power calculation within the TES tank. Figure 16 shows the references of the TES charging/discharging cooling powersQ TES andQ TES,sec , when some uncertainty in T lat pcm is considered. This temperature has been increased/decreased in 1 K from the nominal value in the prediction model. Given the variations of the intermediate temperature for the nominal case shown in Figure 15a, the considered increment/decrement is high enough to cause significant uncertainty in the prediction model with respect to the accurate model used for the plant simulation. Furthermore, Figure 17 compares the TES charge ratio for the three cases.
It is shown in Figure 17 that the proposed scheduling strategy provides very similar results for the references of the charging and discharging cooling powers, in spite of the significant uncertainty in the cooling power calculation and thus in the estimated charge ratio. In the case ofQ TES,sec , the references are identical, as shown in Figure 16b, but some appreciable differences appear in Figure 16a regarding the charging power. The latter variations are due to the greater closeness to the upper limit imposed on the charge ratio in the case of the charging process, as represented in Figure 17. Therefore, the proposed scheduling strategy presents a suitable degree of robustness to significant modelling errors, ensuring cooling demand is satisfied and that the limits of the latent zone are met even if significant errors are introduced in the cooling power calculation.

Conclusions and Future Work
In this work, cold-energy management of a novel setup including a thermal energy storage tank based on phase change materials that complements a refrigeration facility has been addressed. The design of the upgrading process undertaken on the refrigeration facility to include the storage tank has been summarily described, along with the main features of the storage unit. Concerning the modelling, from a very detailed dynamic model previously developed, a simplified model focused on the slower dynamics related to heat transfer within the storage tank has been proposed to act as the prediction model in the scheduling optimization.
A layered scheduling and control strategy has been proposed, where a non-linear predictive scheduler computes the references of the main cooling powers involved, whereas the low-level controller ensures the achievement of the required cooling powers. The predictive scheduling problem is posed as a non-linear optimization with constraints due to cooling demand satisfaction, power feasibility, and latency limits, while considering economic criteria in the objective function.
A case study has been analysed, where a challenging load forecast that requires scheduling of different operating modes must be satisfied. The proposed strategy has been shown to ensure the cooling demand satisfaction and meeting of constraints while reducing the daily operating cost by up to 28% when compared to the refrigeration cycle without TES. Moreover, the latter fails in satisfying the cooling load during the entire day. A sensitivity analysis of the proposed strategy has shown it to provide satisfactory performance even when significant uncertainty in the prediction model is considered.
As future work, the application of the proposed strategy to the experimental plant is scheduled to be performed as soon as the upgrading process undertaken on the facility is finished. Furthermore, the development of an economic NMPC strategy is devised [33], as well as the stability analysis of the proposed NMPC-based scheduling strategy.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: