# Nonlinear MIMO Control of a Continuous Cooling Crystallizer

Authors: Pedro Alberto Quintana-Hernández; Raúl Ocampo-Pérez; Salvador Tututi-Avila; Salvador Hernández-Castro
Journal: Modelling and Simulation in Engineering
 (2012)

---

## 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 (μ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.


---

## Body

## 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)dMsdt=d(VLC)dt=d(εVC)dt-∑kεkqkCk-R′,
where Ms is the solute mass solved in the solvent, VL 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(2)ρVcpdTdt=∑kqkHk-ΔHcR′-Hext,
where ρ is the slurry density, cp is the specific heat capacity, T is the crystallizer temperature, Hk is the enthalpy of the kth stream, ΔHc is the heat of crystallization, and Hext 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:
(3)∂n∂t+∂(nG)∂L=-∑kqknkV+Bδ(L-L0),
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-L0) is the Dirac delta function acting at L=L0. The boundary and initial conditions for (3) are given as follows:
(4)BC:n(0,t)=0n(∞,t)=0t≥0,IC:n(L,0)=f(L)0≤L≤∞.
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 kg, kb, 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
(5)G=kgSrgNrh,B=kbSrbMToNrp.
The relative supersaturation Sr is evaluated as (C-Csat)/Csat, Nr represents the agitation rate, and MT 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.
Equation
Parameter
Value
Units
k
g
9.090× 10−4
(cm/min) (min/Rev)h
Growth
g
1.500
Dimensionless
h
1.337
Dimensionless
Nucleation
k
b
180
#/min cm3 (cm3/g)o (min/Rev)p
b
0.562
Dimensionless
o
0.001
Dimensionless
p
0.050
DimensionlessDifferent 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
(6)dμ0dt=μ0eτ-μ0τ+kbSrbMToNrp,(7)dμjdt=μjeτ-μjτ+kbSrgNrhjμj-1,j=1,2,…,5,(8)dCdt=(εC)e-C+Cμ3e-3τkgSrgNrhμ2(ρc-C)τ(1-μ3),(9)ρVcpdTdt=ρqcp(Te-T)-ΔHc3ρckvVkgSrgNrhμ2-UA(T-Tc),(10)ρwVccpwdTcdt=ρwqwcpw(Tce-Tc)+UA(T-Tc)-U0A0(Tc-T∞).
The mean crystal size and standard deviation can be evaluated using the following relationships of the moments in order to obtain the complete CSD:
(11)L43=μ4μ3,SD4,3=μ5μ3-(μ4μ3)2.
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, cmz, for the ammonium sulfate-water system as function of temperature. The lower limit (saturation concentration, cs) 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
(12)cmz=1.3179-7T3+2.6136-6T2+2.6280-3T+0.70758,(13)cs=4.0×10-5T2+2.0×10-4T+0.73.


## 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:
(14)x˙=f(x)+∑i=1mgi(x)ui,yi=hi(x),1≤i≤p.
The f(x) and gi(x) are smooth vector fields, h1(x),…,hp(x) are smooth scalar fields, ui are the manipulated variables, x˙ is state vector, and yi are the output variables. The mathematical model for the crystallization process may be obtained making the following changes into (6)–(10), μi=xi, C=x6,T=x7,Tc=x8, and reordering the vector fields as follows: (15)f(x)=(ρqcp(Te-x7)-UA(x7-x8)ρVcpμ0eτ-x0τμ1eτ-x1τμ2eτ-x2τμ3eτ-x3τμ4eτ-x4τμ5eτ-x5τ(εC)e-C+Cμ3eτ(1-x3)ρqcp(Te-x7)-UA(x7-x8)ρVcpUA(x7-x8)-U0A0(x8-T∞)ρwVccpwUA(x7-x8)-U0A0(x8-T∞)ρwVccpw)T,g1(x)=(-3τkgSrgNrh-1x2(ρc-x6)τ(1-x3)kbSrbMToNrp-1kgSrgx0Nrh-12kgSrgx1Nrh-13kgSrgx2Nrh-14kgSrgx3Nrh-15kgSrgx4Nrh-1-3τkgSrgNrh-1x2(ρc-x6)τ(1-x3)-3ΔHcρckvVkgSrgNrh-1x2ρVcp0-3τkgSrgNrh-1x2(ρc-x6)τ(1-x3))T,g2(x)=(00000000Tce-x8Vc)T,y1=h1(x)=x3,y2=h2(x)=x7,wherex∈R9 and u∈R2.


### 3.1. Relative Degree

The relative degree of a multivariable nonlinear system may be calculated with the Lie derivatives ofhi(x) along the vector fields f(x) and gi(x). It is defined as the smaller integer ri such that LgiLfri-1h(x)≠0 for all 1≤i≤m in a neighborhood of x0 and an m×m nonsingular matrix defined as
(16)C(x)=(Lg1Lf0h1(x)…LgmLf0h1(x)⋮⋱⋮Lg1Lfmhm(x)⋯LgmLfmhm(x)).
Considering h1(x)=x3 and h2(x)=x7 as the control objectives, u1=Nr and u2=qw as the manipulated variables, and assuming a control problem with direct relationship between u1→x3 and u2→x7, the relative degree vector [r1,r2] is calculated from Lg1Lfr1-1h1(x) and Lg2Lfr2-1h2(x) different from  zero. The evaluation gives relative degrees r1=1 and r2=2. The system overall relative degree is three assuming that (Tce-x8) and x2 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 Tce≠x8. This fact ensures the existence of both the inverse of C(x) and the control law
(17)Lg1Lf0h1(x)=3kgSrgx2Nrh-1≠0,(18)Lg2Lf0h1(x)=0,(19)Lg1Lf1h2(x)=3(ρqcp+UA)ΔHcρckvkgx2SrgNrh-1ρ2Vcp2,(20)Lg2Lf1h2(x)=UAρVcp(Tce-x8)Vc,(21)detC(x)=3kgSrgx2Nr(h-1)UA(Tce-x8)ρVcpVc.


### 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 findn-rm additional functions ϕn-rm to complete the map ϕ(x) in such a way that the Jacobian of ϕ(x) is nonsingular at x0. 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/ρVcp≠0(22)ϕ4=x0,ϕ5=x1,ϕ6=x2,ϕ7=x4,ϕ8=x5,ϕ9=x6.
Finally, the complete map can be represented as follows:(23)ϕ(x)=(z1=h1(x)z2=h2(x)z3=Lfh2(x)z4z5z6z7z8z9)T=(x3x7ρqcp(Te-x7)-UA(x7-x8)ρVcpx0x1x2x4x5x6)T.The system given by (14) may be represented in its new coordinates set as follows:
(24)z˙1=Lfh1(x)+Lg1h1(x)u1,z˙2=Lfh2(x),z˙3=Lf2h2(x)+Lg1Lfh2(x)u1+Lg2Lfh2(x)u2,z˙4=(x0e-x0)τ+kbSrbMTou1p,z˙5=(x1e-x1)τ+kgSrgx0u1h,z˙6=(x2e-x2)τ+2kgSrgx1u1h,z˙7=(x4e-x4)τ+4kgSrgx3u1h,z˙8=(x5e-x5)τ+5kgSrgx4u1h,z˙9=(εC)e-x6+x6x3e-3τkgSrgx2(ρc-x6)u1hτ(1-x3),
where z1-3∈R3 represent the system controllable and observable parts, while z4-9∈R6 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, x3-x3*=z1=0, x7-x7*=z2=0 and z˙1=z˙2=z˙3=0. The system may be represented as follows:
(25)z˙4=(x0e-x0)τ+kbSrbMTou1p,z˙5=(x1e-x1)τ+kgSrgx0u1h,z˙6=(x2e-x2)τ+2kgSrgx1u1h,z˙7=(x4e-x4)τ+4kgSrgx3u1h,z˙8=(x5e-x5)τ+5kgSrgx4u1h,z˙9=(εC)e-x6+x6x3e-3τkgSrgx2(ρc-x6)u1hτ(1-x3).
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: z˙=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 ofNr and qw and plotted together with the upper (12) and lower limits (13) of the metastable zone. It was found that combination of values of Nr smaller than 300 rpm and cooling flow rates bigger than 2200 cm3/min produces crystallizer concentrations outside the metastable zone (see Figure 1). Therefore, it was selected a combination of Nr = 420 rpm and qw = 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 cm3/min and agitation rate of 420 rpm. These values were used to generate the nonsingular matrix A shown in the following:(26)A=(-0.38460-41.6617001.042E30.0049-0.3846-14.624100000.0098-0.765800000-7.6E-4-0.38460000-5E-50.0247-0.3846000000-0.3846).Table 2
Equilibrium state variables calculated at operation conditionsqw = 800 cm3/min and Nr = 420 rpm.
x
0
x
1
x
2
58.6272
0.7646
0.01984
x
3
x
4
x
5
0.0007682
3.955
E
-
5
2.545
E
-
6
x
6
x
7
x
8
0.7706
24.943
19.917Figure 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 valuesx3* and x7* under normal operating conditions. The proposed pairing u1→x3 and u2→x7 do not consider interactions between the other input and output variables; for example,  input variable vi modifies only the output variable yi. The control law may be expressed as follows:
(27)yiri=vi.
The control actions may be calculated solving the system given by (28), where C represents the uncoupled system matrix given previously in (16):
(28)U=C-1[v1-Lfr1h1(x)v2-Lfr2h2(x)],(29)u1=Nr=(v1τ-x3e+x33τkgSrgx2)1/h.
Using a proportional control action only, v1=-k0(x3-x3sp), (29) becomes
(30)u1=Nr=(-k0(x3-x3sp)τ-x3e+x33τkgSrgx2)1/h.
The cooling flow rate, qw, might be evaluated as follows:(31)u2=qw=ΔHcρckvVc(ρqcp+UA)(v1τ-x3e+x3)ρcpUAτ(x8-Tce)-VccpVρ(x8-Tce)UA×{v2-(-ρqcp-UA)[ρqcp(Te-x7)-UA(x7-x8)]ρ2V2cp2-UA[UA(x7-x8)-U0A0(x8-Tinf)]ρVcpρwVccpw},(32)u2=qa=Srgx2C1u1(x8-Tce)-C2(x8-Tce)×[v2-C3(Te-x7)+(x7-x8)(C4-C5)+C6(x8-Tinf)],
where
(33)C1=3ΔHcρckvVckg(ρqcp+UA)ρcpUA,C2=VccpVρUA,C3=(-ρqcp-UA)ρqcpρ2V2cp2,C4=(-ρqcp-UA)UAρ2V2cp2,C5=(UA)2ρVcpρwVccpw,C6=(UA)U0A0ρVCpρwVccpw,v2=-k2y˙2-k1y2=-k2x˙7-k1(x7-x7sp).


## 3.1. Relative Degree

The relative degree of a multivariable nonlinear system may be calculated with the Lie derivatives ofhi(x) along the vector fields f(x) and gi(x). It is defined as the smaller integer ri such that LgiLfri-1h(x)≠0 for all 1≤i≤m in a neighborhood of x0 and an m×m nonsingular matrix defined as
(16)C(x)=(Lg1Lf0h1(x)…LgmLf0h1(x)⋮⋱⋮Lg1Lfmhm(x)⋯LgmLfmhm(x)).
Considering h1(x)=x3 and h2(x)=x7 as the control objectives, u1=Nr and u2=qw as the manipulated variables, and assuming a control problem with direct relationship between u1→x3 and u2→x7, the relative degree vector [r1,r2] is calculated from Lg1Lfr1-1h1(x) and Lg2Lfr2-1h2(x) different from  zero. The evaluation gives relative degrees r1=1 and r2=2. The system overall relative degree is three assuming that (Tce-x8) and x2 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 Tce≠x8. This fact ensures the existence of both the inverse of C(x) and the control law
(17)Lg1Lf0h1(x)=3kgSrgx2Nrh-1≠0,(18)Lg2Lf0h1(x)=0,(19)Lg1Lf1h2(x)=3(ρqcp+UA)ΔHcρckvkgx2SrgNrh-1ρ2Vcp2,(20)Lg2Lf1h2(x)=UAρVcp(Tce-x8)Vc,(21)detC(x)=3kgSrgx2Nr(h-1)UA(Tce-x8)ρVcpVc.


## 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 findn-rm additional functions ϕn-rm to complete the map ϕ(x) in such a way that the Jacobian of ϕ(x) is nonsingular at x0. 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/ρVcp≠0(22)ϕ4=x0,ϕ5=x1,ϕ6=x2,ϕ7=x4,ϕ8=x5,ϕ9=x6.
Finally, the complete map can be represented as follows:(23)ϕ(x)=(z1=h1(x)z2=h2(x)z3=Lfh2(x)z4z5z6z7z8z9)T=(x3x7ρqcp(Te-x7)-UA(x7-x8)ρVcpx0x1x2x4x5x6)T.The system given by (14) may be represented in its new coordinates set as follows:
(24)z˙1=Lfh1(x)+Lg1h1(x)u1,z˙2=Lfh2(x),z˙3=Lf2h2(x)+Lg1Lfh2(x)u1+Lg2Lfh2(x)u2,z˙4=(x0e-x0)τ+kbSrbMTou1p,z˙5=(x1e-x1)τ+kgSrgx0u1h,z˙6=(x2e-x2)τ+2kgSrgx1u1h,z˙7=(x4e-x4)τ+4kgSrgx3u1h,z˙8=(x5e-x5)τ+5kgSrgx4u1h,z˙9=(εC)e-x6+x6x3e-3τkgSrgx2(ρc-x6)u1hτ(1-x3),
where z1-3∈R3 represent the system controllable and observable parts, while z4-9∈R6 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, x3-x3*=z1=0, x7-x7*=z2=0 and z˙1=z˙2=z˙3=0. The system may be represented as follows:
(25)z˙4=(x0e-x0)τ+kbSrbMTou1p,z˙5=(x1e-x1)τ+kgSrgx0u1h,z˙6=(x2e-x2)τ+2kgSrgx1u1h,z˙7=(x4e-x4)τ+4kgSrgx3u1h,z˙8=(x5e-x5)τ+5kgSrgx4u1h,z˙9=(εC)e-x6+x6x3e-3τkgSrgx2(ρc-x6)u1hτ(1-x3).
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: z˙=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 ofNr and qw and plotted together with the upper (12) and lower limits (13) of the metastable zone. It was found that combination of values of Nr smaller than 300 rpm and cooling flow rates bigger than 2200 cm3/min produces crystallizer concentrations outside the metastable zone (see Figure 1). Therefore, it was selected a combination of Nr = 420 rpm and qw = 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 cm3/min and agitation rate of 420 rpm. These values were used to generate the nonsingular matrix A shown in the following:(26)A=(-0.38460-41.6617001.042E30.0049-0.3846-14.624100000.0098-0.765800000-7.6E-4-0.38460000-5E-50.0247-0.3846000000-0.3846).Table 2
Equilibrium state variables calculated at operation conditionsqw = 800 cm3/min and Nr = 420 rpm.
x
0
x
1
x
2
58.6272
0.7646
0.01984
x
3
x
4
x
5
0.0007682
3.955
E
-
5
2.545
E
-
6
x
6
x
7
x
8
0.7706
24.943
19.917Figure 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 valuesx3* and x7* under normal operating conditions. The proposed pairing u1→x3 and u2→x7 do not consider interactions between the other input and output variables; for example,  input variable vi modifies only the output variable yi. The control law may be expressed as follows:
(27)yiri=vi.
The control actions may be calculated solving the system given by (28), where C represents the uncoupled system matrix given previously in (16):
(28)U=C-1[v1-Lfr1h1(x)v2-Lfr2h2(x)],(29)u1=Nr=(v1τ-x3e+x33τkgSrgx2)1/h.
Using a proportional control action only, v1=-k0(x3-x3sp), (29) becomes
(30)u1=Nr=(-k0(x3-x3sp)τ-x3e+x33τkgSrgx2)1/h.
The cooling flow rate, qw, might be evaluated as follows:(31)u2=qw=ΔHcρckvVc(ρqcp+UA)(v1τ-x3e+x3)ρcpUAτ(x8-Tce)-VccpVρ(x8-Tce)UA×{v2-(-ρqcp-UA)[ρqcp(Te-x7)-UA(x7-x8)]ρ2V2cp2-UA[UA(x7-x8)-U0A0(x8-Tinf)]ρVcpρwVccpw},(32)u2=qa=Srgx2C1u1(x8-Tce)-C2(x8-Tce)×[v2-C3(Te-x7)+(x7-x8)(C4-C5)+C6(x8-Tinf)],
where
(33)C1=3ΔHcρckvVckg(ρqcp+UA)ρcpUA,C2=VccpVρUA,C3=(-ρqcp-UA)ρqcpρ2V2cp2,C4=(-ρqcp-UA)UAρ2V2cp2,C5=(UA)2ρVcpρwVccpw,C6=(UA)U0A0ρVCpρwVccpw,v2=-k2y˙2-k1y2=-k2x˙7-k1(x7-x7sp).


## 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: Nr∈ [300   800] rpm and qw∈ [100   8000] cm3. 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.
Symbol
Value
Units
τ
2.6
min
ρ
c
1.769
gr/cm3
ρ
1.244
gr/cm3
c
p
0.6816
cal/gr°C
Δ
H
c
−12.0354
cal/gr
V
3000
cm3
V
c
820
cm3
C
p
w
1
cal/gr°C
T
∞
25
°C
T
c
e
13
°C
U
A
1000
cal/min°C
U
0
A
0
100
cal/min°C
μ
0
e
0.9772
μ
1
e
1.004× 10−2
μ
2
e
1.8213× 10−4
μ
3
e
2.3407× 10−5
μ
4
e
3.0159× 10−8
μ
5
e
3.8891× 10−10Simulations were performed using a seeded inlet stream with crystals of the same sizeL43=129 μm at Te of 30°C and a flow rate of 1153.84 cm3/min. The initial conditions for the state variables used are shown on Table 2. An initial equilibrium operation was assumed at Nr = 420 rpm and qw = 800 cm3. 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 k0 = 1.063, k1 = 2.001, and k2 = 1.065. The PID parameters for the temperature controller are Kc = 98.872, τI = 995.651, and τD = 10.223, and the PID parameters for third moment controller are Kc = 9993.881, τI = 0.991, and τD = 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. Figures2(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 μ3 and T 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.Closed-loop dynamic responses for the MIMO system when step changes inμ3 are introduced.
(a)
(b)
(c)
(d)Figure3 illustrates the metastable zone limit (Cmz), the solubility profile (Cs), 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.Concentration profiles inside the crystallizer when step changes are imposed to the output variableμ3.
(a)
(b)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. Figure4 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).Closed-loop dynamic responses for the MIMO system when step changes in the crystallizer temperature are introduced.
(a)
(b)
(c)
(d)Concentration profiles inside the crystallizer when step changes are imposed to the output variableT.
(a)
(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 Figure6. 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.Table 4
IAE values generated when NLC and PID control strategies were used.
NLC
PID
Run
T
μ
3
× 103
T
μ
3
× 103
1
2.6055
9.60
5.5421
12.80
2
3.8599
6.10
7.9254
9.80
3
3.1019
6.10
9.5816
9.60Closed-loop dynamic responses for the MIMO system when a disturbance change of ±1°C were used in the feed temperature.
(a)
(b)
(c)
(d)


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



---

*Source: 912071-2012-10-22.xml*