Article pubs.acs.org/IECR
Multibatch Model Predictive Control for Repetitive Batch Operation with Input−Output Linearization Junghui Chen* and Yu-Hung Lin Department of Chemical Engineering, Chung-Yuan Christian University, Chung-Li, Taiwan 320, Republic of China ABSTRACT: Because of the repetitive nature of the batch process, in conventional batch process control, the start of the subsequent batch is often assumed to be the same as that of the prior batch. Nevertheless, in real operation, under the influence of various exterior factors, most batch processes involve gradual changes over batches and, thus, the end-point condition would affect the start of the subsequent process batch. The initial condition at each batch may not be reset to be the desired initial condition. In order to ensure an optimum operation, it is essential to taking into account the effect of the previous batch. In this paper, a multibatch dynamics model is derived. It takes into account this batchwise dynamics. The derived process model can provide the benefits of better prediction and performance. In order to overcome the problem of nonlinearity in the batch process, the nonlinear process is transformed to a linear process through input−output feedback linearization. The identification method of the multibatch model is proposed for such batch processes, and the model will be able to reflect the particular type of batch process, which can be transformed to a linear input−output relation. Based on this model, a model predictive control scheme is derived for the case of unconstrained control as well as constrained control. The applications are discussed through a fed-batch reactor problem for the biosynthesis of penicillin to demonstrate the advantages of the proposed method, in comparison with the conventional methods.
1. INTRODUCTION Batch processing has its past applications in fine and specialty chemicals, but has newly found applications in emerging industries, like biotechnology and material science. Batch processes have historically lagged continuous processes in terms of development and deployment of advanced optimization and control tools. They also present unique challenges like nonstationary operating points and variations in the initial charge conditions. Using tools widely adopted for continuous industrial process control to solve the batch control problem is not easy, because the conventional controller design for the continuous process may not be able to achieve a specified profile response and process raw materials into products in finite duration. From the viewpoint of control system design for batch processes, the control problem frequently follows a time-varying reference profile over a fixed time interval. Amann et al. proposed quadratic criterion-based iterative learning control (Q-ILC), which made use of a quadratic optimization criterion based on a model.1 The input trajectory was successively refined using previous batch data, and it could be formulated under different settings. As Q-ILC lacked real-time correction capability, Lee et al. proposed BatchMPC (BMPC).2 Their use of model predictive control (MPC) in the design concept provided effective tracking. It was also able to handle constraints. Based on a time-varying linear system model, the technique utilized not only the incoming measurements from the ongoing batch but also the information from the past batches. Chin et al. then devised a two-stage algorithm by modifying and combining existing QILC and BMPC techniques to handle the effects of real-time disturbance on the performance of ILC.3 In the control design for the nonlinear processes, while the control technique can be applied directly to nonlinear processes, for a certain class of processes that can be expressed in © 2012 American Chemical Society
the canonical form, it is possible to convert the process through a linearization transform and result in a process with linear characteristics. Then, linear control can be directly applied. Brockett introduced the feedback linearization technique to transform a class of nonlinear processes to give a linear relation between the input and the output.4 Their method was a transformation of the nonlinear system into an equivalent linear system and the conventional linear control strategy can be easily applied. Kravaris and Chung have applied the input− output linearization to temperature control of a batch process.5 They were able to successfully control the nonlinear process using a linear control method. Because of the repetitive nature of the batch process, the endpoint condition in real operation would affect the start of the subsequent process batch. Under the influence of various exterior factors, most batch processes involve gradual changes over batches. For example, in penicillin fermentation production, the initial culture volume tends to gradually increase over batches, because the remainder of the culture from the previous batch is added to the coming batch.6 In the polymerization process, the reactor surfaces are extensively fouled with polymer deposits, so it is extremely difficult to clean the polymer deposits, which can affect the initial conditions.7 Thus, multibatch (MB) dynamics can be presented in a batch process. For model-based control, these relations should be taken into account. The derived process model could provide the benefits of better prediction and performance. In this paper, a multibatch model predictive control (MB-MPC) strategy will be constructed. Through input−output feedback Received: Revised: Accepted: Published: 9598
September 4, 2011 March 1, 2012 June 10, 2012 June 11, 2012 dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
where LgLr−1 f h(x) ≠ 0 and r is the order of the differentiation. The input−output system is linearized. The transformed input vi is mapped to the original input u via
linearization, the nonlinear process is transformed to a linear process. The identification method of an MB model is proposed for batch processes whose previous batches have an effect on the subsequent batches. The MB model is able to reflect the particular type of batch process whose dynamics is inevitable. Based on this MB model, a model predictive control (MB-MPC) method is derived for the cases of unconstrained control as well as constrained control. The rest of the paper is organized as follows. Section 2 gives a brief introduction to input−output feedback linearization and formulates the multibatch model. MB-MPC is derived in Section 3. Section 4 then details a case study that simulates a process to assess its performance and applicability before some concluding remarks.
ui(t ) =
y(t ) = h(x(t ))
Lg Lrf − 1h(x i(t ))
= ψ (r)(x i(t ), vi(t )) (5)
Now the input vi and the output yi have a linear relationship:
yi(r)(t ) = vi(t )
(6)
Figure 1 illustrates this transformation. vi and yi have a linear relationship, because the nonlinearity of the process between ui
2. MULTIBATCH MODEL DERIVATION USING NONLINEAR FEEDBACK LINEARIZATION The nonlinear nature means that linear control techniques cannot be directly applied. For a system of companion form, the process can be linearized. Based on the linearized system, because of the cleaning or remainder from the previous batch, a model that takes into account the batchwise dynamics resulting from a gradual change over batches is developed. Before the multibatch dynamic model is proposed, a brief overview of the nonlinear state feedback linearization related to the ground for the proposed method is presented first. 2.1. Overview of Nonlinear State Feedback Linearization. Global linearization involves finding global nonlinear transformation on the states and/or the manipulated inputs so that the transformed system possesses certain linearity characteristics. The feedback linearization approach is a transformation of the nonlinear system into an equivalent linear system through a change of variables and a suitable control input. It can be applied to a system of so-called companion form or controllability canonical form. The input−output feedback linearization is applied to a batch system of x ̇(t ) = f (x(t )) + g (x(t ))u(t )
vi(t ) − Lrf h(x i(t ))
Figure 1. Feedback linearization of a nonlinear batch system.
and yi has been “removed”. Note that if g is the function of the state and the input, one must determine a proper state and input transformation so that the nonlinear system dynamics is transformed to a linear input−output relationship and then a controller based on linear control is formulated. 2.2. Multibatch Dynamic Model. Now the nonlinear dynamics (eq 1) is algebraically transformed into linear one (eq 6). The discrete state space model derived from the input− output relationship (eq 6) then can be described by zi(k) = Azi(k − 1) + bvi(k − 1) yi (k) = q Tzi(k) + d
(7)
where A, b, and q are r × r, r × 1, and 1 × r paramedic matrices, respectively, and d is the model bias. Note that eq 7 is for one batch run only. z is the transformed state vector z∈Rr. By repetitive substitution, the process can be represented by
t = [0 K ] (1)
where x is the state vector (x ∈ Rn), y is the output (y ∈ R) and u is the input (u ∈ R). f(x), g(x), and h(x) are the nonlinear functions of the state. f(x) and g(x) are the vector fields on Rn and h(x) is a scale field on Rn. Almost all the SISO chemical processes can be described by a model form of eq 1.5,8 The righthand side of eq 1 is a linear function of the manipulated input. The process is operated at the fixed time of K and is repeated exactly after each batch. By continuously differentiating the output of eq 1,
yi = s + H 0pvi + d
(8)
where T
s = ⎡⎣q Tzi(0) q TAzi(0) ··· q TAK zi(0)⎤⎦ d = [ d d ··· d ]T T
yi (̇ t ) = Lf h(x i(t )) + Lg h(x i(t ))ui(t )
vi = ⎡⎣vi(0) vi(1) ··· vi(K − 1)⎤⎦
(2)
T yi = ⎡⎣ yi (1) yi (2) ··· yi (K )⎤⎦
where Lf h(x(t )) ≡
∂h(x i(t )) f ( x i (t ) ) ∂x i(t )
Lf g (x i(t )) ≡
∂h(x i(t )) g ( x i (t ) ) ∂x i(t )
⎡ 0 ⎢ T ⎢ q b ⎢ H 0p = ⎢ ⋮ ⎢ ⋮ ⎢ ⎢⎣ q TAK − 1b
(3)
The transformation is given as vi(t ) = Lrf h(x i(t )) + Lg Lrf − 1h(x i(t ))ui(t )
0 ⎤ ⎥ ⋱ ⋮ ⎥ ⎥ ⋱ ⋱ ⋮ ⎥ ⋱ 0 ⎥ ⎥ ··· ··· q Tb ⎥⎦ ··· ···
(9)
In the batch process shown in Figure 2, the subsequent output and input data from the entire run at all sampling times are
(4) 9599
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
different weight matrices. Thus, the method is needed because it can automatically select the number of lags L = l that would be included in the multibatch dynamic model for a better multibatch model. This will be discussed in Section 2.3. Note that eq 12 represents the initial state of the current batch run affected by the final state of the previous batch run. Thus, the dynamic batch model is improper if the other factors occur, such as the structure change which is different from eq 1. 2.3. Multibatch Model Identification. The dynamic behavior of the multibatch process model would be identified using historical data. With the data obtained from the process operation, the model for the process can then be identified. To eliminate the computation load for the entire batchwise dynamic structure, the model parameters at each time point can be sequentially estimated. First, assume that the output is related to the current batch only; thus, the initial model is
Figure 2. The repetitive operations of the batch process and the effect of the previous batch.
denoted as yi+1 and vi+1, respectively. The repetitive nature of the batch process implies that the output yi+1 for the current run can be affected not only by the control input vi+1 but also by an initial condition at the beginning of a batch (i+1) that is the result from zi(K), which is the end state of the previous batch (i), as shown in Figure 2. In this case, the final state of the previous run will affect the subsequent batch, so the batchwise relation is zi(0) = Czi − 1(K )
yi = H 0pvi + d
(10)
where C is the batchwise coefficient matrix. Equation 10 is substituted into the state equation repeatedly until time point k = 0. The linear relation (eq 10) is assumed here for easy explanations. The linearization approximation can be applied if the relationship between the initial condition of a current batch and the end condition of the previous batch is not linear. Assume that the last previous L batches, zi−l(0) (l < L) has a significant effect on the ith batch. Thus, zi(0) can be represented by
where each row of the vector can be written as yi (0) = d yi (1) = h1p0vi(0) + d yi (2) = h2p0vi(0) + h1p0vi(1) + d ⋮
zi(0) = CAK − 1bvi − 1(0) + ... + Cbvi − 1(K − 1)
yi (k) =
hkp0vi(0)
yi (K ) =
hKp0vi(0)
+ CAK Cbvi − 2(K − 1) + ...
+ ... + h1p0vi(K − 1) + d
where d can be directly gotten at the initial time point. Using the historical data at the first time point, a predicted 0 ĥpk can be obtained. With successive evaluation of each row, the value of Ĥ p0 is obtained. The model will be validated against the actual data. If the model (eq 15) is incorrect, the value of l is subsequently increased by one batch at each iteration until the predicted output conforms to the process outputs and as such the model becomes
(11)
Equation 11 is substituted into eq 10 to yield (12)
with ⎡ q T(CAK )l − 1CAK − 1b ··· ··· q T(CAK )l − 1Cb ⎤ ⎢ ⎥ ⎢ T ⎥ K l−1 K−1 T K l−1 q A(CA ) Cb ⎥ ⎢ q A(CA ) CA ⎢ b ⎥ ⎥ H Tp = ⎢ ⋮ ⋮ ⎢ ⎥ ⎢ ⎥ ⋮ ⋮ ⎢ ⎥ ⎢ q TAK (CAK )l − 1C ··· ··· q TAK (CAK )l − 1Cb ⎥ ⎢⎣ ⎥⎦ AK − 1b l = 1, ..., L
+
hKp0− 1vi(1)
+ ... + h1p0vi(k − 1) + d
(15)
+ (CAK )L − 1CAK − 1bvi − L(0) + ...
yi = H 0pvi + H1pvi − 1 + H 2pvi − 2 + ... + HLpvi − L + d
+
hkp‐01vi(1) ⋮
+ CAK CAK − 1bvi − 2(0) + CAK CAK − 2bvi − 2(1) + ...
+ (CAK )L − 1Cbvi − L(K − 1)
(14)
yi = H 0pvi + H1pvi − 1 + ... + Hlpvi − l + d
(16)
3. MULTIBATCH MODEL PREDICTIVE CONTROL The proposed control scheme is based on the model developed from the data of the process. The linearized input− output system based on the feedback linearization can be used to form the basis of the model for predictive control. Similar to MPC in the continuous process, the output predictions over a prediction horizon is calculated at each batch, based on a manipulated variable over a control horizon illustrated in Figure 3. Figure 4 shows the corresponding control structure for MB-MPC. The output (y) is compared to the setpoint (yspt), and the tracking error is used by MB-MPC control to predict the future batch input vector (v̂), which contains the input for the entire batch time. Within the batch, at each time k, the input v̂(k) is transformed to ûi+1(k) before it is applied to the process. Using the multibatch model prediction, control for the future batch can be calculated to create the desired operation. For a
(13)
The initial condition of batch processes may not remain the same, but the batch model parameters (eq 12) can be identified using the current operation batches. The batchwise coefficient matrix (C) is embedded into the dynamic batch model parameters (Hpl ) (eq 13) and the parameters can be identified using the operation data. Like the conventional dynamic processes described by a step response model, the serial correlation is taken into account in eq 12 by augmenting the current batch by the previous batches with 9600
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
linearly transformed input v, the control actions for the future M batches are vi = vi − 1 + Δvi vi + 1 = vi − 1 + Δvi + Δvi + 1 ⋮ i+M−1
vi + M − 1 = vi − 1 +
∑
Δvj (17)
j=i
and the predicted model is ŷ t = Hat Δv̂ t + Ht v 0 + Hpt v pt + dt
(18)
where Hat = Ht + Ht F + ··· + Ht F M − 1 ⎡ H 0p ⎢ ⎢ H1p ⎢ Ht = ⎢ ⋮ ⎢ ⋮ ⎢ p ⎢⎣ HM −1 ⎡0 ⎢ ⎢I F = ⎢0 ⎢ ⎢⋮ ⎢⎣ 0 ⎡ Hp ⎢ 1 ⎢ H 2p ⎢ pt H =⎢ ⋮ ⎢ ⋮ ⎢ p ⎢⎣ HM −1
0 ⎤ ⎥ ⋱ ⋱ ⋮⎥ ⎥ ⋱ ⋱ ⋱ ⋮⎥ ⋱ ⋱ 0 ⎥ ⎥ ··· ··· H1p H 0p ⎥⎦ 0 ··· ···
0 ··· ··· 0 ⎤ ⎥ ⋱ ⋱ ⋮⎥ ⋱ ⋱ ⋱ ⋮⎥ ⎥ ⋱ ⋱ 0⎥ ··· 0 I 0 ⎥⎦
p p⎤ H 2p ··· ··· HM − 1 ··· HL ⎥ ⋰ ⋰ ⋰ ⋰ 0⎥ ⎥ ⋰ ⋰ ⋰ ⋰ ⋮⎥ ⋰ ⋰ ⋰ ⋰ ⋮⎥ ⎥ ··· HLp 0 ··· ··· 0 ⎥⎦
Figure 3. The repetition batch operations: (a) controlled and predicted outputs of each batch run and (b) manipulated input of each batch run.
Like conventional MPC, which is characterized in part by its ability to incorporate constraints, it is also possible to incorporate constraints into the proposed control method within the multibatch MPC structure. For the case of constrained control, the bound is applied to the input as well as the outputs. Therefore, the control method can be derived for the case of unconstrained control as well as that of constrained control. 3.1. Unconstrained Control. For the desired trajectories from batch i to batch i + M − 1,
⎡ Δvi ⎤ ⎡ vi − vi − 1 ⎤ ⎢ ⎥ ⎢ ⎥ vi + 1 − vi ⎢ Δvi + 1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⋮ Δv̂ t = ⎢ ⋮ ⎥ = ⎢ ⎢ ⎥ ⎢ ⋮ ⎥ ⎢ ⋮ ⎥ ⎢ ⎥ ⎢⎣ Δvi + M − 1⎥⎦ ⎢⎣ vi + M − 1 − vi + M − 2 ⎥⎦
⎤T ··· y ispt y spt = ⎡⎣ y ispt y ispt +1 + M − 1⎦
⎡ yi ⎤ ⎡ vi − 1⎤ ⎡ vi − 1 ⎤ ⎢ ⎥ ⎢ ⎥ ⎢v ⎥ v ⎢ yi + 1 ⎥ ⎢ i − 1⎥ ⎢ i−2 ⎥ ŷ t = ⎢ ⋮ ⎥ , v 0 = ⎢ ⋮ ⎥ , v pt = ⎢ ⋮ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⋮ ⎥ ⎢ ⋮ ⎥ ⎢ ⋮ ⎥ ⎢y ⎥ ⎢⎣ vi − 1⎥⎦ ⎢⎣ vi − L ⎥⎦ ⎣ i + M − 1⎦ ⎡d ⎤ ⎢ ⎥ ⎢⋮⎥ t d = ⎢⋮⎥ ⎢ ⎥ ⎢⋮⎥ ⎢⎣ d ⎥⎦
(20)
the difference between the prediction and the future outputs is given as
The control objective is to minimize êct so that the prediction can be close to the desired value. It should also decrease
(19) 9601
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
variability in Δv̂t. The control objective can be formulated as a quadratic minimization function mint J ̂ = min( êct )T (ect̂ ) + (Δv̂ t )T λΔv̂ t t Δû
where
(22)
Δû
where λ is the weighting matrix. In order to reach the required target with adjustment to Δv̂t, the objective must be minimized. The derivations are similar to these of the conventional MPC.9,10 Thus, the evaluated control action (Δv̂t) is −1
Δv̂ t = ⎣⎡(Hat )T Hat + λ⎦⎤ (Hat )T et = Lbtbet
(23)
where Lbtb is the controller gain. −1
Lbtb = ⎡⎣(Hat )T Hat + λ⎤⎦ (Hat )T
(24)
and et = y spt − Hpt v pt − Ht v 0 − dt
(25)
Δv̂ is calculated with Lbtb. Based on Δv̂i and the previous batch v̂i, the control action, v̂i+1, at the (i+1)th batch is calculated. As shown in Figure 4, it is transformed to ûi+1(k) and is applied to the process. 3.2. Constrained Control. For constrained control, the upper and lower bounds can be imposed on the input by defining the control input constrain vector, t
T
uub = ⎡⎣uub(1) ··· uub(k) ··· uub(K − 1)⎤⎦ T
ulb = ⎡⎣ulb(1) ··· ulb(k) ··· ulb(K − 1)⎤⎦
(26)
where uub and ulb are the upper and lower limits of the input, respectively. The inherent structure presented in the feedback linearization models leads to a very intuitive control strategy. The predictive control method is based on the linear model between v (input) and y (output). Thus, the constraints (eq 26) should be mapped when the control problem, in terms of the variables v, is designed. Since u and v do not have a linear relationship, before the optimization under the constraints of the feedback linearization models is conducted, the bounds of input (u) must be mapped into constraints on the transformed variables (v) first. The transformation must be performed at each sampling point, because the mapping is state-dependent. For the batch process, the nonlinear mapping of the constraints at each time k is given as vlb(k) = min Lrf h(x(k)) + Lg Lrf − 1h(x(k))u(k)
( u(k)
(
mint J ̂ = min( êct )T (êct ) + (Δv̂ t )T λΔv̂ t t Δv̂
)
(30)
ΩΔv̂ t ≥ ψ
Since (et)Tet is unaffected by Δv̂t, the objective further is simplified: ⎛1 ⎞ mint ⎜ (Δv̂ t )T SΔv̂ t + q Tet ⎟ ⎝ ⎠ Δv̂ 2
(31)
subject to ΩΔv̂ t ≥ ψ
where S = [(Hat)THat + λ] and qT = −(Δût)T (Hat)T. The QP problem is solved using the MatLab optimization toolbox, and readers are referred to the 2004 work of Boyd and Vandenberghe for full detail.12 The QP algorithm, which is subject to the constraints, is employed to search for the optimal control moves. After a range of future control moves is calculated at each sampled instant, control is implemented in a receding horizon manner. Like most of the QP problems, the active-set algorithm for QP requires an initial feasible point. If the initial feasible point does not exist, the feasible set of QP is empty. The control action repeats the previous one until the next time point. While the above derivation is for batch-to-batch control, it can be modified for within-batch control. In such a case, at time k of batch i, the control action from time 0 to k − 1 would be completed, so the past designed control inputs would not be
(27)
and is subject to the constraints (eq 26). Solving the optimization problem is trivial, because the state is known and the objective function is affine in u. For more explanation, the reader is referred to the paper written by Kurtz and Henson.11 Setting the constraints similarly for the outputs, the overall constraint becomes ⎡ Ω1 ⎤ ⎡ ψ1 ⎤ ⎢ ⎥Δv̂ t ≥ ⎢ ⎥ ⎢⎣ Ω 2 ⎥⎦ ⎣ ψ2 ⎦ ΩΔv̂ t ≥ ψ
Δv̂
subject to
)
vub(k) = max Lrf h(x(k)) + Lg Lrf − 1h(x(k))u(k) u(k)
Incorporating the constraints into the objective, the minimization problem becomes
(28) 9602
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
Figure 4. The MB-MPC structure.
changed. The first k row would be removed from yspt (eq 20), leaving only the row and the column from k + 1 to the end time ⎡ ⎡ spt ⎤ ⎤ ⎢ ⎢ yi (k) ⎥ ⎥ ⎢⎢ ⋮ ⎥⎥ ⎢⎢ ⎥⎥ ⎢ ⎢ y spt (K )⎥ ⎥ ⎣ ⎦⎥ i y spt (k) = ⎢ ⎢ ⎥ spt ⎢ y i+1 ⎥ ⎢ ⎥ ⋮ ⎢ ⎥ ⎢ y spt ⎥ ⎣ i+M−1 ⎦
the biosynthesis of penicillin.14 The MB-MPC scheme for the process is as shown in Figure 5. The control method makes use of the previous data batch for prediction control of the current batch. The substrate flow (U) is fed to maintain the substrate concentration (S), which, in turn, affects the growth of biomass (X) that produces the final product penicillin (P). According to Montague et al., the subtract flow is manipulated to control the biomass concentration.13 Control of the biomass concentration under the specific growth rate in fermentation with sufficient supplies of dissolved oxygen enables the maximization of antibiotic production. The reference biomass concentration trajectory is based on the work of Cuthrell and Biegler.14 The process is summarized in Table 1 and the corresponding parameters, in Table 2, with their definitions in the nomenclature. Based on the nonlinear state feedback linearization discussed in Section 2.1, using the transformation
(32)
The same procedure is done to the corresponding disturbance and bias term, yielding
⎛ X ⎞ ⎟⎟U = v μX − ⎜⎜ ⎝ Sf V ⎠
(35)
Similar to eq 6, the linear model, dX =v dt
and the objective function is modified as ct
T
ct
t
T
t
min J (̂ k) = min (ê (k)) (ê (k)) + (Δv̂ (k)) λΔv̂ (k) t
Δv̂ t (k)
(36)
can be used to represent the process. In this fermentation process, the control aims to maintain a biomass concentration at a desired profile by manipulating the subtract flow (U). The process performs the above prescribed tasks repeatedly. However, in the iterative operation, the initial condition at each batch may not be reset to the desired initial condition or near the desired initial condition, because the small amount of components in the previous batch run remains in the tank. In this process, the batchwise relation is
Δv̂ (k)
(34)
With the modified objective function, the control adjustment at the time point of the immediate batch can be formulated in the same manner as that done previously.
4. CASE STUDY To demonstrate the applicability of the proposed control method, a simulation study is carried out on a batch reactor. The simulation demonstrates the effectiveness of the multibatch model in prediction and its effect when used in control. A comparison with the conventional MPC method will also be demonstrated. The process is a fed-batch reactor problem for
mi(0) = αmi − 1(K )
m = X , P, S
(37)
where α = 0.02 is a constant that represents 2% of the biomass, product, and subtract at the ending time (k = K) of the previous batch brought over to the initial time (k = 0) of the 9603
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
Figure 5. Schematic of penicillin fermenter MB-MPC.
Table 1. Process Equations for the Fermentation Process ⎛X ⎞ ⎜ ⎟U ⎝ Sf V ⎠
dX dt
= μX −
dP dt
= ρX − KdegrP −
dSs dt
=−
dV dt
=
μX Yx / s
by a small random perturbation sequence input (U) around the nominal input condition are first collected. Figures 6a and 6b show the subtract flow and the corresponding output biomasses in the five batches from the historical operation data, respectively. The converted variable (v) is obtained from eq 35. Its value has a variation of 0.2 to −0.2 (shown in Figure 6c). The historical data for the initial biomass concentration, which has a slight variation at each batch, are also plotted in Figure 7. The data obtained is used to identify the model of the process. Model identification using these data with the consideration of the multibatch dynamics is carried out. With historical data, it is first assumed that there is only current batch effect, the model parameter Hp0 is obtained and the result is as shown in Figure 8a. The model is still incorrect, indicating that the past batch effect has not been correctly identified. At subsequent iteration, this past effect is increased by one and Hp1 is obtained. The result is as shown in Figure 7b, in which the predicted value is close to the actual value. 4.2. Multibatch Model Predictive Control. The biomass concentration profile of Cuthrell and Biegler11 has been used as the reference trajectory for the purpose of control, as shown in Figure 9. The control method aims to maintain the desired
−
ρX YP / s
⎛ P ⎞ ⎜ ⎟U ⎝ Sf V ⎠
−
msS X Km + S
⎛ + ⎜1 − ⎝
Ss ⎞ U ⎟ Sf ⎠ V
U Sf
μ = μmax
(
S KxX + S
)
⎧ ⎫ S ⎬ ρ = ρmax ⎨ ⎪ ⎪ ⎩ K p + S[1 + (S /K in)] ⎭ ⎪
⎪
next batch. The above models generate the simulation data. The actual values of the batchwise coefficient matrix would not be known in advance, because the proposed control design is based on the identified multibatch model. 4.1. Multibatch Dynamic Model. To build up the initial multibatch dynamic model, the biomass concentration (X) driven Table 2. Parameters for the Fermentation Process parameter μmax ρmax KX KP Kin
value −1
0.11 h 0.0055 g P g−1 X h−1 0.006 g S g−1 X 0.0001 g S l−1 0.1 g S l−1
parameter
value −1
ms Yx/s Yp/s Sf Kdeg r
0.029 g S g X h 0.47 g X g−1 S 1.2 g X g−1 S 500 g S l−1 0.01 h−1 9604
−1
parameter
value
Km X(0) P(0) S(0) V(0)
0.0001 g S l−1 1.5 g l−1 0.0 g l−1 0.0 g l−1 7l
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
Figure 7. Initial biomass concentration at each batch run.
Figure 6. Historical five batch runs: (a) input feed rate of substrate, (b) output biomass, and (c) converted input variable.
biomass concentration profile that the process should follow for a fixed batch time by manipulating the subtract flow. Based on the identified model, control of the process is simulated. In order to show the effectiveness of the proposed MB-MPC, a comparison between MB-MPC and conventional MPC is made. Figure 10a shows the control results from the first batch. Since this is the first batch, the effect of the batchwise dynamics has no bearing on the output and, hence, it is found that
Figure 8. Validation of the linearized model (a) with the current batch only and (b) with one previous batch.
the outputs of both methods are identical. Nevertheless, when the results from the future batch are compared, the control 9605
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
Figure 9. Reference trajectory of biomass concentration.
performance differs greatly. At the third batch (Figure 10b), because the past batch has come into effect and MPC does not take this effect into account, the resultant biomass concentration profile deviates from the reference trajectory. This deviation starts to occur after ∼30 h. The quality of the learned model with the single batch is not accurate and it is certain that the performance of predictive control is inferior. In contrast, MB-MPC takes this effect into account and the tracking is satisfactory. Figure 10c shows that the control results from the fifth batch. It is found that control using the MB-MPC is able to track the desired trajectory accurately, whereas MPC is not very effective. This outcome is not unexpected as the effect of the previous batch exhibits a great influence on the accuracies of the prediction model. 4.3. Comparison of Feedback. The control law (eq 23) is applied without considering the model error. In order to compensate the difference between the measured outputs and the model prediction, the tracking error is corrected via feedback. It is used to correct the control by compensating for the error at the next prediction; thus, Δv̂ t = Lbtb(y spt − Hpt u pt − Ht u 0 − dt − λf etf )
(38)
where the tracking error is T etf = ⎡⎣(ef )T (ef )T ··· (ef )T ⎤⎦ KM × 1
ef = yi − 1 − yî − 1
(39)
Figure 11 shows the result of control from the first, to the third batch, and to the fifth batch. In this first batch, there is no effect of the previous batch, so both control methods are feasible. Further control results show that control at the third batch using MB-MPC has achieved good tracking of the trajectory, but the MPC feedback cannot achieve similar control performance until the fifth batch. This difference in the convergence speed can be explained by the fact that the correction in the proposed MB-MPC method takes into account a prediction that is an inherent characteristic of the process while the feedback in MPC is simply an empirical correction of the process and it justifies the construction of such a model. In the event of disturbance, even when the model is correct, the control performance will not be
Figure 10. Controlled outputs based on MPC and the proposed MB-MPC at (a) batch run 1, (b) batch run 3, and (c) batch run 5.
accurate. Nevertheless, using the multibatch model, with a correction term similar to eq 38, the disturbance could be rejected. 9606
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
Article
degradation, performance assessment is required to determine and measure the capability of control systems and to improve the degradation performance. The current control performance (J), which is defined as T
J(Δu t ) = (y spt − y t ) (y spt − y t ) + (Δu t )T λ(Δu t )T (40)
should be regularly checked to determine if it is close to the benchmark performance (eq 22). Equation 40 is for the real batch operation objective. If J is significantly different from the benchmark, the multibatch dynamic model should be reidentified. The new benchmark then would be used to keep monitoring the next coming operation batch system. Note that the constraint for the case study is 0 ≤ U (t ) ≤ 50
(41)
The boundary of the constraint region is of no particular significance for finding the minimum J, which is subject to eq 41, because the solution to the unconstrained optimization problem in this case is still in the feasible region for the constraint. Thus, the solution to the constrained optimization problem is the same as the solution to an unconstrained optimization problem. Thus, in this case, the benchmark performance of the constraint region using eq 31 can still be based on the computed benchmark performance using eq 22.
5. CONCLUSION A multibatch model has been proposed in this study. Using input−output linearization, the nonlinear process can be modeled by a linear model. This model accounts for the varying batchwise relation, because of the varying condition as a result of the previous batch operation. Based on the multibatch model, a model predictive control method is derived for the cases of constrained and unconstrained control. The proposed method has been applied to a simulation of a batch fermentation of penicillin. The results illustrate the benefits and the requirement of the multibatch model. It shows that the process can be effectively controlled using the proposed method. A comparison between MB-MPC and MPC shows that a faster convergence speed can be achieved as the correction in the proposed MB-MPC method takes into account a prediction that is an inherent characteristic of the process. Several multistage/multiphase statistical models for batch processes have been developed.15 The 2D window data are used to train the dynamic behavior among batches. Such dynamics exist not only within a particular batch but also from batch to batch. The proposed MB-MPC modeling method and the multistage/multiphase models are developed on different bases, but both of their motivations are essentially rooted in the same soil. The main difference is that the multistage/ multiphase models do not consider the remaining effect of the previous batch. The systematical discussions and comparisons would be an interesting topic in our future work. Also, the proposed MB-MPC method is applied to a single-input singleoutput process and the design is based on the batch process, which can be globally linearized via the nonlinear state feedback linearization. If the dynamic batch model described by eq 1 is improper, the more general methodology for the repetitive batch operation will be discussed in our future study. Moreover, many chemical processes have multivariables. For future work, the current framework will be extended to a multivariable process for control purposes.
Figure 11. Controlled outputs based on MPC with feedback tracking errors and the proposed MB-MPC at (a) batch run 1, (b) batch run 3, and (c) batch run 5.
Note that the essential of successful predictive control counts on the accurate model to describe the operating batch processes. Because of the presence of disturbance and process 9607
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608
Industrial & Engineering Chemistry Research
■
Article
(3) Chin, I.; Qin, S. J.; Lee, K. S.; Cho, M. A Two-Stage ILC Technique Combined with Real-Time Feedback for Independent Disturbance Rejection. Automatica 2004, 40, 1913. (4) Brockett, R. W. Feedback Invariants for Nonlinear Systems. In Proceedings of the Seventh Triennial World Congress, Helinski, Finland, June 12−16, 1978; pp 1115−1120 (5) Kravaris, C.; Chung, C. B. Nonlinear State Feedback Synthesis by Global Input/Output Linearization. AIChE J. 1987, 33, 592. (6) Zhao, C.; Wang, F.; Gao, F.; Zhang, Y. Enhanced Process Comprehension and Statistical Analysis for Slow-Varying Batch Processes. Ind. Eng. Chem. Res. 2008, 47, 9996. (7) Park, H. S.; Ryan, W. D.; Paine, H. L. Cleaning Reactors Contaminated with Carboxyl Containing Polymers. U.S. Patent 4,345,949, 1982. (8) Slotine, J. J. E.; Li, W. Applied Nonlinear Control; Prentice−Hall: Englewood Cliffs, NJ, 1991. (9) Camacho, E.; Borbons, C. Model Predictive Control; Springer− Verlag: London, 1999. (10) Ogunnaike, B. A.; Ray, W. H. Process Dynamics Modeling and Control; Oxford University Press: New York, 1994. (11) Kurtz, M. J.; Henson, M. A. Input−Output Linearizing Control of Constrained Nonlinear Processes. J. Process Control 1997, 7, 3. (12) Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press: New York, 2004. (13) Montague, G. A.; Morris, A. J.; Wright, A. R.; Aynsley, M.; Ward, A. Modelling and Adaptive Control of Fed-Batch Penicillin Fermentation. Can. J. Chem. Eng. 1986, 64, 567. (14) Cuthrell, J. E.; Biegler, L. T. Simultaneous Optimization and Solution Methods for Batch Reactor Control Profiles. Comput. Chem. Eng. 1989, 13, 49. (15) Yao, Y.; Gao, F. A Survey on Multistage/Multiphase Statistical Modeling Methods for Batch Processes. Ann. Rev. Control 2009, 33, 172.
AUTHOR INFORMATION
Corresponding Author
*Fax: +886-3-2654199. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■ ■
ACKNOWLEDGMENTS The authors acknowledge financial support for this work from National Science Council, R.O.C. NOMENCLATURE X = biomass concentration (g L−1) P = substrate concentration (g L−1) X = product (penicillin) concentration (g L−1) V = culture volume (L) μ = specific growth rate of biomass (h−1) ρ = specific production rate of penicillin (g P g−1 X h−1) μmax = maximum specific growth rate of biomass (h−1) ρmax = maximum specific production rate (g P g−1 X h−1) KX = Contois saturation constant for substrate limitation of biomass production (g S g−1 X) KP = Monod saturation constant for substrate limitation of produce production (g S L−1) Kin = substrate inhibition constant for product formation (g S L−1) Kdegr = first-order decay rate constant for product (h−1) Km = maintenance constant (g S L−1) ms = maintenance requirement of substrate (g S g−1 X h−1) YX/S = yield of biomass on substrate (g X g−1 S) YP/S = yield of product on substrate (g P g−1 S) yspt = reference trajectory SF = substrate concentration in the feed stream (g L−1) U = substrate flow rate (L h−1) z = transformed state vector k = unit time point l = model horizon K = batch end time A, b, q = state space model coefficients C = batchwise correlation constant x = state vector y = process output u = process input v = linearized input Lbtb = multibatch controller gain matrix hp = model coefficient Hp = model coefficient matrix d = bias dt = bias vector e = tracking error e = tracking error vector
Subscripts
i = batch number ub = upper bound lb = lower bound
■
REFERENCES
(1) Amann, N.; Owens, D. H.; Rogers, E. Iterative Learning Control for Discrete-Time System with Exponential Rate of Convergence. IEEE Proc. Control Theory Appl. 1996, 143, 217. (2) Lee, K. S.; Chin, I. S.; Lee, H. J.; Lee, J. H. Model Predictive Control Technique Combined with Iterative Learning Control for Batch Processes. AIChE J. 1999, 45, 2175. 9608
dx.doi.org/10.1021/ie2020125 | Ind. Eng. Chem. Res. 2012, 51, 9598−9608