Article pubs.acs.org/IECR
A Max−Min Control Problem Arising in Gradient Elution Chromatography Qinqin Chai,*,†,‡ Ryan Loxton,‡ Kok Lay Teo,‡ and Chunhua Yang† †
School of Information Science & Engineering, Central South University, Changsha, China Department of Mathematics & Statistics, Curtin University, Perth, Australia
‡
ABSTRACT: Gradient elution chromatography is an industrial process used to separate and purify multi-component chemical mixtures. In this article, we consider an optimal control problem in which manipulative variables in the chromatographic process need to be determined to maximize separation efficiency. This problem has two nonstandard characteristics: (i) the objective function is nonsmooth, and (ii) each state variable is defined over a different time horizon. The final time for each state variable, the so-called retention time, is not fixed and actually depends on the control variables. To solve this optimal control problem, we first introduce a set of auxiliary decision variables to govern the ordering of the retention times. Then, we approximate the control by a piecewise-constant function and apply a novel time-scaling transformation to map the retention times and control switching times to fixed points in a new time horizon. The retention times and control switching times become decision variables in the new time horizon. On this basis, the optimal control problem is reduced to an approximate nonlinear optimization problem that can be solved using a recently developed exact penalty method. Numerical results show that our approach is both accurate and efficient. The time transformation method used by Jennings et al.1 is based on the substitution t = stf, where t ∈ [0, tf] is the original time variable, s ∈ [0, 1] is the new time variable, and tf is the final time of the process. The same transformation has been successfully applied to solve time-optimal control problems.3,4 When applied to the chromatography optimal control problem, this transformation does not map the retention times to fixed points; it maps only the final time. Furthermore, the retention times do not necessarily coincide with the control switching times in the new time horizon. In fact, the retention times remain variable under this transformation, which makes them very difficult to compute numerically. This is a major disadvantage, as studies show that the productivity of a chromatographic process is highly sensitive to the retention times.5,6 Thus, although difficult, accurate determination of the retention times is crucial for industrial applications.7 In the chromatography optimal control problem, an equality state constraint is imposed at each retention time. Such constraints are called characteristic-time constraints in the optimal control literature.8,9 Computational methods for solving optimal control problems with characteristic-time constraints were developed by Martin8 and Loxton et al.9 These methods, however, assume that the ordering of the characteristic times is fixed and known. In the chromatography optimal control problem, the characteristic times are the retention times, and their ordering depends on the control variables. One approach that can be used to deal with this situation is the bilevel optimization approach,10 whereby the ordering is optimized in the outer level and the control
1. INTRODUCTION Chromatography plays an important role as a separation and purification process in many industrial settings, especially in the preparation of biochemical and pharmaceutical products. A typical chromatographic system consists of a column containing an absorbent, called the stationary phase, and a liquid that flows through the column, called the mobile phase. The absorbent is fixed in the column. The mixture to be separated is injected into the mobile phase and flows through the column. Because the different components in the mixture are attracted to the stationary phase to varying degrees, they travel through the column at different speeds, and thus they exit the column at different times (called peak times or retention times). Therefore, the mixture is gradually separated while moving through the column. The separated components are analyzed by a detector at the outlet of the column. In practice, chromatography is time-consuming, and the purity of the final product must satisfy strict conditions. Thus, it is essential that the chromatographic process be controlled judiciously by varying mobile-phase conditions such as pH value, ionic strength, and flow rate. To achieve a high-quality separation and improve productivity, the minimum duration between successive retention times should be maximized. The problem of choosing a chromatography control scheme that achieves optimum separation was formulated as a max− min optimal control problem by Jennings et al.1 A computational method for solving this problem, based on the control parametrization technique,2 was proposed. This method involves reformulating the max−min objective function into a more convenient form and then approximating the control by a piecewise-constant function, before finally transforming the time horizon into the fixed interval [0, 1]. This yields an approximate mathematical programming problem that can be solved using existing optimization algorithms. © 2012 American Chemical Society
Received: Revised: Accepted: Published: 6137
October 29, 2011 February 28, 2012 March 29, 2012 March 29, 2012 dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research
Article
where τf is the terminal time for the chromatographic process, defined by
variables and retention times are optimized in the inner level. Another possible approach is the grid search algorithm,7 whereby the ordering of the retention times is determined after the controls are chosen. However, both the bilevel approach and the grid search algorithm are inefficient because they involve solving a computationally intensive discrete optimization problem. In this article, we consider the same chromatography optimal control problem as formulated by Jennings et al.1 We propose a new method for reformulating this problem that facilitates accurate determination of the retention times. First, a set of auxiliary decision variables are introduced to govern the ordering of the retention times. Then, after the control variables have been approximated by piecewise-constant functions, a novel time-scaling transformation is used to map the retention times to fixed points in a new time horizon. We then show that the max−min optimal control problem under consideration is equivalent to a minimization problem subject to additional inequality constraints. This minimization problem can be solved using an exact penalty method.11 We conclude the article with two numerical examples.
τf = max {τk} k = 1,..., N
We now introduce the following optimal control problem. Problem P. Given the system defined by eqs 1 and 2, choose the control u: [0, ∞) → p and the corresponding retention times τk, k = 1, ..., N, to maximize objective function in eq 5 subject to the terminal constraints in eq 3 and the control constraints in eq 4. Compared with standard optimal control problems, problem P has several unusual characteristics: (i) The max−min objective function is nonsmooth. (ii) Each state variable is defined over a different time horizon. (iii) The retention times and the ordering of the retention times are not fixed, but instead depend on the control through eq 3. Thus, problem P is a highly nonstandard optimal control problem and cannot be solved directly using standard methods such as the Pontryagin minimum principle12 or state discretization.13,14
2. PROBLEM STATEMENT A gradient elution chromatographic process with N components can be described by the following dynamical system1 xk̇ (t ) = fk [t , xk(t ), u(t )],
t > 0; k = 1, ..., N
3. PROBLEM TRANSFORMATION Let τ = [τ1, ..., τN]T be a vector containing the retention times. Furthermore, let vij, i = 1, ..., N; j = 1, ..., N, be a set of auxiliary decision variables controlling the order of the retention times, where
(1)
with initial conditions xk(0) = 0,
k = 1, ..., N
(2)
⎧1, if component j has the i th earliest retention time vij = ⎨ ⎩ 0, otherwise ⎪
where xk is the chromatography signal corresponding to the kth component, f k is the signal velocity corresponding to the kth component, and u is a vector representing mobile-phase conditions such as pH value, ionic strength, temperature, and flow rate. In the language of control theory, x = [x1, ..., xN]T ∈ N is called the state, and u = [u1, ..., up]T ∈ p is called the control. Each f k: × × p → is assumed to be a given continuously differentiable function. For each k = 1, ..., N, the retention time τk for the kth component is defined by the equality constraint xk(τk) =
∫0
τk
fk [t , xk(t ), u(t )]dt = Lk , k = 1, ..., N
⎪
Then, clearly N
∑ vij = 1,
t ≥ 0; j = 1, ..., p
(7)
and N
∑ vij = 1,
j = 1, ..., N (8)
i=1
(3)
[v1T,
T T
We collect the auxiliary variables into a vector v = ..., vN ] NN T ∈ , where vi = [vi1, ..., viN] . As an example, consider a three-component mixture in which component 2 is the first component to exit the chromatography column, component 1 is the second, and component 3 is the last. Then, τ2 < τ1 < τ3, and thus, v1 = [0, 1, 0]T, v2 = [1, 0, 0]T, and v3 = [0, 0, 1]T. Standard gradient-based optimization techniques cannot handle binary constraints such as eq 6. Thus, we replace the 0−1 constraints on vij by the constraints
(4)
where aj and bj are the lower and upper bounds, respectively, of the jth control variable. Given a piecewise-continuous function u: [0, ∞) → p satisfying eq 4, we can solve the system defined by eqs 1 and 2 to yield a corresponding state trajectory. To ensure that the gradient elution process is effective, the control variables should be chosen so that the minimum duration between successive retention times is maximized. Hence, we introduce the objective function ⎡ (τ − τ )2 ⎤ j i ⎥ J(u) = min⎢ i≠j ⎢ ⎥⎦ τ f ⎣
i = 1, ..., N
j=1
where Lk is the peak height of the chromatography signal xk. Thus, the state variable xk is defined on the time horizon [0, τk]. We assume that the p control variables in the chromatographic process are bounded aj ≤ uj(t ) ≤ bj ,
(6)
⎡N ⎤2 ⎛2 ⎛ ⎞⎥ 1⎞ ⎢ 1 ∑ vij⎜⎝j − j + ⎟⎠ − ⎢∑ vij⎝⎜j − ⎠⎟⎥ = 1 , 3 2 ⎦ 12 j=1 ⎣ j=1 N
i = 1, ..., N
(9)
and 0 ≤ vij ≤ 1,
i = 1, ..., N ; j = 1, ..., N
(10)
The following result shows that eqs 9 and 10 imply vij ∈ {0, 1}.
(5) 6138
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research
Article
Theorem 1. Suppose that vij, i = 1, ..., N; j = 1, ..., N, satisfy eqs 7 and 10. Then, for each i = 1, ..., N, eq 9 holds if and only if there exists a k ∈ {1, ..., N} such that vik = 1 and vij = 0 for all j ≠ k. Proof. Let i ∈ {1, ..., N} be fixed but arbitrary, and assume that vik = 1 and vij = 0 for all j ≠ k. Then
g (y ) = y 2 − y +
⎧ ⎪ = y 2 − ⎨1 + ⎪ ⎩
⎡N ⎤2 ⎛2 ⎛ ⎞⎥ 1⎞ ⎢ 1 ∑ vij⎜⎝j − j + ⎟⎠ − ⎢∑ vij⎜⎝j − ⎟⎠⎥ 3 2 ⎦ j=1 ⎣ j=1
⎡ 1 +⎢ + ⎢⎣ 3
N
= k2 − k + =
⎛ 1 1⎞ − ⎜k − ⎟ ⎝ 3 2⎠
(11)
To prove the opposite implication, we use arguments similar to those used in the proof of lemma 3.1 in the article by Lee et al.15 First, consider the following optimization problem, which we call problem Q ⎡N ⎤2 ⎛2 ⎛ ⎞⎥ 1⎞ ⎢ 1 min ∑ vij⎜j − j + ⎟ − ∑ vij⎜j − ⎟ ⎝ vi1,..., viN 3 ⎠ ⎢⎣ j = 1 ⎝ 2 ⎠⎥⎦ j=1 N
(12)
N
s.t. ∑ vij = 1 (13)
j=1
vij ≥ 0,
j = 1, ..., N
(14)
Note that problem Q has a continuous objective function and a compact feasible region. Hence, it admits at least one optimal solution. The Lagrangian for problem Q is defined by
N
∑ ηjvij j=1
where ρ is the Lagrange multiplier for the equality constraint and ηj ≥ 0, j = 1, ..., N, are the Lagrange multipliers for the inequality constraints. Because the linear independence constraint qualification is clearly satisfied in problem Q, any optimal solution must satisfy the following Kuhn−Tucker conditions16
(15)
and ηk vik = 0,
k = 1, ..., N
(16)
From eq 15, we obtain g (k) − ηk = 0,
j=1
⎤ 1 ⎞⎟ − ρ⎥ ⎥⎦ 2⎠
where the last inequality follows from 0 ≤ vij ≤ 1. Clearly, from cases i and ii, the minimum value of the cost function in problem Q is 1/12. Furthermore, this minimum is achieved only when there exists a k ∈ {1, ..., N} such that vik = 1 and vij = 0, j ≠ k. This completes the proof. On the basis of theorem 1, we can replace the binary constraints in eq 6 with the non-discrete constraints in eqs 7, 9, and 10. As demonstrated later, this reformulation enables the optimal retention time ordering to be determined using an exact penalty method. To proceed, we now use the control parametrization technique2 to approximate problem P by a finite-dimensional optimization problem. This is done by approximating the control u by a piecewise-constant function that switches value at each retention time and at q − 1 times between each pair of
⎡N ⎤ ⎛ ∂3 1 1⎞ ⎛ 1⎞ = k 2 − k + − 2⎢∑ vij⎜j − ⎟⎥⎜k − ⎟ ⎝ ⎠ ⎝ ⎢⎣ j = 1 ∂vik 3 2 ⎥⎦ 2⎠ k = 1, ..., N
⎜
⎡N ⎤2 ⎛2 ⎛ ⎞⎥ 1⎞ ⎢ 1 ∑ vij⎜⎝j − j + ⎟⎠ − ⎢∑ vij⎝⎜j − ⎠⎟⎥ 3 2 ⎦ j=1 ⎣ j=1 ⎛ ⎛ 1⎞ 1⎞ = vik1⎜k12 − k1 + ⎟ + vik 2⎜k 2 2 − k 2 + ⎟ ⎝ ⎝ 3⎠ 3⎠ 2 ⎡ ⎛ ⎞ ⎛ ⎞⎤ 1 1 − ⎢vik1⎜k1 − ⎟ + vik 2⎜k 2 − ⎟⎥ ⎝ ⎣ ⎝ 2⎠ 2 ⎠⎦ ⎛ ⎛ 2 1⎞ 1⎞ = vik1⎜k1 − k1 + ⎟ + (1 − vik1)⎜k 2 2 − k 2 + ⎟ ⎝ ⎠ ⎝ 3 3⎠ 2 ⎡ ⎛ ⎛ 1⎞ 1 ⎞⎤ − ⎢vik1⎜k1 − ⎟ + (1 − vik1)⎜k 2 − ⎟⎥ ⎝ ⎣ ⎝ 2⎠ 2 ⎠⎦ 1 1 2 2 = + (vik1 − vik1 )(k1 − k 2) ≥ 12 12
⎡N ⎤2 ⎛2 ⎛ ⎞⎥ 1⎞ ⎢ 1 3 = ∑ vij⎜j − j + ⎟ − ∑ vij⎜j − ⎟ ⎝ 3 ⎠ ⎢⎣ j = 1 ⎝ 2 ⎠⎥⎦ j=1
− ρ − ηk = 0,
⎛
∑ vij⎝j −
N
N
⎞ ⎛N ⎜ − ρ⎜∑ vij − 1⎟⎟ − ⎠ ⎝ j=1
N
⎤⎫ ⎡N ⎛ 1⎞ ⎪ 2⎢∑ vij⎜j − ⎟⎥⎬y ⎝ ⎢⎣ 2 ⎠⎥⎦⎪ j=1 ⎭
Clearly g(k) = 0 if and only if ηk = 0. Hence, because g(y) is a quadratic function with at most two real roots, no more than two of the multipliers for the inequality constraints in problem Q are zero. Suppose that ηk > 0 for all k = 1, ..., N. Then, by eq 16, vik = 0 for all k = 1, ..., N, contradicting eq 13 in problem Q. Thus, it suffices to consider the following two cases: (i) There exists an integer k1 such that ηk1 = 0 and ηk > 0 for all k ≠ k1 (exactly one of the multipliers is zero). (ii) There exists integers k1 and k2 such that ηk1 = ηk2 = 0 and ηk > 0 for all k ≠ k1, k2 (exactly two of the multipliers are zero). Consider case i. In this case, eq 16 and the equality constraint eq 13 in problem Q imply that vik1 = 1 and vik = 0 for all k ≠ k1. Thus, as shown in eq 11, the optimal cost of problem Q is equal to 1/12. We now consider case ii. In this case, eq 16 implies that vik = 0 for all k ≠ k1, k2. Hence, from the equality constraint eq 13 in problem Q, we have vik2 = 1 − vik1. Thus, the optimal cost of problem Q is
2
1 12
⎤ ⎡N ⎛ 1 1⎞ ⎛ 1⎞ − 2⎢∑ vij⎜j − ⎟⎥⎜y − ⎟ − ρ ⎝ ⎠⎥⎝ ⎢⎣ 3 2 2⎠ ⎦ j=1
k = 1, ..., N
where g: → is a quadratic function defined by 6139
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research
Article qN
successive retention times. Thus, the time horizon is divided into qN subintervals, with q subintervals between each pair of successive retention times, and the control is approximated by a constant value on each subinterval. For each j = 1, ..., p, the control uj is approximated as
uj̃ (s) = uj[t(s)] =
uj(t ) =
∑ i=1
dxk̃ (s) = ds
⎪
⎪
(18)
with initial conditions xk̃ (0) = 0,
J (̃ σ , θ , v) =
=
2⎫ ⎧ ⎪ (t(i + 1)q − tiq) ⎪ ⎬ min ⎨ ⎪ i = 1,..., N − 1⎪ tqN ⎩ ⎭ 2⎫ ⎧ ⎪ (θiq + 1 + θiq + 2 + ··· + θiq + q) ⎪ ⎬ min ⎨ ⎪ i = 1,..., N − 1⎪ θ1 + θ2 + ··· + θqN ⎩ ⎭
For each i = 1, ..., N, exactly one component in the chromatography system has its retention time at t = tiq. Hence, we have the following interior-point constraints N
∑ vik[xk̃ (iq) − Lk] = 0,
where θi = ti − ti − 1 ≥ 0,
i = 1, ..., qN
(21)
It follows from eq 19 that, for s ∈ [l − 1, l], we have l−1
s
ω(η)dη =
(27)
In view of our discussion above, problem P can be approximated by the following optimization problem. Problem P1. Given the system defined by eqs 24 and 25, find vectors σ ∈ qNp, θ ∈ qN , and v ∈ NN , such that the objective function in eq 26 is maximized subject to eqs 7−10, 21, and 23 and the interior-point constraints in eq 27. Although problem P1 is a finite-dimensional optimization problem, it is still difficult to solve because the objective function in eq 26 is nonsmooth. Thus, we introduce a new decision parameter ξ, where
(20)
i=1
i = 1, ..., N
k=1
qN
∑ θχi [i− 1,i) (s)
∑ θj + θl(s − l + 1)
(θiq + 1 + θiq + 2 + ··· + θiq + q)2
j=1
ξ=
Furthermore, for each l = 1, ..., qN
min
i = 1,..., N − 1
l
θ1 + θ2 + ··· + θqN
Clearly, for each i = 1, ..., N − 1, the following inequality is satisfied
∑ θj = ∑ (t j − t j − 1) = tl j=1
(25)
The objective function in eq 5 becomes
(19)
l
k = 1, ..., N
t(0) = 0
where ω: [0, qN] → [0, ∞) is a piecewise-constant function with switching points at the fixed locations s = i, i = 1, ..., qN − 1. We express ω mathematically as
t (l ) =
k = 1, ..., N
i=1
(26)
t(0) = 0
∫0
qN
∑ θifk [t(s), xk̃ (s), σ i]χ[i− 1,i) (s),
(24)
dt ( s) = ω(s) ds
t (s ) =
(23)
dt ( s) = ω(s) ds
Here, 0 = t0 ≤ t1 ≤ ··· ≤ tqN = τf. Furthermore, for each k = 1, ..., N, the switching time tkq (the right end-point of subinterval [tkq−1,tkq]) coincides with one of the retention times. Our objective is to choose the control heights and control switching times in eq 17 appropriately so that the objective function in eq 5 is maximized. Note that the control approximation scheme used in eq 17 is more flexible than the one developed by Jennings et al.,1 which does not allow the switching times to be determined optimally (instead they are prefixed). It is well-known that treating the control switching times as decision variables causes major problems in numerical computation.4,9 Hence, we use the so-called time-scaling transformation to map the switching times to fixed points in a new time horizon. This involves introducing a new time variable s ∈ [0, qN], and then relating s to t through the differential equation
ω(s) =
j = 1, ..., p ; i = 1, ..., qN
Define σ = [(σ1)T, ..., (σqN)T]T ∈ qNp, where σi = [σi1, ..., σip]T ∈ p, i = 1, ..., qN. Furthermore, let θ = [θ1, ..., θqN]T. Under the time-scaling transformation, the system defined by eqs 1 and 2 becomes
(17)
where ti, i = 0, ..., qN, are the control switching times; σji , i = 1, ..., qN; j = 1, ..., p, are the control heights; and χ[ti−1,ti) is the indicator function for the subinterval [ti−1,ti) defined by ⎧1, if t ∈ [ti − 1 , ti) χ[t , t ) (t ) = ⎨ i−1 i ⎩ 0, otherwise
(22)
where the control heights satisfy the constraints aj ≤ σji ≤ bj ,
t ∈ [0, τf ]
j = 1, ..., p
i=1
qN
σjiχ[t , t ) (t ), i−1 i
∑ σjiχ[i− 1,i) (s),
j=1
(θiq + 1 + θiq + 2 + ··· + θiq + q)2
This shows that the time-scaling transformation defined by eqs 19 and 20 maps the control switching times to fixed integers. After the time-scaling transformation has been applied, the control variables defined in eq 17 are written as
θ1 + θ2 + ··· + θqN
≥ξ
Because θ1 + ··· + θqN > 0, this inequality can be rewritten as 6140
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research ⎛ iq + q ⎞2 ξ ∑ θk − ⎜⎜ ∑ θl ⎟⎟ ≤ 0, ⎝ l = iq + 1 ⎠ k=1
Article
penalty parameter η is large, the final term in eq 30 forces ε to be small, which, in turn, causes the middle term in eq 30 to penalize constraint violations very severely. Hence, minimizing the penalty function will likely lead to feasible points of problem P2. With this in mind, we introduce the following penalty problem for problem P2. Problem P3. Given the system defined by eqs 24 and 25, find vectors σ ∈ qNp, θ ∈ qN , v ∈ NN and parameters ξ and ε such that the penalty function in eq 30 is minimized subject to the bound constraints in eqs 10, 21, and 23 and the interiorpoint constraints in eq 27. Problem P3 is an optimal parameter selection problem with decision parameters σ, θ, v, ξ, and ε and the interior-point constraints in eq 27. Unlike the nonsmooth approximate problem derived by Jennings et al.,1 problem P3 can be solved using standard computational methods such as those in Chapter 5 of Teo et al.2 These computational methods were implemented in the optimal control software package MISER.18 To solve problem P3, MISER requires gradient formulas for both the objective function and the constraint functions. The gradients of Jη̂ are given by
qN
i = 1, ..., N − 1 (28)
Thus, problem P1 is equivalent to the following problem. Problem P2. Given the system defined by eqs 24 and 25, find vectors σ ∈ qNp, θ ∈ qN , and v ∈ NN and the parameter ξ to minimize the cost function J (̂ σ , θ , v, ξ) = −ξ
(29)
subject to eqs 7−10, 21, 23, 27, and 28.
4. COMPUTATIONAL METHOD Problem P2 is a finite-dimensional optimization problem with nonlinear equality and nonlinear inequality constraints. The major difficulty with solving this problem is that the constraints in eqs 7−10 force vij to be binary decision variables and, thus, the feasible region for problem P2 is disjoint. Standard gradient-based optimization methods usually fail miserably when applied to problems with disjoint feasible regions. In this section, we show how to solve problem P2 using an exact penalty method. First, let ⎛ mq + q ⎞2 gm(θ , ξ) = ξ ∑ θk − ⎜⎜ ∑ θl ⎟⎟ , ⎝ l = mq + 1 ⎠ k=1
∂Jη̂ ∂θi
qN
N−1
= 2ε−α
∑ [max{0, gm(θ , ξ) − ε γ κm}] m=1
⎛ mq + q ⎞ × [ξ − 2⎜⎜ ∑ θk ⎟⎟χ[mq + 1, mq + q] (i)] ⎝ k = mq + 1 ⎠
m = 1, ..., N − 1
Furthermore, for each m = 1, ..., N, let ∂Jη̂
N
hm(v) =
∑ vmj − 1
∂vij
j=1 N
hm̂ (v) =
∑ ⎢hm(v)
⎣ m=1 ⎢
+ hm̂ (v)
∑ vim − 1
⎡
N
= 2ε−α
∂hm(v) ∂vij
∂h ̃ (v) ⎤ ∂hm̂ (v) + hm̃ (v) m ⎥ ∂vij ∂vij ⎥⎦
i=1
∂Jη̂
⎡N ⎤2 ⎛ ⎞ ⎛ ⎞ 1 1 1 2 ⎢ hm̃ (v) = ∑ vmj⎜j − j + ⎟ − ∑ vmj⎜j − ⎟⎥ − ⎝ ⎠⎥ ⎝ ⎠ ⎢ 3 2 12 ⎣ j=1 ⎦ j=1 N
∂ξ
k=1
∂Jη̂ ∂ε
∂Jη̂
(30)
∂σji
∑ {−γκmε γ− 1[max{0, gm(θ , ξ) − ε γ κm}]} m=1 −α − 1
[Δ1(θ , ξ , ε) + Δ2 (v)] + βηε β − 1
=0
where
N−1
∑ [max{0, gm(θ , ξ) − ε γ κm}]2
⎧ ⎤ ⎡N ⎛ ⎛ 1⎞ 1⎞ 1⎞ ⎪⎛ 2 ∂hm̃ (v) ⎪ ⎝⎜j − j + ⎠⎟ − 2⎜⎝j − ⎟⎠⎢∑ vij⎜⎝j − ⎟⎠⎥, if i = m 3 2 ⎢⎣ j = 1 2 ⎥⎦ =⎨ ∂vij ⎪ ⎪ 0, if i ≠ m ⎩
m=1 N
∑ {[hm(v)]2
N−1
= 2ε−α − αε
where ε > 0 is a new decision variable and
Δ2 (v) =
m=1
× ∑ θk}
Jη̂ (σ , θ , v, ξ , ε) = −ξ + ε−α[Δ1(θ , ξ , ε) + Δ2 (v)]
Δ1(θ , ξ , ε) =
∑ {[max{0, gm(θ , ξ) − ε γ κm}]
qN
We use an exact penalty approach, recently developed by Yu et al.,11 to handle the equality constraints in eqs 7−9 and the inequality constraints in eq 28. This approach has previously proved very effective at solving semi-infinite and discrete optimization problems.17,20 Consider the following penalty function
+ ηε β
N−1
= −1 + 2ε−α
+ [hm̂ (v)]2 + [hm̃ (v)]2 }
m=1
and
Here, α > 0; β > 2; γ > 0; and κm ∈ (0, 1), m = 1, ..., N − 1, are fixed constants, and η > 0 is the penalty parameter. Note that Δ1(θ,ξ,ε) measures violations in eq 28, whereas Δ2(v) measures violations in eqs 7−9. The idea is that, when the
∂hm(v) ⎧1, if i = m =⎨ ∂vij ⎩ 0, if i ≠ m 6141
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research
Article
∂hm̂ (v) ⎧1, if j = m =⎨ ∂vij ⎩ 0, if j ≠ m
solution of problem P2 as the penalty parameter approaches infinity.11,20 Hence, one can obtain an approximate solution of problem P2 by solving problem P3 for large η. A corresponding suboptimal control for problem P can then be constructed according to eq 17. In the next section, we use this approach to solve two examples.
⎪ ⎪
We now consider the gradient of the interior-point constraints in eq 27. First, define N
Φm[x(̃ mq), v] =
∑ vmk[xk̃ (mq) − Lk],
m = 1, ..., N
5. NUMERICAL EXAMPLES 5.1. Example 1. Jennings et al.1 have described a linear gradient elution chromatographic process for separating four protein solutes. The retention time for each protein solute is controlled by adjusting the ionic strength of the mobile-phase composition. The rate of change of each chromatography signal is given by the following system functions
k=1
That is, Φm is the left-hand side of the interior-point constraints in eq 27 (with index i replaced by m). The partial derivatives of Φm with respect to ξ, ε, and vij can be computed directly as ∂Φm ∂Φm =0 = ∂ε ∂ξ
and ⎧ if i = m ⎪ xj̃ (mq) − Lj , ∂Φm =⎨ ⎪ ∂vij if i ≠ m ⎩ 0,
The partial derivatives with respect to σ and θ, however, are more difficult because σ and θ influence Φm implicitly through the dynamic system defined by eqs 24 and 25. These partial derivatives can be computed using the gradient formulas reported in Chapter 5 of Teo et al.2 First, define the Hamiltonian function as
∑ λkm(s)ω(s)fk [t(s), xk̃ (s), ũ(s)] + λNm+ 1(s)ω(s) k=1
where the costate functions λmk : → satisfy the dynamic systems ∂H m λ̇ k (s) = − m , ∂xk̃
λkm(mq)
f2 [u(t )] =
u(t ) 1.4u(t ) + 5.4
f3 [u(t )] =
u(t ) 1.4u(t ) + 7.2
f4 [u(t )] =
[u(t )]2 1.4[u(t )]2 + 9
xk(τk) = 1,
s ∈ [0, mq]; k = 1, ..., N
k = 1, 2, 3, 4
(31)
We consider the problem of maximizing the objective function in eq 5 with N = 4 subject to the interior-point constraints in eq 31. This problem is in the form of problem P. We use the method outlined in sections 3 and 4 with q = 1 to approximate problem P by problem P3. Problem P3, the penalty problem, can be solved using the software package MISER. The parameters in the penalty function for problem P3 are set as follows:
= vmk
and ∂H m λ̇ N + 1(s) = − m , ∂t
[u(t )]2 1.4[u(t )]2 + 3.6
where the control is subject to the bound constraints 0.5 ≤ u(t) ≤ 2. A chromatography signal of 1 indicates that the protein solute has left the chromatography column. Thus, the state variables must satisfy the following interior-point constraints at the retention times
N
Hm =
f1 [u(t )] =
s ∈ [0, mq]
λNm+ 1(mq) = 0
α = 1.5,
Then, it follows from theorem 5.2.1 in Teo et al.2 that ⎧ i ∂Hm ds , if i ≤ mq ⎪ ∂Φm ⎨ i − 1 ∂uj̃ i = ∂σj ⎪ if i > mq ⎩ 0,
β = 2.2,
γ = 3,
κ1 = κ2 = κ3 = 0.3
We start by solving problem P3 using MISER with η = 10. Then, we increase η and re-solve problem P3 using the previous solution as the initial guess. We continue to increase η until η = 30, at which point the objective function value is 0.8205 with ε = 0.0024. The objective function value is slightly better than the result reported by Jennings et al.1 The largest violation of the interior-point constraints in eq 27 is 6 × 10−4, which is smaller than the violation of 1.07 × 10−3 reported in Jennings et al.1 We next consider problem P3 with q = 3. Thus, the time horizon is divided into 12 subintervals. We again use MISER to solve problem P3. The optimal trajectories and optimal control computed by MISER are shown in Figure 1. The optimal control values and the corresponding time-scaling transformation parameters are reported in Tables 1 and 2, respectively. The optimal objective function value is 0.8276 when ε = 0.001, η = 300, and τ = [2.7507, 5.9795, 9.1807, 12.3820]. As expected, the optimal value for q = 3 is higher than the optimal value for q = 1.
∫
and ⎧ i ∂Hm ⎪ ds , if i ≤ mq ∂Φm = ⎨ i − 1 ∂ω ∂θi ⎪ if i > mq ⎩ 0,
∫
The gradient formulas given above for the penalty and constraint functions can be combined with a standard optimization algorithmfor example, a conjugate gradient method or sequential quadratic programming19to solve problem P3 as a nonlinear programming problem. The optimal control software MISER does this automatically. It can be shown that a local solution of problem P3 converges to a local 6142
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research
Article
K 2[I(t ), z(t )] = 5.48 × 106(0.1867I −1)−26.9 + 13.8z − 1.15z K3[I(t ), z(t )] = 2.21 × 106(0.1867I −1)−24.7 + 12.9z − 1.1z
2
2
The controls are bounded by 0.405 ≤ I(t ) ≤ 0.495,
4.8 ≤ z(t ) ≤ 5.2
A protein species is considered to be separated when its corresponding chromatography signal reaches 1. That is, the state variables must satisfy the interior-point constraints xk(τk) = 1,
k = 1, 2, 3
Our aim is to solve problem P, which involves maximizing the objective function in eq 5 for the system with N = 3 subject to the interior-point constraints given above. Using the process outlined in sections 3 and 4 with q = 3, we approximate problem P by problem P3. The parameters for the exact penalty function are the same as in example 1. Using MISER to solve problem P3, we obtained the optimal trajectories and optimal controls shown in Figure 2. The
Figure 1. Optimal controls and states for example 1.
Table 1. Optimal Control Values for Example 1 parameter
value
parameter
value
σ1 σ2 σ3 σ4 σ5 σ6
1.6314 1.6324 1.6326 0.8729 0.8730 0.7901
σ7 σ8 σ9 σ10 σ11 σ12
0.5000 0.5000 0.5000 0.5644 0.5645 0.5645
Table 2. Optimal Interval Durations for Example 1 parameter
value
parameter
value
θ1 θ2 θ3 θ4 θ5 θ6
0.9347 0.9153 0.9007 1.0740 1.0777 1.0770
θ7 θ8 θ9 θ10 θ11 θ12
1.0801 1.0608 1.0603 1.0541 1.0728 1.0743
Figure 2. Optimal results for example 2.
optimal objective function value is 1.8673, when η = 200, ε = 2.24895 × 10−4, and τ = [5.3029, 10.8293, 16.3557]. The largest violation of the interior-point constraints in eq 27 is 1 × 10−5. All of the other constraints in problem P2 are satisfied.
The equality and inequality constraints in eqs 7−10, 21, 23, and 28 in problem P2 are all satisfied. Compared with the results obtained by Jennings et al.,1 our optimal objective function value is better. Moreover, the largest violation of the interior-point constraints in eq 27 is 8 × 10−4, which is slightly smaller than the violation in Jennings et al.1 5.2. Example 2. We now apply our optimal control method to a more complicated system with multiple inputs. Consider the ion-exchange chromatographic process for separating three protein species described by Kaltenbrunner et al.21 The peak time for each protein solute is controlled by adjusting the ionic strength (I) and the pH value (z) of the mobile-phase composition. The rate of change of each chromatography signal is given by the following system functions
6. CONCLUSIONS In this article, we have presented a new computational method for solving a challenging max−min optimal control problem arising in gradient elution chromatography. The state variables in this problem represent the chromatography signals of the mixture components, and each signal is required to reach its desired target at a different retention time. Our method is based on an approximation scheme whereby the control is approximated by a piecewise-constant function. We use a novel time-scaling transformation to map the retention times to fixed points and another transformation to convert the max− min objective function into a smooth function. An exact penalty method is then used to solve the resulting approximate optimization problem. The numerical simulations in section 5 show that our new method gives better results than the previous method proposed by Jennings et al. 1 More importantly, our new method is capable of accurately computing the retention times, and it can also be applied to more general problems of larger dimension.
fk [I(t ), z(t )] ⎡ ⎤−1 (1 − εb)εm (1 − εb)(1 − εm) Kk[I(t ), z(t )]⎥ , = ⎢1 + + εb εb ⎣ ⎦ k = 1, 2, 3
where εb = 0.25, εm = 0.8, and the affinity distribution coefficients are defined by K1[I(t ), z(t )] = 0.49 × 106(0.1867I −1)−18.5 + 9.9z − 0.77z
2
6143
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
Industrial & Engineering Chemistry Research
■
Article
(20) Yu, C.; Li, B.; Loxton, R.; Teo, K. L. Optimal Discrete-valued Control Computation. J. Global Optim., Published online Feb. 8, 2012. (21) Kaltenbrunner, O.; Giaverini, O.; Woehle, D.; Asenjo, J. A. Application of Chromatographic Theory for Process Characterization Towards Validation of an Ion-Exchange Operation. Biotechnol. Bioeng. 2007, 98, 201−209.
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was partially supported by a grant from the Australian Research Council and Grant 61025015 from the China National Science Fund for Distinguished Young Scholars.
■
REFERENCES
(1) Jennings, L. S.; Teo, K. L.; Wang, F. Y.; Yu, Q. Optimal Protein Separation. Comput. Chem. Eng. 1995, 19, 567−573. (2) Teo, K. L.; Goh, C. J.; Wong, K. H. A Unified Computational Approach to Optimal Control Problems; Longman Scientific and Technical: Essex, U.K., 1991. (3) Woinaroschy, A.; Isopescu, R. Time-Optimal Control of Dividing-Wall Distillation Columns. Ind. Eng. Chem. Res. 2010, 49, 9195−9208. (4) Lee, H. W. J.; Teo, K. L.; Jennings, L. S.; Rehbock, V. Control Parametrization Enhancing Technique for Time Optimal Control Problems. Dyn. Syst. Appl. 1997, 6, 243−261. (5) Toumi, A.; Hanisch, F.; Engell, S. Optimal Operation of Continuous Chromatographic Processes: Mathematical Optimization of the VARICOL Process. Ind. Eng. Chem. Res. 2002, 41, 4328−4337. (6) García, M. G.; Balsa-Canto, E.; Banga, J. R. Dynamic Optimization of a Simulated Moving Bed (SMB) Chromatographic Separation Process. Ind. Eng. Chem. Res. 2006, 45, 9033−9041. (7) Yin, X. Y.; Yao, W. F.; Hu, Y. Z. Optimization of Gradient Separation for Chromatographic Fingerprint of Herbal Medicine: A Predictive Model Combined with Grid Search. J. Chromatogr. Sci. 2008, 46, 722−729. (8) Martin, R. B. Optimal Control Drug Scheduling of Cancer Chemotherapy. Automatica 1992, 28, 1113−1123. (9) Loxton, R. C.; Teo, K. L.; Rehbock, V. Optimal Control Problems with Multiple Characteristic Time Points in the Objective and Constraints. Automatica 2008, 44, 2923−2929. (10) Farhadinia, B.; Teo, K. L.; Loxton, R. C. A Computational Method for a Class of Non-standard Time Optimal Control Problems Involving Multiple Time Horizons. Math. Comput. Model. 2009, 49, 1682−1691. (11) Yu, C.; Teo, K. L.; Zhang, L.; Bai, Y. A New Exact Penalty Function Method for Continuous Inequality Constrained Optimization Problems. J. Ind. Manage. Optim. 2010, 6, 895−910. (12) Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.; Mishchenko, E. F. The Mathematical Theory of Optimal Processes; Gordon and Breach Science Publishers: New York, 1986. (13) Hager, W. W. Runge-Kutta Methods in Optimal Control and the Transformed Adjoint System. Numer. Math. 2000, 87, 247−282. (14) Kaya, C. Y.; Martínez, J. M. Euler Discretization and Inexact Restoration for Optimal Control. J. Optim. Theory Appl. 2007, 134, 191−206. (15) Lee, H. W. J.; Teo, K. L.; Cai, X. Q. An Optimal Control Approach to Nonlinear Mixed Integer Programming Problems. Comput. Math. Appl. 1998, 36, 87−105. (16) Nocedal, J.; Wright, S. Numerical Optimization, 2nd ed.; Springer: New York, 2006. (17) Yu, C.; Teo, K. L.; Bai, Y. An Exact Penalty Function Method for Nonlinear Mixed Discrete Programming Problems. Optim. Lett., published online Aug. 30, 2011. DOI: 10.1007/s11590-011-0391-2. (18) Jennings, L. S.; Fisher, M. E.; Teo, K. L.; Goh, C. J. MISER3 Optimal Control Software: Theory and User Manual; University of Western Australia: Perth, Australia, 2004. (19) Luenberger, D. G.; Ye, Y. Linear and Nonlinear Programming, 3rd ed.; Springer: New York, 2008. 6144
dx.doi.org/10.1021/ie202475p | Ind. Eng. Chem. Res. 2012, 51, 6137−6144
