**1. Introduction**

Generation of cold utility for industrial processes consumes about 9% of the net electricity in Germany, with similar trends in other countries, and an overall chiller energy efficiency ratio (EER) of less than two [1]. Cold utility systems encompass cooling towers (evaporative coolers), air coolers, and refrigeration. Air coolers incur near-zero running costs (i.e., "free cooling chiller"), cooling towers require water make-up and chemicals, while refrigeration is normally driven by electricity with substantial running expenses that are approximately one order of magnitude higher than a cooling tower system. Maximizing the use of the lowest-cost utility saves running costs and high-quality energy, reducing environmental burdens. The ongoing digitization of many processes in the industry enables the potential to automate and optimize the operation of cold utility units. Coupled with an efficient and flexible utility system, a site can achieve substantially highier levels of system efficiency. Critical to this vision is linking sectors—e.g., process heating and cooling and electricity grid—combined with technologies, processes, and control systems to collect and analyze the data and implement complex energy optimization methods. Prerequisite for an efficient cooling system is the knowledge and the optimal combination of different operating and performance conditions of individual chillers. The performance of cooling systems depends on their part-load performance and their condensing temperature, which are often overlooked and not continuously measured. Manufacturer performance data often substantially di ffer from in-plant performance [2].

The optimized operation of refrigeration units in combination with a complex control algorithm can improve the energy e fficiency of cooling systems when implemented. Clausen et al. [3] show that the coupling of optimization and simulation of parcel sorting plants lead to e fficiency increases. Dynamic simulations support both the validation of e fficiency concepts prior to implementation and model-based optimization. Augenstein [4] discusses that a linear optimization of the chiller deployment strategy reduces power consumption by 11%. At the same time, there is also optimization potential for specifically adapted expert control strategies in the sequencing of refrigeration systems. Rule-based controls show a good performance for specific systems under normal conditions [5]. Changes in boundary conditions lead to nonlinear behavior of equipment, and equipment degradation can shift the optimum operating set-points and control of a plant over its life cycle [6].

Several approaches are reported in the literature to optimize cooling system operations. Hovgaard et al. [7] use the example of a supermarket refrigeration supply to show that electrical energy savings of 9%–32% are possible with the help of model predictive control (MPC) with a thermal storage system, predicting the changing load and taking into account varying energy prices. Braun [8] reduced the energy consumption with a component-based and system-based optimization algorithm and correlated the power consumption of the chillers, condenser pump, and cooling tower fans with a quadratic function. Olson et al. [9] describe a mathematical approach to the distribution of loads over several chillers to minimize energy costs. Many di fferent methods for optimal chiller sequencing control have been investigated. Some examples are the branch-and-bound method [10] dynamic programming [11], gradient method [12], genetic algorithms [13], simulated annealing [14], and the particle swarm optimization [15]. These methods generally perform well in a simulation environment.

One of the challenges with practical implementation is system dynamics leading to exponentially increased computational demand. To reduce computing time for dynamic optimization, Powell et al. split up the problem solving into a dynamic one for the total load at each time step and a static one for the optimal chiller loading. For a 24 h time horizon, the energy savings are around 9.4% and the energy costs savings are 17.4%. [16] A literature review reveals that the energy-saving potential varies between 5% [17] and 40% [18] depending on the case study.

Another important factor is the availability of cold storage for a cooling system. Kapoor et al. [19] show that implementing a 15,000 m<sup>3</sup> cold-water storage tank could result in a cost saving of 13%. Cold energy storage presents an opportunity to reduce costs by shifting loads. Cole et al. [20] use an MPC with thermal energy storage and variable electricity prices for peak and o ff-peak hours to reduce the electrical energy consumption and the operating costs for the cooling system of one building. In the best case, the annual costs are reduced by around 42% while the total energy usage increases by more than 6%. Kircher and Zhang [21] investigate the cooling strategy of an o ffice building with ice and chilled water storages using an MPC and demand response incentives. The optimization horizon is limited to one day, and the non-linear part-load performance of the chillers is neglected. Ma et al. [22] show that an MPC implementation for a building cooling system considering a storage tank, cooling towers, condensing, and evaporation temperature increases the plant coe fficient of performance (COP) by around 19%. A direct charging and discharging of the storage by the cooling tower is not analyzed. With a similar idea to the current study, Soler et al. [23] present a performance optimization of a district cooling network with eight chillers and a thermal storage tank for a one day period, assuming the use of chillers according to their e fficiency rating and best performance in full load. For minimizing shut-downs and start-ups, no-increase/no-decrease constraints are added.

Many production sites use cold-water storage tanks to smooth the cooling load profile but not for load managemen<sup>t</sup> and load shifting. A sprinkler tank contains a high volume of water for use in emergencies such as extinguishing fires. There is an opportunity to utilize the sprinkler tank for cold-water storage with load shifting while maintaining a minimum water level for emergency operations. Using a free cooling chiller (i.e., cooling by air or cooling water from a cooling tower) offers the opportunity to provide cold water with a higher energy efficiency compared to a compression chilling machine. However, if ambient temperatures equal or exceed the cold water set point temperature, the free cooling chiller is no longer effective. Extending the operation hours of the free cooling chiller and reducing compression chiller's could save a significant amount of electrical energy. This needs a complex pre-charge and discharge control for the cold-water storages. To this end, a mixed integer linear optimization (MILP) algorithm is advantageous.

In practical applications, the performance is highly dependent on the accuracy of the performance models and the effectiveness of the convergence of the optimization algorithm. As demonstrated in the literature review, a current gap is a lack of properly accounting for in-plant chiller performance and considerations for how performance changes through its lifetime. This weakness can be compensated by implementing continuous self-learning performance models and linking the optimization algorithm and dynamic simulations including part-load (off-design).

The aim of this study is to develop and apply a predictive optimization algorithm for multiple compression chillers, a free cooling chiller, and a sprinkler tank as a cold-water storage tank considering part-load system operation under a dynamic simulation to improve overall energy cost efficiency. The algorithm extends the work of Peesel et al. [2], who introduced a predictive simulation-based optimization for a food processing factory, calculating electricity savings of up to 23% per year. The novelty of this study is the evaluation of the interdependency between part-load ratio, variable electricity prices, and a sprinkler tank using data from an industrial energy monitoring system. Furthermore, the influence of the prediction horizon and the ambient temperature on the optimization result is evaluated. The modeling of the utilities, the load profile of cold water, and the framework conditions are based on the energy data of a plastics processing company. For the performance curves of the utility, a continuous self-learning model is used. Additionally, the initial states of the optimization are reset by the simulation results of the previous timestep. The combination of the linear optimization and the dynamic simulation gains the advantages of both methods. This results in more realistic and accurate results for the energy and costs-saving potential of the predictive simulation-based optimization. The developed optimization algorithm is applied to case studies of plastic manufacturers in Germany and Spain to demonstrate its potential.

The remainder of the article is structured as follows. Section 2 presents the general approach of the optimization, the data pre-processing, modeling of the performance curve, and the simulation model. The case study of a plastic processing company is introduced in Section 3. This is followed by the explanation of the results, including the influences of the prediction time horizon, ambient temperature, and variable electricity price, in Section 4. The conclusion is presented in Section 5 together with an outlook for further research.

#### **2. Method and Simulation Model**

This section describes the overall method and simulation modeling approach of the systems using a model-based predictive optimization considering the individual operating performance of the refrigeration supply system by continuously self-learning performance curve models. Figure 1 shows the general optimization method with the energy data from an industrial plant and the self-learning performance curve models. The scheme of the simulation-based optimization is also visualised. The method extends the work of Peesel et al. [24] and comprises four steps. The first step includes collection and pre-processing of the data. The input conditions include the ambient temperature, the residual load, and the electricity price for the production site. At this point, a forecasting model for the weather or the cooling load can be integrated. Weather forecasts are common, whereas load forecasting has been attempted, for example, by statistical methods [25], artificial neural networks [26], and support-vector machines [27]. The energy monitoring system collects the operating data of the cold utility system for all production states. The second step is the modeling of the continuous self-learning performance models for the utilities. This is followed by the predictive optimization for the sequencing of the different chillers. The last step is the simulation of the complete utility system. The different steps are explained in more details in the subsequent subsections.

**Figure 1.** Simulation-based optimization scheme.

#### *2.1. Data Pre-Processing*

Energy and performance data of an energy monitoring system often contain some obvious outliers and incorrect values. This requires pre-processing the data to remove these unusable values and extract the reliable data, such as for cold water temperature (*Tcold*), condensing temperature (*Tcond*), part-load ratio (PLR), energy efficiency ratio (EER), electrical power (*Pel*), and thermal energy flow ( . *Qcool*). The definition of the EER and the PLR are defined by Equations (1) and (2).

$$\text{EER} = \frac{\dot{Q}\_{\text{cool}}}{P\_{el,\text{CC}}} \tag{1}$$

$$\text{PLR} = \frac{\dot{Q}\_{\text{cool}}}{\dot{Q}\_{\text{cool},\text{max}}(T\_{\text{cond}})} \tag{2}$$

The data pre-processing is split into two steps—data preparation based on system behavior and the application of statistical methods. The binary operating state is the variable that helps filter data for times, such as when plants are shut-down. These are up to 50% of all values, depending on the time span of the data such as weekly, monthly, or annual data records. Data during plant shut-downs measure electrical power values <1 kW for the chillers in standby, while cooling load data are often erroneous due to thermal inertia in the fluid. By including these data in the characteristic curve generation, an extremely high EER value results that misrepresents the actual situation. A physical method to remove outliers of the EER is the comparison of the ERR calculated from the measurement data with half of the ideal Carnot efficiency. Other important logical criteria are both the chillers' total electrical power consumption and their transferred cooling capacity. Measured values for these variables must not assume negative values and be lower than the nominal electrical power and cooling capacity at the lowest condensing temperature stated by the manufacturer.

Statistical methods are used to identify incorrect values that pass through the analysis of system behavior. The filtering of values is performed for values smaller than a percentile (1% or 5%) or greater than a given percentile from 95% or 99%. This method is also applied to the rate of change of the recorded measurements to filter out extreme changes, especially for values of the electrical power. A special method is a classification from the minimum to the maximum condensing temperature into one Kelvin class. For this transformation, a matrix approach is used. The direction vector with the largest standard deviation is determined for each classified data matrix for PLR and EER. Outliers are determined in linear or non-linear data point distributions. The vector of the largest scatter indicates a transformed coordinate system for the size pairs of EER and PLR. In this coordinate system, further outliers are identified by the interquartile distance (IQR). In statistical studies, 1.5 times the IQR is often used to identify outliers. The analysis of energy monitoring data of six different compression chillers in three sites by the authors shows that this must be selected much smaller to filter out data. In this work, 0.5 times the IQR is used. The distribution of the data points in the respective classes is decisive for filtering out the values. The various methods are run through one after the other, starting with the method that takes the least time to calculate and filters out the most values. The order is important to minimize the total computation time. After the data pre-processing, the scattering and the number of large outliers for different measured variables were considerably reduced, shifting the median mainly due to the removal of the measured data in the switched off state. The change in the measured values after data pre-processing is shown in the following Figure 2.

**Figure 2.** (**a**) Data overview before and after data pre-processing for the electrical load of chiller 1; (**b**) data overview before and after data pre-processing for the thermal load of chiller 1.

#### *2.2. Continuous Self-Learning Performance Modeling*

After the data pre-processing, the modeling of the load performance curve models is performed. For the free cooling chiller, a second order linear regression model is used. Since the data pre-processing shows that measurement data are available in sufficient quantity and for all operating states, only interpolations are necessary. The following Equation (3) with the regression coefficients *Bx* defines the model:

$$\text{EER}\_{\text{FC}} = B\_0 + \, B\_1 \cdot \text{PLR}\_{\text{FC}} + \, B\_2 \cdot \, T\_{amb} + \, B\_3 \cdot \text{PLR}\_{\text{FC}}^2 + B\_4 \cdot \text{PLR}\_{\text{FC}} \cdot T\_{amb} + \, B\_5 \cdot T\_{amb}^2 \tag{3}$$

The *EERFC* of the free cooling chiller is dependent on the PLR and the ambient temperature (*Tamb* ) for a fixed cold-water supply temperature. For the optimization, the semi-empirical regression chiller model—the Gordon–Ng universal model [28]—based on the evaporator outlet water temperature is used. Equation (4) presents the Gordon–Ng model,

$$\frac{1}{EER} = -1 + \frac{T\_{cmd}}{T\_{cold}} + \frac{-A\_0 + A\_1 T\_{cond} - A\_2 \left(\frac{T\_{cmd}}{T\_{cold}}\right)}{\dot{Q}\_{cool}}\tag{4}$$

where *A*0, *A*1, *A*2 are regression coefficients and *Qcool* is the cooling load of the chiller. In this study, the cold-water temperature (*Tcold*) is kept constant. *Tcond* is the ingoing condensing temperature of the compression chiller.

.

Former studies show that even two structurally identical compression chillers can perform differently in the same location [29]. This results in a need for empirical models for each utility operation. In operating areas with a significant quantity of energy data, empirical models describe the system performance. In other areas, extrapolating semi-empirical models is more accurate. The linking of physical and empirical modeling ensures an acceptable accuracy in both fields of interpolating and extrapolating. For the modeling of the utility, manufacturer and performance data are combined and continuously adapted to form part-load performance curve models. The existing and new data points are weighted differently in the continuous self-learning model. In considering changes in utility performance, recent data are weighted higher. If the characteristic curves are regularly adjusted over a longer period, the deviation between the previous and new characteristic curve decreases. The disadvantage of this method is that an irregularity in the measurement system or in the utilities leads to significant deviations.

For the validation of the Gordon–Ng universal model, the dataset of 1652 data points is split into a training set (82%) and a validation set (18%). The overall deviation of the EER is 8.3% and the root-mean-square error (RMSE) is 0.355. This is because of a low quantity of data for low and high PLR values. Lee et al. [30] show slightly lower values for the validation of the Gordon–Ng universal model.

#### *2.3. Mathematical Model*

The optimization is formulated as a mixed-integer linear programming (MILP) problem. Since the control parameters of the chillers are discrete, the solutions must be integers. The goal of the optimization is to satisfy the cold-water demand by multiple chilling systems for the next NT timesteps at the lowest energy costs. The sprinkler tank offers the opportunity to store cold water. This potential is used to show the load shifting potential of the cooling system by using a variable electricity price as an input parameter of the optimization. The following equations describe the minimization problem.

$$\min \sum\_{t=1}^{N\_T} \sum\_{n=1}^{N\_{\text{ct}}} \sum\_{st=1}^{N\_{\text{St}}} (\text{costs}\_{t,n,\text{ st}} \cdot X\_{t,n,\text{st}}) + \text{SC}\_{t,n} \cdot Z\_{t,n} \tag{5}$$

$$\text{Isub joint to } X\_{t,u,st} \in \{0,1\}, Z\_{t,u} \in \{0,1\} \tag{6}$$

The different part-load operating states (*St*) are between 0% and 100%, while the smallest step size is 1%. The active chilling machines set the lower bound of the operating states and the step size. The cost term (*costst*,*n*,*st*) includes the costs for each part-load operating state (*St*) of the chillers. The costs are calculated by multiplying the electrical power consumption of the chiller by the electricity costs for the current time step. The electrical power consumption is based on the regression models—Equation (3) or (4). The number of chillers is defined by *Ncc*. The number of timesteps (*NT*) describes the prediction horizon of the optimization in hours. For example, an *NT* of 3 by a step size of 3600 s means that the optimal load distribution is found for the sum of loads of the current hour and the following two hours. In this paper, the prediction horizon is varied between 6 and 48 timesteps. Additionally, the starting costs (SC) consider the extra energy demand when starting a machine. SC also helps to prevent an unnecessary starting and stopping of the machines if the gain in efficiency is very low. The active states of the chillers are described by *Xt,n,st*, while *Zt,n* is the binary on/off-variable.

$$\left(T\_{\rm res} \cdot \sum\_{t=2}^{N\_{\rm T}} \sum\_{n=1}^{N\_{\rm St}} \sum\_{st=1}^{N\_{\rm St}} (\mathbf{C}C\_{t,n,st} \mathbf{X}\_{t,n,st}) \right) \le Q\_{\rm pral} + Q\_{\rm plus} \tag{7}$$

$$\text{T}\_{\text{res}} \cdot \sum\_{t=2}^{N\_{\text{T}}} \sum\_{n=1}^{N\_{\text{ct}}} \sum\_{st=1}^{N\_{\text{St}}} (\text{CC}\_{t,n,st} \cdot \text{X}\_{t,n,st}) \ge \text{Q}\_{\text{prod}} - \text{Q}\_{\text{minus}} \tag{8}$$

$$\mathcal{Q}\_{\text{plus}} = \mathbf{m}\_{\text{Tetra}m} \cdot \mathbf{c}\_{p,\text{water}^\*} \left( T\_{\text{return}} - T\_{\text{cold}} \right) \tag{9}$$

$$Q\_{\rm minus} = \mathbf{m}\_{Tcold} \cdot c\_{p,water^\*} \left( T\_{cold} - T\_{return} \right) \tag{10}$$

Equations (7) and (8) describe the constraints for the lower and upper bounds of the optimization. The total cooling supply by the chillers must be lower or equal to the total cooling demand of all timesteps (*Qpred*) and the cooling capacity of the sprinkler tank (*Qplus*). *CC* is the cooling capacity of the chillers. The cooling capacity of the sprinkler tank *Qplus* is defined by the product of the mass of water ( *m*) at the cold water return temperature *Treturn*, the specific heat capacity of water *cp*,*water*, and the temperature di fference of the cold water *Tcold* and returning water *Treturn*. The current potential cooling energy of sprinkler tank *Qminus* is calculated analogously. By varying the mass of water at the di fferent temperatures, the cooling capacity of the thermal storage can be changed. Additionally, the total supply cannot be lower than the total cooling demand minus the current cooling power of the sprinkler tank. For solving the mixed-integer problem, the Gurobi Optimization 8.1 software is used [31].

#### *2.4. Simulation Model*

The results of the optimization are the control parameters for each chilling machine in the current timestep. These parameters represent the operating schedule for the simulation. A dynamic simulation continuously validates the set points in the optimization based on the response and feedback effects from the system. These include thermal inertia of the cooling system, start-up characteristics, as well as the physically modeled heat losses of the storage tank. There is also a di fference in the start-up characteristics for a compression chilling machine in the optimization and the simulation. In the optimization, the full cooling capacity of the chiller is immediately available due to the linear programming, whereas, in the simulation, the chiller needs a few minutes to ramp to full capacity. This is modeled by a first-order delay element. Comprehensive measurements have shown that the ramp time can vary from one to five minutes or even longer for larger machines. During starting time, the cooling capacity and the energy e fficiency of the compression chillers is reduced according to the part-load of the load performance curve. For part-load values below 20% of the full capacity, the EER is below two. This compensates overestimated EER by the linear programming and integrates lower efficiency in starting and shutting down times in the simulation. Therefore, the simulation compensates the limitation of the linear optimization and verifies the operating schedule. Additionally, the complex heat loss mechanisms of a storage tank are modeled more accurately in a simulation than in the linear model. The applied storage model uses the Multi-Node model of the Matlab/Simulink ®CARNOT Blockset. [32] For the calculation of the Multi-Node storage model, the storage volume is divided along its height into *n* nodes. In each layer, an ideal mixing is assumed, i.e., the temperature is constant in the layer element in both the radial and vertical directions. For each layer and each simulation time step, the balance equation, Equation (11), of the energy transport is solved according to the first law of thermodynamics.

$$\frac{d\mathcal{U}}{dt} = \dot{H} + \dot{Q}\_{\text{vertical}} + \dot{Q}\_{\text{loss}} \tag{11}$$

The change in the internal energy, *U*, describes the change in the temperature, *T*, of the individual layers as a function of entering and leaving enthalpy streams . *H* caused by the piston flow in the tank and charging mass flows, direct vertical heat conduction e ffects . *Qvertical* through the layers, and the heat losses . *Qloss* through the container wall. To reduce computational time, the nodes in the model are limited to two—one for the cold-water supply temperature and one for the cold water return temperature.

This method eliminates unrealistic starting points. The initial states of the optimization are reset by the simulation results of the previous timestep. The combination of the linear optimization and the dynamic simulation enables the advantages of both methods. This results in more realistic and accurate results for the energy and costs saving potential of the predictive simulation-based optimization.

#### **3. Plastic Processing Case Study**

The evaluation of the optimization and the simulation is done by a case study on two plastic processing companies that produce injection-molded parts for the food sector. The purpose of the case study is to establish the level of possible energy cost savings for the considered site. The case study can provide motivation to the site to implement the new optimization approach as well as others to follow a similar procedure. This requires high hygiene standards for the production hall. One is located in Germany and the other in Spain. The variation of the location quantifies the influence of the ambient temperature on the results of the optimization. The injection-molding machines are cooled continuously by two di fferent cooling circuits. The mold-cooling circuit is set at a temperature level of about 14 ◦C and the machine cooling at about 30 ◦C. The relatively high process cooling temperature of 14 ◦C o ffers the opportunity to integrate free cooling chillers as well as multiple compression chillers. The refrigeration of the molding circuit is carried out by two air-cooled compression chiller machines and a free cooling chiller, which is also used for winter relief. Cooling towers provide the cooling to machine cooling of the plastic processing machines. Moreover, the company has a 900 m<sup>3</sup> sprinkler tank that can be used as a cold-water storage tank at a temperature level of 14 ◦C. The plastic processing companies often have large warehouses. Fire safety regulations require quick access to a high volume of extinguishing water. Because of the multiple operational states, the predictive simulations-based optimization is well-suited for the mold cooling circuit in the plastic processing industry.

In this study, the prediction horizon is varied between 6 and 48 h and the step size is set to 1 h. The *NCC* is 3 and the *NST* varies from 4 for the compression chillers to 10 for the free cooling chiller. The cooling demand and the ambient temperature are based on the measured data for the year 2017. For this case study, it is assumed that the forecast for the weather and the cooling demand is perfectly known in each time step. The cooling demand varies between 0 and 190 kW. The simulation time (*Tres*) is equal to the step size. The maximum cooling energy stored in the cooling tank is fixed by the maximum temperature di fference of 3 K and the 900 m<sup>3</sup> of water, which results in 3140 kWh. The costs for the optimization are calculated by the product of the current energy consumption of the chillers and the current spot market price for electricity. The prices vary between 75 and 200 €/MWh. The fluctuating price for electricity is also compared to a fixed price scenario. The fixed electricity price is 160 €/MWh. Table 1 gives an overview of the installed cooling utilities and Figure 3 shows the scheme of the cooling system of the plastic processing company.


**Table 1.** Utility system of the plastic processing company.

**Figure 3.** Scheme of cooling system plastic processing company.

The following Figure 4 visualizes the partial load performance curve of the air-cooled chiller, according to the Gordon–Ng model, after data processing. The lower the ambient temperature, the higher the EER because of the lower condensing temperature of the chiller and less energy consumed by the fan of the condenser.

**Figure 4.** Part-load performance curve for air-cooled compression chiller.
