Nonlinear MIMO Control of a Continuous Cooling Crystallizer

Abstract

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 (𝜇₃) and the crystallizer temperature (𝑇); 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.

1. 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–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–14] has proven being popular. It has been applied successfully to several processes such as reactors [15, 16], fermenters [17–19], distillation [19], continuous isothermal crystallization [7, 20–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 model-based 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.

2. 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 𝑑𝑀_𝑠=𝑑𝑉𝑑𝑡^𝐿𝐶=𝑑()𝑑𝑡𝜀𝑉𝐶𝑑𝑡−_𝑘𝜀_𝑘𝑞_𝑘𝐶_𝑘−𝑅^′, (1) where 𝑀_𝑠 is the solute mass solved in the solvent, 𝑉_𝐿 is the suspension solids-free volume, 𝑉 is the total suspension volume, 𝐶 is the solute concentration in the continuous phase, 𝜀 is the fraction of the volume free of solids, 𝑞 is the volumetric flow into and out of the system, and 𝑅^′ 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 𝜌𝑉𝑐_𝑝𝑑𝑇𝑑𝑡=_𝑘𝑞_𝑘𝐻_𝑘−Δ𝐻_𝑐𝑅^′−𝐻_ext, (2) where 𝜌 is the slurry density, 𝑐_𝑝 is the specific heat capacity, 𝑇 is the crystallizer temperature, 𝐻_𝑘 is the enthalpy of the 𝑘th stream, Δ𝐻_𝑐 is the heat of crystallization, and 𝐻_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, 𝑛(𝐿)). 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 𝑛(𝐿) as function of time and crystal size: 𝜕𝑛+𝜕()𝜕𝑡𝑛𝐺𝜕𝐿=−_𝑘𝑞_𝑘𝑛_𝑘𝑉+𝐵𝛿𝐿−𝐿₀, (3) in which 𝑛 is the population density, 𝐺 is the crystal growth rate, 𝐿 is the internal coordinate (characteristic crystal length), 𝐵 is the nucleation rate density, and 𝛿(𝐿−𝐿₀) is the Dirac delta function acting at 𝐿=𝐿₀. The boundary and initial conditions for (3) are given as follows: ()()()(𝐿)BC∶𝑛0,𝑡=0𝑛∞,𝑡=0𝑡≥0,IC∶𝑛𝐿,0=𝑓0≤𝐿≤∞. (4) 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 𝑘_𝑔, 𝑘_𝑏, 𝑔, ℎ, 𝑜 and 𝑝 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 𝐺=𝑘_𝑔𝑆^𝑔_𝑟𝑁^ℎ_𝑟,𝐵=𝑘_𝑏𝑆^𝑏_𝑟𝑀^𝑜_𝑇𝑁^𝑝_𝑟. (5) The relative supersaturation 𝑆_𝑟 is evaluated as (𝐶−𝐶_sat)/𝐶_sat, 𝑁_𝑟 represents the agitation rate, and 𝑀_𝑇 is the total crystal mass per volume. Table 1 shows the kinetic parameters included in (5).

Table 1 

Adjusted parameters for constitutive equations Salcedo-Estrada.

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 𝑑𝜇₀=𝜇𝑑𝑡𝜏^𝑒0−𝜇₀𝜏+𝑘_𝑏𝑆^𝑏_𝑟𝑀^𝑜_𝑇𝑁^𝑝_𝑟,𝑑𝜇_𝑗=𝜇𝑑𝑡𝜏^𝑒𝑗−𝜇_𝑗𝜏+𝑘_𝑏𝑆^𝑔_𝑟𝑁^ℎ_𝑟𝑗𝜇_𝑗−1,𝑗=1,2,…,5,𝑑𝐶=()𝑑𝑡𝜀𝐶_𝑒−𝐶+𝐶𝜇^𝑒₃−3𝜏𝑘_𝑔𝑆^𝑔_𝑟𝑁^ℎ_𝑟𝜇₂𝜌_𝑐−𝐶𝜏1−𝜇₃,𝜌𝑉𝑐_𝑝𝑑𝑇𝑑𝑡=𝜌𝑞𝑐_𝑝𝑇_𝑒−𝑇−Δ𝐻_𝑐3𝜌_𝑐𝑘_𝑣𝑉𝑘_𝑔𝑆^𝑔_𝑟𝑁^ℎ_𝑟𝜇₂−𝑈𝐴𝑇−𝑇_𝑐,𝜌_𝑤𝑉_𝑐𝑐_𝑝𝑤𝑑𝑇_𝑐𝑑𝑡=𝜌_𝑤𝑞_𝑤𝑐_𝑝𝑤𝑇_𝑐𝑒−𝑇_𝑐+𝑈𝐴𝑇−𝑇_𝑐−𝑈₀𝐴₀𝑇_𝑐−𝑇_∞. (6)(7)(8)(9)(10) The mean crystal size and standard deviation can be evaluated using the following relationships of the moments in order to obtain the complete CSD: 𝐿₄₃=𝜇₄𝜇₃,SD_4,3=𝜇₅𝜇₃−𝜇₄𝜇₃². (11) 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, 𝑐_mz, for the ammonium sulfate-water system as function of temperature. The lower limit (saturation concentration, 𝑐_𝑠) was evaluated using the polynomial reported by Perry et al. [33]. Both (12) and (13) are valid in the range from 20 to 50°C 𝑐_mz=1.3179⁻⁷𝑇³+2.6136⁻⁶𝑇²+2.6280⁻³𝑐𝑇+0.70758,_𝑠=4.0×10⁻⁵𝑇²+2.0×10⁻⁴𝑇+0.73. (12)(13)

3. 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: 𝑥̇=𝑓(𝑥)+^𝑚_𝑖=1𝑔_𝑖(𝑥)𝑢_𝑖,𝑦_𝑖=ℎ_𝑖(𝑥),1≤𝑖≤𝑝. (14) The 𝑓(𝑥) and 𝑔_𝑖(𝑥) are smooth vector fields, ℎ₁(𝑥),…,ℎ_𝑝(𝑥) are smooth scalar fields, 𝑢_𝑖 are the manipulated variables, 𝑥̇ is state vector, and 𝑦_𝑖 are the output variables. The mathematical model for the crystallization process may be obtained making the following changes into (6)–(10), 𝜇_𝑖=𝑥_𝑖, 𝐶=𝑥₆,𝑇=𝑥₇,𝑇_𝑐=𝑥₈, and reordering the vector fields as follows: 𝑓(𝑥)=𝜇^𝑒₀𝜏−𝑥₀𝜏𝜇^𝑒₁𝜏−𝑥₁𝜏𝜇^𝑒₂𝜏−𝑥₂𝜏𝜇^𝑒₃𝜏−𝑥₃𝜏𝜇^𝑒₄𝜏−𝑥₄𝜏𝜇^𝑒₅𝜏−𝑥₅𝜏(𝜀𝐶)_𝑒−𝐶+𝐶𝜇^𝑒₃𝜏1−𝑥₃𝜌𝑞𝑐_𝑝𝑇_𝑒−𝑥₇−𝑈𝐴𝑥₇−𝑥₈𝜌𝑉𝑐_𝑝𝑈𝐴𝑥₇−𝑥₈−𝑈₀𝐴₀𝑥₈−𝑇_∞𝜌_𝑤𝑉_𝑐𝑐_𝑝𝑤^𝑇,𝑔₁(𝑥)=𝑘_𝑏𝑆𝑟^𝑏𝑀^𝑜_𝑇𝑁_𝑟^𝑝−1𝑘_𝑔𝑆𝑟^𝑔𝑥₀𝑁_𝑟^ℎ−12𝑘_𝑔𝑆𝑟^𝑔𝑥₁𝑁_𝑟^ℎ−13𝑘_𝑔𝑆𝑟^𝑔𝑥₂𝑁_𝑟^ℎ−14𝑘_𝑔𝑆𝑟^𝑔𝑥₃𝑁_𝑟^ℎ−15𝑘_𝑔𝑆𝑟^𝑔𝑥₄𝑁_𝑟^ℎ−1−3𝜏𝑘_𝑔𝑆𝑟^𝑔𝑁_𝑟^ℎ−1𝑥₂𝜌_𝑐−𝑥₆𝜏1−𝑥₃−3Δ𝐻_𝑐𝜌_𝑐𝑘_𝑣𝑉𝑘_𝑔𝑆𝑟^𝑔𝑁_𝑟^ℎ−1𝑥₂𝜌𝑉𝑐_𝑝0^𝑇,𝑔₂(𝑥)=𝑇00000000^𝑐𝑒−𝑥₈𝑉_𝑐^𝑇,𝑦₁=ℎ₁(𝑥)=𝑥₃,𝑦₂=ℎ₂(𝑥)=𝑥₇, (15)

where 𝑥∈𝑅⁹ and 𝑢∈𝑅².

3.1. Relative Degree

The relative degree of a multivariable nonlinear system may be calculated with the Lie derivatives of ℎ_𝑖(𝑥) along the vector fields 𝑓(𝑥) and 𝑔_𝑖(𝑥). It is defined as the smaller integer 𝑟_𝑖 such that 𝐿_𝑔𝑖𝐿^𝑟𝑖_𝑓⁻¹ℎ(𝑥)≠0 for all 1≤𝑖≤𝑚 in a neighborhood of 𝑥⁰ and an 𝑚×𝑚 nonsingular matrix defined as 𝐶(𝑥)=𝐿_𝑔1𝐿⁰_𝑓ℎ₁(𝑥)…𝐿_𝑔𝑚𝐿⁰_𝑓ℎ₁(𝑥)𝐿⋮⋱⋮_𝑔1𝐿^𝑚_𝑓ℎ_𝑚(𝑥)⋯𝐿_𝑔𝑚𝐿^𝑚_𝑓ℎ_𝑚(𝑥). (16) Considering ℎ₁(𝑥)=𝑥₃ and ℎ₂(𝑥)=𝑥₇ as the control objectives, 𝑢₁=𝑁_𝑟 and 𝑢₂=𝑞_𝑤 as the manipulated variables, and assuming a control problem with direct relationship between 𝑢₁→𝑥₃ and 𝑢₂→𝑥₇, the relative degree vector [𝑟₁,𝑟₂] is calculated from 𝐿_𝑔1𝐿^𝑟1_𝑓⁻¹ℎ₁(𝑥) and 𝐿_𝑔2𝐿^𝑟2_𝑓⁻¹ℎ₂(𝑥) different from zero. The evaluation gives relative degrees 𝑟₁=1 and 𝑟₂=2. The system overall relative degree is three assuming that (𝑇_𝑐𝑒−𝑥₈) and 𝑥₂ are different from zero. The first condition is fulfilled at normal cooling operation and the second in a seeded crystallization process. Matrix 𝐶(𝑥) can be constructed with the values given in (17)–(20). The evaluation of the determinant of 𝐶(𝑥) gives a value different from zero (nonsingular) when 𝑇_𝑐𝑒≠𝑥₈. This fact ensures the existence of both the inverse of 𝐶(𝑥) and the control law 𝐿_𝑔1𝐿⁰_𝑓ℎ₁(𝑥)=3𝑘_𝑔𝑆𝑟^𝑔𝑥₂𝑁_𝑟^ℎ−1𝐿≠0,_𝑔2𝐿⁰_𝑓ℎ₁(𝑥)𝐿=0,_𝑔1𝐿¹_𝑓ℎ₂(𝑥)=3𝜌𝑞𝑐_𝑝+𝑈𝐴Δ𝐻_𝑐𝜌_𝑐𝑘_𝑣𝑘_𝑔𝑥₂𝑆𝑟^𝑔𝑁_𝑟^ℎ−1𝜌²𝑉𝑐_𝑝²,𝐿_𝑔2𝐿¹_𝑓ℎ₂(𝑥)=𝑈𝐴𝜌𝑉𝑐_𝑝𝑇_𝑐𝑒−𝑥₈𝑉_𝑐,det𝐶(𝑥)=3𝑘_𝑔𝑆𝑟^𝑔𝑥₂𝑁_𝑟^(ℎ−1)𝑈𝐴𝑇_𝑐𝑒−𝑥₈𝜌𝑉𝑐_𝑝𝑉_𝑐. (17)(18)(19)(20)(21)

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 𝑛−𝑟_𝑚 additional functions 𝜙_{𝑛−𝑟𝑚} to complete the map 𝜙(𝑥) in such a way that the Jacobian of 𝜙(𝑥) is nonsingular at 𝑥⁰. 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(𝐽_𝜙(𝑥))=−𝑈𝐴/𝜌𝑉𝑐_𝑝≠0𝜙₄=𝑥₀,𝜙₅=𝑥₁,𝜙₆=𝑥₂,𝜙₇=𝑥₄,𝜙₈=𝑥₅,𝜙₉=𝑥₆. (22) Finally, the complete map can be represented as follows:𝜙(𝑥)=𝑧₁=ℎ₁(𝑥)𝑧₂=ℎ₂(𝑥)𝑧₃=𝐿_𝑓ℎ₂(𝑥)𝑧₄𝑧₅𝑧₆𝑧₇𝑧₈𝑧₉^𝑇=𝑥₃𝑥₇𝜌𝑞𝑐_𝑝𝑇_𝑒−𝑥₇−𝑈𝐴𝑥₇−𝑥₈𝜌𝑉𝑐_𝑝𝑥₀𝑥₁𝑥₂𝑥₄𝑥₅𝑥₆^𝑇. (23)

The system given by (14) may be represented in its new coordinates set as follows: 𝑧̇₁=𝐿_𝑓ℎ₁(𝑥)+𝐿_𝑔1ℎ₁(𝑥)𝑢₁,𝑧̇₂=𝐿_𝑓ℎ₂(𝑥),𝑧̇₃=𝐿²_𝑓ℎ₂(𝑥)+𝐿_𝑔1𝐿_𝑓ℎ₂(𝑥)𝑢₁+𝐿_𝑔2𝐿_𝑓ℎ₂(𝑥)𝑢₂,𝑧̇₄=𝑥^𝑒₀−𝑥₀𝜏+𝑘_𝑏𝑆𝑟^𝑏𝑀^𝑜_𝑇𝑢^𝑝₁,𝑧̇₅=𝑥^𝑒₁−𝑥₁𝜏+𝑘_𝑔𝑆𝑟^𝑔𝑥₀𝑢^ℎ₁,𝑧̇₆=𝑥^𝑒₂−𝑥₂𝜏+2𝑘_𝑔𝑆𝑟^𝑔𝑥₁𝑢^ℎ₁,𝑧̇₇=𝑥^𝑒₄−𝑥₄𝜏+4𝑘_𝑔𝑆𝑟^𝑔𝑥₃𝑢^ℎ₁,𝑧̇₈=𝑥^𝑒₅−𝑥₅𝜏+5𝑘_𝑔𝑆𝑟^𝑔𝑥₄𝑢^ℎ₁,𝑧̇₉=(𝜀𝐶)_𝑒−𝑥₆+𝑥₆𝑥^𝑒₃−3𝜏𝑘_𝑔𝑆𝑟^𝑔𝑥₂𝜌_𝑐−𝑥₆𝑢^ℎ₁𝜏1−𝑥₃, (24) where 𝑧₁₋₃∈𝑅³ represent the system controllable and observable parts, while 𝑧₄₋₉∈𝑅⁶ represent the nonobservable and uncontrollable parts of the system (internal dynamic) (Slotine and Li [34]).

3.3. Zero Dynamics

The zero dynamics of the system described by (24) is obtained when the output errors are equal to zero, for example, 𝑥₃−𝑥^∗₃=𝑧₁=0, 𝑥₇−𝑥^∗₇=𝑧₂=0 and 𝑧̇₁=𝑧̇₂=𝑧̇₃=0. The system may be represented as follows: 𝑧̇₄=𝑥^𝑒₀−𝑥₀𝜏+𝑘_𝑏𝑆𝑟^𝑏𝑀^𝑜_𝑇𝑢^𝑝₁,𝑧̇₅=𝑥^𝑒₁−𝑥₁𝜏+𝑘_𝑔𝑆𝑟^𝑔𝑥₀𝑢^ℎ₁,𝑧̇₆=𝑥^𝑒₂−𝑥₂𝜏+2𝑘_𝑔𝑆𝑟^𝑔𝑥₁𝑢^ℎ₁,𝑧̇₇=𝑥^𝑒₄−𝑥₄𝜏+4𝑘_𝑔𝑆𝑟^𝑔𝑥₃𝑢^ℎ₁,𝑧̇₈=𝑥^𝑒₅−𝑥₅𝜏+5𝑘_𝑔𝑆𝑟^𝑔𝑥₄𝑢^ℎ₁,𝑧̇₉=(𝜀𝐶)_𝑒−𝑥₆+𝑥₆𝑥^𝑒₃−3𝜏𝑘_𝑔𝑆𝑟^𝑔𝑥₂𝜌_𝑐−𝑥₆𝑢^ℎ₁𝜏1−𝑥₃. (25) 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: 𝑧̇=𝐴𝑥, where 𝐴 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 𝑁_𝑟 and 𝑞_𝑤 and plotted together with the upper (12) and lower limits (13) of the metastable zone. It was found that combination of values of 𝑁_𝑟 smaller than 300 rpm and cooling flow rates bigger than 2200 cm³/min produces crystallizer concentrations outside the metastable zone (see Figure 1). Therefore, it was selected a combination of 𝑁_𝑟 = 420 rpm and 𝑞_𝑤 = 800 as the steady state starting operation conditions that ensure the equilibrium operation with a dynamic closed-loop response asymptotically stable. The eigenvalues of 𝐴 at the selected conditions were −𝟎.𝟑𝟖𝟒𝟔,−𝟎.𝟑𝟗𝟗𝟐,−𝟎.𝟓𝟔𝟕𝟗+𝟎.𝟑𝟐𝟓𝟐𝐢, and −𝟎.𝟓𝟔𝟕𝟗−𝟎.𝟑𝟐𝟓𝟐𝐢. 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³/min and agitation rate of 420 rpm. These values were used to generate the nonsingular matrix 𝐴 shown in the following:𝐴=−0.38460−41.6617001.042𝐸30.0049−0.3846−14.624100000.0098−0.765800000−7.6𝐸−4−0.38460000−5𝐸−50.0247−0.3846000000−0.3846. (26)

Table 2 

Equilibrium state variables calculated at operation conditions 𝑞_𝑤 = 800 cm³/min and 𝑁_𝑟 = 420 rpm.

Figure 1 

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

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 𝑥^∗₃ and 𝑥^∗₇ under normal operating conditions. The proposed pairing 𝑢₁→𝑥₃ and 𝑢₂→𝑥₇ do not consider interactions between the other input and output variables; for example, input variable 𝑣_𝑖 modifies only the output variable 𝑦_𝑖. The control law may be expressed as follows: 𝑦^𝑟𝑖_𝑖=𝑣_𝑖. (27) The control actions may be calculated solving the system given by (28), where 𝐶 represents the uncoupled system matrix given previously in (16): 𝑈=𝐶⁻¹𝑣₁−𝐿^𝑟1_𝑓ℎ₁(𝑥)𝑣₂−𝐿^𝑟2_𝑓ℎ₂(𝑥)𝑢₁,=𝑁_𝑟=𝑣₁𝜏−𝑥^𝑒₃+𝑥₃3𝜏𝑘_𝑔𝑆𝑟^𝑔𝑥₂^1/ℎ. (28)(29) Using a proportional control action only, 𝑣₁=−𝑘₀(𝑥₃−𝑥_3𝑠𝑝), (29) becomes 𝑢₁=𝑁_𝑟=−𝑘₀𝑥₃−𝑥_3𝑠𝑝𝜏−𝑥^𝑒₃+𝑥₃3𝜏𝑘_𝑔𝑆𝑟^𝑔𝑥₂^1/ℎ. (30) The cooling flow rate, 𝑞_𝑤, might be evaluated as follows:𝑢₂=𝑞_𝑤=Δ𝐻_𝑐𝜌_𝑐𝑘_𝑣𝑉_𝑐𝜌𝑞𝑐_𝑝+𝑈𝐴𝑣₁𝜏−𝑥^𝑒₃+𝑥₃𝜌𝑐_𝑝𝑈𝐴𝜏𝑥₈−𝑇_𝑐𝑒−𝑉_𝑐𝑐_𝑝𝑉𝜌𝑥₈−𝑇_𝑐𝑒×𝑣𝑈𝐴₂−−𝜌𝑞𝑐_𝑝−𝑈𝐴𝜌𝑞𝑐_𝑝𝑇_𝑒−𝑥₇−𝑈𝐴𝑥₇−𝑥₈𝜌²𝑉²𝑐²_𝑝−𝑈𝐴𝑈𝐴𝑥₇−𝑥₈−𝑈₀𝐴₀𝑥₈−𝑇_inf𝜌𝑉𝑐_𝑝𝜌_𝑤𝑉_𝑐𝑐_𝑝𝑤, (31)𝑢₂=𝑞_𝑎=𝑆^𝑔_𝑟𝑥₂𝐶₁𝑢₁𝑥₈−𝑇_𝑐𝑒−𝐶₂𝑥₈−𝑇_𝑐𝑒×𝑣₂−𝐶₃𝑇_𝑒−𝑥₇+𝑥₇−𝑥₈𝐶₄−𝐶₅+𝐶₆𝑥₈−𝑇_inf, (32) where 𝐶₁=3Δ𝐻_𝑐𝜌_𝑐𝑘_𝑣𝑉_𝑐𝑘_𝑔𝜌𝑞𝑐_𝑝+𝑈𝐴𝜌𝑐_𝑝𝑈𝐴,𝐶₂=𝑉_𝑐𝑐_𝑝𝑉𝜌,𝐶𝑈𝐴₃=−𝜌𝑞𝑐_𝑝−𝑈𝐴𝜌𝑞𝑐_𝑝𝜌²𝑉²𝑐²_𝑝,𝐶₄=−𝜌𝑞𝑐_𝑝−𝑈𝐴𝑈𝐴𝜌²𝑉²𝑐²_𝑝,𝐶₅=(𝑈𝐴)²𝜌𝑉𝑐_𝑝𝜌_𝑤𝑉_𝑐𝑐_𝑝𝑤,𝐶₆=(𝑈𝐴)𝑈₀𝐴₀𝜌𝑉𝐶_𝑝𝜌_𝑤𝑉_𝑐𝑐_𝑝𝑤,𝑣₂=−𝑘₂𝑦̇₂−𝑘₁𝑦₂=−𝑘₂𝑥̇₇−𝑘₁𝑥₇−𝑥_7𝑠𝑝. (33)

4. 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: 𝑁_𝑟∈ [300   800] rpm and 𝑞_𝑤∈ [100   8000] cm³. 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.

Table 3 

Operating conditions and physical constants used in the model.

Simulations were performed using a seeded inlet stream with crystals of the same size 𝐿₄₃=129 μm at 𝑇_𝑒 of 30°C and a flow rate of 1153.84 cm³/min. The initial conditions for the state variables used are shown on Table 2. An initial equilibrium operation was assumed at 𝑁_𝑟 = 420 rpm and 𝑞_𝑤 = 800 cm³. 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 𝑘₀ = 1.063, 𝑘₁ = 2.001, and 𝑘₂ = 1.065. The PID parameters for the temperature controller are 𝐾_𝑐 = 98.872, 𝜏_𝐼 = 995.651, and 𝜏_𝐷 = 10.223, and the PID parameters for third moment controller are 𝐾_𝑐 = 9993.881, 𝜏_𝐼 = 0.991, and 𝜏_𝐷 = 0. 011.

Run 1
In order to evaluate the servo control for the third moment (steady-state value of 0.0040), two step changes (−0.0015 and +  0.0015) were introduced at times 20 and 60. The reference value for the crystallizer temperature was kept constant at 23°C. Figures 2(a) and 2(c) show the closed-loop dynamic responses obtained when traditional PID and nonlinear control (NLC) were implemented. For both controllers, the system reaches the set point values, but the IAEs for 𝜇₃ and 𝑇 are smaller when the NLC algorithm is implemented. Figures 2(b) and 2(d) show the manipulated variables effort for both types of controllers, and no appreciable difference was found. The stirring rate has a large effect on the distribution moments and a small effect on the crystallizer temperature.

(a)

(b)

(c)

(d)

Figure 2 

Closed-loop dynamic responses for the MIMO system when step changes in 𝜇3 are introduced.

Figure 3 illustrates the metastable zone limit (𝐶_mz), the solubility profile (𝐶_𝑠), 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 𝜇₃.

(a)

(b)

Figure 3 

Concentration profiles inside the crystallizer when step changes are imposed to the output variable 𝜇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 𝜇₃ and 𝑇 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).

(a)

(b)

(c)

(d)

Figure 4 

Closed-loop dynamic responses for the MIMO system when step changes in the crystallizer temperature are introduced.

(a)

(b)

Figure 5 

Concentration profiles inside the crystallizer when step changes are imposed to the output variable 𝑇.

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 𝜇₃ and 𝑇 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.

Table 4 

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

(a)

(b)

(c)

(d)

Figure 6 

Closed-loop dynamic responses for the MIMO system when a disturbance change of ±1°C were used in the feed temperature.

5. 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.

Acknowledgment

The authors want to thank DGEST (Project 2600.09-P) for providing financial support.