Environmental Health Oriented Optimal Temperature Control for Refrigeration Systems Based on a Fruit Fly Intelligent Algorithm

The recent decades have witnessed refrigeration systems playing an important role in the life of human beings, with wide applications in various fields, including building comfort, food storage, food transportation and the medical special care units. However, if the temperature is not controlled well, it will lead to many harmful public health effects, such as the human being catching colds, food spoilage and harm to the recovering patients. Besides, refrigeration systems consume a significant portion of the whole society’s electricity usage, which consequently contributes a considerable amount of carbon emissions into the public environment. In order to protect human health and improve the energy efficiency, an optimal control strategy is designed in this paper with the following steps: (1) identifying the refrigeration system model based on a least squares method; (2) tuning an initial group of parameters of the proportional-integral-derivative (PID) controller via the pidTuner Toolbox of Matlab; (3) using an intelligent algorithm, namely fruit fly optimization (FOA), to further optimize the parameters of the PID controller. By comparing the optimal PID controller and the controller provided in the reference, the simulation results demonstrate that the proposed optimal PID controller can produce a more controllable temperature, with less tacking overshoot, less settling time, and more stable performance under a constant set-point.


Introduction
Nowadays, refrigeration systems play an essential part in the daily life of human beings. Controlling the temperature by refrigeration techniques is involved in various areas such as human comfort, food storage, food transportation and the environment [1]. However, the systems have to work in the manner of moving the heat from a cold reservoir to a hot reservoir [2], requiring high energy consumption. With the development of the society and the urbanization process, refrigeration techniques are applied everywhere, which causes an acceleration of the growth of the carbon emissions around the world [3]. According to the surveys, almost 30% of the energy consumed all over the world is utilized for Heating, Ventilating, and Air Conditioning (HVAC) [4], while refrigerators account for about 28% of all the energy consumption of a typical US family [5].
Owing to the fact that refrigeration systems are closed cycles, their elements are connected with diverse valves and pipes, which leads to a strong nonlinearity [2], which adds to the dynamic modeling difficulties. Plus, the coupling characteristics increase the difficulty of the design of a controller for the system. All these characteristics of refrigeration systems make their temperature control challenging, and inefficient control strategies may even cause the problem of temperature fluctuation. When the temperature of the environment where people live varies from time to time, it's very likely for those people to get ill, and this is bad for public health.
To solve the problem of saving energy and improving the refrigeration effect, efficient control strategies are of great importance [6]. Several methods have been used for the control of refrigeration systems. In Ma's and Bayram's studies [7,8], fuzzy logic control was applied to control the temperature of a refrigeration system, while in Pedersen's study [9], a neural network is combined with a gain scheduling-based PI controller to control the overheating of a refrigeration system. In addition, in Yin's and Schalbart's studies [10,11], MPC controllers are utilized to control refrigeration systems, and a L-Band SBQP-Based MPC control scheme has also been applied to control two different devices in a supermarket refrigeration system [12]. Although all these control strategies achieve a satisfactory control performance, they have the same disadvantage, which is complexity. For example, to implement a MPC scheme, a great amount of calculation is required, which needs to be performed by a high-performance computer, and this makes the application of MPC schemes hard to realize [13].
Proportional-Integral-Derivative (PID) controllers have been widely used in industrial applications for a long time for their simple structures, accuracy and degree of stability performance [14], thus they can also be used for the control of refrigeration systems. Generally speaking, there are two ways to tune a PID controller, which are analytically and numerically [15]. However, due to the coupling of the refrigeration system, the controller consists of two PID controllers and this constitutes a MIMO process. As a result, there will be six parameters that need to be determined, which makes the controller hard to tune.
The fruit fly optimization algorithm (FOA), proposed by Pan in 2011 [16], is a kind of stochastic optimization algorithm that selects a result with certain rules. With the merits of fast convergence and easy programmability, it is widely used to solve optimization problems [17][18][19]. To this end, a fruit fly optimization algorithm is applied to tune and optimize the parameters of a refrigeration controller in this paper. The position of each fruit fly stands for a set of parameters of the controller, and with iteration, the fly swarm will finally arrive at a location with the best smell concentration [20], which represents the set of parameters possessing the best control performance.
To summarize, this paper: (1) identified the transfer function the refrigeration system model; (2) uses FOA to optimize the parameters of the PID controller based on the identified model; (3) uses the optimal PID controller to control the refrigeration system model. This paper is organized as follows: in Section 2, the model of the refrigeration system is described and the control problems are analyzed, while in Section 3, the transfer function of the refrigeration system is identified, the RGA paring is done to analyze the relationship between the system variables, and the parameters of the PID controller for the identified model are optimized. In Section 4, the optimal PID controller is put into use to control the refrigeration system, and conclusions are drawn in Section 5.

Model Description
The refrigeration system model in this paper, based on vapor compression, is proposed in Bejarano's study [2]. The system consists of a compressor, a condenser, an expansion valve, and an evaporator, and the composition of the system is shown in Figure 1 [2]. The goal of this system is to remove the heat in the secondary flux of the evaporator and then deliver it to the secondary flux of the condenser. The system is based on an inverse Rankine cycle, and it works as follows: firstly, the refrigerant flows through the evaporator under low temperature and low pressure conditions. In this way, the heat in the evaporator secondary flux is removed. Then, the temperature and the pressure of the refrigerant are increased with the help of the compressor, and then it enters the condenser. The temperature of the refrigerant decreases after it flows through the condenser, and during this process, it may become a sub-cooled liquid. Finally the pressure and the temperature of the refrigerant decrease again after flowing through the expansion valve, and the next round of the cycle will start once it enters the evaporator. decrease again after flowing through the expansion valve, and the next round of the cycle will start once it enters the evaporator. As illustrated in Figure 1, this system is a MIMO system, where by manipulating the two variables Av (the opening percentage of the expansion valve) and N (the speed of the compressor), another two variables Te,sec,out (the outlet temperature of the evaporator secondary flux) and TSH (the degree of superheating) are controlled. All the other variables are the disturbances of the system. The range and the initial values of the variables in the system are displayed in Table 1, according to [2].

Expansion
Valve N A v m c,sec T c,sec,in m e,sec T e,sec,in   [21]. With different conditions of the fluid in the heat exchanger, the model can be classified into different modes. As for the evaporator, the modes are classified by the amount of superheated steam as is shown in Figure 2, and the conditions of the condenser are categorized into five modes as illustrated in Figure 3 [22]. As illustrated in Figure 1, this system is a MIMO system, where by manipulating the two variables A v (the opening percentage of the expansion valve) and N (the speed of the compressor), another two variables T e,sec,out (the outlet temperature of the evaporator secondary flux) and T SH (the degree of superheating) are controlled. All the other variables are the disturbances of the system. The range and the initial values of the variables in the system are displayed in Table 1, according to [2]. The model of the system is established by the switched moving boundary (SMB) method, which was first proposed by McKinley and Alleyne in 2008 [21]. With different conditions of the fluid in the heat exchanger, the model can be classified into different modes. As for the evaporator, the modes are classified by the amount of superheated steam as is shown in Figure 2, and the conditions of the condenser are categorized into five modes as illustrated in Figure 3 [22].

Mode 2
Two-phase Figure 2. Modes of the evaporator.

Mode 5
Superheated Figure 3. Modes of the condenser.

Control Problems
The goal of the control design is to maintain the outlet temperature of the evaporator secondary flux Te,sec,out and the degree of superheating TSH following their reference as accurately as possible. The problems to designing the control strategy are listed below: • Strong nonlinearity: Owing to the fact that the refrigeration system is a closed cycle, its elements are connected with diverse valves and pipes, this leads to the result of strong nonlinearity, which adds to the difficulty of dynamic modeling.

•
High coupling: This makes the design of controller for this system complicated and challenging. • Frequent disturbance: This requires the controller to have high robustness and be able to control the system efficiently and accurately to restrain the effect of the disturbance.

•
Constrained control variables: The control variables in this paper is the condenser speed and the valve opening, and they are constrained between 30 Hz~50 Hz and 10~100%, respectively, and this may cause the problem of controller saturation.

Control Design
Owing to the fact that the refrigeration system is a nonlinear system, it's hard to tune the parameters of the controller with the original model. Thus, in this section, the refrigeration system is

Mode 2
Two-phase Figure 2. Modes of the evaporator.

Mode 5
Superheated Figure 3. Modes of the condenser.

Control Problems
The goal of the control design is to maintain the outlet temperature of the evaporator secondary flux Te,sec,out and the degree of superheating TSH following their reference as accurately as possible. The problems to designing the control strategy are listed below: • Strong nonlinearity: Owing to the fact that the refrigeration system is a closed cycle, its elements are connected with diverse valves and pipes, this leads to the result of strong nonlinearity, which adds to the difficulty of dynamic modeling.

•
High coupling: This makes the design of controller for this system complicated and challenging. • Frequent disturbance: This requires the controller to have high robustness and be able to control the system efficiently and accurately to restrain the effect of the disturbance.

•
Constrained control variables: The control variables in this paper is the condenser speed and the valve opening, and they are constrained between 30 Hz~50 Hz and 10~100%, respectively, and this may cause the problem of controller saturation.

Control Design
Owing to the fact that the refrigeration system is a nonlinear system, it's hard to tune the parameters of the controller with the original model. Thus, in this section, the refrigeration system is

Control Problems
The goal of the control design is to maintain the outlet temperature of the evaporator secondary flux T e,sec,out and the degree of superheating T SH following their reference as accurately as possible. The problems to designing the control strategy are listed below: • Strong nonlinearity: Owing to the fact that the refrigeration system is a closed cycle, its elements are connected with diverse valves and pipes, this leads to the result of strong nonlinearity, which adds to the difficulty of dynamic modeling.

•
High coupling: This makes the design of controller for this system complicated and challenging. • Frequent disturbance: This requires the controller to have high robustness and be able to control the system efficiently and accurately to restrain the effect of the disturbance. • Constrained control variables: The control variables in this paper is the condenser speed and the valve opening, and they are constrained between 30 Hz~50 Hz and 10~100%, respectively, and this may cause the problem of controller saturation.

Control Design
Owing to the fact that the refrigeration system is a nonlinear system, it's hard to tune the parameters of the controller with the original model. Thus, in this section, the refrigeration system is identified as a linear model, and two PID controllers will be designed based on the identified model. Fruit fly optimization algorithm (FOA) is used to optimize the parameters of the controller.

Transfer Function Identification
As illustrated in Figure 1, this system is a 2 × 2 MIMO system, where by manipulating the two variables A v (the opening percentage of the expansion valve) and N (the speed of the compressor), another two variables T e,sec,out (the outlet temperature of the evaporator secondary flux) and T SH (the degree of superheating) are controlled. As a result, it can be described as a model with two inputs, two outputs, and four transfer functions. The structure of the identified model is shown in Figure 4. G 11 presents the transfer function from A v to T e,sec,out , G 21 represents the transfer function from A v to T SH , G 12 presents the transfer function from N to T e,sec,out , G 22 represents the transfer function from N to T SH . The process of system identification is as follows: Firstly, a step signal is set on A v , while N is kept as a constant, and the step response of T e,sec,out and T SH are obtained. Then, a step signal is set on N, while keeping A v as a constant to obtain the system response of T e,sec,out and T SH . Finally, the step response curves are analyzed by the System Identification Toolbox of Matlab to identify the system. The comparison of the step response curves is drawn in Figure 5, and we get the transfer functions as follows: identified as a linear model, and two PID controllers will be designed based on the identified model. Fruit fly optimization algorithm (FOA) is used to optimize the parameters of the controller.

Transfer Function Identification
As illustrated in Figure 1, this system is a 2 × 2 MIMO system, where by manipulating the two variables Av (the opening percentage of the expansion valve) and N (the speed of the compressor), another two variables Te,sec,out (the outlet temperature of the evaporator secondary flux) and TSH (the degree of superheating) are controlled. As a result, it can be described as a model with two inputs, two outputs, and four transfer functions. The structure of the identified model is shown in Figure 4. G11 presents the transfer function from Av to Te,sec,out, G21 represents the transfer function from Av to TSH, G12 presents the transfer function from N to Te,sec,out, G22 represents the transfer function from N to TSH. The process of system identification is as follows: Firstly, a step signal is set on Av, while N is kept as a constant, and the step response of Te,sec,out and TSH are obtained. Then, a step signal is set on N, while keeping Av as a constant to obtain the system response of Te,sec,out and TSH. Finally, the step response curves are analyzed by the System Identification Toolbox of Matlab to identify the system. The comparison of the step response curves is drawn in Figure 5, and we get the transfer functions as follows: For all the identified transfer function, the fitness values are all more than 99%. As a result, the identified model can be an ideal linearized model for the refrigeration system.   For all the identified transfer function, the fitness values are all more than 99%. As a result, the identified model can be an ideal linearized model for the refrigeration system.

RGA Paring
The relative gain array (RGA) which is recommended in [23], is a helpful tool for analyzing the interaction of the variables in a system. Setting s to 0, we get the steady state matrix of G 11 , G 21 , G 12 and G 22 in (2): (2) and the RGA matrix is calculated by (3): For the identified system, the RGA matrix is given as (4):

RGA Paring
The relative gain array (RGA) which is recommended in [23], is a helpful tool for analyzing the interaction of the variables in a system. Setting s to 0, we get the steady state matrix of G11, G21, G12 and G22 in (2): and the RGA matrix is calculated by (3): For the identified system, the RGA matrix is given as (4): According to [23,24], the bigger the RGA matrix element is, the stronger the relationship between the system input and output. As is shown in (4), the RGA elements λ11 and λ22 are much larger than λ12 and λ21, this indicates that there is a stronger relationship between Av and Te,sec,out and between N and TSH. As a result, Te,sec,out can be controlled by Av easily, and TSH can be controlled by N easily. According to [23,24], the bigger the RGA matrix element is, the stronger the relationship between the system input and output. As is shown in (4), the RGA elements λ 11 and λ 22 are much larger than λ 12 and λ 21 , this indicates that there is a stronger relationship between A v and T e,sec,out and between N and T SH . As a result, T e,sec,out can be controlled by A v easily, and T SH can be controlled by N easily.

Controller Design
With the merits of simplicity and reliability, PID controllers are still widely used for industrial process control [25][26][27]. The control equation of a PID controller is shown in (5): where u(t) is the control action, K p , K i , K d are the proportional, integral and derivative gain respectively, and e(t) is the tracking error.
For the high coupling of the refrigeration system, two PID controllers are applied to control the evaporator secondary flux T e,sec,out and the degree of superheating T SH , respectively, and the structure of the controller is shown in Figure 6.  Figure 6. The structure of the controller.
Av,ini and Nini are two constants, and they are set at 50 and 40 respectively, which are at the middle of the range of Av and N. Among the controller, there are six parameters we need to decide: the proportional, integral and derivative gain Kp1, Ki1, Kd1 for PID1, and the proportional, integral and derivative gain Kp2, Ki2, Kd2 for PID2.
We use pidTuner Toolbox of Matlab to tune the controller initially, and the initial values of the parameters of the controllers are shown as follows:

Introduction of Fruit Fly Optimization Algorithm (FOA)
Fruit fly optimization algorithm (FOA) is a new method for solving optimization problems which was proposed by Pan [16]. This algorithm is based on the behavior of fruit fly foraging, and the parameters we aim to optimize are set as the position of the fruit fly swarm. The algorithm is executed as follows: Firstly, the position of the fruit fly swarm is initialized randomly. After that, the direction and distance of the movement of each fly is assigned, and the position of the fly swarm is updated. Then, with the calculation result of the judgment function, the fly with the best position is determined, and all the other flies will gather together in this place. Afterwards, we assign the direction and distance of the movement of each fly again, and repeat the process, and the fly swarm will finally reach to the place nearest to the food, in other words, the parameters are optimized.

Tuning of PID Controllers Based on FOA
In this paper, Kp1, Ki1, Kd1, Kp2, Ki2 and Kd2 are set as the parameters of the position of an individual fly, and the block diagram of FOA applied in this paper is drawn in Figure 7. We set swarm size P = 20, and the number of iterations N = 100, respectively. To achieve this scheme, a P × 6 matrix X in (7) is applied to represent the position of the fly swarm, and each position of the flies represents a candidate solution of the PID controller parameters. A v,ini and N ini are two constants, and they are set at 50 and 40 respectively, which are at the middle of the range of A v and N. Among the controller, there are six parameters we need to decide: the proportional, integral and derivative gain K p1 , K i1 , K d1 for PID1, and the proportional, integral and derivative gain K p2 , K i2 , K d2 for PID2.
We use pidTuner Toolbox of Matlab to tune the controller initially, and the initial values of the parameters of the controllers are shown as follows:

Introduction of Fruit Fly Optimization Algorithm (FOA)
Fruit fly optimization algorithm (FOA) is a new method for solving optimization problems which was proposed by Pan [16]. This algorithm is based on the behavior of fruit fly foraging, and the parameters we aim to optimize are set as the position of the fruit fly swarm. The algorithm is executed as follows: Firstly, the position of the fruit fly swarm is initialized randomly. After that, the direction and distance of the movement of each fly is assigned, and the position of the fly swarm is updated. Then, with the calculation result of the judgment function, the fly with the best position is determined, and all the other flies will gather together in this place. Afterwards, we assign the direction and distance of the movement of each fly again, and repeat the process, and the fly swarm will finally reach to the place nearest to the food, in other words, the parameters are optimized.

Tuning of PID Controllers Based on FOA
In this paper, K p1 , K i1 , K d1 , K p2 , K i2 and K d2 are set as the parameters of the position of an individual fly, and the block diagram of FOA applied in this paper is drawn in Figure 7. We set swarm size P = 20, and the number of iterations N = 100, respectively. To achieve this scheme, a P × 6 matrix X in (7) is applied to represent the position of the fly swarm, and each position of the flies represents a candidate solution of the PID controller parameters.
The initial position of the fly swarm is set as K p1 = −20.03, K i1 = −15.21, K d1 = −0.05, K p2 = 5.14, K i2 = 0.21, K d2 = 0.13. In each round, the position of each fly will change in a random direction and distance. The process can be achieved by the calculation of matrix in (8): where R is a P × 6 matrix indicating the radius of the range each parameter is able to change, and rand() is a random number between 0 and 1. In order to make the movable range larger at the beginning to enable the flies to get to the best position earlier, and smaller in the end to guarantee the accuracy, we make R vary with the iteration time. The calculation of R is shown in (9): where R 0 is a P × 6 matrix of the initial number of the radius, τ is the radios adjustment factor which is between 0 and 1, and i is the iteration time. Here we set τ = 0.97, and the value of R 0 is shown in (9). With the position of each fly as the parameters of the controller, we are able to get the dynamic performance of the controller, and calculate the judgment function. We take the control accuracy, the overshoot, and the saturation time into consideration to value the control performance, and we define the judgment function as follows: where e 1 (t) and e 2 (t) are the input of PID1 and PID2, respectively, and u 1 (t) and u 2 (t) are the output of PID1 and PID2, respectively. ω 1 , ω 2 , ω 3 , ω 4 , ω 5 and ω 6 are the combined weights, and α is the adjustment factor. We set ω 1 = 0.7, ω 2 = 0.01, ω 3 = 300, ω 4 = 1.4, ω 5 = 0.01 and ω 6 = 600, respectively, and α = 0.02.
Among all the results, we pick out the position vector with the lowest judgment value, and the assignment of the positions in the next round will be based on this. After hundreds of iteration, the fly swarm will reach to the best position, and the parameters will be tuned.
Among all the results, we pick out the position vector with the lowest judgment value, and the assignment of the positions in the next round will be based on this. After hundreds of iteration, the fly swarm will reach to the best position, and the parameters will be tuned.

Optimization Result
We applied the identified model in 3.1 to execute FOA, and we put step signals to make Te,sec,out change from −22.15 °C to −22.65 °C at 120 s, and make Tsh change from 14.65 °C to 7.2 °C at 120s, respectively. After 100 times of iteration, we finally get the optimized parameters in (11): The comparison of the step response curves is shown in Figure 8, and the trend curves of the change of judgment value and each parameter are shown in Figures 9 and 10, respectively.

Optimization Result
We applied the identified model in 3.1 to execute FOA, and we put step signals to make T e,sec,out change from −22.15 • C to −22.65 • C at 120 s, and make T sh change from 14.65 • C to 7.2 • C at 120 s, respectively. After 100 times of iteration, we finally get the optimized parameters in (11): The comparison of the step response curves is shown in Figure 8, and the trend curves of the change of judgment value and each parameter are shown in Figures 9 and 10    As we can tell from the result, the optimization accelerates the control action of the controller, and reduces the settling time. The judgment value decreases from 765 to 220, and converges at about 220 after 85 iterations, Kp1 swings lower and settles at about −28.5 after 80 iterations, and Kp2 swings higher and settles at about 6.25 in the end. Ki1 keeps decreasing and finally reaches about −45.5, while Ki2 keeps increasing to reach about 1.37 at the end. Kd1 decreases from about 0.05 to 0.063, while Kd2 decreases from about 0.1 to 0.067.

Nonlinear Simulation
In this section, to show the robustness of the proposed controller, simulation of step response based on the proposed controller and the original controller in [2] is implemented, and the result is compared and analyzed.

Simulation Result
During the step response process, Te,sec,out changes from −22.15 °C to −22.65 °C at 2 min, while Te,sec,out changes from 14.65 °C to 7.2 °C at 2 min, from 7.2 °C to 22.2 °C at 9 min, and from 22.2 °C to 11.65 °C at 16 min, respectively. Controller 1 represents the original controller, while controller 2 represents the optimal PID controller in Section 3.3. The simulation result is shown in Figure 11. As we can tell from the result, the optimization accelerates the control action of the controller, and reduces the settling time. The judgment value decreases from 765 to 220, and converges at about 220 after 85 iterations, K p1 swings lower and settles at about −28.5 after 80 iterations, and K p2 swings higher and settles at about 6.25 in the end. K i1 keeps decreasing and finally reaches about −45.5, while K i2 keeps increasing to reach about 1.37 at the end. K d1 decreases from about 0.05 to 0.063, while K d2 decreases from about 0.1 to 0.067.

Nonlinear Simulation
In this section, to show the robustness of the proposed controller, simulation of step response based on the proposed controller and the original controller in [2] is implemented, and the result is compared and analyzed.

Simulation Result
During the step response process, T e,sec,out changes from −22.15 • C to −22.65 • C at 2 min, while T e,sec,out changes from 14.65 • C to 7.2 • C at 2 min, from 7.2 • C to 22.2 • C at 9 min, and from 22.2 • C to 11.65 • C at 16 min, respectively. Controller 1 represents the original controller, while controller 2 represents the optimal PID controller in Section 3.3. The simulation result is shown in Figure 11.

Discussion
From the simulation result, we can tell that the optimal PID controller has a better control performance. The absolute value of the overshoot of the original controller varies from 0 to 11.24, while the overshoot of the optimal PID controller keeps at 0 for most of the conditions. Plus, the settling time of the optimal PID controller (4.20 s~120.18 s) is much shorter than the original controller (49.02 s~169.20 s). In Figure 11b, we can tell that the saturation time of the optimal PID controller is also much shorter than the original controller. In addition, in Figure 11c, the compressor efficiency and the coefficient of performance (COP) of the two controllers are about the same in the steady state, and the optimal PID controller obtains a faster dynamic response.
In conclusion, the proposed controller has the merits of less overshoot, less settling time, less saturation time, and faster dynamic response. As a result, the optimized controller is better than the original controller.

Conclusions
Refrigeration systems are critical for public health and carbon emissions. To deal with the characteristics of nonlinearity and high coupling, this paper solves this problem by applying an intelligent FOA algorithm to optimize the parameters of PID controllers to improve the control performance. This method not only is easy to code and realize, but also retains simplicity and reliability advantages of PID controllers. In this paper, a linear 2 × 2 MIMO system is identified based on the refrigeration system model, and FOA is employed to optimize the parameters of the PID controller for this system. The optimal PID control is finally put into use to control the nonlinear refrigeration system, and the simulation results illustrate that the optimal PID controller is able to regulate the overshoot, reduce the settling time, and restrain the controller saturation phenomenon. As a result, the optimal PID controller has a more accurate control performance and can help the refrigeration system emit less carbon dioxide during the process of dynamic regulation. With these advantages, the optimal PID controller based on FOA is both good for the environment, and public health. Therefore, it can be an ideal controller for refrigeration systems Author Contributions: All authors collectively conceived the research and carried out the analysis. Q.Y. led the simulation and paper writing with contributions and guidance from L.S. and Q.H.