Model Predictive Control Strategies for Batch Sugar Crystallization Process

1. Introduction

The industrial processes are governed generally by general principles of the physics and chemistry. With the aid of data acquisition systems supported in microprocessor it is possible to obtain real data of the industrial process, that it characterizes in detail his dynamics and input-output dependency. Several methods of identification allow, from these data, to obtain linear and nonlinear models of these processes (Rossiter, 2003; Morari, 1994); which are the base to predict the process behaviour within all the family of the model based predictive controllers (MPC).

Diverse algorithms MPC have demonstrated its effectiveness in those control loops characterized by strong nonlinearities, difficult dynamic, inverse answers and great delay; that they are generally those of greater influence in the final product quality and the process efficiency (Allgöwer et al., 2004; Qin & Badgwell 2003).

One of the most important steps in the implementation of a MPC is just the obtaining of the model that can predict with reliability the future behaviour of the controlled variable, like answer to a predefined optimized control action (Rawlings 2000). This work applies two kind of MPC: (i) Classical Model-Based Predictive Control and (ii) Neural Network Model Predictive Control (NNMPC).

The classical MPC strategy uses a discrete model obtained from general phenomenological model of the feed-batch crystallization process, consisting of mass, energy and population balance. The NNMPC strategy uses to obtain a neural network, the training algorithms proposed in the Neural Network Toolbox of MatLab (version 7.04) (Bemporad et al., 2005).

In this particular case it is analyzed a fed-batch sugar crystallization process, in this process there is abundant information, detailed mathematical models and real industrial data. (Chorão, 1995; Feyo de Azevedo & Gonçalves 1988; Georgieva et al., 2003). This fact motivated the use of the neural networks to model the process and to propose a neural network MPC (NNMPC) that considers the process like a gray box, of which has input-output information and the historical experience of he process behaviour.

2. Batch sugar crystallization process

2.1. General description

The operation of crystallization is applied in the sugar industry to obtain the sucrose dissolved in the extracted juice of the sugar cane or the sugar beet basically.

Typical industrial fed-batch evaporative sugar crystallization is performed in a vacuum pan crystallizer. The reactor has a cylindrical form with volume that can vary between 20-60 m3. The feed system is usually equipped with an extra water input to dilute the sugar solution if necessary. The heat transfer system is a calandria type, to permit the heat interchange between steam and suspension. The vacuum pressure in the pan is generated by the contact barometric condenser and the pan is equipped with a mechanical agitator to keep the suspension homogeneous. The operation is conducted in a fed-batch mode with an average duration of a cycle about 90 minutes.

Sugar crystallization occurs through the mechanisms of nucleation, growth and agglomeration. In the course of production, the crystallization phenomenon is driven by two mechanisms (Jancic & Grootscholten, 1984): i) mass transfer from dissolved sucrose to crystal surface and ii) heat transfer in the calandria. Shortly before the grain setting and continuing during the beginning of the crystallization phase, the available crystalline surface to deposit the molecule of sucrose is much smaller than the mass of dissolved sucrose. During this period the evaporation rate is high, the crystal area/mass of crystallized sucrose rate is very low, therefore the process is driven by the mass transfer. The supersaturation tends to increase and if not controlled, it often achieves the undesirable zone of secondary crystal nucleation. Later on, when the total crystal area and the crystallization capacity increases, the crystal area/mass of crystallized sucrose rate gets high and the process is driven by the heat transfer.

The process objectives are to maximize the speed of crystal growth, keeping high the produced sugar quality and minimizing the costs and losses. These objectives must be fulfilled without occurrence of secondary nucleation or agglomeration. The sugar quality is evaluated by the particle size distribution (PSD) at the end of the process which is quantified by two parameters - the final average (in mass) particle size (MA) and the final coefficient of particle variation (CV). The main challenge of the sugar production is the large batch to batch variation of the final PSD. This lack of process repeatability is caused mainly by improper control policy and results in product recycling and loss increase. The sugar production is heuristically operated, and while the traditionally applied PI(D) controllers are still the preferred solutions they usually lead to energy and material loss that can easily be reduced if an optimized operation policy is implemented. These problems constitute the main motivation for the operation strategy formulated in the next section.

2.2. Crystallization model

The general phenomenological model of the fed-batch crystallization process consists of mass, energy and population balances, including the relevant kinetic rates for nucleation, linear growth and agglomeration (Simoglou et al., 2005). While the mass and energy balances are common expressions in many chemical process models, the population balance is related with the crystallization phenomenon, which is still an open modelling problem. The Appendix A shows a detailed phenomenological model for crystallization process.

2.3. Problem formulation

The final values of the crystal size distribution function (CSD) parameters: mass averaged crystal size (MA) and coefficient of variation (CV) are the best indicators of the quality and efficiency of the crystallization process. The direct measurement and control of these parameters are very difficult to make actually, in fact there are no references of its industrial implementation. The most used solution in the sugar industry consists of establishing a strategy that manipulate other variables; which allows to arrive at the end of the process with acceptable values in the CSD parameters.

The batch operation imposes to the process frequent operational changes that depend of: the quality of the raw material, disturbances in the work conditions and market demand changes. The previous problem, the nonlinearities and the restrictions imposed to the process motivated the use a nonlinear MPC (NMPC).

When a NMPC algorithm is applied, the first challenge consists of obtaining the model to use, which must be viable and trustworthy. Although the sugar crystallization process has been studied in depth and efficient mathematical models exist to represent it, these must be validated and be fit before its application in a NMPC algorithm, which will cause frequent updates if the process is batch.

Like an alternative, in this work it is tried to demonstrate the efficiency that has the use of the neuronal networks in a NMPC, where the neural networks could be trained from industrial data with the input-output answer of the process.

3. Problem solution

Sugar production is characterized by strongly non-linear and non-stationary dynamics and goes naturally through a sequence of relatively independent stages: charging, concentration, seeding, setting the grain, crystallization (the main phase), tightening and discharge (Georgieva et al., 2003). Therefore the operation strategy is formulated as a cascade of individual control loops for each of the stages (Fig. 1). The feedback control policy is based on measurements of the flowrate, the temperature, the pressure, the stirrer power and the supersaturation (by a refractometer). Measurements of these variables are usually available for a conventional crystallizer.

3.1. Operation strategy

Sugar production is still a very heuristically operated process, with classical proportional integral and eventually derivative (PID) controllers being the most typical solution. The different phases of the sugar production are comparatively independent and moved by distinct driving forces, thus a single controller can hardly be effective for the complete process. Instead, individual controllers for each stage where it seems appropriate, was the adopted framework (Fig. 1). See Table 1 for more details on the formulated operation strategy.

In the present study, the control actions are performed by manipulating the valves of the liquor/syrup feed flowrates (Ff) and the steam flowrate (Fs), while the volume of massecuite (Vm), the supersaturation (S) and the current of the agitator (IA) are the controlled variables. This choice is completely inspired by the industrial practice in several refineries.

Charging (stage 1): During the first stage the crystallizer is fed with liquor until it covers approximately 40 % of the vessel height. The process starts with vacuum pressure of around 1 bar (equal to the atmospheric pressure) and reduces it up to 0.23 bar. When the vacuum pressure reaches 0.5 bar, the feed valve is completely open such that the feed flowrate is kept at its maximum value. When the liquor covers 40 % of the vessel height, the feed valve is closed and the vacuum pressure needs some time to stabilize around the value of 0.23 bar before the concentration stage starts.

Figure 1.

Cascade MPC control - strategy

Concentration (stage 2): Once the vacuum pressure stabilizes, the stirrer is switched on and the concentration begins. In order to guarantee unperturbed operation of the barometric condenser and the steam production boiler, the steam flowrate must increase slowly (from 0 to 2 kg/s, in two minutes approximately). The concentration of the dissolved sucrose by evaporation under vacuum results in volume reduction. However, for technological reasons, the minimum suspension level of the pan must be above the calandria. Therefore a feed flowrate action is required to control the level (the volume) of the pan and this constitutes the first control loop. In this stage, the supersaturation increases rapidly (at about a rate of 0.025 per min.). When it reaches a value of 1.06, the feeding is stopped and the steam flowrate is reduced slowly to 1.4 kg/s, with the same speed as it was increased. The concentration stage is over when the supersaturation reaches the value of 1.11.

Seeding (stage 3): At this moment seed crystals are introduced into the pan to provoke crystallization. This stage is rather unstable and to prevent seed crystals from dissolution in the liquor, the feed valve must be closed and the steam flowrate kept at its minimum for a short period (about 2 min.). Keeping these conditions unchanged contributes to the formation of the grain and is also important for the final crystal size distribution. The supersaturation continues naturally to increase but usually no control action is required.

Crystallization with liquor (stage 4): During this stage the supersaturation is first controlled by a proper feeding to be around a set point of 1.15. This constitutes the second control loop. At the beginning of this stage, the mass transfer is the driving crystallization force, the crystallization rate increases and the controller usually reduces the feed flowrate

Stage	Action	Control
Charge	The steam valve is closed and the stirrer is off.; The vacuum pressure changes from 1 to 0.23 bar.; The vacuum pressure reaches 0.5 bar, feeding starts with max rate.; Liquor covers 40 % of the vessel height.	No control; The feed valve is completely open
Concentration	The vacuum pressure stabilizes around 0.23 bar.; The stirrer is on.; The volume is kept constant.; The steam flowrate increases to 2 kg/s; The supersaturation reaches 1.06, the feeding is closed, the steam flowrate is reduced to 1.4 kg/s	Control loop 1; Controlled variable: Volume; Manipulated variable: liquor feed flowrate
Seeding and setting the grain	The supersaturation reaches 1.11.; Seed crystals are introduced.; The steam flowrate is kept at the minimum for two minutes.	No control; The feed valve is closed
Crystallization with liquor; (phase 1)	The steam flowrate is kept around 1.4 kg/s.; The supersaturation is controlled at the set point 1.15.	Control loop 2 Controlled variable: supersaturation Manipulated variable: liquor feed flowrate
Crystallization with liquor; (phase 2)	The volume of crystallizer reaches ≈ 22 m^3.; The feed valve is closed.; The supersaturation is controlled at the set point 1.15.; The stirrer power reaches 20.5 A.	Control loop 3 Controlled variable: supersaturation Manipulated variable: steam flowrate
Crystallization with syrup	The steam flowrate is kept around the maximum of 2.75 kg/s. (hard constraint).; The volume fraction of crystals is kept at the set point 0.45.; The volume reaches its maximum value (30 m^3); The feed valve is close.	Control loop 4 Controlled variable: volume fraction of crystals. Manipulated variable: syrup feed flowrate
Tightening	The stirrer power reaches the maximum value of 50 A (hard constraint).; The steam valve is closed.; The stirrer and the barometric condenser are stopped.	No control

Table 1.

Summary of the sugar crystallization operation strategy.

to maintain the reference value of the supersaturation. When all liquor quantity is introduced, the feeding is stopped and the supersaturation is now kept at the same set point of 1.15 by the steam flowrate as the manipulated variable. This constitutes the third control loop. The heat transfer is now the driving crystallization force. A typical problem of this control loop is that at the end of this stage the steam flowrate achieves its maximum value of 2.75 kg/s but it is not sufficient to keep the supersaturation at the same reference value therefore a reduction of the set point is required. The stage is over when the stirrer power reaches the value 20.5 A.

Crystallization with syrup (stage 5): A stirrer power of 20.5A corresponds to a volume fraction of crystals equal to 0.4. At this moment the feed valve is reopened, but now a juice with less purity (termed syrup) is introduced into the pan until the maximum volume (30 m³) is reached. The control objective is to maintain the volume fraction of crystals around the set point of 0.45 by a proper syrup feeding. This constitutes the fourth control loop.

Tightening (stage 6): Once the pan is full the feeding is closed. The tightening stage consists principally in waiting until the suspension reaches the reference consistency, which corresponds to a volume fraction of crystals equal to 0.5. The supersaturation is not a controlled variable at this stage because due to the current conditions in the crystallizer, the crystallization rate is high and it prevents the supersaturation of going out of the metastable zone. The stage is over when the stirrer power reaches the maximum value of 50 A. The steam valve is closed, the water pump of the barometric condenser and the stirrer are turned off. Now the suspension is ready to be unloaded and centrifuged.

4. Model based predictive control

The term model-based predictive control (MPC) does not refer to a particular control method, instead it corresponds to a general control approach (Rossiter, 2003). The MPC concept, introduced in late seventies, nowadays has evolved to a mature level and became an attractive control strategy implemented in a variety of process industries (Camacho & Bordons, 2004). The main difference between the MPC configurations is the model used to predict the future behavior of the process or the implemented optimization procedure. First the MPC based on linear models gained popularity (Morari, 1994) as an industrial alternative to the classical proportional-integral-derivative (PID) control and later on nonlinear cases as reactive distillation columns (Balasubramhanya & Doyle, 2000) and polymerization reactors (Seki et al., 2001) were reported as successfully MPC controlled processes.

4.1. Classical model based predictive control

The main difference between MPC configurations is the model used to predict the future behaviour of the process and the optimization procedure. Nonlinear model predictive control (NMPC) is an optimisation-based multivariable constrained control technique that uses a nonlinear dynamic model for the prediction of the process outputs (Qin & Badgwell, 2003). At each sampling time k the model predicts future process responses to potential control signals over the prediction horizon (H_p). The predictions are supplied to an optimization procedure, to determine the values of the control action over a specified control horizon (H_c) that minimizes the following performance index:

$$\underset{u_{\min }\le \left[u_c(k),u_c(k+1)\mathrm{,.....}{u}_c(H_c)\right]\le u_{\max }}{\min }J={\lambda }_1{\displaystyle \sum _{k=1}^{H_p}{\left(y_r(k)-\widehat{y}(k)\right)}^2}-{\lambda }_2{\displaystyle \sum _{k=1}^{H_c}{\left(u_c(k-1)-u_c(k-2)\right)}^2}$$E1

Subject to the following constrains

$$u_{\min }\le u_c\le u_{\max }$$E2

$$\Delta u_{\min }\le \Delta u\le \Delta u_{\max }$$E3

$$y_{\min }\le y_p\le y_{\max }$$E4

Where $$u_{\min }$$and $$u_{\max }$$ are the limits of the control inputs, $$\Delta u_{\min }$$and $$\Delta u_{\max }$$ are the minimum and the maximum values of the rate-of-change of the inputs and $$y_{\min }$$and $$y_{\max }$$ are the minimum and maximum values of the process outputs.

H_p is the number of time steps over which the prediction errors are minimized and the control horizon H_c is the number of time steps over which the control increments are minimized, $$y_r$$is the desired response (the reference) and $$\widehat{y}$$ is the predicted process output (Diehl et al., 2002). $$u_c^{}(k),u_c(k+1),u_c(H_c)$$are tentative future values of the control input, which are parameterized as peace wise constant. The length of the prediction horizon is crucial for achieving tracking and stability. For small values of H_p the tracking deteriorates but for high H_p values the bang-bang behavior of the process input may be a real problem. The MPC controller requires a significant amount of on-line computation, since the optimization (1) is performed at each sample time to compute the optimal control input. At each step only the first control action is implemented to the process, the prediction horizon is shifted or shrunk by usually one sampling time into the future, and the previous steps are repeated (Rossiter, 2003). $${\lambda }_1$$and $${\lambda }_2$$are the output and the input weights respectively, which determine the contribution of each of the components of the performance index (1).

4.2. Neural network model predictive control

The need for neural networks arises when dealing with non-linear systems for which the linear controllers and models do not satisfy. Two main achievements contributed to the increasing popularity of the NNs: (i) The proof of their universal approximation properties and the development of suitable algorithms for NN training as the backpropagation and (ii) The adaptation of the Levenberg-Marquard algorithm for NN optimization.

The most used NN structures are Feedforward networks (FFNN) and Recurrent (RNN) ones. The RNNs offer a better suited tool for nonlinear system modelling and is implemented in this work (Fig.2). The Levenberg-Marquard (LM) algorithm was preferred as the training method due to its advantages in terms of execution time and robustness. Since the LM algorithm requires a lot of memory, a powerful (in terms of memory) computer is the main condition for successful training. In order to solve the problem of several local minima, that is typical for all derivative based optimization algorithms (including the LM method), we have repeated several time the optimization specifying different starting points.

The individual stages of the crystallization process are approximated by different RNNs of the type shown in Fig. 2. Tangent sigmoid hyperbolic activation functions are the hidden computational nodes (Layer 1) and a linear function is located at the output (Layer 2). Each NN has two vector inputs (r and p) formed by past values of the process input and the NN output respectively. The architecture of the NN models trained to represent different process stages is summarized as follows:

Figure 2.

Neural network architecture

$$u_{NN}=\left[r,p\right]=\left[u_c(k-1),u_c(k-2),y_{NN}(k-1),y_{NN}(k-2)\right]$$E5

$$x=W_{11}r+W_{12}p+b_1$$E6

$$n_1=\left(e^x-e^{-x}\right)/\left(e^x+e^{-x}\right)$$E7

$$n_2=w_{21}{n}_1+b_2$$E8

Where$$W_{11}\in R^{m\times 2}$$, $$W_{12}\in R^{m\times 2}$$, $$w_{21}\in R^{1\times m}$$, $$b_1\in R^{m\times 1}$$, $$b_2\in R$$are the network weights (in matrix form) to be adjusted during the NN training, m is the number of nodes in the hidden layer.

Since the objective is to study the influence of the NNs on the controller performance, a number of NN models is considered based on different training data sheets.

1. Case 1 (Generated data): Randomly generated bounded inputs ($$u_i$$) are introduced to a simulator of a general evaporative sugar crystallization process introduced in Georgieva et al., 2003. It is a system of nonlinear differential equations for the mass and energy balances with the operation parameters computed based on empirical relations (for no stationary parameters) or keeping constant values (for stationary parameters). The simulator responses are recorded ($$y_{\text{i}}$$) and the respective mean values are computed ($$u_{\text{i}\text{,mean}}$$,$$y_{\text{i}\text{,mean}}$$). Then the NN is trained supplying as inputs $$u_{\text{i}}-u_{\text{i}\text{,mean}}$$ and as target outputs$$y_{\text{i}}-y_{\text{i}\text{,mean}}$$.

2. Case 2: Industrial data: The NN is trained with real industrial data. In order to extract the underlying nonlinear process dynamics a prepossessing of the initial industrial data was performed. From the complete time series corresponding to the input signal of one stage only the portion that really excites the process output of the same stage is extracted. Hence, long periods of constant (steady-state) behavior are discarded. Since, the steady-state periods for normal operation are usually preceded by transient intervals, the data base constructed consists (in average) of 60-70% of transient period data. A number of sub cases are considered. Case 2.1: Industrial data of two batches is used for NN training. Case 2.2: Industrial data of four batches is used for NN training. Case 2.3: Industrial data of six batches is used for NN training.

Figure 3.

Case1: NN data generation

4.3. Selection of MPC parameters: H_p, H_c, ₂

The choice of $$H_p$$ is related with the sampling period ($$\Delta t$$) of the digital control implementation, which in its turn is a function of the settling time t_s (the time before entering into the 5% around the set-point) of the closed loop system. As a rule of thumb, it is suggested $$\Delta t$$ to be chosen at least 10 times smaller than t_s, (Soeterboek, 1992). Hence, the prediction horizon can be chosen as $$H_p$$= round-to-integer(t_s/$$\Delta t$$). It is well known that the smaller the sampling time, the better can a reference trajectory be tracked or a disturbance rejected. However, choosing a small sampling time yields a large prediction horizon. In order to compute the optimal control input, the optimization (1) is performed at each sampling time, therefore MPC controller requires a significant amount of on-line computation. This can cause problems related with large amount of computer memory required and additional numerical problems due to the large prediction horizon. The introduction of the ET MPC as in (7) serves as a compromise between these conflicting issues and reduces significantly the computational efforts.

Parameters $${\lambda }_1$$ and $${\lambda }_2$$ determine the contribution (the weight) of each term of the performance index, the output error (e) and the control increments ($$\Delta u$$). In this work the parameter $${\lambda }_1$$ is set to the normalized value of 1, while the choice of $${\lambda }_2$$ is based on the following empirical expression:

$${\left(u_{\max }-u_{\min }\right)}^2\cdot {\lambda }_2=e_{\max }\cdot P/100$$E9

where P defines the desired contribution of the second term in (1) (0% ≤ P ≤100%) and

$$e_{\max }=\max \left({\left(ref-y_{\max }\right)}^2,{\left(ref-y_{\min }\right)}^2\right)$$E10

The intuition behind (9-10) is to make the two terms of (1) compatible when they are not normalized and to overcome the problem of different numerical ranges for the two terms. Table 2 summarize the set of MPC parameters used in the four control loops define in the section 3.

Control; loop (CL)	ts (s); settling; time	Δt (s); sampling; period	Hp; prediction; horizon	Hc; control; horizon	λt; weight	Controlled; variable	Set-point
CL1	40	4	10	2	1000	Volume	12.15
CL2	40	4	10	2	0.1	Supersaturation	1.15
CL3	60	4	15	2	0.01	Supersaturation	1.15
CL4	80	4	20	2	10000	Fraction; of crystals	0.43

Table 2.

MPC design parameters for the control loops define in Table 1

5. PID controllers

The PID parameters were tuned, where$$k_p$$,$${\tau }_i$$ ,$${\tau }_d$$ are related with the general PID terminology as follows (Aström & Hägglund, 1995):

$$u(t+k)=K_p\left[e(t+k)+\frac{\Delta t}{{\tau }_i}\cdot {\displaystyle \sum _{i=0}^{k}e(t+i)}+\frac{{\tau }_d}{\Delta t}\cdot \left(e(t+k)-e(t+k-1)\right)\right]$$E11

Since the process is nonlinear, classical (linear) tuning procedures were substituted by a numerical optimization of the integral (or sum in the discrete version) of the absolute error (IAE):

$$IAE={\displaystyle \sum _{k=1}^N\left|ref(t+k)-y_p(t+k)\right|}$$E12

Equation (12) was minimized in a closed loop framework between the discrete process model and the PID controller. For each parameter an interval of possible values was defined based on empirical knowledge and the process operator expertise. A number of gradient (Newton-like) optimization methods were employed to compute the final values of each controllers summarized in Table 1. All methods concluded that the derivative part of the controller is not necessary. Hence, PI controllers were analyzed in the next tests.

	Control loop 1	Control loop 2	Control loop 3	Control loop 4
kp	0.05	-0.5	20	-0.01
τi	30	40	10	70
τp	0	0	0	0

Table 3.

Optimized PID parameters for the control loops define in Table 1

6. Discussion of results

The operation strategy, summarized in Table 1 and implemented by a sequence of Classical- MPC, NNMPC or PI controllers is comparatively tested in Matlab environment. The output predictions are provided either by a simplified discrete model (with the main operation parameters kept constant) or by a trained ANN model (5-8). A process simulator was developed based on a detailed phenomenological model (Georgieva et al., 2003). Realistic disturbances and noise are introduced substituting the analytical expressions for the vacuum pressure, brix and temperature of the feed flow, pressure and temperature of the steam with original industrial data (without any preprocessing(Scenario-2)). The test is implemented for two different scenarios of work.

1. Scenario - 1: The simulation uses, like process, the set of equations differentials proposed in (Georgieva et al. 2003) with empirical operation parameters.

2. Scenario - 2: The simulation uses, like process, the set of equations differentials proposed in (Georgieva et al. 2003), but are used like operation parameter e real industrial data batch not used in neural network training.

Time trajectories of the controlled and the manipulated variables for the control loop 1, 2 and 4 of one batch (Batch 1) are depicted in Figs. 4-6. The three controllers guarantee good set point tracking. However, the quality of the produced sugar is evaluated only at the process end by the crystal size distribution (CSD) parameters, namely AM and CV. The results are summarized in Table 4 and both classical and NNPMC outperform the PI. Our general conclusion is that the main benefits of the MPC strategy are with respect to the batch end point performance.

Figure 4.

Controlled (Volume of massecuite) and control variables (F_f- feed flowrate) over time for the 1^st control loop.

Figure 5.

Controlled (Supersaturation) and control variables (F_f- feed flowrate) over time for the 2^nd control loop.

Figure 6.

Controlled (Volume fraction of crystals) and control variables (F_f- feed flowrate) over time for the 4^th control loop.

Performance measures	Classical MPC	NN-MPC	PI
AM (mm) (reference 0.56)	0.586	0.584	0.590
CV (%)	32.17	31.13	32.96

Table 4.

Batch end point performance measures (Batch - 1)

Performance measures	Classical MPC	NN-MPC	PI
AM (mm) (reference 0.56)	0.615	0.609	0.613
CV (%)	29.39	30.28	31.14

Table 5.

Batch end point performance measures (Batch - 2)

Performance measures	Classical MPC	NN-MPC	PI
AM (mm) (reference 0.56)	0.636	0.631	0.639
CV (%)	28.74	29.42	29.23

Table 6.

Batch end point performance measures (Batch - 3)

7. Appendix A. Crystallization model

Sugar crystallization occurs through the mechanisms of nucleation, growth and agglomeration. The general phenomenological model of the fed-batch crystallization process consists of mass, energy and population balances, including the relevant kinetic rates for nucleation, linear growth and agglomeration [Ilchmann, et al., 1994]. While the mass and energy balances are common expressions in many chemical process models, the population balance is related with the crystallization phenomenon, which is still an open modeling problem.

Mass balance

The mass of all participating solid and dissolved substances are included in a set of conservation mass balance equations:

$$\dot{M}=f_1(M(t),F(t),S_1(t)),t_0\le t\le t{}_f,M(0)=M_0$$E13

where $$M(t)\in {\Re }^q$$ and $$F(t)\in {\Re }^m$$ are the mass and the flow rate vectors, with q and m dimensions respectively, and $$t_f$$is the final batch time. $$S_1(t)\in {\Re }^{r1}$$is the vector of physical time dependent parameters as density, viscosity, purity, etc. For the process in hand, the detailed form of the macro-model (13) is as follows

$$M_{sol}=M_a+M_i+M_w$$E14

$$M_m=M_{sol}+M_c$$E15

$$\frac{d{M}_w}{dt}=F_f{\rho }_f\left(1-B_f\right)+F_w{\rho }_w-J_{vap}$$E16

$$\frac{d{M}_i}{dt}=F_f\cdot {\rho }_f\cdot B_f\cdot \left(1-Pu{r}_f\right)$$E17

$$\frac{d{M}_a}{dt}=F_f\cdot {\rho }_f\cdot B_f\cdot Pu{r}_f-J_{cris}$$E18

$$\frac{d{M}_c}{dt}=J_{cris}$$E19

$$V_m=\frac{M_c+M_{sol}}{{\rho }_{sol}}$$E20

$$J_{vap}=\frac{W+Q}{{\lambda }_{vap}}+K_{vap}\cdot \left(T_m-T_{w(vac)}-BPE\right)$$E21

Energy balance

The general energy balance model is

$$\frac{d{T}_m}{dt}=a{J}_{cris}+b{F}_f+c{J}_{vap}+d$$E22

where parameters a, b, c and d incorporate the enthalpy terms and specific heat capacities derived as time dependent functions of physical and thermodynamic properties asfollows

$$a=\frac{H_{sol}-H_c+(1-B_{sol})\frac{d{H}_{sol}}{d{B}_{sol}}+\frac{1-Pu{r}_{sol}}{B_{sol}}\cdot \frac{d{H}_{sol}}{dPu{r}_{sol}}}{M_{sol}\cdot C{p}_{sol}+M_c\cdot C{p}_c}$$E23

$$b=\frac{{\rho }_f\left(H_f-H_{sol}+(B_f-B_{sol})\frac{d{H}_{sol}}{d{B}_{sol}}+\frac{B_f\left(Pu{r}_f-Pu{r}_{sol}\right)}{B_{sol}}\cdot \frac{d{H}_{sol}}{dPu{r}_{sol}}\right)}{M_{sol}\cdot C{p}_{sol}+M_c\cdot C{p}_c}$$E24

$$c=\frac{H_{sol}-H_{vap}-B_{sol}\cdot \frac{d{H}_{sol}}{d{B}_{sol}}}{M_{sol}\cdot C{p}_{sol}+M_c\cdot C{p}_c}$$E25

$$d=\frac{W+Q+F_w{\rho }_w(H_w-H_{sol}+B_{sol})\frac{d{H}_{sol}}{d{B}_{sol}}}{M_{sol}\cdot C{p}_{sol}+M_c\cdot C{p}_c}$$E26

$$\frac{d{H}_{sol}}{d{B}_{sol}}=-29.7T_m+4.6Pu{r}_{sol}{T}_m+0.075T_m^2$$E27

$$\frac{d{H}_{sol}}{dPu{r}_{sol}}=4.61B{x}_{sol}{T}_m$$E28

Population balance

Mathematical representation of the crystallization rate can be achieved through basic mass transfer considerations or by writing a population balance represented by its moment equations. Employing a population balance is generally preferred since it allows to take into account initial experimental distributions and, most significantly, to consider complex mechanisms such as those of size dispersion and/or particle agglomeration/aggregation. The basic moments of the number-volume distribution function are

$$\frac{d{\tilde{\mu }}_0}{dt}={\tilde{B}}_0-\frac{1}{2}\cdot \beta \text{'}\cdot {\tilde{\mu }}_0^2$$E29

$$\frac{d{\tilde{\mu }}_1}{dt}=G_v\cdot {\tilde{\mu }}_0$$E30

$$\frac{d{\tilde{\mu }}_2}{dt}=2\cdot G_v\cdot {\tilde{\mu }}_1+\beta \text{'}\cdot {\tilde{\mu }}_1^2$$E31

$$\frac{d{\tilde{\mu }}_3}{dt}=3\cdot G_v\cdot {\tilde{\mu }}_2+3\cdot \beta \text{'}\cdot {\tilde{\mu }}_2^2$$E32

$$J_{cris}={\rho }_c\cdot \frac{d{\tilde{\mu }}_1}{dt}$$E33

,

where$${\tilde{B}}_0$$, G and $${\beta }^{\text{'}}$$ are the kinetic variables nucleation rate, linear growth rate and the agglomeration kernel, respectively with the following mathematical descriptions

$${\tilde{B}}_0=K_n\cdot 2.894\cdot {10}^{12}\cdot G^{0.51}{\left(\frac{{\tilde{\mu }}_1}{k_v\cdot V_m}\right)}^{0.53}\cdot V_m$$E34

$$\beta \text{'}=\frac{K_{ag}\cdot G\cdot {\tilde{\mu }}_1}{V_m^2}$$E35

$$G=K_g\cdot \exp \left(-\frac{57000}{R\left(T_m+273\right)}\right)\cdot \left(S-1\right)\cdot \exp \left(-13.863\left(1-P_{sol}\right)\right)\cdot \left(1+2\cdot \frac{v}{V_m}\right)$$E36

$$G_v=3\cdot k_v{\left(\frac{v}{{\tilde{\mu }}_0}\right)}^{2/3}\cdot G$$E37

.

The crystallization quality is evaluated by the particle size distribution (PSD) at the end of the process which is quantified by two parameters - the final average (in mass) particle size (AM) and the final coefficient of particle variation (CV) with the following definitions:

$$AM=\overline{L}$$E38

$$CV=\frac{\sigma }{\overline{L}}$$E39

Where $$\sigma$$ and $$\overline{L}$$are computed from:

$$\overline{L}={\left(\frac{{\eta }_3}{1+3\cdot {\left(\frac{\sigma }{\overline{L}}\right)}^2}\right)}^{1/3}$$E40

$$15\cdot {\eta }_3^2\cdot {\left(\frac{\sigma }{\overline{L}}\right)}^6+\left(45\cdot {\eta }_3^2-9\cdot {\eta }_6\right){\left(\frac{\sigma }{\overline{L}}\right)}^4+\left(15\cdot {\eta }_3^2-6\cdot {\eta }_6\right){\left(\frac{\sigma }{\overline{L}}\right)}^2+{\eta }_3^2-{\eta }_6=0$$E41

In (41, 42), $${\eta }_j$$represent moments of mass-size distribution functions, that are related to the moments of the number-volume distribution functions $${\mu }_j$$ by the following relationships:

$${\eta }_3=\frac{{\mu }_2}{k_v\cdot {\mu }_1}$$E42

,

and

$${\eta }_6=\frac{{\mu }_3}{k_v{}^2\cdot {\mu }_1}$$E43

Correlations for physical properties

$$Q={\alpha }_s\cdot F_s\cdot \Delta H_s$$E44

$${\rho }_f=\left(1000+\frac{B{x}_f\cdot \left(200+B{x}_f\right)}{54}\right)\cdot \left(1-0.036\cdot \frac{T_f-20}{160-T_f}\right)$$E45

$$C{p}_f=4186.8-29.7\cdot B{x}_f+4.61\cdot B{x}_f\cdot Pu{r}_f+0.075\cdot B{x}_f\cdot T_f$$E46

$$H_f=C{p}_f\cdot T_f$$E47

$${\rho }^{*}{}_{sol}=\left(1000+\frac{B{x}_{sol}\cdot \left(200+B{x}_{sol}\right)}{54}\right)\cdot \left(1-0.036\cdot \frac{T_m-20}{160-T_m}\right)$$E48

$${\rho }_{sol}={\rho }^{*}{}_{sol}+1000\cdot \left(-1+\exp \left[\left(-6.927\cdot {10}^{-6}\cdot B{x}_{sol}{}^2-1.164\cdot {10}^{-4}\cdot B{x}_{sol}\right)\cdot \left(Pu{r}_{sol}-1\right)\right]\right)$$E49

$$C{p}_{sol}=4186.8-29.7\cdot B{x}_{sol}+4.61\cdot B{x}_{sol}\cdot Pu{r}_{sol}+0.075\cdot B{x}_{sol}\cdot T_m$$E50

$$H_{sol}=C{p}_{sol}\cdot T_m$$E51

$${\rho }_m=\frac{{\rho }_{sol}\cdot {\rho }_c}{{\rho }_c-w_c\cdot \left({\rho }_c-{\rho }_{sol}\right)}$$E52

$$Pu{r}_{sol}=\frac{M_a}{M_a+M_i}$$E53

$$B_{sol}=\frac{M_a+M_i}{M_{sol}}$$E54

$$B{x}_{sol}=100\cdot B_{sol}$$E55

$$B{x}_{sat}=64.447+8.222\cdot {10}^{-2}\cdot T_m+1.66169\cdot {10}^{-3}\cdot T_m{}^2-1.558\cdot {10}^{-6}\cdot T_m{}^3-4.63\cdot {10}^{-8}\cdot T_m{}^4$$E56

$$S^{*}=1.129-0.284\cdot \left(1-Pu{r}_{sol}\right)+\left(2.333-0.0709\cdot \left(T_m-60\right)\right)\cdot {\left(1-Pu{r}_{sol}\right)}^2$$E57

$$S=\frac{\frac{B{x}_{sol}}{100-B{x}_{sol}}}{\frac{B{x}_{sat}}{100-B{x}_{sat}}\cdot C_{sat}}$$E58

$$C_{sat}=0.1\cdot \frac{B{x}_{sol}}{100-B{x}_{sol}}\cdot \left(1-Pu{r}_{sol}\right)+0.4+0.6\cdot \exp \left(-0.24\cdot \frac{B{x}_{sol}}{100-B{x}_{sol}}\cdot \left(1-Pu{r}_{sol}\right)\right)$$E59

$$v=\frac{M_c}{{\rho }_c}$$E60

$$w_c=\frac{M_c}{M_c+M_{sol}}$$E61

$$C{p}_c=1163.2+3.488\cdot T_m$$E62

$$H_c=C{p}_c\cdot T_w$$E63

$${\rho }_w=1016.7-0.57\cdot T_w$$E64

$$T_{w(vac)}=122.551\cdot \exp \left(-0.246\cdot P_{vac}\right)\cdot {\left(P_{vac}\right)}^{0.413}$$E65

$$T_{w(s)}=100.884\cdot \exp \left(-1.203\cdot {10}^{-2}\cdot P_s\right)\cdot {\left(P_s\right)}^{0.288}$$E66

$${\lambda }_{w(vac)}=2263.28-58.21\cdot \ln \left(P_{vac}\right)$$E67

$${\lambda }_s=2257.51-85.95\cdot \ln \left(P_s\right)$$E68

$$H_w=2323.3+4106.7\cdot T_w+T_w{}^2$$E69

$$H_{w(s)}=2323.3+4106.7\cdot T_{w(s)}+T_{w(s)}{}^2$$E70

$$H_s=2491860-13270\cdot P_s+\left(1946.5+37.9\cdot P_s\right)\cdot T_s$$E71

$$H_{vac}=2499980-24186\cdot P_{vac}+\left(1891.1+106.1\cdot P_{vac}\right)\cdot T_m$$E72

$$\Delta H_s=H_s+H_{w(s)}$$E73

$$BPE=\left(0.03-0.018\cdot Pu{r}_{sol}\right)\cdot \left(T_{w(vac)}+84\right)\cdot \left(\frac{B{x}_{sol}}{100-B{x}_{sol}}\right)$$E74

For more detailed presentation of the process model, refer to [Georgieva et al., 2003].

8. Conclusion

With the results obtained in this work it has been demonstrated that algorithm NNMPC is a viable solution to control nonlinear complexes processes, still in the case that only exists input-output information of the process.

An aspect very important to obtain successful results with NNMPC is the representative quality of the available data, which was demonstrated with the results obtained in the third control loop analyzed.

The weighting factor $${\lambda }_2$$ has a crucial paper in the good NNMPC performance. A constrain very hard can impose that the control signal can not follow the dynamics of the process, but a very soft constrain can cause instability in the control signal, when the model is not precise.

Acknowledgments

Several institutions contributed for this study: 1) Foundation of Science and Technology of Portugal, which financed the scholarship of investigation of doctorate SFR/16175/2004; 2) Laboratory for Process, Environmental and Energy Engineering (LEPAE), Department of Chemical Engineering, University of Porto; 3) The Institute of Electronic Engineering and Telematics of Aveiro (IEETA); 4) Sugar refinery RAR, Portugal; The authors are thankful to all of them.