Homotopy with Second-Order Correction Based Backtracking Method

Article pubs.acs.org/IECR

Homotopy with Second-Order Correction Based Backtracking Method for Chemical Process Simulation Weifeng Chen,† Zhijiang Shao,*,‡ Lingyu Zhu,§ and Xi Chen‡ †

College of Information Engineering, Zhejiang University of Technology, Hangzhou 310023, China State Key Laboratory of Industrial Control Technology, Department of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China § College of Chemical Engineering and Materials Science, Zhejiang University of Technology, Hangzhou 310014, China ‡

ABSTRACT: The design and management of chemical manufacturing processes, including performance assessment and product quality evaluation, always involve frequent simulation. The equation-oriented approach is one of the most competitive strategies in process simulation, but this approach is difficult to implement because of the complexity and nonlinearity of chemical processes. In this study, a homotopy-based backtracking method (HBM) for simulation is improved with a secondorder correction, which can be calculated with low computational cost. The theoretical foundation for the proposed homotopy with second-order correction based backtracking method (SOC-HBM) is also established. Simulation runs of production of secreted protein in a fed-batch reactor, a penicillin fermentation process, and a binary distillation column system are demonstrated using SOC-HBM.

1. INTRODUCTION Process simulation is vital to the design and management of chemical manufacturing processes. Through simulation, plant design and operation for productivity, performance, and product quality can be evaluated in a safe, reliable, and economical manner. The equation-oriented approach is one of the most competitive strategies in process simulation because of its unique advantages, such as direct description and advanced simultaneous solving strategy. The equation-oriented approach assembles all the variables and equations describing the process model together and then solves the system of equations. Given the complexity and nonlinearity of chemical processes, process simulation using the equation-oriented approach is a challenging task. The essence of the equation-oriented approach is to solve a system of nonlinear equations. The strategies for solving a system of nonlinear equations could be classified into the following two groups: the derivative-free and the derivativebased. Mousa et al.1 transformed the system of nonlinear equations into a nonlinear programming problem (NLP) with an additional parameter that was used to define the initial precision of the system, after which a coevolutionary algorithm is implemented for solving the resulting NLP. Hirsch et al.2 proposed an enhanced continuous greedy randomized adaptive search procedure (GRASP) for finding all the roots of a system of nonlinear equations. Ren et al.3 introduced an efficient genetic algorithm using symmetric and harmonious individuals for solving a system of nonlinear equations. These heuristic rules-based derivative-free approaches have good global convergence, but are not suitable for large-scale process design and simulation, and their solving efficiencies may not meet the needs of practical applications. The well-known derivativebased strategy for solving system of nonlinear equations is Newton’s method, which has local Q-quadratic convergence. The assumption that the starting point should be sufficiently © 2014 American Chemical Society

close to the solution restricts the widespread application of Newton’s method. To reduce the dependence of the starting point, line search and trust region strategies were introduced to improve the global convergence of Newton’s method.4−6 However, the above-mentioned strategies suffer from one shortcoming, wherein these strategies may converge to a local minimum of the merit function instead of a solution of the nonlinear system.7 Different from a line search or trust region based Newton’s method, the homotopy approach does not rely on the descent of a merit function and is more likely to lead to solutions in difficult cases.8 Jalali et al.9 propounded the homotopy method in the complex domain based on bifurcation theory for solving the nonlinear equations of phase equilibrium and stability analysis. Malinen et al.10 utilized homotopy parameter bounding to improve the robustness of the homotopy method in finding multiple steady state solutions of chemical processes. Considering that the homotopy method has good global convergence under weak assumptions and retains the advantages of Newton’s method, certain problems are still encountered when it is used in practical applications. The critical issues are how to design a proper homotopy function and define a suitable sequence of homotopy parameter values to guarantee that the solution for each parametric system of nonlinear equations can be found and the whole procedure can converge rapidly. Zhu et al.11 presented a homotopy-based backtracking method (HBM), which heads toward the target and backtracks at failure. Chen et al.12 proposed a sensitivity embedded homotopy-based backtracking method (S-HBM) by introducing sensitivity analysis for a system of nonlinear Received: Revised: Accepted: Published: 15080

March 17, 2014 August 27, 2014 September 8, 2014 September 8, 2014 dx.doi.org/10.1021/ie501139g | Ind. Eng. Chem. Res. 2014, 53, 15080−15088

Industrial & Engineering Chemistry Research

Article

are still utilized to represent the homotopy based backtracking strategy and sensitivity embedded homotopy based backtracking method with homotopy (3) in this paper.] Different from the traditional homotopy path tracking method that increases the homotopy parameter gradually from 0 to 1, the backtracking strategy makes the solving procedure always lead toward the target and backtracks at failure, and can be described as follows. Backtracking Strategy. The backtracking strategy includes five steps: Step 1. Set θc = 0 (current homotopy parameter), θt = 1 (target homotopy parameter), and x0 = x0 (given initial guess). Step 2. Solve h(x,θt) = 0 using Newton-based method, with x0 as the starting point. Step 3. If step 2 fails, proceed to step 4. Otherwise, proceed to step 5. Step 4. If θt − θc ≥ η (a predefined number), calculate θt = θc + (θt − θc)/λ and return to step 2. λ is a backtracking constant. Otherwise, the backtracking strategy fails. Stop. Step 5. If θt is not equal to 1, θc = θt, x0 = x(θt), and θt = 1 and return to step 2. Otherwise, the backtracking strategy is successful. Stop. The following assumptions are constructed to obtain the properties of backtracking strategy. Assumption A. Assumption A consists of three parts.

equations, which can provide derivative information on the roots with very little extra computational costs. As long as the homotopy parameter does not reach the value of 1, S-HBM corrects the solution of current instance of homotopy function with first order derivative information and then utilizes it as the starting point for the instance of a homotopy function whose homotopy parameter is 1, i.e., the original system of nonlinear equations. When the curvature of the solution curve with respect to the homotopy parameter is large, first-order correction is not sufficient. In the present work, the secondorder derivative with respect to the homotopy parameter is approximated. Homotopy with a second-order correction (SOC) based backtracking strategy is employed. The paper is organized as follows. Section 2 analyzes the sensitivity of the homotopy function. The homotopy with second-order correction based backtracking method is developed by applying the sensitivity calculation, and the theory foundation is established. Section 3 demonstrates two case studies to illustrate the performance of the proposed method, including simulations of two continuous stirred tank reactors in series and an argon side arm column in the cryogenic air separation unit. Section 4 concludes this paper.

2. HOMOTOPY WITH SOC BASED BACKTRACKING METHOD The chemical processes described with an equation oriented mode for a steady state situation usually could be defined as follows: f (x ) = 0

A1. All the pairs (x(θ),θ), θ ∈ [0,1], belong to a bound and closed set D. A2. h(x,θ) is three times differentiable with respect to x and θ in set D. A3. ∇xh(x(θ),θ) is nonsingular, and the matrix norm of [∇xh(x(θ),θ)]−1 is less than τ for θ ∈ [0,1], where τ is a positive constant. Lemma 1. Under assumptions A1 and A2, the matrix norms of ∂h(x(θ),θ)/∂θ, ∂2h(x,(θ),θ)/(∂θ ∂x), and ∂2h(x,(θ),θ/∂x2 are bound by M > 0 for all pairs (x(θ),θ), θ ∈ [0,1]. Proof. Under assumptions A1 and A2, since D ⊂ Rn+1 (finite dimension space) is bound and closed; D is a compact set. ∂h(x,θ)/∂θ, ∂2h(x,θ)/(∂θ ∂x), and ∂2h(x,θ)/∂x2 are continuous functions, so the images of the functions in D are also compact sets. Hence, they are bounded. ■ Under assumptions A2 and A3, x(θ) is three times differentiable with respect to θ according to the Classical Implicit Function Theorem and then it could be expanded at θ = θc by Taylor expansion, resulting in the following equation:

(1)

with the unknown vector x ∈ Rn, which consists of temperature, pressure, flow rate, composition variables, and so on. f:Rn → Rn is a continuously differentiable function vector, which consists of the so-called MESH equations including mass balances, equilibrium, summation of compositions, and heat balances. 2.1. Homotopy-Based Backtracking Strategy. Homotopy is defined according to the characteristics of process simulation under variable operation conditions as follows:11,12 h(x , θ ) = f (x , α(θ )) = 0

(2)

where α(0) is the base operation condition, and α(1) is the target operation condition. The homotopy parameter represents the relative distance between the base operation condition and the target operation condition. The issue of (2) is that the simulation solution under the base operation condition cannot be always obtained directly. In this paper, homotopy is still considered as the convex linear combination of real functions with respect to the homotopy parameter θ as shown in the following equation: h(x , θ ) = θf (x) + (1 − θ )g (x) = 0

x(θ ) = x(θc) +

dx 1 d2x (θc) (θ − θc) + (θc) (θ − θc)2 dθ 2 dθ 2

+ O(|θ − θc|3 )

(3)

(4)

Differentiating the homotopy (2) with respect to homotopy parameter θ to obtain the following equation

where f(x) is the system of nonlinear equations (1) to be solved, and g(x) is the auxiliary function for which a solution is known or is easily found. The most widely cited g(x) functions are the fixed-point function g(x) = (x − x0), affine function g(x) = f ′(x0)(x − x0), and Newton function g(x) = f(x) − f(x0),where x0 is an initial guess.13 Here, h(x,0) = 0 has the same solution as g(x) = 0, and h(x,1) = 0 is equivalent to f(x) = 0. When the homotopy parameter is gradually varied, a solution path x(θ) satisfying h(x(θ),θ) = 0 is tracked from the starting point to the solution of (1). Backtracking strategy is stated based on homotopy (3). [The abbreviations HBM and S-HBM

∇x h(x(θ ), θ )

dx (θ ) + ∇θ h(x(θ ), θ ) = 0 dθ

dx (θ ) = −(∇x h(x(θ ), θ ))−1 ∇θ h(x(θ ), θ ) dθ

(5)

(6)

From formula 4, there exists a neighborhood N(θc,r), r > 0, when θ ∈ N(θc,r); the following formulation is satisfied according to Lemma 1. 15081

dx.doi.org/10.1021/ie501139g | Ind. Eng. Chem. Res. 2014, 53, 15080−15088

Industrial & Engineering Chemistry Research x(θ ) − x(θc) ≤ (∇x h(x(θ ), θ ))−1

Article

∇θ h(x(θ ), θ )

x0 = x(θc) +

|θ − θc| + o(|θ − θc|) < (τM + 1)|θ − θc|

Step 4. Solve h(x,θt) = 0 using Newton-based method, with x0 as the starting point. Step 5. If step 4 fails, proceed to step 6. Otherwise, proceed to step 7. Step 6. If θt − θc ≥ η (a predefined number), calculate θt = θc + (θt − θc)/λ and return to step 4; λ is a backtracking constant. Otherwise, the backtracking strategy fails. Stop. Step 7. If θt is not equal to 1, θc = θt, x(θc) = x(θt), θt = 1, and return to step 2. Otherwise, the backtracking strategy is successful. Stop. In the following, the properties of SOC-HBM are analyzed. The error between x(θ) and the approximation x̃(θ) is considered. Calculating B̅ (θc) using formula 10, ∇xh(x,θ ̅ c + δ) should be nonsingular, and this can be guaranteed by Theorem 1. Lemma 2. G ∈ Rn×n, I ∈ Rn×n is a identity matrix. If the matrix norm of G is less than 1, then I − G is nonsingular.14 Theorem 1. Under Assumption A, for an arbitrary given ρ, 0 < ρ < 1, there exists a constant δ0 > 0. If 0 < δ < δ0, then −1 ∇xh(x,θ ̅ c + δ) is nonsingular and ∥[∇xh(x,θ ̅ c + δ)] ∥ < τ/(1 − ρ). Proof. According to assumption A2, h(x,θ) is three times differentiable with respect to x. ∇xh(x,θc + δ) can be expanded at (x(θc + δ), θc + δ).

(7)

The backtracking strategy can make x(θc) sufficiently close to the solution of (3) by decreasing |θ − θc|, and then Newtonbased method can always solve (3) successfully with x(θc) as the initial guess. Hence, backtracking strategy can implement the solving procedure step by step in theory under Assumption A. However, the gradient and even the curvature of x(θ) with respect to θ may become very large, and the backtracking strategy should make |θ − θc| small enough to guarantee the success of the Newton-based method. The efficiency of backtracking strategy becomes low and unacceptable. 2.2. Homotopy with SOC Based Backtracking Method (SOC-HBM). In backtracking strategy, the solution of (3) corresponding to the current homotopy parameter θc is utilized as the starting point directly. Formula 7 shows that x(θc) is not a suitable starting point for Newton’s method when |θ − θc| is large. However, formula 7 could be improved using the gradient and curvature information according to formula 4. Using formula 6, (dx/dθ(θc) could be calculated easily. However, for the second order derivative information (d2x/ dθ2)(θc), further transformation of formula 5 shows that it is hard to evaluate exactly because of its requirement for the exact Hessian matrix of homotopy 3, which is often not desirable. In this work, the second order derivative information is approximated by finite difference, as shown in the following equation: ⎡ dx ⎤ d2x dx (θ ) ≈ B(θc) = ⎢ (θc + δ) − (θc)⎥ /δ 2 c ⎣ ⎦ dθ dθ dθ

∇x h(x ̅ , θc + δ) = ∇x h(x(θc + δ), θc + δ) +

∂ 2h (x(θc + δ), θc + δ) ∂x 2

(x ̅ − x(θc + δ)) + o( x ̅ − x(θc + δ) )

(8)

(12)

where δ is a small positive value. According to eq 6, to be able to use formula 8 to approximate (d2x/dθ2)(θc), x(θc + δ) should be obtained first. x(θc + δ) could be approximated by the following formula: dx x(θc + δ) ≈ x ̅ = x(θc) + (θc) δ dθ

dx 1 (θc) (θt − θc) + B̅ (θc) (θt − θc)2 dθ 2

Hence, according to Lemma 1, there exists a constant ω > 0, when ∥x̅ − x(θc + δ)∥ < ω; the following formulation is satisfied. ∇x h(x ̅ , θc + δ) − ∇x h(x(θc + δ), θc + δ) < (M + 1) x ̅ − x(θc + δ)

(9)

2

Considering that d x(θ)/dθ is continuous in the closed interval [0,1], then d2x(θ)/dθ2 is bound by the positive constant σ. Formulas 4 and 9 show that

Hence, we have the following formula: B(θc) ≈ B̅ (θc) ⎡ ⎤ dx = ⎢ −(∇x h(x ̅ , θc + δ))−1 ∇θ h(x ̅ , θc + δ) − (θc)⎥ ⎣ ⎦ dθ /δ

x ̅ − x(θc + δ) < σδ 2

(14)

According to assumption A3, ∇xh(x(θc + δ), θc + δ) is nonsingular, and the inverse matrix is bound.

(10)

[∇x h(x(θc + δ), θc + δ)]−1 < τ

x(θ ) ≈ x(̃ θ ) = x(θc) +

(13)

2

(15)

Combining formulas 13−15, we obtain the following formula:

dx 1 (θc) (θ − θc) + B̅ (θc) (θ − θc)2 dθ 2

K = [∇x h(x(θc + δ), θc + δ)]−1 [∇x h(x ̅ , θc + δ)

(11)

Based on formulas 6, 9, and 10, the SOC-based backtracking strategy can be described as follows. SOC-HBM. SOC-HBM includes seven steps. Step 1. Set θc = 0, θt = 1, and x(θc) = x(0) = x0 (given initial guess). Step 2. Calculate (dx/dθ)(θc) and B̅ (θc) using formulas 6, 9, and 10. Step 3. Calculate

− ∇x h(x(θc + δ), θc + δ)] < (M + 1)τσδ 2

(16)

Let ⎧ ω , δ0 = min⎨ ⎩ σ 15082

⎫ ρ ⎬ (M + 1)τσ ⎭



Article

When δ < δ0, we have ∥K∥ < ρ < 1. According to Lemma 2, ∥I − K∥ is nonsingular; i.e., [∇xh(x(θc + δ), θc + δ)]−1∇xh(x,̅ θc + δ) is nonsingular. Hence, ∇xh(x,̅ θc + δ) is nonsingular. When δ < δ0, ∥K∥ < ρ < 1 is satisfied, and we obtain the following formulas:

B(θc) − B̅ (θc) = [(∇x h(x(θc + δ), θc + δ))−1∇θ h(x(θc + δ), θc + δ) − (∇x h(x ̅ , θc + δ))−1∇θ h(x ̅ , θc + δ)] /δ

n n

lim (I − K )∑ K = lim (I − K

n →∞

n+1

n →∞

i=0

)=I

= [(∇x h(x(θc + δ), θc + δ))−1∇θ

(17a)

h(x(θc + δ), θc + δ)

(I − K )−1 = [∇x h(x ̅ , θc + δ)]−1 ∇x h(x(θc + δ), θc + δ)

− (∇x h(x ̅ , θc + δ))−1∇θ

∞

=

∑ Kn


(17b)

i=0

+ (∇x h(x ̅ , θc + δ))−1∇θ

∞

[∇x h(x ̅ , θc + δ)]−1 ≤ (∑ K n )


i=0 −1

[∇x h(x(θc + δ), θc + δ)]

− (∇x h(x ̅ , θc + δ))−1∇θ h(x ̅ , θc + δ)]

τ

τ < < 1− K 1−ρ

/δ ⎤ ⎡ τ 2σ(M + 1) 0, when 0 < δ < δ0, we have the following formula:

From formula 8, we obtain the following formula: d2x (θc) − B(θc) = O(δ) dθ 2

d2x (θc) − B̅ (θc) = O(δ) dθ 2

= (∇x h(x(θc + δ), θc + δ))−1[∇x h(x ̅ , θc + δ)

(19)

According to assumption A2, ∇θh(x, θc + δ) can be also expanded at (x(θc + δ), θc + δ). ∇θ h(x ̅ , θc + δ) = ∇θ h(x(θc + δ), θc + δ) +

∂ 2h (x(θc + δ), θc + δ) ∂θ ∂x

(x ̅ − x(θc + δ)) + o( x ̅ − x(θc + δ) )

3. NUMERICAL RESULTS In this study, the performance of SOC-HBM was tested with simulations of production of secreted protein in a fed-batch reactor, a penicillin fermentation process, and a binary distillation column system. Comparisons among SOC-HBM,

(20)

Hence, according to Lemma 1, there exists constant υ > 0, when ∥x̅ − x(θc + δ)∥ < υ; the following formula is satisfied: ∇θ h(x ̅ , θc + δ) − ∇θ h(x(θc + δ), θc + δ) < (M + 1) x ̅ − x(θc + δ)

(21)

Table 1. Comparisons among HBM, S-HBM, and SOCHBM for Simulating Production of Secreted Protein in a Fed-Batch Reactor

Combining formula 21 with formula 14, when δ < (υ/σ) , we obtain the following formula: 1/2

∇θ h(x ̅ , θc + δ) − ∇θ h(x(θc + δ), θc + δ) < (M + 1)σδ

2

(25)

Hence, ∥x(θ) − x̃(θ)∥ = O(δ)|θ − θc|2 + O(|θ − θc|3), and the desired result is obtained. ■ From formulas 24 and 25, we can see that the approximation B̅ (θc) has the same effect as B(θc). For a given homotopy parameter θg, consider a neighborhood N(x(θg),ε), 0 < ε < 1, in which the Newton-based method can solve h(x,θg) = 0 successfully with the starting point x0 ∈ N(x(θg),ε). According to Theorem 2, for SOCHBM, if δ reaches the value of 0, the largest allowable distance between θc and θg can tend to O(ε1/3), and the Newton-based method can solve h(x,θg) = 0 successfully with x0 = x(θc). By contrast, for HBM, the largest allowable distance between θc and θg is O(ε), and the Newton-based method can succeed with x0 = x(θc).

− ∇x h(x(θc + δ), θc + δ)](∇x h(x ̅ , θc + δ))−1 τ < τ(M + 1)σδ 2 1−ρ τ 2σ(M + 1) 2 δ 1−ρ

(24)

Combining formulas 23 and 24, the following formulation is satisfied:

(∇x h(x(θc + δ), θc + δ))−1 − (∇x h(x ̅ , θc + δ))−1

=

(23)

(22)

In the following discussion, the error between B(θc) and B̅ (θc) is considered. 15083

method

sequence of successful homotopy params

total no. iterations

time cost (s)

flag

HBM S-HBM SOC-HBM

0, 0.5, 1.0 0, 1.0 0, 1.0

1048 34 28

134.66 4.58 3.86

successful successful successful



Article

Table 2. Euclidean Distances among x(θ), x(̃ θ), and x(θc) for Production of Secreted Protein in a Fed-Batch Reactora θc

norm of B̅ (θc) − 1.8956 1.7990 1.6985 1.6102 1.5269 1.4453 1.3742 1.3027 1.2373 1.1771

1.00 0.99 0.98 0.97 0.96 0.95 0.94 0.93 0.92 0.91 0.90

× × × × × × × × × ×

|θ − θc|

∥x(θ) − x(θc)∥

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

103 103 103 103 103 103 103 103 103 103

0 2.4066 4.7401 7.0009 9.1903 1.1310 1.3362 1.5349 1.7272 1.9136 2.0941

× × × × × ×

normlzn of ∥x(θ) − x(θc)∥ 0 1.1492 2.2635 3.3432 4.3887 5.4009 6.3808 7.3296 8.2479 9.1381 1.0000

101 101 101 101 101 101

× × × × × × × × × ×

10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−1

∥x(θ) − x̃(θ)∥ 0 1.3169 9.4573 2.8265 6.0101 1.0247 1.5537 2.1761 2.8474 3.5879 4.3680

× × × ×

10−2 10−2 10−1 10−1

normlzn of ∥x(θ) − x̃(θ)∥ 0 3.0149 2.1651 6.4709 1.3759 2.3459 3.5570 4.9819 6.5188 8.2141 1.0000

× × × × × × × × × ×

10−6 10−5 10−5 10−4 10−4 10−4 10−4 10−4 10−4 10−3

a In order to facilitate comparison, ∥x(θ) − x(θc)∥ is normalized by dividing each element with 2.0941 × 102 and ∥x(θ) − x̃(θ)∥ is normalized by dividing each element with 4.3680 × 103.

Figure 1. Profiles of normalization of ∥x(θ) − x(θc)∥ and cube root of the normalization of ∥x(θ) − x̃(θ)∥ for production of secreted protein in a fed-batch reactor.

3.1. Production of Secreted Protein in a Fed-Batch Reactor. The production of secreted protein in a fed-batch reactor is considered as originally posed by Park and Ramirez.15 Here the feed flow rate q (L/h) in the reactor is the only control action. The host-cell growth, gene expression, and the secretion of expressed polypeptides can be described by the following set of equations.15 q PṀ = A(S)(PT − PM) − PM , PM(t0) = PM0 (26a) V q PṪ = B(S)X − PT , PT(t0) = PT0 (26b) V q Ẋ = C(S)X − X , X(t0) = X 0 (26c) V q S ̇ = −YC(S)X + (m − S), S(t0) = S0 (26d) V

Table 3. Model Parameters and Initial Conditions name

value

unit

name

value

unit

name

value

unit

μm Km Ki Yx

0.02 0.05 5 0.5

L/h g/L g/L −

Yp v Sin X0

1.2 0.004 200 1

− L/h g/L g/L

S0 P0 V0 tf

0.5 0 150 150

g/L g/L L h

Table 4. Comparisons among HBM, S-HBM, and SOCHBM for Simulating Penicillin Fermentation Process method



time cost (s)

flag

HBM S-HBM SOC-HBM

0, 0.25, 1.0 0, 1.0 0, 1.0

2182 65 18

348.04 11.13 3.22


S-HBM, and HBM were also presented. Newton’s function is utilized as an auxiliary function in formula 3, and the solver fsolve is used as the Newton-based method for SOC-HBM, SHBM, and HBM. The small positive constant δ for finite difference is set as 10−4. The backtracking constant, λ, is set as 2, and the predefined positive value η is set as 2−8. The maximum number of iterations and the maximum number of function evaluations are set as 1000 and 5000, respectively. The following numerical experiments were carried out on a SONY VGN-SR58F running Windows 7 with Intel Core 2 Duo CPU P8800, 2.66 GHz, and 4.0 GB RAM.

V̇ = q ,

V (t0) = V0

(26e)

with 4.75C(S) S e−5.0S , B (S ) = , 0.12 + C(S) 0.1 + S 21.87S C(S) = (S + 0.4)(S + 62.5)

A (S ) =

(27)

where PM (arbitrary units/L) is the concentration of the secreted protein, PT (arbitrary units/L) is the concentration of 15084



Article

Table 5. Euclidean Distances among x(θ), x(̃ θ), and x(θc) for Penicillin Fermentation Processa θc 1.00 0.99 0.98 0.97 0.96 0.95 0.94 0.93 0.92 0.91 0.90

norm of B̅ (θc) − 3.2436 3.2863 3.3294 3.3726 3.4161 3.4600 3.5040 3.5484 3.5930 3.6380

× × × × × × × × × ×

|θ − θc|

∥x(θ) − x(θc)∥

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

102 102 102 102 102 102 102 102 102 102

0 2.3238 4.6455 6.9652 9.2828 1.1598 1.3912 1.6223 1.8531 2.0838 2.3143

× × × × × × × × × ×

normlzn of ∥x(θ) − x(θc)∥ 0 1.0041 2.0073 3.0096 4.0111 5.0115 6.0113 7.0099 8.0072 9.0040 1.0000

101 101 101 101 102 102 102 102 102 102

× × × × × × × × × ×

∥x(θ) − x̃(θ)∥ 0 7.1482 5.7788 1.9617 4.6722 9.1661 1.5908 2.5372 3.8042 5.4410 7.4982

10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−1

× × × × × × × × × ×

normlzn of ∥x(θ) − x̃(θ)∥ 0 9.5332 7.7069 2.6162 6.2311 1.2224 2.1216 3.3837 5.0735 7.2564 1.0000

10−5 10−4 10−3 10−3 10−3 10−2 10−2 10−2 10−2 10−2

× × × × × × × × × ×

10−7 10−6 10−5 10−5 10−4 10−4 10−4 10−4 10−4 10−3

a In order to facilitate comparison, ∥x(θ) − x(θc)∥ is normalized by dividing each element with 2.3143 × 103 and ∥x(θ) − x̃(θ)∥ is normalized by dividing each element with 7.4982 × 101.

Figure 2. Profiles of normalization of ∥x(θ) − x(θc)∥ and cube root of the normalization of ∥x(θ) − x̃(θ)∥ for penicillin fermentation process.

the secretion rate with respect to S, B(S) is the function form of the expression rate with respect to S, and C(S) is the function form of the specific growth rate with respect to S. The initial conditions of state variables and other system parameters are shown as follows.

Table 6. Comparisons among HBM, S-HBM, and SOCHBM for Simulating Binary Distillation Column method



time cost (s)

flag

HBM S-HBM SOC-HBM

0, 0.5, 0.75, 0.875, 1.0 0, 0.5, 0.75, 1.0 0, 0.5, 1.0

93 59 44

138.40 88.26 68.64


PM0 = 0.0, V0 = 1.0,

the total protein, X (g/L) is the culture cell density, S (g/L) is the culture glucose level, V (L) is the culture volume, Y (=7.3) is the yield of glucose (g)/cell mass (g), m (g/L) is the glucose concentration of the feed stream, A(S) is the function form of

PT0 = 0.0,

X 0 = 1.0,

S0 = 5.0,

m = 20.0

The goal is to evaluate the production of secreted protein in a fed-batch reactor after 10 h (tf = 10) with the feed flow rate set to 0.5 L/h. A simultaneous approach is utilized to simulate

Table 7. Euclidean Distances among x(θ), x(̃ θ), and x(θc) for Binary Distillation Columna θc 1.00 0.99 0.98 0.97 0.96 0.95 0.94 0.93 0.92 0.91 0.90

norm of B̅ (θc) − 3.4749 1.7155 1.0228 6.8515 4.9479 3.7640 2.9743 2.4194 2.0133 1.7065

× × × × × × × × × ×

106 106 106 105 105 105 105 105 105 105

|θ − θc| 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

∥x(θ) − x(θc)∥ 0 1.2605 2.0319 2.6258 3.1153 3.5354 3.9056 4.2379 4.5404 4.8186 5.0766

× × × × × × × × × ×

103 103 103 103 103 103 103 103 103 103

normlzn of ∥x(θ) − x(θc)∥ 0 2.4830 4.0025 5.1724 6.1366 6.9641 7.6933 8.3479 8.9438 9.4918 1.0000

× × × × × × × × × ×

10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−1

∥x(θ) − x̃(θ)∥ 0 2.1584 3.7679 5.7592 7.7346 9.6056 1.1359 1.3000 1.4541 1.5993 1.7365

× × × × × × × × × ×

102 102 102 102 102 103 103 103 103 103

normlzn of ∥x(θ) − x̃(θ)∥ 0 1.2430 2.1698 3.3166 4.4541 5.5316 6.5413 7.4863 8.3737 9.2099 1.0000

× × × × × × × × × ×

10−4 10−4 10−4 10−4 10−4 10−4 10−4 10−4 10−4 10−3

In order to facilitate comparison, ∥x(θ) − x(θc)∥ is normalized by dividing each element with 5.0766 × 104 and ∥x(θ) − x̃(θ)∥ is normalized by dividing each element with 1.7365 × 106. a

15085



Article

Figure 3. Profiles of normalization of ∥x(θ) − x(θc)∥ and cube root of the normalization of ∥x(θ) − x̃(θ)∥ for binary distillation column.

this problem.16,17 The time domain is equally divided into 50 finite elements. The dynamic system of eqs 26a−26e is discretized with Radau collocation (three collocation points) in each finite element. There are 1445 variables in the resulting nonlinear algebraic equations. The initial guess is obtained in this way that all the discretized variables are set to their corresponding initial conditions. HBM, S-HBM, and SOC-HBM were used to simulate the resulting nonlinear algebraic equations, respectively. The results are shown in Table 1. From Table 1, we can see that HBM, S-HBM, and SOCHBM are all successful. SOC-HBM needs the least CPU time cost. In addition, the computational cost for formulas 6, 9, 10, and 11 in SOC-HBM is 0.12 s, which accounts for just 3.11% of the total time cost. According to Theorem 2, ∥x(θ) − x̃(θ)∥ ≈ O(|θ − θc|3) and ∥x(θ) − x(θc)∥ = O(|θ − θc|) are satisfied if θ is sufficiently close to θc, i.e., ∥x(θ) − x̃(θ)∥ ≈ O(∥x(θ) − x(θc)∥3). In order to measure the satisfaction of the relationship in the neighborhood N(θ = 1, ε = 0.1), here the Euclidean distances among x(θ), x̃(θ), and x(θc) are evaluated. The results are shown in Table 2, The cube root of the normalization of ∥x(θ) − x̃(θ)∥ is shown as a dashed line and the normalization of ∥x(θ) − x(θc)∥ is shown as a solid line in Figure 1. From Figure 1, it can be seen that ∥x(θ) − x̃(θ)∥ ≈ O(∥x(θ) − x(θc)∥3) is nearly satisfied. The satisfaction of the relationship can be measured by the error between the cube root of the normalization of ∥x(θ) − x̃(θ)∥ and the normalization of ∥x(θ) − x(θc)∥, which can be indicated with the area between the solid line and the dashed line. It is 4.9833 × 10−4 in this case. The reason for the appearance of the area is that the absolute average rate of change of B̅ (θc) with respect to |θ − θc| is not small (7.9833 × 103), which can be calculated from the second column of Table 2. 3.2. Penicillin Fermentation Process. The penicillin fermentation process is considered next as posed by Srinivasan et al.18 Here the feed rate of substrate u (L/h) is the only control action. The process can be described by the following set of equations. u Ẋ = μ(S)X − X , X(0) = X 0 (28a) V Ṡ = −

μ(S)X vX u − + (Sin − S), Yx Yp V

P ̇ = vX −

V̇ = u ,

u P, V

P(0) = P0

V (0) = V0

(28c) (28d)

with μ(S) = μm S /(K m + S + S2/K i)

(29)

where S (g/L) is the concentration of the substrate; X (g/L) is the concentration of the biomass; P (g/L) is the concentration of the product; V (L) is the volume; Sin (g/L) is the inlet substrate concentration; μm, Km, Ki, and v are kinetic parameters; and Yx and Yp are yield coefficients. The model parameters and initial conditions are given in Table 3. The penicillin fermentation process is simulated in the time domain [0, 150] (h) with the feed rate of substrate set to 1.0 L/ h. The time domain is equally divided into 100 finite elements. The dynamic system of eqs 28a−28d is discretized with Radau collocation (three collocation points) in each finite element. There are 1896 variables in the resulting nonlinear algebraic equations. The initial guess is obtained the same way as in section 3.1. HBM, S-HBM, and SOC-HBM were used to simulate the resulting nonlinear algebraic equations, respectively. The results are shown in Table 4. From Table 4, we can see that HBM, S-HBM, and SOCHBM are all successful. SOC-HBM needs the least CPU time cost. In addition, the computational cost for formulas 6, 9, 10, and 11 in SOC-HBM is 0.15 s, which accounts for just 4.66% of the total time cost. The Euclidean distances among x(θ), x̃(θ), and x(θc) are evaluated in the neighborhood N(θ = 1, ε = 0.1); the results are shown in Table 5. The cube root of the normalization of ∥x(θ) − x̃(θ)∥ is shown as a dashed line and the normalization of ∥x(θ) − x(θc)∥ is shown as a solid line in Figure 2. From Figure 2, it can be seen that ∥x(θ) − x̃(θ)∥ ≈ O(∥x(θ) − x(θc)∥3) is well satisfied. The error between the cube root of the normalization of ∥x(θ) − x̃(θ)∥ and the normalization of ∥x(θ) − x(θc)∥ is 3.2603 × 10−5 in this case. The absolute average rate of change of B̅ (θc) with respect to |θ−θc| is 4.3822 × 102, which is much less than that for the production of secreted protein in a fed-batch reactor. 3.3. Binary Distillation Column System. The binary distillation column separates components cyclohexane and heptane. The dynamic model is based on the Raoult’s law

S(0) = S0 (28b) 15086



Article

resulting nonlinear algebraic equations, respectively. The results are shown in Table 6. From Table 6, we can see that HBM, S-HBM, and SOCHBM are all successful. SOC-HBM needs the least CPU time cost. In addition, the computational cost for formulas 6, 9, 10, and 11 in SOC-HBM is 4.01 s, which accounts for just 5.84% of the total time cost. The Euclidean distances among x(θ), x̃(θ), and x(θc) are evaluated in the neighborhood N(θ = 1, ε = 0.1); the results are shown in Table 7, The cube root of the normalization of ∥x(θ) − x̃(θ)∥ is shown as a dashed line and the normalization of ∥x(θ) − x(θc)∥ is shown as a solid line in Figure 3. From Figure 3, it can be seen that the error between the cube root of the normalization of ∥x(θ) − x̃(θ)∥ and the normalization of ∥x(θ) − x(θc)∥ is the largest compared with the first two cases, which is 1.1448 × 10−3. The reason is that the absolute average rate of change of B̅ (θc) with respect to |θ − θc| is 3.6714 × 107, which is also the largest.

assumption (constant relative volatility), with vapor−liquid equilibrium described by the Wilson equation. The reflux ratio r is the only control variable in the model. This problem was originally described by Hahn and Edgar19 and is paraphrased below. x1̇ =

1 V (y2 − x1), Mcond

xj̇ =

1 L(xj − 1 − xj) − V (yj − yj + 1), M tray

x1(0) = x1,0

(30a)

xj(0) = xj ,0 ,

j = 2, ..., 16

(30b)

1 (Fxf + Lx16) − (F + L)x17 − V (y17 − y18 ), M tray

x17 ̇ =

x17(0) = x17,0 xj̇ =

(30c)

1 (F + L)(xj − 1 − xj) − V (yj − yj + 1), M tray

xj(0) = xj ,0 , x32 ̇ =

j = 18, ..., 31

4. CONCLUSION In this study, the SOC-HBM strategy for process simulation is proposed for improving HBM by introducing second-order correction. The theoretical foundation for SOC-HBM is established based on some assumptions. Simulation runs on the production of secreted protein in a fed-batch reactor, a penicillin fermentation process, and a binary distillation column system were used as test problems. The numerical results show that SOC-HBM exhibits the best performance compared with HBM and S-HBM. In addition, the computational cost of the second-order correction in SOC-HBM accounts for a few percent of the total time cost. They are 3.11, 4.66, and 5.84% respectively in the test problems. It is much less than the computational cost saved by reducing the number of iterations in the solving procedure. That is why SOC-HBM is feasible. The theoretical analysis shows that ∥x(θ) − x̃(θ)∥ ≈ O(∥x(θ) − x(θc)∥3) is satisfied if θ is sufficiently close to θc. In the numerical tests, the satisfaction is indicated by the error between the cube root of the normalization of ∥x(θ) − x̃(θ)∥ and the normalization of ∥x(θ) − x(θc)∥ in the neighborhood N(θ = 1, ε = 0.1). They are 4.9833 × 10−4, 3.2603 × 10−5, and 1.1448 × 10−3 respectively in the test problems. The satisfaction is best in the penicillin fermentation process simulation, while it is worst in the binary distillation column simulation, whose absolute average rate of change of B̅ (θc) with respect to |θ − θc| is 3.6714 × 107. This implies that the fourth item of the Taylor expansion plays an important role. The only method to improve the satisfaction is to reduce the radius of the neighborhood.

(30d)

1 (F + L)x31 − (F − D)x32 − Vy32 , M reb

x32(0) = x32,0

(30e)

sat sat P = xjγA, jPA, j + (1 − xj)γB, jPB, j ,

sat yj P = xjγA, jPA, j,

j = 1, ..., 32 sat sat PA, j = PA (Tj),

(30f)

PB,satj = PBsat(Tj),

j = 1, ..., 32 (30g)

γA, j = γA(xj , Tj),

γB, j = γB(xj , Tj),

j = 1, ..., 32 (30h)

D = 0.5F ,

L = rD ,

V=L+D

(30i)

where x1 and y1 are the reflux drum liquid and vapor mole fractions of cyclohexane; xj and yj, j = 2, ..., 31, are the liquid and vapor mole fractions of cyclohexane on tray j − 1; x32 and y32 are the reboiler liquid mole fractions of cyclohexane; Tj, j = 1, ..., 32 (K), is the temperature at each stage; F (mol/min) is the feed flow rate; xf is the mole fraction of feed; D is the distillate flow rate (mol/min); L (mol/min) is the flow rate of the liquid in the rectification section; V (mol/min) is the vapor flow rate in the column; Mtray is the total molar holdup on each tray; Mcond is the total molar holdup in the condenser; Mreb is the total molar holdup in the reboiler; P (Pa) is the total pressure in column; γA,j and γB,j, j = 1, ..., 32, refer to activity sat coefficients defined by the Wilson equation; and Psat A,j and PB,j, j = 1, ..., 32, are the vapor pressures. Additional problem data and sat expressions for γA,j, γB,j, Psat A,j, and PB,j, j = 1, ..., 32, can be found in http://www.hedengren.net/research/models.htm. The binary distillation column is simulated in the time domain [0, 240] (min) with the reflux ratio set to 10. The time domain is equally divided into 20 finite elements. The dynamic system of eqs 30a−30i is discretized with Radau collocation (three collocation points) in each finite element. There are 6368 variables in the resulting nonlinear algebraic equations. The initial guess is obtained the same way as in section 3.1. HBM, S-HBM, and SOC-HBM were used to simulate the

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the 973 Program (No. 2009CB320603) and National Natural Science Foundation of China (No. 21206149, No. 61203132). 15087



■

Article

REFERENCES

(1) Mousa, A. A.; El-Desoky, I. M. GENLS: Co-evolutionary Algorithm for Nonlinear System of Equations. Appl. Math. Comput. 2008, 197 (2), 633. (2) Hirsch, M. J.; Pardalos, P. M.; Resende, M. G. C. Solving Systems of Nonlinear Equations with Continuous GRASP. Nonlinear Anal.: Real World Appl. 2009, 10 (4), 2000. (3) Ren, H.; Wu, L.; Bi, W.; Argyros, I. K. Solving Nonlinear Equations System Via an Efficient Genetic Algorithm with Symmetric and Harmonious Individuals. Appl. Math. Comput. 2013, 219 (23), 10967. (4) Dennis, J. E., Jr.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1996. (5) Zhang, J.; Wang, Y. A New Trust Region Method for Nonlinear Equations. Math. Methods Oper. Res. 2003, 58 (3), 283. (6) Fan, J.; Pan, J. An Improved Trust Region Algorithm for Nonlinear Equations. Comput. Optim. Appl. 2011, 48 (1), 59. (7) Nocedal, J.; Wright, S. Numerical Optimization; Springer: New York, 1999. (8) Billups, S. C.; Watson, L. T. A Probability-one Homotopy Algorithm for Nonsmooth Equations and Mixed Complementarity Problems. SIAM J. Optim. 2002, 12 (3), 606. (9) Jalali, F.; Seader, J. D.; Khaleghi, S. Global Solution Approaches in Equilibrium and Stability Analysis Using Homotopy Continuation in the Complex Domain. Comput. Chem. Eng. 2008, 32 (10), 2333. (10) Malinen, I.; Tanskanen, J. Homotopy Parameter Bounding in Increasing the Robustness of Homotopy Continuation Methods in Multiplicity Studies. Comput. Chem. Eng. 2010, 34 (11), 1761. (11) Zhu, L.; Chen, Z.; Chen, X.; Shao, Z.; Qian, J. Simulation and Optimization of Cryogenic Air Separation Units Using a Homotopybased Backtracking Method. Sep. Purif. Technol. 2009, 67 (3), 262. (12) Chen, W.; Zhu, L.; Chen, X.; Shao, Z. Sensitivity Embedded Homotopy-based Backtracking Method for Chemical Process Simulation. 10th IEEE International Conference on Control and Automation; IEEE Press: Piscataway, NJ, 2013. (13) Wayburn, T. L.; Seader, J. D. Homotopy Continuation Methods for Computer-Aided Process Design. Comput. Chem. Eng. 1987, 11 (1), 7. (14) Horn, R. A.; Johnson, C. R. Matrix Analysis; Cambridge University Press: London, 1985. (15) Park, S.; Ramirez, W. R. Optimal Production of Secreted Protein in Fed-Batch Reactors. AIChE J. 1988, 34 (9), 1550. (16) Biegler, L. T.; Cervantes, A. M.; Wachter, A. Advances in Simultaneous Strategies for Dynamic Process Optimization. Chem. Eng. Sci. 2002, 57 (4), 575. (17) Kameswaran, S.; Biegler, L. T. Simultaneous Dynamic Optimization Strategies: Recent Advances and Challenges. Comput. Chem. Eng. 2006, 30 (10−12), 1560. (18) Srinivasan, B.; Bonvin, D.; Visser, E.; Palanki, S. Dynamic Optimization of Batch Process II. Roles of Measurements in Handling Uncertainty. Comput. Chem. Eng. 2003, 27 (1), 27. (19) Hahn, J.; Edgar, T. F. An Improved Method for Nonlinear Model Reduction Using Balancing of Empirical Gramians. Comput. Chem. Eng. 2002, 26 (10), 1379.

15088


Homotopy with Second-Order Correction Based Backtracking Method

Recommend Documents