Product Dynamic Transitions Using a Derivative-Free Optimization

Jul 15, 2016 - *(A.F.-T.) E-mail: [email protected]. .... In this work we use a variant of the DFO trust-region algorithm, ..... This is one o...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/IECR

Product Dynamic Transitions Using a Derivative-Free Optimization Trust-Region Approach Israel Negrellos-Ortiz Department of Chemical Engineering, Universidad Iberoamericana, Prolongación Paseo de la Reforma 880, México D.F. 01219, México

Antonio Flores-Tlacuahuac* Escuela de Ingenieria y Ciencias, Tecnologico de Monterrey, Campus Monterrey, Ave. Eugenio Garza Sada 2501, Monterrey, Nuevo León, México 64849

Miguel Angel Gutiérrez-Limón Department of Energy, Universidad Autónoma Metropolitana-Azcapotzalco, Av. San Pablo 180, C.P. 02200 México D.F., México ABSTRACT: Traditionally the optimization of processing systems has relied on the availability of an explicit model together with the corresponding gradient information. However, there are some practical scenarios such as (a) nondifferentiable systems, (b) physical experimental systems, (c) simulation environments, and (d) reduced order systems where such a model and its gradient are not available. Under these scenarios the deployment of derivative-free optimization strategies provides an alternative manner to cope with the optimization of such systems. In particular, in this work we deploy a derivative-free optimization trust region approach to deal with the product dynamic optimization problem of processing systems. To this aim, we use a closed-loop model predictive control strategy where the system to be optimized is embedded in a black-box dynamic simulation environment. The results demonstrate that black-box dynamic models can be dynamically optimized assuming that the number of decision variables is not large. The firstprinciples dynamic model of a binary distillation column embedded in the ASPEN dynamic simulation environment was deployed as our black-box dynamic model, to demonstrate the advantages of solving product dynamic transition problems when an explicit model of the dynamic model and/or its gradient information are not available.

1. INTRODUCTION During start-up, shut-down, or normal operation, optimal dynamic transitions among products are normally required to achieve an operation target in the best possible way. Such a target can be related to desired values of processing variables such as temperature, composition, flow rates, etc. Therefore, in these sort of problems the aim consists in computing the time domain values of the control variables u(t) such that the system response y(t) attains a desired value embedded in an objective function Ω(y, u) which can be either maximized or minimized. Optimal product transitions can be systematically calculated by setting the product transition problem as a dynamic optimization problem.1−3 Several optimization techniques for handling the solution of dynamic optimization problems can be found elsewhere.4 All these optimization algorithms assume that the gradient of the objective function and related constraints are somehow available. When this is the case, the user has access to a wide variety of efficient and powerful algorithms © XXXX American Chemical Society

which compute either local or global optimal solutions. Nowadays, when the explicit form of the optimization formulation is available, accurate derivatives are normally computed by automatic differentiation techniques,5 although derivatives can also still be approximated by finite differences approaches. However, this gradient evaluation approach is becoming obsolete in commercial and academic optimization software due to its disadvantages especially when dealing with noisy and expensive to evaluate functions. On the other hand, there are some practical optimization problems where derivatives are not available or are not reliable making the use of gradient-based optimization techniques impossible or undesirable. Such situations have to do either Received: January 19, 2016 Revised: July 10, 2016 Accepted: July 15, 2016

A

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

gradient information is not available or reliable. Specifically, the DFO model predictive control of a benzene−toluene conventional distillation column is addressed. We consider the onepoint composition control of the distillate and bottoms streams. Moreover, the more challenging dual composition control of both distillate and bottoms streams is also considered. The results are compared against a PI control structure. Presently, one important drawback of derivative free techniques is the limitation on the number of optimization variables that can be handled. The outline of this work is as follows. In Section 2 an updated literature review about optimal product transitions and applications of DFO techniques is provided; in Section 3 we explain the ways in which optimal product dynamic transitions were addressed; Section 4 contains the definition of the problem solved in this work; and some case studies are given in Section 5. Conclusions and future work are highlighted in Section 6.

with real industrial or lab-scale processes or with huge models embedded in simulation platforms. As far as processing equipment is concerned, gradient information will not be available, and for safety, cost, or operational reasons such equipment normally is not used for this aim. The situation is similar when dealing with simulation codes. Depending upon the software, the user has access only to the binary code or to limited portions of the code. Therefore, most of the time finite differences are the only way to obtain gradient approximation. However, expensive computer simulations and noisy functions seem to be the two major reasons to avoid gradient approximation along this way. Moreover, gradient approximation by finite differences techniques is not a straightforward task since it requires making complicated choices such as the proper size of perturbations for reliable derivative approximations. When access to the simulation code is granted, there exists the chance of obtaining gradient approximation by adjoint techniques.6,7 Some of the earliest optimization algorithms did not require gradient information for determining an optimal solution.8,9 However, these algorithms lacked a theoretical basis and were intended for small scale systems. Most of these algorithms were abandoned in favor of gradient-based algorithms, since the use of derivative information about the objective function and constraints with respect to the decision variables can dramatically improve the efficiency of optimization algorithms. However, in engineering practice, it is common to find optimization problems where derivatives are not available or are not reliable. Even under these circumstances the optimization of such systems can be desirable. Derivative-free optimization (DFO) algorithms have been proposed to deal with optimization problems with no gradient information.10,11 DFO methods can be broadly classified into deterministic and stochastic algorithms. Some of the deterministic DFO methods are based on well-known unconstrained optimization techniques such as line search and trust region methods.12,13 However, some other variations of deterministic DFO methods have also been proposed such as spatial search, genetic algorithms, and hit-and-run methods.14 On the other hand, stochastic DFO methods lack theoretical basis but tend to be easier to use. Such methods include genetic and particle swarm algorithms.15,16 It should be remarked that there are some theoretical convergence proofs for deterministic DFO methods.11 In this work we use a variant of the DFO trust-region algorithm, as implemented in the BOBYQA code,17 to handle the dynamic optimization of chemical processes during product transitions. Optimal product transitions imply the computation of control actions which lead to efficient and optimal transitions among plant products.2 Commonly, during product transitions, minimum transition time and/or minimum off-spec product should be achieved. However, other objective functions or conflicting objective functions can also be considered.18 To compare the performance of DFO algorithms we run the dynamic optimization of some systems when a dynamic mathematical model is available. Then, the same problem is solved (a) using a penalty function approach to move the nonlinear constraints into the objective function and (b) assuming that the dynamic model is not available. Moreover, using a nonlinear model predictive control approach, we also solve a product dynamic transition problem where the dynamic model is embedded into a process simulation software. This is the practical scenario where probably DFO methods should find their most interesting and valuable application: the optimization of processing systems embedded in large or complex computer codes where

2. LITERATURE REVIEW In this section we provide a literature review summary on the two main topics addressed in the present work: product transitions and applications of derivative-free optimization techniques. Product Transitions. In this work product transitions refer to a kind of optimization problem where the aim is to compute the time value of the manipulated variables to drive a dynamic system from an initial to a target point in an optimal manner.19 Such problems are common in the processing industry. In particular, in the polymerization industry they are referred to as grade transition problems.2 Such product transition problems can be off-line computed, requiring then a control system to perform the closed-loop tracking of the optimal value of the manipulated variables.20 On the other hand, product transition can also be online computed using a control strategy such as Model Predictive Control21 to take care of the impact of disturbances and to a some extent to address the effect of model uncertainty on the transition problem. In this section we provide a short review of some of the more relevant works addressing product transition problems. McAuley and McGregor22 deployed an open-loop optimal control approach to address the grade transition problem in polymerization reaction systems. The authors used density and melt index measurements to quantify the quality of the product. The impact of modeling errors on the transition trajectories was also addressed. Using heuristic transition trajectories as operating policies Debling et al.23 analyzed the open-loop dynamic response of a polyolefin production system. Moreover, since the work only addressed the calculation of open-loop policies, their closed-loop implementation remained as an issue. Cervantes et al.24 were among the first authors to consider the grade transition problem of a complex large-scale polymerization system. The authors addressed a polyethylene tubular reactor optimal grade transition problem. The dynamic model turns out to be highly nonlinear. Moreover, the model comprises a distributed system. For addressing the efficient solution of the optimization problem, finite differences were used for handling the discretization of the spatial behavior, while orthogonal collocation on finite elements took care of the discretization of the time domain behavior. Moreover, to improve the numerical solution of the underlying optimization problem, the authors also exploited the mathematical structure of the problem. Chatzidoukas et al.25 addressed the problem of determining the optimal transition trajectories and the structure of the control system used for enforcing such optimal control trajectories. To solve this problem the authors formulated and solved a B

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Because the Quasi-Newton approach is based on the finitedifferences numerical approximation of the gradient, the authors concluded that this optimization approach is highly sensitive to numerical noise leading to slow convergence. Of course, the selection of step size for the proper computation of the gradient turns out be an additional problem. Very small step sizes can give rise to poor performance or convergence problems since some elements of the gradient can be null. Moreover, the number of simulation runs for computing the gradient can be large. To cope with issue, the authors used a DFO approach which turned out to be a simple and efficient way of addressing the optimization of the approached energy systems. Although both optimization approaches achieved the same optimal solution, after proper tuning of each approach, they concluded that the CPU time required by the DFO approach was smaller than the CPU demanded by the Quasi-Newton optimization approach. Schafer et al.33 and Ugur et al.34 considered the optimization of the configuration of a stirrer for mixing purposes. They used a rigorous computational fluid dynamics (CFD) approach for the modeling and simulation of a tank where mixing operations took place. Using the CFD results as a black box system, they compared the performance of two derivative free optimization methods (DFO and CONDOR). The objective function consisted in the minimization of the dimensionless Newton number whereas the decision variables were the disk thickness, the bottom clearance, and the baffle length. The authors concluded that in one of the addressed examples both optimizers led to different optimal local solutions. The reason for this behavior has to do with the way both solvers build the quadratic model representing the objective function: DFO deploys Newton polynomials whereas CONDOR uses Lagrange polynomials. The authors also stated that for larger number of decision variables (20 or more) the computational load required for both solvers tend to be excessive. Korkel et al.35 compared the performance of derivative based vs DFO methods for handling nonlinear optimum experimental design problems using a simulated experiment to get experimental information. They were interested in the design of experiments for computing control actions which led to statistically reliable parameter estimations with the aim of identifying dynamic mathematical models. First, they solved the problem assuming that model derivatives were available using as objective function the minimization of a function closely related to minimal statistical uncertainty (i.e., function of the variancecovariance matrix). Then, in an inner optimization loop, they considered the parameter estimation problem. The problem was then solved using a standard successive quadratic programming (SQP) approach. On the other hand, the same conceptual problem (i.e., reliable estimation of model parameters) was formulated and solved using a DFO methodology. To this end, the authors used a multidirectional search method proposed by Torczon.36 This is a DFO method which can handle only simple bounds on decision variables in addition to the objective function. To compare the performance of both approaches the same objective function (i.e., a given function of the variancecovariance matrix) was used. Although the authors reported similar results using both approaches the computational load for solving the DFO problem was significantly larger since in this case the number of function evaluations was large. The perceived advantage of the DFO approach was the simplicity of its use. It remains as a challenge to apply the DFO framework

mixed-integer dynamic optimization problem for the computation of both continuous and binary decision variables in a dynamic processing framework. Nystrom et al.26 were among the first authors to consider the interplay between scheduling and dynamic optimization problems in polymerization reaction systems. To address this problem the authors proposed a decomposition approach based on the construction of a master and primal problem. The iterative solution of those problems provided the optimal production schedule and optimal dynamic transition trajectories. Asteasuain et al.27 addressed the problem related to the optimal process design of a styrene polymerization system together with the design of the control system. The reactor design problem included the equipment design, selection of the initiator, and computation of the steady-state conditions and transition trajectories, while the design of the control system included pairings of the control structure and the tuning of the PI controllers. Since optimal transition policies were determined minimizing transition time and offspecification polymer and because these two objectives are in conflict, a multiobjective optimization framework was used to address the optimal solution of the underling mixed-integer dynamic optimization problem. Prata et al.28 also considered the production scheduling and dynamic optimization problem. They proposed a logical mixed-integer dynamic optimization approach to deal with this problem and applied their approach to an industrial polymerization system. Terrazas-Moreno et al.29 addressed the problem of the simultaneous design, scheduling, and optimal control of a continuous methyl-methacrylate polymerization reactor. The resulting optimization problem was cast as a mixed-integer dynamic optimization problem which was discretized and solved as a mixed-integer optimization problem. Moreover, this was one of the first works to consider process uncertainty and dealing with this problem as a stochastic optimization problem using a two-stage optimization formulation. Recently, Weng et al.30 considered the dynamic optimization of a high-density polyethylene process using the molecular weight distribution as the quality index for driving optimal control calculations. Dynamic optimization based on the full dynamic molecular weight distribution turns out to be a computationally expensive problem. To deal with this issue the authors proposed a combination of a small scale problem based on the moments method for representing the molecular weight distribution and a steady-state molecular weight distribution method. However, it remains as a challenge to consider the full dynamic molecular weight distribution as a quality index in grade transition problems. Patil et al.31 proposed a novel methodology for assessing design, scheduling, and control problems under process disturbances and uncertainty conditions in multiproduct production systems. Process disturbances were characterized as sinusoidal signals, whereas uncertainty was incorporated using a multiscenario approach and assigning a probability to each scenario. The proposed strategy was applied to two reaction systems. Applications of Derivative-Free Optimization. Van der Lee et al.32 carried out the optimization of energy production in Rankine-like power cycles. They considered as the objective function the thermodynamic efficiency of the cycle and deployed several energy efficiency case studies which resulted in different numbers and types of decision variables. Simple lower and upper box constraints on the decision variables were considered. To perform the optimization of the energy system they compared the performance of a Quasi-Newton gradientbased approach and a DFO method proposed by Conn et.al.10 C

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

objective function, different sets of decision variables were used. Simple bounds on the decision variables were also included. For solving the DFO problem the BOBYQA algorithm17 was deployed. The authors reported a quick convergence rate of the algorithm in this experimental setting and improved results of the addressed compounds related to the target properties. This work clearly demonstrates that a combination of scientific knowledge and expertise with numerical DFO algorithms provides an efficient, systematic, and reliable way to improve the chemical synthesis of compounds by just testing a small set of synthesis conditions. Forouzanfar and Reynolds42 addressed the optimal oil wellplacement problem using DFO methods. Modeling the behavior of oil wells turns out to be a complex and computationally demanding problem. To this aim, the authors deployed a commercial simulator. In their approach the authors used as objective function the maximization the net present value (NPV) of production in a reservoir. The decision variables had to do with good placement of parameters such as location, length, and trajectory of each well. Simple bounds on the decision variables were also included. For solving the optimization problem the authors considered their own implementation of the BOBYQA algorithm as reported by Powell17 and compared both the quality of the solution and the performance of that algorithm against a genetic algorithm. They found that, for the addressed examples, DFO methods required fewer computer simulations to achieve the optimum solution. However, the solutions obtained using both optimization algorithms tended to be similar but not equal. The authors also reported that different optimal solutions were obtained by starting the DFO algorithm from different initial points. This behavior reflects the fact that the problem features several local optimal solutions and that when converged DFO methods can only achieve such a local optimal solution. Asadollahi et al43 compared the performance and reliability of four DFO strategies while approaching the production optimization of oil wells. The objective function consisted of the maximization of the NPV while water-cuts for different producers were taken as the decision variables. A commercial simulator was used as a black box to evaluate the production of the wells and to provide input information for solving the optimization problems. The addressed DFO methods were as follows: pattern search Hooke−Jeeves, generalized pattern search, Nelder−Mead, and a Line search method. They also compared the performance of those DFO methods against a sequential quadratic gradient based method (SQP), where numerical derivatives were calculated by finite differences. The convergence rate strongly depended upon the initial guesses of the decision variables. In addition, in the case of the SQP method, the results also depended upon the discretization size for the numerical evaluation of derivatives. Globally, the DFO Line search method provided the best results in terms of the number of simulations used to converge the optimization problems. Although, in most of the cases, all the addressed algorithms were eventually able to achieve similar values of the objective function but at larger computational cost. Ciaurri et al.44 deployed a DFO approach for handling the optimization of production wells for generally nonlinear constrained optimization problems. This is one of the first reported attempts for including more general nonlinear constraints in a DFO setting where normally simple box bounds are dealt with. They reported the use of the penalty function and filter methods13 for efficient constraint handling. The authors compared the performance of two DFO methods (generalized pattern search36

for optimum nonlinear experiment design to a real or laboratory scale system. Kramer et al.37 provide a quick overview of both computer algorithms and practical applications of derivative-free optimization methods. They classify derivative-free optimization methods as based on pattern search methods (i.e., generalized pattern search or direct search), methods which approximate the objective function and then use trust-region methods for handling the solution of the optimization problem, and optimization methods based on heuristics and stochastic components such as evolutionary, genetic, and particle-swarm methods. Some practical recommendations about how to handle the presence of constraints in a derivative-free optimization approach are also provided. The authors also mention that future work on the theory and applications of DFO in dynamic optimization, mixed-integer nonlinear programming, and the optimization of uncertain systems is required. Giuliani and Camponogara38 used DFO methods to address the optimization of gas-lifted oil fields. The authors deployed an augmented Lagrangian approach for dealing with constraints which were originally stated as bounds on decision variables. However, for dealing which constraints of the type Ax − b ≤ 0, they just stated that only solutions which made feasible this set of linear constraints would be accepted. After inequality constraints were embedded into the Lagrangian function, the optimization problem was solved using a trust-region algorithm proposed by Conn et al.39 The authors also compared their DFO trust-region method results against a mixed-integer liner programming (MILP) formulation which deployed piece-wise linear approximation of the optimization model. They concluded that the DFO trust-region method led to better optimal solutions and reported infeasible MILP solution for two of the addressed production scenarios. This has to do with the fact that in those cases the piece-wise linear model was not able to represent the model behavior. This behavior was somehow expected, since the DFO trust-region method builds a local approximated model at each iteration. Botelho et.al.40 considered the parameter estimation problem of an industrial polymerization reactor. Prior to solving the parameter estimation problem the authors used an identifiability analysis to decide the set of parameters to be identified. Once that decision was taken, a DFO nonlinear parameter estimation problem was solved. For handling the numerical solution of the DFO the authors deployed the BOBYQA optimization software.17 Rios and Sahinidis14 published a full and exhaustive comparison about the performance of free-derivative optimization algorithms. They compared 22 DFO algorithms using 502 problems. This is practically the first systematic work aiming to compare the performance of DFO algorithms to gain insight about the merits of different DFO algorithms. They came to the conclusion that although important advances (especially in model based DFO algorithms and in global convergence proofs) have been reported both in the theoretical and numerical aspects, more work is needed in DFO especially in extending DFO methods to deal with generally nonlinear constrained problems. Schilling et al.41 provided one of the few reported results of application of DFO to an experimental system aimed to the synthesis of new inorganic−organic compounds with target properties. The authors deployed a set of objective functions related to increasing crystallinity, inhibition of crystallization, increasing crystal size, and tuning of the particle size of different inorganic−organic compounds. Depending upon the deployed D

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research and Hooke−Jeeves direct search8) against a heuristic genetic optimization and a gradient-based sequential quadratic method. They reported good performance of the DFO methods when dealing with constraint handling and came to the conclusion that the generalized pattern method augmented with a filter provided the best results when derivatives are not available. However, nowadays one of the main drawbacks of DFO methods lies in the fact that normally only small scale problems can be dealt with. Koller and Ulbrich45 addressed the production of metal sheets by the hydroforming process as an optimal control problem. They used a commercial simulator to provide a dynamic model for describing the spatial and time variation of the shape of the metal sheets when subject to variation in the fluid pressure which was considered to be the control variable, whereas the objective function was taken as the volume difference between the sheet metal and the die. They compared the performance of the Nealder−Mead method and BOBYQA17 approaches for solving the underlying optimal control problem and came to the conclusion that the BOBYQA solver outperformed the Nealder−Mead for solving the addressed optimal control problem in terms of the value of the objective function and the number of simulation runs required to solve the optimization problem. Dowling et al.46 approached the dynamic optimization of a complex pressure−swing absorption (PSA) system aimed to CO2 capture. The dynamic model consisted of a series of mass and energy conservation equations and some constitutive relationships representing both time and one spatial coordinate system behavior. High CO2 capture was one the objective functions, and several control variables were also selected. Because of the complexity of the mathematical model, a DFO approach17 was used to address the underlying optimal control problem. The authors also used a gradient-based optimization algorithm, embedded in the IPOPT software,47 to deal with the optimization problem. Because they used the solution of the mathematical model as a black box, proper gradient information was calculated by an adjoint approach. Although the authors concluded that the gradient-based optimization algorithm outperformed the DFO strategy, from Tables 10, 11, and 12 of ref 46 it is clear that this was done at the cost of larger CPU times. Fisher and Jiang48 deployed four different DFO methods for the evaluation of kinetic parameters in a combustion system. The methods considered by the authors were as follows: adaptive random search and a genetic algorithm and two deterministic algorithms coded in the programs CONDOR49 and BOBYQA.17 The authors concluded that BOBYQA was the optimization algorithm which provided the best model parameter fitting to reproduce the provided pseudoexperimental information.

problem in systems where a dynamic model is available. Then, the typical constrained dynamic optimization problem was transformed into an unconstrained optimization problem by moving all the sets of constraints, comprising the dynamic model and related additional constraints, into the objective function and penalizing such a set of constraints. For doing so, the set of constraints related to the dynamic model were approximated using the transcription approach.4 The resulting unconstrained penalized optimization problems were then solved using (a) typical nonlinear optimization solvers, available in the GAMS optimization environment,50 and (b) the BOBYQA implementation of a DFO solver.17 As previously stated, we clearly recognize that this is not the real scenario where DFO strategies should be used, since when the dynamic mathematical model and related derivatives are available powerful and efficient gradient optimization techniques13 should be deployed. Our aim was just to compare the performance of a typical DFO algorithm implementation (i.e., BOBYQA) against a straightforward solution using a nonlinear programming solver. Finally, it should be remarked that product transitions computed deploying the above-described procedure are open-loop or off-line control policies. Afterward, a control systems would be needed to implement such control policies in a closed-loop scenario.20 • Model predictive control using black box dynamic modeling. DFO techniques have been used for approaching either linear or nonlinear model predictive control of processing systems.51 Model predictive control (MPC) is a closed-loop control technique well suited for the control of nonlinear systems using a receding horizon approach.21 In this control environment a set of values of the input or manipulated variables (i.e., u1, ..., uN) are computed such that the output or controlled variables (i.e., x1, ..., xM) are forced toward target or desired values, where N and M stand for the length of the control and prediction horizons, respectively. Commonly, only the first control action u1 is applied to the controlled system. When a new measurement becomes available, the control scenario is moved ahead by a time equivalent to the sampling time and the procedure is repeated until the target values are obtained. When a dynamic model or derivatives of the underlying systems are not available, a kind of DFO algorithms, the so-called genetic algorithms, have been applied for addressing the model predictive control problem of such systems.52,53 In this work we deploy a kind of DFO algorithm based on using trust-region optimization techniques, which are numerical techniques aimed to the efficient optimization of unconstrained systems12,13 featuring only a given objective function. In the case of DFO trust-region methods the objective function is approximated by a quadratic function.11 Such an approximation of the objective function is built using only input and output information. Hence, using DFO techniques no information, in the form of an explicit dynamic model or their related derivatives, is required for the dynamic optimization of processing systems. We only require (a) input/output information on the controlled system and (b) solving the dynamic optimization problem by approximating the objective function using DFO trust-region

3. DFO APPROACHES FOR PRODUCT TRANSITIONS In this section we provide details about the method of using DFO techniques for handling MPC problems and a short explanation about the particular DFO algorithm used in this work. • Unconstrained optimization problem with penalized dynamic model. Even when DFO strategies are aimed at systems where no derivative information is available, in a first attempt to try to optimize such a kind of systems and to realize the quality of DFO solutions obtained from typical software implementations, we have assessed the product transition E

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

to the target operating point in an optimal manner. Since this normally implies using an explicit dynamic model for computing such control actions and we assume such a model is not available, we need to implement an MPC algorithm in a different manner when a derivative-free optimization approach is deployed. This means that at any time a number of control actions (given by the prediction horizon) were computed, but in the spirit of the MPC strategy, only the first one of those control actions was implemented, discarding the remaining ones. When a new measurement became available the solution of the closed-loop optimal control problem represented by eqs 4−8 was repeated and so on until target conditions were reached. We would like to highlight that nonlinear dynamic process simulation models of the addressed systems (in particular of the distillation column) were used to perform the optimization calculations. Particularly, for the fourth case study, a DFO nonlinear model predictive control approach was used for the online computation of the control actions. • Implementation of the DFO trust-region method. BOBYQA (Bound Optimization BY Quadratic Approximation) is a derivative-free optimization code, proposed by Powell (2009),17 which approximates the behavior of the objective function deploying a quadratic interpolation polynomial. This optimization approach finds the minimum of an objective function J(u) subject to simple constrains:

optimization techniques. Assuming that the system response would be available in the form of a black-box model, the remaining issue consists in formulating an objective function which mimics the model predictive control approach used when a dynamic model is available. After testing several forms of objective function for product transitions, we came to the conclusion that the following form of objective function is suitable for performing dynamic product transitions: N

Ω=

∑ [αx(xi − x f )2 + αu(ui − u f )2 ] i=1 N

N

+ λx ∑ [Δxi] + λu ∑ [Δui]2 2

i=1

(1)

i=1

where the superscript f stands for final, desired or target values, αx and αu are weights deployed to penalyze for deviations from the target values, and λx and λu are penalties impossed on variations of the magnitude of the states and the manipulated variables, respectively. Actually, the form of the objective function defined in eq 1 was first proposed in ref 51. The first two terms in eq 1 are used to drive the states x and manipulated variable u toward their target values. The aim of the last two terms in eq 1 is to penalize for large changes in both the states and manipulated variable values. This reflects our desire of obtaining smooth control actions and in consequence smooth system response as well. Note that those changes are defined as follows: Δxi = xi − xi − 1 (2) Δui = ui − ui − 1

min J(u)

where u is an n-dimensional decision vector. At each iteration k, a quadratic polynomial is used to approximate the behavior of the objective function:

(3)

therefore at the beginning of the product transition (i.e., i = 1) x0 and u0 stand for the initial values of the states and manipulated variables, respectively. It should be remarked that in eq 1 we have made equivalent the size of the control and prediction horizons. However, the Ω objective function can clearly be cast for different values of N and M as well. In summary, the MPC computation of optimal product dynamic transitions using a DFO approach can be cast as follows:

Q (uik) = J(uik), i = 1, ..., m

The quadratic approximation model Qk (u) can assume the following form: Q k(u) = ak + (u − u0k)T f k +

min Ω = x, u

∑ [αx(xi − x

note that in this form ak is a scalar, f k is the Jacobian, and Fk stands for the Hessian matrix. Also note that the term uk0 stands for the approximation of the decision vector at iteration k. To solve the derivative-free optimization problem, a truncated gradient trust region method is deployed. The optimization problem is run to optimize the quadratic approximation model of the objective function. At each iteration the trust region radius is computed. The optimization problem ends when the trust-region radius becomes smaller than a given target value. More details about the implementation of the BOBYQA algorithm can be found elsewhere.17

f 2

) + αu(ui − u ) ]

i=1 N

N

+ λx ∑ [Δxi] + λu ∑ [Δui]2 2

i=1

i=1

s.t. x l ≤ xi ≤ x u ,

i = 1, ..., N (5)

ul ≤ ui ≤ uu ,

i = 1, ..., N (6)

(4)

Δx l ≤ Δxi ≤ Δx u , i = 1, ..., N (7) Δul ≤ Δui ≤ Δuu , i = 1, ..., N (8)

1 (u − u0k)T 2

F k(u − u0k)

N f 2

such that ulow ≤ u ≤ uup

(5)

where the l and u superscripts stand for lower and upper bounds, respectively. Because this is a derivative-free optimization approach, no equations related to the deployment of the MPC algorithm should be used. Solving an MPC problem requires the online or feedback computation of control actions able to drive the system

4. PROBLEM DEFINITION The problem to be solved in this work can be described as follows. Given • An objective function representing the target grade or product plant transitions. F

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research • A black box system (either industrial or laboratory system or computer simulation) from which only input and output information about the dynamic performance of such a system is recorded. • A set of simple box bound constraints on the system states. • We assume that numerical values of the system derivatives of the objective function are not available. Then the problem consists of computing the time value of the set of control actions u(t) such that the objective function attains its optimal value, meaning that the optimal solution of the dynamic product transition problem can take the system from an initial x0(t) to a final steady-state xf(t) in the best possible way.

the advantages offered by DFO for the dynamic optimization of systems when an explicit dynamic model is not available. Fundamental Differential Equation. The first case study deals with the fundamental differential equation which reads as follows: dy = −ay(t ) + bu(t ) dt

(9)

where y is the state variable, t is the independent variable, u is a control variable, and a, b are constant with a = b = 1. Hence, the penalized state transition objective function, which includes the discretization of the dynamic mathematical model, reads min Ω =

5. CASE STUDIES In this section we carried out the optimal dynamic product transition of several processing systems under different modeling assumptions. The first three examples deal with dynamic systems where a dynamic model is available, whereas a binary distillation column, deployed as the fourth example, is used to demonstrate

1 2

∫0

1

(y 2 + u 2) dt + μF(y , y )̇

(10)

where F(y, ẏ) stands for a function which represents the discretized dynamic model as a function of the states (y) and states derivates (ẏ) and μ is a penalty parameter. Figure 1 displays the optimal dynamic transtions as computed by the CONOPT solver embedded in the GAMS optimization environment50 and using a Fortran implementation of the BOBYQA

Figure 1. Comparison of optimal dynamic optimization results for the fundamental ordinary differential equation. The optimization problem was solved using the transcription approach for system discretization and moving the discretized dynamic system into the objective function using the CONOPT and DFO BOBYQA solvers. (a and b) For 10 finite elements and μ = 30 and (c and d) for 5 finite elements and μ = 100. G

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research DFO solver.17 The vectors of initial [y0, u0]and target [yf, uf ] transition conditions are given by [1, −0.375] and [0.3, 0], respectively. As shown in Figure 1, the dynamic transition problem in BOBYQA was solved using two different values of the penalty parameter μ. We note an oscillation of the manipulated variable response for small μ values reflected as deviation of the state variable with respect to the solution obtained using the dynamic model as an explicit constraint of the optimization problem (see Figure 1a,b). As expected, as the value of μ is increased the oscillation in the manipulated variables practically disappears and the two optimal dynamic transition solutions look similar (see Figure 1c,d). The values of the objective functions using GAMS and BOBYQA are 0.193 and 1.879 × 10−1, respectively. When an explicit dynamic mathematical model is used for product transitions, the set of model constraints must be met. Moreover, in the transcription approach such a set of constraints must be approximated by discretization of the original continuous model. In this work, the function F(y, ẏ), which appears in eq 10 and alike, represents such a discretized dynamic mathematical model. Although our aim is to solve dynamic product transition problems using a derivative-free optimization approach, we also included the optimal solution of such a problem, when gradient information was available like in the first three case studies, using the transcription approach and moving the discretized dynamic model F(y, ẏ) into the objective function as a penalty term. This would allow us to compare both gradient and derivative-free optimization approaches for solving product transition problems. Of course, when gradient information is available it can be more effective to solve the problem as a fully constrained optimization problem instead of penalizing the objective function. Hicks Reactor. In order to compute product transition trajectories, the CSTR model as proposed by Hicks and Ray19 was used. Because the original parameters set used by these authors did not lead to multiple steady states, some of the parameters values were modified in order to end up with a multiplicity map. In dimensionless form the model is given by dy1 dt dy2 dt

= =

1 − y1

− k10e−N / y2y1

θ yf − y2 θ

− k10e−N / y2y1 − αU (y2 − yc )

Table 1. Parameters Used for the Multiple Steady States Hicks Reactor

∫0

t

description

θ J cf α Tf k10 Tc N

20 100 27.6 1.95 × 10−4 300 300 290 5

residence time −ΔH/(ρCp) feed concentration dimensionless heat transfer area feed temperature preexponential factor coolant temperature E1/(RJcf)

representing the discretized dynamic model and μ = 3 × 107 is a penalty parameter. As shown in Figure 2, we assume that a product transition from steady-state B to the steady-state A will be carried out. Optimal product transitions are shown in Figure 3. As shown there, the optimal dynamic transition trajectory using an explicit dynamic model turns out to be different from the solution computed when no dynamic model and related derivatives are used. In fact, from Figure 3a the optimal solution found by BOBYQA is just a linear ramp between the initial and final value of the manipulated variables. On the other hand, the optimal solution found by using an explicit dynamic model is more involved initially implying closing the control valve and then opening it up to an upper bound. Figure 3b,c displays the states dynamic response as a result of imposing the control actions shown in Figure 3a. We note that independently of the dynamic optimization approach, in both cases the minimum transition time turns out to be the same; i.e., product transitions are achieved deploying similar transition time. We can conclude that the DFO optimal dynamic transition policy does not look aggressive in the sense that it avoids the overshooting behavior in the reactor temperature as shown in Figure 3c. The value of the objective functions using GAMS and BOBYQA are 292.411 and 1.223 × 10−3, respectively. Isothermal Biochemical Reactor. In this case study we address the optimal product grade transition problem of an isothermal biochemical reactor proposed by Agrawal et al.54 The dimensionless model reads as follows:

(11)

(12)

{α1(y1(t ) − y1d )2 + α2(y2 (t ) − y2d )2

+ α3(u(t ) − ud)2 } dt + μF(y, y )̇

value

Figure 2. Steady-state multiplicity map of the Hicks reactor. The continuous line stands for stable steady states, whereas the dashed line stands for unstable steady states.

where y1 stands for dimensionless concentration (c/cf), y2 is the dimensionless temperature (Tr/Jcf), yc is the dimensionless coolant temperature (Tc/Jcf), yf is the dimensionless feed temperature (Tf/Jcf), u is the cooling flow rate, and c and Tr are the concentration and reactor temperature, respectively. Table 1 contains the numerical values of the parameters used in this work. In Figure 2 a steady-state multiplicity map is shown. As the objective function the requirement of minimum transition time between the initial steady state and the final desired one is imposed: min Ω =

parameter

(13)

where the superscript d stands for desired or target values and α1 = 1 × 107, α2 = 5 × 103, and α3 = 1 × 10−1 are weights used for stressing the importance of the y1, y2, and u decision variables. Similar to the past case study, the term F(y, y)̇ is a function

dx1 = −x1 + Da(1 − x 2) exp(x 2/γ )x1 dt H

(14)

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 3. Comparison of the product optimal dynamic transition results for the Hicks reactor using an explicit dynamic model solved using the CONOPT and the BOBYQA solvers. Ten finite elements with three internal collocation points were used for the discretization of the dynamic model.

As in the past cases, the term F(x, x)̇ is a function representing the discretized dynamic model and μ = 5 is a penalty parameter. We assume that a product transition from the initial steady state A to the final steady state B (see Figure 4) will be carried out. Optimal product transitions are shown in Figure 5. As shown, this time the optimal product dynamic transition solutions, using the discretized explicit model and without using the dynamic model and related derivatives, are not so different. Initially, as shown in Figure 5a, the DFO optimal solution has a larger undershoot behavior. However, this initial behavior seems to be helpful to avoid the overshooting behavior displayed by the optimal solution which uses the full explicit dynamic model. It should be also noted that the two optimal product dynamic transition strategies feature essentially the same value of the transition time. The value of the objective functions using GAMS and BOBYQA are 8 × 10−3 and 1.48 × 10−1, respectively. Product Transitions Using an ASPEN Simulation Black Box Model. In this section we provide an example of the typical environment in which DFO methods represent an advantage over gradient-based optimization strategies: an explicit dynamic model of the underlying system and or its gradient are not available. We take advantage of the fact that there are some commercial process simulators which have embedded powerful first-principles mathematical models. Commonly, such process

dx 2 = −x 2 + Da(1 − x 2) exp(x 2/γ )(1 + β) dt /(1 + β − x 2)

(15)

where x1 is the normalized cell concentration, x2 is the substrate conversion, and Da is the Damköhler number. Using β = 0.1 and γ = 0.4, we obtain the bifurcation diagram depicted in Figure 4 for the x2 state, while the nominal steady states are shown in Table 2. It should be noted that in this case x1 and x2 are the output variables, whereas Da is the manipulated variable. For comparing the quality of the optimal product transition policies we carry out a product transition from the initial steady state A to the final steady state B. In this case the following minimum transition time objective function between the initial steady state and the final desired one is imposed: min Ω =

∫0

t

{αx1(x1(t ) − x1d)2 + αx2(x 2(t ) − x 2d)2

+ α Da(Da(t ) − Dad)2 } dt + μF(x , x)̇

(16)

where the superscript d stands for desired or target values. For this problem αx1 = 1 × 101, αx2 = 1 × 102, and αDa = 1 × 101. I

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 4. Multiple steady states diagram of the biochemical reactor.

Figure 6 displays the flowsheet of the steady-state simulation of the column deployed to obtain a dynamic model of the binary distillation column using the Pressure Driven mode available in ASPEN PLUS.55 To test the performance of the DFO MPC control scheme the following closed-loop scenarios were examined: (a) one point composition control of the distillate stream, (b) one point composition control of the bottoms stream, and (c) simultaneous dual composition control of both distillate and bottoms streams. To control the benzene mole fraction in the distillate stream we deployed the mass reflux flow rate as the manipulated variable, whereas the toluene mole fraction in the bottoms stream was regulated using the reboiler thermal duty as the manipulated variable. It should be stressed that proportional−integral (PI) primary control loops were also included. Accordingly, the column pressure was regulated using the condenser thermal duty, the condenser level was controlled using the distillate flow rate, and the reboiler level was controlled deploying the bottoms flow rate. Table 3 shows the values of the Ziegler and Nichols gains and integral times for each one of the primary control loops. Moreover, to compare the performance of the DFO MPC control system, two PI composition controllers were also tuned. Such controller tuning values were found deploying the Autotune tool available in Aspen Dynamics. As previously discussed in Section 3, for approaching the closed-loop calculation of optimal dynamic product transitions we use a DFO MPC formulation for the objective function. The formulation of the objective function for one point composition control reads as follows:

Table 2. Nominal Steady States for the Bioreactor System x1 x2 Da

A

B

C

D

0.1683 0.8926 1

0.2685 0.6343 0.56

0.2335 0.3363 0.65

0.1354 0.1581 0.8

simulators have limited optimization capabilities being most of them related to steady-state process simulation scenarios. In these simulation environments dynamic optimization algorithms are not available and are difficult to implement or are limited to a kind of linear MPC technique. This is true especially for sequential process simulators, since in equationoriented simulation environments, because of the simulator structure, dynamic optimization facilities can be available. In any case, in this case study we will assume that a typical black-box dynamic first-principles model is available. In this black box system we provide values of the manipulated variables and record the system responses. Using this black box scenario we build and solve a DFO product transition problem capable of taking the dynamic system from an initial to a final target point in an optimal manner using a nonlinear MPC optimization environment. A nonlinear dynamic model of a benzene/toluene binary distillation column was built using the ASPEN Dynamics simulation environment. The equimolar feedstream is composed of a total flow rate of 1000 kmol/h at 100 °C and 5 bar. The distillation column has 12 trays and a total condenser, and the number of the feed tray is 6 counted down from top to bottom. Moreover, at the condenser the column operates at 1 bar with a 0.1 bar drop pressure per tray. We used the rigorous vapor−liquid equilibrium RadFrac module available in Aspen Plus. In addition, the Chao-Seader method was deployed for vapor−liquid equilibrium calculations. In addition, the values of the reflux ratio and reboiler thermal duty were set at 6.3457 and 26.568 GJ, respectively. Under these processing conditions the benzene mole fraction in the distillate stream is 0.8, whereas the toluene mole fraction in the bottoms stream is 0.8 as well.

N

min Ωi =

∑ [α x(xij − xid)2 + α u(uij − uid)2 ] j=1

N

N

+ α x ∑ [Δxij]2 + α u ∑ [Δuij]2 , j=1

J

j=1

i = [b , t ] (17)

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 5. Comparison of the product optimal dynamic transition results for the biochemical reactor using an explicit dynamic model solved using the CONOPT and the BOBYQA solvers. Twenty finite elements with three internal collocation points were used for the discretization of the dynamic model.

Figure 6. Flowsheet of the benzene−toluene distillation column. Discharge pressure in PUMP-D and PUMP-B pumps is 6 bar. Output pressure in CVD and CVB control valves is 3 bar. Moreover, a 5 min residence time, for sizing the accumulator and sump tanks, was set.

where x stands for the mole fraction which is the controlled variable, u stands for the manipulated variable, N is the prediction and control horizons, and the superscript d stands for

desired or target values. Moreover, the index i stands for the number of control loops, which in this case, according to Table 3, can only take the following values: i = b, t (i.e., i = b means control K

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Moreover, the following values of the weighting functions were deployed: αxb = 1, αxt = 1, αub = 1, and αut = 0.1 One-Point Composition Control of the Distillate Stream. In this case study we carried out two product dynamic transitions: (a) from the initial steady-state given by [xb0, R0] = [0.8,26.0923] to the target steady state given by [xbf , Rf ] = [0.9,36.700], and (b) then the product transition in the opposite direction was done. Here, xb stands for benzene mole fraction and R is the mass reflux flow rate [Ton/h], whereas the subscripts 0 and f stand for initial and final values, respectively. We have assumed that online measurements of the benzene mole fraction in the distillate stream will be available every 30 min and that the length of the prediction/control horizon N is 5. In addition, the bottoms composition was not put under control. For implementing both product dynamic transitions we proceeded as follows: (1) first, we initialize the simulation with the steady state conditions, (2) then change to dynamic mode, (3) run the dynamic simulation for 30 min, (4) stop the simulation, and (5) use the values of mass reflux flow rate and benzene concentration in the DFO BOBYQA solver as decision variables and solve the DFO problem; BOBYQA yields the new calculated value of the mass reflux flow rate to be implemented in the simulation environment. Repeat steps 3−5 when new measurements become available until the target steady state is reached. In Figure 7 the optimal product transitions using a nonlinear MPC approach and the BOBYQA DFO trust-region solver are shown. As a note, the DFO BOBYQA solver is capable of performing the requested product dynamic transitions. One of the things to notice from the results displayed in Figure 7 is that product transitions are not symmetric meaning that the dynamic behavior of the benzene/toluene distillation column comprises a nonlinear system. Moreover, a relatively small number of control actions were required to reach the target steady states. In addition, it should be stressed that the dynamic optimization results were obtained without using an explicit dynamic model of the distillation column and/or its gradient, which represents an advantage when a dynamic model is not available or nondifferentiable. The value of the objective function is 5.67 × 10−4.

Table 3. Ziegler−Nichols PI Controller Parameter Values for the Primary (three first control loops) and Secondary Composition Control Loops of the Benzene−Toluene Distillation Column System loop

controlled variable

manipulated variable

1

column pressure

2

condenser level

3

reboiler level

b

benzene mole fraction toluene mole fraction

condenser thermal duty distillate flow rate bottoms flow rate mass reflux flow rate reboiler thermal duty

t

integral time

action

20

12

reverse

20

1000

direct

20

20

direct

gain

24.93 6.135

9.5

reverse

7.5

reverse

of the benzene mole fraction, whereas i = t means control of the toluene mole fraction). In addition, Δx and Δu have similar definitions like in eqs 2 and 3. Finally, αx and αu are weighting functions for the composition and manipulated variable, respectively. Similarly, for the dual composition control of the distillation column the following form of the objective function was deployed: t

min Ω =

N

∑ ∑ [αix(xij − xid)2 + αiu(uij − uid)2 ] i=b j=1

t

+

N

t

∑ αix ∑ [Δxij]2 i=b

+

N

∑ αiu ∑ [Δuij]2

j=1

i=b

j=1

(18)

Independently of the type of composition control, the following bounds on the controlled and manipulated variables were enforced: 0 ≤ xb ≤ 1

(benzene concentration)

0 ≤ xt ≤ 1

(toluene concentration)

26.0923 ≤ ub ≤ 62.957 t

26.56 ≤ u ≤ 41.86

(mass reflux flow rate [Ton/h]) (reboiler duty [GJ/h])

Figure 7. One-point control of the benzene mole fraction in the distillate stream using the DFO-MPC approach. The sampling time was set at 30 min. L

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 8. One-point control of the toluene mole fraction in the bottoms stream using the DFO-MPC approach. The sampling time was set at 30 min.

Figure 9. Two-point control of the benzene mole fraction in the distillate stream and of the toluene mole fraction in the bottoms stream using the DFO-MPC approach. The sampling time was set at 30 min. M

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 10. Comparison of one-point control of the benzene mole fraction in the distillate stream versus a Ziegler−Nichols PI controller. The sampling time was set at 30 min.

Comparison against Proportional−Integral Controllers. In this section we carried out a comparison between the performance of the DFO-MPC controllers and the corresponding closed-loop behavior obtained from simple PI control structures. Intuitively, we expected better closed-loop response from DFO-MPC controllers because they use more information regarding the controlled system. Moreover, the presence of the relatively long measurement delay will create serious problems when deploying PI control systems. Accordingly, in Figure 10 closed-loop results are compared using the DFO-MPC control system and a Ziegler−Nichols PI control system when the one-point composition control of the distillate stream was attempted. It is clear that even for a single-point composition control the PI controller is not able to perform the requested product transition. The main reason for the poor performance of the PI controller has to do with the length of the sampling time. Since the remaining control scenarios displayed similar closed-loop response when deploying PI control system, only the comparison of the closed-loop control of the distillate stream is shown.

One-Point Composition Control of the Bottoms Stream. In this case the mole fraction of toluene in the bottoms stream was closed-loop controlled using the reboiler thermal duty as the manipulated variable. It should be stressed that the mole fraction of the benzene in the distillate stream was not under closed-loop control. Two product transitions were addressed in this case: (a) from the initial steady-state given by [xt0, Q0] = [0.8,26.5681] to the target steady state given by [xtf , Qf ] = [0.9,29.99] and (b) the product transition in the opposite direction. Similarly to the past case, xt stands for toluene mole fraction, and Q is the reboiler thermal duty [GJ/h], whereas the subscrips 0 and f stand for initial and final values, respectively. In Figure 8 the results for the controlled and manipulated variables are displayed. As it can be seen, the results are acceptable in terms of achieving the product transitions in a relatively smooth manner. Globally speaking and as expected, one-point control was not a major issue using the DFO MPC approach. The value of the objective function is 4.08 × 10−4. Dual Composition Control of the Distillate and Bottoms Streams. The more challenging and interesting problem in conventional distillation control refers to the simultaneous control of both distillate and bottoms stream mole fraction. The problem is hard to deal with because of the commonly strong interactions between such control loops and large measurement delays. Because of its multivariable nature and its ability to deal with nonlinear systems, nonlinear MPC controllers are a good option to tackle the dual composition control issue in conventional distillation columns. In Figure 9 closed-loop results for both distillate and bottoms stream are displayed. The initial and final steady states are given as follows: [xb0, xt0, R0, Q0] = [0.8, 0.8, 26.092, 26.5681] and [xbf , xtf , Rf, Qf ] = [0.9, 0.9, 62.657, 41.73], respectively. Similarly to the two past one-point control problems, two product transitions were attempted: (a) from their initial steady-state to the target steady-state and (b) in the opposite direction. As noted from Figure 9, the DFO MPC controller is able to perform the requested product transitions in a relatively small number of control steps. Moreover, the closed-loop response is smooth enough in any transition direction. The value of the objective function is 5.44 × 10−1.

6. CONCLUSIONS There are some practical science and engineering situations where the deployment of optimization tools can bring important benefits measured in terms of economic, health, and environment indicators. Commonly, the development of a good mathematical model is a first requisite to apply formal gradient-based optimization tools. There are some situations where such models can be difficult to obtain because of the involved complexity in building experimental equipment for model verification purposes or reliable model parameters lacking. In other practical cases, the mathematical model can be only available in the form of a large computer program or in a simulation environment. In any of these cases, no explicit model would be available for optimization purposes. Under these scenarios, DFO techniques can provide an effective and efficient way of solving optimization problems for which no gradient information is available. In particular, in this work we have deployed DFO trust-region techniques for addressing the solutions of such problems. When the optimization problem N

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

(23) Debling, J.; Han, G.; Kuijpers, J.; VerBurg, J.; Zacca, J.; Ray, W. Dynamic Modeling of Product Grade Transitions for Olefin Polymerization Processes. AIChE J. 1994, 40, 506−520. (24) Cervantes, A.; Tonelli, S.; Brandolin, A.; Bandoni, J.; Biegler, L. Large-Scale Dynamic Optimization for Grade Transitions in a Low Density Polyethylene Plant. Comput. Chem. Eng. 2002, 26, 227−237. (25) Chatzidoukas, C.; Perkins, J.; Pistikopoulos, E.; Kiparissides, C. Optimal Grade Transition and Selection of Closed-Loop Controllers in a Gas-Phase Olefin Polymerization Fluidized Bed Reactor. Chem. Eng. Sci. 2003, 58, 3643−3658. (26) Nystrom, R. H.; Franke, R.; Harjunkoski, I.; Kroll, A. Production Campaign Planning Including Grade Transition Sequencing and Dynamic Optimization. Comput. Chem. Eng. 2005, 29, 2163−2179. (27) Asteasuain, M.; Bandoni, A.; Sarmoria, C.; Brandolin, A. Simultaneous Process and Control System Design for Grade Transition in Styrene Polymerization. Chem. Eng. Sci. 2006, 61, 3362−3378. (28) Prata, A.; Oldenburg, J.; Kroll, A.; Marquardt, W. Integrated Scheduling and Dynamic Optimization of Grade Transitions for a Continuous Polymerization Reactor. Comput. Chem. Eng. 2008, 32, 463−476. (29) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. Simultaneous Design, Scheduling, and Optimal Control of a MethylMethacrylate Continuous Polymerization Reactor. AIChE J. 2008, 54, 3160−3170. (30) Weng, J.; Shao, Z.; Chen, X.; Gu, X.; Yao, Z.; Feng, L.; Biegler, L. A Novel Strategy for Dynamic Optimization of Grade Transition Processes Based on Molecular Weight Distribution. AIChE J. 2014, 60, 2498−2512. (31) Patil, B.; Maia, E.; Ricardez-Sandoval, L. Integration of Scheduling, Design, and Control of Multi product Chemical Processes under Uncertainty. AIChE J. 2015, 61, 2456−2470. (32) Van der Lee, P.; Terlaky, T.; Woudstra, T. A New Approach to Optimizing Energy Systems. Comput. Methods Appl. Mech. Engrg. 2001, 190, 5297−5310. (33) Schafer, M.; Karasozen, B.; Uludag, Y.; Yapici, K.; Ugur, O. Numerical Method for Optimizing Stirrer Configurations. Comput. Chem. Eng. 2005, 30, 183−190. (34) Ugur, O.; Karasozen, B.; Schafer, M.; Yapici, K. Derivative Free Optimisation Methods for Optimising Stirrer Configurations. European Journal of Operational Research 2008, 191, 855−863. (35) Korkel, S.; Qu, H.; Rucker, G.; Sager, S. Derivative Based vs. Derivative Free Optimization Methods for Nonlinear Optimum Experimental Design. In Current Trends in High Performance Computing and Its Applications; Zhang, W., Tong, W., Chen, Z., Glowinski, R., Eds.; Springer: Berlin, Heidelberg, 2005; pp 339−344. (36) Torczon, V. On the Convergence of Pattern Search Algorithms. SIAM Journal on Optimization 1997, 7, 1−25. (37) Kramer, O.; Echeverria, D.; Koziel, S. Derivative-Free Optimization. In Computational Optimization, Methods and Algorithms; Koziel, S., Yang, X.-S., Eds.; Springer: Berlin, Heidelberg, 2011; pp 61−83. (38) Guilani, C.; Camponogara, E. Derivative-Free Methods Applied to Daily Production Optimization of Gas-Lifted Oil Fields. Comput. Chem. Eng. 2015, 75, 60−64. (39) Conn, A.; Gould, N.; Sartenaer, A.; Toint, P. Convergence Properties of an Augmented Lagrangian Algorithm for Optimization with a Combination of General Equality and l Linear Constraints. SIAM, J. Optim. 1996, 6, 674−703. (40) Botelho, V.; de Souza Moreira, I.; Trierweiler, J. O.; Neumann, G. A.; Farenzena, M. Estimation of Kinetic Parameters of a Polymerization Reactor using Real Data. IFAC Proceedings Volumes 2012, 45, 685−690. (41) Schilling, L.; Niekiel, F.; Stock, N.; Hartke, B. ComputerAssisted Synthesis Optimisation of Inorganic-Organic Hybrid Compounds using the Local Optimisation Algorithm BOBYQA. ChemPlusChem 2014, 79, 863−871. (42) Forouzanfar, F.; Reynolds, A. Well-Placement Optimization using Derivative-Free Method. J. Pet. Sci. Eng. 2013, 109, 96−116.

involves a few decision variables, DFO algorithms perform well. We have shown that, using an MPC framework, the product optimal dynamic transition problem can be efficiently solved using a dynamic model embedded in the ASPEN process simulator. In future work we will address the following problems: (a) processing systems involving uncertainty, (b) nondifferentiable systems, and (c) design scenarios involving discrete decision variables.



AUTHOR INFORMATION

Corresponding Author

*(A.F.-T.) E-mail: antonio.fl[email protected]. Notes

The authors declare no competing financial interest.



REFERENCES

(1) Flores-Tlacuahuac, A.; Moreno, S. T.; Biegler, L. T. Global Optimization of Highly Nonlinear Dynamic Systems. Ind. Eng. Chem. Res. 2008, 47, 2643−2655. (2) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldívar-Guerra, E. Optimal Grade Transitions in the HIPS Polymerization Process. Ind. Eng. Chem. Res. 2006, 45, 6175−6189. (3) Biegler, L. T. An Overview of Simultaneous Strategies for Dynamic Optimization. Chem. Eng. Process. 2007, 46, 1043−1053. (4) Biegler, L. Nonlinear Programming; Society for Industrial and Applied Mathematics: 2010. (5) Griewank, A.; Walther, A. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation; SIAM: 2008. (6) Cao, Y.; Li, S.; Petzold, L.; Serban, R. Adjoint Sensitivity Analysis for Differential-Algebraic Equations: The Adjoint DAE System and Its Numerical Solution. SIAM Journal on Scientific Computing 2003, 24, 1076−1089. (7) Strang, G. Computational Science and Engineering; WellesleyCambridge Press: Wellesley, MA, 2007. (8) Hooke, R.; Jeeves, T. A. Direct Search Solution of Numerical and Statistical Problems. J. Assoc. Comput. Mach. 1961, 8, 212−229. (9) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Computer Journal 1965, 7, 308−313. (10) Conn, A.; Scheinberg, K.; Toint, P. Recent Progress in Unconstrained Nonlinear Optimization without Derivatives. Mathematical Programming 1997, 79, 397−414. (11) Conn, A.; Scheinberg, K.; Vicente, L. Introduction to DerivativeFree Optimization; SIAM: 2009. (12) Dennis, J.; Schnabel, B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; SIAM: 1987. (13) Nocedal, J.; Wrigth, S. Numerical Optimization; Springer: 2006. (14) Rios, L.; Sahinidis, N. Derivative-Free Optimization: A Review of Algorithms and Comparison of Software Implementations. J. Glob. Optim. 2013, 56, 1247−1293. (15) Goldberg, D. Genetic Algorithms in Search, Optimization and Machine Learning; Addison-Wesley: 1989. (16) Clerc, M. Particle Swarm Optimization; Wiley: 2006. (17) Powell, M. J. D. The BOBYQA Algorithm for Bound Constrained Optimisation without Derivatives; Report No. DAMTP 2009/NA06; University of Cambridge: 2009. (18) Morales, P.; Flores-Tlacuahuac, A. Muti-Objective Nonlinear Model Predictive Control of Semibatch Polymerization Reactors. Macromol. React. Eng. 2012, 6, 252−264. (19) Hicks, G. A.; Ray, W. H. Approximation Methods for Optimal Control Synthesis. Can. J. Chem. Eng. 1971, 49, 522−528. (20) Flores-Tlacuahuac, A.; Alvarez, J.; Saldívar-Guerra, E.; Oaxaca, G. Optimal Transition and Robust Control Design for Exothermic Continuous Reactors. AIChE J. 2005, 51, 895−908. (21) Henson, M. A.; Seborg, D. E. Nonlinear Process Control; Prentice-Hall, 1997. (22) McAuley, K.; MacGregor, J. Optimal Grade Transition in a Gas Phase Polyethylene Reactor. AIChE J. 1992, 38, 1564−1576. O

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (43) Asadollahi, M.; Naevdal, G.; Dadashpour, M.; Kleppe, J. Production Optimization using Derivative Free Methods Applied to Brugge Field Case. J. Pet. Sci. Eng. 2014, 114, 22−37. (44) Ciaurri, D. E.; Isebor, O.; Durlofsky, L. Application of Derivative-Free Methodologies to Generally Constrained Oil Production Optimisation Problems. Int. J. Mathematical Modelling and Numerical Optimization 2011, 2, 134−161. (45) Koller, D.; Ulbrich, S. Optimal Control of Hydroforming Processes. PAMM 2011, 11, 795−796. (46) Dowling, A.; Vetukuri, S.; Biegler, L. Optimization Strategies for Pressure Swing Adsorption Cycle Synthesis. AIChE J. 2012, 58, 3777− 3791. (47) Waechter, A. An Interior Point Algorithm for Large-Scale Nonlinear Optimization with Applications in Process Engineering. Ph.D. thesis, Carnegie Mellon University, 2002. (48) Fischer, M.; Jiang, X. Numerical Optimisation for Model Evaluation in Combustion Kinetics. Appl. Energy 2015, 156, 793−803. (49) Vanden Berghen, F.; Bersini, H. CONDOR, a New Parallel, Constrained Extension of Powell’s UOBYQA Algorithm: Experimental Results and Comparison with the DFO Algorithm. Journal of Computational and Applied Mathematics 2005, 181, 157−175. (50) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMS: A User’s Guide; GAMS Development Corporation: 1998; http://www. gams.com. (51) Daehlen, J. S.; Eikrem, G. O.; Johansen, T. A. Nonlinear Model Predictive Control using Trust-Region Derivative-Free Optimization. J. Process Control 2014, 24, 1106−1120. (52) Onnen, C.; Babuška, R.; Kaymak, U.; Sousa, J. M.; Verbruggen, H. B.; Isermann, R. Genetic Algorithms for Optimization in Predictive Control. Control Engineering Practice 1997, 5, 1363−1372. (53) Sarimveis, H.; Bafas, G. Fuzzy Model Predictive Control of Non-Linear Processes using Genetic Algorithms. Fuzzy Sets and Systems 2003, 139, 59−80. (54) Agrawal, P.; Lee, C.; Lim, H. C.; Ramkrishna, D. Theoretical Investigations of Dynamic Behavior of Isothermal Continuous Stirred Tank Biological Reactors. Chem. Eng. Sci. 1982, 37, 453−462. (55) Jana, A. Process Simulation and Control Using Aspen; PHI Learning: 2012.

P

DOI: 10.1021/acs.iecr.6b00268 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX