Stochastic optimization of simulated moving bed process

Computers and Chemical Engineering 29 (2005) 1577–1589

Stochastic optimization of simulated moving bed process Sensitivity analysis for isocratic and gradient operation Grzegorz Ziomek, Dorota Antos∗ Faculty of Chemistry, Chemical Engineering Department, Rzesz´ow University of Technology, Poland Received 9 March 2004; received in revised form 12 September 2004; accepted 13 December 2004 Available online 20 January 2005

Abstract A random search strategy has been used for designing and optimization of simulated moving bed process (SMB) under isocratic as well as under solvent gradient conditions. The effectiveness of both the process modes has been compared. For predictions of the objective functions, i.e., the minimum of eluent consumption and/or the maximum of the process productivity a mathematical model of the process dynamics has been employed and implemented in the optimization procedure. Four-dimensional space of decision variables corresponding to the flowrates in the SMB zones has been searched in order to find the optimal set of the process parameters. The optimization was constrained to the purity demand in the outlet streams withdrawn in the raffinate and the extract port. The obtained set of random decision variables fulfilling purity constraints was used to construct the operating window of parameters guaranteeing successful separation. For feasible operating points the sensitivity of the purity constraints with respect to the operating parameters has been calculated. The results of calculations indicated that operating conditions, which ensure similar process efficiency, could correspond to different sensitivity of the process constraints. Such an analysis was found to be useful for the selection of process conditions, for which the best trade-off between the process efficiency and its robustness can be achieved. This appears to be particularly important for designing the gradient SMB process, for which robustness of the operating conditions is a factor of a major importance. In order to improve the efficiency of calculations a modification of the original random search procedure based on the Luus–Jaakola algorithm has been proposed. © 2004 Elsevier Ltd. All rights reserved. Keywords: Stochastic optimization; Simulated moving bed chromatography; Gradient elution

1. Introduction In the last decade the simulated moving bed (SMB) process has been successfully implemented as a separation technique in the petrochemical, biochemical and fine chemical industries. The SMB is a well-established separation technology basing on continuous chromatography process. The SMB unit was designed as a practical realization of a true moving bed (TMB) where the solid and fluid phases move countercurrently. The general concept of a classical four-zone SMB unit is illustrated schematically in Fig. 1. There are two incoming streams: the feed mixture to be separated (V˙ F ) and desorbent or eluent (V˙ D ). Two streams leave the unit, ∗

Corresponding author. Tel.: +48 17 8651730; fax: +48 17 8651730. E-mail address: dorota [email protected] (D. Antos).

0098-1354/$ – see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2004.12.007

one enriched with the less adsorbable raffinate (V˙ R ), and one enriched with the more adsorbable extract (V˙ E ). The four streams divide the unit into four zones (I–IV). Each of these zones contains at least one fixed bed (column) and has to fulfill distinct tasks, i.e., in the zones II and III separations takes place, while in zones I and IV the solid and the fluid phases are regenerated, respectively. The movement of the solid bed is simulated by switching of ports or columns in certain time intervals (Ruthven & Ching, 1989). The SMB technique is implemented in an industrial scale using the same solvent to prepare the feed and to perform the adsorbent regeneration. This so-called isocratic process is presently well understood (Mazzotti, Storti, & Morbidelli, 1994, 1996, 1997; Migliorini, Mazzotti, & Morbidelli, 1998, 1999; Storti, Masi, Carra, Mazzotti, & Morbidelli, 1989).

1578

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

Nomenclature a, b linear coefficients in Eqs. (28)–(30) C concentration in the fluid phase (vol.%) EC eluent consumption (vol. eluent/(vol. product)) H Henry’s constant km lumped mass transport coefficient (min−1 ) Keq equilibrium constant (vol.−1 ) m flowrate ratio m vector of decision variables N total number of iterations in the each pass NB number of fixed beads (columns) NC number of components Nend number of points in the time grid NK number of cells NEQ total number of external passes NGP total number of evaluations NIQ total number of internal passes pH , pKeq parameters in Eqs. (10) and (11) PR productivity (vol. product/(vol. solid phase min)) Pu purity q solid phase concentration (vol.%) q˜ i solid phase concentration in equilibrium with the local fluid phase concentration (vol.%) rH , rKeq parameters in Eqs. (10) and (11) R number of points generated in each iteration S column cross area (cm2 ) t time (min) t* switching time (min) u superficial velocity (cm/s) V column volume (cm3 ) ˙ V volumetric flowrate (cm3 /min) Greek letters β control parameter for reducing the size of search region γ contraction factor for external passes δ parameter defining the size of search region for variable m εt total void fraction τ residence time Subscripts 0 initial values D desorbent E extract f final search region F feed i component index j space interval index mod modifier neq neq th pass of optimization algorithm niq niq th pass of optimization algorithm opt current best value R raffinate

Superscripts K zone index N time interval index

Recently, the idea of modulating of solvent strength in order to increase productivity was introduced to the liquid SMB separation (Abel, Mazzotti, & Morbidelli, 2002, 2004; Antos and Seidel-Morgenstern, 2001, 2002a; Houwig, van Hateren, Billiet, & van der Wielen, 2002; Jensen, Reijns, Billiet, & van der Wielen, 2000). The solvent strength can be modulated in two steps using different solvents at the two inlet ports. The feed is dosed continuously in a relatively weak solvent, whereas as the desorbent a stronger solvent is used (see Fig. 1). Thus, in the classical four zones closed-loop SMB process two distinct levels of internal solvent composition exist which are separated by the two inlet positions. These two characteristic levels of solvent strength can be adjusted by using different amounts of a suitable modifier in the two feed streams. As a result, the components to be separated are more retained in the solvent regeneration zone of the SMB process (zone IV) and more easily eluted in the sorbent regeneration zone (zone I). Recent results of studying this type of twostep gradient SMB process in the closed-loop arrangement demonstrated its potential to reduce significantly the solvent consumption and, thus, to increase product concentrations (Abel et al., 2002, 2004; Antos & Seidel-Morgenstern, 2001, 2002a; Ziomek, Kaspereit, Je˙zowski, & Seidel-Morgenstern, 2005). Design and optimization of such a complex process in an industrial scale can be done by the use of an optimization procedure coupled with an adequate mathematical model of process dynamics. Recently, stochastic optimization algorithms have been successfully applied for the optimization of isocratic SMB process, e.g., a genetic algorithm was implemented in the multiobjective optimization of a reactive SMB process

Fig. 1. Scheme of a four-zone SMB unit allowing the implementation of a gradient operation.

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

(Subramani, Hidajat, & Ray, 2003; Zhang, Hidajat, & Ray, 2002a), SMB and the Varicol process (Zhang, Hidajat, Ray, & Morbidelli, 2002b; Zhang, Mazzotti, & Morbidelli, 2003), a random search algorithm of Luus and Jaakola (1973) was used for the optimization of mobile phase composition in the isocratic and gradient SMB process (Ziomek et al., 2005). In this work the multi-pass Luus and Jaakola (1973), Luus (2001) MPLJ algorithm was modified and utilized to construct the window of the operating parameters for isocratic as well as for gradient SMB process. The effectiveness of both the process modes has been compared. Since the application of the solvent gradient was found to be very promising with respect to the productivity and the eluent consumption, main efforts have been focused on the developing of a reliable numerical tool for the designing of the gradient SMB operation, i.e., selecting conditions favorable for both the process effectiveness and its robustness. The optimization involved four decision variables corresponding to the flowrates in the SMB zones and it was constrained to the purity demand in the outlet streams withdrawn in the raffinate and the extract port. For the operating points generated during the calculations sensitivity functions of the purity with respect to the operating parameters have been calculated. The sensitivity analysis indicated that operating conditions, which guarantee similar process effectiveness can correspond to different values of the sensitivity functions. Low absolute values of sensitivity are related to high robustness, i.e., ability of the unit to be retained in the state, in which it can perform purity demand required. Hence, for the process realization a set of the operating parameters favorable for both effectiveness and robustness should be selected. This appears to be particularly important for the designing of the solvent gradient SMB process, for which small deviations from the operating conditions can results in failure of the separation.

2. Design of SMB process Instructive information concerning possible process performances under isocratic conditions is offered by the socalled triangle theory (Mazzotti et al., 1997). For the ideal situation of an infinite number of equilibrium stages regions of acceptable flowrates leading to pure raffinate and extract streams can be calculated analytically using the equilibrium theory (Mazzotti et al., 1997). The regions are conventionally plotted in the plane (mIII –mII ) of dimensionless flowrates:

mK =

V˙ K t ∗ − Vεt , V (1 − εt )

K = I, . . . , IV

(1)

1579

The switching time t* is defined as t∗ =

V (1 − εt ) V˙ s

(2)

where εt is the total porosity, V˙ K the volumetric flowrate in zone (or fixed bed) K, V˙ s the solid flowrate in the equivalent TMB unit and V the volume of a single fixed bed. If the conditions for the complete regeneration of the solid and the mobile phase are fulfilled, i.e., the lower bound for mI and the upper one for mIV is satisfied, the region of complete separation has a triangular shape limited by the adsorption isotherm. The vertex of the triangle corresponds to the maximal throughput, hence, to the maximal productivity and the minimum solvent consumption (Mazzotti et al., 1994, 1996, 1997; Storti, Mazzotti, Morbidelli, & Carra, 1993). The bounds of the regions calculated by the use of the triangle theory are related to 100% of purity of both the outlet streams at infinite system efficiency. For a real system with a finite column efficiency and reduced purity demands the equilibrium theory overestimates or underestimates the sizes of regions. The size of region diminishes with decreasing the system efficiency and widens with reduction of the purity demands (Biressi, Ludemann-Hombourger, Mazzotti, Nicoud, & Morbidelli, 2000; Migliorini et al., 1999). The accurate position of the region bounds is available only by the numerical solving of a mathematical model of process dynamics. In order to construct the operating window all the four dimensionless flowrates should be considered; a simple, two-parameter screening of the (mII –mIII ) plane is not sufficient to determine the operating conditions. A systematic screening of the four-dimensional space appears to be excessively time consuming. Hence, an efficient optimization routine is necessary for the process designing. Therefore, in this work the attempts have been made in order to select an adequate numerical algorithm and improve its effectiveness. The determining of the operating conditions guaranteeing successful separation is much more difficult for the closedloop gradient SMB process, which involves changes of adsorption equilibria according to the local modifier distribution in the SMB unit. Due to the process complexity and a number of operating parameters involved the optimization of the gradient process should be preceded by the optimization of the operating conditions for the isocratic mode. The effectiveness of both the processes with respect to the process productivity and its robustness should be compared.

3. Equilibrium stage model For the modeling of the SMB process under isocratic as well as gradient conditions the kinetic stage model has been employed (Antos & Seidel-Morgenstern, 2001; Ray, Carr, & Aris, 1994). In this model each fixed bed in the Kth zone of

1580

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

volume V is divided in a cascade of NK cells: K dCi,j

dt K dqi,j

dt

+

K K K − Ci,j Ci,j−1 1 − εt dqi,j = εt dt τK

work by the use of the following finite difference scheme: (3)


+

=

K K km (˜qi,j (CK j ) − qi,j ),

(4) where τ K stands for the liquid phase residence time in a cell: τ

εt V εt V = K K = K K , N u S N V˙

K = I, . . . , IV

(5)

K designates the conIn the above equations t is the time, Ci,j K the solid centration in the fluid phase in a cell j of zone K, qi,j K K phase concentration of component i, q˜ i,j (Cj ) the solid phase concentration at the equilibrium with the vector of fluid phase concentrations Cj , km the lumped mass transport coefficient, uK the superficial velocity in K zone, and uK = V˙ K /S, S is the column cross area. The number of cells NK corresponds to the term HETPK DL in the van Deemter formula: K HETPK = HETPK DL + HETPkm = a +

b + cuK uK

t 1 − εt εt



K,n+1 K,n − qi,j qi,j

t

K,n+1 K,n − qi,j qi,j

i = 1, . . . , N C , mod; j = 2, . . . , N K , K = I, . . . , IV

K

K,n+1 K,n − Ci,j Ci,j

(6)

K represents the effect of axwhere HETPK DL = a + b/u ial dispersion as a sum of two terms of eddy diffusivity and molecular diffusion, respectively, while the value of K HETPK km = cu indicates the effect of mass transport resistances correlated with the lumped mass transport coefficient km . Since the contributions of axial dispersion and mass transport resistances to the column efficiency are additive, the NK value can be assumed to be constant and the effect of the velocity on the overall system efficiency can be simulated by the proper adjusting of the values of km coefficient. Note that the coefficient km strongly depends on the mobile phase composition and, less significantly, on the concentration of components to be separated. For systems with a very low, limited by kinetic effects efficiency both these dependencies should be accounted for; a detailed discussion of this issue can be found elsewhere (Antos, Kaczmarski, Pi˛atkowski, & Seidel-Morgenstern, 2003). For typical systems of moderate or high column efficiency all kinetic and dispersion effects can be lumped and accounted for indirectly by the number of cells NK . The km value can be assumed to be sufficiently high to mimic a permanent equilibrium between the solid and the liquid phases. In this case the kinetic equation (4) is kept for numerical reasons only (Antos & Seidel-Morgenstern, 2001). The same high k value is used for all components. Such a model coupled with the proper boundary conditions was solved in this

t

+

K,n K,n − Ci,j−1 Ci,j

τK

=0

(7)

K,n K,n K,n = km (˜qi,j (Cj ) − qi,j ),

i = 1, . . . , N C , mod; j = 2, . . . , N K , K = I, . . . , IV, n = 2, . . . , N end

(8)

where t is the time increment, which size has to be sufficiently small to reduce the influence of numerical diffusion on the accuracy of the solution (see details in Antos & SeidelMorgenstern, 2001). The scheme (7) and (8) has to be repeatedly solved for the components to be separated and optionally (for the gradient mode) for the modifier until a cyclic steady state is achieved. The criterion adopted to decide if this state is achieved was based on a check of the mass balance consistency between the masses of component i collected in the outlet ports and introduced into the feed during the switching period.

4. Adsorption isotherms The competitive adsorption between components 1 and 2 was described by the competitive Langmuir model: q˜ i (C1 , C2 ) =

Hi (Cmod )Ci , 2 1 + l=1 Keql (Cmod )Cl

i = 1, 2

(9)

where the values of the isotherm coefficients: the Henry’s constant Hi and the equilibrium constant Keq depend on the modifier concentration. To quantify the effect of the modifier concentration on the isotherm coefficients in Eq. (9) a power dependency was used: Hi (Cmod ) = (pHi Cmod )−rHi ,

i = 1, 2

−rKeqi

Keqi (Cmod ) = (pKeqi Cmod )

,

i = 1, 2

(10) (11)

For the modifier the single component Langmuir model was assumed: q˜ mod (Cmod ) =

Hmod Cmod 1 + Keq,mod Cmod

(12)

Because the concentration of the modifier typically exceeds markedly the concentration of the component to be separated, the competition between the sample and the modifier can be neglected. The parameters Hmod and Keq,mod as well as all pH , rH , pKeq , rKeq have to be determined experimentally.

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

The details of the adsorption model can be found elsewhere (Antos & Seidel-Morgenstern, 2001, 2002b; Antos, Pi˛atkowski, & Kaczmarski, 2000).

5. Adaptative random search optimization 5.1. Problem formulation In this work the following optimization problems have been solved: 1.

Min(EC) = f (mI , mII , mIII , mIV )

(13)

subject to the purity constraints: PuEx ≥ PuEx,min

(14)

PuRa ≥ PuRa,min

(15)

2.

In addition the following overall mass balances have to be met: mD = mI − mIV > 0

(21)

mEx = mI − mII > 0

(22)

mRa = mIII − mIV > 0

(23)

mF = mIII − mII > 0

(24)

where mEx , mRa and mF , mD correspond to the dimensionless flowrates of out- (mEx , mRa ) or in-coming (mF , mD ) streams (see Fig. 1). In order to simplify analysis a typical separation system has been considered, for which the pressure drop in the SMB unit was a limiting factor with respect to the process performance. The maximal dimensional flowrate max(V˙ K ) was set equal to unity: max(V˙ K ) = 1

or Max(Pu) = f (mI , mII , mIII , mIV )

(16)

In the above equations EC is the eluent consumption defined as EC =

V˙ D + V˙ F C2Ex V˙ Ex + C1Ra V˙ Ra

(17)

Pu is the purity of the out-coming streams: Pu = PuEx + PuRa

(18a)

where PuEx =

C2Ex C2Ex + C1Ex

(18b)

PuRa =

C1Ra C1Ra + C2Ra

(18c)

1581

(25)

The pressure constraint (25) was assumed as a scaling factor, which can be adjusted to pressure constraints of a real system. All the dimensional streams can be recalculated maintaining the optimal dimensionless decision variables. If optimal dimensionless mK are known the constraint max(V˙ K ) determines the value of switching time and the simulated solid phase flowrate (Eqs. (1) and (2)). Hence, the four dimensionless flowrates in combination with the pressure drop constraint determine all the streams in the SMB unit. If the column efficiency is a limiting factor (see Zhang et al., 2003) more detailed approach including plate equation (Eq. (6)) can be easily adopted. 5.2. Optimization algorithm

(20)

5.2.1. Multi-pass optimization algorithm In this work to optimize SMB chromatographic process the so-called multi-pass algorithm (MPLJ) proposed by Luus and Jaakola (1973), Luus (2001) was used. The algorithm is efficient and straightforward; the numerical code contains no more than 200 hundreds lines. Nevertheless, in general, random search algorithms require a number of repetitive solutions of a dynamic model for randomly generated decision variables. Hence, the attempts have been made to improve its effectiveness. The modification of the multi-pass algorithm proposed here allowed decreasing of the computation time compared to the ordinary Luus–Jaakola algorithm tested in our previous work (Ziomek et al., 2005).

Solving both the problems Min(EC) or Max(PR) leads to the similar set of the optimal decision variables, however, the employing of EC as an objective function was found to be more effective in the optimization algorithm.

5.2.2. Modification of multi-pass Luus–Jaakola algorithm The procedure used in the work is briefly summarized in the following:

with N end Ci,port =

n n=2 Ci,port t , t∗

i = 1, 2, port = Ex, Ra (19)

n Ci,port is calculated from the cell model (Eqs. (7) and (8)) in the position corresponding to the extract and raffinate port (the end of the section K = I and the section K = III, respectively). Additionally the value of the process productivity has been observed:

PR =

C2Ex V˙ Ex + C1Ra V˙ Ra V (1 − εt )N B

1582

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

(1) Set the number of iterations (N) in one pass, the number of generated points in each iteration (R), a total number of internal passes (NIQ), total number of external passes (NEQ) and the contraction factor for external passes (γ K ). I II III IV Set the starting point mK 0 = [m , m , m , m ], the K size of the initial search region (δ0 ) and the size of the final search region for each pass (δK f ) for each K variable. K , where mK is currently the best solution. Set mK = m opt opt 0 (2) Set the internal passes counter niq equal to 1 (niq = 1), and the external passes counter neq equal to 1 (neq = 1). Set the iterations counter n equal to 1 (n = 1). (3) Find random R values for each variable and check the decision constraints (Eqs. (21)–(24)). If the decision constraints are met, solve the model (7) and (8), calculate value of the objective function (e.g., Eq. (13)) and check the purity constraints (14) and (15). Accept the points, for which purity constraints are met. (4) Compare the values of the goal function, if better value K has been found, set mK opt = mn . (5) Reduce the size of the search region: K K δK n = β δn−1

(26)

The size of the search region βK is calculated according to the formula:  1/N δK f (27) βK = δK 0 (6) Increase the iteration counter n by 1, go to step 3. Repeat steps 3–5 until the number of current iterations is equal to N. (7) Increase the NIQ counter niq by 1. Calculate new size of the search region for a next internal pass from the linear equation: K K δK 0,niq = a niq N + b

(28)

where aK =

K δK f − δ0 NIQ × N

bK = δK 0

(29) (30)

βK was calculated from formula (27), which was suggested by Je˙zowski and Bochenek (2002) as a more effective in the search procedure. • In the original MPLJ algorithm the size of a new search region for a subsequent pass is calculated from the changes of decision variables in a current pass: K K δK 0,niq +1 = |mopt,niq −1 − mopt,niq |

(32)

In the optimization procedure for SMB process the search region calculated in such a way is likely to shrink drastically to very low values. Such a rapid narrowing of the search space brings a risk that the region containing the global optimum is cut off. Linear increase of the size of subsequent search regions in each new internal pass (see Eqs. (28)–(30)) allowed reducing the risk of trapping the calculations in a local optimum. In the presented modification the search bounds for each variable were varying dynamically after each of iterations. The best solution from a current iteration was set as center of the search region in a subsequent iteration. Typical illustration of changes in the size of search region in subsequent iterations is presented in Fig. 2. The modification proposed allowed increasing probability of locating a global minimum. The effectiveness of the procedure was preliminarily verified by comparing the results of calculations to those obtained with the original MPLJ algorithm for a typical test problem. Namely, a multi-modal function was selected and the optimization problem was solved for few starting points differing in the distance to the known global optimum. As a criterion for the comparison of the algorithms the success rate s.r. was adopted: s.r. =

NSGP NGP

(33)

where NGP is the number of all generated points (evaluations) for all the selected starting points and NSGP the number of successful runs (locating a global optimum) from among NGP. Typical results are presented in Section 7.1. Then, the effectiveness of the algorithm was verified for the SMB process. The optimization problem (Eqs. (13)–(15)) was solved with both the original and the modified algorithm

Go to step 3, and continue for the given number of NIQ. (8) Increase the NEQ counter neq by 1. Calculate the new size of the search region for a subsequent external pass: K K δK 0,neq = δ0 γ

(31)

The contraction coefficient (γ K ) could be dropped if NEQ is set to be equal to 1. (9) Go to step 2. Repeat points 2–8 until neq = NEQ. Two modifications has been introduced into the original MPLJ algorithm: • In the original MPLJ algorithm βK is a constant selected within the range 0.6–0.99, while in the modified algorithm

Fig. 2. Size of the search region in each of iterations for two external and four internal passes (NIQ = 4, NEQ = 2, N = 10).

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

1583

with the same number evaluations and the values of the objective function achieved with both the algorithms were compared (see Section 7.1). 5.2.3. Multi-objective optimization The optimization method presented can be employed for a single-objective function, for which the most effective direction for reducing the search region can be found. However, the random nature of the procedure allows easy switching between different objective functions. After each single-objective optimization the randomly generated points are analyzed; as the starting point for the second optimization the most favorable conditions for the subsequent objective function can be selected. Then, the optimization is switched into another direction driving solution rapidly to the next optimum.

6. Sensitivity analysis In addition to the optimization, the average sensitivities of the purity in the extract and in the raffinate port were evaluated under steady state conditions. According to Eqs. (18b) and (18c) the purity is a function of the concentrations of the components to be separated in the raffinate and the extract streams. The concentration profile is calculated by the use of the discrete dynamic model (e.g., Eqs. (7) and (8)) and depends on the model parameters, i.e., the flowrates, bed porosity, equilibrium parameters. If the equilibrium parameters and bed porosity can be assumed to be constant during the process (isocratic mode), the flowrates are the crucial parameters, which assure maintainability of the purities on the desired level. The relative sensitivity function of the purity with respect to flowrates can be described as V˙

SPuport =

∂Puport V˙ port , ∂V˙ port

port = Ex, Ra

(34)

Since all the flowrates in the SMB unit are mutually depenV˙

dant (see Eqs. (21)–(24)) the value of SPuport represent sensitivity of the purity to changes in one of the flowrates in the unit. The relative sensitivity of the purity with respect to the flowrate in the extract and the raffinate port (V˙ E and V˙ R ), can be determined by differentiating Eqs. (18b) and (18c): ∂PuEx ˙ VEx = ∂V˙ Ex (∂C1Ex /∂V˙ Ex )C2Ex − (∂C2Ex /∂V˙ Ex )C1Ex

V˙ R SPu,Ex

=−

(C1Ex + C2Ex )2

∂PuRa ˙ VRa ∂V˙ Ra (∂C1Ra /∂V˙ Ra )C2Ra − (∂C2Ra /∂V˙ Ra )C1Ra

V˙ Ex

(35)

VR SPu,Ra = ˙

=

(C1Ra + C2Ra )2

V˙ Ra

(36)

Fig. 3. Profiles  of the test functions (Michniev et al., 2000). Function de scription: F = ni=1 5j=1 j cos[(j + 1)xi + j], with bounds (0, 10), with two optimized parameters.

The ∂Ci,port /∂V˙ port are calculated by differentiating (analytically or numerically) in each time step n during the switching n+1 interval t* the finite difference solution Ci,j=port (Eqs. (7) and (8)) with respect to V˙ in position j corresponding to the extract and the raffinate port. The flowrate V˙ in Eqs. (7) and (8) is accounted for through the residence time (see Eq. (5)). The average sensitivity S¯ Pu,port in the switching period can be defined as follows: N steps V˙ S t ¯SPu,port = n=2 Pu,port , port = Ex, Ra (37) ∗ t Additionally, for the gradient mode, the sensitivity of the purity with respect to the modifier concentration Cmod was calculated ∂Puport mod SPu,port = (38) Cmod,port ∂Cmod,port 7. Procedures 7.1. Numerical test of the modified MPLJ algorithm For a preliminary test a multi-modal and hard to optimize using stochastic optimization methods function F has been selected (see Fig. 3). The test function was taken from Michniev, Stoyanov, and Keil (2000). Tests were performed with final size δf = 0.001. The number of the external passes (NEQ) was set to be equal to 2; the number of internal passes (NIQ) was a control parameter. For three different starting points the procedure was ran 2000 times. The success rate was calculated from Eq. (33). The results, which differed

1584

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

Table 1 Efficiency of the modified MPLJ algorithm Optimization algorithm

Original MPLJ Modified MPLJ a b

SMBa (goal function, EC)

Success rate (%) (test function, F) 4b

6b

8b

4

6

8

50.62 64.65

50.05 72.02

51.17 74.55

1.6981 1.3706

1.6637 1.2703

1.5353 1.1344

Calculated for the conditions CF = 3 vol.%, Cmod = 15 vol.% (see Fig. 4). Total pass numbers (NIQ × NEQ).

from the global optimum less than 0.01%, were accepted as successful runs. The global optimum was the real value of the goal function (Michniev et al., 2000). The results of calculations were compared to those obtained by the use of the original MPLJ algorithm. The comparison of the efficiency presented in Table 1 indicated visible improvement in the calculations efficiency, i.e., the probability of locating the global optimum was increased ca. 20–40%. The modified multi-pass Luus–Jaakola algorithm was found to be a robust tool for optimization of a multi-modal mathematical function. For several optimization tests performed for “easier” mathematical functions effectiveness of both the algorithms was found to be comparable. In the next step the SMB optimization was performed by the use of both the algorithms. The non-linear constrained problem (Eqs. (13)–(15)) was solved. For both the algorithms the same number of evaluations was used. The modified algorithm was found to be much more effective, offering better values of the objective function (see Table 1). 7.2. Searching procedure 7.2.1. Isocratic mode The optimization of isocratic mode was performed for three different feed concentration CF and constant modifier concentration. In order to improve effectiveness of the optimization for the starting point a set of dimensionless flowrates, for which the purity constraints were met, i.e., a “feasible point”, was selected. The search was organized in the following way: (1) Searching for the starting point: • To calculate value of the dimensionless flowrates m the triangle theory was used. • For the adequate values of mI , mIII , mIV the mII parameter was screened discretely by the use of the following equation: mII = mIII − y where y is the constant set to be initially equal to 0.1mII and decreased by the step of 0.01mII until the purity constraints were met. The last point, for which the purity constraints were respected was selected as a starting point. In order to find a feasible point several screening trials were sufficient.

(2) Optimization procedure: The optimization procedure was ran with the following control parameters: NEQ = 2, NIQ = 4, N = 10 and R = 10. Final sizes for each variable were set as δK f = 0.01, the initial size of the first pass was equal to δK = 0.4H 2 and contraction 0 factor γ K = 0.7. The total number of evaluations in the optimization procedure can be calculated as follows: NGP = N × R × NEQ × NIQ The number of evaluations depends on the difficulty of the separation, e.g., to optimize separation under strongly non-liner conditions in the isocratic mode (Cmod = 15 vol.%, CF = 3 vol.%) 1200 evaluations were required, with NEQ set to be equal to 3. The remaining optimized problems were solved with 800 generated points with (NEQ = 2). For each point the sensitivity of the purity in the extract and raffinate stream with respect to the flowrates was calculated. 7.2.2. Gradient mode To optimize gradient SMB process the following search procedure was applied: ¯ mod in the (1) The average concentration of the modifier C unit was set. This corresponds to the average modifier concentrations in the SMB unit under gradient conditions, whereby the desorption zones I and II operate at higher modifier concentrations, while the adsorption zones III and IV operate at lower modifier contents. (2) The best choice for this value is the optimal mobile phase composition for the isocratic process (Ziomek et al., 2005). ¯ mod the values of mI , mIV were calcu(3) For the assumed C lated from the triangle theory. (4) For the obtained mI , mIV parameters the optimization of purity in extract and raffinate streams was performed with mII and mIII as the decision variables (ca. 100 evaluations were required). (5) The point, for which the purity constraints were met was used as a starting point in the main optimization routine involving all the four decision variables. The optimization was performed for 800 generated points. For each point the sensitivity of the purity in the extract and raffinate streams with respect to the flowrates and the modifier concentration were calculated.

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

8. Result and discussion 8.1. Isocratic mode The calculations were performed for the selected modifier concentration Cmod = 15 vol.%, and the feed concentrations: CF = 0.3, 1, 3 vol.%. In order to increase the speed of calculations the number of cells in each zone was assumed as low as NK = 50. For sake of simplicity the influence of the flowrate variations on the system efficiency was neglected. Reduced purity demands were assumed, i.e., PuEx,Ra ≥ 90%. The simplest configuration of the unit was assumed: 1/1/1/1, i.e., one column per each K zone. The following column characteristics were assumed: length and column diameters L = 15 cm and d = 0.46 cm, respectively, total porosity εt = 0.75. The isotherm parameters Hmod and Keq,mod as well as all pH , rH , pKeq , rKeq were assumed the same as obtained for the separation of the two cycloketones mixture in a normal phase system analyzed by Antos and Seidel-Morgenstern (2001, 2002a). Typical values of isotherm parameters for different mobile phase composition are summarized in Table 2. The points generated in the optimization procedure, which met purity constraints are depicted in Fig. 4; the regions predicted by the triangle theory are overlaid on the same plot. It is evident that the location and shape of regions obtained by Table 2 Typical isotherm coefficients for the component to be separated Cmod

H1

H2

Keq1

Keq2

15 20 30

4.777 3.479 2.225

7.195 5.420 3.635

0.174 0.126 0.079

0.369 0.281 0.191

For the modifier Hmod = 0.92, Keq,mod = 0.05 vol.−1 .

Fig. 4. Regions of the (mIII –mII ) plane with points generated in the optimization procedure, for which the purity constraints were met (PuEx,Ra ≥ 90%). Isocratic mode, Cmod = 15 vol.%, CF = 0.3, 1, 3 vol.%. Solid lines: the predictions by the use of the triangle theory; symbols: the points generated by the random search procedure.

1585

the use of the optimization procedure and the triangle theory are similar and follow the same trend — with an increase of overloading (non-linearity of the isotherm) the operating window shrinks and orientates closer to the diagonal. However, as mentioned above, the triangle method provides only a rough estimation of the optimal conditions, e.g., it underestimates the length of the separation region. Feasible points situated close to the vertex of the triangle correspond to the maximal throughput and to the optimal conditions. V˙ V˙ The sensitivity values S¯ Pu,Ex and S¯ Pu,Ra versus purity of the outlet streams are depicted in Fig. 5a and b. The analysis corresponds to the region of CF = 1 vol.% shown in Fig. 4. As it can be expected, the vertex position of the optimal point in the operating region is related to high absolute values of the sensitivity in the extract port as well in the raffinate port. However, the position in the operating region concerns only two operating parameters: mIII , mII and cannot provide general information concerning maintainability of the purity constraints. If the flowrates of the regenerating streams: mI , mIV do not assure complete regeneration of the solid and fluid phases, which is allowed for reduced purity constraints, all the decision variables influence on the sensitivity values. Therefore, the sensitivity of the purity can be also high for points situated inside the separation region. The plots presented in Fig. 5a and b indicate that generV˙ V˙ ally the sensitivity S¯ Pu,Ex and S¯ Pu,Ra decrease with increasing purity of the extract and the raffinate streams. The points corresponding to the highest purity are also related to the highest maintainability of the purity at the desired level. Nevertheless, it is apparent, particularly for the raffinate stream, that points corresponding to the highest purities are characterized by different sensitivity of the purity constraints. For the process realization the points corresponding to relatively low values of sensitivity could be selected. Note that the sensitivity of the raffinate purity is higher than that for the extract. This results from the fact that unit operate under non-linear conditions, for which self-sharpened fronts are formed. The raffinate stream is contaminated by the self-sharpened front of the extract, while the extract stream is contaminated by the dispersed front of the raffinate. The positions of self-sharpened fronts are obviously more difficult to control. Typical concentration profiles corresponding to the optimum are presented in Fig. 6. V˙ V˙ In Fig. 7a and b the plots of S¯ Pu,Ex and S¯ Pu,Ra versus the objective function EC are shown. Note that operating under optimal conditions related to the highest process efficiency does not lead to the highest sensitivity values. For a number of points situated inside the region the sensitivity is higher. This means that the locating of operating points between the bounds of the region does not guarantee improving robustness of the process. Such a geometrical “safety margin” is illusory for a real process realized in systems with a limited efficiency. Hence, the sensitivity analysis should be considered as an additional performance index providing information, which concern maintainability of the purity at the desired level.

1586

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

Fig. 5. Sensitivity of the purity in (a) the extract stream and (b) the raffinate stream with respect to the flowrates vs. the purity for the feasible points (PuEx,Ra ≥ 90%) generated during the optimization. Isocratic mode, Cmod = 15 vol.%, CF = 1 vol.%.

Fig. 6. Internal concentration profiles for the optimum under steady state conditions at the end of the switching time. Solid line: CF = 0.3 vol.%; dashed line: CF = 1 vol.% (see Fig. 4).

8.2. Gradient mode The optimization for the gradient process was performed for CF = 0.3 and 1 vol.%. The average modifier concentration in the units was assumed according to the optimal modifier

concentration in the isocratic mode, e.g., for CF = 0.3 vol.%, ¯ mod = 20 vol.%; for CF = 1 vol.%, C ¯ mod = 30 vol.% (see C details in Ziomek et al., 2005). The concentration of the modifier in the feed and the desorbent port was set to Cmod,F = 0 vol.% and Cmod,D = 100 vol.%, respectively. This corresponded to the maximal “driving force” for the gradient. However, the improvement of the efficiency reported in the following does not represent a general tendency; the efficiency of the gradient mode depends mainly on the dependence of the Henry’s constants on the modifier content. The particular system being studied concerned compounds exhibiting relatively weak retention dependency. Moreover, for high overloadings the whole “driving force” for the gradient (i.e., Cmod,F = 0 vol.% and Cmod,D = 100 vol.%) cannot be always exploited (Ziomek et al., 2005). Two regions of the operating parameters calculated for CF = 0.3 and 1 vol.% are depicted in Fig. 8. The optimum of the objective function is significantly improved compared to that obtained for the same feed concentration in the isocratic mode (compare the data in Table 3). V˙ V˙ The sensitivity S¯ Pu,Ex and S¯ Pu,Ra versus the objective function EC are plotted in Fig. 9. The optimization results for the case of the separation illustrated (CF = 1 vol.%) is particu-

Fig. 7. Sensitivity of the purity in (a) the extract stream and (b) the raffinate stream with respect to the flowrates vs. the objective function EC for the feasible points generated during the optimization. Isocratic mode, Cmod = 15 vol.%, CF = 1 vol.%.

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

1587

Table 3 Comparison of the optimization results obtained for isocratic and gradient modes CF (vol.%) Gradient Isocratic, Cmod = 15 vol.% (Cmod = opt)a Gradient Isocratic, Cmod = 15 vol.% (Cmod = opt)a a

0.3

1

PuEx

PuRa

PR

EC

u SPu,Ex

u SPu,Ra

mod SPu,Ex

mod SPu,Ra

0.967 0.900 (0.900) 0.900 0.900 (0.901)

0.957 0.900 (0.900) 0.970 0.901 (0.900)

0.113 0.042 (0.044) 0.271 0.103 (0.118)

1.999 4.347 (4.398) 0.723 1.615 (1.470)

0.0383 0.079 (0.084) 0.086 0.081 (0.060)

0.711 −0.112 (−0.113) −0.025 −0.132 (−0.154)

0.740 –

−5.342 –

1.186 –

−0.574 –

For the optimized mobile phase compositions, i.e. for CF = 0.3 vol.%, Cmod = 20 vol.%; for CF = 1 vol.% Cmod = 45 vol.%.

Fig. 8. Regions of the (mIII –mII ) plane with points, which met the purity constraints (PuEx,Ra ≥ 90%). Gradient mode, CF = 0.3, 1 vol.%.

larly favorable for the raffinate purity. High raffinate purities can be achieved at comparable or lower sensitivity values compared to the isocratic mode and at the significantly improved value of the objective function (compare in Table 3). Moreover, the purity for the gradient mode can be improved

Fig. 9. Sensitivity of the purity in the extract (left scale) and the raffinate (right scale) streams with respect to the flowrate vs. the objective function EC for feasible points generated during the optimization. Gradient mode, CF = 1 vol.%.

almost without worsening the values of the objective function while both these performance indexes are conflicting for the isocratic mode (see Fig. 10). The plot of the conflicting performance indexes, which is typical for multiobjective optimization, is generated directly during the random search calculations. However, for a single-objective optimization the complete plot containing the optima of both objectives cannot be created. For this, the optimization have to be performed in the two steps: (1) the optimum for the EC function is found; (2) among the randomly generated points in the EC optimization the set of decision variables, which is the most favorable for the objective Pu is selected and used as a starting point in the subsequent optimization. The number of evaluation necessary to identify the global optimum in the second optimization was equal only ca. 150 (e.g., NEQ = 1, NIQ = 3, N = 7 and R = 7). The trend illustrated in Figs. 9 and 10 is not general. The improvement in the process efficiency for the gradient mode depends on the isotherm non-linearity and it is particularly enhanced for relatively low feed concentration (see Table 3) while diminishes slowly with an increase of overloading. The gradient process can be expected to be sensitive to the concentration profile of the modifier. Therefore, the sensitivity of the purity constraints with respect to the modifier concentration was calculated. The results can be observed in

Fig. 10. Comparison of the optimization results for isocratic and gradient mode for CF = 1 vol.%. Solid squares: gradient mode; open squares: isocratic mode. Left upper corner: full scale of EC for the isocratic mode (created in the multiobjective optimization).

1588

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589

References

Fig. 11. Sensitivity of the purity in the extract and the raffinate stream with respect to the modifier concentration vs. the EC function for feasible points generated during the optimization. Gradient mode, CF = 1 vol.%. mod Fig. 11. The values of sensitivity S¯ Pu,Port are approximately V˙ one order of magnitude higher compared to that of S¯ Pu,port . This indicates that in the gradient operation maintaining of the modifier concentration at the level predicted is crucial for successful separation.

9. Conclusions In this work an efficient numerical tools has been proposed for optimization of SMB process. The main efforts have been focused on the designing of the solvent gradient SMB operation, which was proved to outperform isocratic mode with respect to the eluent consumption and the process productivity. The multi-pass Luus–Jaakola algorithm has been modified and used for the optimization of the separation process. For predictions of the objective functions, i.e., the minimum of eluent consumption the equilibrium stage model accounting for real system efficiency was employed and implemented in the optimization procedure. Such an optimization tool allowed the construction of the operating window of parameters guaranteeing successful separation. For feasible operating points generated during the random search calculations several performance indexes were simultaneously observed: eluent consumption (the goal function), extract and raffinate port purities (constraints), productivity, and sensitivity of the purity with respect to such operating parameters as the flowrates and the modifier concentration. The results of calculations indicated that a number of operating parameters correspond to the similar process efficiency but different maintainability of the purity at the desired level. The sensitivity analysis was found to be helpful for balancing high efficiency of the process with its high robustness. The results of the sensitivity analysis revealed that in the gradient operation, which performance is superior compared to the isocratic mode, maintaining of the modifier concentration at the level predicted is a factor of major importance.

Abel, S., Mazzotti, M., & Morbidelli, M. (2002). Solvent gradient operation of simulated moving beds. Journal of Chromatography A, 944, 23. Abel, S., Mazzotti, M., & Morbidelli, M. (2004). Solvent gradient operation of simulated moving beds. 2. Langmuir isotherms. Journal of Chromatography A, 1026, 47. Antos, D., Pi˛atkowski, W., & Kaczmarski, K. (2000). Determination of mobile phase effect on single-component adsorption isotherm by use of numerical estimation. Journal of Chromatography A, 874, 1. Antos, D., & Seidel-Morgenstern, A. (2001). Application of gradients in the simulated moving bed process. Chemical Engineering Science, 56, 6667. Antos, D., & Seidel-Morgenstern, A. (2002a). Two-step gradient in simulated moving bed chromatography. Journal of Chromatography A, 944, 77. Antos, D., & Seidel-Morgenstern, A. (2002b). Continuous step gradient elution for preparative separations. Separation Science and Technology, 37, 1469. Antos, D., Kaczmarski, K., Pi˛atkowski, W., & Seidel-Morgenstern, A. (2003). Journal of Chromatography A, 1006, 51. Biressi, G., Ludemann-Hombourger, O., Mazzotti, M., Nicoud, R., & Morbidelli, M. (2000). Design and optimisation of a simulated moving bed unit: Role of deviations from equilibrium theory. Journal of Chromatography A, 876, 3. Houwig, J., van Hateren, S. H., Billiet, H., & van der Wielen, L. A. M. (2002). Effect of salt gradients on the separation of dilute mixtures of proteins by ion-exchange in simulated moving beds. Journal of Chromatography A, 952, 85. Jensen, T. B., Reijns, T. G. P., Billiet, H. A. H., & van der Wielen, L. A. M. (2000). Novel simulated moving-bad method for reduced solvent consumption. Journal of Chromatography A, 873, 149. Je˙zowski, J., & Bochenek, R. (2002). Experiences with the use of Luus–Jaakola algorithm and its modifications in optimization of process engineering problems. In R. Luus (Ed.), Recent developments in optimization and optimal control in chemical engineering. Trivandrum, India: Research Signpost. Luus, R. (2001). Use of Luus–Jaakola optimization procedure for singular optimal control problems. Nonlinear Analysis, 47, 5647. Luus, R., & Jaakola, T. H. I. (1973). Optimization by direct search and systematic reduction in the size of search region. AIChE Journal, 19, 760. Mazzotti, M., Storti, G., & Morbidelli, M. (1994). Robust design of countercurrent adsorption separation processes. 2. Multicomponent systems. AIChE Journal, 40, 1825. Mazzotti, M., Storti, G., & Morbidelli, M. (1996). Robust design of countercurrent adsorption separation. 3. Nonstoichiometric systems. AIChE Journal, 42(10), 2784. Mazzotti, M., Storti, G., & Morbidelli, M. (1997). Optimal operation of simulated moving bed units for non-linear chromatographic separations. Journal of Chromatography A, 769, 3. Michniev, I., Stoyanov, S., & Keil, F. J. (2000). Investigation and modification of the Luus–Jaakola global optimization algorithm. Hungarian Journal of Industrial Chemistry, 28, 231. Migliorini, C., Mazzotti, M., & Morbidelli, M. (1998). Continuous chromatographic separation through simulated moving beds under nonlinear conditions. Journal of Chromatography A, 827, 161. Migliorini, C., Mazzotti, M., & Morbidelli, M. (1999). Design of simulated moving-bed units under nonideal conditions. Industrial and Engineering Chemistry Research, 38, 2400. Ray, A. K., Carr, R., & Aris, R. (1994). The simulated countercurrent moving bed chromatographic reactor: A novel reactor-separator. Chemical Engineering Science, 49, 469. Ruthven, D. M., & Ching, C. B. (1989). Countercurrent and simulated countercurrent adsorption separation process. Chemical Engineering Science, 44, 1011–1038.

G. Ziomek, D. Antos / Computers and Chemical Engineering 29 (2005) 1577–1589 Storti, G., Masi, M., Carra, S., Mazzotti, M., & Morbidelli, M. (1989). Optimal design of multicomponent countercurrent adsorption separation process involving nonlinear equilibria. Chemical Engineering Science, 44, 1329. Storti, G., Mazzotti, M., Morbidelli, M., & Carra, S. (1993). Robust design of binary countercurrent adsorption separation processes. AIChE Journal, 39, 471. Subramani, H. J., Hidajat, K., & Ray, A. K. (2003). Optimization of reactive and Varicol systems. Computers and Chemical Engineering, 27, 1883. Zhang, Z., Hidajat, K., & Ray, A. K. (2002). Multiobjective optimization of simulated countercurrent moving bed chromatographic reactor

1589

(SCMCR) for MTBE synthesis. Industrial and Engineering Chemistry Research, 41, 3213. Zhang, Z., Hidajat, K., Ray, A. K., & Morbidelli, M. (2002). Multiobjective optimization of simulated moving bed and Varicol process for chiral separation. AIChE Journal, 48(12), 2800. Zhang, Z., Mazzotti, M., & Morbidelli, M. (2003). Multiobjective optimization of simulated moving bed and Varicol processes using genetic algorithm. Journal of Chromatography A, 989, 95. Ziomek, G., Kaspereit, K., Je˙zowski, J., Seidel-Morgenstern, & Antos, D. (2005). Effect of mobile phase composition on the SMB process efficiency. Stochastic optimization of isocratic and gradient process. Journal of Chromatography A, manuscript submitted for publication.