Article pubs.acs.org/IECR
Controllability of Semibatch Nonisothermal Antisolvent Crystallization Processes G. Cogoni,† S. Tronci,† R. Baratti,*,† and Jose A. Romagnoli‡ †
Dipartimento di Ingegneria Meccanica, Chimica e dei Materiali, Università di Cagliari, Cagliari, Italy Department of Chemical Engineering, Louisiana State University, Baton Rouge, Louisiana 70803, United States
‡
ABSTRACT: In this article the problems of input multiplicities arising in the nonisothermal crystallization process is first considered. Using controllability tools it was found that the system, along a certain trajectory is ill conditioned. Furthermore, the singularity loci of the gain matrix, calculated at asymptotic condition, were obtained. The ill-conditioned nature of the corresponding control problem may imply limitations and affect the responsiveness of the system, and thus the design of a proper control strategy is also addressed in this paper. The proposed control strategy consists of two-stages. In the first stage a feedforward control is used, where the antisolvent flow rate and temperature trajectories are the results of a static optimization procedure. At the second stage, a feedback control action is introduced, and it is used to eliminate the possible offsets due to modeling errors. Results are illustrated for the simulated NaCl−water−ethanol antisolvent crystallization system. the case of linear growth term13 and in asymptotic CSD for nonlinear growth term.14 In order to use the model over the whole operational range, Cogoni et al.15 has recently proposed a relationships for describing the dependence of model parameters from the two process variables, that are antisolvent flow rate and temperature. This procedure led to discover that the same asymptotic distribution, in terms of the first two moments, can be obtained at different operating conditions, thus, implying an input multiplicities of the system, which is consistent with the opposite effects that antisolvent flow rate and temperature have on the crystallization process. The presence of input multiplicities should be properly taken into account when designing the process controller, since the use of standard feedback controller may lead to undesired system behavior.16,17 For this reason a comprehensive analysis of the input-output relationships of the crystallization system and its controllability has been taken into account to design a proper control strategy. The controller has been inspired by the two-stages approach proposed by Liptack et al.18 and it is based on the following actions: (i) the open-loop constrained optimization problem is solved to minimize the mismatch between desired and calculated mean and variance trajectories, (ii) the trajectories of the manipulated variables are implemented to drive the system close to targets, and (iii) a SISO feedback control, based on the linear FPE, is then used to eventually reduce the offset of the mean size of crystals. The choice of applying feedback control, only for the mean (excluding the variance), mainly depends on two characteristics of the crystallization system. The first issue is that mean and variance of the CSD are not independent15 and for a given mean value there exist a region on feasible variance values.
1. INTRODUCTION Antisolvent-aided crystallization is an advantageous separation technique when the solute is highly soluble or heat sensitive. The driving force in crystal formation is the supersaturation that establishes the thermodynamic equilibrium for the solid− liquid separation. The main technologies used to obtain the supersaturation are cooling of the solution and antisolvent addition, and a proper combination of them could improve the quality of the product, as recently demonstrated for few systems.1,2 The organic systems, paracetamol and acetylsalicylic acid, used in the aforementioned two articles have solubility that change significantly with temperature, therefore, to incorporate cooling with antisolvent crystallization can significantly increase crystallization yield. It has also been shown recently3 that even for systems with solubility weakly dependent on temperature; it is possible to impart significantly improved control over both the distribution mean size and coefficient of variation by manipulating temperature together with antisolvent feed rate. This strategy would allow us to add a second degree of freedom (cooling) to be used to control the crystallization process. The control of the crystal size and the crystal size distribution is an important and challenging problem and several factors can affect the size and the widening of the size distribution. The development of effective mathematical models describing the crystal growth dynamics is, therefore, a crucial issue toward finding the optimal process performance and to control the characteristics of the product. Antisolvent crystallization has been modeled for many systems using the traditional population balance modeling approach.4−9 As an alternative, it was recently shown10−13 that it is possible to describe a crystallization process by means of a stochastic approach, which allows the obtainment of the crystal size distribution (CSD) evolution with respect to time using the Fokker−Planck equation. Furthermore, this type of representation leads to simpler models and allows for the analytical solution to describe the CSD over time in antisolvent crystallization operations in © 2014 American Chemical Society
Received: Revised: Accepted: Published: 7056
November 26, 2013 April 3, 2014 April 5, 2014 April 6, 2014 dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
Article
Therefore, if the chosen set-points violate this constraint, no controller system is able to eliminate the offsets on one of the controlled outputs. Furthermore, the estimation of the mean is more efficient with respect to the variance, as recently demonstrated in a previous paper of the authors.14 The article is organized as follows. The Stochastic Modeling Approach section is a short description of the stochastic approach for crystal size distribution, followed by the global models with the direct dependencies with the manipulated variables such as antisolvent flow rate and temperature. In the Input-output Controllability Analysis, a detailed study is described on the relationships between manipulated and controlled variable along with input multiplicities on the CSDs. The section Control Strategy describes an original procedure to control the semibatch crystallization process, based on a two-stage control strategy and finally Conclusions.
ψ (y , t0) =
(1)
In eq 1 g(L;θ) represents the deterministic contribution to the growth rate of L, θ is the vector parameter defined in the model, Lη(t) is the random component linearly depending on crystal size (Geometric Brownian Motion), and η(t) is the Langevin force. It is further assumed that E[η(t )] = 0 E[η(t )η(t ′)] = 2Λδ(t − t ′)
(2)
⎛ dμ(t ) μ(t ) ⎞ = f1 (t ; q , T ) = r2⎜1 − ⎟ with μ(t0) = μ0 dt K2 ⎠ ⎝
where Λ is the intensity of the Langevin force, which has been assumed as constant and equal to √D, with D being the diffusivity parameter. Dividing both terms of eq 1 by L, the Langevin equation becomes dy = g ( y ; θ ) + η( t ) dt
(6)
where μ0 and σ0 are, respectively, the experimental initial mean size of crystals (logarithmic scale) and the standard deviation at t = t0. The two growth terms chosen to describe the crystallization system (eqs 4a,4b) can be regarded as the phenomenological model, which provide the main features of a generic growth process, tending toward an equilibrium value. The first law, eq 4a, was proposed in Grosso et al.,11,12 and it allowed a very good description of the experimental data. The second eq (4a,4bb) was later introduced in Cogoni et al.13 as a simplification of the growth model. The use of a linear term allows the analytical solution of the FPE (eq 5), that can be more easily applied for model-based control applications. In the present paper, the FPE with the logistic term (eq 4a) has been used to simulate the crystallization process. On the other hand, the FPE with the linear growth model (eq 4a) has been assumed to be the available model of the process and has been used for the optimization problem, the derivation of the controller and the selection of the desired final values (mean and variance) of the CSD. When the eq4a is considered, the analytical solution of eq 5 is not available; therefore the CSD is obtained by numerical integration. Conversely, when the growth term is described by the linear expression 4b, the PDF/CSD will preserve its Gaussian shape at any time,14 thus the distribution is described by its first two moments, the mean, μ, and the variance, σ2. The first moment then follows the deterministic Gompertz equation in logarithmic scale (eq 4b), and therefore, its dynamic and analytical solutions are given by
2. STOCHASTIC MODELING APPROACH The crystallization process has been modeled using the approaches proposed previously,11−13 where the growth of crystal size L is assumed to be governed by a deterministic and a stochastic term: dL = Lg (L ; θ ) + Lη(t ) dt
⎡ (y − μ )2 ⎤ 1 0 ⎥ exp⎢ − ⎢⎣ σ0 2π 2σ02 ⎥⎦
⎞ ⎛ r μ(t ) = (μ0 − K 2) exp⎜ − 2 (t − t0)⎟ + K 2 ⎠ ⎝ K2
(3)
(7a,7b)
where y is defined as y = ln L and g(y;θ) can assume the following structures:
The dynamics of the second moment is described by the following differential equation,20 followed by its analytical solution:
⎛ y⎞ g (y ; θ1) = r1y⎜1 − ⎟ K1 ⎠ ⎝ ⎛ y ⎞ g (y ; θ2) = r2⎜1 − ⎟ K2 ⎠ ⎝
r dσ 2(t ) = f2 (t ; q , T ) = 2D2 − 2 2 σ 2(t) dt K2 with σ 2(t0) = σ02 (4a,4b)
⎞ ⎛ ⎛ DK ⎞ r DK σ 2(t ) = ⎜σ02 − 2 2 ⎟ exp⎜ −2 2 (t − t0)⎟ + 2 2 r2 ⎠ r2 ⎠ ⎝ K2 ⎝
The probability density function, PDF, related to the y random variable can then be described by the following Fokker−Planck equation,19FPE:
(8a,8b)
In order to use the model over the whole operating condition window, a global model with explicit functionality of the parameter’s vector θi has been developed.15 The relationship between the manipulated variables q and T for the FPE parameters has been chosen with simple polynomial relations where the polynomial degree could be only equal to 1 or 2. The relationships used, equal for both models, are described as follow:
∂ψ (y , t ) ∂ 2ψ (y , t ) ∂ =D − [g (y , t ; θi)ψ (y , t )] 2 ∂t ∂y ∂y t ≥ 0, y ∈ 9, with i = 1, 2
(5)
where g(y;θi) can assume one of the two expressions (eqs 4a,4b) and the initial condition for t = t0 is assumed as a normal distribution: 7057
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
Article
deep analysis on the input−output relationships is of paramount importance. Preliminary information can be obtained by considering the derivative of the process outputs (mean and variance) with respect to the manipulated variables (q and T)
ri(q , T ) = γi ,0r + γi ,1rq + γi ,2rT K i(q , T ) = γi ,0K + γi ,1K q2 + γi ,2K T 2
i = 1, 2
Di(q , T ) = γi ,0D + γi ,1Dq2 + γi ,2D T
(9)
The coefficients of expressions of eq 9 are reported in Table 1, where it is possible to notice that the behavior of the parameters estimated for every operating condition is maintained for both proposed models.
Kp = ij
nonlinear FPE (i = 1)
linear FPE (i = 2)
γi,0r γi,1r γi,2r γi,0K γi,1K γi,2K γi,0D γi,1D γi,2D
0.5264 0.5983 −6.4588 × 10−4 4.9176 −0.0238 1.7139 × 10−4 0.2134 0.0277 −0.0019
5.7639 2.8342 −0.1584 4.8593 −0.0244 2.0 × 10−4 0.3864 0.0287 −0.0094
∂mj
,
i , j = 1, 2
(y1 , y2 ) = (μy , σy2)
Table 1. Values of the Model Parameters Describing the Dependence of (r, K, D) on Antisolvent Flow Rate and Temperature coefficient
∂yi
(m1 , m2) = (q , T )
(10)
The Kpij variables would be equivalent to the process gain if the system were linear and time independent, but being the crystallization a nonlinear process conducted in semibatch mode, they give punctual information on the effects of inputs on output response depending on the prescribed input trajectories. Different simulations have been run in order to have a broader analysis on system behavior, considering both the linear and nonlinear drift models. For sake of brevity only the results obtained in a specific case are shown, in the understanding that similar behavior can be obtained at different conditions. The time evolution of the four gains obtained with input trajectories reported in Figure 1 are represented in Figure 2, where it is evident that mean and variance gains with respect to q (Kp11 and Kp21) change sign as the system evolves, according to the behavior of system exhibiting input multiplicities. This behavior may indicate that the input−output relationships degenerate at some conditions, and the determinant of the process matrix Kp(t) (eq 11) may approach zero.17
3. INPUT−OUTPUT CONTROLLABILITY ANALYSIS In recent works it was shown that, representing the asymptotic mean and variance on an antisolvent (q) vs temperature (T) plane, input multiplicities were obtained15,21 for a certain window of operating conditions. This means that it is possible to obtain an asymptotic CSD, characterized by its mean and variance, by using two different sets of input values. This particular behavior has been corroborated experimentally,15 and it is due to the presence of competing effects of the process inputs; for example, high antisolvent flow rate and low temperature favor crystal growth rate, while low flow rate and high temperature lead to higher asymptotic mean crystal size. From a control perspective, input multiplicities imply that the steady-state gains of the system may change sign as the system moves from one operating point to another; therefore, if a conventional feedback control is applied, it can lead to undesired output values, unstable or oscillatory responses.16 With the objective of designing a robust and efficient control for the crystallizer it is important to eliminate the difficulties related to the occurrence of input multiplicities; therefore, a
⎛ ∂μ ∂μ ⎞ y y ⎛ K p (t ) K p (t )⎞ ⎜⎜ ∂q ∂T ⎟⎟ 11 12 ⎟= K p(t ) = ⎜ ⎜ K (t ) K (t )⎟ ⎜⎜ ∂σ 2 ∂σ 2 ⎟⎟ p22 y y ⎝ p21 ⎠ ⎜ ⎟ ∂ ∂ q T ⎝ ⎠t
(11)
As a rigorous framework for a more reliable indication of the conditioning of the process, the minimum singular value and condition number of Kp matrix have been also calculated. The SVD has been applied after a proper scaling of the variables, using the following relationships (eq 12):
Figure 1. Antisolvent flow rate (left panel) and temperature (right panel) trajectory. 7058
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
Article
Figure 2. Time evolution of the derivatives of process output with respect to the manipulated variables (q: left panel, T: right panel), calculated for the trajectory reported in Figure 1
Figure 3. Time evolution of the minimum singular value (left panel) and condition number of the process matrix Kp(t) with respect to the manipulated variables (q and T), calculated for the trajectory reported in Figure 1
⎧ ∂μ ̅ ⎪ y ⎪ ∂q ̅ ⎪ ⎪ ∂μy̅ ⎪ ⎪ ⎪ ∂T̅ ⎨ ⎪ ∂σ ̅y2 ⎪ ⎪ ∂q ̅ ⎪ ⎪ ∂σy̅ 2 ⎪ ⎪ ∂T̅ ⎩
=
∂μy (q(t ) − q(t )) 0 f ∂q (μy (t 0) − μy (tf ))
=
∂μy (T (t ) − T (t )) 0 f ∂T (μy (t0) − μy (tf ))
indicating sensitivity to model errors and possible difficulties for achieving acceptable control performance. It has been also verified that the occurrence of the maximum value of the condition number (or lesser value of the minimum singular value) corresponds to zero determinant for the given trajectory. At this condition the input−output relationship degenerates, and controlling the system is practically impossible. The difficulty of designing an effective feedback control strategy has been also assessed by analyzing first the rank of the Kalman matrix22 (eq 13) and, since the system presents nonlinearities, the determinant of the Gramian matrix23 (eq 15) for both asymptotic and dynamic conditions. In order to analyze the controllability of the system, eqs 7a,7b−8a,8b have been linearized with respect to the manipulated variables q and T.
∂σy2 (q(t 0) − q(tf )) = ∂q (σy2(t0) − σy2(tf )) =
∂σy2 (T (t0) − T (tf )) ∂T (σy2(t0) − σy2(tf ))
(12)
where t0 and tf are, respectively, the initial and the final sampling time. The obtained values are reported as a function of time in Figure 3, using the input trajectory shown in Figure 1. The system along the entire trajectory is constantly illconditioned with the condition number always greater than 10,
K = [ BAB]|t = t *
(13)
Where the matrixes A and B are defined as follows: 7059
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research ⎡ ∂f ∂f1 ⎤ ⎡ ∂f ⎢ 1 ⎥ ⎢ 1 2 ⎢ ∂μy ∂σy ⎥ ⎢ ∂q ⎥ and B = ⎢ A=⎢ ⎢ ∂f2 ∂f2 ⎥ ⎢ ∂f2 ⎢ 2⎥ ⎢⎣ ∂q ⎢⎣ ∂μy ∂σy ⎥⎦ G (t ) =
∫t
∂f1 ⎤ ⎥ ∂T ⎥ ⎥ ∂f2 ⎥ ∂T ⎥⎦
Article
(14a,14b)
t
exp( At ) BB′exp(A′ t) dt 0
(15)
Performing the controllability tests for different input trajectories, it is possible to verify that the rank of eq 13 is always equal to 2, but the controllability Gramian (eq 15) is always ill-conditioned. The Gramian condition numbers in fact vary from tens (milder ill-conditioning) to thousands (stronger ill-conditioning) along the prescribed system trajectories. This result implies that, in principle, it should be possible to control CSD mean and variance, but the ill-conditioned nature of the present control problem may imply limitations and affect the responsiveness of the system, according to the result obtained with analyzing the Kp matrix. System Behavior at Asymptotic Conditions. In the definition of the final product, in terms of mean and variance of CSD, it is important to define the region of input multiplicities at the asymptotic conditions and then avoid the targets close to the ill-conditioned region. When considering the linear model (eqs 7a,7b and 8a,8b), it is possible to use an explicit relationship between the outputs and the manipulated variables: ⎧ μy ,asymp (q , T ) = K (q , T ) ⎪ ⎪ ⎨ D(q , T )K (q , T ) ⎪ σy2,asymp(q , T ) = ⎪ r (q , T ) ⎩
Figure 4. Degeneracy locus of the nonlinear FPE (thick solid black line) and linear FPE (dashed black line) along with the iso-level curves representing respectively the variances (solid red lines) and the mean sizes (solid blue lines), in linear scale; a) FPE with eq 4b as drift term; b) FPE with eq 4a as drift term. +∞ ⎧ y Ψasymp(y ; q , T ) dy ⎪ μy ,asymp (q , T ) = −∞ ⎪ ⎪ 2 ⎨ σy ,asymp(q , T ) +∞ ⎪ = (y − μy ,asymp (q , T ))2 Ψasymp ⎪ −∞ ⎪ (y ; q , T ) d y ⎩
∫
∫
(16a,16b)
The singularity loci of the gain matrix Kp calculated at asymptotic conditions can be obtained substituting the eqs 16a,16b in eq 11 and finding the values of q and T leading to zero determinant. The obtained values are reported in Figure 4a (dotted black line), along with the iso-level curves representing respectively the variances (solid red lines) and the mean sizes (solid blue lines) in linear scale in case of linear growth term. From a control point of view, where the gain array becomes singular, there is no direct guidance on the loop pairing between the manipulated and the controlled variables, and a conventional control system may be inadequate if the desired asymptotic conditions are lying on (or close to) the singularity loci. It is worth noticing that zero determinant conditions imply also minimum variance for a given mean crystal size; therefore, obtaining a narrow CSD manipulating q and T in a feedback control framework may be a demanding task. The same procedure to individuate the singularity at asymptotic conditions can be applied to the global model obtained through the nonlinear FPE, in which the drift term is defined as a logistic model (eq 4a). The definition of the system in asymptotic conditions is not straightforward as for the linear drift term. The system obtained using the nonlinear FPE is defined as an implicit system of integral equations with respect to the manipulated variables:
(17a,17b)
where the asymptotic CSD distribution is given by the following relationship:14 ⎡ r(q , T ) exp⎢ D(q , T ) ⎣ Ψasymp(y ; q , T ) = ⎡ r(q , T ) +∞ ∫−∞ exp⎢⎣ D(q , T )
y2 2
y3 3K (q , T )
( − ( − y2 2
)⎤⎥⎦ )⎤⎥⎦ dy
y3 3K (q , T )
(18)
Given a desired mean size, it is possible to solve iteratively eqs 17a,17b using eq 18, within the operability window and then obtain all the possible variances for the mean size of crystals considered. According to the results previously obtained, the region with the minimum variance for a given mean size of crystals corresponds to the frontier line where the sign of the asymptotic gain matrix changes and, thus, where the system is ill-conditioned. Repeating this approach for the range of mean sizes included in the operability window and, for each obtaining the minimum variance, it is possible to define numerically the zero determinant conditions (for the nonlinear global model) shown in Figure 4b (solid black line). Analyzing the two operative maps for the linear (Figure 4a) and nonlinear (Figure 4b) models, it is possible to observe that temperature has a small effect on asymptotic variance, 7060
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
Article
Figure 5. Upper panel: Mean size and variance behavior (continuous line) and desired set-point (dotted line) when applying the optimal trajectory; Bottom panel: optimal antisolvent flow rate trajectory (left) and optimal feasible temperature trajectory (right).
operational challenges discussed by Bonvin25 and Juba and Hamer,26 the most remarkable is the Time-varying characteristics, which means that the transition in the controlled variables may be large compared to typical excursions for continuous systems, implying that in case a standard linear transfer function model (gain and time constant) is used, it may need to be timevarying and, batch characteristics can change from run-to-run i.e.: growth velocity (r), diffusivity constant (D), and asymptotic dimension of crystals (K). It was demonstrated by Liptak18 and Bonvin23 that a combination of an on−off controller and a PID controller can faster achieve the desired conditions for a batch (or semibatch) crystallizer. This combination was made up for a heated batch reactor, where the maximum heating is applied until just before or within ±5% of the desired set-point; then the on−off controller is switched off and the PID controller activated to bring the controlled variable smoothly to the desired value. This strategy could be used for the considered crystallization problem, with some modification aimed to eliminate the difficulties encountered by the presence of input multiplicities, ill-conditioning, and possible violation of output constraints, leading to selection of unfeasible set-points. In this framework a two-stage control procedure for the crystallizer has been developed. In the first stage feedforward control is used, where the antisolvent flow rate and temperature trajectories are the results of a constrained optimization procedure, aimed to minimize the area between the desired mean size and variance of crystals and the calculated outputs. In order to mimic a real situation, where a mismatch generally exists between the process model and the plant, the linear model (eq 4b) is used as nominal model for calculating the optimal trajectories and select the desired output values by means of the asymptotic operative map (Figure 4a), while the nonlinear model (eq 4a) is used to simulate the crystallization process.
confirming results obtained with the analysis of the controllability matrix. Furthermore, comparison of the two maps evidences that the interval of feasibility variance values for a given mean is different for the two models. This means that the use of a nominal model with uncertainties to select the setpoints for the controlled outputs may lead to unfeasible conditions and then may determine oscillating response, close to the singularity loci, or may cause the manipulated variables to reach restriction values (saturation). As an illustrative example the asymptotic conditions (μ,σ2) = (126.7 μm, 3429.6 μm2) have been indicated in Figure 4a with a black circle. By inspection of Figure 4b, it is possible to observe that the two iso-lines for the asymptotic conditions selected using the linear model do not intersect for q and T values belonging to the operative region. A significant difference is also detected for the singularity loci reported in a and b of Figure 4, where it is possible to observe that the region obtained for the nonlinear FPE is wider with respect to using the linear FPE. Again, an uncertain estimation of the ill-conditioned region may lead to undesired behavior of the controller, if process conditions have been selected not sufficiently far from the actual critical zone.16,17 The above observations evidence the importance of a good characterization of the system.
4. CONTROL STRATEGY The development of an efficient and robust control for the crystallization system should take into account the following issues: (i) occurrence of input multiplicities, (ii) ill-conditioning of the input−output relationship, and (iii) semibatch conduction mode. The first two aspects evidence that a conventional feedback control could be inadequate for the problem at hand, while batch or semibatch processes are challenging to track the set-point since there are not steadystate operating points and there are wide operating ranges,24 as for the case studied for antisolvent crystallization. Among the 7061
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
Article
system over the entire operating window, a 2 × 2 feedback control structure has demonstrated its inability to work properly for each possible process condition. The proposed solution is a two-stage controller, where a feedforward action is used to bring the system close to the desired conditions and a feedback action is then switched on to eliminate the offset of the mean crystal size due to modeling errors. The choice of renouncing to a strict control on variance and accepting the offset resulting from the application of the feedforward action is related to the high uncertainty on variance measurement (or estimation). It has been assessed that the mismatch between actual and desired variance using the feedforward control is mostly within the experimental error; therefore, the adjustment with feedback controller is not applied, avoiding problems related to selection of unfeasible variance with the assumed linear model and system ill-conditioning. The choice of the manipulated variable for the feedback control on the mean crystal size is based on the input effects on the controlled variable, in terms of gain and response time. From the analysis of the response, the antisolvent flow rate has been shown to better fit the requirements. The design of the feedback controller has been based on the Internal Model Controller, IMC,27 approach, using the linear FPE as the process model. In this case, the mean crystal size dynamics is ruled by eqs 7a,7ba, which are nonlinear with respect to the input variables. In order to obtain the control law, a linear approximation of the mean dynamics is obtained28 by considering a small deviation (δμ, δq) from a reference trajectory (μ*(t), q*(t)):
The second stage of the developed control strategy is based on a feedback control action, and it is used to eliminate the possible offsets due to modeling errors. The control, based on an IMC algorithm, is switched when the mean size has reached a certain percentage within its set-point. Because of poor conditioning of the input−output relationships evidenced in the previous section, and the presence of uncertainty in the nominal model, which may lead to select unfeasible desired conditions, a single loop is used for achieving the desired mean crystal size. The variance is only feed-forwardly controlled, because its measure is not as robust as the mean,14 and it eventually will present a slight offset due to model errors. Optimal Input Trajectories. The optimal trajectories are obtained by minimizing the residual area between the desired mean (μy,d) and the variance (σ2yd) values of the CSD and their time evolution trajectories obtained by the FPE with the linear growth term. The targets cannot be chosen independently, and they have been determined through the operative map (black circle in Figure 4a). The optimization problem was subjected to a maximum and a minimum constraint for the temperature and antisolvent flow rate and a maximum constraint to exclude the region above the singularity locus (Figure 4a). +∞ tf ⎧ |μy , d − y Ψ(y , t ; q , T ) dy|dt ⎪ min −∞ ⎪ q(t ), T(t ) t 0 ⎪ +∞ tf ⎪ |σy2, d − min (y − μy , d )2 Ψ −∞ ⎪ q(t ), T(t ) t0 ⎪ (y , t ; q , T ) dy|dt ⎨ ⎪ subject to: ⎪ ⎪ 0.7 mL/min ≤q ≤ 3.0 mL/min ⎪ 10°C ≤ T ≤ 30°C ⎪ ⎪ T < g (q)|ill − conditioning frontier ⎩
∫
∫
∫
∫
μ(t ) = δμ + μ*(t) q(t ) = δq + q*(t )
(20a,20b)
Applying eqs 20a,20b to the system eq 7a and performing a Laplace transform, the transfer function of the process in a canonical form becomes
(19)
The final time (tf) has been chosen sufficiently large to ensure that mean and variance of the crystal size distribution will not undergo further changes (asymptotic conditions). The resulting optimal input trajectories are step functions, meaning that the minimum residual area can be obtained if the inputs are directly moved to the values leading to the desired asymptotic mean and variance. Considering the maximum heat exchangeable by the cooling/heating bath, the temperature trajectory has been rearranged from a step response to a ramp transition between the initial and the final temperature. The behavior of the mean and standard deviation of the CSD for the given suboptimal trajectory is reported in Figure 5 (upper panel) along with the input trajectories (bottom panel). The presence of an offset for both mean and standard deviation is due to modeling errors; in fact the manipulated inputs for the given desired conditions have been calculated using the linear model (eq 4b), while the plant has been simulated using the nonlinear growth term (eq 4a). It is worth noticing that the greatest difference between the two models is for the variance of the CSD, where the offset is equal to 80 μm2. The mismatch for the mean calculated with the linear and the nonlinear model is ∼1.7 μm, indicating that the linear model gives a good approximation for the prediction of the first moment of the CSD, but the use of the nonlinear growth term allows for a better description of variance.13 Two-Stage Controller. According to the time-varying behavior of the semibatch system and ill-conditioning of the
Gp(s) =
K p(q*, T *) τ(q*, T *)s + 1
(21)
where ⎡⎛ μ ⎞ ∂r rμ ∂K ⎤ + 2 K p(q , T ) = τ(q*, T *)⎢⎜1 − ⎟ ⎥ K ⎠ ∂q K ∂q ⎦ ⎣⎝ τp(q , T ) =
μ*, q*, T *
K (q*, T *) r(q*, T *) (22)
The transfer function eq 21 can be used to design IMC, getting to Gq (s) =
τ(q*, T *)s + 1 K p(q*, T *)[εs + 1]
(23)
where ε is the IMC filter (tuning parameter). The parameters τp and Kp depend on the state of the system (eq 22); therefore, there are two possibilities: (i) the controller gain (Kc = τp/ Kp/ε) and integral time (τi = τp) are updated at each sampling time obtaining an adaptive controller, (ii) or they are kept constant. In the latter case, it is important to properly select the parameter values that guarantee the desired controller performance even if the system characteristics change with time. For the system considered, the best performance has 7062
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
Article
Figure 6. Upper panel: Response of mean (left) and variance (right) for the two-stage controller. Bottom panel: antisolvent flow rate (right) and temperature (left) trajectory.
Figure 7. Set-points selected with the perfect model. Upper panel: Response of mean (left) and variance (right) for the two-stage controller. Bottom panel: antisolvent flow rate (right) and temperature (left) trajectory.
been obtained by keeping constant the integral time and equal to the characteristic time of the process evaluated at the switching time and adjusted by a proportional factor of 1.5, in order to avoid oscillations on the process response. The controller gain has been updated each sampling time, where the IMC filter has been defined proportionally to the characteristic time of the process of a factor of 1.5. The controller parameters are summarized as follows:
Kc(q , T ) =
1 αK p(q , T )
τi = βτ(tswitch)
(24)
where α and β are both equal to 1.5. The switching from the feedforward controller to the SISO feedback controller for the mean size has been realized, 7063
dx.doi.org/10.1021/ie404003j | Ind. Eng. Chem. Res. 2014, 53, 7056−7065
Industrial & Engineering Chemistry Research
■
considering the absolute value of the first derivative of the mean size of crystals with respect to the time. The value of the switch condition has been considered to be less than the tuning parameter ξ, which has been set equal to 0.05: dμy dt
Article
AUTHOR INFORMATION
Corresponding Author
*Telephone: +39 070 6755056. Fax: +39 070 6755067. E-mail:
[email protected]. Notes
