Optimal Grade Transition for Polyethylene Reactors Based on

Article pubs.acs.org/IECR

Optimal Grade Transition for Polyethylene Reactors Based on Simultaneous Strategies and Trust Region Particle Swarm Optimization Wenxing Xu,†,|| Zhiqiang Geng,† Qunxiong Zhu,*,† and Xiangbai Gu†,‡ †

College of Information Science and Technology, Beijing University of Chemical Technology, Beijing, 100029 Sinopec Engineering, Beijing, 100029

‡

ABSTRACT: Various grades of polyethylene are required to be produced in the same industrial equipment by changing the operating conditions. However, transitions between different grades are rather slow and result in the production of a considerable amount of off-specification polymer. To improve grade transition is viewed here as a dynamic optimization problem, for which simultaneous strategies are especially beneficial. This paper first used a simultaneous strategy to fully discretize the state and control variables, converting the mathematical optimization model into a large-scale NLP problem. Second, a hybrid particle swarm optimization (PSO) integrated with trust region method (TRPSO) was proposed. The simulation results on some multimodal global optimizations showed that the algorithm is far more effective than linear decreasing weight PSO (LDWPSO) and sequential quadratic programming (SQP) to search the global optimum. Later, simultaneous strategies and the TRPSO were combined to solve the slurry high density polyethylene grade transition in a real plant. The detailed comparisons were derived from TRPSO, LDWPSO, iteration dynamic programming (IDP), and the real historical transition data of a plant demonstrated the effectiveness and practicability of the proposed method.

1. INTRODUCTION Polyethylene is widely used today in a multitude of products, each of which requires specific polymers with some designated specifications. To meet specific needs and reduce cost, polymer companies often have to produce various grades in the same equipment by changing the operating conditions. Normally, a considerable amount of off-specification polymer is produced during grade transition processes,1 so the grade transition optimization problems are of great interest. This problem has been studied by a number of researchers since the mid-1990s. McAuley and MacGregor,2 for the first time, proposed optimal grade transition strategies in a gas phase polyethylene reactor by solving a dynamic model-based optimization problem. Debling et al.1 tested different grade transition operations using the simulation package POLYRED. Takeda and Ray3 and Wang et al.4 applied the control vector parametrization (CVP) method to calculate optimal grade transition strategies for a slurry-phase loop reactor and a gas-phase fluidized-bed reactor, respectively. BenAmor et al.5 applied an industrial real-time optimization package (ROMeo) to the nonlinear model predictive control (NLMPC) of a simulated polymer grade transitions. The NLMPC algorithm was formulated by combining orthogonal collocation and sequential quadratic programming. Bonvin et al.6 used a measurement-based approach to track the necessary conditions of optimality for grade transition. These methods are relatively easy to construct and to apply and have indeed improved process operation, but they may become computationally expensive for large scale problems due to repeated numerical integration of the DAE model. Flores-Tlacuahuac and Biegler et al.7−11 cast this type of problem as a mixed-integer dynamic optimization (MIDO) formulation, exploiting a structure with simultaneous approaches to solve them. Wang and Yang et al.12,13 introduced sequential © 2013 American Chemical Society

quadratic programming (SQP) and iteration dynamic programming (IDP) into grade transition optimization of fluidized bed gas-phase and slurry polyethylene to improve industrial operations. Lee et al.14 established a simple melt index controller with only proportional-integral mode based on an artificial neural networks model for grade transition and significantly reduced the grade transition time in comparison with the base strategy by step-changing of the operating recipe in simulation. Generally, either target function methods or dynamic simulation software have been used to simulate and optimize the polymerization process of ethylene and obtain the optimized trajectory of operation variables and polymer properties. However, compared with dynamic simulation, it is more likely to achieve more accurate optimization by solving objective function models. Furthermore, in the target function method for grade transition optimization, a dynamic optimization problem is to be solved. For this kind of problem, there are three basic types of solving methods: the variational approach, sequential strategy, and simultaneous strategy. As proved by FloresTlacuahuac and Biegler et al.,9−11 a better way to solve the dynamic optimization model of grade transition is to convert the model using a simultaneous strategy into a NLP problem and then to use software or algorithms specific for solving NLP problems. Simultaneous approaches usually lead to fast solution of dynamic optimization problems and can be applied to handle open loop unstable systems and path constraints. Until now, the algorithms used for the derived large scale NLP after changing from the original dynamic optimization problems are Received: Revised: Accepted: Published: 3363

March 16, 2012 December 24, 2012 January 22, 2013 January 22, 2013 dx.doi.org/10.1021/ie300712e | Ind. Eng. Chem. Res. 2013, 52, 3363−3372

Industrial & Engineering Chemistry Research

Article

Figure 1. Process flow diagram of slurry HDPE polymerization process.

material, hexane as solvent, butene as comonomer, and hydrogen as chain terminator. The ethylene is passed into in gas phase, dissolved in the solvent, and spread to the surface of the catalyst particles to react, which makes a gas−liquid−solid threephase reaction. In the real production process, the mass transfer between phases is strengthened by selecting the stirring speed and gas entrance location. The post-treatment section includes separation, granulation, and other processes, in which the stabilizer is to be added to granulation process and has a comparatively great impact on the polymer properties. 2.2. Process Modeling. For many chemical engineering problems, process modeling always involves nonlinear differential algebraic equations (DAEs). The general formulation for optimization of a dynamic system can be written as follows:

basically nonlinear optimization methods, such as SQP and IDP. However, the SQP algorithm has a risk to fall into the local optimal solution depending on the selection of initial solution; the IDP needs more iteration steps with the increasing desired accuracy and thus has the disadvantage of slow convergence. With the development of swarm intelligence, heuristic intelligent methods, which are relatively easy to achieve and do eventually locate the desired solution, have also been tried to solve large-scale NLP problems.15−20 In this paper, we try to combine the advantages of PSO and traditional gradient method together for solving the derived large scale NLP in grade transition. Furthermore, correlation models with the variation of melt index (MI), density, temperatures, and gas concentrations in the reactor are utilized to make optimization to be solved easily.2,12,13,21 In this work, the optimal grade transition will be considered as a minimal off-specification products problem. After establishing the objective and constraint functions for grade transition in a slurry polyethylene reactor, the dynamic optimization problem is first converted into a NLP problem by a simultaneous strategy. Then, a novel TRPSO algorithm, which is a hybrid of PSO and trust region method, is proposed and applied to solve this NLP problem. The performance of the proposed method is illustrated and evaluated with an actual slurry high density polyethylene grade transition process. Based on some real industrial processing data, the evaluation and research of the simultaneous strategy and TRPSO algorithm presented in this paper are carried into execution. The results obtained are then discussed, and remarks about the method are also presented. The organization of this paper is as follows. Process description and modeling are introduced in section 2, followed by description of our novel nonlinear optimization approach in section 3 and the application, results, and discussion in section 4. Finally, we draw the conclusions in section 5.

min

z(t ), y(t ), u(t )

s.t.

J = Φ(z(t ), y(t ), u(t ))

dz (t ) = f (z(t ), y(t ), u(t )) dt

t ∈ [t0 , t f ] z(t0) = z 0

(1)

(2)

g (z(t ), y(t ), u(t )) = 0

(3)

u low ≤ u(t ) ≤ u up

(4)

ylow ≤ y(t ) ≤ yup

(5)

z low ≤ z(t ) ≤ zup

(6)

where u(t), y(t), and z(t) are vector functions of control, algebraic, and differential state variables, respectively, J = Φ(z(t),y(t),u(t)) is the objective function, and eqs 2 and 3 and boundary inequalities 4−6 are the constraints. An optimizer is to choose functions of the scalar “time” parameter to minimize the objective function while satisfying the constraint conditions, which include a set of nonlinear DAEs. A simultaneous approach fully discretizes the DAE optimization problem into an NLP by approximating state and control profiles by a family of polynomials on finite elements (t0 < t1 < ··· < tNE = tf). These polynomials can be represented as power series, sums of orthogonal polynomials, or in Lagrange form. Here, we use Radau collocation points because they allow constraints to be set at the end of each element and to stabilize the system more efficiently if high index DAEs are present. Monomial basis representation for the differential profiles is as follows:

2. PROCESS DESCRIPTION AND MODELING 2.1. Process Description. In this study, the grade transition process of industrial high density polyethylene (HDPE) continuous slurry polymerization is studied. The schematic diagram of the reactor system is given in Figure 1.13 The productive technology consists of three major parts: the pretreatment section of raw materials and catalyst, reaction section, and post-treatment section. The reaction section involves two reactors. These two reactors can be switched between in series and in parallel by adjusting the pipeline, while the polymer grades can be controlled by varying the operating parameters and compositions of the polymers in two reactors. During the production process, both temperature and pressure are low. The main domestic catalyst, BCH catalyst, and imported cocatalyst, AT cocatalyst, are added to the reactors in certain flow rates. The reaction process is carried forward with ethylene of high purity as the main raw

Ncp ⎛ t − t l − 1 ⎞ dz z(t ) = zl − 1 + hl ∑ Ωq⎜ ⎟ ⎝ hl ⎠ dtl , q q=1

3364

(7)

dx.doi.org/10.1021/ie300712e | Ind. Eng. Chem. Res. 2013, 52, 3363−3372


Article

u low ≤ ul , q ≤ u up ,

Here zl−1 is the value of the differential variable at the beginning of element l, hl is the length of element l, dz/dtl,q is the value of its first derivative in element l at the collocation point q, and Ωq is a polynomial of order Ncp, satisfying eq 8. ⎧ q = 1, ···, Ncp ⎪ Ωq(0) = 0 ⎨ ⎪ ′ ⎩ Ωq(ρr ) = δq , r q , r = 1, ···, Ncp

(20)

As discussed in the literature, for the grade transition process of slurry high density polyethylene here, the melt index (MI) and density (ρ) in the two reactors and pelletizer are the main polymer properties, which are also known as output variables of the process, and the reaction temperature (T), the ratio of hydrogen to ethylene ([H2]/[C2]) and the ratio of comonomer to ethylene ([C4]/[C2]) are the main manipulating parameters, which are also called input variables of the process. Therefore, vectors of state and control variables are represented as follows:

where ρr = (tl,r − tl−1)/hl is the location of the rth collocation point within each element, δq,r is sign function: ⎪

⎪

(9)

z = (z1 , z 2 , ···, z NZ)T

To satisfy eqs 8 and 9, Ωq is constructed as the integration of Lagrange polynomial: t

Ωq(t ) =

Ncp

∫0 ∏

r = 1, r ≠ q

= (ln MI1,c , ρ1,c , ln MI 2,c , ρ2,c , ln MI p)T

(22)

(10) ⎛ u = (u1, u 2 , ···, uNU) = ⎜⎜T1, ⎝

The continuity of the differential profiles is enforced by Ncp q=1

dz dt l , q

T

⎛ t − tl − 1 ⎞ y(t ) = yl − 1 + hl ∑ Ψq⎜ ⎟y ⎝ hl ⎠ l , q q=1

(12)

Ncp ⎛ t − tl − 1 ⎞ u(t ) = ul − 1 + hl ∑ Ψq⎜ ⎟u l , q ⎝ hl ⎠ q=1

(13)

specifications produced during transition from grade A to grade B, the objective function is converted to eq 24.

Here, yl,q and ul,q represent the vector of the algebraic and control variables, respectively, in element l at collocation point q. Ψq is the Lagrange polynomial of degree Ncp − 1 satisfying

J=

(14)

min

s.t.

J = Φ(zl , q , yl , q , ul , q)

dz = f (zl , q , yl , q , ul , q) dtl,q Ncp

zl , q = zl − 1 + hl

∑ Ωq′(ρq ) q ′= 1

t ∈ [t0 , tf ]

z(t0) = z 0

dz dt l , q

g(zl , q , yl , q , ul , q) = 0

tf

w(x − xsp)2 dt

(24)

NZ + NY Nfe Ncp

J ≈ tf

∑ ∑ ∑ whi l(xi ,l ,p − xisp)2 i=1

(15)

l=1 p=1

(25)

Suppose the hybrid condition of gas−liquid−solid threephase in two reactors is good and ignore the effect of diffusion process on polymerization. The instantaneous quality index can be predicted using the steady state model of objective grade according to the literature.13

(16)

ln MI i = A 0 +

(17)

⎡ ⎡ [H ] ⎤ ⎡ [C ] ⎤⎤ A1 + A 2 ln⎢A3 + A4 ⎢ 2 ⎥ + A5⎢ 4 ⎥⎥ ⎢⎣ T ⎣ [C2] ⎦ ⎣ [C2] ⎦⎥⎦ (26)

(18)

0.5 ⎤ ⎡⎡ ⎡ [H ] ⎤0.5 [C ] ⎤ ρi = B0 + B1⎢ 2 ⎥ + B2 ln MI + B3⎢⎢ 4 ⎥ × ln MI⎥ ⎥⎦ ⎣ [C2] ⎦ ⎣⎢⎣ [C2] ⎦ (27)

Ncp

dz zl = zl − 1 + hl ∑ Ωq(1) dt l , q q=1

∫0

Using the simultaneous strategy described above to fully discretize the DAE system shown above over finite elements, the objective function is further approximately transformed to eq 25.

From eq 7, the differential variables are required to be continuous throughout the time horizon, while the control and algebraic variables are allowed to have discontinuities at the boundaries of the elements. Substitution of eqs 7−14 into eqs 1−6 leads to the following NLP: dz / dl , qyl , q , ul , q

T ⎡ [C4 ] ⎤ ⎞ ⎥⎟ ⎢ ⎣ [C2] ⎦2 ⎟⎠

were MI1,c, ρ1,c, MI1,i, ρ1,i, MI1,sp, ρ1,sp and MI2,c, ρ2,c, MI2,i, ρ2,i, MI2,sp, ρ2,sp are the cumulative, instantaneous, and specified values of polymer properties in reactor one and two. MIp and MIp,sp are the instantaneous and objective melt indexes of products in the pelletizer. t0 is the start time, and tf is the end time of grade transition. NZ is the total number of differential state variables. NY is the total number of algebraic variables; NU is the total number of manipulated variables. ⎡z⎤ Suppose x = ⎢ y ⎥ and our objective is to minimize the off⎣ ⎦

Ncp

q , r = 1, ···, Ncp

⎡ [H ] ⎤ ⎡ [H 2] ⎤ ⎥ , T2 , ⎢ 2 ⎥ , ⎢ ⎣ [C2] ⎦2 ⎣ [C2] ⎦1

(23)

(11)

In addition, the control and algebraic profiles are approximated using a Lagrange basis representation, which takes the form:

Ψq(ρr ) = δq , r

(21)

y = (y1 , y2 , ···, yNY )T = (ln MI1,i , ρ1,i , ln MI 2,i , ρ2,i )T

(z − zr ) dz (zq − zr )

zl = zl − 1 + hl ∑ Ωq(1)

z low ≤ zl , q ≤ zup

13

(8)

⎧ 0, q ≠ r δq , r = ⎨ ⎩1, q = r

ylow ≤ yl , q ≤ yup ,

(19) 3365



Article

3. NONLINEAR OPTIMIZATION 3.1. Trust Region (TR) Method. As one of the traditional numerical methods for solving NLP problem, the trust region method determines the direction and step length simultaneously so as to avoid linear search.22 Considering the unconstrained minimization problem,

So part of the state variable constraints can be represented as: yi , j , q = gĩ (ul , q),

i = 1, ···, NY,

l = 1, ···, Nfe ,

q = 1, ···, Ncp

(28)

where g̃i(ul) are regression functions of eqs 26 and 27 according to historical data. The relationship between instantaneous and cumulative properties at the exit of the reactor is as shown in eq 29, dPc(t ) = Fr[Pi(t ) − Pc(t )] dt

min f (x)

An approximate function within the neighborhood of current iteration point is determined to derive the trust region subproblem min qk (x)

(29)

where k is the iteration number. In a standard trust region algorithm, the quadratic approximation is used. Thus, the trust region subproblem is in the form of eq 38.

d zi = f (zi , l , q , yi , l , q , ul , q) = Fr(yi , l , q − zi , l , q), dtl,q

s.t. (30)

ln MI p = D0 + D1 ln MI1, i + D2 ln MI 2,i + D3 ln MI 2,c + D4 ρ2,c

(31)

The last state variable constraint is represented as dz NZ = f (zl , q , yl , q , ul , q) = D0 + D1y1, l , q + D2y3, l , q dtl,q + D3z 3, l , q + D4 z4, l , q

(32)

Combining eqs 28, 30 and 32, and the industrial operation rules, they are represented as manipulating bounds, so the constraints for all state variables are derived. Therefore, after using a simultaneous strategy to complete the full discretization, the objective function for a slurry polyethylene grade transition is considered here as an NLP problem shown by eqs 25, 30, 32, 17, 28, 19, and 20. Since in real plant, only manipulated variable profile u(t) can be adjusted, and it can be manipulated once during an element, which means ul(t), l = 1, ···, Nfe, stays the same at each collocation point during an element. This NLP can be rewritten as eqs 33−35: min f (x) c(x) = 0

(39)

Presk = q(0) − q(sk)

(40)

Aresk Presk

(41)

Good algorithms exist for solving eq 38. Preconditioned conjugate gradient (PCG) method is one of them, and it restricts the trust-region subproblem to a two-dimensional subspace spanned by the gradient direction g and d, where d is either an approximate Newton direction (i.e., a solution to eq 42) or a direction of negative curvature computed by eq 43. This method overcomes the limitation that the Hessian matrix is guaranteed to be positive definite only in the neighborhood of a strong minimizer and requires trivial memory space for computation. Thus, it is appropriate for a large scale optimization.

(34)

xlow ≤ x ≤ xup

Aresk = f (xk) − f (xk + sk)

rk =

(33)

x ∈ Rfe

(38)

where, s = x − xk is a trial step, gk = ∇f(xk) is the gradient of f(x) at the current point xk, Bk is the Hessian matrix (the symmetric matrix of second derivatives) of f(x), Δk ≥ 0 is the radius of the subregion, ∥·∥ is a certain norm, for which usually 2-norm is adopted. After solving the subproblem, the trial step sk is derived. Aresk and Presk computed by eqs 39 and 40 are defined as the actual and predicted decline amount of the kth step respectively. Then, the ratio of the two amounts rk = Aresk/Presk is applied to determine whether this trial is to be adopted and the way to adjust the value of Δk.

As the MI in pelletizer is modeled as follows:

(35)

Bk d = −g

(42)

dTBk d < 0

(43)

A sketch of box-constrained minimization using the trustregion reflective algorithm is given as follows:

where x = {ul|l = 1, ···,Nfe}, f: R →R, and c: R →R . The initial mathematical optimization model is composed of eqs 21−24, 26, 27, 29, and 31. The results in the literature13 have shown that since the factor of pelletizer was added into the objective function, the sharp change of controlled variable is effectively restrained and the security of the grade transition process is ensured. Thus, the NLP model of eqs 33−35 derived after full discretization may reflect the practical reaction process satisfactorily. fe

1 T s Bk s + gkT s 2 || s || ≤ Δk

min qk (s) =

i = 1, ···, NZ − 1

fe

(37)

x∈N

where P is the properties of polymer, such as ln(MI) and ρ, and Fr is the fraction of polymer produced in the reactor per unit interval. Thus, eq 16 is in the form

s.t.

(36)

x ∈ Rn

m

Step 1: Use a scaled modified Newton step to formulate the two-dimensional trust-region subproblem with the aid of PCG algorithm. Step 2: Solve eq 38 to determine the trial step sk. Step 3: Calculate the value of rk. Step 4: If rk meets the conditions, of which often f(xk + sk) < f(xk) is used, then xk+1 = xk + sk. 3366



Article

Table 1. Benchmark Functions function name

function

Rosenbrock

searching space [−30, 30]n

n−1

f1 (x) =

∑ [100(xi+ 1 − xi2)2 + (xi − 1)2 ] i=1

Rastrigin

[−5.12, 5.12]n

n

f2 (x) =

∑ [xi 2 − 10cos(2πxi) + 10] i=1

Girewank

f3 (x) =

1 4000

n

n i=1

i=1

3.3. Hybrid Trust Region and Particle Swarm Optimization (TRPSO) Algorithm. The TRPSO algorithm takes use of global search ability of PSO and fast convergence of the TR method. In every iteration of PSO, each particle is searching its trust region to get the local optimum of its surrounding region with the help of gradient information. Thus, the whole particle swarm is more likely to search the whole space based on fast local convergence. Not only the accuracy but also efficiency of the algorithm is promoted. Meanwhile, the PCG method, which needs small memory space, is used solving the trust region subproblem. This algorithm is suitable for a large scale optimization. TRPSO algorithm procedure: Step 1: Randomly initialize the particle swarm within each variable’s domain. Set the problem’s dimension N. Empirically set swarm size M and maximum function evaluations Max_Fes. Uniformly, randomly initialize each particle’s position within the box constraint bounds as xi = (xi1, xi2, ···, xiN)T, velocity vi = (vi1, vi2, ···, viN)T. Step 2: Starting from each xi, use the TR algorithm to search a local optimum for initializing the particle’s individual optimum pi and global optimum pg = minpopsize pi. i=1 Set total iteration count j = 1; update FEs according to eq 46:

Step 5: If the convergence is received, the iteration will be terminated. Then, return the result. Otherwise, go to Step 6. Step 6: Adjust Δk and go to Step 1. Herein, box constraints are handled by reflective Newton step and linear equality constraints are dealt with a sparse leastsquares step and reduced preconditioned conjugate gradients, as described in the literature.23 For nonlinear equality constraints, only a large constant penalty factor 1 × 106 is introduced here. Despite a high convergence rate, the TR algorithm should be used very carefully since it is very sensitive to initial points. Generally it is not guaranteed to attain the global optimum. 3.2. PSO Algorithm. Particle Swarm Optimization (PSO) is a swarm intelligence heuristic search technique proposed by Kennedy and Eberhart in 1995.15 In order to emulate the foraging behavior of birds, in PSO, every solution is a “bird” in the search space, which is also called a “particle”. A swarm is initialized with random positions xi = (xi1, xi2, ···, xiN)T and velocities vi = (vi1, vi2, ···, viN)T of each particle. Then, while flying through the problem search space, each particle adjusts its velocity to find a better solution (position) by applying its own flying experience (personal best position found in its earlier flights) pi = (pi1, pi2, ···, piN)T and experience of neighboring particles (global best position found of the population) pg = (pg1, pg2, ···, pgN)T. In standard PSO, the new velocity and position for each particle are updated by the following equations:

M

FEs = FEs +

xidk+ 1 = xidk + vidk+ 1

∑ FEs_TR i i=1

vid k + 1 = wvid k + c1 × rand() × (pid − xid k) + c 2 × rand() × (pgd − xid k)

[−600, 600]n

⎛ xi ⎞ ⎟+1 i⎠

∑ xi 2 − ∏ cos⎜⎝

(46)

where FEs_TRi represents the function evaluation number consumed by the ith particle using TR for the local search in the jth iteration. Step 3: If FEs < Max_Fes, go on; otherwise, terminate and return optimization result f(pg). Step 4: Update particle velocity vi and position xi according to eqs 44 and 45 of all M particles. Step 5: Taking each xi as a starting point, calculate local optimum position Pi′ for each particle within its neighborhood by TR algorithm. j = j + 1. Update FEs according to eq 46. Step 6: Update individual optimum pi for each particle respectively according to eq 47, and calculate the global optimum pg. Turn to Step 3.

(44) (45)

where i = 1, 2, ···, M and M is the swarm size. d = 1, 2, ···, N and N is the dimensionality of the search space; the superscript k indicates the iteration count, rand() is a uniform random value in the range [0,1], scaling parameters c1 and c2 are positive constants for regulating the maximum step sizes for the particles to fly toward pi and pg, respectively. The dimensional speed of each particle is clamped to the range [vdmin, vdmax] to control excessive roaming of particles outside the search space. In the standard PSO, w is a constant during the searching iteration. If the parameter w linearly decreases, as recommended from Shi and Eberhart,16 the algorithm evolves into the so-called Linear Decreasing Weight Particle Swarm Optimization (LDWPSO). As a random search strategy, the PSO algorithm, which is easy to achieve, has not only the advantages of low memory space and high efficiency but also the flexibility of balancing global and local searches. However, although it converges quickly in the beginning of the search, around the optimum its local search slows down quickly.

pi = min{pi , pi′}

(47)

Three commonly used benchmark functions shown in Table 1 are used to test the performance of the algorithm described (TRPSO) with comparison to that of LDWPSO and SQP.24 Algorithm parameters are set the same as commonly used in literature without tuning beforehand, which are as follows: dimension N = 10/50/100; population size M = 20; maximum number of function evaluations Max_FEs = 1000(N + 10); 3367



Article

Table 2. Results Comparison among TRPSO, LDWPSO, and SQP AvgBest function Rosenbrock

Rastrigin

Griewank

Dim 10 50 100 10 50 100 10 50 100

TRPSO 1.55 2.93 1.03 6.65 9.91 2.01 0.00 0.00 0.00

−8

× 10 × 10−15 × 10−24 × × × × ×

101 102 1000 1000 1000

SR

LDWPSO 3.58 7.18 1.61 1.26 1.30 2.26 4.58 7.01 2.39

× × × × × × ×

SQP

2

10 104 105 101 102 102 10−1

1.12 2.09 9.66 7.04 3.20 5.86 3.53 2.21 2.90

× 101

wmax = 0.9, wmin = 0.4; c1 = c2 = 2; the upper and lower bounds for velocities are vmax = (χmax − χmin)/2 and vmin = −vmax. Fifty independent runs have been executed for each test function. Values that are less than 1.00 × 10−100 are recorded as 0.00 × 1000. The simulation results are shown in Table 2, where AvgBest means the average of optimal function values found by the algorithm in all independent runs, VTR represents values to reach, and SR is the success rate. A run of an optimization algorithm is considered to be successful with respect to the VTR if and only if the optimizer has discovered a candidate solution for which the objective function takes on a value below or equal to the VTR before FEs exceeds Max_Fes. If SRN and TRN are the number of successful runs and total number of runs respectively, SR is computed as given in eq 29. SR (%) =

SRN × 100 TRN

× × × × × × × ×

101 101 101 102 102 10−1 10−8 10−8

VTR 1.00 1.00 1.00 2.00 1.20 2.50 1.00 1.00 1.00

× × × × × × × × ×

−6

10 10−6 10−6 101 102 102 10−6 10−6 10−6

TRPSO

LDWPSO

SQP

50/50 50/50 50/50 50/50 50/50 50/50 50/50 50/50 50/50

0/50 0/50 0/50 42/50 19/50 31/50 4/50 3/50 3/50

36/50 9/50 0/50 0/50 0/50 0/50 14/50 50/50 50/50

Figure 2. Logarithmic plots of mean function value versus function evaluations of TRPSO and LDWPSO on Rosenbrock of 50 dimensions.

(29)

Table 2 shows that, for some benchmarks, such as Rosenbrock and Griewank, SQP can get more accurate results with larger success rate that LDWPSO, and for some benchmarks, such as Rastrigin, SQP can hardly get satisfactory results without skillfully setting initial search point. Moreover, the effects of both LDWPSO and SQP are dependent on the dimension of the problem. However, for all tested benchmarks of all dimensions, the average optimal function values found by TRPSO are the closest to the theoretical optima (zero). That is to say that the proposed TRPSO method outperforms the other two methods in terms of accuracy. Since for the VTR predefined TRPSO can derive successful rates of 100% for all tested benchmarks despite of dimensions, which is far better than that derived by LDWPSO and SQP, it can also be concluded that TRPSO has the strongest robustness of the three. Furthermore, comparisons of the algorithmic convergence speed between TRPSO and LDWPSO on the three benchmarks of 50 dimensions are illustrated in Figures 2−4. It is shown that the fitness values decreases much faster in the TRPSO algorithm than in the LDWPSO algorithm. For function Rosenbrock, in the TRPSO algorithm, the fitness values fly to a satisfactory accuracy after less than 200 function evaluations, while in the LDWPSO algorithm the fitness values decrease slowly after 5000 function evaluations and never reach a satisfactory accuracy until the iteration ends. For function Rastrigin, the convergence speed of TRPSO exceeds that of LDWPSO after about 8000 function evaluations and TRPSO can get a result only a little more accurate than that derived by LDWPSO. That is because there are too many local optima in function Rastrigin and the improvement

Figure 3. Logarithmic plots of mean function value versus function evaluations of TRPSO and LDWPSO on Rastrigin of 50 dimensions.

Figure 4. Logarithmic plots of mean function value versus function evaluations of TRPSO and LDWPSO on Griewank of 50 dimensions and amplified comparison plot of convergence for the beginning dozens of steps.

of TR method on PSO for this function is not so evident. Also, it is clear from Table 2 that another classic gradient method 3368



Article

size M = 30; max function evaluations Max_Fes = 600; wmax = 0.9, wmin = 0.4, c1= 2, c2 = 2, vmax = (χmax − χmin)/2, vmin = −vmax. Fifty independent runs have been executed for each algorithm. Figure 5 shows the comparison of fitness convergence between

SQP can hardly solve this benchmark. From Figure 3, it is seen that in the TRPSO algorithm, the downward trend is more obvious after about 6000 function evaluations, which means larger probability in getting a more accurate result. For function Griewank, which can be efficiently solved by gradient methods, in the TRPSO algorithm, the fitness values decrease especially quickly in the beginning several iterations and run into the satisfactory really fast. In order to describe the convergence process more clearly, the convergence plot of the beginning 50 function evaluations is shown here in Figure 4. It is obvious that TRPSO provides more accurate results with less computing time for NLP problems.

4. APPLICATION AND RESULTS DISCUSSION In a real plant, the polymerization reaction of high density polyethylene occurs mainly in the two reactors in serial configuration. A model for grade transition process has been reported by Wang et al.13 The two commercial polyethylene grades A and grade B and their production condition of the steady-state process are shown in Tables 3 and 4. Our proposed algorithm is applied to the Table 3. Properties of Grade A and B MI (g·(10 min)−1)

Figure 5. Plots of mean fitness value versus function evaluations of TRPSO and LDWPSO.

ρ (g·cm−3)

variable

first reactor

second reactor

pelletizer

first reactor

second reactor

grade A grade B upper limit lower limit

1000 140 +10% −10%

0.06 0.95 +10% −10%

0.225 0.65 +10% −10%

0.949 0.960 +5% −5%

0.953 0.964 +5% −5%

TRPSO and LDWPSO while optimizing the grade transition model. The TR stage in the red box is magnified to show this convergence process more clearly. It can be concluded that both searching accuracy and efficiency have been greatly improved for TRPSO in comparison with LDWPSO. The optimal trajectories calculated by TRPSO are shown in Figures 6−8, in which parts a−d are the trajectories of the instantaneous melt index, the cumulative melt index, the instantaneous density, and the cumulative density changing over time, respectively, and parts e−g are the manipulated variable profiles. The horizontal line in the graphs represents the steady state operation. In the first reactor, the transition involves a sharp decrease of MI from 1000 g·(10 min)−1 to 140 g·(10 min)−1 and an upward change of density from 0.949 g·cm−3 to 0.960 g·cm−3. The transition is accomplished by adjusting the temperature and feeding flow rate of hydrogen. The temperature is decreased immediately and then gradually increased to 83 °C at the steady-state point. The feeding flow rate of hydrogen is reduced to decrease the value of [H2]/[C2], which is first decreased to 3.33 for about 5 h and then gradually increased to 3.7 at the steady-state point. In the operation, the instantaneous melt index of polymer quickly reduces to 20 g·(10 min)−1 or so, and then back up to the target of 140 g·(10 min)−1. The cumulative melt index is more smoothly transferred to the target, as shown in Figure 6a and b. In the second reactor, the transition involves a sharp increase of MI from 0.06 g·(10 min)−1 to 0.95 g·(10 min)−1 and an upward change of density from 0.953 g·cm−3 to 0.964 g·cm−3. The transition is accomplished by adjusting the temperature and feeding flow rate of hydrogen and comonomer butylene. The temperature is increased immediately to 85 °C and then gradually decreased to 80 °C at the steady-state point. The feeding flow rate of hydrogen is improved to increase the value of [H2]/[C2], which is first increased to 0.308 for about 2 h and then gradually decreased to 0.28 at the steady-state point. The feeding flow rate of butylene is reduced to decrease the value of [C4]/[C2], which is first decreased to 0.007 and then gradually increased to 0.00778 at the steady-state point. In the operation,

Table 4. Manipulated Variables of Grade A and B first reactor

second reactor

variable

T (°C)

[H2]/[C2]

T (°C)

[H2]/[C2]

[C4]/[C2]

grade A grade B upper limit lower limit

85 83 90 75

5.8 3.7 +10% −10%

78 80 85 70

0.06 0.95 +10% −10%

0.0225 0.00778 +10% −10%

same process, and the results are compared with that in literature and in plant. The goal of this case is to minimize the amount of offspecification products while the polymers produced change from grades A to grade B to complete the transition process. For the NLP model derived in section 2 after using the simultaneous strategy, parameters A0−A5, B0−B3, and D0−D4 take the fitted values as shown in the literature, and Fr and w take the experience value of 0.22 and [0.55, 0.022, 0.1375, 0.0055, 0.45, 0.018, 0.1125, 0.0045,1.0]T respectively.13 In the proposed model, there are five control variables and nine state variables. Since in a real plant the transition lasts for 20 h, set t0 = 0 h, tf = 30 h. Since the manipulated variables can be adjusted every hour, divide the transition time span into Nfe = 30 equal pieces, which means 30 finite elements. The simultaneous approach is applied to the model with three-point Radau collocation with the location of ρ1 = 0.155, ρ2 = 0.645, ρ3 = 1. As shown in Table 3, all of the cumulative state variables have their upper and lower limits, which are called state constraints. The state constraints are all dealt with by the penalty function method to turn the constrained nonlinear optimization problem into a bound constrained nonlinear programming problem. Then, the algorithm parameters for TRPSO and LDWPSO are set as follows: dimension N = 5 × 30; swarm 3369



Article

Figure 6. Optimal variable profiles for grade transition from A to B in the first reactor.

Figure 7. Optimal variable profiles for grade transition from A to B in the second reactor.

transition in the pelletizer lasts from 0.12 h to 2.88 h. While the transition process is terminated, all polymer products in two reactors and pelletizer are required to achieve the requirements of the target grade. Therefore, the transition time calculated from solving the objective function is 6.23 h. In the commercial plant, the operations were taken according to experiences. There was rarely any overregulation for each of the operational variable. All operational variables were tuned gradually to the objective value step by step. This ensured a

the instantaneous melt index of polymer quickly rises to 3.8g·(10 min)−1 or so, and then back down to the target of 0.95g·(10 min)−1. The cumulative melt index is more smoothly transferred to the target, as shown in Figure 7a and b. Since the two grades have similar density, only the transition of the melt index is considered to calculate the transition time. It can be seen from Figures 6b, 7b, and 8 that the transition in the first reactor starts at 0.18 h and ends at 5.4 h, the transition in the second reactor starts at 0.07 h and ends at 6.3 h, and the 3370



Article

Table 5. Comparisons of Time and Quantity of Off-Spec Product of Grade Transition TRPSO optimal

IDP optimal

commercial operation

transition strategy

first reactor

second reactor

pelletizer

first reactor

second reactor

pelletizer

first reactor

second reactor

pelletizer

part transition time (h) total transition time (h) off-spec product (kg)

5.22

6.23 6.23

2.76

7.9

4.43 7.9

4.96

20

14 20

8

6.6 × 104

8.9 × 104

2.5 × 105

with the help of derived trajectories of grade properties and manipulated variables, it is much easier to use process measurements to adjust input variables for optimal grade transition implementation.

5. CONCLUSIONS A simultaneous strategy and TRPSO-based methodology to optimize grade transition process was developed and tested. The simultaneous strategy was used to fully discretize both state and control variables in the objective function model established for the grade transition process, which leads to a large-scale NLP problem. Then, the TRPSO, which uses PSO as the global optimizer while using TR to accelerate convergence of particles in local field, was proposed to solve NLP. Three optimization algorithms for NLP, including TRPSO, LDWPSO, and SQP were tested on three multimodal benchmarks and the comparison results showed that TRPSO outperforms the other two algorithms on both convergence precision and robustness. At last, the combination of simultaneous strategy and TRPSO was used to solve a commercial grade transition optimization in a real plant. The comparison results with the existing IDP optimization and actual data before optimization proved the validity of the proposed combination of simultaneous strategy and TRPSO algorithm. It will have a promising potential for practical use in grade transition optimization.

Figure 8. Optimal MI profile for grade transition from A to B in pelletizer.

smooth transition to the properties of the objective product. However, it took a really long time for the transition to complete. A more reasonable operation profile was proposed from IDP.13 It has been proved that [H2]/[C2] is the primary variable used to control MI. So, the feeding flow rate of hydrogen in the first reactor is reduced to decrease the value of [H2]/[C2], which is first decreased to 3.2 for about 2 h and then gradually increased to 3.5 at the steady-state point. In the second reactor, the feeding flow rate of hydrogen is improved to increase the value of [H2]/[C2], which is first increased to 0.6 for about 3 h and then gradually decreased to 0.25 at the steady-state point. In both reactors, the feeding flow rate of butylene is reduced properly to decrease the value of [C4]/[C2]. Thus, the density increases gradually. In the first reactor, the instantaneous melt index of polymer started to deviate from that of grade A after 0.145 h of transition. The production of off-specifications started. After 8.135 h of transition, acceptable products of grade B started to be produced. The transition time of reactor one was supposed to be 7.9 h. Since the transition progress in reactor one was much slower than that in the second reactor and in the pelletizer, the whole transition time was also 7.9 h, which was much less than the actual transition time in commercial plant, but more than that derived by TRPSO. Therefore, we can conclude that the profiles derived by TRPSO are superior to that by IDP and commercial plant. Table 5 compares the optimization results with that derived by IDP13 and the transition time and off-specifications produced in actual commercial operation. It is concluded that the transition time required is least if the production process is executed according to trajectory calculated using TRPSO. Since the off-specifications produced varies directly as transition time, according to the formula given in the literature,13 based on simultaneous strategy and TRPSO described, the off-specifications produced during transition are 6.6 × 104 kg, 73.6% less than 2.5 × 105 kg, which is the actual situation in commercial operation, and 25.8% less than 8.9 × 104 kg, which is the optimization result of IDP. Thus, it is proven that using the method proposed has a promising potential for notable economic benefits. In a real plant, the rigorous models are insufficient to predict product quality due to the presence of uncertainty in the form of model mismatch and process disturbances. However,

■

AUTHOR INFORMATION

Corresponding Author

*Tel: +86-10-64426960. Fax: +86-10-64437805. E-mail: [email protected]. Present Address ||

College of Information Engineering, Beijing Institute of Petrochemical Technology, Beijing, 102617

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (Grant No. 61074153), the National High Technology Research and Development Program (863, Grant No. 2006AA04Z184), and the Fundamental Research Funds for the Central Universities (Grant No. zz1136); their support is hereby acknowledged.

■

NOMENCLATUREINDICES i = 1, ···, NZ differential state variables j = 1, ···, NY algebraic variables l = 1, ···, Nfe finite elements q, q′, r = 1, ···, Ncp collocation points Parameters

NZ NY NU Nfe 3371

number number number number

of of of of

differential state variables algebraic variables manipulated variables finite elements


Industrial & Engineering Chemistry Research Ncp Ai, Bi [H2/C2] [C4/C2] Fr J MI ρ P T t u z wi

Article

(16) Shi, Y.; Eberhart, R. C. A modified particle swarm optimizer. Proc. IEEE Congr. EVol. Comput.; IEEE Press: Piscataway, NJ, 1998; pp 69−73. (17) Geng, Z. Q.; Zhu, Q. X.; Gu, X. B.; Lin, X. Y. Optimal control of cracking depth based on multi-swarm competitive PSORBFNN for ethylene cracking furnace. CIESC J. 2010, 61 (8), 1942−1948. (18) Tseng, L. Y.; Chen, C. Multiple trajectory search for large scale global optimization. Proc. IEEE Congr. EVol. Comput.; IEEE Press: Piscataway, NJ, 2008; pp 3052−3059. (19) Zhao, S. Z.; Liang, J. J.; Suganthan, P. N.; Tasgetiren, M. F. Dynamic multi-swarm particle swarm optimizer with local search for large scale global optimization. Proc. IEEE Congr. EVol. Comput.; IEEE Press: Piscataway, NJ, 2008; pp 3845−3852. (20) Hsieh, S. T.; Sun, T. Y.; Liu, C. C.; Tsai, S. J. Solving large scale global optimization using improved particle swarm optimizer. Proc. IEEE Congr. EVol. Comput; IEEE Press: Piscataway, NJ, 2008; pp 1777−1784. (21) Li, J. B.; Liu, X. G. Melt index prediction by RBF neural network optimized with an MPSO-SA hybrid algorithm. Neurocomputing 2011, 74 (5), 735−740. (22) Moré, J. J.; Sorensen, D. C. Computing a trust region step. SIAM J. Sci. Stat. Comput. 1983, 3, 553−572. (23) Coleman, T. F.; Verma, A. A preconditioned conjugate gradient approach to linear equality constrained minimization. Comput. Optim. Appl. 2001, 20 (1), 61−72. (24) Xu, W. X.; Geng, Z. Q.; Zhu, Q. X.; Gu, X. B. A piecewise linear chaotic map and sequential quadratic programming based robust hybrid particle swarm optimization. Inf. Sci. 2013, 218 (1), 85−102.

number of collocation points parameter of function ratio of hydrogen to ethylene ratio of comonomer to ethylene fraction of polymer produced in the reactor before (t + Δt) objective function melt index, g·(10 min)−1 density, g·cm−3 properties of polymer temperature of polymerization, °C time, h trajectory of manipulated variable trajectory of state variable weighting factor in objective function

Subscripts

c i low up

■

cumulative property instantaneous property lower limit of manipulated variables upper limit of manipulated variables

REFERENCES

(1) Debling, J. A; Hun, G. C.; Kuijpers, F.; VerBurg, J.; Zacca, J.; Ray, W. H. Dynamic modeling of product grade transitions for olefin polymerization process. AIChE J. 1994, 40 (3), 506−620. (2) McAuley, K. B.; MacGregor, J. F. Optimal grade transitions in a gas-phase polyethylene reactor. AIChE J. 1992, 38 (10), 1564−1576. (3) Takeda, M.; Ray, W. H. Optimal grade transition strategies for multistage polyolefin reactors. AIChE J. 1999, 45 (8), 1776−1793. (4) Wang, Y.; Seki, H.; Ohyama, S.; Akamatsu, K.; Ogawa, M.; Ohshima, M. Optimal grade transition control for polymerization reactors. Comput. Chem. Eng. 2000, 24 (2−7), 1555−1561. (5) BenAmor, S.; Doyle, F. J.; McFarlane, R. Polymer grade transition control using advanced real-time optimization software. J. Process. Control 2004, 14 (4), 349−364. (6) Bonvin, D.; Bodizs, L.; Srinivasan, B. Optimal grade transition for polyethylene reactors via NCO tracking. Chem. Eng. Res. Des. 2005, 83 (A6), 692−697. (7) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldivar-Guerra, E. Dynamic optimization of HIPS open-loop unstable polymerization reactors. Ind. Eng. Chem. Res. 2005, 44 (8), 2659−2674. (8) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldivar-Guerra, E. Optimal grade transitions in the HIPS polymerization process. Ind. Eng. Chem. Res. 2006, 45 (18), 6175−6189. (9) Kameswaran, S.; Biegler, L. T. Simultaneous dynamic optimization strategies: Recent advances and challenges. Comput. Chem. Eng. 2006, 30 (10−12), 1560−1575. (10) Flores-Tlacuahuac, A.; Biegler, L. T. Simultaneous mixed-integer dynamic optimization for integrated design and control. Comput. Chem. Eng. 2007, 31 (5−6), 588−600. (11) Flores-Tlacuahuac, A.; Biegler, L. T. Integrated control and process design during optimal polymer grade transition operations. Comput. Chem. Eng. 2008, 32 (11), 2823−2837. (12) Wang, J. D.; Yang, Y. R. Study on optimal strategy of grade transition in industrial fluidized bed gas phase polyethylene production process. Chin. J. Chem. Eng. 2003, 11 (1), 1−8. (13) Wang, P.; Wang, J. D.; Yang, Y. R. Optimization of grade transition for high-density polyethylene continuous slurry process. J. Chem. Ind. Eng. (China) 2006, 57 (11), 2657−2663. (14) Lee, H. Y.; Yang, T. H.; Chien, I. L.; Huang, H. P. Grade transition using dynamic neural networks for an industrial highpressure ethylene−vinyl acetate (EVA) copolymerization process. Comput. Chem. Eng. 2009, 33 (8), 1371−1378. (15) Kennedy, J.; Eberhart, R C. Particle swarm optimization. Proc. IEEE Congr. EVol. Comput.; IEEE Press: Piscataway, NJ, 1995; pp 1942−1948. 3372


Optimal Grade Transition for Polyethylene Reactors Based on

Recommend Documents