HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 33(1-2). pp. 97–104 (2005) MODEL PREDICTIVE CONTROL OF CONTINUOUS CRYSTALLIZERS N. MOLDOVÁNYI1 and B.G. LAKATOS2 1Honeywell Process Solutions, H-1139 Budapest Petneházy u. 2-4 HUNGARY, 2Department of Process Engineering, University of Veszprém, H-8201 Veszprém, P.O. Box 158, HUNGARY The problem of model predictive control of continuous isothermal crystallizers, using a detailed moment equation model is analysed. The mean size of the crystalline product and the variance of crystal size are the controlled variables, while the manipulated variables are the input concentration of the solute and the flow-rate. The controllability and observability, as well as the coupling between the inputs and the outputs are analyzed by simulation using the linearised model. The crystallizer has proved to be a nonlinear multi-input multi-output system with strong coupling between the state variables. It is shown that the mean crystal size and the variance of can be controlled nearly separately by the residence time and the inlet solute concentration, respectively. By seeding, the controllability of the crystallizer increases significantly. The linear model predictive controller synthesized using the moment equation model appears to be an efficient controller for continuous crystallizers. Keywords: Model predictive control; continuous isothermal crystallizer, computer simulation Introduction Model predictive control (MPC) refers to a class of com- puter control algorithms that utilize an explicit process model to predict the future response of the plant. At each control interval an MPC algorithm attempts to optimize future plant behaviour by computing a sequence of future manipulated variable adjustments. The first input in the optimal sequence is then sent into the plant, and the entire calculation is repeated at subsequent control intervals [1]. Originally developed to meet the specialized control needs of power plants and petroleum refineries, MPC technology can now be found in a wide variety of ap- plication areas including chemicals (Honeywell has MPC at polypropylene units at TVK, Hungary), food pro- cessing (Honeywell is right now working on a dairy pro- duct unit at the UK), automotive and aerospace applica- tions. The presented work is an opening to an another new application, the MPC control of continuous crystal- lizers. Crystallization is a widely used cleaning, separation and grain producing technique in the chemical industry, particularly at the pharmaceutical works. From the point of view of controlling a crystallizer the main quality crite- ria are the properties of the produced crystals, first of all the size-distribution and the mean size. Crystallization is a multi-variable system, with multi input and multi output (MIMO) often with strong coupl- ing. Thus a good, up-to-date control is possible using a model-based MIMO control system. There are only few examples in the literature for this [2-5]. Since one can do nothing to change the size distribution of the crystals in the system once crystals have grown beyond a stable nucleus size. Therefore a predictive type of control would be better than the corrective type. One of the main prob- lems is that for a proper model-based control of the size- distribution, because of the mentioned properties of the population balance equation, high-order control is re- quired, which means serious difficulties. But the crystal- lizers are dissipative systems [6], so that a crystallizer as a dynamical system possesses finite dimensional global attractors [7] that create an adequate basis for the syn- thesis and usage of good quality, low-order model-based control systems. At the same time it means that for the synthesis of the model-based control system of crystal- lizers the moment equation model, generated from popu- lation balance equation governing the crystal size distri- bution can be used. Chui and Cristofides [8] applied this property to design a nonlinear SISO controller. In this paper a model-based MIMO control system of a continuous isothermal crystallizer is presented. For the synthesis of the control system a multi-variable state- space model is composed. Linear controllability and ob- servability analysis is presented, and the coupling of the inputs and the outputs are analysed. The efficiency of the developed model predictive controller is demonstrated by simulation. 98 Concept of MPC In model predictive control, the control action is pro- vided after solving – on-line at each sampling instant – an optimization problem, and the first element in the op- timized control sequence is applied to the process (re- ceding horizon control). The “moving horizon” concept of MPC is a key feature that distinguishes it from clas- sical controllers, where a pre-computed control law is employed. The major factor of the success of predictive control is its applicability to problems where analytic control law is difficult, or even impossible to obtain. The methodology of all the controllers belonging to the MPC family is characterized by the following strat- egy, represented in Fig.1 (y is the output, w is the set- point and u is the input): FUTUREPAST Contr ol horizon Predicti on horizon u(.) w(.) y(.) k k+1k-1k-2 k+2 Fig 1. MPC horizons Prediction horizon (Hp) represents the number of samples taken from the future over which MPC com- putes the predicted process variable profile and mini- mizes the predicted error. The control signals change only inside the control horizon, Hc remaining constant afterwards 1,...,),1()( −=−+=+ pcc HHjHkujku (1) The basic steps: 1. As it is shown, in the MPC future outputs for a de- termined prediction horizon Hp are predicted at each instant k using a prediction model. These pre- dicted outputs pHjkjky ,...1),(ˆ =+ (means the value at the instant k+j, calculated at instant k) depend on the known values up to instant k (past inputs and outputs) and the future control signals 1,...0),( −=+ pHjkjku , which are those to be sent to the system and to be calculated. 2. The set of control signals is calculated by op- timizing a cost function in order to keep the process as close as possible to the reference trajectory . This criterion usually takes the form of a quadratic function of the errors between the predicted output signal and the reference tra- jectory. The control effort is included in the ob- jective function in most of the cases. An explicit solution can be obtained if the criterion is quadratic, the model is linear and there are no constraints, otherwise an iterative optimization method has to be used. pHjjkw ,...1),( =+ 3. The control signal )( kku is sent to the process whilst the next control signals calculated are rejected, be- cause at the next sampling instant y(k+1) is already known and step 1 is repeated with this new value and all the sequences are brought up to date. Thus the )11( ++ kku is calculated (which in principle will be different to the )1( kku + because of the new information available) using a receding ho- rizon concept. In order to implement this strategy, the basic struc- ture shown is Fig.2 is used. A model is used to predict the future plant outputs, based on past and current values and on the proposed optimal future control ac- tions. These actions are calculated by the optimizer taking into account the cost function (where the future tracking error is considered) as well as the constraints. Optimizer Model Past Inputs and Outputs Future Inputs Cost Function Constraints Future Errors Predicted Outputs Reference Trajectory - + Fig.2: Basic structure of MPC The process model plays, in consequence, a decisive role in the controller. The chosen model must be cap- able of capturing the process dynamics so as to pre- cisely predict the future outputs as well as being simple to implement and to understand. As MPC is not a unique technique but a set of different methodologies, there are many types of models used in various formula- tions. Honeywell uses mostly black-box models at the refineries, getting them by stepping the plant. The new tendency is using chemical engineering, so called “grey- box” models. The presented case study clearly fills this requirement. The optimizer is another fundamental part of the strategy as it provides the control actions. If the cost function is quadratic, its minimum can be obtained as an explicit function (linear) of past inputs and outputs and the future reference trajectory. In the presence of in- equality constraints the solution has to be obtained by more computationally taxing algorithms. The size of op- timization problems depends on the number of variables and on the prediction horizon used and usually turn out to be relatively modest optimization problems which do not require sophisticated computer codes to be solved. However the amount of time needed for the constrained and robust cases can be various orders of magnitude higher than that needed for the unconstrained case and 99 the bandwidth of the process to which constrained MPC can be applied is considerably reduced. For a continuous-time model (the cost function is discrete), the MPC problem can be represented as [ ,)(),(min )( kUkYJ kU Ψ= ] (2) [ ] tHttttutxftx p∆+≤≤= *,*)(*),(*)(ˆ (3) (3) [ ] tHtttHttHtutu pcc ∆+≤≤∆−+∆−+= *)1(,)1(*)( [ ],)(),(),(0 kUkYkXΦ= (5) )(0 kDU≥ (6) )(0 kDY≥ (7) where ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + + ≡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + + ≡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ −+ + ≡ )( )1( )( )( , )( )1( )( )(, )1( )1( )( )( kHkx kkx kkx kX kHky kky kky kY kHku kku kku kU p pc M MM . Here, )( kku is the input calculated from informa- tion available at time k, )(ku )( kky is the output cal- culated from information available at time k, Hc is the control horizon and Hp is the prediction horizon, while x denotes the state variable. Constraint (3) corresponds to satisfaction of the continuous-time model equations over the prediction horizon, while (4) enforces the re- quirement that all inputs beyond the control horizon are held constant. Algebraic equation (5) represents con- straints for the model, and for the sake of completeness Eqs (6) and (7) correspond to the constraints on the in- put and output variables, respectively. )(ky The process model is assumed to have the following discrete-time representation, (8) [ ),(),()1( kukxFkx =+ ] ] (9) [ ,)()( kxhky = where x is the n-dimensional vector of state variables, u is the m-dimensional vector of manipulated input vari- ables, and y is the p-dimensional vector of controlled output variables. Such a model can be obtained by dis- cretizing a continuous-time, state-space model or by de- riving a state-space realization of a discrete-time, input- output model. It is important to note that time delays can be handled by augmenting the state vector such that the resulting state-space model has no delays. The optimization problem for the prototypical MPC formulation is [9]: [ ] [ ,)(),(),( )(min 1 0 )1(),...1(),( ∑ − = −++ +∆++ ++= p c H j p kHkukkukku kjkukjkukjkyL kHkyJ φ ] (10) where )1()()( kjkukjkukjku −+−+=+∆ , φ and L are (possibly) (non)linear functions of their arguments. The optimization problem is solved to the constraints discussed below. The functions φ and L can be chosen to satisfy a wide variety of objectives, including minim- ization of overall process cost. However, economic op- timization may be performed by a higher-level system which determines appropriate setpoints for the MPC controller. In this case it is meaningful to consider quad- ratic functions of the following form: [ ] [ ] [ ] [ ] ),j(k)j(k )()ju(k)()ju(k )()jy(k)()jy(k kuSku kukRkuk kykQkykL T s T s s T s +∆+∆+ −+−++ −+−+= (11) [ ] [ ])()H(k)()H(k pp kykyQkyky sTs −+−+=φ (12) where and are steady-state targets for u and y, respectively, and Q, R, S are positive definite weighting matrices. The principal controller tuning parameters are Hc, Hp, Q, R, S and the sampling period ∆t. )(kus )(kys The prediction outputs are obtained from the model (8-9). Successive iterations of the model equations yield [ ] [ ][ ] [ ] [ ] [ ][ ] [ ] [ ],)1(),...1(),(),()( ,)1(),(),( ,)1(,)(),( ,)1(),1()2( ,)(),( ,)(),()1()1( 2 1 1 1 kjkukkukkukxGkjky kkukkukxG kkukkukkxFG kkukkxGkky kkukxG kkukkxFhkkxhkky j −++=+ +≡ +≡ ++=+ ≡ =+=+ M (13) where )()( kxkkx = is a vector of current state vari- ables. If the control horizon (Hc) is less than the pre- diction horizon (Hp), the output predictions are ge- nerated by setting inputs beyond the control horizon equal to the last computed value: .,)1()( pcc HjHkHkukjku ≤≤−+=+ Note that the prediction )( kjky + depends on the current stable variables, as well as on the calculated input sequence. Therefore, MPC requires measurements or estimates of the state variables. Solution of the MPC problems yields the input se- quence { })1(),...,1(),( kHkukkukku c −++ Only the first input vector in the sequence is actually imple- mented: )()( kkuku = . Then the prediction horizon is moved forward one time step, and the problem is re- solved using new process measurements. This receding horizon formulation yields improved closed-loop per- formance in the presence of unmeasured disturbances and modelling errors. 100 Case Study: Continuous crystallizer Moment equation model The mathematical model of a continuous MSMPR crystallizer consists of the population balance equation for crystals, of the balance equations for sol-vent and crystallizing substance, and of the equations describing the variation of the equilibrium saturation concentration. In the present analysis, the crystallizer is assumed to be isothermal, thus the equilibrium saturation concentration c* is constant during the course of the process. It is assumed that the following conditions are satisfied: (1) the volumetric feed and withdrawal rates of the crystallizer are constant and equal, thus the working volume is constant during the course of the operation; (2) the crystals can be characterized by a linear di- mension L; (3) all new crystals are formed at a nominal size Ln≅0 so that we assume Ln=0; (4) crystal breakage and agglomeration are neg- ligible; (5) no growth rate fluctuations occur; (6) the overall linear growth rate of crystals G is size-dependent and has the form of the power law ex- pression of supersaturation ; (14) ( aLcckG gg +−= 1*)( ) (7) the primary nucleation rate Bp is described by the Volmer model ⎟⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −= * ln exp 2 c c k kB epp ε (15) the secondary nucleation rate Bs is described by the power law relation (16) jbbb cckB 3*)( µ−= where µ3 is the third of the ordinary moments of the population density function n, which are defined as (17) ...3,2,1,0,),( 0 == ∫ ∞ mdLtLnLmmµ With these assumptions the population balance equa- tion governing the crystal size dynamics becomes: [ ] [ 0,0 ),(),( ),(*),,(),( >> −=⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Lt tLntLnq L tLncLcG t tLn V in∂ ∂ ∂ ∂ ] (18) subject to the following initial and boundary conditions: (19) 0),()0,( 0 ≥= LLnLn 0,*),,(),(*),,(lim 0 ≥== → tbpccBtLncLcG L νν (20) (21) 0,0),(lim ≥= ∞→ ttLn L Here, n(L,t)dL expresses the number of crystals having sizes in the range L to L+dL at time t in a unit volume of suspension. The mass balance of the crystallizing substance has the form [ ] [ cinc cqqcdt cd V ρεε ρεε )1( )1( −+−= −+ ] (22) with the initial condition (23) 0)0( cc = where the voidage of suspension ε is related to n and L by (24) ∫ ∞ −=−= 0 3 3 ),(11 dLtLnLkk VV µε Finally, the mass balance of the solvent is written in the form ( ) svsvin sv cqqc dt cd V ε ε −= (25) with the initial condition . (26) 0)0( svsv cc = Therefore, the state at time t≥0 of the continuous isothermal MSMPR crystallizer is given by the triple [c(t),csv(t),n(t)], and its dynamics is described by the dis- tributed parameter model formed by the mixed set of partial and ordinary differential eqs (18), (22) and (25), subject to the initial and boundary conditions (19-21) (23) and (26). The evolution in time of this system oc- curs in the state space R2×N that is the Descartes product of the vector space R2 of concentrations and of the function space N of the population density func- tions. Consideration of dynamical problems of crystal- lizers in this product space, however, seems to be quite complex and not constructive. In the present study, we concentrate on a reduced case, considering the problem in a finite dimensional state space model based on the moments of the population density function instead of the distributed parameter system (18)-(21). Since the overall crystal growth rate (14) is a linear function of size L, the population balance Eq. (18) can be converted into an infinite set of recursive ordinary differential equations for the moments of population density function: ( ) bpVBq dt d V in ,,00 0 =+−= νµµ µ ν (27) ( ) )(*)( 1min mmggmm accVmkqdt d V µµµµ µ +−+−= − m=1,2,3 (28) which can be closed by Eq.(22), describing the mass balance of the crystallizing substance, at the equation for the third order moment. Then Eq.(22) takes the form )(*))((3 1 )( 32 µµρεε acccVkkcc q dt dc V gcgVin +−−−−= (29) while Eq.(12) can be rewritten as )(*)(3 1 )( 32 µµεε accVckkcc q dt dc V gsvgVsvsvin sv +−−−= (30) 101 where csv stands for the concentration of solvent. Here, because of the selective withdrawal, the voidage in the crystallizer and that in the outlet stream are not equal. Therefore, the first four moment equations from the system (27-28) with Eqs (29-30) provide a closed mo- ment equations model of the crystallizer. Dimensionless equations. Scaling We introduce the following set of dimensionless vari- ables svincsvinsvcsvincin cmmmt csycsyccsy ccsymsxts ==−= −==== ,*),( *),(,3,2,1,0,, µξ into Eqs (27)-(30), where st, sc and sm, m=0,1,2,3, are scale factors defined as *}max{ 1 : cc s in c − = , , ( gintgV ccskks 3330 *}max{6: −= − ) ( ) gintgV ccskks 2221 *}max{6: −= − ( )gintgV ccskks *}max{3: 12 −= − Vks =:3 and max{cin} denotes the maximal value of inlet con- centration, as well as the set of dimensionless para- meters q V sts tt ==:τ , ( 1*}max{*)(: −−−= ccc inρα ) ( )gintg ccask *}max{: 1 −= −β ( ) gintgVpap ccskkkD 343 *}max{6: −= − ( ) bgintgjVbab ccskkkD +−− −= 3431 *}max{6: *}max{ * *: cc c cs in c − ==γ Then the dimensionless governing equations take the form: bp xx d dx in ,,000 =Θ+ − = ν τξ ν (31) 3,2,1),( 1 min =++ − = − mxmxy xx d dx mm gmm β τξ (32) ( ) 3 32 3 1 )3()( 1 x xxyy x yy d dy gin − +− − − − = βα τξ (33) ( ) 3 32 3 1 )3( 1 x xxyy x yy d dy gsvsvsvinsv − + + − − = β τξ (34) subject to the initial conditions 000 )0(,)0(,3,2,1,0,)0( svsvmm yyyymxx ==== where ( ) ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + −−= γ γ Θ y k xD eapp 2 3 ln exp1 (36) and jb abb xyD 3=Θ . (37) It follows from physical reasoning that the phys- ically admissible solutions to Eqs (31-34) should sat- isfy the constraints 110,0 ,0,0,0 minmax33 3 1 022 3 2 01100 <−=<≤≤≤ ≤≤≤≤≤≤ εxxxAx xAxxxyy m mmin (38) where x0m denotes the maximal value of the zero order moment, while ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = ρ ε * 1 min c V Vsv (39) where Vsv is the volume of solvent in the crystallizer. The parameters, which in the case of primary and secondary nucleation form the vectors of real numbers pp=(τ,α,g,β,Dap,ke,γ) and pb=(τ,α,g,β,Dab,b,j), res- pectively, are also bounded: 0,0,0,0 0,0,,0,0,0 min ≥≥≥≥ ≥≥≥≥≥≥ γ ββατ jbk DDg e abap (40) As a consequence, the state of crystallizer (31)-(34) is represented by the vector of variables (x0,x1,x2,x3,y, ysv), and its time evolution occurs in the feasible region of solutions (38) of the six-dimensional state space R6. The behaviour of crystallizer in the neighbourhood of a stationary state may be deduced by examining the eigenvalues of the Jacobian matrix of Eqs (36-38) at this state which becomes ( ) ( ) ( ) ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − + − − − − − − − − − − − − − − − − − − − − − )1( 1 )1( )( 1 3 1 00 0 1 )(3 1 )( 00 0 1 300 00 1 20 000 1 000 1 3 3 333 55 33 33 22 11 1514 S S SS svSsvin S g SsvS S g SsvS S S g SS S g SS S inSg S g S S inSg S g S S inSg S g S SS x x xy yyg x yy x yy j x yy x yy y xx gyy y xx gyy y xx gyy jj τ β αβα ττ β ττ β ττ β τ νν (41) where in the case of primary nucleation ( )S Sin pS x xx j 3 00 14 1− − = τ and ( ) ( ) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + + − = γ γ γτ SS inSe pS y y xxk j 3 00 15 ln 2 (42) ( )[ ] ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − −−− +− − = )( )(1 1 1 3 55 SS SSSin S S yy ygyyy x j ατ α τ while for secondary nucleation S S sS x x jj 3 0 14 = and S S sS y sj 015 = x (43) In order to formulate the state-space model, we define ( ) ( )6543213210 ,,,,,,,,,, uuuuuuwyxxxx ininininin ==u (44) where τττ Ω ==== Vqtq Vs ss V qsw . (45) Now: 102 bp xx w d dx in ,,000 =Θ+ Ω − = ν ξ ν (46) )( 1 min mm gmm xmxy xx w d dx β Ωξ ++ − = − m=1,2,3 (47) ( ) 3 32 3 1 )3()( 1 x xxyy x yy w d dy gin − +− − −Ω − = βα ξ (48) ( ) 3 32 3 1 )3( 1 x xxyy x yy w d dy gsvsvsvinsv − + − −Ω − = β ξ (49) subject to the initial conditions 000 )0(,)0(,3,2,1,0,)0( svsvmm yyyymxx ==== (50) i.e. bp xu u d dx ,,016 0 =Θ+ Ω − = ν ξ ν (51) )( 16 mm gmmm xmxy xu u d dx β Ωξ ++ − = − , m=1,2,3 (52) ( ) 3 32 3 5 6 1 )3()( 1 x xxyy x yu u d dy g − +− − −Ω − = βα ξ (53) ( ) 3 32 3 5 6 1 )3( 1 )( x xxyy x yuy u d dy gsvsvsvinsv − + − −Ω − = β ξ (54) where ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −=⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −= ρ ρ ρ ρ 5 u s y sy csv in csvsvin (55) To summarize as a control engineering problem: vector of state-variables is x=(x0,x1,x2,x3,,y,ysv), its changes rep- resented by a nonlinear state space model (51)-(55); the input vector of is u=(x0in,x1in,x2in,x3in,yin,w) and the out- put is defined as y=x. Analysis of the model Stability and bifurcation In linear dynamics, one seeks the fundamental solutions from which one can build all other solutions. In non- linear dynamics, the main questions are: What is the qualitative behaviour of the system? Which and how many non-wandering sets (i.e. a fixed point, a limit cycle, a quasi-periodic or chaotic orbit) occur? Which of them are stable? How does the number of non- wandering sets change while changing a parameter of the system (called control parameter)? The appearance and disappearance of a non-wandering set is called a bifurcation. Change of stability and bifurcation always coincide. The number of attractors in a nonlinear dynamical system can change when a system parameter is changed. This change is called bifurcation. It is accompanied by a change of the stability of an attractor. In a bifurcation point, at least one eigenvalue (λ) of the Jacobian matrix gets a zero real part. There are three generic types of so- called co-dimension-one bifurcations (the term co- dimension counts the number of control parameters for which fine tuning is necessary to get such a bifurcation). Back to the crystallizer, changing the value of ke, the parameter of primary nucleation rate and observing the supersaturation, Hopf bifurcation occurs as it shown in Fig.3. 0,45 0,4 ymin ymax 0,35 0,3 Fig. 3: Bifurcation diagram ke-yS of the crystallizer For further studies of controlling the crystallizer an operating point has been chosen from the region of stable steady states exhibiting only damped oscillations, that is at ke=0.01. Controllability and observability There are two basic problems we need to consider. The first one is the coupling between the input and the state: Can any state be controlled by the input? This is a con- trollability problem. Another is the relationship between the state and the output: Can all the information about the state be observed from the output? This is an ob- servability problem. For the controllability and observability test a lin- earized model (at the operating point) of the nonlinear system was used. xcy BuAxx T= +=& (56) where the state transition matrix (A) is the Jacobi matrix of the system, the input matrix (B) can also be derived from the model and the output matrix (cT) is a diagonal matrix. For a MIMO system the necessary and sufficient conditions for the system to have completely control- lability is nRank =⎥⎦ ⎤ ⎢⎣ ⎡ − BABAABB 1n2 L (57) For the general system the necessary and sufficient condition of a linear system for complete observability is nRank T nTTT =⎥⎦ ⎤ ⎢⎣ ⎡ − cAcAcAc 12 )()( L (58) The results of the calculation is that the linearized system is completely controllable and observable at a certain operating point. 0 0,0 0,1 0,15 0,2 0,25 bifurcation point yS damped oscillations 0,001 0,01 0,1 limit-cycle oscillations 1ke 103 Relative –gain array )1()1( ))(ˆ)(())(ˆ)(( ),,,,( 1 21 2 1 −+−++ +−++−+ = ∑ ∑ = = jkjk jkjkjkjk HHHJ c p p H j T T H Hj cpp ∆uR∆u ywQyw QR (61) The relative-gain array provides exactly a methodology, whereby we select pairs of input and output variables in order to minimize the amount of iteration among the resulting loops. It was first proposed by Bristol and today is a very popular tool for the selection of control loops. where denotes the predicted process outputs, and are the minimum and the maximum prediction horizons, is the control horizon, Q and R are positive definite weighting matrices. )(ˆ jk +y 1pH 2pH cH In our crystallization system the control variables can be the mean size of crystals, the variance of the crystal size (σ2) and the productivity, i.e. the total vol- ume of the crystals: 33 2 0 1 0 2 2 0 1 1 ,, xx x x x x x =⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −== υυυ (59) The manipulated variables are the input concentra- tion of the solute and the flow-rate: . (60) wuyu in ==== 6251 , ϑϑ An optimization algorithm will be applied to com- pute a sequence of future control signals that minimizes the cost function. For unconstrained control based on linear process model(s) and quadratic cost function the control sequence can be analytically calculated. After tuning the , the and 3=cH 11 =pH 52 =pH . By seeding, the controllability of the crystallizer in- crease, the overshoots and the oscillation are smaller. The results of the controlling study have shown that the linear MPC is an adaptable and feasible controller as it illustrated by Fig.4. Here, the first two rows are the out- puts (solid lines) with the corresponding setpoints (dashed lines), while the third and fourth rows present the time variations of the inputs of the crystallizer. Note that since the volume of the crystal suspension was kept constant the mean residence time was varied by changing the volumetric feed. For the crystallizer we have two outputs and two inputs, there are three possible pairs of control variables, so three different relative-gain arrays can be formed and computed (The value 0 and 1 are rounded.): ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ 2 1 2 1 21.12821.127 21.12721.128 ϑ ϑ υ υ , ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ 2 1 3 1 10 01 ϑ ϑ υ υ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ 2 1 3 2 10 01 ϑ ϑ υ υ The results show that controlling the mean size and the variance together would be very difficult. However, by putting crystal grains to the input (seeding), the control of the variance also becomes possible. The new different relative-gain array: Fig.4. Performance of the MPC of the continuous isothermal MSMPR crystallizer ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∆ ∆ 2 1 2 1 24.124.0 24.024.1 ϑ ϑ υ υ In this case, the mean-size and the variance can be nearly separately controlled. For the further experiments these two outputs will be selected. The system is very sensible to the quality and the quantity of the seeding. It is assumed to be fixed to a suitable operating point. Results of simulations The cost function is chosen to satisfy a wide variety of objectives, including minimization of overall process costs. However, economic optimization may be per- formed by a higher-level system which determinates the appropriate setpoints for the controller. In this case cost function is formulated reflecting the reference tracking error and the control action. The general expression of such an objective function is 104 Next steps References 1. QIN S. J. and BADGWELL T.A Control Engineering Practice, 2003, 11, 733-764 The presented model is going to be developed in UniSim, the simulation software of Honeywell, to make possible to connect it to the Profit Controller. the MPC of Honeywell. The final step is to control a real con- tinuous crystallizer at a plant. 2. MYERSON A.S. et al. Proc. 10th Symposium. Industrial Crystallization. Academia, Praha, 1987 3. JAGER, J., et al. Powder Technology, 1992, 69, 11 4. MILLER, S.M. and RAWLINGS, J.B., AIChE Journal, 1994, 40, 1312 Conclusions 5. ROHANI, S. et al. Computers and Chemical Engin- eering, 1999, 23, 279. 6. LAKATOS, B.G. and SAPUNDZHIEV, Ts.J., Bulgarian Chemical Communications, 1997, 29, 28 The moment equation model of a continuous isothermal MSMPR crystallizer was presented and the model was analyzed. Better control of the variance of the crystal size was possible by introducing some seeding into the crystallizer. By seeding, the controllability of the crys- tallizer increase, the overshoots and the oscillation are smaller. The results of the controlling study have shown that the linear MPC is an adaptable and feasible con- troller for continuous crystallizers. 7. TEMAM, R., Infinite-Dimensional Dynamical Sys- tems in Mechanics and Physics. Springer-Verlag, New York, 1988 8. CHIU, T. and CHRISTOFIDES, P.D., AIChE Journal, 1999, 45, 1279 9. MEADOWS, E.S. and RAWLINGS, J.B, Model predictive control, Englewoods Cliffs. NJ: Prentice- Hall, 1997 Nonlinear process control, Chapter 5. (233-310).