Dynamic Study and PI Control of Milk Cooling Process

— The background of this research is to understand the operation process, which is the main goal of developing the process model. This model is often used for operator training, process design, safety system analysis, or control system design. The dynamic model of the milk cooling process from 36˚C to 4˚C using chilled water available at 2˚C was performed. Chilled water was maintained at a constant temperature by using a refrigerant unit. The process being investigated was a Packo brand milk cooling tank belonging to KUD SAE Pujon (Malang - Indonesia). A fundamental heat balance method was used to derive the model, leading to a first-order transfer function process. For a 2-hr cooling process, the gain and time constant values are 1.00 and 42.3548 mins, respectively. Heat balance was then extended to continuous processes so that its transfer function could also be obtained. This study simulated and investigated the behavior of batch and continuous processes. Process Identification via input-output data was also introduced for continuous process. The process model obtained from the system identification toolbox was very useful in control, such as for determining tuning parameters via the Ziegler-Nichol method for Proportional-Integral control. However, a small delay was required to be introduced to the system as the first order process without time Ziegler Nichol method cannot be implemented. Further research may include other system identification methods, such as ARX, ARMAX, Output-Error, Box Jenkins etc., or implementing advanced process control for milk cooling.


I. INTRODUCTION
Milk is a white liquid food rich in nutrients mammalian milk glands produce.This is the main source of nutrition for baby mammals, including humans, before they can digest other types of food [1].Early lactation milk contains colostrum [2], which carries maternal antibodies to its young and can reduce the risk of many diseases [3].It contains many other nutrients, including protein and lactose [4], [5].Consumption of inter-species milk is common, particularly among humans, who consume other mammalian milk [6].Milk is a good medium for the growth of microorganisms [7], especially pathogenic bacteria such as Staphylococcus aureus, Escherichia coli, and Salmonella sp [8], [9].Microorganisms in raw milk can thrive if stored at room temperature [10], [11].These microorganisms will be dangerous if consumed, so raw milk that has been milked must be stored at a certain temperature [12].It aims to inhibit the growth of microorganisms that can reduce the quality of fresh milk.
Methods of milk collection differ across countries, and there is no single, universally applicable scheme.Milk is highly perishable due to bacterial activity, and therefore, it must be collected and cooled soon to avoid damage due to deterioration and to maintain the standard.Since milk production is often remote from markets and treatment facilities, in Indonesia, the milk cooling center at Koperasi Unit Desa (KUD) provides the facilities to maintain quality through cooling and hygienic storage before further transportation to the milk factory [13].
The chilling of raw fresh milk to below 4 °C as soon as possible after milking and within two to four hours at the most -is recognized as the best way of maintaining standard and reducing milk damage.Therefore, it was done in Cooperatives (Koperasi Unit Desa) as a batch cooling process and this process is performed twice a day (in the morning and afternoon).However, in a milk factory, milk cooling is usually done in a continuous process.This cooler needs to be installed with a controller to ensure the temperature of the milk comes out 4 °C [14].
The main problem that often occurs in dairy factories is the operator's lack of understanding of the process, which creates difficulties in controlling the process or in designing appropriate equipment when needed.Therefore, this paper aims to present the dynamic model of the batch and steady-state milk cooling process.The model in Laplace transform is important for controlling the process, and understanding the process operation is a major overall objective for developing a dynamic process model.This model is often used for operator training, process design, safety system analysis or design, or control system design.
Milk cooling is one of the largest consumers of electricity [15]- [17]; therefore, milk production is an exhaustive energy process.Milk cooling can be done in two steps: precooling using water and refrigerant cooling to 4°C.Cooling costs at a Milk Cooling Center can be decreased by precooling fresh warm milk using water from the main supply, surface, well, or groundwater.Precooling reduces the refrigeration burden and hence reduces energy needs and costs.When the temperature difference between the water and the fresh raw milk is significant, refrigeration costs can be reduced by up to 64 percent [13].However, in tropical countries like Indonesia, we normally skip the first step since the water temperature and fresh raw milk are about the same.The energy required for milk cooling depends mainly on the efficiency of the cooling equipment and the temperature difference in milk between the initial and end of the process [15].
Chilled water can be used to cool the milk in the bulk tank cautiously.Water can be chilled using a chiller or ice bank equipment, which usually consists of an insulated water tank with an evaporator and a holding vessel.The water is chilled to approximately 2˚C, and the milk is cooled to 4˚C.Another method of chilling water is circulation through an ice bank.Ice builds up around the copper tubes in cylindrical formation.Water is circulated in a closed loop through the bulk tank and back to the ice bank [15].When milk is stored above 30˚C, it can spoil after 4 h, depending on its total bacteria count.Storing at 20˚C, 10˚C and 4˚C can partially inhibit bacterial growth for the first 8, 16, and 24 h after milking, respectively [18].
Murphy et al. [15] studied feed forward and feedback controller for milk cooling system for two steps cooling (precooling and refrigerant cooling) and tuned using Ziegler Nichols ultimate gain method [19].However, no information for the plant models was provided.Coughanowr and LeBlanc [19] presented a fundamental method for the obtaia mercury thermometer's first-order transfer function.This paper adopted the method for derivation of batch milk cooling transfer function, which is similar in nature.
Investigations of waste heat recovery from bulk milk coolers [20] have also been conducted.Heat dissipated into the atmosphere through the condenser was recovered to improve the plant's energy efficiency.The waste heat can be utilized to heat the water which was used to clean the milk processing equipment.Another modeling work that focuses on an ice bank for milk cooling after milking has also been performed [21].Agustriyanto et al. [22], [23] provided a transfer function model for continuous milk cooling, while this paper also derived a batch cooling system.Operators in dairy factories often do not understand the process that occurs, so they cannot control or design the equipment according to their needs.This is a major problem in dairy factories.
The models addressed in this paper were based on fundamental theories or laws, such as the conservation of energy [24].Simulink was also used to simulate the model.Simulink is a block diagram environment for multi-domain simulation and Model-Based Design.It supports system-level design, simulation, automatic code generation, and continuous test and verification of embedded systems.Simulink provides a graphical editor, customizable block libraries, and solvers for modelling and simulating dynamic systems.It is integrated with MATLAB, enabling us to incorporate MATLAB algorithms into models and export simulation results to MATLAB for further analysis.Therefore, it is easier for us to generate and analyze the data.

II. MATERIALS AND METHOD
This section describes the method for mathematical modeling of batch and continuous milk cooling process and Proportional Integral (PI) control of continuous milk cooling because this study was a simulation.However, Pujon Malang's survey was conducted to collect equipment data at KUD SAE.

A. Mathematical Method for Batch Process
The system being studied was the milk cooling facilities belonging to KUD SAE Pujon, as shown in Figure 1.It is a rectangular milk cooling tank with a half-cylinder tank contained therein.A rectangular tank is a place for cooling water, while a half-cylinder tank is a container for milk.The scheme of the milk cooling tank is shown in Figure 2. Table 1 shows batch and parameter data for simulation.
The tank capacity (V) is 2500 L. Therefore, the mass of milk in the tank can be calculated since the milk density is known.The volume of chilled water is unknown, but its temperature is kept constant at 2˚C (T1) by using refrigerant.If we set the cooling time equal to 2 hours, then the total heat that must be released from milk to cooling water can be calculated as follows: Log mean temperature different (LMTD) for parallel flow then can be calculated [25], [26], [27], [28] : Batch milk cooling process can be modelled based on energy balance as follows [19]: Unsteady state energy balance: @ 5 @ EF "GGH HI$'JK% When Qin = 0 , then: Equation ( 7)-( 8): Taking Laplace Transform of both sides in Equation ( 9): When To(0)=0, then: The resulting Equation ( 11) is a first order process transfer function [22] [29]with: Hence:

B. Mathematical Method for Continuous Process
The same tank was then operated continuously.It was assumed that the volume of milk or the mass of milk in the tank was constant so that the mass balance could be neglected.Table 2 shows the parameters at steady state conditions used in the mathematical model equation.Standard procedure could be applied to the continuous system to derive Laplace Transfer Function [22], [23]

C. Proportional-Integral Control of Continuous Milk Cooling
The system and block diagram for closed loop systems are shown in Figures 3 and 4. For simplicity, the Ziegler-Nichol tuning method would be used.Ziegler Nichol formula for Proportional-Integral controller are as follows [19], [30], [31], [32].The ultimate period, Pu, is defined as the period of sustained cycling that would occur if a proportional controller with gain Ku were used.Ku and Pu were easily determined from the Bode diagram [34], [35] ωco is cross-over frequency [36].

A. Dynamic Study of Batch Process
The milk cooling process in Cooperatives (Koperasi Unit Desa) is generally carried out as a batch process.The batch process only depends on time.When the milk temperature has reached the expected, then the milk is released.Assuming that the overall heat transfer coefficient (U) is constant, the result of the mathematical model is shown in Equation 14.
From the mathematical model of the batch system above, the simulation results are shown in Figure 5. Based on Figure 5, it can be seen that the batch system only affects cooling time.The milk temperature comes in at 36˚C and the temperature reaches 4˚C with a cooling time of 120 minutes.Since the gain is 1 while the time constant is 42.3548 minutes, the process will reach 2 ̊ C at about more than 4 times the time constant [21].The time constant is the time required for the first-order process to reach 63.2% to the new steady state value when a step disturbance is given to the input variable [19], [37].Now, suppose that cooling water temperature varies as shown in Figure 6.The milk temperature profile in the tank can be investigated and shown in Figure 7. Various results for different cooling water temperatures can also be obtained easily using Simulink.

B. Dynamic Study of Continuous Process
The mathematical model in Equation 17was simulated in Simulink to obtain the input-output data.The input and output results were then identified using the System Identification Toolbox (SIT) in Matlab.Suppose we did not derive Equation ( 17), but using the input-output data collected directly from the plant (or, in this study, simulated using Equation 17), we still can identify the system transfer function [38].

C. Proportional-Integral Control of Continuous Process
Proportional Integral controller settings using the Ziegler-Nichol tuning method are shown in Table 3.Here we assumed several values of time delay since the process itself is a firstorder process without time delay, and Cross-over frequency cannot be determined for a first-order process without time delay.Therefore, a small value of time delay should be added to the process so that we can have cross-over frequency, ultimate gain and ultimate period.
The Bode Diagram was used in this method to calculate Kc and τI.The following is an example of a calculation assuming Time Delay = 1 min.
From the Bode diagram (Figure 11), the cross-over frequency value was 3.283 rad/time, and the magnitude was 0.00003061 abs, so that the ultimate gain (Ku) and ultimate period (Pu) could be determined:   Figure 12,13,14 shows the output profile (To) and its set point for using the assumption of time delay from 1, 2, and 3 mins, respectively.It can be seen that the higher the time delay assumption, the greater the error will be.

IV. CONCLUSION
The dynamic model of the batch milk cooling process has been derived and simulated using Simulink.The overall heat transfer coefficient (U) was assumed constant and could be calculated as the batch process finished in 2 hours.The milk temperature profile shows reasonable dynamics during the process.The dynamic model of continuous milk cooling has also been presented.Suppose we did not derive Equation ( 17), but we can still identify the system transfer function g the input-output data collected directly from the plant (or, in this study, simulated using Equation 17).Several methods for system identification are provided in MATLAB.These include transfer function models, state space models, polynomial models (i.e., ARX, ARMAX, Output-Error, Box Jenkins) etc.Those models can be converted to Laplace Transfer Function easily.The system identification can be applied for the continuous system and the batch milk cooling process.The mathematical model is particularly useful for control purposes.Here we demonstrated Proportional-Integral control of the continuous milk cooling process, i.e., using the Ziegler-Nichol tuning method for determining controller parameters.Also, we introduced a small-time delay to the system to simplify the tuning process.Other tuning methods are also available such as Tyreus-Luybean, Hagglund-Astrom, Ciancone, Direct Synthesis etc. [43], [44], [45].Further research may include system identification for different sizes of equipment or different type of cooling processes, and advanced control strategies for the milk cooling process.

Fig. 2
Fig. 2 Stirred milk cooling tank system seen for the side (courtesy of koperasi sae pujon)

Fig. 3 Fig. 4
Fig. 3 Continuous milk cooling system controlled by PI controller

Fig. 6
Fig. 6 Temperature profile of inlet chilled water

Fig. 8
Fig. 8 Comparison of output between original model and SIT with Ti perturbation

Fig. 12
Fig. 12 Output profile and setpoint for 1 min time delay

Figure 15
Figure15and 16 show profile of output for regulatory problem (disturbance rejection) for controller setting using 4and 5-mins assumption.Again, it can be seen that the smaller the time delay assumption the better.

Fig. 15
Fig. 15 Output profile (To) and its setpoint for 4 min time delay with changes in temperature of inlet milk (Ti)

TABLE I BATCH
AND PARAMETER DATA

TABLE III PI
CONTROLLER SETTING BASED ON THE ZIEGLER-NICHOL METHOD