Bilevel Optimizing Control Structure for a Simulated Moving Bed

Ind. Eng. Chem. Res. 2010, 49, 3689–3699

3689

Bilevel Optimizing Control Structure for a Simulated Moving Bed Process Based on a Reduced-Order Model Using the Cubic Spline Collocation Method Kiwoong Kim,† Kwang Soon Lee,*,† and Jay H. Lee‡ Department of Chemical and Biomolecular Engineering, Sogang UniVersity, 1 Shinsoodong, Mapogu, Seoul 121-742, Korea, and School of Chemical and Biomolecular Engineering, Georgia Institute of Technology, 778 Atlantic DriVe, Atlanta, Georgia 30332-0100

A bilevel strategy for optimizing control of a simulated moving bed (SMB) process is proposed. In the lower level, repetitive model predictive control (RMPC) is used to regulate product purities; in the upper level, optimal feed/desorbent flow rates and the switching period are determined. Both levels employ a fundamental SMB model reduced to a set of nonlinear discrete-time dynamic equations using the cubic spline collocation method and exact discretization. For RMPC, the SMB model is linearized successively along the operating trajectories seen in the previous switching period. It is assumed that the flow rates can be varied within a switching period and the average product purities over each switching period can be measured albeit with a significant analysis delay. Numerical studies using linear isotherms showed that the proposed strategy is successful at driving the process to the intended optimum and maintaining it there while robustly regulating the product purities despite various uncertainties. Introduction The simulated moving bed (SMB) is a real-world approximation of the countercurrent true moving bed (TMB), which is an idealized continuous chromatography process. In SMB, continuous feeding and product withdrawal are enabled by periodic switching of the inlet and outlet ports in the direction of the desorbent flow.1 SMB was first commercialized in the petrochemical industry for the separation of xylene isomers and other hydrocarbon mixtures2 and, later, has been used in the food industry for the separation of fructose from glucose. Today it is widely employed in the pharmaceutical and biochemical industries for various difficult separations, e.g., those involving mixtures of enantiomers or of heat-sensitive products. As the importance of SMB grew, nonconventional SMB operation techniques, e.g., VariCol,3 solvent gradient,4 ModiCon,5 and PowerFeed,6,28 have been developed as a way to enhance the productivity and to reduce the desorbent consumption. A typical SMB process is run under a cyclic steady state (CSS) where process variables vary with time along their respective periodic trajectories due to a periodic switching of the columns. Hence, an important objective in SMB operation can be stated as “seeking and settling on a CSS where desired product purities are obtained with maximum productivity and minimum desorbent consumption”. For the initial determination of operating condition, the celebrated Mazzotti’s triangle theory7 can be used for the case of constant flow rates. Alternatively, numerical search can be performed based on a mathematical model of the SMB process. So, the determined operating condition is revised further after the implementation in order to compensate for the effects of various uncertainties. Control of SMB poses many challenges. First and foremost, due to its cyclic dynamics, conventional feedback control methods intended for a process with a fixed steady state may perform poorly. The availability of product purity measurements varies from case to case, but it can definitely cause additional * To whom correspondence should be addressed. E-mail: kslee@ sogang.ac.kr. Phone: +82-2-705-8477. Fax: +82-2-3272-0319. † Sogang University. ‡ Georgia Institute of Technology.

complications in the control system design. In the case that the measurements are infrequent and delayed, the controller has to rely heavily on the SMB model and model accuracy becomes a major issue. The adsorption isotherm, dead volume, axial dispersion, and bed porosity are major sources of uncertainties to overcome in this case. Once the controller is designed that can ensure on-spec separation, maximizing the product throughput while minimizing the desorbent consumption becomes the main goal in the SMB operation. The process control research community has attempted to address the above-mentioned difficulties. Kloppenburg and Gilles8 proposed a nonlinear control method using asymptotically exact input/output linearization and nonlinear state estimation on the basis of a hypothetical TMB model. Klatt et al.9 exploited a two-layer control architecture where an offline dynamic optimizer based on a rigorous SMB model calculates the optimal flow rates and switching period in the upper level and multiple SISO controllers regulate the process in the lower level. The optimizer minimizes the desorbent usage at a CSS, imposing inequality constraints on the extract and raffinate purities. The controllers take action once every switching period to maintain the so-called front positions at the reference values provided by the optimizer. Toumi and Engell10 presented a nonlinear model predictive control scheme for a reactive SMB process. The controller minimizes the production cost online by enforcing the purity specifications using inequality constraints. They considered the switching period as one of the manipulated variables (MV’s) and relied on a detailed process model and the successive quadratic programming (SQP) technique for the calculation of the MV’s. In the above work, the MVs are updated only once every switching period. As evidenced by the new SMB operations like Modicon5 and PowerFeed;6 however, the productivity and desorbent curtailment can be enhanced further when the MVs are allowed to change within a switching period. In such a case, the SMB operation becomes repetitive and the so-called repetitive control (RC) technique11,12 can greatly improve the operation. The RC technique, based on the internal model principle, can completely eliminate the tracking offset despite model uncertainty and repeated disturbances. Lee et al.13

10.1021/ie901121y  2010 American Chemical Society Published on Web 03/23/2010

3690

Ind. Eng. Chem. Res., Vol. 49, No. 8, 2010

proposed a time-domain RC technique called RMPC (repetitive model predictive control), which conducts period-to-period feedback together with real-time feedback within the framework of MPC. They applied it to the start-up problem of a numerical SMB process.14 The idea of RMPC was exploited and extended by Erdem et al.15 and Abel et al.16 in their construction of a receding horizon optimizing controller for an SMB process. They considered the same optimizing control objective as in the work of Toumi and Engell10 but derived their RMPC law by allowing multiple control actions within a switching period and using a time-varying linear model derived from the linearization of a discrete-time nonlinear model (obtained from applying the finite difference method (FDM)) along the nominal periodic trajectories. The optimizing control method was successfully implemented in an experimental SMB system with continuous purity measurement17 and also discrete purity measurement averaged over a switching period.18 In the present contribution, a novel bilevel optimizing control structure for the SMB process is proposed where an RMPC layer is cascaded by a CSS optimizer. The proposed scheme has some novel aspects compared to the previously published ones. First, the nominal SMB model was derived using the cubic spline collocation method (CSCM) and exact discretization instead of the finite difference approximation. The linear time-varying model for the RMPC design is obtained by linearizing the nominal model and is updated after each switching operation by relinearizing the model around the operating trajectories seen in the previous operation period. The optimizer performs the CSS optimization on the basis of the nominal SMB model with switching period and feed/desorbent flow rates as decision variables, whereas the RMPC conducts the regulation of product purities averaged over a switching period. Product purities are assumed to be measured once every period with some analysis delay. The same measurement situation has also been considered by other researchers to reflect the industrial practice.18,19 All the flow rates are allowed to vary within a switching period as in PowerFeed and ModiCon operations. In the proposed scheme, only the RMPC is run continuously and the optimizer is invoked only when a new optimum is thought to be needed, which is in contrast with the existing methods.10,16 The performance of the proposed technique is demonstrated through simulation of various perturbation scenarios for the case of linear adsorption isotherms.

Figure 1. Schematic diagram of a four-zone SMB.

hypothetical adsorbed phase concentration equilibrated with the bulk phase concentration, D, k, H, S, and ε denote the axial dispersion coefficient, mass transfer coefficient, equilibrium constant, column cross-sectional area, and bed porosity, respectively, and tˆ and zˆ represent time and axial distance, respectively. Let L and T represent the column length and the duration of a switching operation, respectively. Using the normalized variables, t } tˆ/T, z } zˆ/L, η } (1 - ε)/ε, uj } (QjT)/(εV), λi } DiT/L2, and κi } kiT, where V ) SL, eqs 1-3 can be rewritten as ∂cij ∂wij ∂cij ∂2cij + uj ) λi 2 , +η ∂t ∂t ∂z ∂z

t, z ∈ [0.1]

dwij ) κi(Hicij - wij) dt

(4)

(5)

The initial and boundary conditions (IC and BC’s) are given as IC: BC:

cij(0, z) ) wij(0, z) ) 0

cij(t, 0) ) cin ij (t),

∂cij ∂z

|

∞

(6)

|

∂cij )0 ∂z zL for a sufficiently large zL (7)

≈

Numerical SMB Models Governing Equations. We consider a four-zone SMB with identical columns as shown in Figure 1, where the feed stream contains binary components, A and B, and volumetric flow rates for each stream are represented by Qj. The column number, I to IV, is the logical index and rotates along the physical columns with the port switching. Assuming the linear adsorption isotherms, the mass balance for component i ∈ {A, B} in the jth column can be written as

( )

2

∂cij Qj ∂cij ∂ cij (1 - ε) ∂wij + + ) Di 2 , ε εS ∂zˆ ∂tˆ ∂tˆ ∂zˆ i ) A, B, j ) I, ..., IV

(1)

dwij ) ki(w*ij - wij) dtˆ

(2)

w*ij ) Hicij

(3)

where c and w represent the concentrations in the bulk phase and adsorbent in moles per cubic centimeter, w* represents the

The second BC in eq 7 was introduced by Guiochon et al.20 and called the far-side BC (FSBC) by Yoon et al.21 It is a physically more reasonable representation of the normally considered BC of ∂cij/∂z|z)1 ) 0. In addition to the individual column model, the following mass balances should be satisfied at the nodal points for inlet and outlet ports: QI ) QIV + QD, ciI(t, 0)QI ) ciIV(t, 1)QIV QII ) QI - QE, ciII(t, 0) ) ciI(t, 1) ) ciE QIII ) QII + QF, ciIII(t, 0)QIII ) ciII(t, 1)QII + ciFQF QIV ) QIII - QR, ciIV(t, 0) ) ciIII(t, 1) ) ciR, i ) A, B

(8) Recall that t in eqs 4 and 5 is normalized with respect to T. It can be shown that the column performance remains the same despite changes in T as long as uj values are kept constant when λi and κi are sufficiently small and large, respectively. Cubic Spline Collocation Method (CSCM).21-23 The CSCM was employed to reduce the partial differential equation


(PDE) model in eqs 4 and 5 to a set of ODE’s since the polynomial-based orthogonal collocation method (OCM) is prone to generate spurious ripples in the resulting model for high Peclet number columns.21 FDM was not considered since the order of an FDM-based ODE model is usually much higher than that of a collocation-based model, and moreover, the FSBC, when applied, increases the model order significantly. In the CSCM, the interpolation function is given as a concatenation of piecewise cubic polynomials defined over each subdomain such that the adjacent polynomials satisfy some smoothness conditions, called the spline conditions, at the connecting points. Hence, the CSCM can yield a much smoother profile than the OCM. Indeed, the CSCM can be considered as a special case of the OCM or the Galerkin’s method (GM) on finite elements (FE), which have been adopted for accurate numerical modeling of the SMB system.24 Normally, the OCMFE and GM-FE require that the function value and first-order derivative are continuous at the nodal points between the subdomains. If continuation of second-order derivative is additionally enforced and a cubic polynomial is chosen as the trial function for each subdomain, all the coefficients can be determined by the interpolation condition at the nodal points and both methods converge to the CSCM. Let [0, zL] be divided into nc + 1 intervals with the collocation points at z0 () 0), z1, ..., znc () 1), and znc+1 () zL), respectively. The cubic spline function pk(z) for the kth interval [zk-1, zk] represents a cubic polynomial that satisfies pk(zk) )

pk+1(zk), p(1) k (zk)

)

(1) pk+1 (zk), p(2) k (zk)

) k ) 1, ..., nc (2) pk+1 (zk),

[ ] (1)

˜ } A

T

m (z0) l m(1)(znc+1)T

[ ] (2)

and

B˜ }

(nc+2)×(nc+2)

3691

T

m (z0) l m(2)(znc+1)T

(nc+2)×(nc+2)

Reduced-Order Ordinary Differential Equation (ODE) Models. Let c˜ij(t, z) be the interpolation function for cij(t, z) and let cij } [c˜ij(t, z0) · · · c˜ij(t, znc) c˜ij(t, znc + 1)]T. Let us define ˜ ij similarly. The method of collocation dictates w˜ij(t, z) and w that eqs 4 and 5 are satisfied by the interpolation functions at the collocation points. This leads to, after dropping the subscript ij for notational simplicity, dc˜ dw ˜ ˜ )c˜ ) (λB˜ - uA +η dt dt

(14)

dw ˜ ) κ(Hc˜ - w ˜) dt

(15)

The PDE’s hold in the domain (open set) and the BC’s hold at the boundary points. Accordingly, the equations at z0 and zn+1 should be excluded from eqs 14 and 15, and the BC’s are applied to the remaining ODE’s. After rearrangement, the ODE column model is derived as follows: dc ) Epc + ηκw + dpc(t, 0) dt

(16)

dw ) κ(Hc - w) dt

(17)

where the superscript (n) denotes the nth order derivative. Each cubic polynomial has four coefficients, and we have 4(nc + 1) coefficients to determine. The spline condition in eq 9 specifies 3nc relations for the coefficients. The interpolation condition at the nc internal collocation points and the two BC’s provide nc + 2 conditions. Therefore, two more conditions are needed for the cubic polynomials to be uniquely determined. For the additional conditions, the followings are generally considered:

where c(t) } [c(t, z1) · · · c(t, znc)]T and w is defined similarly. Definitions of Ep and dp are given in Appendix A with a brief description of the derivation procedure. The numerical process used for the simulation is given as eqs 16 and 17 with 25 equally spaced internal collocation points (nc ) 25) at {0, 0.04, ..., 0.96, 1, 5}. The last component “5” is the point zL in eq 7 for the far-side boundary condition.21 For the nominal model used for the control calculation, 10 equally spaced internal collocation points (nc ) 10) were chosen and also adsorption equilibrium

(2) p(2) 1 (z0) ) pnc+1(znc+1) ) 0

wij ) Hicij

(9)

(10)

The interpolation function is a concatenation of the cubic polynomials for each zone. Using Sk(z) defined by Sk(z) ) 1 for [zk-1, zk] and 0 otherwise, the interpolation function is represented by

(18)

was assumed instead of the rate process to lower the computational burden. This simplifies the column equation as dc ) Emc + dmc(t, 0) dt

(19)

nc+1

y˜(z) }

∑ θ p (z)S (z) k k

k

(11)

k)1

Suppose that y(z) is interpolated at zk, k ) 0, 1, ..., nc + 1 by y˜(z). After some straightforward manipulations, y˜(z) can be represented as y˜(z) ) m(z)Ty ) m(z)Ty˜

(12)

where m(z) is an (nc + 2) × 1 vector function of z; y } [y(z0) · · · y(znc) y(znc+1)] and y˜ is defined similarly. Now, the first and second order derivatives of y˜ at the collocation points can be written as dy˜ ˜ y˜ )A dz where

(and)

d2y˜ ) B˜y˜ dz2

Equations 16 and 17 (likewise eq 19) represents the model equations for a single column. The entire SMB models for the numerical process and also the nominal model are constituted by stacking the respective equations of four individual columns and combining them with the nodal point mass balances in eq 8 and a port rotation relationship. We assume the following generic form to express the resulting models:

(13)

dz j (u, v)z + G j (u, v) )F dt

(20)

In the above, the state z for the numerical process is composed of the fluid and adsorbed phase concentrations of A and B at the collocation points in the four columns; z for the nominal model is similarly defined but does not contain adsorbed phase concentrations; u and v represent the MV’s for RMPC and the decision variables for the optimizer, respectively, which were chosen as

3692


u}

[

]

QIT QIIT T , εV εV

v}

[

QFT QDT εV εV

]

T

(21)

Indeed, the optimizer determines optimum u and v, but only v is implemented in the process, the reason for which will be explained later. It is assumed that the products are separately collected for each switching period, and the average purities over a period are measured after each period with some analysis delay. The nominal ODE model is transformed to a discretetime equivalent model and linearized as described subsequently. Discrete-Time Nominal Model. The ODE model in eq 20 is linear in z, and the exact discrete-time equivalent model can be readily derived if zero-order hold is assumed for u and v.25 If we let t and h represent the discrete-time index and sampling interval, respectively, the discrete-time model for the kth switching period can be obtained as zk(t + 1) ) F(uk(t), vk(t))zk(t) + G(uk(t), vk(t)), zk(0) ) zk-1(N)

(22) j (t)h) and G(t) } F j (t)-1[exp(F j (t)h) where N } T/h; F(t) } exp(F j (t); F j (t) is the short-hand notation of F j (u(t), v(t)), similarly - I]G j (t), and G(t). for F(t), G The average purity over a switching period is defined as p j (N) where p¯(t) }

(

[ ] pjAE(t) pjBR(t)

)

sAE(t) sAE(t) + sBE(t) sAE(t) ) sAE(t - 1) + QE(t - 1)cAE(t - 1) sBE(t) ) sBE(t - 1) + QE(t - 1)cBE(t - 1) and similarly for pjR(t), sAR(t) and sBR(t)

pjE(t) )

(23)

(24)

Figure 2. Overall structure of the optimizing control system.

determines QF, QD, QI, QII, and T simply based on the SMB model but implements only QF, QD, and T leaving QI and QII to be recalculated by RMPC. In this respect, the performance of the optimizer is completely dependent on the model quality whereas the performance of RMPC is not due to the cyclewise feedback. Linearized Model for Controller Construction. RMPC is a linear model-based technique whereas the SMB model in eq 25 is nonlinear. To reflect the nonlinear aspects of the process in the controller design as much as possible, the SMB model is repeatedly linearized around the operating trajectories seen during the previous switching period and used for the controller design. By representing the model error, process disturbance, and measurement error by random processes, the linearized model for the kth period can be readily derived as ∆xk(t + 1) ) Ak-1(t)∆xk(t) + Bk-1(t)∆uk(t) + w(t) yk(t) ) ∆p¯k(t) ) yk-1(t) + Hk-1∆xk(t) m yk+1(d) ) yk(N) + n

(26) and the subscripts E and R represent the extract and raffinate, respectively. Incorporating eqs 23 and 24 into eq 22 and augmenting the state as x } [zT sAE sBE sAR sBR]T yields

where ∆xk(t) } xk(t) - xk-1(t), ∆uk(t) } uk(t) - uk-1(t) ∂A ∂(A + B) ∂H , Bk-1 } , Hk-1 } T Ak-1(t) } T T k-1 ∂x k-1 ∂u ∂x

|

xk(t + 1) ) A(uk(t), vk(t), xk(t)) + B(uk(t), vk(t)), xk(0) ) [zk-1(N)T 0 0 0 0 ]T yk(t) ) p¯k(t) ) H (xk(t)), zk+1(0) ) zk(N) m (d) ) yk(N) + n yk+1

(25) where ym and n represent the purity measurement and zeromean noise, respectively. It is assumed that the purity measurement requires d-steps of analysis delay with d < N, and thus, the average purity for the kth period is available at t ) d in the k + 1th period. If N < d < 2N, the measurement can be available m in eq 25 should be modified to at the k + 2th period and yk+1 m yk+2, and other variables and equations associated with this change should be modified accordingly. In the above, N was fixed at 30 irrespective of T which can be varied by the optimizer. Optimizing Control System Overall Structure. The overall structure of the proposed optimizing control system is depicted in Figure 2. It is centered around the nonlinear discrete-time model in eq 25. The model is updated offline when needed and is used for optimization, construction of a linearized model for RMPC, and generation of state estimate trajectory for linearization. The optimizer

|

|

k-1

(27)

Detailed formulas of A, B, and H are given in Appendix B. In the above, w(t) and n are assumed to be zero-mean white noises with covariance matrices Rw and Rn, respectively. Error from linearization as well as the effect of other unmodeled disturbances are represented by the sum of the error carried over from the previous switching period plus a new random error that occurs during the kth switching period.26 In the differenced model, the former gets canceled through differencing and only the latter part remains, which is reasonably modeled by the zeromean white noise sequence of w(t). For linearization in eq 27, the state estimates should be available, and thus, a nonlinear observer was built on eq 25. Because the purity measurement is available only once during an entire switching period with a measurement delay, the state estimation for the (k - 1)th period is performed offline before the start of the kth period according to the following equation and the estimated trajectories are used for the linearization: m j (yk-1 xk-1|k-1(0|d) ) xk-1|k-2(0|d) + K (d) - H (xk-1|k-2(0|d))) xk-1|k-1(t + 1|d) ) A(uk(t), vk(t), xk-1|k-1(t|d)) + B(uk(t), vk(t)),

t ) 0, ..., N - 1

(28)


where is the measurement of the product purity at t ) 0 of the (k - 1)th period, which is obtained at t ) d due to the analysis delay. xk1|k2(t1|t2) denotes an estimate of xk1(t1) based on the information available up to (k2, t2). RMPC. In the present SMB process, it is assumed that only one measurement is taken for each of extract and raffinate streams, respectively, after each period of operation and the input sequence for the kth period is determined based on the average purities at the end of the k - 1th period. Because of the analysis delay, however, it is assumed that the input sequence for the kth period can be computed at t ) d in the kth period, after the measurement result is available, and hence, the part of the input sequence corresponding to t ) d to t ) N - 1 is implemented and the rest that corresponds to t ) 0 to t ) d - 1 is carried over to the (k + 1)th period. The overall computation steps in RMPC can be described as follows: (1) At (k, d), the measurement ymk (d) is available. On the basis of this value, obtain ∆xk|k(0|d) using the Kalman filter (actually a smoother) in eq 31. (2) Using ∆xk|k(0|d) and the output predictor built with eq 31, compute the control input sequence {∆uk(0), ..., ∆uk(N - 1)} by solving eq 32. (3) Apply {∆uk(d), ..., ∆uk(N - 1)} to the process in the present kth period and carry over {∆uk(0), ..., ∆uk(d - 1)} to the (k + 1)th period as {∆uk+1(0), ..., ∆uk+1(d - 1)}. For construction of the RMPC law, we first perform “lifting” on eq 26 over a switching period and obtain the following model that describes the state transition from the present period to the next: m yk-1 (d)

∆xk+1(0) ) Φk∆xk(0) + Γk∆Uk + ξk yk(N) ) yk-1(N) + Πk∆xk(0) + Gk∆Uk + Hk-1ξk m yk+1 (d) ) yk(N) + n (29) where ∆Uk } [∆uk(0)T · · · ∆uk(N - 1)T]T. Definitions of associated matrices and noise terms in eq 29 are given in Appendix C. Let r be the purity set point and define ek(N) } yk(N) - r. Then, eq 29 can be rearranged to a standard stochastic state space form:

[

] [ ][

∆xk+1(0) ek(N)

] [ ] [ ]

[]

Φk0 ∆xk(0) Γk I + ∆Uk + ξ Gk ΠkI ek-1(N) H k ∆xk(0) ek(N) ) [Πk I ] + Gk∆Uk + Hξk ek-1(N) m ek+1(d) ) ek(N) + n (30) )

The Kalman filter to estimate ∆xk(0|d) is given by the following equation:

.

[

[

∆xk|k(0|d) ek|k-1(N|d)

∆xk+1|k(0|d) ek|k(N|d)

] [ ] ] [ ][ ] [ ] )

∆xk|k-1(0|d) + K(em k (d) - ek-1|k-1(N|d)) ek-1|k-1(N|d)

)

Φk0 ΠkI

∆xk|k(0|d) Γk + ∆Uk Gk ek|k-1(N|d)

(31)

The input sequence ∆Uk is determined at (k, d) by solving min ∆Uk

1 2

{∑

mp-1

i)0

mc-1

||ek+i|k(N|d)||Q2 +

∑ ||∆U

k+j||R

j)0

2

}

(32)

subject to input inequality constraints imposed by the pressure drop limitation and minimum flow rate for column wetting. In

3693

the above, mp and mc are prediction and control horizons, respectively; ek+i|k(N|d) can be represented as a linear function of ∆xk|k(0|d) and ∆Uk+j, j < i using eq 30. The RMPC law in the above has been derived by assuming that the computation time for eq 32 is shorter than one sampling interval. If the computation time is significant (in the order of one switching period), ∆Uk has to be computed m (d) (most SMB after the (k - 1)th period on the basis of yk-1 processes can have a pause period after a column switching). In this case, the algorithms in eqs 27, 31, and 32 need to be modified accordingly. Input Blocking. The number of sample units in each period, N, should be large enough for precise estimate of the average product purities. A large N, however, results in a high dimension for ∆Uk and consequently increases the computational burden. In reality, it is not necessary to change the MV’s at every sampling instant. The number of input changes, to be denoted as Nc, can be restricted to a smaller value than N using a blocking matrix.27 For demonstration, if N ) 3 and Nc ) 2 and we want to change u(t) at the first and third sample times but not at the second, the blocking matrix is designed as

[ ]

[ ]

1 0 1 0 ¯ k ) BU ¯ k f ∆Uk ) 0 0 ∆U ¯k Uk ) 1 0 U 0 1 0 1

(33)

Inserting the above into eqs 30-32 modifies the minimization j k instead of ∆Uk. problem to determine ∆U Model Parameter Estimator. The model parameter estimator updates key process parameters so that the sum of squared error of the purity estimate over a window of cycles DK is minimized such that min θ

∑ ||y

m k+1(d)

- H (xk|k(N|d;θ))||2

(34)

k∈DK

The purity estimate is obtained using the SMB model in eq 25 and the state estimator in eq 28. θ denotes the parameter vector to estimate, which preferably includes the adsorption equilibrium constants. Since RMPC has a capability to overcome model uncertainty to some extent, it is enough to invoke the parameter estimator only when there is a significant process change and/ or before the optimizer is activated. If there is no sudden and significant process change that cannot be overcome by RMPC, optimization will be conducted intermittently with a long interval and a sufficiently large amount of measurements can be collected before each optimization. Under this situation, the estimator can provide accurate parameter estimate unless the measurements are not informative and/ or insensitive to parameter changes. For this reason, the parameter estimator was not employed in the numerical study under the assumption that the estimator can work accurately. Optimizer. The existing SMB optimizing control methods10,15–18 are designed to carry out a minimization of an objective function composed of the desorbent consumption minus the product yield while enforcing the purity requirement by imposing inequality constraints on the (predicted) product purities. As the optimization is repeated online, the cyclic steady state (CSS) optimum can be approached unless the optimizing controller is unstable or there are persistent disturbances that cannot be eliminated completely with the existing actuators. In this research, an optimizer that conducts a CSS optimization, as in the work of Klatt et al.,9 on the basis of the nonlinear SMB model in eq 25, is considered. The most important difference from the Klatt et al.’s approach9 is that

3694


in our formulation the flow rates are allowed to vary within a switching period. The optimization problem is formulated as N-1

min

xm(0),T,s,v(τ),u(τ),τ)t1,...tNo

subject to

∑

1 J) (a1QD(t) - a2QF(t)) + sTMs 2 t)0 xm+1(0) ) xm(0) H (x(N)) g r - s, s g 0 eq 25 Qj,min e Qj e Qj,max, j ) I, ..., IV, F, D

(35)

Table 1. Parameters for SMB Process and Controller for Numerical Study Model Parameters and Nominal Operating Conditions L ) 19.7 cm, S ) 4.91 cm2, ε ) 0.47 collocation points: {0, 1/nc, 2/nc, ..., 1, 5} nc ) 25 for numerical process, nc ) 10 for nominal model DA, DB ) 1 cm2/min HA ) 1.286, HB ) 0.5667, kA ) 1 min-1, kB ) 5 min-1 cAF ) cBF ) 0.03 g/cm3 QF ) 0.7 cm/min, QD ) 2 cm/min switching period, T ) 30 min, analysis delay ) 25 min SMB configuration 1-1-1-1 (4 columns) Nominal Values of Controller and Optimizer Parameters N ) 30, h ) T/N ) 1 min, Nc (for RMPC) ) 6 Q ) I, R ) 0.02I, Rw ) 10-4I, Rn ) 10-7I, M ) 106I 0 e QI, QII, QIII, QIV e 10 a1 ) 0.25, a2 ) 0.5

In the above, the subscript m denotes a dummy index for the switching period and hence xm+1 ) xm expresses the CSS condition, u and v represent [ (QIT/εV) (QIIT/εV)]T, which is also used as an MV for RMPC, and [ (QFT/εV) (QDT/εV)]T, respectively (see eq 25), No is the number of flow rate changes within a switching period, s contains the slack variables for the purity requirement, M is a positive definite matrix penalizing the magnitude of the slack variables. The number of sampling steps for analysis delay is adjusted whenever T is updated by the optimizer. In finding x(0) that satisfies the CSS condition, we did not employ any active search method but just repeated the operation until a CSS is reached, which has been referred to as the nested/single discretization technique in Kawajiri and Biegler.28 The optimizer in this study relies on the nonlinear SMB model, the parameters of which are assumed to be updated based on the process measurements. It is sufficient to invoke the optimizer only when an appreciable change occurs in the SMB model. Once the optimizer determines the decision variables, only T and v(t) are implemented in the process, whereas u(t) is left to be manipulated by the lower level controller. Regarding the above way of implementing the decision variables, the following question may arise: Does CSS u(t) by RMPC coincide with u(t) from the optimizer? Before answering to this question, one needs to note from eq 35 that only v(t) and T determine minimum J whereas u(t) satisfies the constraints. Therefore, if RMPC attains the purity regulation objective with the optimal v(t) and T implemented in the process, the optimization objective is fulfilled as well. The optimization was solved using fmincon.m in the optimization toolbox in MATLAB. Numerical Study In Table 1, parameter values for the process and controller and the nominal operating conditions used in the numerical study are given. All the nominal operating conditions except cA,F and cB,F were allowed to be modified during the optimization. It is first shown how reasonably the CSCM with the FSBC can represent the dynamics of a single SMB column in comparison with the OCM. Next, the tracking and regulation performance as well as the robustness of the proposed RMPC technique are examined for various disturbance and set point change scenarios. To demonstrate the improvement from using the successively linearized model, a fixed linear time invariant model-based RMPC technique was also applied to

some selected scenarios and the performance is compared. Finally, the combined optimization and RMPC was simulated. The optimization was performed for various cases of decision variables from {v(t), u(t)} with No ) 1, i.e., no flow rate change within a switching operation, to {v(t), u(t), T} with No ) 3, i.e., up to three instances of flow rate changes within a switching operation. Emphasis was given on examining how much the optimal value of J is affected by the degree of freedom in the decision variables and also how closely RMPC achieves the required purity regulation under optimization. Results and Discussion Performance of the CSCM/FSBC column model. In Figure 3a and b, we exhibit the response of the ODE column model derived using the CSCM and FSBC with nc ) 6. For comparison, we also show in Figure 3c and d a similar response of the column model derived using the OCM. It is assumed that only the desorbent flows initially through the column and cA are increased stepwise at the column inlet at t ) 0. It can be observed that the CSCM model produces physically more reasonable responses than the OCM model. For example, the unrealistic ripples in the response curves are effectively suppressed and the artificial zero-slope at the column outlet in Figure 3d is removed. Open-Loop Behavior of the SMB Model. Open-loop responses of the numerical SMB process are demonstrated in Figure 4a by changing the flow rates for three different cases, which are also located in the Mazzotti’s m2-m3 and m1-m4 planes7 as in Figure 4b marked with the pure separation and complete regeneration regions. mj is defined as mj } (QjT - Vε)/(V(1 - ε))

(36)

One can see that the purities become higher and higher as the flow rates approach and eventually enter the region of complete separation and regeneration. However, perfect separation is not realized with the numerical process due to the dispersion effect and column switching that are not considered in the Mazzotti’s analysis. Effectiveness of the Successively Linearized Model-Based RMPC Method. This time, we assess the performance of the successively linearized model-based RMPC method relative to that of the fixed linear model-based RMPC method for the case of set point changes. The results are depicted in Figure 5a and b. Figure 5a and b shows the tracking performance of the fixed linear model-based RMPC method. Here, linearized SMB models are obtained at the respective initial steady states. As


3695

Figure 3. Step responses of ODE column models by (a and b) CSCM with FSBC (nc ) 6, zL ) 5) and (c and d) OCM (nc ) 7). For both models, u ) 3.52, κ ) 200, η ) 1.13, λ ) 0.103, and HA ) 1.286.

Figure 4. Open-loop responses of average purity for different flow rates: (a) flow rates and average purities; (b) flow rate changes in the m2-m3 and m1-m4 planes.

can be observed, the set point tracking succeeds when the set point change is small but fails when it is large. This suggests that the region within which a linearized SMB model is valid (and therefore RMPC is effective) is limited. Figure 5c shows the tracking performance of the successively linearized model-based RMPC technique. The linearization is performed after each switching operation and thus, the linearized model can represent the process dynamics over the next

switching period with a reasonable accuracy. Consequently, the tracking of a large set point change that was unsuccessful by the fixed model-based RMPC method is successfully achieved as is shown in Figure 5c. Regulation against the Changes in Equilibrium Constants and Bed Porosity. It is generally agreed that the adsorption isotherm is the most critical and sensitive part of an SMB model. The existing isotherm equations have a fundamental limitation in correctly representing the reality, and the associated parameters also can vary with column aging and bed temperature. Bed porosity is another important factor that varies with time as the packing state of the bed gradually changes. Such changes can be traced using an online parameter estimation technique. Despite this, it is important that the controller possess the robustness needed to cope with such parameter changes. Figure 6 shows the RMPC performance when both HA and HB increase stepwise by 15% from their nominal values at k ) 75. When the system is left uncontrolled at the change, the average product purities drop from 99.5% to around 91% for both extract and raffinate. Under RMPC, the purities are recovered to 99.5% after some transient. Figure 7 exhibits a result of RMPC when the bed porosity decreases from 0.47 to 0.37 stepwise at k ) 200. After a somewhat long transient, RMPC helps the purities recover to their original values. In this study, it is assumed that the purity measurements are available only once for each switching period. If more frequent measurements were assumed, the transient periods under RMPC in the above simulations would have been reduced further. Performance of Combined Optimization and RMPC. The optimizer relies only on the SMB model, the parameters of which may be updated using the process measurements. It is enough to invoke the optimizer only at the process startup and when there is a change in the process, purity specifications, and/ or cost factors.

3696


Figure 5. Set point tracking by RMPC constructed based on (a and b) fixed linear models and (c) a successively linearized model.

Figure 6. Performance of RMPC against a change in adsorption equilibrium constants: HA ) 1.286 f 1.479, HB ) 0.5667 f 0.651.

Figure 7. Performance of RMPC against a change in bed porosity; ε ) 0.47 f 0.37.

We first considered the case where the decision variables for the optimization are time-invariant u and v, and x(0). From the optimal solution, only v is implemented on the process and u is left to be redetermined by RMPC. The results are summarized in Figure 8. One can see that the cost function J is unchanged under RMPC before the optimization, but evidently decreased after the optimization. The purity perturbation by the change in v is gracefully overcome by RMPC. In Figure 8b, the flow rate changes are shown on the m2-m3 and m1-m4 planes. For this, the time-varying flow rates were averaged over each switching period and converted to mi values. It can be seen that the flow rates enter the complete separation and regeneration regions (p1 f p2) as RMPC starts to work and tend to the vertices (p2 f p3),

recognized as the optimum condition of an ideal SMB process, after the optimization. In addition to the simplest optimization case shown in Figure 8, we considered five more cases by increasing the degree of freedom in the decision variables. The cases included timeinvariant v and u plus T, time-varying v(t) and u(t) with No ) 2 and 3 with and without T, respectively, where No refers to the number of flow rate changes within a switching period. Figure 9a exhibits the objective function values for the different choices of the decision variables. One can obviously see that the time-varying flow rates, i.e., PowerFeed operation, reduces J more than the time-invariant flow rates. In addition, it appears that J can be reduced considerably further when the switching period is optimized together with the flow rates. In Figure 9b, optimized v(t) with No ) 3 and also T are exhibited. One can see that the optimized T is 17.35 min which is smaller than the analysis delay d ) 25 min. As we implemented the new T, therefore, the period indices for ym and em in eqs 26 and 28-31 were modified to reflect this change. Computation for optimization of cases a and f in Figure 9 required approximately 3 and 20 min, respectively, on a PC with Intel 6600 quad-core processor. Indeed, the computation time of this order should not cause any practical problem unless it is too long since the optimization is intermittently carried out offline but could be shortened if a CSS x(0) were actively sought using the methods of Kawajiri and Biegler28 or others. Use of a Nominal Model with Bias Correction. For comparison with the proposed optimizing control scheme, optimizing control was also attempted using an SMB model with parameter error but with bias correction to compensate for the prediction error of the purity as in the standard MPC approach. We revisited the situation in Figure 6 where HA and HB of the SMB process vary from 1.286 and 0.5667 to 1.479 and 0.651, respectively. It was assumed that the nominal model is given as xk(t + 1) ) A(uk(t), vk(t), xk(t);θold) + B(uk(t), vk(t);θold) yk(N) ) H (xk(N);θold) + b (37)


3697

Figure 8. Performance of combined optimization and RMPC, where the optimizer searches for optimum time-invariant flow rates: (a) profiles of average purity, cost function, and [QD QF] under the nominal state, RMPC, and combined optimization and RMPC; (b) flow rate plots in the in the m2-m3 and m1-m4 planes.

Figure 9. Performance of combined optimization and RMPC depending on the choice of decision variables: (a) improvement of objective function; (b) optimized time-varying flow rates with No ) 3 and duration of switching period.

where θ )[HA HB]T and θold ) [1.268 0.5667]T; b denotes the bias correction term obtained as an average value of the prediction error of the purity over the period up to 110 after the parameter change in Figure 6.

In conclusion, the above optimizing control method was found to either succeed or fail depending on the operating point whereas the proposed optimizing control method with accurate parameter values always succeeded in all the test

3698


Acknowledgment This work was supported by Energy Resource R&D program(200931022) under the Ministry of Knowledge Economy, Republic of Korea and also by KOSEF through the Center of Advanced Bioseparation Technology in Inha University. Appendix A: Derivation of Equation 16 Figure 10. Result of combined optimization and RMPC based on a nominal model with parameter error but with bias correction instead. T is not optimized and No ) 1 is chosen.

cases. It was observed that the possibility of failure increases as the purity target as well as a1, the unit cost of QD, in eq 35 are increased. Figure 10 shows a case of failure. For the same case, however, optimizing control correctly worked when the purity set point was changed to a lower value than around 0.99. The reason for the failure is simply the inappropriateness of the bias term at higher purities than where it was obtained. Conclusions Through this study, a new bilevel optimizing control structure for the SMB process has been proposed. The proposed structure is composed of an offline optimizer that determines the optimal time-varying flow rates and duration of switching period under a CSS condition and a repetitive controller called RMPC that conducts the regulation of extract and raffinate purities. The optimizer relies only on a nonlinear SMB model, which is assumed to be continuously updated using the process measurements. Hence, only RMPC is run in real-time and the optimizer is invoked only when there is a change in the estimated model, the purity specifications, and/or the cost factors in the objective function. Both the controller and the optimizer were constructed on the basis of a discrete-time nonlinear SMB model that was reduced from the PDE SMB model by applying the CSCM with the FSBC and exact discretization. For RMPC, the model was successively linearized after each switching period along the operating trajectories of the previous period. Numerical studies showed that a fixed model-based RMPC fails to track a large set point change, which suggests strong nonlinearity of the SMB process even with linear adsorption isotherms. On the other hand, the proposed successive linearization-based RMPC copes well with the process nonlinearity and successfully fulfilled the tracking of a large set point change, an important capability especially during a start-up operation. The optimizer was designed to minimize an objective function that consists of desorbent consumption and product yield under a CSS condition. Various choices of decision variables were considered, from time-invariant flow rates to time-varying flow rates plus the switching period. The optimization was attempted for up to three instances of flow rate changes per period plus the switching period itself and we found the objective function decreased noticeably as the degree of freedom in the decision variables was increased. It is believed that the proposed optimizing control scheme has advantages over the existing ones in terms of both the optimizer and the controller performances and also the way it is implemented, real-time control, and offline optimization. In future works, the proposed scheme will be extended to nonlinear isotherm systems and tested on an experimental SMB process.

Numerical Process j and B j , which are (n + 2) × (n + 2) Let us partition A matrices, such that

c r c where a0c, an+1 , an+1 , b0c, and bn+1 are n × 1 vectors; h0 and hn+1 are scalars; and A and B are n × n matrices. Rewriting eqs 14 and 15 after discarding the ODE’s at the boundary points gives

dc dw ) (λB - uA)c + (λbc0 - uac0)c(z0) + +η dt dt c c (λbn+1 - uan+1 )c(zn+1) (39) dw ) κ(Hc - w) dt

(40)

The far-side BC results in 0)

dc(zn+1) rT ) h0c(z0) + an+1 c + hn+1c(zn+1) dz

(41)

Eliminating c(zn+1) from eq 39 using eq 41 and substituting eq 40 into eq 39 yields dc ) Epc + ηκw + dpc(z0) dt

(42)

where Ep } λB - uA - ηκHI dp } λbc0 - uac0 -

1 c rT (λbc - uan+1 )an+1 hn+1 n+1

h0 c (λbc - uan+1 ) hn+1 n+1

(43) Nominal Model In the nominal model, adsorption equilibrium is assumed, thus eq 15 is replaced by w j ) Hcj. The model equation in eq 14 then becomes dc¯ 1 ¯ )c¯ (λB¯ - uA ) dt 1 + ηH

(44)

By the same procedure as in eq 38 to eq 41, we have dc ) Emc + dmc(z0) dt where

(45)


(

1 1 c rT λB - uA (λbc - uan+1 )an+1 1 + ηH hn+1 n+1 h0 1 c dm } λbc0 - uac0 (λbc - uan+1 ) 1 + ηH hn+1 n+1

Em }

(

)

[ ]

) (46)

Appendix B: A, B, and H in Equation 27

F 0 0 ∂QEcE Fz I 0 ∂A ∂ QEcE ) ∂zT A) T ) ∂x ∂[zTsETsRT] QRcR ∂QRcR 0 I ∂zT

[ ]

(47)

In ∂QEcE/∂zT ∈ R4×8nc and ∂QRcR/∂zT ∈ R4×8nc, only the following four terms are nonzero: ∂QEcBE ∂QEcAE ) QE ) Q E, ∂cAE ∂cBE ∂QRcBR ∂QRcAR ) Q R, ) QR ∂cAR ∂cBR ∂(A + B) ∂ B) ) T ∂uT ∂u

[ ]

Fz + G 0 ) 0

(48)

[ ] ∂(Fz + G) ∂uT 0 0

(49)

where ∂(Fz + G)/∂uT is computed numerically. In H ) ∂H /∂xT ∈ R2×(4nc+4), only the following terms are nonzero: ∂pjE ∂pjE sAE sBE )) , 2 ∂sAE ∂sBE (sAE + sBE)2 (sAE + sBE) sBR ∂pjR sAR ∂pjR ), ) 2 ∂sAR ∂sBR (sAR + sBR) (sAR + sBR)2

(50) Appendix C: Definition of Matrices in Equation 30 N-1

Φk } P . Γk } P[

∏A

k-1(t)

) PΦ* k

t)0

N-1

∏

N-1

Ak-1(t)Bk-1(0)

t)1

ξk } P[

∏A

N-1

k-1(t)

t)1

Literature Cited

k-1(t)Bk-1(1)...Bk-1(N

- 1)] ) PΓ* k

t)2

. Πk } Hk-1Φ* k-1 . Gk } Hk-1Γ* k N-1

∏A

∏ t)2

[ ]

w(0) w(1) Ak-1(t) · · · I] l w(N - 1)

(1) Ruthven, D. M.; Ching, C. B. Conunter current and simulated counter current adsorption separation processes. Chem. Eng. Sci. 1989, 44, 1011– 1038. (2) Broughton, D.; Gerhold, C. Continuous sorption process employing fixed bed of sorbent and moVing inlets and outlets. U.S. Patent 2,985,589, 1961. (3) Ludemann-Hombourger, O.; Nicoud, R. M.; Bailly, M. The “VARICOL” process: A new multicolumn continuous chromatographic process. Sep. Sci. Technol. 2000, 35, 1829–1862. (4) Abel, S.; Mazzotti, M.; Morbidelli, M. Solvent gradient operation of simulated moving beds: I. Linear isotherms. J. Chromatogr. A 2002, 944, 23–39.

3699

(5) Schramm, H.; Kaspereit, M.; Kienle, A.; Seidel-Morgenstern, A. Simulated moving bed process with cyclic modulation of the feed concentration. J. Chromatogr. A 2003, 1006, 77–86. (6) Zhang, Z.; Mazzotti, M.; Morbidelli, M. PowerFeed operation of simulated moving bed units: Changing flow-rates during the switching interval. J. Chromatogr. A 2003, 1006, 87–99. (7) Mazzotti, M.; Storti, G.; Morbidelli, M. Optimal operation of simulated moving bed units for nonlinear chromatographic separations. J. Chromatogr. A 1997, 769, 3–24. (8) Kloppenburg, E.; Gilles, E. Automatic control of the simulated moving bed process for C8 aromatics separation using asymptotically exact input/output linearization. J. Process Control 1999, 9, 41–50. (9) Klatt, K.-U.; Hanisch, F.; Du¨nnebier, G. Model-based control of a simulated moving bed chromatographic process for the separation of fructose and glucose. J. Process Control 2002, 12, 203–219. (10) Toumi, A.; Engell, S. Optimization-based control of a reactive simulated moving bed process for glucose isomerization. Chem. Eng. Sci. 2004, 59, 3777–3792. (11) Hara, S.; Yamamoto, Y.; Omata, T.; Nakano, N. Repetitive control system: A new type servo system for periodic exogeneous signals. IEEE Trans. A. C. 1998, 33, 659–668. (12) Tomizuka, M.; Tsao, T.; Chew, K. Analysis and synthesis of discrete-time repetitive controllers. Trans. ASME, J. Dyn. Sys. Meas. Control 1988, 111, 353–358. (13) Lee, J. H.; Natarajan, S.; Lee, K. S. A model-based predictive control approach to repetitive control of continuous processes with periodic operations. J. Process Control 2001, 11, 195–207. (14) Natarajan, S.; Lee, J. H. Repetitive model predictive control applied to a simulated moving bed chromatography system. Comput. Chem. Eng. 2000, 24, 1127–1133. (15) Erdem, G.; Abel, S.; Morari, M.; Mazzotti, M.; Morbidelli, M.; Lee, J. H. Automatic control of simulated moving beds. Ind. Eng. Chem. Res. 2004, 43, 405–421. (16) Abel, S.; Mazzotti, M.; Erdem, G.; Morari, M.; Morbidelli, M. Optimizing control of simulated moving beds: linear isotherm. J. Chromatogr. A 2004, 1033, 229–239. (17) Abel, S.; Erdem, G.; Mazzotti, M.; Morari, M.; Morbidelli, M. Optimizing control of simulated moving beds - experimental implementation. J. Chromatogr. A 2005, 1092, 2–16. (18) Amanullah, M.; Grossmann, C.; Mazzotti, M.; Morari, M.; Morbidelli, M. Experimental implementation of automatic ‘cycle to cycle’ control of a chiral simulated moving bed separation. J. Chromatogr. A. 2007, 1165, 100–108. (19) Grossmann, C.; Erdem, G.; Morari, M.; Amanullah, M.; Mazzotti, M.; Morbidelli, M. ‘Cycle to cycle’ optimizing control of simulated moving beds. AIChE J. 2008, 54, 194–208. (20) Guiochon, G.; Shirazi, S. G.; Katti, A. Fundamentals of preparatiVe and nonlinear chromatography; Academic Press: Boston, 1994. (21) Yun, W.; Lee, K. S. The use of cubic spline and far-side boundary condition for the collocation solution of a transient convection-diffusion problem. Korean J. Chem. Eng. 2007, 24, 204–208. (22) Choi, S. Y.; Kim, H. S.; Lee, K. S.; Yoo, K. P.; Lee, W. H. Reduced-order distillation model using collocation method with cubic splines. Korean J. Chem. Eng. 1991, 8, 44–52. (23) Manickam, A. V.; Moudgalya, K. M.; Pani, A. K. Second-order splitting combined with orthogonal cubic spline collocation method for the Kuramoto-Sivashinsky equation. Comp. Math. Appl. 1998, 35, 5–25. (24) Du¨nnebier, G.; Klatt, K.-U. Modelling and simulation of nonlinear chromatographic separation processes: a comparison of different modelling approaches. Chem. Eng. Sci. 2000, 55, 373–380. (25) Åstro¨m, K. J.; Wittenmark, B. Computer-Controlled Systems: Theory and Design, 3rd ed.; Prentice-Hall: New York, 1996. (26) Lee, J. H.; Lee, K. S.; Kim, W. C. Model-based iterative learning control with a quadratic criterion for time-varying linear systems. Automatica 2000, 36, 641–657. (27) Lee, J. H.; Chikkula, Y.; Yu, Z.; Kantor, J. Improving computational efficiency of a model predictive control algorithm for multi-time-scale systems using wavelet transformation. Int. J. Control 1995, 61, 859–883. (28) Kawajiri, Y.; Biegler, L. T. Optimization strategies for simulated moving bed and powerfeed processes. AIChE J. 2006, 52, 1343–1350.

ReceiVed for reView July 12, 2009 ReVised manuscript receiVed February 20, 2010 Accepted February 24, 2010 IE901121Y

Bilevel Optimizing Control Structure for a Simulated Moving Bed

Recommend Documents