NMPC of an industrial crystallization process using model-based observers

Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

Contents lists available at ScienceDirect

Journal of Industrial and Engineering Chemistry journal homepage: www.elsevier.com/locate/jiec

NMPC of an industrial crystallization process using model-based observers Ce´dric Damour a, Michel Benne a,*, Lionel Boillereaux b, Brigitte Grondin-Perez a, Jean-Pierre Chabriat a a b

Lab. of Energetic, Electronics and Processes, University of La Reunion, 15 Av. Cassin, BP 7151, 97715 Saint-Denis, France GEPEA, UMR CNRS 6144, ENITIAA, rue de la Ge´raudie`re, BP 82225, F-44322 Nantes cedex 3, France

A R T I C L E I N F O

A B S T R A C T

Article history: Received 29 January 2010 Accepted 4 April 2010

This paper illustrates the benefits of a nonlinear model-based predictive control (NMPC) approach applied to an industrial crystallization process. This relevant approach proposes a setpoint tracking of the crystal mass. The controlled variable, unavailable, is obtained using an extended Luenberger observer. A neural network model is used as internal model to predict process outputs. An optimization problem is solved to compute future control actions taking into account real-time control objectives. The performances of this strategy are demonstrated via simulation in cases of setpoint tracking and disturbance rejection. The results reveal a significant improvement in terms of robustness and energy efficiency. ß 2010 The Korean Society of Industrial and Engineering Chemistry. Published by Elsevier B.V. All rights reserved.

Keywords: Crystallization Nonlinear model predictive control Industrial control Extended Luenberger observer Artificial neural network

1. Introduction Crystallization is a unit operation for achieving the extraction of a solute from a saturated solution. The presence of a continuous phase and a dispersed phase gives rise to physicochemical phenomena such as nucleation, growth, agglomeration and dissolution. Despite of the long history and widespread application of crystallization processes, due to the strong nonlinearity of this process, its monitoring and control remain an interesting challenge in terms of quality and global efficiency improvement. Software sensors have proved to be useful to improve crystallization monitoring. Especially, due to their robustness, observers are used to estimate the purity, the crystal size distribution (CSD) or the solute concentration [1–3]. At the same time, advanced control algorithms are developed to improve the process control (nonlinear modelbased predictive control, generic model control, H1-control, etc.) for both continuous and batch crystallization processes [4–13]. In cane sugar industry, in order to improve exhaustion, crystallization is achieved grade wise through multiple stages, with a decrease in purity from the first to the last stage. In the first stage, the main objective is a good product quality, such as a good CSD. In this stage, the concentration of impurities being relatively small, a wide variety of mathematical models can be found in the scientific literature to describe the sugar crystallization phenomena, both steady state and dynamic models [14–17]. The most widely approach to estimate the crystallization rate involves the population balance [18,19] that gives information

* Corresponding author. Tel.: +33 262 262 938 223; fax: +33 262 262 938 673. E-mail address: [email protected] (M. Benne).

about CSD such as the mean particle size and deviation, or the number of particles. The widely used dynamic model for first stage is based on nine ordinary differential equations. In the second and more especially the last stages of extraction, the control objective consists in a maximal sucrose exhaustion of the solution. The main information remains the mass of crystals, without consideration of the CSD. More, the large presence of impurities prevent from using the widely employed momentum equations. In all stages the control strategy is designed to maintain the supersaturation s of the solution in a specific range (1.00 < s < 1.35) [20]. In practice, the electrical conductivity of the solution k is used to represent the supersaturation s and the process control is based on an electrical conductivity setpoint tracking using a proportional integral derivative (PID) controller [21,22]. Yet, experimental observations emphasize the shortcomings of this control strategy, especially in the presence of a high level of impurities. Indeed, k is strongly influenced by the variability of the experimental conditions and the presence of solid phases as well as various ionic compounds in the raw material [23,24]. Besides, the current PID control strategy has proved to be inadequate regarding the strong nonlinearities of the crystallization process [25]. This statement of fact leads to the proposition of alternative controlled variable and a more suitable control algorithm to improve the process control and the global process efficiency. In this paper, due to its robustness and its successful implementation in numerous industries [26,27], an NMPC approach is presented to improve the process control, based on a setpoint tracking control. The crystal mass is used as controlled variable and the nonlinear predictive approach is used as control algorithm instead of the linear PID controller. Due to the

1226-086X/$ – see front matter ß 2010 The Korean Society of Industrial and Engineering Chemistry. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jiec.2010.07.014

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

literature, the main process nonlinearities are including in the crystallization rate, which is usually solved from the population balance. Regarding the exhaustion objective of the last crystallization stage, the main information remains the crystal mass, without consideration of the CSD. Assuming nucleation, dissolution and agglomeration negligible during the growing phase [28,29] in supersaturated conditions, this point allows considering a new expression of the crystallization rate. In this study, the proposed model considers an explicit conversion of dissolved sucrose to crystallized sucrose thanks to a crystallization conversion rate parameter acryst.

Nomenclature Thermodynamic variables Brix (mass percentage of dry substance) (%) Bx% cc crystal content (%) Cp specific heat capacity (J K kg1) F flow rate (m3 s1) h specific enthalpy (J kg1) L specific latent heat (J kg1) m mass (kg) m˙ mass flow rate (kg s1) % purity (mass percentage of sugar) (%) Pte Pty purity (mass fraction of sugar) Q˙ heating power (W) T temperature (8C) W˙ stirring power (W) s supersaturation k conductivity (%) s concentration a adjusted parameter r density (kg m3)

2.1. Energy and mass balance The last crystallization step is performed through a semi-batch vacuum pan, operating as continuously stirred tank reactors (CSTR, Fig. 1). The dynamic model proposed in this paper is derived from four ordinary differential equations (1)–(4) that represent the mass balance for dissolved sucrose ms, crystals mc, water mc and emitted vapor mvap. Bx%f Pte%f dms dmc ¼ rfF f  dt 100 100 dt

(1)

dmc ¼ acryst ms dt

(2)

Subscripts i impurities c, Crys crystal s dissolved sucrose w water cw condensate (from heating steam) f feed hs heating steam mc mass of crystals ml mother liquor mg magma (mother liquor and crystals) vap emitted vapor

unavailability of on-line crystal mass measurement the proposed control strategy is based on observer. An artificial neural network model based on the estimate mass of crystal is used as internal model to predict the process output. The first step in observer design is the model formulation. In this study, a model dedicated to last stage crystallization is presented. In this model, the sugar is either dissolved or as crystals, without consideration of the CSD, which leads to a simplification compared to the well-known momentum approach. The crystallization is thus represented by a kinetics of conversion of the sugar between solution and crystals, considered as irreversible. This model reduces the problem formulation to four ordinary differential equations. An extended Luenberger observer is designed based on this model to achieve our objectives of on-line monitoring and control. This paper is organized as follow: the second section focuses on the model design and validation. The extended Luenberger observer design, including observability analysis and experimental validation, is presented in section three. Section four deals with the NMPC design. Finally, in fifth section, the benefits of the proposed control strategy are presented and discussed via simulation results.

709

Bx%f dmw ¼ rfF f 1  dt 100

! þ rw F w 

dmvap dt

Q˙  r f F f h f  rw F w hw  acryst ms ðhc  hml Þ=ðhml  hvap Þ dmvap ¼ 1 þ ððhml þ hc Þ=ðhml  hvap ÞÞ dt ðhml þ hc Þðr f F f ðBx%f =100ÞðPte%f =100Þ þ ðdmi ðtÞ=dtÞ ... þ

þr f F f ð1  ðBx%f =100ÞÞ þ rw F w Þ=ðhml  hvap Þ 1 þ ððhml þ hc Þ=ðhml  hvap ÞÞ

[(Fig._1)TD$IG]

(4)

2. Model design In this section, a model dedicated to last stage crystallization monitoring and control is developed. In the usual dedicated

(3)

Fig. 1. CSTR used for sugar crystallization.

710

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

With, mi ðtÞ ¼ mi ð0Þ þ r f

Bx%f 100

1

Pte%f 100

! F ft

and acryst [s1] an adjustable parameter that represents the crystallization conversion rate and which is a function of the 3 2 x˙ 1 acryst 6 x˙ 7 6 a 6 2 7 6 cryst 6 7¼6 4 x˙ 3 5 4 k6 k5 2

x˙ 4

k6 k5

32

0

0

0

0

0

3 2 k1 x1 7 6 6 0 0 76 x 2 7 7 6 76 7 þ 6 0 54 x3 5 4 k2 ð1 þ k4 Þ þ k6 ðr f h f  k3 Þ

0

0

0

0

x4

k2 k4 ð1 þ k4 Þ þ k6 ðk3  r f h f Þ

Proof. The energy balance takes into account both power supplies and specific enthalpy terms associated to mass transport and phase change: ˙ heating flux  Q,  hf and hw, feeding and water dilution specific enthalpy terms, respectively  (hvap  hml) and (hc  hml), evaporation and crystallization specific latent heat, respectively  (hml + hc), specific enthalpy of the solution It results that the energy balance can be described by Eq. (5), where Tmg represents the temperature of the sucrose solution:

Q˙ þ r f F f h f þ rw F w hw  ðdmvap =dtÞðhvap  hml Þ þðdmc =dtÞðhc  hml Þ  ðdmmg =dtÞðhml þ hc Þ mmg C pmg

Q˙ þ r f F f h f þ rw F w hw  ðdmvap =dtÞðhvap  hml Þ þðdmc =dtÞðhc  hml Þ  ððdms =dtÞ þ ðdmc =dtÞ dT mg þðdmw =dtÞ þ ðdmi ðtÞ=dtÞÞðhml þ hc Þ ¼ mmg C pmg dt

(5)

(6)

¼

Q˙ þ r f F f h f þ rw F w hw þ acryst ms ðhc  hml Þ ... ðhvap  hml Þ ðr f F f ðBx%f =100ÞðPte%f =100Þ þ r f F f ð1  ðBx%f =100ÞÞ



þrw F w  ðdmvap =dtÞ þ ðdmi ðtÞ=dtÞÞðhml þ hc Þ ðhvap  hml Þ (7)

Finally, the derivative of the emitted vapor with respect to time can be described by Eq. (4). 2.2. State space representation The nonlinear state space representation deduced from Section 2.1 is given by:  X˙ ¼ AX þ BU Y ¼ gðXÞ

0 0

rw ðð1 þ k4 Þ  hw k6 Þ rw ðk4 ð1 þ k4 Þ þ hw k6 Þ

0

3

2 3 u1 0 7 76 7 74 u 2 5 k6 5 u3 k6

with, k1 ¼ r f

Bx%f Pte%f 100 100

k2 ¼ r f 1 

Bx%f

!

100

k3 ¼ ðhml þ hc Þ r f

k4 ¼

Bx%f

!

100

hml þ hc hml  hvap

k5 ¼ acryst ðhc  hml Þ k6 ¼

1 þ k4 hml  hvap

In this case, three on-line industrial measurements are available: the total dry mater mass percentage, Bx% mg (given by a microwave sensor), the dissolved dry mater mass percentage, Bx% (given by a lm refractometer sensor) and the mass of emitted vapor mvap . The crystal content cc is obtained by an indirect measurement using the following expression: 100  cc ¼

Assuming a constant temperature of the solution and using the corresponding expression for dmw =dt, dms/dt and dmc/dt lead to Eq. (7). dmvap dt

the state and input matrix are expressed as follow:

0

supersaturated state and the temperature of the sucrose solution. The concentration of the solution is deduced from the estimate state by the ratio of the mass of dissolved sucrose to the mass of water. It can be notice that the emitted vapor mvap is deduced from the energy balance (5).

dT mg ¼ dt

According to the following state and input vectors: 2 3 2 3 2 3 2 3 ms x1 u1 Ff 6 mc 7 6 x2 7 4 6 6 7 7 X¼4 ¼ 4 5; U ¼ F w 5 ¼ 4 u 2 5 5 mw x3 u3 Q˙ mVa p x4

% Bx% mg  Bxlm

100  Bx% lm

This approach leads to the measured output vector Y and the nonlinear output function g(x): " #     x2 cc y1 ¼ ; gðxÞ ¼ x1 þ x2 þ x3 þ mi ðtÞ Y¼ mVa p y2 x4 with cc ¼

mc : ms þ mc þ mw þ mi ðt Þ

In the considered process, it illustrated that the use of a ‘‘pseudo-output variable’’ added to the measured outputs is required to fill the observability condition. Let us consider the ‘‘pseudo-output’’ variable ms+c that represents the entire mass of sucrose contained in the solution (dissolved sucrose ms plus crystal mc). Due to the perfect knowing of initial conditions ms+c is a reliable measurement. With, msþc ðtÞ ¼ msþc ð0Þ þ r f

Bx%f Pte%f 100 100

u1 t

This approach leads to the new measured outputs vector Y and the nonlinear output function g(x): 2 3 2 3 x2 cc 6 7 x þ x þ x þ m ðtÞ 2 3 i Y ¼ 4 msþc 5; gðxÞ ¼ 4 1 5 x1 þ x2 mVa p x4

[(Fig._3)TD$IG]

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

711

In this approach, only the output function of the model is nonlinear that allows considering an extended Luenberger observer. 2.3. Model identification and validation The model identification is performed using a set of industrial databases. The crystallization conversion rate parameter acryst is estimated to make the model fit the data. An iterative fitting procedure has been implemented. The algorithm is based on the minimization of a quadratic criterion J expressing the mean squared difference between simulated and experimental crystal content cc. In this case, the optimization problem can be described by the following expression: ! 1X i 2 i min J ¼ min ðccmeasure  ccmodel Þ acryst 2 i Using 5 databases this optimization procedure on acryst leads to the following values:

acryst [s1]

Database 1

Database 2

Database 3

Database 4

Database 5

1.2333  1004

1.2087  1004

1.1224  1004

1.2090  1004

1.1534  1004

Fig. 3. Simulated compared with industrial measurement (dotted line) for cc, ms+c and mvap (database 7).

3. Observer design 3.1. Brief recall about observers Let us consider the following nonlinear system:

With a standard deviation of 4.5791  1004 s1 it can be noticed that the value of the crystallization conversion rate is almost the same for the 5 databases. That could be explained by the fact that, during the crystal growth phase, experimental conditions (temperature, concentration, supersaturation range) are nearly the same from one to another. Therefore, in the following, the crystallization conversion rate acryst will be assumed to be constant and equal to the mean value 1.1854  1004 s1. The experimental validation step is performed on cc, cs+c and cvap using industrial data. Prior to each run, the initial conditions for the state variables are determined from off-line measurements for Brix, purity and level at time t = 0 s. Figs. 2 and 3 show simulated versus measured results for cc, ms+c and mvap for database 6 and 7, respectively. Comparisons between simulated and industrial measurements show that the model estimates with a good accuracy the crystal content, the entire mass of sucrose in solution and the mass of emitted vapor. However, let us remind that the crystallization conversion rate is taken as a constant and this assumption is true only in the appropriate temperature (60 8C < Tmg < 70 8C), concentration (2.5 < concentration < 4) and supersaturation range (1 < s < 1.35).

[(Fig._2)TD$IG]

dx ¼ f ðxðtÞ; uðtÞÞ dt y ¼ gðxðtÞÞ

(8)

where x 2 Rn represents the state vector, u 2 Rm the input vector (or manipulated vector) and y 2 Rp the output (or observation) vector. An observer for this system is the following auxiliary dynamical system: dz ¼ f ðzðtÞ; uðtÞ; yðtÞÞ dt ˆ xðtÞ ¼ gðzðtÞ; uðtÞ; yðtÞÞ

(9)

where the input vector is constituted of the input and output ˆ vectors of the system, and the output vector xðtÞ is the estimated state vector. A necessary condition to assign the dynamics of the observer is that the system must be observable, that is, it does not have couples of indistinguishable initial states [30,31]. In the case of non linear systems, the system (8) is locally observable in x if and only if: 9 1 0 8 g > > > > > = C B @ < Lfg >  C B rankB (10)  C¼n .. @@x > > A > . > > > : n1 ; L g f

x

Lfg represents the Lie derivative [32]. 2 3 f  6 f 1 7 @g @g @g 6 2 7 L f g ¼ rg f ¼ ... 6 7 @x1 @x2 @x3 4 ... 5 fn If condition (10) is verified 8x 2 Rn, the system is uniformly observable. 3.2. Observability analysis

3 cc Let us consider the measured outputs vector Y ¼ 4 msþc 5 and mVa p the corresponding nonlinear output function 2 3 x 2

Fig. 2. Simulated compared with industrial measurement (dotted line) for cc, ms+c and mvap (database 6).

6 7 gðxÞ ¼ 4 x1 þ x2x þþx3x þ mi ðtÞ 5 1

2

x4

2

712

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

The observer gain K(t) is calculated at each time t according to the eigenvalues of:   @g  A0 ¼ A  KðtÞ  @x

The necessary and sufficient condition for the system to be uniformly observable is to fill the following condition 8x 2 Rn: 91 0 8  g > > > > > =C B @ C B rankB C ¼ 4 L2f g > A @@x > > > > > : L3 g ; f

t

where A0 is the system matrix of the observer dynamics. Since the system is observable, the matrix K(t) can be chosen in such a way that the observer dynamics matrix A0 has any arbitrarily assigned poles in the left half of the complex plane. However, for convergence of the estimator error, it is mandatory that the matrix A0 is a convergent one. In particular, K(t) can be chosen so that the error converges very quickly to zero. However, the observer becomes very sensitive to the disturbances (for instance, noises of measurement) [37,38]. In this study, a pole placement technique is applied in order to obtain a good trade-off between convergence rate and robustness of the observer.

x

Let us note OM the observability matrix: 9 8  g > > > > > > < @ L f g = OM ¼  ; L2f g > @x > > > > ; : L3 g > f x

it can be demonstrated that the observability matrix is full rank

8x 2 Rn only using the six first lines. Proof. Let us note OM1 the 4  4 matrix composed of line 1, 2, 3 and 6 of the observability matrix.

3.4. Observer validation

A sufficient and necessary condition of the full rank of the observability matrix OM is that OM1 is a full rank matrix: 0

x2 B ðx þ x þ x þ m ðtÞÞ2 B 1 2 3 i OM1 ¼ B 1 B @ 0 k5 k6

x1 þ x3 þ mi ðtÞ

The experimental validation of the extended Luenberger x2

2

ðx1 þ x2 þ x3 þ mi ðtÞÞ 1 0 0

2

ðx1 þ x2 þ x3 þ mi ðtÞÞ 0 0 0

OM1 is a full rank matrix if and only if det(OM1) 6¼ 0, with: x2 detðOM1Þ ¼ k5 k6 ðx1 þ x2 þ x3 þ mi ðtÞÞ2 Besides, 8t, t 2 R+*, (x1, x2, x3, x4, mi(t)) 2 R+*  R+*  R+*  R+*  R+* and k5k6 6¼ 0 So, 8x 2 Rn, det(OM1) 6¼ 0 In other words, OM1 is a full rank matrix 8x 2 Rn that involves that the rank of OM is equal to 4 8x 2 Rn. In conclusion, the system is uniformly observable and so the dynamics of the observer can be assigned. 3.3. Extended Luenberger observer design Different types of state estimators have been developed for nonlinear systems, based on extended Kalman filters estimators (EKF) [33,34], or using extended Lunberger observer (ELO) [35,36]. The main difficulty using the EKF approach is the appropriate choice of the process covariance matrix [12] whereas the ELO requires a not trivial pole placement. In this paper, partially due to the weak nonlinearity of the model, the ELO approach is proposed. Let us consider the following nonlinear model: 2 3 2 8 k1 acryst 0 0 0 > > > 6 7 6 > 0 a 0 0 0 cryst > 6 7 6 > xþ4 > > x˙ ¼ 4 k6 k5 k2 ð1 þ k4 Þ þ k6 ðr f h f  k3 Þ 0 0 05 > < k2 k4 ð1 þ k4 Þ þ k6 ðk3  r f h f Þ k6 k5 0 0 0 2 3 x2 > > > > > 6 7 x þ x þ x þ m ðtÞ > 1 2 3 i > y ¼ gðxÞ ¼ 4 5 > x1 þ x2 > : x4

1 0C C 0C C 1A 0

observer is performed on cc, mvap and ms+c using industrial databases. Figs. 4 and 5 show estimated versus measured values of cc, ms+c and mvap for databases 8 and 9, respectively. It can be notice that for the database 9 the extended Luenberger observer is voluntarily initialized far from the experimental initial conditions (initial conditions are determined at t = 0 s using off-line measurements). In fact, to analyze the convergence rate of the extended Luenberger observer the initial conditions (I.C) are set as follow:

Experimental I.C Observer I.C

ˆ with K(t) is the Luenberger observer gain and xðtÞ is the estimate state vector.

mc [kg]

mw [kg]

mvap [kg]

6130 3000

4660 3000

0 100

Fig. 6 is a zoom of Fig. 5 that shows the convergence rate of the estimated values to the measured for mvap. The observer estimate shows good convergence properties and estimated with a sufficient accuracy the crystal content, the mass of sucrose in solution and the mass of emitted vapor to reach our control objectives. 4. NMPC design The performance and stability of process control depend both

0 0 rw ðð1 þ k4 Þ  hw k6 Þ rw ðk4 ð1 þ k4 Þ þ hw k6 Þ

An extended Luenberger observer is designed, based on the linearization of the output function around the trajectory at each time t to determine the gain of the observer: ˆ˙ ˆ xðtÞ ¼ AxðtÞ þ BuðtÞ þ KðtÞ½yðtÞ  gðxðtÞÞ

ms [kg] 18,000 15,000

3 0 2 3 u1 0 7 74 u 2 5 k6 5 u3 k6

on past and present changes, sometimes due to irregular operating conditions: the reduction of available heating steam, the variation of the quality of feeding juices and molasses, etc. If PID controllers have proved to be suitable in a wide range of operating conditions, advanced strategies may be required to meet steadiness and quality demands, in spite of strong irregularities.

[(Fig._4)TD$IG]

[(Fig._7)TD$IG]

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

713

Fig. 7. NMPC control scheme.

Fig. 4. Estimated compared with industrial measurement (dotted line) for cc, ms+c and mvap (database 8).

[(Fig._5)TD$IG]

(described in Section 4.2) and a nonlinear optimizer for minimizing the mean squared difference between predicted outputs and target values. The aim of the optimization problem is to find the optimal solution of a nonlinear cost function with respect to some predefined performance criterion [39,40]. In the present case, the optimization problem can be described by the following equation: 0 1 Ny NX u 1 X 2 2 ˆ þ i=kÞÞ þ min JðkÞ ¼ min@ g ðsR ðk þ iÞ  yðk b Du ðk þ iÞA i

i¼N1

Fig. 5. Estimated compared with industrial measurement (dotted line) for cc, ms+c and mvap (database 9).

In the present study, a NMPC strategy is proposed to overcome the limitation of PID controllers. It comprises a predictive model of the process and an optimization problem.

i

i¼0

with, SR: target value defined by the reference trajectory to make ˆ predicted process output the process output fit the setpoint w; y: (mass of crystals in solution); u: controlled variable (mass flow rate of the feed); Du(k + i) = u(k + i)  u(k + i  1); g and b: weighting parameters. The prediction horizon Ny corresponds to the future time interval used to compute predictions with the predictor model. The control horizon Nu corresponds to the time interval when present and future control actions are computed. The optimization parameter N1 determines, together with the prediction horizon Ny the coincidence horizon. Identification leading to Nu < Ny, after Nu iterations the controller repeats the last computed signal u(k + Nu) over the prediction horizon Ny. Fig. 8 illustrates an example of the predictive and control horizons for predictive control strategy and Fig. 9 summarizes the predictive control principle. 4.2. Internal predictive model design and validation Process modeling is a key stage in a predictive control strategy. The [(Fig._8)TD$IG] model must be able to predict the process behavior several

4.1. NMPC control scheme The control scheme shown in Fig. 7 comprises the NMPC controller and the process. The adopted NMPC strategy comprises [(Fig._6)TD$IG] an artificial neural network model to predict process outputs

Fig. 6. Convergence of estimated to industrial measurement (cross) for mvap (database 9).

Fig. 8. Prediction and control horizons.

[(Fig._9)TD$IG]

714

[(Fig._10)TD$IG]

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

Fig. 10. ANN closed-loop validation – predicted (bold) compared with estimated (sampling time fixed to 30 s).

ahead prediction (from the beginning to the end), with an error about 5%. Those results are considered satisfactory for achieving NMPC objectives. 5. Control of the industrial crystallization process Fig. 9. Model-based predictive control scheme.

steps ahead, over a future horizon. And it must be simple enough to make sure the associate optimization problem can be solved in reasonable computation time according to the sampling time. In the following, an artificial neural network (ANN) architecture is considered. The model input is the mass flow rate of the feed (u) ˆ and the predicted output is the mass of crystals in solution ðyÞ, which leads to a SISO topology: nu = 1 and ny = 1. To identify a multistep ahead predictor, yˆ is sent back on the input cell, which leads to a closed-loop ANN. Once the structure is established, the following choice that has to be made is the number of past input– output signals used as regressors, i.e. the model order. Thus the regression vector at each instant k is defined as follow: T ˆ  1Þ . . . yðk ˆ  na Þ; uðk  1  nk Þ . . . uðk  nb  nk Þ ’ðkÞ ¼ ½yðk

with: nk: time delay, fixed to 1 in the following; na: number of past estimated outputs; nb: number of past inputs. The closed-loop topology comprises an input layer with nunb + nyna cells, nc hidden neurons with hyperbolic tangent function as activation function and an output layer with 1 linear neuron. The selection of model order is done using the Lipschitz criterion, introduced by He and Asada [41]. Numbers of past inputs and estimated outputs are both fixed to the value four (na = nb = 4). The identification methodology leads to an ANN model with 4 hidden neurons. The identification of the ANN predictor is performed using the estimated mc. The estimate crystal mass is obtained using the extended Luenberger observer presented in Section 3 applied to industrial data. To be closer to the industrial context, a white Gaussian noise with a signal to noise ratio of 1% is added to the estimated mc for both identification and validation. The validation of the closed-loop ANN predictor is done using two additional databases. Fig. 10 shows the predicted versus estimated (by observer) values for mc with a sampling time of thirty second. The closed-loop topology allows a 200 samples

In this section, the benefits of NMPC strategy based on a new controlled variable are presented. Applicability of the NMPC strategy is tested via simulation, using industrial data. Simulation results demonstrate good performances in setpoint tracking even in the presence of simulated vacuum accidents in the data (Section 5.1). Besides, comparisons between simulation results and industrial data (from a PID-controlled process) show the benefits of the new strategy in term of energy and time saving (Section 5.2). In the following, the predictive horizon is fixed to 10, and the control horizon is fixed to 2, according to Clarke’s criterion [42,43]. The sampling time is fixed to 30 s. 5.1. Control performance This subsection considers two different simulation tests to demonstrate the control performance in setpoint tracking (test 1) and in disturbance rejection in setpoint tracking (test 2). The test 1 has two main objectives. The first one is the validation of the predictive control regarding the setpoint tracking (the setpoint is chosen as a linear function), whereas the second one is the determination of an optimal setpoint in order to improve the energy efficiency. This optimal setpoint is chosen in such a way that following this trajectory ensures to be in the metastable zone. Fig. 11 shows the setpoint tracking performance of the NMPC scheme with a slope of 230 kg min1. Series of simulation runs have to be conducted to determine the optimal slope according to the initial conditions and the process kinetics. This optimal slope is determined using a try-error method. The slope is changed to an upper value test after test until a limit value is reached. Over this limit value, the setpoint tracking could not be guaranteed, due to process kinetics, manipulated variable limitations and/or experimental conditions. The application of this iterative procedure to experimental data leads to a mean value of 240 kg min1 for the optimal slope. Test 2 is conducted to demonstrate the robustness of the process control strategy. In this purpose, the response to disturbances in setpoint tracking is examined.

[(Fig._1)TD$IG]

[(Fig._13)TD$IG]

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

Fig. 11. Setpoint tracking (slope of 230 kg min1).

First, several vacuum accidents are simulated during the growth phase: 5 min long vacuum accidents randomly repeated at different time instants. Each time, the reduced pressure inside the crystallizer is increased up to atmospheric pressure (1 bar) instead of 200 mbar. Second, to be closer to experimental conditions, a white Gaussian noise with a signal to noise ratio of 1% is added to the simulated mc. Fig. 12 shows the setpoint tracking with a slope of 240 kg min1 in presence of disturbances. A vacuum accident is simulated from time t = 29 min to time t = 34 min. In industrial context, the electrical conductivity traditionally used as controlled variable is very sensitive to variations in reduced pressure in the pan. In case of vacuum accidents, this sensitivity significantly decreases the process control performance. This result illustrates the robustness of the proposed NMPC strategy and the interest of mc as controlled variable to improve the process control.

715

Fig. 13. NMPC strategy versus PID controller.

[(Fig._14)TD$IG]

Fig. 14. Simulated concentration for controlled mc shown in Fig. 13.

Besides, using mc instead of the electrical conductivity as controlled variable enhances the controller robustness and ensures to be in the appropriate supersaturated state (Fig. 14).

5.2. Benefits of the NMPC strategy proposed 6. Conclusion and prospects A third test is conducted to show the benefits of using NMPC instead of PID controller, and mc instead of the electrical conductivity as controlled variable. In this purpose, a comparison is made between experimental mc obtained with the PID controller and the simulated mc obtained with NMPC controller. The optimal setpoint, with a slope of 240 kg min1, is used for NMPC simulation. Considering the same experimental conditions and the same heating power, Fig. 13 shows the efficiency of the NMPC strategy in terms of process time saving. Indeed, using the same volume of liquor supply to reach almost the same final mass of crystals in solution, the NMPC strategy leads to significant time saving, that involves energy saving: the simulated mc reaches the final value of 24,400 kg after about 63 min. According to the measured mc, about 71 min were necessary for the PID-controlled process to reach this final value.

[(Fig._12)TD$IG]

In this paper, a nonlinear predictive control strategy has been developed to improve the process supervision and control of an industrial process. In a first stage, an alternative controlled variable, the crystal mass, has been proposed to replace the widely used electrical conductivity of the solution. Due to the unavailability of on-line crystal mass measurements, the proposed control strategy is based on an extended Luenberger observer. In a second stage, a more suitable control algorithm, an NMPC approach, has been developed to overcome the limitation of the classical PID controller. The control scheme comprises an artificial neural networks predictive controller and an optimization problem. In terms of process control, simulation results show that the proposed control strategy leads to good performance in setpoint tracking and disturbance rejection. Besides, the robustness of the control strategy has been tested and revealed that the mass of crystals is less sensitive to vacuum accidents than the electrical conductivity. In terms of global efficiency improvement, the proposed control strategy leads to a significant time saving involving a significant energy saving. Acknowledgments

Fig. 12. Setpoint tracking (slope of 240 kg min1).

This study is the fruit of an industrial partnership between BoisRouge sugar mill, eRcane research institute and University of La Reunion. The authors gratefully acknowledge the help of Mr. JeanClaude Pony, Dir. of the BR sugar mill (www.bois-rouge.fr), and its technical staff, and Mr. Jean-Paul Dijoux, networks engineer (www.ercane.re).

716

C. Damour et al. / Journal of Industrial and Engineering Chemistry 16 (2010) 708–716

References [1] T. Bakir, S. Othman, G. Fevotte, H. Hammouri, AIChE J. 52 (2006) 2188. [2] S. Motz, S. Mannal, E.D. Gilles, J. Process Control 18 (2008) 361. [3] A. Mesbah, A.N. Kalbasenka, A.E.M. Huesman, H.J.M. Kramer, P.M.J. Van den Hof, in: Proceedings of the 17th World Congress IFAC, 2008, p. 3246. [4] U. Vollmer, J. Raisch, Control Eng. Pract. 9 (2001) 837. [5] N. Moldova´nyi, B. Lakatos, F. Szeifert, J. Crystal Growth 275 (2005) 1349. [6] Q. Hu, S. Rohani, D.X. Wang, A. Jutan, Powder Technol. 156 (2005) 170. [7] P. Georgieva, S. Feyo de Azevedo, Comput. Intell. 3 (2006) 224. [8] W. Paengjuntuek, A. Arpornwichanop, P. Kittisupakorn, Chem. Eng. J. 139 (2008) 344. [9] M. Sheikhzadeh, M. Trifkovic, S. Rohani, Chem. Eng. Sci. 63 (2008) 829. [10] Z.K. Nagy, Comput. Chem. Eng. 33 (2009) 1685. [11] Z.K. Nagy, B. Mahn, R. Franke, F. Allgo¨wer, Control Eng. Pract. 15 (2007) 839. [12] Z.K. Nagy, R.D. Braatz, AIChE J. 49 (2003) 1776. [13] W. Paengjuntuek, P. Kittisupakorn, A. Arpornwichanop, J. Ind. Eng. Chem. 14 (2008) 442. [14] S. Feyo de Azevedo, J. Chora˜o, M.J. Gonc¸alves, L. Bento, Int. Sugar J. 95 (1993) 483. [15] S. Feyo de Azevedo, J. Chora˜o, M.J. Gonc¸alves, L. Bento, Int. Sugar J. 96 (1994) 18. [16] D. Devogelaere, M. Rijckaert, O.G. Leon, G.C. Lemus, VIIth Brazilian Symp. on Neural Networks, 2002, p. 2. [17] A. Simoglou, P. Georgieva, E.B. Martin, A.J. Morris, S. Feyo de Azevedo, Comput. Chem. Eng. 29 (2005) 1411. [18] H.M. Hulburt, S. Katz, Chem. Eng. Sci. 19 (1964) 555. [19] D. Ramkrishna, Rev. Chem. Eng. 3 (1985) 49. [20] A. Mesbah, J. Landlust, A.E.M. Huesman, H.J.M. Kramer, P.J. Jansens, P.M.J. Van den Hof, Chem. Eng. Res. Design (2009). [21] E. Hugo, Sucrecrie de canne, third ed., Tech & Doc Lavoisier, 1979. [22] P.W. van der Poel, H. Schiweck, T. Schwartz, Sugar Technology: Beet and Cane Sugar Manufacture, Verlag Dr. Albert Bartens KG, Berlin, 1998. [23] B. Grondin-Perez, M. Benne, C. Bonnecaze, J.-P. Chabriat, J. Food Eng. 66 (2005) 361.

[24] B. Grondin-Perez, M. Benne, J.-P. Chabriat, J. Food Eng. 76 (2006) 639. [25] S. Beyou, Ph.D. Thesis, Proposition d’un algorithme PID a` parame`tres variables pour l’ame´lioration de la conduite d’un proce´de´ de cristallisation industriel, University of La Re´union, France, 2008. [26] A. Arpornwichanop, P. Kittisupakorn, Y. Patcharavorachot, I.M. Mujtaba, J. Ind. Eng. Chem. 14 (2007) 175. [27] W. Weerachaipichasgul, P. Kittisupakorn, A. Saengchan, K. Konakom, I.M. Mujtaba, J. Ind. Eng. Chem. 16 (2010) 305. [28] C. Pautrat, J. Ge´notelle, M. Mathlouthi, Sucrose crystal growth: effect of supersaturation, size and macromolecular impurities, in Sucrose crystallization, science and technology, VanHook et al., Bartens, 1997, p. 57. [29] S. Barth, Ph.D. thesis, Utilization of FBRM in the Control of CSD in a Batch Cooled Crystallizer, Georgia Institute of Technology, 2006. [30] G. Bornard, F. Celle-Couenne, G. Gilles, in: A.J. Fossard, D. Normand-Cyrot (Eds.), Observabilite´ et observateurs, Masson, Inc., Paris, 1993, p. 177. [31] E. Sontag, Mathematical Control Theory, Springer Verlag, 1990. [32] R. Marino, P. Tomei, Nonlinear Control Design, Prentice Hall, 1995. [33] P. Bogaerts, Bioprocess Eng. 20 (1999) 249. [34] M. Boutayeb, D. Aubry, IEEE Trans. Autom. Control 44 (1999) 1550. [35] D.G. Luenberger, IEEE Trans. Autom. Control 16 (1971) 596. [36] E. Quintero-Ma´rmol, W.L. Luyben, C. Georgakis, Ind. Eng. Chem. Res. 30 (1991) 1870. [37] V. Sundarapandian, Math. Comput. Model. 43 (2006) 42. [38] C. Pichardo-Almarza, A. Rahmani, G. Dauphin-Tanguy, Simul. Model. Pract. Theory 13 (2005) 179. [39] J. Richalet, Techniques de l’inge´nieur. Informatique industrielle, ISSN 1632-3831 S2: R7423.1-R7423.17 (1997). [40] J. Richalet, A. Rault, J.L. Testud, J. Papon, Automatica 14 (1978) 413. [41] X. He, H. Asada, ACC proc. 2 (1993) 2520. [42] D.W. Clarke, C. Mohtadi, P.S. Tuffs, Automatica 23 (1987) 137. [43] D.W. Clarke, C. Mohtadi, P.S. Tuffs, Automatica 23 (1987) 149.