Automatic Control of Simulated Moving Beds - ACS Publications

Ind. Eng. Chem. Res. 2004, 43, 405-421

405

Automatic Control of Simulated Moving Beds Gu 1 ltekin Erdem,† Stefanie Abel,‡ Manfred Morari,*,† Marco Mazzotti,‡ Massimo Morbidelli,§ and Jay H. Lee| Automatic Control Laboratory, Physikstrasse 3, and Institute of Process Engineering, Sonneggstrasse 3, ETH-Zu¨ rich, 8092 Zu¨ rich, Switzerland, Institute for Chemical and Bio-Engineering, ETH-Zu¨ rich, 8093 Zu¨ rich, Switzerland, and Integrated Sensing System Identification and Control Laboratory, Georgia Institute of Technology, 778 Atlantic Drive, Atlanta, Georgia 30332-0100

In recent years, simulated moving bed (SMB) technology has become increasingly attractive for complex separation tasks and emerged as the leading chromatographic separation technique in the fields of pharmaceuticals, fine chemicals, and biotechnology. Extensive usage of SMB technology has created a need for robust automatic control techniques to exploit its full economic potential. Automatic control of SMB units is a challenging task not only because of its nonsteady-state, nonlinear, mixed discrete, and continuous nature, but also because of long delays involved in exhibiting the effect of disturbances. This work proposes a new optimization-based control strategy, in which a linearized time-varying reduced-order model that accounts for the periodic nature of the SMB process is used for on-line optimization and control of SMBs applied to systems characterized by linear adsorption isotherms. Four internal flow rates, which can be adjusted via external flow rates, are used as the manipulated variables. On-line concentration measurements at the product outlets together with a periodic Kalman filter are used to reduce the effect of model errors. Optimal input adjustments that allow for the achievement of process specifications and optimal performance are calculated on the basis of the predicted evolution of the plant. The realization and implementation of the concept on a virtual eight-column SMB unit is described, and the capability of the controller to run the SMB plant at its economic optimum regardless of possible disturbances and model uncertainties is assessed. 1. Introduction Simulated moving bed (SMB) chromatography is a powerful technique for complex separation tasks. It has become an attractive technology especially in the areas of biotechnology, fine chemicals, and pharmaceuticals1-3 because it provides significant advantages over conventional separation techniques such as batch chromatography. SMB technology overcomes the limitations of batch chromatography and provides high productivity through overloaded operating conditions. Moreover, it reduces solvent consumption and product dilution, which leads to an easier and cheaper product recovery step. SMB technology attracts interest particularly in the area of chiral separations, e.g., for enantiopure drug development, where the resolution of enantiomers is usually a binary separation problem with low selectivities and high costs of solvent and chiral stationary phases. Market competition requires the pharmaceutical companies to work with short drug development times, where they require only a few grams of the chiral drug for preliminary biological tests, and to look for costefficient conditions for large-scale production. SMB technology satisfies the needs both for the early testing and production stages.4 The underlying principle of SMB technology is a continuous simulated countercurrent chromatographic * To whom correspondence should be addressed. E-Mail: [email protected]. Tel.: +41 1 632 76 26. Fax: +41 1 632 12 11. † Automatic Control Laboratory, ETH-Zu¨rich. ‡ Institute of Process Engineering, ETH-Zu¨rich. § Institute for Chemical and Bio-Engineering, ETH-Zu¨rich. | Integrated Sensing System Identification and Control Laboratory, Georgia Institute of Technology.

Figure 1. Scheme of a simulated moving bed (SMB) unit. The dash arrow indicates the inlet/outlet positions after the first switch.

process, and the key idea of the SMB is to simulate the solid-phase motion of the corresponding true moving bed (TMB) unit by using standard fixed-bed chromatographic columns and periodically switching the inlet and outlet ports of the unit in the same direction as the fluid flow (see Figure 1). As a result, SMB technology profits from the countercurrent contact of the stationary and mobile phases without facing the difficulties associated with the movement of the solid phase. The most widely adopted SMB configuration includes four sections, with each section performing a specific task in the separation process. With reference to the separation of a binary

10.1021/ie030377o CCC: $27.50 © 2004 American Chemical Society Published on Web 11/12/2003

406

Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004

Figure 2. Concentration profile of the light component (B) in the raffinate outlet of an SMB unit.

mixture, the feed stream, carrying a mixture of the more adsorbable component A and the less adsorbable component B, is introduced between sections II and III. The separation is carried out in the two central sections, where component A is adsorbed and carried with the solid phase to the extract outlet and component B is desorbed and carried with the liquid phase to the raffinate outlet. Component B is adsorbed in section IV so that the eluent is regenerated before being mixed with fresh eluent and recycled to section I, where the adsorbent is regenerated through the desorption of component A. It is common practice to divide each section of an SMB unit into several subsections in order to closely approximate true countercurrent movement of the solid phase (see Figure 1). Detailed descriptions of the principles of SMB units can be found elsewhere.1 Essential for an understanding of the challenges in the automatic control problem is that the SMB process will never reach a steady state with constant profiles of all of the process variables. The stationary regime of an SMB unit is a cyclic steady state in which the concentration profile along the unit moves with the liquid flow and is taken back periodically with each switch. As a result, under the assumption that all of the columns comprising the unit are identical, SMB units exhibit the same time-dependent behavior during each time period between two successive switches of the inlet and outlet ports. Therefore, modeling of the stationary regime of an SMB unit requires a time-dependent model. Figure 2 shows a typical time-varying concentration profile of component B at the raffinate outlet of an SMB unit. The modeling of SMB units has been widely studied by several research groups, and models with varying complexity are available for the design and optimization of SMB separations.5 An important contribution was made by the introduction of a shortcut design method, the so-called triangle theory,6 developed in the frame of equilibrium theory, that accounts for competitive adsorption thermodynamics but neglects mass-transfer resistance and axial dispersion effects. Triangle theory provides several advantages for the design and operation of SMB units. It allows for the selection of nearoptimal and robust operating conditions after a few algebraic computations for Langmuir adsorption equilibrium isotherms, i.e., those of greatest practical interest. Moreover, this approach provides a thorough understanding of how changes in the feed composition and flow rate affect the separation performance. As a result,

triangle theory is widely employed by SMB practitioners. On the other hand, detailed SMB models that account for column efficiency effects as well can be very valuable for detailed process design and optimization.7,8 The performance of both detailed and shortcut SMB design approaches depends on the quality of the physical data provided to the model, particularly the adsorption isotherms, which are difficult to measure with high precision and characterize the separation performance. Optimal operation of SMB units suffers from the fact that SMB separations are highly sensitive to disturbances, e.g., temperature or feed composition changes, and changes in operating parameters, e.g., chemical and mechanical aging of the stationary phase, which can result in suboptimal operating conditions or off-specification production. The natural solution to the above problems is the development and implementation of an automatic control strategy. The control strategy should exploit the full economic power of SMB technology by both operating the SMB plant at its optimal operating conditions and adapting conditions in an evolutionary manner to achieve the best possible separation regardless of the possible disturbances and changes in the operating parameters. This control task is challenging not only because of the non-steady-state, nonlinear, mixed discrete (inlet-outlet port switches), and continuous nature of the SMB units but also because of the long delays exhibited by the effects of disturbances. Recently, several approaches to the automatic control of SMBs have been introduced. These approaches can be grouped into two main categories. Control strategies within the first category make use of a TMB model to approximate the SMB dynamics in order to overcome the problems associated with the hybrid and non-steadystate nature of the SMB unit. These two underlying characteristics of the SMB unit are hidden from the controller.9,10 In the second category, an optimal operating trajectory of the process is calculated by off-line optimization, e.g., on the basis of detailed SMB models, and local SISO controllers are used to keep the process on the optimal operating trajectory, i.e., to maintain the position of the internal profiles, despite the disturbances. Changes in the operating parameters necessitate off-line reoptimization of the process to update the optimal trajectory and to resynthesize the local controllers according to the updated optimal trajectory.11 At the heart of all of these control strategies is the availability of accurate physical data on the system and/ or the assumption that the optimal operating conditions of the SMB are known a priori. However, there usually is no optimal operating trajectory fixed a priori. Optimal profiles depend on the operating parameters, e.g., feed concentration and throughput, and the physical parameters of the system, e.g., adsorption isotherm and singlecolumn properties such as column void fraction and column efficiency, which are subject to change. As a consequence, the automatic control of SMB units goes far beyond a simple regulation or tracking problem. It is posed here as a general dynamic constrained optimization problem. Model predictive control (MPC) has proven to be a powerful control technique in the past two decades. The main advantage of MPC is its ability to handle multivariable systems with time delays and to incorporate hard constraints on both the inputs and the states explicitly through on-line optimization.12,13 Recently, a new model predictive control method formulated for

Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 407

periodic process control problems has been introduced.14,15 This approach, repetitive model predictive control (RMPC), makes use of the key concepts of repetitive control (RC) and combines them with the main ideas of model predictive control. RMPC can be easily used for tracking or regulation problems, but its main significance comes from its flexible structure, which allows any economic objective to be used in the control problem.14,16,17 RMPC provides a quite general control strategy for SMB plants, and we believe that it should be considered as an alternative to other existing approaches. 2. SMB Unit One of the most widely adopted SMB schemes includes four sections consisting of chromatographic columns, each of which has a volume of V ) LAcr, where L is the length and Acr is the cross-sectional area. The cl volume of each section is Vj ) ncl j V, where nj is the number of columns in section j (j ) I, ..., IV), i.e., Vj ) 2V for all j for a 2-2-2-2 configuration. The SMB system considered in this study is a closed-loop foursection unit arranged in a 2-2-2-2 configuration (see Figure 1). For the sake of clarity, the virtual simulated plant will be referred to as the SMB plant throughout this work. A dynamic model of the SMB process is obtained by interconnecting the dynamic models of the single chromatographic columns. Single-column dynamic models of varying complexity are available in the literature.5 The equilibrium dispersive model (EDM) is regarded as a good compromise between model accuracy and computational efficiency.7 Therefore, a linear adsorption equilibrium together with the EDM is used to model the single-column dynamics in this study (please refer to the notation section for the meanings of the symbols)

h

/ ∂qi,h ∂ci,h ∂ci,h ∂2ci,h + (1 - h) + vh ) hDi ∂t ∂t ∂z ∂z2 / ) Hici,h (i ) A, B; h ) 1, ..., 8) qi,h

(1) (2)

In the above equations, vh represents the superficial velocity, which is the same for all columns belonging to the same section, i.e., for the given configuration in Figure 1, v1 ) v2 ) QI/Acr, v3 ) v4 ) QII/Acr, v5 ) v6 ) QIII/Acr, and v7 ) v8 ) QIV/Acr. Proper boundary conditions and proper initial conditions (enforcing the inlet/ outlet switching mechanism), together with the following node balances, complete the mathematical model

v1 ) v8 +

QI/O 1 Acr

in out ) v8ci,8 + v1ci,1 I/O Qh+1 vh+1 ) vh + Acr in vh+1ci,h+1

)

out vhci,h

(3)

QI/O 1 cI/O Acr i,1

for h ) 1, ..., 7

I/O Qh+1 + cI/O Acr i,h+1

(4)

(5)

for h ) 1, ..., 7 (6)

Here, QI/O identifies the flow rate of the inlet or the h

Figure 3. Scheme of the optimizing control concept.

outlet stream entering or leaving the SMB loop just in before column h, and cout i,h and ci,h+1 are the concentrations of component i at the outlet of column h and at the inlet of column h + 1, respectively. The SMB configuration in Figure 1 corresponds to the following I/O ) QD; specifications: QI/O h ) 0 for h ) 2, 4, 6, 8; Q1 I/O I/O I/O E F R Q3 ) -Q ; Q5 ) Q ; Q7 ) -Q . Moreover, we have I/O I/O out I/O F I/O out ci,1 ) cD i , ci,3 ) ci,2 , ci,5 ) ci , and ci,7 ) ci,6 . The switching operation is explicitly considered by shifting the boundary conditions, i.e., QI/O h ) 0 for h ) 1, 3, 5, 7; D; QI/O ) -QE; QI/O ) QF; QI/O ) -QR; and ) Q QI/O 2 4 6 8 I/O I/O out I/O F I/O out ci,2 ) cD i , ci,4 ) ci,3 , ci,6 ) ci , and ci,8 ) ci,7 after the first switching operation. The equilibrium dispersive model can be solved efficiently by performing a finite difference approximation of the PDE model and replacing the effect of the apparent axial dispersion term on the right-hand side of eq 1 with numerical dispersion as discussed elsewhere.7,18 The following physical parameters complete the characterization of the SMB: Each chromatographic column has a length and diameter of 10 and 1 cm, respectively. Differences in the single-column properties due to packing differences and velocity variations are neglected, and the nominal porosity value is ) 0.7. Apparent axial dispersion is defined such that each column has 100 theoretical stages with respect to each solute. The Henry’s constants in eq 2 characterizing the retention behavior of the two components to be separated are HA ) 4 and HB ) 2, which leads to a selectivity of S ) HA/HB ) 2, a measure of the difficulty of the separation task. Finally, the switching time, t*, is fixed a priori at 480 s. 3. Control Concept The control strategy proposed in this work is based on on-line optimization. The conceptual structure of the controller is depicted in Figure 3. An explicit model of the SMB unit to be controlled is employed to predict the future evolution of the plant. The performance of the plant is optimized over a future horizon according to the current state of the plant. The process and product specifications are considered explicitly as constraints in the optimization problem. The optimization yields an optimal control action sequence that is applied according to a receding horizon strategy, i.e., the first element of the calculated optimal control action sequence corresponding to the current time t is implemented, and the remaining calculated optimal inputs are discarded. A new optimization problem is solved at time t + 1 on the basis of the information coming from the new measurements collected from the plant. In this approach, the switch time, t*, is fixed a priori, and the manipulated variables are the four internal flow rates, i.e., QI, QII, QIII, and QIV, which are adjusted by acting

408


on QI itself and on the external flow rates QE, QF, and QR. The key variables characterizing the SMB performance in the frame of triangle theory are the flow rate ratios in the four sections, which are defined as:

mj )

Qjt* - V V(1 - )

(j ) I, II, III, IV)

(7)

Note that the flow ratios mj are linear functions of the internal flow rates Qj. Two on-line measurement systems are installed at the raffinate and extract outlets of the plant to monitor the compositions of the product streams that constitute the measured output variables. One of the most critical components of the control scheme is the SMB model that serves for the prediction of the future behavior of the plant. On-line optimization of the process on the basis of detailed dynamic SMB models, which lead to nonlinear programming problems, is difficult because of the needed modeling effort and the excessive computation time requirement. Therefore, a simplified SMB model that captures the key characteristics of the process is needed. The underlying idea of RMPC is that possible model prediction errors and the effect of periodwise constant disturbances are likely to repeat from cycle to cycle and can be estimated by using the available measurements of the plant output, which is then used to correct for the possible model prediction errors in the future cycles. A periodic Kalman filter is used to combine the measurement information and the model estimation in an optimal sense to correct for the model errors recursively.19 3.1. Simplified SMB Model. The SMB process exhibits periodic steady-state dynamics where the concentration profiles in the unit show a time-dependent but period-invariant behavior. Even though the periodicity of the concentration profiles can conceptually be defined over successive switches, in real applications, the profiles are periodic only from cycle to cycle, i.e., one cycle is completed after eight switches for the given system in Figure 1, because of the unavoidable variations in single-column properties and dead volumes in the system. A linear time-varying (LTV) model that captures the periodic nature of the process is employed. The proposed LTV model is extracted from the dynamic SMB model described in section 2, which comprises 16 partial differential equations (PDEs), i.e., one for each species and column given by eq 1 and the corresponding node balances (eqs 3-6). The finite difference method is used to convert the PDEs into ordinary differential equations (ODEs). Ten grid points along each column are defined, i.e., 80 grid points along the unit, for the finite difference approximation. This yields 160 ODEs, i.e., one for each component per grid point, describing the internal concentration profiles of the two species A and B as a function of time and position along the SMB unit. A cycle comprises eight successive switches where each switch corresponds to a different port (inputoutput) configuration as described in section 2. As a result, the set of equations, i.e., the ODEs together with the node balances, describing the system changes according to the eight different port configurations. The resultant ODEs have some nonlinear terms, i.e., the convective terms in eqs 1, 4, and 6, where concentration and flow rate are multiplied, and they can be linearized at their cyclic steady-state values. Because the periodic steady-state profile is time-varying, the duration of a

Figure 4. Reference concentration profiles of species A and B inside the SMB unit at different sample points, i.e., t ) 0, t ) 3 and t ) 7. 0 and O show the corresponding concentration values of components A and B, respectively, at grid points defined along the unit.

cycle is divided into N sample points (N ) 64 in our case), and the nonlinear equations are linearized at each sample point using the corresponding cyclic steady-state concentration values and flow rates. The reference profile used for linearization at time t, where t ) 0, ..., N - 1, corresponds to the internal concentration values cref i,g (t), where i ) A, B, and g ) 1, ..., 80, are the component and space index, respectively. The reference ) 4.0, mref operating point is chosen as mref I II ) 2.1, ref ref mIII ) 3.9, and mIV ) 2.1, and the cyclic steady-state profiles are obtained by simulating the detailed SMB dynamics with the model described in section 2. The propagation of the reference internal concentration profiles, i.e., the profile at sample point t, t + 1, ..., is shown in Figure 4. The linearized equations can be written in state-space form, and the discrete-time statespace equations can be obtained by using a zero-order hold on the inputs and a sample time equal to the duration between two successive linearization sample points

xk(n+1) ) A(n) xk(n) + B(n) uk(n) yk(n) ) C(n) xk(n)

for n ) 0, ..., N - 1

The transition from one cycle to the next can be written as

xk+1(0) ) xk(N)

(8)

where k and n are the cycle and time indices, respectively; x is the state vector with a dimension of 160, which is composed of the concentration values along the unit, i.e., ci,g(n) for i ) A, B and g ) 1, ..., 80; and u is the manipulated variable vector consisting of the internal flow rates in the four sections, i.e., Qj for j ) I, ..., IV. The concentrations of the two components in the raffinate, i.e., cRA(n), cRB(n), and the extract outlet, i.e., cEA(n), cEB(n), constitute the measured output vector y. Here, all of the variables, i.e., x, u, and y, are defined in terms of deviation variables from the reference values. For instance, the state vector consisting of the internal


concentration values is defined as x(n) ) ci,g(n) cref i,g (n). Similarly, the manipulated variable vector is defined as u(n) ) Q(n) - Qref(n), where Q(n) and Qref(n) are vectors comprising the internal flow rates, i.e., Q ) [QI QII QIII QIV]T. The resulting linear time-varying state-space model captures both the time-dependent cyclic steady state and the hybrid nature of the process explicitly through the time-varying state-space matrices and the cycle-to-cycle transition operator. A more detailed description of the modeling concept can be found elsewhere.16 3.1.1. Reduction of the Model Order. The order of the model leads to a significant computational load for on-line applications, but it can be reduced with available model reduction techniques. These techniques are applicable to time-invariant systems. Therefore, the LTV model is lifted by grouping the input and output values for one cycle, i.e., for time interval n ) 0-N, to obtain a time-invariant cycle-to-cycle description of the system as given below

xk+1(0) ) Φxk(0) + ΓUk

(9)

Yk ) Πxk(0) + GUk

(10)

Φ, Γ, Π, and G can be interpreted as the state transition, reversed controllability, observability, and impulse response coefficient matrices of the system, respectively. Equation 9 describes the transition from cycle k to k + 1.

Φ ) A(N - 1) A(N - 2) ‚‚‚ A(0)

[

C(0) C(1) A(0) Π) ·· · C(N - 1) A(N - 2) · · · A(0)

]

(11)

(13)

Yk and Uk are defined as

Yk } [ykT(0) ‚‚‚ ykT(N - 1)]T ∈ R(N×ny)

(15)

Uk } [ukT(0) ‚‚‚ ukT(N - 1)]T ∈ R(N×nu)

(16)

where ny and nu are the dimensions of the output and manipulated variable vectors, respectively, i.e., ny ) nu ) 4 in our case. The order of the model is reduced from 160 to 32 via balanced model reduction.20-22 Figure 5 shows a typical comparison of the reduced-order model output and the nonlinear detailed simulation model output.

Figure 5. Comparison of reduced-order model output and the nonlinear detailed dynamic simulation model output. O, Reduced model output; -, plant output. cEA (g/L) and cRB (g/L) are the concentrations of species A in the extract and species B in the raffinate, respectively. The operating point is mI ) 4.0, mII ) 2.08, mIII ) 3.93, and mIV ) 2.1.

3.1.2. Disturbance Modeling. One of the main tasks of the controller is to reduce the effect of disturbances on the process performance. Therefore, prediction of the combined overall effect of all possible disturbances on the plant output is a critical issue affecting the control performance. Generally, a description of the disturbance as white noise added to the deterministic models gives poor performance. One possibility is to characterize the disturbances first through available first-principles process models and then incorporate them into the deterministic part of the process model. On the other hand, if the physical sources for disturbances cannot be identified, their overall effect on the output can be modeled as a stochastic process. For SMB units, two types of disturbances can act on the system. The first type of disturbance is periodwise-persisting, e.g., changes in the feed concentration that are generally experienced for different feed batches, and their effects on the plant output tend to repeat from cycle to cycle. The second type of disturbance is specific to each period, and hence, their effects will be specific to a cycle and will not propagate to the subsequent cycles. RMPC is based on the idea that, because the possible model prediction errors or the overall effects of periodwise-persisting disturbances on the plant output repeat themselves, they can be modeled directly on the output and compensated by using the information from the previous cycles. The reduced-order lifted model together with the residual term, hk, which comprises the unknown model error and the influences of the persisting as well as the random disturbances, follows directly from eqs 9 and 10

x˜ k+1(0) ) Φ ˜ x˜ k(0) + Γ h Uk

(17)

Yk ) Π ˜ x˜ k(0) + GUk + hk

(18)

where x˜ is the state vector of the reduced-order lifted model and Φ ˜,Γ ˜, Π ˜ , and G are the state-space matrices. The term hk is modeled with the stochastic difference equation

h k + wk h h k+1 ) h

(19)

˜ k + vk hk ) h

(20)

410


where

hk } [hkT(0) ‚‚‚ hkT(N - 1)]T ∈ R(N×ny)

(21)

In the above equation, h h k represents the model prediction error and the effect of periodwise-persisting disturbances that repeat themselves, whereas vk can be interpreted as the part that represents the effect of random disturbances. Both wk and vk are zero-mean white-noise variable sequences with respect to the period index. 3.1.3. Unlifting (Time-to-Time Transition Model). The resulting lifted formulation (eqs 17 and 18) can be used to develop a period-to-period control algorithm where the future optimal control sequence for a complete cycle is calculated at the beginning of each cycle. An obvious drawback of this approach is its periodwise compensation for the disturbances, which will mostly be effective for the persisting disturbances but not for the instantaneous ones. Therefore, we need to introduce a real-time feedback formulation, i.e., a time-to-time transition formulation that utilizes the measurements as they are available, that will overcome this limitation and address the disturbances as they take place. A brief explanation of the unlifting procedure to obtain a timeto-time transition model is given below; a more detailed description can be found elsewhere.15,23 Notice that the cycle-to-cycle description of the system (eqs 17 and 18) can be expressed in terms of incremental changes by differencing the model for two successive cycles as follows

˜ ∆x˜ k(0) + Γ ˜ ∆Uk ∆x˜ k+1(0) ) Φ

(22)

where

The transition from cycle to cycle completes the timevarying description of the system.

The periodic time-varying system of eqs 28, 29, and 31 can be written in the following compact form

zjk(n+1) ) A h zjk(n) + B h (n) ∆uk(n) h (n) zjk(n) + H(n)vk yk(n) ) C

n ) 0, ..., N - 1 (33)

zjk+1(0) ) Ψ h zjk(N) + [0 I]Twk Yk ) Y h k + vk

(24)

Here, ∆ indicates the backward difference with respect to the cycle index, i.e., ∆Uk ) Uk - Uk-1 and ∆uk(n) ) uk(n) - uk-1(n). Note that Y h k involves the repeating part of the disturbance effects and Yk incorporates the effect of random disturbances. Let us define

δk(n) } ∆x˜ k+1(0)

and

y j k(n) } Y hk with ∆uk(j) ) 0 for j g n (25)

This implies that the same input as for cycle k - 1 is implemented starting at time n of the cycle k and leads to

˜ ∆x˜ k(0) + δk(n) ) Φ ˜ n-1][∆ukT(0) ‚‚‚ ∆ukT(n-1)]T (26) [Γ ˜ 0 ‚‚‚ Γ h k-1 + Π ˜ ∆x˜ k(0) + y j k(n) ) Y [G0 ‚‚‚ Gn-1][∆ukT(0) ‚‚‚ ∆ukT(n-1)]T + wk-1 (27) where Γ ˜ n and Gn are the nth columns of the matrices corresponding to the input at time n, as shown in eqs 12 and 14. Note that δk(N) ) ∆x˜ k+1(0) and y j k(N) ) Y hk by definition (eq 25). The model equations (eqs 26 and 27) can be written for two consecutive time steps, i.e., for times n and n + 1, and taking the difference, we obtain

(32)

(34)

3.2. Kalman Filtering. The Kalman filter, which is based on the minimization of the variance of the estimation error, has been accepted as the standard linear optimal state estimator.19,24 A periodically timevarying Kalman filter is found to best suit our purposes. The one-step-ahead prediction equations for the timevarying system (eqs 32-34) are given by

h zjk(n|n-1) + B h (n) ∆uk(n) + zjk(n+1|n) ) A h (n) zjk(n|n-1)] Kk(n)[ykmeas(n) - C zjk+1(0|-1) ) Ψ h zjk(N|N - 1)

n ) 0, ..., N - 1 (35)

where zjk(n+1|n) denotes the prediction of zjk(n+1) based on measurements available up to time n. Kk(n) is the time-varying Kalman filter gain matrix. Let Rv and Rw be the covariance matrices of wk and vk, respectively; then, the filter gain is given by

Kk(n) ) h T(n)[H(n) RvHT(n) + C h (n) Pk(n) C h T(n)]-1 (36) Pk(n) C where Pk(n) is the covariance matrix of the estimate

Pk(n+1) ) Pk(n) - Pk(n) C h T(n)[H(n) RvHT(n) + h T(n)]-1C h (n) Pk(n) C h (n) Pk(n) C h Pk(N)Ψ h T + [0 I]TRw[0 I] Pk+1(0) )Ψ for n ) 0, ..., N - 1 (37)


Figure 7. Outlet purities for the controlled and uncontrolled SMB plant. A step change in the Henry’s constants take place at cycle 70, and the constants attain the values HA ) 4.4 and HB ) 2.4. The set purity is 99% for both outlets.

input moves, respectively.

Y nk p(n+2|n) )

[

∆Unk c(n+1) )

Figure 6. Operating parameter space in terms of the flow rate ratios (mj, j ) I, ..., IV).

Instead of using eqs 36 and 37 at each time step, one can, once and for all, iterate on the periodically timevarying Riccati difference equation (eq 37) until it converges to a periodic steady-state solution, i.e., Pk(n) f P∞(n) with n ) 0, ..., N - 1, and obtain the periodic steady-state gain matrices, i.e., K∞(0), ..., K∞(N - 1), according to the filter gain estimation equation (eq 36). A more detailed description of periodically time-varying Kalman filter equations can be found elsewhere.22,25 3.3. Prediction. The np-step-ahead prediction of the plant output based on the measurements up to time n and the nc-step future input moves can be obtained from the time-varying process model given by eqs 32-34

Y nk p(n+2|n) ) Szj(n+1) zjk(n+1|n) +

and

]

(39)

n+np

E QEk (i) cA,k (i) ∑ i)n+1

k ) 1, 2, 3, ...

zjk(n+1|n) is the state estimate obtained at time n from the observer (eq 35). Y nk p(n+2|n) and ∆Unk c(n+1) consist of the future system outputs and the future

[

∆uk(n+1) ·· · ∆uk(n+nc)

]

The matrices S zj(n+1) and Su(n+1) are constructed by successive substitution of model eqs 32-34 for future outputs over np steps as a function of nc future input movements. It is important to notice that the time index n is running as a subindex of the cycle index k. For instance, let j ) n + nc; if j g N, then j ) j - N and k ) k + 1. Also note that this interpretation keeps the sizes of the prediction and control horizons fixed. 3.4 Optimization of the Process Performance. The optimization of SMB processes has been addressed by several researchers, and it has been shown that the potential gain from the optimal operation of SMB units can be significant.26 Regardless of the complexity of the optimization schemes used, the SMB practitioner eventually faces the uncertainties in the system and the challenge of the robust operation of the plant at the calculated optimal operating conditions. Therefore, integration of on-line optimization and automatic control is required. For a given plant layout and switching period, t*, fixed a priori, the goal of optimal operation of an SMB could be to maximize the throughput and to minimize the solvent consumption, provided that the process and product specifications are fulfilled. These specifications can be formulated as constraints in the optimization problem. Because the final purity of the collected product is of main interest, purity constraints are defined on the average purities of the raffinate and extract outlets over the prediction horizon

Su(n+1) ∆Unk c(n+1) (38) for n ) 0, ..., N - 1

yk(n+2|n) ·· · yk(n+np+1|n)

Pave E (n+1) )

n+np

E E (i) + cB,k (i)] ∑ QEk (i)[cA,k

i)n+1

g Pmin (40) E

412


Figure 8. Controller action in terms of flow rate ratios, i.e., mj values. The mj values are calculated by using the average internal flow rates over a cycle. n+np

∑

Pave R (n+1)

R QRk (i) CB,k (i)

i)n+1

)

g Pmin (41) R

n+np

R R (i) + cB,k (i)] ∑ QRk (i)[cA,k

i)n+1

Here, QEk and QRk represent the extract and raffinate flow rates, respectively, and can be written in terms of internal flow rates, i.e., QE ) QI - QII and QR ) QIII QIV. Because these constraints are nonlinear, they are successively linearized at each time step around steadystate values and introduced into the optimization problem as linear inequality constraints. The output concentration profile and flow rate sequence of cycle k - 1 can be used to linearize the purity equations at cycle k, and moreover, the corresponding values can be updated at each time step by the new outlet measurements and the implemented flow rates. The linearized purity constraints subject to the system by eq 38 take the following form T

E nc Pave E (n+1) ) J1 (n+1) ∆Uk (n+1) + T

JE2 (n+1) zjk(n+1|n) + bE(n+1) - s1 g Pmin E

with s1 g 0

(42)

T

R nc Pave R (n+1) ) J1 (n+1) ∆Uk (n+1) + T

JR2 (n+1) zjk(n+1|n) + bR(n+1) - s2 with s2 g 0 g Pmin R

(43)

where JE1 , JE2 , JR1 , and JR2 are gradient vectors and bE and bR are scalars that result from the linearization. If

the purity constraints are infeasible, i.e., if there exists no input sequence that satisfies the purity constraints, then this would clearly cause problems for on-line applications. Therefore, the purity constraints are defined as soft constraints by introducing the slack variables s1 and s2 that guarantee feasibility for sufficiently large s1 and s2 values. The slack variables are included in the cost function and highly penalized to avoid excessive softening of the constraints. The predicted output concentration values have to be nonnegative to calculate the product purities given by eqs 42 and 43 correctly. Lower bounds on the predicted outputs can be introduced with the following inequality

Y nk p(n+2|n) g Y low k (n+2) - s3 where

Y

low k (n+2)

[

and

ylow k (n+2) ·· ) · low yk (n+np+1)

s3 g 0 (44)

]

(45)

is the vector of the output constraint trajectory over the prediction horizon. Certain constraints on the internal as well as on the external flow rates have to be included not only to be consistent with the physics of the process, but also to ensure a smooth operation of the plant with the available hardware. The implemented constraints on the internal and external flow rates are briefly mentioned here; their detailed description is provided in the Appendix. (i) The external flow rates must be nonnegative. Note that constraints on the external flow rates can be easily transformed into constraints on the internal flow rates (refer to section A.1 of the Appendix).


(ii) There is a maximum allowable pressure drop in the chromatographic columns. Because sections I and III are subjected to the highest flow rates in the unit, their internal flow rates must be limited. Moreover, nonnegativity of the internal flow rates has to be imposed as well (refer to section A.2 of the Appendix). (iii) Other constraints on the plant inputs can be defined to ensure a smooth and safe operation. For instance, a maximum allowable change in internal flow rates can be imposed to avoid sudden and significant pressure changes in the columns (refer to section A.3 of the Appendix). (iv) Ideally, SMB units have periodic dynamics over successive switches. Even though modeling of the process periodicity over successive switches would lead to a smaller model with fewer states and inputs, and therefore a smaller problem to be solved on-line, defining the period of the process as a complete cycle allows for the correction of differences in the individual column characteristics. Nevertheless, the controller has been restricted in such a way that the calculated control input sequence has identical values for successive switches. This allows the number of degrees of freedom to be decreased and leads to a smoother operation (refer to section A.4 of the Appendix). The defined inequality constraints on the outputs and inputs can be combined into the following linear inequality

C uk (n+1) ∆Unk c(n+1) g C k(n+1)

(46)

The defined equality constraints can be grouped in a similar way (refer to section A.5 of the Appendix) nc eq C u,eq k (n+1) ∆Uk (n+1) ) C k (n+1)

(47)

The performance objective for the SMB unit can be defined as maximization of the feed flow rate and minimization of the solvent consumption over the control horizon n+nc

min ∆Unk c(n+1),sk

[

[λ1(i) QD(i) - λ2(i) QF(i)] + λ3(n+1)s] ∑ i)n+1

Figure 9. Controller action represented in operating parameters spaces, i.e., the (mII, mIII) and (mIV, mI) planes. p1 and p3 are steady-state operating points for the nominal plant. p2 and p4 are steady-state operating points after the disturbance. The mj values are calculated by using the average internal flow rates over a cycle.

(48)

Here, λ1, λ2, and λ3 are possibly time-varying weights. The elements of the vector λ3 are chosen to be large relative to the other weights so that the slack variables are kept as small as possible. The minimization (eq 48), together with the inequality and equality constraints on the input and output values, given by eqs 46 and 47, respectively, constitutes a standard linear programming (LP) problem to be solved on-line. 4. Assessment of the Controller Performance on a Virtual Plant Because the flow rate ratios mj, given by eq 7, are widely used to characterize SMB performance, they will be used in the following discussion instead of the internal flow rates to present our results. Unless indicated otherwise, all controller actions are given in terms of mj values calculated as averages over one cycle. In the frame of equilibrium theory, which is based on the equivalence of the TMB and SMB configurations under the assumptions that the axial dispersion and the mass-transfer resistance effects are negligible, the suf-

Figure 10. Outlet purities for the controlled SMB plant.The isotherm parameters of the actual plant are HA ) 3.35 and HB ) 2.35, whereas the model parameters are HA ) 4 and HB ) 2. The set purity is 99% for both outlets.

ficient and necessary conditions to achieve complete separation and complete regeneration are given by the following inequalities

414


Figure 11. Controller action in terms of flow rate ratios, i.e., mj values. The mj values are calculated by using the average internal flow rates over a cycle.

HA < mI < ∞

(49)

HB < mII < mIII < HA

(50)

0 < mIV < HB

(51)

for nonporous particles as in this case.6 If the flow rate ratios in sections I and IV fulfill the necessary conditions (eqs 49 and 51), the operating parameters mII and mIII define a plane that is divided into several regions characterizing the separation performance (Figure 6). First, the separation can take place if and only if mII < HA and mIII > HB, because violation of the former or the latter constraint leads to flooding of the extract or raffinate outlet with the solvent, respectively. The triangular area defined by inequality 50 is the region where the products at both the raffinate and the extract outlets are pure. For the operating points constituting the region defined by mII < HB and HB < mIII < HA, the light component B is carried to the extract outlet, and the purity for the extract is low, whereas the purity of raffinate is not affected. On the other hand, the region where mIII > HA and HB < mII < HA defines the operating points where the raffinate outlet is polluted, and the extract purity is not affected. All product streams have poor purity for the operating points corresponding to the region mIII > HA and mII < HB. The vertex of the complete separation region W1 represents the optimal operating conditions in terms of solvent consumption and productivity per unit mass of stationary phase. As mentioned, the separation performance is characterized by the position of the operating point in the (mII, mIII) plane, as well as in the (mIV, mI) plane (Figure 6). The violation of the constraint on the operating variable mI, i.e., HA > mI, leads to incomplete regeneration of

the solid phase and results in poor raffinate purity. On the other hand, violation of the constraint on the operating variable mIV, i.e., mIV > HB, leads to pollution of the extract stream as a result of incomplete regeneration of the mobile phase. The minimum possible solvent consumption corresponds to the vertex of the rectangular area, i.e., W2. The equilibrium theory is validated by comparison with experimental data in several previous works.6,27 In all simulations, the plant is operated for one complete cycle at a fixed operating point, namely, mI ) 4.2, mII ) 2.2, mIII ) 3.8, mIV ) 1.9, and the controller is switched on at the beginning of the second cycle. This is necessary because the control algorithm needs the input and output sequences of the first cycle for initialization. All of the simulations include zero-mean white noise added to the output measurements. A standard deviation of 2% of the measured concentration value is assumed. The lowest acceptable purity for the extract and the raffinate is defined as 99%, i.e., Pmin ) Pmin ) E R 99%. The examples below address different types of disturbances that can act on the system and affect the operation performance negatively. The adsorption isotherm parameters, i.e., Henry’s constants, are the most significant physical data of the system characterizing the separation of the species. Hence, the operation is highly sensitive to disturbances that cause changes in the isotherm parameters. For some of the examples, the evolution of the product purities for the open-loop operation of the SMB unit, i.e., uncontrolled operation, is also given to illustrate the effects of the disturbances. It is worth mentioning that open-loop simulations aim only at illustrating the significance of the studied disturbances; instead of choosing proper operating conditions for each case, all open-loop simulations are performed under the same operating conditions, namely


Figure 12. Controller action represented in operating parameters spaces, i.e., the (mII, mIII) and (mIV, mI) planes. p5 and p7 are the startup operating points. p6 and p8 are the steady-state operating points. The mj values are calculated by using the average internal flow rates over a cycle.

those corresponding to the reference operating paramref ref ref eters, i.e., mref I ) 4, mII ) 2.1, mIII ) 3.9, and mIV ) 2.1. 4.1. Step Disturbance. The first example shows the performance of the controller when the isotherm parameters undergo a step change, which might occur as the result of sudden temperature changes in the system. With reference to Figure 7, the plant is started with the operating parameters given above. Because the columns constituting the SMB unit are initially saturated with the mobile phase only, i.e., there are no solutes, it takes at least a couple of switches for the solutes to reach the outlet positions, i.e., A to reach the extract and B the raffinate. It is worth mentioning that it takes more time, depending on the operating conditions, for the impurities, i.e., B in the extract and A in the raffinate, to reach the outlet positions. This behavior can be observed in the purities of the uncontrolled plant. In this particular case, the chosen mII and mIII values are correct, whereas mI and mIV lead to poor regeneration in sections I and IV, respectively. As a result, a significant decrease in the extract and raffinate purities starts after three cycles and continues until the unit reaches the cyclic steady state.

The controller is switched on at cycle two, and the operating conditions are adapted by the controller such that the purity requirements are fulfilled and the operating performance is optimized. Note that, in the first part of the operation, i.e., until the disturbance takes place, the plant and model have the same physical parameters, i.e., HA ) 4, HB ) 2, and the only plantmodel mismatch is due to linearization and model reduction. The disturbance takes place at cycle 70 where the Henry’s constants attain the new values HA ) 4.4 and HB ) 2.4, i.e., corresponding to 10 and 20% changes in HA and HB, respectively. As a consequence, the selectivity decreases to S ) 1.83, i.e., 8% less than the nominal value, which implies that the separation task becomes more difficult. The significance of the illustrated disturbance can be observed from the evolution of the product purities for the uncontrolled operation of the plant. Such a disturbance leads to a dramatic reduction in the product purities for the uncontrolled plant (Figure 7). On the other hand, the controller adapts the operating conditions and fulfills the product specification within five and eight cycles for the extract and raffinate outlets, respectively. With reference to Figure 7, it is worth discussing the and Pmin difference between the set purities, i.e., Pmin E R (dotted line in the figure), and the values of the plant outputs under the control action (solid and dashed lines in the figure). The latter are calculated by numerical integration of the continuous plant outputs, thus accurately representing the average purities, i.e., the physically meaningful purities that have to fulfill process specifications. However, these values are not accessible to the controller, because the controller is based on information about the plant output sampled only at discrete times. In fact, the average purities used by the optimizer are obtained through the formulas given in eqs 40 and 41, hence representing a rather rough approximation of the real values. Note that this is the case also in the absence of offset between the instantaneous purity values calculated by the model and those measured at the plant. The magnitude, as well as the sign, of the differences observed in the Figure 7 depends on the shape of the composition profiles of the product outputs and decreases as the profiles get smoother. It is worth mentioning that, in practice, the output concentration profiles are generally smoother and more symmetric because of the unavoidable dead volume in the line from the column outlets to the detectors.28 Therefore, one might expect that the difference between the set average purity and the true purity values will be insignificant for practical applications. The controller action, in terms of the mj values, for the whole operation is given in Figure 8. Let us interpret the controller action on the (mII, mIII) and (mIV, mI) planes. With reference to Figure 9, the solid lines indicate the complete separation and regeneration regions for the model and for the nominal plant, whereas the areas outlined by dashed lines represent the regions for the actual plant with the new isotherm parameters attained after the disturbance. The points p1 and p3 indicate the steady-state operating conditions attained during the first part of the operation (nominal plant). It is observed that the controller adapts the operating conditions in such a way that the plant is operated very close to the vertex of the triangular area

416


Figure 13. Controller action for a complete cycle in terms of external flow rates, i.e., desorbent, extract, feed, and raffinate flow rates, at steady-state operating conditions. The period between two successive switches, i.e., t*, is 480 s, and a complete cycle comprises eight switching periods.

in the (mII, mIII) plane, i.e., the optimal operating point according to the equilibrium theory, and close to the boundary of the complete regeneration region in the (mIV, mI) plane. After the disturbance, the complete separation and regeneration regions for the plant are significantly different. The controller reacts to the changes and drives the operating points into the new separation and regeneration regions in order to fulfill the product specifications as soon as possible. After satisfying the purity requirements, the operating conditions are adjusted to maximize throughput (the operating point is moved to the vertex of the new complete separation region) and to minimize solvent consumption. The points p2 and p4 are the new steady-state operating conditions attained after the disturbance. 4.2. Plant-Model Mismatch. Parameter uncertainty concerning the adsorption isotherm of the species to be separated or the column characteristics, i.e., column efficiency, length, or porosity, is typically faced in practice. These parameters are difficult to determine and can change during the operation of the unit. The particular example presented here addresses significant uncertainty in isotherm parameters, whereas model bias due to differences in column characteristics is discussed elsewhere.16 Let us consider a plant with the isotherm parameters HA ) 3.35 and HB ) 2.35, i.e., -16 and 17% different from the corresponding model values, respectively. Note that these Henry’s constants lead to a selectivity of S ) 1.43, which implies that the actual separation task is much more difficult than the nominal case. Same startup operating conditions as in the previous example are used for the first cycle. The chosen startup operating parameters correspond to the “no pure outlet” region in the (mII, mIII) plane (Figure 6) where both the raffinate and extract purities are low.

The evolution of the product purities for the controlled plant is illustrated in Figure 10. The controller is able to fulfill the product specifications within eight cycles despite the significant parameter uncertainty. The controller action is given in Figures 11 and 12. Referring to Figure 12, the complete separation and regeneration regions for the model are given by solid lines, whereas the regions indicated by dashed lines correspond to the actual plant. One can see that the complete separation region applicable to the plant is smaller than that corresponding to the model, i.e., the plant is operating with a lower selectivity. On the other hand, the complete regeneration area is such that the separation task requires less-than-nominal fresh solvent. The points indicated as p5 and p7 show the startup operating points. It can be observed that the controller drives the operating points into the correct separation region as soon as the purities decrease below their set values and then operates the unit close to the vertex of the triangle in the (mII, mIII) plane. At the same time, it reduces the solvent consumption by adapting the operating conditions in sections I and IV, i.e., the operating point is driven close to the borders of the complete regeneration area in the (mIV, mI) plane. The points indicated as p6 and p8 correspond to the steady-state operating points. The controller action in terms of external flow rates for a complete cycle, comprising eight switches, at steady-state operating conditions is illustrated in Figure 13. The inlet flows (QD, QF) and outlet flows (QE, QR) are given by dashed and solid lines, respectively. In our case, the time period t* between two successive switches is 480 s and is divided into eight segments, each of which comprise 64 sample points for a complete cycle, i.e., N ) 64. It can be seen from Figure 13 that the controller is implementing similar flow rate sequences for successive switches.


Figure 15. Outlet purities for the controlled and uncontrolled SMB plant. A ramp change in the Henry’s constants take place between cycle 70 and 110. The final values of the Henry’s constants are HA ) 4.4 and HB ) 2.3. The set purity is 99% for both outlets.

Figure 14. Controller action for a switch period in terms of flow rate ratios, i.e., mj values, at steady-state operating conditions. The moving directions of the operating points are indicated by arrows.

Figure 14 shows the steady-state flow rate sequence implemented between two successive port switches in the (mII, mIII) and (mIV, mI) planes. With reference to the (mII, mIII) operating plane, the first operating point of the switch period is located within the “complete separation” region. Then, the controller moves the operating point into the no pure outlet region to increase the throughput and drives it back into the complete separation region at the end of the switch period (the direction of the movement is indicated by the arrow). Note that the controller is operating the plant such that the throughput is 19% higher in the second half of the switch period than it is in the first half. It has been shown recently29 that such operating modes can improve the productivity of SMB units. The average operating point corresponds to the steady-state operating point indicated as p6 in Figure 12. One can see that there are six different operating points shown in the figure, because two of them are repeated within one switch period. The operating conditions in the (mIV, mI) plane are also changing in time. The controller is increasing the amount of fresh desorbent fed into section I during the first half of the switch period and decreasing it during the second half (the direction of the movement is indicated by the arrow). The average operating point

corresponds to the steady-state operating point indicated as p8 in Figure 12. 4.3. Gradual Changes in the System Characteristics. Common in SMB operation are slow changes in the characteristics of the system because of chemical and mechanical aging of the solid phase, leading to gradual changes in the retention behavior of the species to be separated. This situation can be simulated by implementing a ramp change in the Henry’s constants. Figure 15 illustrates a situation in which the isotherm parameters are changing from their initial values, HA ) 4 and HB ) 2, to their final values, HA ) 4.4 and HB ) 2.3, over a time interval of 40 cycles. The first part of the operation corresponds to the nominal case in which the plant and the model have identical isotherm parameters. The change in the Henry’s constants starts at cycle 71 and continues over 40 cycles, i.e., until the end of cycle 110. This results in a continuous decrease in the selectivity, i.e., from S ) 2 to S ) 1.91. It can be seen that the product purities of the uncontrolled plant exhibit a gradual decrease over the period where the isotherm parameters undergo the ramp change. On the other hand, the controller handles the changes in the system characteristics and fulfills the product specifications throughout the operation. Figure 16 illustrates the gradual shift in the complete separation and regeneration regions because of the changes in the retention behavior of the species. The areas indicated by solid lines correspond to the nominal plant, which is applicable until the isotherm parameters start to change, i.e., until the end of cycle 70. The regions start to shift because of the changes in the Henry’s constants and reach their final position at cycle 110 (depicted by dotted lines). The areas given by dashed lines correspond to the separation and regeneration regions at cycle 95. The controller is adapting the operating conditions to the changes in the characteristics of the system and moves the operating points so as to fulfill the necessary conditions for high-purity separation. The symbols indicate the operating points for cycle 70, i.e., the steady-state operating conditions for the nominal plant, and for cycles 95, 110, and 130, i.e., the new steady-state operating conditions after the Henry’s constants attain their final values. This example shows the advantage of a real-time feedback formulation that utilizes the measurements as they become available.

418


Figure 17. Number of cycles that the plant produces offspecification products versus the duration of the ramp change in the isotherm parameters. The Henry’s constants HA and HB undergo changes of +10 and +15%, respectively, i.e., their final values are HA ) 4.4 and HB ) 2.3.

Figure 16. Trajectory of the operating point and the corresponding complete separation and regeneration regions. O, Steady-state operating conditions for the nominal plant; 4, operating points at cycle 95; b, operating points at cycle 110; 0, steady-state operating point after the Henry’s constants attain their final values. The mj values are calculated by using the average internal flow rates over a cycle.

Figures 7 and 15 illustrate two extreme cases, one in which the controller reacts to an instantaneous step disturbance and the products are off-specification during a few cycles and one in which the disturbance is gradual and the controller is able to avoid any off-specification production (Figure 15). In Figure 17, the transition between these two types of behavior is investigated. The number of cycles during which either the extract or the raffinate product is off-specification, as a function of the duration of the ramp disturbance (with the overall change of Henry’s constants always the same, i.e., HA from 4 to 4.4, and HB from 2 to 2.3), is shown. For instance, zero duration corresponds to a step change, and a ramp duration of 40 cycles corresponds to the case reported in Figures 15 and 16. It can be seen that, for a step change in the isotherm parameters, which might occur as a result of sudden temperature changes, the controller can fulfill the product specifications within six and eight cycles for the extract and the raffinate outlets, respectively. If the changes take place over a longer time, the number of off-specification cycles first

increases and then decreases, reaching a maximum of about 15 cycles for both extract and raffinate, but at different ramp durations. If the system characteristics undergo a very slow change, which might occur as a result of slow aging of the solid phase, then the controller can adjust the operation so that the product specifications are always fulfilled. The presence of a maximum in the diagrams of Figure 17 indicates a competition between two opposite trends. On one hand, increasing the time period over which the same disturbance occurs makes such a disturbance milder and easier to handle by the controller. On the other hand, a sharper disturbance, corresponding to a short ramp duration, yields larger deviations from the set product specifications and activates a faster response from the controller. No special effort has been made to tune the controller. One can optimize the behavior of the controller to obtain that most appropriate for the disturbance dynamics typical of the system under consideration. The performance of the controller is encouraging enough to extend its application for separation tasks characterized by nonlinear adsorption isotherms. It is also important to devise a strategy to enforce constraints on the way the flow rates are changed, so as to allow for smoother variations if necessary. Another issue to be addressed in future work is the optimal control of SMB units with fewer columns, which are more sensitive to disturbances but also more cost-efficient given that the packing material, particularly for chiral separations, constitutes a significant part of the fixed cost for SMB units. In addition, experimental verification of the controller performance is in progress. Acknowledgment The authors are grateful to Dr. Giancarlo FerrariTrecate, Dr. Francesco Borrelli, and Marc Lawrence for helpful discussions. The support of ETH Zurich through Grant TH-23′/00-1 is gratefully acknowledged. Notation A, B, C ) state-space matrices A h, B h, C h , H, Ψ h ) state-space matrices of the time-to-time transition model Acr ) column cross section, cm2 c ) concentration, g/L

Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 419 D ) apparent axial dispersion coefficient, cm2/s δQmax ) maximum allowable flow rate change ∆Unc ) vector of nc-step future inputs ) bed void fraction Φ, Γ, Π, G ) state-space matrices of the lifted model Φ ˜, Γ ˜, Π ˜ , G ) state-space matrices of reduced-order lifted model h ) residual term in the model equations HA ) Henry’s constant of component A HB ) Henry’s constant of component B K ) Kalman filter gain matrix k ) cycle index L ) column length, cm λ ) weighting factor in the cost function m ) flow rate ratio N ) number of time steps per cycle n ) time index within a cycle nc ) control horizon ncl j ) number of columns in section j np ) prediction horizon Ntst ) number of time steps during a switching period P ) covariance matrix of the estimate PE ) purity of extract outlet PR ) purity of raffinate outlet Q ) volumetric fluid flow rate, mL/s Q ) vector comprising internal flow rates q* ) adsorbed-phase concentration, g/L Rv ) covariance matrix of the noise sequence v Rw ) covariance matrix of the noise sequence w s ) slack variable t ) time, s t* ) switch time, s u ) vector of manipulated variables, mL/s U ) vector grouping the input values for one cycle v ) superficial velocity, cm/s V ) volume of one column, mL v, w ) zero-mean white-noise sequences x ) state vector x˜ ) state vector of the reduced-order model y ) vector of output concentrations Y ) vector grouping the output values for one cycle Y low ) vector of the output constraint trajectory over the prediction horizon Y np ) vector of np-step future outputs z ) axial coordinate, cm zj ) state vector of the time-to-time transition model zj(n+1|n) ) one-step-ahead prediction of the states based on the measurements at time n Subscripts and Superscripts ave ) average D ) desorbent E ) extract F ) feed g ) space index (grid point number) h ) column position index, h ) 1, ..., 8 i ) component index, i ) A, B I/O ) inlet/outlet stream in ) column inlet j ) section index, j ) I, ..., IV max ) maximum meas ) measurement min ) minimum out ) column outlet R ) raffinate ref ) reference value used for linearization

Appendix A.1. Constraints on External Flow Rates. Nonnegativity of external flow rates (QE, QR, QF) can be imposed through internal flow rates.

QE ) QI - QII g 0

(A.1)

QR ) QIII - QIV g 0

(A.2)

QF ) QIII - QII g 0

(A.3)

Let V be the (mnc) × (nnc) block diagonal matrix of the form

[ ]

v 0 V) · ·· 0

0 v ·· · 0

··· ··· ··· ···

0 0 ·· · v

(A.4)

where m and n are the numbers of rows and columns of the matrix v, respectively. nc is the length of the control horizon, and the matrix v is defined as

[

1 -1 0 0 v ) 0 0 1 -1 0 -1 1 0

]

(A.5)

The nonnegativity of the feed, extract, and raffinate flow rates for the whole control horizon can be enforced by nc V∆Unk c(n+1) g -V[Qnrefc + Uk-1 (n+1)]

(A.6)

Qnrefc is a vector consisting of the reference internal flow nc (n+1) rates grouped for the whole control horizon. Uk-1 is the input sequence of the previous cycle defined in the same way as Unk c(n+1) given by eq 39. A.2. Lower and Upper Bounds on the Internal Flow Rates. Lower and upper bounds on the control inputs, i.e., on the internal flow rates, can be defined with the following inequalities nc c + Uk-1 (n+1) -∆Unk c(n+1) g -Unmax

(A.7)

nc c - Uk-1 (n+1) ∆Unk c(n+1) g Unmin

(A.8)

where the upper and lower bounds are defined in terms of deviation variables c c ) Qnmax - Qnrefc Unmax

(A.9)

c c ) Qnmin - Qnrefc Unmin

(A.10)

nc c and Qnmin are the vectors consisting of the maxiQmax mum and minimum allowable internal flow rates grouped for the whole control horizon, respectively. A.3. Constraints on the Maximum Internal Flow Rate Changes. Let δ indicate the difference with respect to the time index, i.e., δQk(n+1) ) Qk(n+1) Qk(n). The rate of input change can be bounded as

|δQk(n+1)| e δQmax

(A.11)

Implementation for the whole control horizon follows nc c -J3∆Unk c(n+1) g -δQnmax + J3Uk-1 (n+1) - J4uk(n) (A.12)

420


nc c J3∆Unk c(n+1) g -δQnmax - J3Uk-1 (n+1) + J4uk(n) (A.13)

where

[

I -I J3 ) 0 ·· · 0

0 I -I ··· 0

··· ··· ··· ··· ···

0 0 ·· · ··· -I

0 0 0 ·· · I

]

and

[]

I 0 J4 ) · ·· 0

[

where

I s -I s 0 s 0 -I s J5 ) I· ·· ··· ·· · s 0 0 I

··· ··· ···

0 0 l · · · -I s

]

C u,eq k (n+1) ) J5

where

C uk (n+1) )

[] I -I V -J3 J3

[

and

Literature Cited

(A.15)

(A.16)

(A.17)

(A.18)

u

and

(A.20)

nc C eq k (n+1) ) -J5Uk-1(n+1) (A.21)

(A.14)

Note that the identity matrix denoted by I s has a dimension of (nuNtst) × (nuNtst), where Ntst is the number of sampling times between two successive switches. A.5. Grouping Inequality and Equality Constraints. The defined inequality constraints can be put into the form

C uk (n+1) ∆Unk c(n+1) g C k(n+1)

nc eq C u,eq k (n+1) ∆Uk (n+1) ) C k (n+1)

where

nc is the vector consisting of the maximum allowδQmax able flow rate change for the whole control horizon, and uk(n) is the control input implemented at time n. A.4. Identical Control Input Sequence for Successive Switches. Identical control inputs for successive switches can be implemented in different ways, one of which is given here

nc (n+1) J5∆Unk c(n+1) ) -J5Uk-1

The defined equality constraints can be written as

S (n+1) JE1 (n+1) JR1 (n+1)

C k(n+1) ) nc c Unmin - Uk-1 (n+1) nc nc -Umax + Uk-1(n+1) nc (n+1)) -V(Qnrefc + Uk-1 nc c + J3Uk-1 (n+1) - J4uk(n) -δQnmax nc c - J3Uk-1 (n+1) + J4uk(n) -δQnmax zj low Y k (n+2) - S (n+1) zjk(n+1|n) - s3 Pmin - JE2 (n+1)zjk(n+1|n) - bE(n+1) - s1 E Pmin - JR2 (n+1) zjk(n+1|n) - bR(n+1) - s2 R

]

(A.19)

(1) Juza, M.; Mazzotti, M.; Morbidelli, M. Simulated movingbed chromatography and its application to chirotechnology. Trends Biotechnol. 2000, 18, 108. (2) Nicoud, R. M. Simulated Moving Bed: Some possible applications of biotechnology. In Bioseparation and Bioprocessing; Subramanian, G., Ed.; Vol. I, Part I; Wiley-VCH: New York, 1998. (3) Francotte, E.; Richert, J. Application of simulated moving bed chromatography to the separation of the enantiomers of chiral drugs. J. Chromatogr. A 1997, 769, 101. (4) Imamoglu, S. Simulated Moving Bed Chromatography (SMB) for Application in Bioseparation. Adv. Biochem. Eng. Biotechnol. 2002, 76, 211. (5) 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. (6) Mazzotti, M.; Storti, G.; Morbidelli, M. Optimal operation of simulated moving bed units for nonlinear chromatographic separations. J. Chromatogr. A 1997, 769, 3. (7) Migliorini, C.; Gentilini, A.; Mazzotti, M.; Morbidelli, M. Design of Simulated Moving Bed Units under Nonideal Conditions. Ind. Eng. Chem. Res. 1999, 38, 2400. (8) Biressi, G.; Ludemann-Hombourger, O.; Mazzotti, M.; Nicoud, R. M.; Morbidelli, M. Design and optimisation of a simulated moving bed unit: Role of deviations from equilibrium theory. J. Chromatogr. A 2000, 876, 3. (9) Kloppenburg, E.; Gilles, E. D. Automatic control of the simulated moving bed process for C8 aromatics separation using asymptotically exact input/output linearization. J. Process Control 1999, 9, 41. (10) Schramm, H.; Gru¨ner, S.; Kienle, A.; Gilles, E. Control of moving bed chromatographic processes. In Proceedings of European Control Conference 2001; Martins de Carvalho, J. L., Fontes, F. A. C. C., de Pinho, M. d. R., Eds.; FEUP Edic¸ o˜es: Porto, Portugal, 2001; p 2528. (11) 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. (12) Bemporad, A.; Morari, M. Model predictive control: A survey. Lect. Notes Control Inf. Sci. 1999, 245, 207. (13) Morari, M.; Lee, J. Model predictive contol: Past, present and future. Comput. Chem. Eng. 1999, 23 (4-5), 667. (14) Natarajan, S.; Lee, J. Repetitive model predictive control applied to a simulated moving bed chromatography system. Comput. Chem. Eng. 2000, 24, 1127. (15) Lee, J.; Natarajan, S.; Lee, K. A model-based predictive control approach to repetitive control of continuous processes with periodic operations. J. Process Control 2001, 11, 195. (16) Abel, S.; Erdem, G.; Mazzotti, M.; Morari, M.; Morbidelli, M. Optimizing Control of Simulated Moving BedssLinear Isotherm. J. Chromatogr. A, manuscript submitted. (17) Abel, S.; Mazzotti, M.; Erdem, G.; Morari, M.; Morbidelli, M. Optimization Based Adaptive Control of Simulated Moving Beds. In Proceedings of the 3rd Pasific Basin Conference on Adsorption Science and Technology Kyongju, Korea; Lee, C. H., Ed.; World Scientific Pub. Co. Inc., ISBN: 9812383492: 2003. (18) Guiochon, G.; Golshan-Shirazi, S.; Katti, A. M. Fundamentals of Preparative and Nonlinear Chromatography; Academic Press: Boston, 1994. (19) Kalman, R. E. A new approach to linear filtering and prediction problems. Trans. ASME, J. Basic Eng. 1960, 82, 35. (20) MATLAB, Control System Toolbox; The MathWorks, Inc.: Natick, MA, 2001. (21) Pernebo, L.; Silverman, M. Model reduction via balanced state space representations. IEEE Trans. Autom. Control 1982, 27, 382.

Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 421 (22) Morari, M.; Lee, J.; Garcia, C. Model Predictive Control, Preprint, to be published by Prentice-Hall in 2003. (23) Lee, K.; Chin, I. S.; Lee, H.; Lee, J. Model predictive control technique combined with iterative learning for batch processes. AIChE J. 1999, 45 (10), 2175. (24) Grewal, M.; Andrews, A. Kalman Filtering: Theory and Practice, Using MATLAB, 2nd ed.; John Wiley and Sons: New York, 2001. (25) Lee, K. S.; Lee, J. H. Implementations of Iterative Learning Control Methodology. In Iterative Learning Control: Analysis, Design, Integration and Applications; Bien, Z., Xu, J., Eds.; Kluwer Academic Publishers: Boston, 1998. (26) Du¨nnebier, G.; Fricke, J.; Klatt, K.-U. Optimal design and operation of simulated moving bed chromatographic reactors. Ind. Eng. Chem. Res. 2000, 39, 2290.

(27) Pedeferri, M.; Zenoni, G.; Mazzotti, M.; Morbidelli, M. Experimental analysis of a chiral separation through simulated moving bed chromatography. Chem. Eng. Sci. 1999, 54, 3735. (28) Zenoni, G.; Pedeferri, M.; Mazzotti, M.; Morbidelli, M. Online monitoring of enantiomer concentration in chiral simulated moving bed chromatography. J. Chromatogr. A 2000, 888, 73. (29) Zhang, Z.; Mazzotti, M.; Morbidelli, M. PowerFeed: An innovation to the simulated moving bed technology. J. Chromatogr. A, in press.

Received for review May 5, 2003 Revised manuscript received September 12, 2003 Accepted September 24, 2003 IE030377O

Automatic Control of Simulated Moving Beds - ACS Publications

Recommend Documents