Nonlinear MIMOControl of a Continuous Cooling Crystallizer

In this work, a feedback control algorithm was developed based on geometric control theory. A nonisothermal seeded continuous crystallizer model was used to test the algorithm. The control objectives were the stabilization of the third moment of the crystal size distribution (μ3) and the crystallizer temperature (T); the manipulated variables were the stirring rate and the coolant flow rate. The nonlinear control (NLC) was tested at operating conditions established within the metastable zone. Step changes of magnitudes ±0.0015 and ±0.5◦C were introduced into the set point values of the third moment and crystallizer temperature, respectively. In addition, a step change of ±1◦C was introduced as a disturbance in the feeding temperature. Closed-loop stability was analyzed by calculating the eigenvalues of the internal dynamics. The system presented a stable dynamic behavior when the operation conditions maintain the crystallizer concentration within the metastable zone. Closed-loop simulations with the NLC were compared with simulations that used a classic PID controller. The PID controllers were tuned by minimizing the integral of the absolute value of the error (IAE) criterion. The results showed that the NLC provided a suitable option for continuous crystallization control. For all analyzed cases, the IAEs obtained with NLC were smaller than those obtained with the PID controller.


Introduction
Crystallization is an operation of great importance in industry due to the wide variety of materials that are sold in crystal form.Well-established examples include the crystallization of salts, sugar, detergents, fertilizers, foodstuffs, pigments, and pharmaceuticals.Its wide use is due to several reasons: a crystal formed from an impure solution is essentially pure, it requires little energy for the separation, and it provides a practical method for obtaining the product in a suitable condition for packaging and storage.Note, however, that all these benefits are possible when a good control strategy is followed and all control objectives are attained.At industrial level, the objective of crystallization is to obtain a product with the specified characteristics of grain size and purity.
The recent efforts on model-based control of crystallization process have been motivated by significant advances in the modeling of crystallization processes.Population balances have provided a framework for the mathematical modeling of crystallization processes [1].Dynamics and stability of isothermal continuous crystallization have received much attention in the literature [2][3][4][5][6].However, the analysis of nonisothermal crystallization systems has not been widely explored because of the difficulty of having accurate measurements of the metastable zone width and the crystallization kinetics.Recent advances on crystallization modeling and new developments in the field of online measurement technology have motivated the implementation of new control strategies such as model-based control [7].
In the chemical industry, most of the processes are controlled by means of PID control strategies.However, not always this type of controllers achieves the desired objectives for instance when the process has an intrinsic unstable behavior with time-varying parameters [8].Continuous crystallization is a highly nonlinear process mainly because it involves two nonlinear distinct kinetic steps: nucleation and growth.In addition, there are interactions between kinetic, fluid dynamics, and crystal size distribution.All those factors might lead to an unstable operation (cycling behavior [9]) that affects the quality of the final product.Therefore, A control strategy that warranties a stable operation is necessary.
NLC [10,11] has become an important field in model-based control techniques.Input/output linearization approach [12][13][14] has proven being popular.It has been applied successfully to several processes such as reactors [15,16], fermenters [17][18][19], distillation [19], continuous isothermal crystallization [7,[20][21][22][23], batch crystallization [15,24], vaccination controls in epidemic models [25], and power systems [26].Significant improvements in the mean size of the final product have been obtained with the application of NLC to crystallization process.Other modelbased control techniques, such as model predictive control, have been also applied successfully to continuous isothermal crystallization process [27,28].Therefore, it is expected that the complex dynamic behavior of the nonisothermal continuous crystallization process might be addressed in a better way by using new control schemes that allow an efficient, stable, and safe operation.
In this paper, we develop a nonlinear feedback controller by using geometric control theory to be applied to a nonisothermal crystallizer.Operating conditions are analyzed based on the combined effects of secondary nucleation and cooling temperature.The algorithm includes two control loops.The first one relates the crystallizer temperature with the coolant flow rate and the second one the third moment of the CDS with the stirring rate.Closed-loop stability is analyzed by calculating the eigenvalues of the internal dynamics.Finally, closed-loop simulations with the NLC are compared with simulations that used a traditional PID controller.The operating ranges are the same as those used for calculating the kinetic equations.These operating values impose additional constraints to the crystallizer controllers.

Nonlinear Model of the Continuous Crystallizer
The overall model for a nonisothermal continuous crystallizer includes mass, energy, and population balances.The mass balance gives the solute concentration in the continuous phase and accounts for flow of solute into and out of the system and the mass transfer to the solid phase via nucleation and growth of crystals.Equation (1) shows the general mass balance where M s is the solute mass solved in the solvent, V L is the suspension solids-free volume, V is the total suspension volume, C is the solute concentration in the continuous phase, ε is the fraction of the volume free of solids, q is the volumetric flow into and out of the system, and R is the global mass transfer of solute.The energy balance accounts for enthalpy differences between the flow streams, the heat of crystallization and heat removal by cooling systems.The following equation represents the energy balance for an incompressible system and constant volume of continuous phase ρV c p dT dt where ρ is the slurry density, c p is the specific heat capacity, T is the crystallizer temperature, H k is the enthalpy of the kth stream, ΔH c is the heat of crystallization, and H ext is the net heat removal.The population balance extended the concept of particle size distribution to include an arbitrary number of internal coordinates necessary to describe the state of the particles (crystal population density, n(L)).Randolph and Larson [29] presented a general crystallization model developed under the following assumptions: continuous operation, perfect mixing, infinitesimal new generated particles' size, and constant volume.The following population balance describes the evolution of n(L) as function of time and crystal size: in which n is the population density, G is the crystal growth rate, L is the internal coordinate (characteristic crystal length), B is the nucleation rate density, and δ(L − L 0 ) is the Dirac delta function acting at L = L 0 .The boundary and initial conditions for (3) are given as follows: In order to complete the modeling of the crystallization process, constitutive relationships are required for growth and nucleation.Salcedo-Estrada [30] reported growth and nucleation equations (5) for the ammonium sulfate-water system.The equations included the following parameters k g , k b , g, h, o and p that were evaluated from experimental data obtained at operating condition located in the metastable zone.Both equations are reported as function of supersaturation and agitation rate ( The relative supersaturation S r is evaluated as (C − C sat )/C sat , N r represents the agitation rate, and M T is the total crystal mass per volume.Table 1 shows the kinetic parameters included in (5).Different approaches to solve the model presented by (1)-( 3) have been shown in the literature.Ramkrishna [31] presented some numerical strategies for accurate computation of the solution.In this work, the population balance is reduced by a process based on the method of moments.Equation ( 3) is transformed into ( 6) and (7) assuming that the dominant dynamic of the system may be represented by a small number of degrees of freedom [20].Equations ( 8) and ( 9) represent the crystallizer mass and energy balances, and (10) corresponds to the jacket energy balance The mean crystal size and standard deviation can be evaluated using the following relationships of the moments in order to obtain the complete CSD: As mentioned previously, this work was developed for a crystallizer operating inside the metastable zone.Therefore, it is important to know the limits of the metastable zone.Lugo-Martinez [32] reported a polynomial function for the upper metastable zone limit, c mz , for the ammonium sulfate-water system as function of temperature.The lower limit (saturation concentration, c s ) was evaluated using the polynomial reported by Perry et al. [33].Both ( 12) and ( 13) are valid in the range from 20 to 50

Geometric Properties
Nonlinear geometric control has been widely analyzed in the literature, and the theoretical basis may be reviewed at Isidori [10] and Slotine and Li [34] among others.The mathematical model for a MIMO system that contains the state variable and controlled output vectors is shown in the following: The f (x) and g i (x) are smooth vector fields, h 1 (x), . . ., h p (x) are smooth scalar fields, u i are the manipulated variables, ẋ is state vector, and y i are the output variables.The mathematical model for the crystallization process may be obtained making the following changes into ( 6)-( 10), μ i = x i , C = x 6 ,T = x 7 ,T c = x 8 , and reordering the vector fields as follows: where x ∈ R 9 and u ∈ R 2 .

Relative Degree.
The relative degree of a multivariable nonlinear system may be calculated with the Lie derivatives of h i (x) along the vector fields f (x) and g i (x).It is defined as the smaller integer r i such that L gi L ri−1 f h(x) / = 0 for all 1 ≤ i ≤ m in a neighborhood of x 0 and an m×m nonsingular matrix defined as Considering h 1 (x) = x 3 and h 2 (x) = x 7 as the control objectives, u 1 = N r and u 2 = q w as the manipulated variables, and assuming a control problem with direct relationship between u 1 → x 3 and u 2 → x 7 , the relative degree vector [r 1 , r 2 ] is calculated from L g1 L r1−1 f h 1 (x) and L g2 L r2−1 f h 2 (x) different from zero.The evaluation gives relative degrees r 1 = 1 and r 2 = 2.The system overall relative degree is three assuming that (T ce − x 8 ) and x 2 are different from zero.The first condition is fulfilled at normal cooling operation and the second in a seeded crystallization process.Matrix C(x) can be constructed with the values given in ( 17)- (20).The evaluation of the determinant of C(x) gives a value different from zero (nonsingular) when T ce / = x 8 .This fact ensures the existence of both the inverse of C(x) and the control law 3.2.Exact Linearization via Feedback.A nonlinear system may be transformed into a linear and controllable system by means of feedback and change of coordinates in the state space.In this work, the overall relative degree is less than the number of state variables, so it is possible to find n − r m additional functions φ n−rm to complete the map φ(x) in such a way that the Jacobian of φ(x) is nonsingular at x 0 .The proposed functions to complete the map are shown in (22).
It should be noted that these additional functions were set arbitrarily and satisfy the nonsingular Jacobian constrain, det(J φ(x) ) = −UA/ρV c p / = 0 Finally, the complete map can be represented as follows: The system given by ( 14) may be represented in its new coordinates set as follows: where z 1−3 ∈ R 3 represent the system controllable and observable parts, while z 4−9 ∈ R 6 represent the nonobservable and uncontrollable parts of the system (internal dynamic) (Slotine and Li [34]).

Zero Dynamics.
The zero dynamics of the system described by ( 24) is obtained when the output errors are equal to zero , for example, x 3 −x * 3 = z 1 = 0, x 7 −x * 7 = z 2 = 0 and ż1 = ż2 = ż3 = 0.The system may be represented as follows: q w ( c m 3 /m in ) Nr (r pm ) If the system given by ( 25) is linearized using an equilibrium set of operation conditions that guarantee a stable performance, the following linear equation system is obtained: ż = Ax, where A represents the Jacobian matrix of the system (25).
In order to find the equilibrium set, steady-state crystallizer concentrations were evaluated for all possible combination of N r and q w and plotted together with the upper (12) and lower limits (13) of the metastable zone.It was found that combination of values of N r smaller than 300 rpm and cooling flow rates bigger than 2200 cm 3 /min produces crystallizer concentrations outside the metastable zone (see Figure 1).Therefore, it was selected a combination of N r = 420 rpm and q w = 800 as the steady state starting operation conditions that ensure the equilibrium operation with a dynamic closed-loop response asymptotically stable.The eigenvalues of A at the selected conditions were −0.3846, −0.3992, −0.5679 + 0.3252i, and −0.5679 − 0.3252i.It is clear that the real part of all eigenvalues is negative; therefore, the system presents a local stability in its internal dynamics.Table 2 shows the steady-state values calculated at cooling flow rate of 800 cm 3 /min and agitation rate of 420 rpm.These values were used to generate the nonsingular matrix A shown in the following: 3.4.MIMO Control Law.If the system relative order is well defined and the internal dynamics is asymptotically stable, then the lineal input/output control law will make the third moment of the CSD and the crystallizer temperature to converge to the reference values x * 3 and x * 7 under normal operating conditions.The proposed pairing u 1 → x 3 and u 2 → x 7 do not consider interactions between the other input and output variables; for example, input variable v i modifies only the output variable y i .The control law may be expressed as follows: The control actions may be calculated solving the system given by (28), where C represents the uncoupled system matrix given previously in ( 16): Using a proportional control action only, v 1 = −k 0 (x 3 −x 3sp ), (29) becomes The cooling flow rate, q w , might be evaluated as follows: where

Results
The crystallizer mathematical model and the control laws given by ( 30) and ( 31) were used to build a Matlab-Simulink algorithm (NLC).In addition, a second algorithm with a PID controller was developed in order to compare the dynamic responses.Table 3 includes the physical parameters of the model used in this work.The manipulated variables were limited to the following intervals: N r ∈ [300 800] rpm and q w ∈ [100 8000] cm 3 .These values were selected from physical operating conditions that maintain the crystallizer concentration within the metastable zone.The lower stirring limit guarantees a homogeneous solution inside the crystallizer, and the upper limit represents the limit before secondary nucleation being considerable due to crystal breakage.The cooling flow rate limits were the same used by Salcedo-Estrada [30] at the experiments where the kinetics parameters were evaluated.Simulations were performed using a seeded inlet stream with crystals of the same size L 43 = 129 μm at T e of 30 • C and a flow rate of 1153.84 cm 3 /min.The initial conditions for the state variables used are shown on Table 2.An initial equilibrium operation was assumed at N r = 420 rpm and q w = 800 cm 3 .The distribution moments were evaluated at these conditions and used in all simulations, and their values are shown on Table 3.The parameters values for the NLC and the PID controllers were obtained minimizing the integral absolute error (IAE) using the routine lsqnonlin of SIMULINK.The NLC parameters are k  has a large effect on the distribution moments and a small effect on the crystallizer temperature.Figure 3 illustrates the metastable zone limit (C mz ), the solubility profile (C s ), and the solution concentration within the crystallizer as function of time and temperature.It is important to observe that concentration remains within the metastable zone regardless the introduction of step changes in μ 3 .
Run 2. The control algorithms were tested introducing step changes in the crystallizer output temperature.Two changes (+0.5 • C and −0.5 • C) were introduced at times 20 and 60 minutes.The set point value for the third moment was kept constant at 0.004.Figure 4 shows the dynamic response of the output and manipulated variables.The step changes in the crystallizer temperature have no effect on the third moment.Final crystallizer temperature is reached faster when the NLC algorithm is used compared to the PID algorithm.The manipulated variables present smooth profiles with larger overshoots for the NLC algorithm.As it is expected, there is a relationship between the crystallizer temperature and agitation rate through (9) and (10).For both controllers, the system reaches the set point values, but the IAE, for μ 3 and T are smaller when the NLC algorithm is implemented.Again, the solution concentrations for these runs are kept inside the metastable zone as shown in Figure 5(b).Run 3. Finally, the control laws were tested introducing disturbances in the feed temperature.
Step changes of magnitude ±1 • C were used.The simulation results are shown in Figure 6.Notice how the cooling flow rate increases considerably when the feed temperature is one degree higher.
The controller is able to keep the system at the desired set point values, and the crystallization process is performed within the metastable zone.Table 4 shows the IAE generated during the full run.For both controllers, the system reaches the set point values, but, again, the IAEs for μ 3 and T are smaller when the NLC algorithm is implemented.
In addition, it is important to point out that the NLC algorithm uses three control parameters only while the two-PID algorithm add six parameters.Therefore, it is easier to tune the three control parameters of the NLC algorithm.

Conclusions
A nonlinear geometric feedback control algorithm was developed and applied for the control of a continuous MIMO crystallizer.Two feedback control loops were established.The first loop controls the crystallizer temperature manipulating the cooling flow rate, and the second loop controls the third distribution moment manipulating the agitation rate.The stability of the closed-loop system was analyzed using the eigenvalues of the internal dynamics showing that if the solution concentration is maintained inside the metastable zone, the system presents an asymptotically stable behavior.
Set point tracking and responses to feed load disturbances in temperature were carried out.In the three runs analyzed, the IAE values generated when NLC strategies were smaller than those generated with traditional PID control strategies.Finally, it is easier to tune the three control parameters of the NLC strategy than the six parameters of the two PID controllers, and it has been shown in this work that NLC provides an excellent potential for the control of nonisothermal crystallization.

Figure 1 :
Figure 1: Steady-state concentration and metastable zone limits at different operating conditions.

Run 1 .Figure 2 :
Figure 2: Closed-loop dynamic responses for the MIMO system when step changes in μ 3 are introduced.

Figure 3 :Figure 4 :Figure 5 :Figure 6 :
Figure 3: Concentration profiles inside the crystallizer when step changes are imposed to the output variable μ 3 .

Table 2 :
Equilibrium state variables calculated at operation conditions q w = 800 cm 3 /min and N r = 420 rpm.

Table 3 :
Operating conditions and physical constants used in the model.

Table 4 :
IAE values generated when NLC and PID control strategies were used.