Yujie Wei,Yongheng Jiang,Dexian Huang*

Process Systems Engineering and Process Safety

A Three-section Algorithm of Dynamic Programming Based on Three-stage Decomposition System Model for Grade Transition Trajectory Optimization Problems☆

Yujie Wei,Yongheng Jiang,Dexian Huang*

National Laboratory for Information Science and Technology,Department of Automation,Tsinghua University,Beijing 100084,China

A R T I C L EI N F O

Article history:

Gradetransition trajectory optimization

Adaptivegrid allocation

Heuristic modif i cations

Three-section dynamic programming

Three-stage decomposition model

This paper introduces a practical solving scheme of gradetransition trajectory optimization(GTTO)problems under typical certi fi cate-checking-updating framework.Due to complicated kinetics of polymerization, differential/algebraic equations(DAEs)always cause great computational burden and system non-linearity usually makes GTTO non-convex bearing multiple optima.Therefore,coupled with the three-stage decomposition model,a three-section algorithm of dynamic programming(TSDP)is proposed based on the general iteration mechanism of iterative programming(IDP)and incorporated with adaptivegrid allocation scheme and heuristic modi fi cations.The algorithm iteratively performs dynamic programming with heuristic modi fi cations under constant calculation loads and adaptively allocates the valued computationalresources to the regions that can further improve the optimality under the guidance of local error estimates.TSDP is fi nally compared with IDP and interior point method(IP)to verify its ef fi ciency of computation.

?2014TheChemicalIndustry andEngineeringSocietyofChina,andChemicalIndustryPress.Allrightsreserved.

1.Introduction

Nowadays,diverse polymer products have been widely used in many standard applications including automobile industries,electronics and packaging materials[1,2].These polymers are in demand with different physical and chemical specif i cations,i.e.,different grades. Therefore,to meet the various requirements and high quality standards from industrial consumers,many different products must be produced and transitions are inevitable.However,during transitions,two different kinds of products are mixed with long reactor residence time (e.g.,2-4 h),thus resulting in a long transition time(e.g.,5-30 h)and a large amount of off-spec polymers,which have to be sold at a much lower price[2,3].Over the last three decades,the need to increase prof i tability in polymer industrial productions has given rise to many attracting research problems such as dynamic modeling of polymerization reactors(or quality index)[4,5],reference trajectory optimization off l ine[6-8],monitoring quality index online and hierarchical integration optimization[9].Among these topics,the dynamic models of processes and trajectory optimization are the basal but challenging problems,which will be the concern in this paper.

For these basic dynamic optimization problems,one common used method is sequential method with non-linear differential/algebraic equations(DAEs),in which only the control variables are discretized. McAuley employed a gradient based nonlinear programming problem (NLP)optimizer updating feasible solutions which are integrated to determine the resulting trajectories of the state variables based on f i rstprinciple models[7].Takeda and Ray proposed the methodology that utilized a dynamic process simulator embedded in the algorithm iterations,and coupled to a sequential quadratic programming(SQP)[8]. Wang and Yang presented a grade transition model with state constraints and chosen SQP to solve it[6].Related studies can be referred to He et al.(2012)[10],Wang and Yang(2003)[6]and Fei et al.[11]. However,Lee et al.pointed out that non-linear and multimodal chemical systems had diff i culties in the convergence to the global optimum and therefore,they proposed a modif i ed differential evolution algorithm[12,13].

Biegler et al.believed that sequential methods with repeated integration of the DAE may become time-consuming[14,15].Therefore, they adopted simultaneous approaches which discretized state and control variables and transformed the original problem into a much larger non-linear programming problem(NLP)which was solved by the interior point strategy[15].Fei et al.indicated that the development of mechanism models based on rigorous theoretical knowledge of physicochemistry was usually quite diff i cult and may be impractical [11].Considering the diff i culties and practicality of modeling,they proposed a new data-based model using ARX(AutoRegressive with eXogenous inputs)combined with neural network under partial least squares framework(ARX-NNPLS),in which less specif i c knowledge of the process was required but the input and output data[11]. The summary of the up-to-now research work is illustrated in Table 1.

As quality demands have gradually increased in recent years,they cause great pressure and motivations to make polymer productions more eff i cientandcost-effective.Thus,theaccuracy ofquality measurements andthef l exibilityandprof i tabilityofgradetransitionshavetobe promoted further[1].

For the expressions of dynamic system models,DAEs are often diff i cult to establishand will causetime-consumingevaluation processes forfeasiblesolutions,whereas,data-drivenmodels as localmodels f i t trainingsetswellbutmaymismatchthetruesystemsoverawiderange. Thereby,ifstructuralknowledgeanddata-drivenmethodscanbeincorporated,it will share advantages of these two kinds of approaches.Due to the physical structural information involved,the model can capture the nature of the wide range of process f l uctuations,while on the otherhand,theunknownnon-linearities,dynamicsofhighlycomplicatedkineticbehavioralongwithunknownnoisescanbeapproximatedby the data-driven methods.A three-stage decomposition model based on structural information to describe the dynamic non-linear behavior of dynamic transitional processes has been developed and the estimation accuracy is validated to be effective in comparison with black-box modeling methods[20].

Foroptimizationalgorithms,dynamicprogramming(DP)isatypical solution method for dynamic optimization problems.The advantages are independence on problem structures and gradient information, and can obtain optimal solutions based on the principle of optimality. However,in most cases,the lack of closed form solutions of Bellman equations always leads to the use of computational methods—such as value iteration and policy iteration—to solve those equations,typically requiring computation time and memory that grow exponentially in the number of state variables,decision variables and random variables that affect the system in each time period,which are“curses of dimensionality”[21].Just these diff i culties render DP infeasible for high-dimensional problems of extensively practical use[21].

In recent years,in order to overcome the“curses of dimensionality”, iterative dynamic programming(IDP)and adaptive grid scheme for DP have been proposed.IDP is f i rst proposed by Luus[22]and has been extensively applied[23,24].The procedure of IDP is performed for a numberofiterationsofDP.Theoptimalcontrolpolicyfromtheprevious iteration gives the central values of the search regions at each stage of the current iteration.The radii of corresponding regions are contracted by a constant to provide a f i ner resolution.Grune and Semmler proposed an adaptive grid scheme for the discretization of state space to fi nd the global solutions of discrete time Bellman equations[25,26]. Local posteriori error estimates were established to guide the computation to the re fi nement on certain critical regions of the whole space.

For IDP,search regions at each time stage are contracted by a constant factor which is assigned in advance.Therefore,the computational resources may not be adaptively adjusted to concentrate on the interesting regions.For an adaptive grid scheme,the task of algorithm is just to solve a fi xed-point equation without a practical trajectory used for operations,therefore neglecting some heuristic improvements at certain stages of the time horizon.Contrarily,if the local posteriori error estimate equations proposed by Grune can be incorporated into IDP,itwillleadustoselectaf i nerregionandconcentrateonthecomputations,and consequently the approximately optimal solutions can be achieved under a certain amount of calculations.Also,if the heuristic modif i cations can be performed at a particular time stage,a better solution in a short time can also be obtained[27].Thereby,on one hand,relying on the principle of optimality of dynamic programming, optimal or approximately optimal solutions can be obtained in principle,while on the other hand,iteration mechanisms of IDP coupled with the adaptive grid scheme and heuristic modif i cations can eff i ciently allocate computational resources and overcome the curse of dimensionality.

In order to increase the algorithmic eff i ciency of grade transition problems,a three-section algorithm of dynamic programming(TSDP) coupled with the three-stage system model is proposed in this paper. Similarly with IDP,this novel algorithm starts from a DP algorithm of coarse discretization with fast computation characteristics.The following idea that improves eff i ciency is to design iteration cycles with heuristic modif i cations under limited calculation loads and to allocate adaptively the valued computational resources to the regions that can further improve the optimality of feasible solutions with the guidance of local posteriori error estimates.

2.Problem Formulations

The grade transition problem in polymerization industries is a dynamic optimization problem.The goal is to f i nd the input trajectories during transitions to have the best dynamic performance under certain objective functions,considering states and outputs of a system in a dynamical process[8].Therefore,system models should be included.

The general continuous grade transition optimization problem can be written as the following formulation.

Table 1The summary of the up-to-now research work

where s(t)is the m-dimensional state vector and u(t)is the n-dimensional input vector.The objective function is the integration form of continuous state and input trajectories from start time 0 to the end time tf,in which v(s(t),u(t))is the current cost function of state s(t)and input u(t)at time t,and α(t)is the time-weight factor.

Themostchallengingproblemfordynamicoptimizationproblemsis the complicated non-linear dynamic system behavior due to the polymerization kinetics and energy effects[16].On one hand,it is diff i cult to accurately describe the unknown non-linearity over a wide range. On the other hand,the non-linearity usually makes grade transition a non-convex optimization problem bearing multiple optima[16].The complexity of algorithm applied up to now that intends to f i nd optimal solutions will always grow exponentially with respect to the size of the input.Therefore,approximation algorithm is an alternative approach that is less ambitious to aim in f i nding exact optimal solutions.

The off l ine optimization strategy research is typically performed using a dynamic model in conjunction with optimization.In the search for feasible solutions,two algorithms are usually designed to f i nish the optimization,which are the certif i cate-checking algorithm around solutions at hand and the updating algorithm to f i nd new solutions as illustrated in Fig.1[28].Many progresses have been made under this generally accepted framework[7,8,13-15].

The certif i cate-checking algorithm around solutions at hand is mainly the scheme dealing with the complex system constraints.The objectives of the algorithm are checking whether the solutions arefeasible or not and evaluating system performances under the decision inputs.Consequently,it is closely related to the expressions of dynamic system models.

The updating algorithm to f i nd new solutions attempts to f i nd out a better solution based on the evaluation of existing solutions and some other information about systems,constraints and objective functions. One common employed is thegradientvalue.Consequently,it is closely related to optimization algorithms.

In this paper,the three-stage system model in our previous work willbeadoptedasthesystemexpressionsandathree-sectionalgorithm of dynamic programming(TSDP)will be proposed as a novel optimization algorithm to update solutions,which is illustrated in Fig.2.

Fig.1.Framework utilized in the gradetransition strategy research.

3.Three-stage Decomposition System Model

A three-stage decomposition model based on the structural analysis hasbeendevelopedandthusitcancapturethenatureofthewiderange of process f l uctuations[20].The most important is that compared with DAEs in which feasible solutions can be evaluated more eff i ciently due to their explicit expressions and applications of data-driven methods.

The three-stage architecture is established on the basis of the structural knowledge as shown in Fig.3[20].Relying on the structural information and intermediate variables(i.e.,reactant concentrations),the integral process can be decomposed into three stages as the“white”structure part of the model in which each stage is much easier to be dealt with[20].That is because diff i culties incurred by dynamic characteristic,insuff i ciency of training data and non-linearity are no longer dealt with simultaneously but sequentially centralized in different stages.Meanwhile,the“black”part due to the unknown dynamic non-linearities and noises can be approximated from data.

In this model,the f i rst stage describes a dynamic balance of gasreactant concentrations maintained by operational variables using the modeling algorithm of adaptive hinging hyperplanes(AHH)[20]. In the second stage,the variables of reaction conditions determine instantaneous quality index based on the mechanism analysis.Finally, thethirdstagecharacterizestherelationofinstantaneousandcumulative quality indices,employing impulse response template(IRT).

Fig.2.Model and algorithm proposed in this paper under the certif i cate-checkingupdating framework.

3.1.First-stage model based on AHH algorithm

The f i rst-stage model is expressed as follows.

where r1,r2,…,rm,rT,rPare the orders of the corresponding variables, Ri(i=1,…,m)represents the concentration of the i-th reactant,T is the temperature of the reaction,P is the pressure of the reactor and ui(i=1,…,m)are the operational variables.

The function of g is approximated in terms of linear combinations of hinging functions and the details of order testing and modeling parameter selection can be referred to the studies.

3.2.Second-stage model using non-linear parameter estimation methods

Ingeneral,thereexists anempiricalrelationshiptodescribehowthe reaction conditions determine instantaneous quality index,which is written as

Fig.3.Three-stage quality index modeling paradigm.

In the feasible region of parameters{K,Tf,τf},the error e(t)= Qc(t)?Qcf(t)between the estimated value and actual value is minimized.

4.Three-section Dynamic Programming with Heuristically Adaptive Grid Scheme

This section will propose the framework and algorithms of threesection dynamic programming incorporating adaptive grid scheme and heuristic modif i cations into the general iteration mechanism of IDP.Under the limitation of certain amounts of calculations,the aim of this novel algorithm is to adaptively allocate the computational resources to the regions that can further improve the optimality of feasible solutions with the guidance of local posteriori error estimates proposed by Grune.

4.1.Framework of three-section dynamic programming with heuristically adaptive grid scheme

The general iteration mechanism of IDP is employed as the basic framework of TSDP,which is illustrated in Fig.4.In this framework, the standard dynamic programming with limited calculations is iteratively performed which is the f i rst-section of the algorithm. After standard DP was performed,heuristic modif i cations are implemented and the objective value can be improved,which is the second-section of the algorithm.

Fig.4.Framework of algorithm of TSDP.

The third-section is concerning the adaptive grid scheme,based on which the next iteration will be performed and a better solution can be obtained.In the algorithm of IDP,the radius of the search region will be contracted by a constant in the current iteration and the number of grids in the search region of each time stage is assigned in advance and f i xed throughout the algorithm.However, in the algorithm of TSDP,two-sequential steps are included:adaptively selecting the radius of the search region and adaptively allocating the number of grids in each search region.The adaptive selection of the search radius is carried out based on the improvement ranges of heuristic modif i cations in the second section.After radiuses are selected,grid allocation scheme will be presented explicitly.The goal of the gird allocation is to minimize the ratio of intervals in the current and last iterations for the same time-stage under the limitation of total calculations.The differences between IDP and TSDP are summarized explicitly in Table 2.

In summary,three elements have been incorporated into the algorithm of TSDP:the general iteration mechanism of IDP with limited calculations,heuristic modif i cations and grid adaptive allocation scheme,which are the reasons that TSDP can be performed with higher eff i ciency.

4.2.First-section standard dynamic programming

In the f i rst section,DP algorithm starts from coarse discretizations with fast computation characteristics.The dynamic programming is developed backward in time as

In this equation,fkn(skn)represents the value of cost-to-go function at the initial state skn,which is the initial state value of the k-th stage in the n-th iteration.Similarly,uknis the input variable of the k-th stage in the n-th iteration as shown in Fig.5.

Thesearchregionatstagekinthen-thiterationis Ukn,which is ref i ned as a proper subset of initialized search space U after grid adaptive allocation.The expression form of Uknis the grid point array as shown in Fig.5.

4.3.Second-section heuristic modif i cations

As the control-vector parameterization is applied,the input value is kept constant in each time interval.As a result of time discretizations and zero-order holding,the operational variable is still kept at a constant that used to be optimal when the output trajectory has exceeded the specif i ed target line.

The occasionsin whichtheobjective can be improvedare illustrated explicitly as an example in Fig.6.Let us now concentrate on the third time-stage in Fig.6(a)and(b).Fig.6(a)is the output trajectory and (b)is the corresponding input trajectory.Fig.6(c)and(d)display partially enlarged details of the third time stage.From the f i gure,we can f i nd that at time of 42,the initial point of the third time-stage, the input value of this stage is suitable,but from the time of 55 to 61 this value is no longer optimal because the output curve has exceeded the target line.If the input value in this tail period could be slightly decreased with some fl exibility,a better output trajectory could be obtained.

From the above analysis,a suf fi cient condition that can judge whether there exist improvements at one time point can be concluded as follows.

Condition:For a gradetransition problem with time length L under N-stage dynamic programming the optimal output trajectory y( t)can be obtained with the optimal input trajectory u( t)as

Table 2Summary of the differences between IDP and TSDP

4.4.Search-radius adaptive selection scheme

The f i rst difference between IDP and TSDP is the search-radius selection as mentioned in Table 2.In IDP,the radius is contracted by a constant from the original radius.But in TSDP,the radius selection is guided by the error gap estimated and adaptively enlarges the searchradius of time stages with larger errors.

The local error estimates equation proposed by Grune is thus introduced to give the gap ηknbetween fkn(skn)and f?(skn).In order to estimatehowgoodfkn(skn)approximatesf?(skn),Grunedef i nedaposteriori error estimates based on local values of the residual as follows[25].

Basedontheaboveerrorestimation,thesearch-radiuscanbeselected as half of the variation range of previous,current and next timestages if the error is low.If the error is large,the radius selected will depend on the variation range of more time-stages.

The search-radius selection rule(SRSR)is as follows.

(1)Iftheerrorofthek-thtime-stageηnk≤θ mkax ηnk,selectthesearchradius as

Grid array along time stages in the n-th iteration

Fig.5.Grid array.

4.5.Adaptive grid allocation scheme

The goal of the adaptive grid allocation is to minimize the ratio of intervals in the current and last iterations for the same time stage under the limitation of total calculations.In other words, it is to make the intervals of two adjacent iterations decrease as far as possible.

Fig.6.Illustration of heuristic modif i cation condition.

The adaptive grid allocation problem of two-dimensional decision variables can be similarly written as

where the variable dn+1,k,irepresents the number of grids on the k-th time-stage(k=1,2,…,N)in the n+1-th iteration for the i-th decision variable ui(i=1,2).The variable δn+1,k,iis the search-radius for the decision variable uion the k-th time-stage of the n+1-th iteration. The variable of σn,k,iis the discretization interval of the continuous input space of uifor the k-th time-stage of the n-the iteration.The relation of

Theobjectiveofthisoptimizationproblemistominimizetheratioof the intervals of two adjacent iterationsto make the intervals for the same time-stage decline.The constraint is the limitation of the total calculations of grid combinations dn+1,k,1×dn+1,k,2along the time-stages.

Let F be the Lagrangian relaxation of this problem

Furthermore,under the Karush-Kuhn-Tucker(KKT)conditions,the adaptive grid decisions can be stated as follows.

Different input variables may have different inf l uences on the performance index and require different minimum intervals.One can initially assign different grid numbers dn,k,ito different input variables on account of the experiential knowledge of different inf l uences.And then the variation of the search-radiuses between δn,k,iand δn+1,k,iwill adaptively lead to the adaptive adjustment from dn,k,ito dn+1,k,i.

Therelativerelationofradiusesalongthetime-stagesinthesameiteration is also considered in these two equations.For the larger radius, the grid allocated to this stage will be much larger and hence much more computational resources will be concentrated on this stage.

Therefore,this allocation scheme considers different inf l uences of different input variables,the radius variation of the same time-stages in adjacent two iterations and the relatives of radiuses along the time stages inthesame iteration,f i nallymakingthe intervalsof twoadjacent iterations decrease as far as possible.

5.Optimization Results for Gas-phase Polyethylene Process

5.1.Gas-phase polyethylene process

In this paper,gas-phase polyethylene process will be considered as an instance and the schematic diagram of the reactor system is shown inFig.7[20].Inthereactor,thecopolymerizationofethyleneandbutene is carried out using the Ziegler-Natta catalyst.

In the system,the fresh feed includes ethylene,butene,hydrogen and nitrogen,and the product of polymerization is discharged from near the base of the reactor[7].Hydrogen/ethylene and butane/ethylene of the fresh feed injected into the reactor are two operational variables mainly considered,which are controlled in the close-loop.The set-points of controllers are the input variables of system models and also decision variables for optimization problems.The melt index and density of products are mainly considered as the output variables,which are also the state variables and output variables of the optimization problems.

5.2.Optimization problems for the grade transition

A polymer grade is typically characterized by various physical properties.In this paper,the melt index and density are mainly considered.As for input operational variables,it is restricted to the set-points of ratio controllers of hydrogen/ethylene σH2/C2and butane/ethylene σC4/C2ofthefreshfeed.Itisbecausewhenoptimalgradetransitionstrategies are installed in an actual plant,it is better to have a simple strategy.The constraint values are summarized in Table 3.

The state variable vector is s(t)=(MI(t),ρ(t))T.The objective sta（te variable vector）is(MIobj,ρobj)T.The input variable vector is u( t)=σH2/C2(t),σC4/C2(t)T.The current cost function of the objective functions in this work is the form of quadratic function,i.e.,

Fig.7.Gas-phase polyethylene reactor system.

and the time-weight factor is αt=βt(0＜β＜1).

A case will be studied in the following subsection using the algorithm of TSDP.In order to test the performance of TSDP,IDP is also applied and put into comparisons.Up to now,the interior point method(IP)incorporated into simultaneous methods has been one of the most widely used methods[14,15,18].The famous algorithm is the IPOPT algorithm for large-scale non-linear programming which follows a primal-dual approach and applies a Newton method to the resulting optimality conditions[16].In the case studies,the IP that follows abarrier approach will also be implemented using the powerful solver of MATLAB's optimization toolbox.This solver has been used extensively by various researchers for the study of scheduling of the grade changes[8,11,16].

The algorithm of IDP was coded based on the MATLAB platform under the guidance of Luus[22].As studied by Bojkov and Luus[22], the region contraction factor should be in the range 0.75 to 0.95 for many typical chemical engineering systems.In this paper,this factor is chosen as 0.9 after some comparisons.The initial search region selected isthe whole range of input space as shown in Table 3.The method of IP is applied using the Optimizer(Fmincon in MATLAB 7.9)and the tolerances used have been reported in Table 4.

Table 3Bound constraints for variables associated with grade transition

5.3.Case study

AgradetransitioninpolyethyleneproductfromgradeAcharacterized by the melt index of 0.5697 and the density of 0.9314 to grade B with the melt index of 1.5314 and density of 0.9233 is considered,as presented in Table 5.

The optimization results of three algorithms are shown in Table 6. In the comparison experiments,it is shown how the discussed algorithm behaves for a grade transition problem with increasing number of time-stages.The number of time stages is 6,10 and 15, respectively.The initial number of grids per time-stage for TSDP is designed as 4,3,and 3.The number of grids for IDP is set to 3,which has been selected after many comparisons.Under these basic parameters sets,three algorithms of TSDP,IDP and IP are performed and the best objective value and CPU time are recorded in Table 6.

FromtheresultsofTable6,TSDPprovidesthebestobjectivefunction value with the highest eff i ciency than IDP and IP.

ForIDPandTSDP,whenthenumberoftime-stagesis6,TSDPobtains theobjectiveof82.4361 withintheCPU timeof 134.52,butIDP canonly achieve the objective of 83.6590 even with more computational time, because IDP is only the considerable iterations of many DPs but TSDP has incorporated heuristic information and adaptive grid scheme into the general framework of IDP.Thus the improvement based on the foundation of IDP can guide the computational resources to the critical regions and saves the solving time.

It was also found that the best objectives of IP and TSDP are very close,especially when the number of time-stages is 10,but the CPU time of IP increases sharply with the increase of the number of time stages.The CPU time of TSDP increases slightly at the same order of magnitude.This is because the total number of calculations of DP is limited and only a few iterations are needed under the guidance of the adaptive allocation scheme.Even if the number of stages is increased, the total calculations may stay constant as a result of sequential mechanisms of TSDPalgorithms.Themechanism of IPis notthe case.As incorporated into simultaneous methods,IP is used to solving NLP problems. With the increase of the number of time stages,the total number of decision variables that need to be solved simultaneously not sequentially, increases too.The enlargement of decision variables makes the NLP more diff i cult to solve and easily fall into local optimum.

Table 5Steady-state operating conditions for Grades A and B

Table 6Comparisons of TSDP,IDP and IP for grades A and B

6.Conclusions

In this paper,a practical solving scheme of grade transition trajectory optimization problems is introduced under the typical certif i catechecking-updating framework.

For certif i cate-checking algorithms,the three-stage system model is adopted as a practical model of dynamic optimization problems.This model not only can accurately describe dynamic non-linear behaviors of the systemsbut also is muchmore eff i cient for evaluations as a result of its explicit expressions.

For updating algorithms,a three-section algorithm of dynamic programming(TSDP)based on the general iteration mechanism of iterative programming(IDP)is proposed incorporated with the adaptive grid allocation scheme and heuristic modif i cations.The algorithm is designed as the iterations of DP with heuristic modif i cations under limited calculation loads and can adaptively allocates the valued computational resources to the regionsthat canfurther improve the optimality under the guidance of local error estimates.TSDP is fi nally compared with IDP and interior point method(IP)to verify its ef fi ciency of computation.

[1]G.Dunnebier,D.van Hessem,J.V.Kadam,K.U.Klatt,M.Schlegel,Optimization and control of polymerization processes,Chem.Eng.Technol.28(5)(2005)575-580.

[2]P.O.Larsson,J.Akesson,N.Carlsson,N.Andersson,Model-based optimization of economical grade changes for the Borealis Borstar(R)polyethylene plant,Comput. Chem.Eng.46(2012)153-166.

[3]C.Chatzidoukas,S.Pistikopoulos,C.Kiparissides,A hierarchical optimization approach to optimal production scheduling in an industrial continuous ole fi n polymerization reactor,Macromol.React.Eng.3(1)(2009)36-46.

[4]K.B.McAuley,J.F.MacGregor,On-line inferenceof polymerpropertiesinan industrial polyethylene reactor,AICHE J.37(6)(1991)825-835.

[5]B.Kou,K.B.McAuley,J.C.C.Hsu,D.W.Bacon,Mathematical model and parameter estimation for gas-phase ethylene/hexene copolymerization with metallocene catalyst,Macromol.Mater.Eng.290(6)(2005)537-557.

[6]J.D.Wang,Y.R.Yang,Study on optimal strategy of grade transition in industrial fl uidized bed gas-phase polyethylene production process,Chin.J.Chem.Eng.11 (1)(2003)1-8.

[7]K.B.McAuley,J.F.MacGregor,Optimal grade transitions in a gas phase polyethylene reactor,AICHE J.38(10)(1992)1564-1576.

[8]M.Takeda,W.H.Ray,Optimal-grade transition strategies for multistage polyolef i n reactors,AICHE J.45(8)(1999)1776-1793.

[9]A.Prata,J.Oldenburg,A.Kroll,W.Marquardt,Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor,Comput. Chem.Eng.32(3)(2008)463-476.

[10]D.HE,T.ZOU,L.YU,Receding horizon control for grade transition of large-scale polypropylene production plants,Chin.J.Chem.Eng.024(2010).

[11]Z.S.Fei,B.Hu,L.B.Ye,J.Liang,ARX-NNPLS model based optimization strategy and its application in polymer grade transition process,Chin.J.Chem.Eng.20(5)(2012) 971-979.

[12]M.H.Lee,C.Han,K.S.Chang,Hierarchical time-optimal control of a continuous copolymerization reactor during start-up or grade change operation using genetic algorithms,Comput.Chem.Eng.21(1997)S1037-S1042(Suppl.).

[13]M.H.Lee,C.Han,K.S.Chang,Dynamic optimization of a continuous polymer reactor using a modif i ed differential evolution algorithm,Ind.Eng.Chem.Res.38(12)(1999) 4825-4831.

[14]L.T.Biegler,V.M.Zavala,Large-scale nonlinear programming using IPOPT: an integrating framework for enterprise-wide dynamic optimization,Comput. Chem.Eng.33(3SI)(2009)575-582.

[15]A.M.Cervantes,S.Tonelli,A.Brandolin,J.A.Bandoni,L.T.Biegler,Large-scaledynamic optimization for grade transitions in a low density polyethylene plant,Comput. Chem.Eng.26(2)(2002)227-237.

[16]N.Padhiyar,S.Bhartiya,R.D.Gudi,Optimal grade transition in polymerization reactors:a comparative case study,Ind.Eng.Chem.Res.45(10) (2006)3583-3592.

[17]H.Pakravesh,A.Shojaei,Optimization of industrial CSTR for vinyl acetate polymerization using novel shuff l ed frog leaping based hybrid algorithms and dynamic modeling,Comput.Chem.Eng.35(11)(2011)2351-2365.

[18]A.Flores-Tlacuahuac,L.T.Biegler,E.Saldivar-Guerra,Dynamic optimization of HIPS open-loop unstable polymerization reactors,Ind.Eng.Chem.Res.44(8) (2005)2659-2674.

[19]A.Flores-Tlacuahuac,L.T.Biegler,E.Saldivar-Guerra,Optimal grade transitions in the high-impact polystyrene polymerization process,Ind.Eng.Chem.Res.45(18) (2006)6175-6189.

[20]Y.Wei,Y.Jiang,F.Yang,D.Huang,Three-stage Decomposition Modeling for Quality of Polymerization Processes Based on Adaptive Hinging Hyperplanes and Impulse Response Template,Ind.Eng.Chem.Res.52(16)(2013)5747-5756.

[21]W.B.Powell,A.George,B.Bouzaiene-Ayari,H.P.Simao,Approximate dynamic programming for high dimensional resource allocation problems,Proceedings of theInternationalJointConferenceonNeuralNetworks(IJCNN)1-5(2005)2989-2994.

[22]R.Luus,Iterative Dynamic Programming,Chapman&Hall/CRC Press,New York, 2000.

[23]M.H.Chang,Y.C.Park,T.Lee,Iterative dynamic programming of optimal control problem using a new global optimization technique,Proceedings of the 8th International Symposium on Process Systems Engineering,2003,pp.416-421.

[24]R.Luus,Piecewise linear continuous optimal control by iterative dynamic programming,Ind.Eng.Chem.Res.32(5)(1993)859-865.

[25]L.Gr?ne,An adaptive grid scheme for the discrete Hamilton-Jacobi-Bellman equation,Numer.Math.75(3)(1997)319-337.

[26]L.Gr?ne,W.Semmler,Using dynamicprogramming with adaptive grid scheme for optimal control problems in economics,J.Econ.Dyn.Control.28(12)(2004)2427-2456. [27]S.Patra,S.K.Goswami,B.Goswami,Fuzzy and simulated annealing based dynamic programming for the unit commitment problem,Expert Syst.Appl.36(3)(2009) 5081-5086.

[28]H.P.Christos,S.Kenneth,Combinatorial Optimization:Algorithms and Complexity, Princeton Press,New Jersey,1982.

29 January 2013

☆Supported by the National Basic Research Program of China(2012CB720500)and the NationalHighTechnologyResearchandDevelopmentProgramofChina(2013AA040702).

*Corresponding author.

E-mail address:huangdx@tsinghua.edu.cn(D.Huang).

http://dx.doi.org/10.1016/j.cjche.2014.09.002

1004-9541/?2014 The Chemical Industry and Engineering Society of China,and Chemical Industry Press.All rights reserved.

Received in revised form 5 June 2013

Accepted 30 June 2013

Available online 6 September 2014