Fed-Batch Sucrose Crystallization Model for the B Massecuite Vacuum Pan, Solution by Deterministic and Heuristic Methods

: Fed-batch crystallization is a crucial step for sugar production. In order to relate parameters that are di ﬃ cult to measure (average diameter of the crystals and total mass formed) to other easier to measure parameters (volume, temperature, and concentration), a model was developed for a B massecuite vacuum pan composed of mass and energy balances together with empirical relations that describe the crystal development inside equipment. The generated system of ordinary di ﬀ erential equations (ODE) had eight parameters which were adjusted through minimization of relative di ﬀ erences between the model results and experimental data. It was solved through the function fmincon, available in MATLAB TM , which is a deterministic and gradient-based optimization method. The objective of this paper is to improve the model obtained and, for this purpose, two metaheuristic functions were used: genetic algorithm and particle swarm. To compare the results, the convergence time of each algorithm was used as well as the resulting quadratic deviation. The particle swarm method was the best option among the three used, presenting a shorter execution time and lower quadratic relative deviation. Author Contributions: Conceptualization, P.E.d.M.G.; methodology, P.E.d.M.G.; software, P.E.d.M.G.; validation, P.E.d.M.G., M.A.d.S.P.J. and J.E.O.; formal analysis, P.E.d.M.G.; investigation, P.E.d.M.G. and M.A.d.S.P.J.; resources, P.E.d.M.G., M.A.d.S.P.J. and J.E.O.; curation, P.E.d.M.G.; writing—original draft and visualization,


Introduction
Crystallization is the main process in the sugar industry for the separation and obtaining of sucrose in its commercial form, sugar, and is one of the crucial steps to increase the yield from the plant [1]. Among all steps of sugar production, fed-batch crystallization is the most complex in operational terms since it requires greater care by the plant operators [2].
The industrial sucrose crystallization process consists of three basic equipment: vacuum pan, crystallizer, and centrifuge. In the vacuum pan, sugar crystals form and grow, at a constant temperature, in a fed-batch process. In crystallizers, the crystals grow due to the cooling of the massecuite, which increases supersaturation. The centrifuges separate the crystals formed from the mother liquor which, once exhausted, is called molasses. In this sequence, about half of the sucrose is recovered as sugar crystal. To maximize this recovery, crystallization schemes are used, inserting new sequences of this equipment that are fed by the molasses from the previous step. These sequences can be repeated once (two massecuites scheme), or twice (three massecuites scheme). It is interesting to note that because each sequence of these equipment recovers around half of the crystals, a two massecuite scheme recovers 75% of the sucrose as crystal, while a three massecuite scheme recovers around 87.5%. Because the cost of establishing a three massecuites plant is 50% higher with a mere 12.5% increase in sucrose recovery and most Brazilian industries have an ethanol plant attached that consumes molasses, a mere 12.5% increase in sucrose recovery and most Brazilian industries have an ethanol plant attached that consumes molasses, the two massecuites scheme is the most widespread in the country [2,3]. Figure 1 shows the flowchart of the two massecuites scheme. Figure 1. Two massecuites scheme flowchart [4], based on [2].
As Brazil is among the largest sugar producers in the world, the scheme of the two massecuites shown in Figure 1 was the one used in the experiments [3]. As local plants also produce electricity with the steam from the boilers, the massecuite cooling in the crystallizers is not used, which eliminates the need to reheat the massecuite to feed the centrifuges (cooling increases the viscosity of the mother liquor) [3,4]. Thus, the crystallization of sucrose occurs only in vacuum pans, making its study the focus of this work.
The instrumentation of massecuite vacuum pans is a challenge for automation companies due to the high values of the sensors and large amount of equipment to be monitored [5,6]. In order to reduce these costs, two models were developed to relate the mean diameter of the crystals D (4,3) and the total formed crystal mass (FCM)-the data of which are difficult to measure-with the equipment concentration, volume, and temperature; one for A massecuite vacuum pan and the other for B massecuite vacuum pan [7]. Both models were solved through the function fmincon, available in MATLAB™, which is a deterministic and gradient-based optimization method [8,9].
With these models, it is possible to develop soft sensors to measure crystal size and quantity in the vacuum pan. It is essential that the program is fast and accurate. Thus, to improve the model obtained for B massecuite vacuum pan, which has eight parameters for adjustment (seven kinetic parameters and a thermodynamic data), two metaheuristic methods were used-genetic algorithm (GA) and particle swarm. To compare the results, the convergence time of each algorithm was used as well as the relative resulting quadratic deviation.

Methodology
The equipment used to obtain the experimental data is installed in a São Paulo state sugarcane industry. During the experiment, the volume, concentration, and temperature data of the fluid were recorded by the sensors installed in the equipment. Each batch occurred in two hours, and samples were taken every fifteen minutes through the equipment probe to evaluate both FCM and D (4,3) As Brazil is among the largest sugar producers in the world, the scheme of the two massecuites shown in Figure 1 was the one used in the experiments [3]. As local plants also produce electricity with the steam from the boilers, the massecuite cooling in the crystallizers is not used, which eliminates the need to reheat the massecuite to feed the centrifuges (cooling increases the viscosity of the mother liquor) [3,4]. Thus, the crystallization of sucrose occurs only in vacuum pans, making its study the focus of this work.
The instrumentation of massecuite vacuum pans is a challenge for automation companies due to the high values of the sensors and large amount of equipment to be monitored [5,6]. In order to reduce these costs, two models were developed to relate the mean diameter of the crystals D (4,3) and the total formed crystal mass (FCM)-the data of which are difficult to measure-with the equipment concentration, volume, and temperature; one for A massecuite vacuum pan and the other for B massecuite vacuum pan [7]. Both models were solved through the function fmincon, available in MATLAB TM , which is a deterministic and gradient-based optimization method [8,9].
With these models, it is possible to develop soft sensors to measure crystal size and quantity in the vacuum pan. It is essential that the program is fast and accurate. Thus, to improve the model obtained for B massecuite vacuum pan, which has eight parameters for adjustment (seven kinetic parameters and a thermodynamic data), two metaheuristic methods were used-genetic algorithm (GA) and particle swarm. To compare the results, the convergence time of each algorithm was used as well as the relative resulting quadratic deviation.

Methodology
The equipment used to obtain the experimental data is installed in a São Paulo state sugarcane industry. During the experiment, the volume, concentration, and temperature data of the fluid were recorded by the sensors installed in the equipment. Each batch occurred in two hours, and samples were taken every fifteen minutes through the equipment probe to evaluate both FCM and D (4,3) through a Nikon Eclipse E200-LED Binocular microscope. Figure 2 shows how the crystals look at the end of the crystallization process in B massecuite vacuum pan and Figure 3 shows the distribution curve of the crystal size at the last point of the first experiment.
Processes 2020, 8, x FOR PEER REVIEW 3 of 12 through a Nikon Eclipse E200-LED Binocular microscope. Figure 2 shows how the crystals look at the end of the crystallization process in B massecuite vacuum pan and Figure 3 shows the distribution curve of the crystal size at the last point of the first experiment.  The mathematical modeling was adapted to a model of a process in pilot scale-batch crystallization with cooling (using evaporation and cooling the solution)-as developed in [10]. The pilot scale process starts with approximately 752 cm 3 of syrup, operating for 40 min at 60 °C and is then cooled for 50 min to 40 °C. During its operation, the massecuite volume decreases since it is not fed with syrup and due to evaporation of part of its contents. The B massecuite vacuum pan operation starts with one-third of the equipment volume (10 m 3 ) filled with syrup, occurring for two hours at a constant temperature of 80 °C, receiving a variable flow of A molasses so that the volume of equipment reaches its end capacity (30 m 3 ) in compensating for evaporation of the solvent. The differences between the two processes are summarized in Table 1. through a Nikon Eclipse E200-LED Binocular microscope. Figure 2 shows how the crystals look at the end of the crystallization process in B massecuite vacuum pan and Figure 3 shows the distribution curve of the crystal size at the last point of the first experiment.  The mathematical modeling was adapted to a model of a process in pilot scale-batch crystallization with cooling (using evaporation and cooling the solution)-as developed in [10]. The pilot scale process starts with approximately 752 cm 3 of syrup, operating for 40 min at 60 °C and is then cooled for 50 min to 40 °C. During its operation, the massecuite volume decreases since it is not fed with syrup and due to evaporation of part of its contents. The B massecuite vacuum pan operation starts with one-third of the equipment volume (10 m 3 ) filled with syrup, occurring for two hours at a constant temperature of 80 °C, receiving a variable flow of A molasses so that the volume of equipment reaches its end capacity (30 m 3 ) in compensating for evaporation of the solvent. The differences between the two processes are summarized in Table 1. The mathematical modeling was adapted to a model of a process in pilot scale-batch crystallization with cooling (using evaporation and cooling the solution)-as developed in [10]. The pilot scale process starts with approximately 752 cm 3 of syrup, operating for 40 min at 60 • C and is then cooled for 50 min to 40 • C. During its operation, the massecuite volume decreases since it is not fed with syrup and due to evaporation of part of its contents. The B massecuite vacuum pan operation starts with one-third of the equipment volume (10 m 3 ) filled with syrup, occurring for two hours at a constant temperature of 80 • C, receiving a variable flow of A molasses so that the volume of equipment reaches its end capacity (30 m 3 ) in compensating for evaporation of the solvent. The differences between the two processes are summarized in Table 1. The population balance of crystals of the model used in this paper is calculated using the following partial differential equation (PDE) [7,10,11]. This PDE is shown in Equation (1): with initial condition n(L, t) = n 0 (L, t) t = 0 and the boundary condition α(L) is known as production-reduction term and it groups the rates of birth and destruction of crystals. Applying the method of moments (MOM) in Equation (1), the following system of ordinary differential equations (ODEs), presented in Equations (3) and (4), is obtained [10]: where V is the mass volume inside the equipment in cm 3 . The number of moments necessary for the model resolution can tend to infinity, but due to the use of the empirical equation shown below (Equation (5)), the moment µ 3 enters the mass balance, causing the four moments (µ 0 , µ 1 , µ 2 , and µ 3 ) to be enough to solve the ODE system [12]. B 0 (N • particles/cm 3 min) and G (cm/min) are the nucleation and growth rates, respectively, that can be modeled by the empirical power equations presented in Equations (5) and (6) [10]: N r is the stirrer speed in rpm. The terms K b (N • part/cm 3 ·min·(g/cm 3 ) j ·(rpm) p ), b, j, p, K g (cm/min·(rpm) q ), g, and q are the kinetic parameters of these equations and were estimated and optimized using the experimental data.
In Equation (5), M T represents the total mass of crystals and can be obtained using Equation (7): where ρ c is the density of the sucrose in g/cm 3 , and K v = π/6 is the characteristic form factor for sugar crystals [10,13]. In Equations (5) and (6), S r represents the relative supersaturation, and is defined by Equation (8): where C is the sucrose concentration and C sat is the saturation concentration point, both in g/cm 3 . The mass percentage of sucrose in the saturated solution (Brix sat ) can be calculated as a function of temperature using Equation (9) [4,7,14].
where T is the temperature in • C. The density of a sucrose solution (ρ) can be obtained by Equation (10) (T in • C) [7,15]: The solution concentration, as well as its saturation concentration, can be obtained by Equation (11): The mass balance of the process is presented in Equation (12) [7]: where m sucrose and m sucrose 0 is the total mass of sucrose (g) in the equipment at time t and 0, respectively, C mA is the concentration of A-molasses (g/cm 3 ), and V mA is the accumulated flow of A-molasses added in equipment (cm 3 ). Equation (13) shows that the total mass of sucrose inside the equipment is partly in solution (CV S ) and partly crystallized (ρ c V c ): Equation (14) shows that the volume of the solution (V S ) plus the volume of the crystals (V c ) is equivalent to the volume of massecuite in the equipment (V). Rearranging: Thus, replacing Equations (13) and (14) in Equation (12): deriving Equation (16) the Equation (17) is obtained: Substituting V c for K v µ 3 (t) and isolating dC/dt leads to Equation (18): with initial condition The energy balance of B massecuite vacuum pan is shown in Equation (20) [2]: where E(t) is the internal energy of the system in J, F in and F out are the mass input and output in g, h in and h out the relative enthalpies in J/g. Replacing the Equation (20) variables with those of the B vacuum pan [7,10]: The subscripts mA and VC represent A molasses and mass in vacuum pan, respectively. Deriving Equation (21): The heat rate can be obtained by Equation (23) [2]: T j is the temperature of the steam in the equipment's calandria ( • C), U l is the global coefficient of thermal exchange (J/( • C·min·cm 2 )), and A l is the heat exchanger area (cm 2 ). The work rate is calculated by Equation (24) [11]: where P is the absolute pressure in J/cm 3 . Replacing Equations (23) and (24) in Equation (22) and explaining the volumes of A molasses and crystals formed, Equation (25) is obtained: with initial condition ∆H c is the crystallization heat in J/g, which can be calculated by means of Equation (27) (T in • C) [10,13]: Also in Equation (25), m ev is the evaporated accumulated mass flow rate and H ev is the vaporization enthalpy in J/g, which can be obtained by Equation (28) (T in • C) [16]: The absolute pressure P, in J/cm 3 , within the equipment can be obtained using the Antoine equation presented in Equation (29) (T in • C) [16]: P = 2.21 × 10 1 e (6.53247− 7173.79 1.8T+421.4747 ) .
The specific heat of sucrose solutions, Cp in J/g • C, can be obtained by means of Equation (30) (T in • C) [7,15,17]: The dynamics of the volume of liquid inside the equipment, V in cm 3 , is described by Equation (31), which is a cubic correlation obtained through the equipment data during the experiment (t in min): V(t) = −1.5896t 3 − 75.299t 2 + 1.9859 × 10 5 t + 10 7 .
(31) Table 2 presents the parameters used in this modeling: Table 2. Process parameters to solve the population, mass, and energy balances.

Parameter Value Unit
In Equations (5) and (6), the proposed mathematical model includes seven kinetic parameters (K b , b, j, p, K g , g, and q) specific for each batch crystallization operating condition [10]. These values need to be adjusted so that the model provides representative results of this process. The global coefficient of thermal exchange (U I ) was also added in the adjustment [7]. Based on this, it sought to reduce the relative quadratic differences between the experimental data and those obtained by the model, according to Equation (32): In the model deterministic solution, the non-linear optimization tool fmincon, available in MATLAB TM , was used in conjunction with the use of the ode23s [7,10]. In this paper Equation (32) was also solved using two metaheuristic tools, one based on the genetic algorithm (ga) and the other on the particle swarm (particleswarm) [18,19]. The results were obtained with the hardware configuration shown in Table 3: Table 3. Hardware configuration used in the modeling.

CPU
i7 4770k RAM 2 × DDR3 1600 MHz-08 GB HD 240 GB SSD SATA 3 The algorithms were executed in order to record the convergence time and the relative quadratic deviation of the obtained results with the experimental data.  Table 4 presents the operational data from experiments 1 and 2: Table 4. Operational data from experiments 1 and 2 (V mA is accumulative). In these conditions, the size of formed crystals and their quantity were evaluated, shown in Table 5: With known data from Tables 2, 4 and 5, the ODE system was solved using the two meta-heuristic methods to reduce the differences between the results of the model obtained in relation to the two experiments (Equation (32)), thereby, obtaining the unknown parameters-the seven kinetic parameters and global coefficient of thermal change (U l ) of system-presented in Tables 6 and 7 along   Table 8 presents the average execution time of each algorithm and the relative quadratic deviation obtained at the end of the execution of the algorithm. It is noted that the particle swarm method was able to obtain a smaller deviation with a shorter execution time compared to the other methods. The genetic algorithm also performed better than the deterministic method used by the fmincon function. It is noteworthy that some point results of the deterministic solution were able to converge in less time than the average of the other methods, but without approaching the deviations obtained by them. This occurs due to the high dependence of the deterministic methods of the initial point, which in most cases leads them to local minimums [20].

Experiment 1 Experiment 2 t (min)
In the graphs shown in Figures 4-6, it is possible to verify that all the algorithms made a good adjustment to the data.  It is noted that the particle swarm method was able to obtain a smaller deviation with a shorter execution time compared to the other methods. The genetic algorithm also performed better than the deterministic method used by the fmincon function. It is noteworthy that some point results of the deterministic solution were able to converge in less time than the average of the other methods, but without approaching the deviations obtained by them. This occurs due to the high dependence of the deterministic methods of the initial point, which in most cases leads them to local minimums [20].
In the graphs shown in Figures 4-6, it is possible to verify that all the algorithms made a good adjustment to the data.    The correlation coefficient (R 2 ) values were obtained through multiple linear regressions with both experiments, and followed the behavior of the relative square deviations shown in Table 8. The graphs obtained for the metaheuristic methods were similar to those obtained for the deterministic method [7].  The correlation coefficient (R 2 ) values were obtained through multiple linear regressions with both experiments, and followed the behavior of the relative square deviations shown in Table 8. The graphs obtained for the metaheuristic methods were similar to those obtained for the deterministic method [7]. The correlation coefficient (R 2 ) values were obtained through multiple linear regressions with both experiments, and followed the behavior of the relative square deviations shown in Table 8. The graphs obtained for the metaheuristic methods were similar to those obtained for the deterministic method [7].

Conclusions
This paper has tried to refine a model obtained for B massecuite vacuum pan by implementing metaheuristic methods in the optimization of the kinetic parameters and global coefficient of thermal change. As shown in Table 8 and Figures 4-6, the particle swarm method was the best option among the three used, presenting a shorter execution time, a lower quadratic relative deviation, and higher correlation coefficients (R 2 ). Nevertheless, when analyzing the data fit, it is noted that all solutions were very close, making the execution time the determining variable for choosing the algorithm to be used, keeping in mind that the development of a soft sensor to measure the size of crystals and their quantity in the vacuum pan is the target of this work. The solution through the deterministic method was highly influenced by the initial point, as expected. As a next step to this work, the refinement of the model of the operation for A massecuite vacuum pan can be carried out, as well as the use of other functions for comparison.