The Combination of Adaptive Pseudospectral Method and Structure

Article pubs.acs.org/IECR

The Combination of Adaptive Pseudospectral Method and Structure Detection Procedure for Solving Dynamic Optimization Problems with Discontinuous Control Profiles Ping Wang,† Chaohe Yang,*,† and Zhihong Yuan*,‡ †

State Key Laboratory of Heavy Oil Processing, China University of Petroleum, Qingdao, 266580, China Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States

‡

ABSTRACT: To cope with computational challenges on dynamic optimization problems with discontinuous control profiles, a methodology combining the adaptive pseudospectral method with the structure detection procedure is presented. First, the adaptive pseudospectral method, which divides time intervals based on an estimation of approximation error, is utilized to obtain potential sequence and type of arcs composed of discontinuous control profiles. Second, the derived information is exploited by the structure detection procedure to reformulate the original nonsmooth problem as a multistage problem (MSP), where adjacent arcs with the same input type are merged into one to remove unnecessary LGR points for the final approximation. Additionally, a stop criterion by monitoring variations of the Hamiltonian function at discrete points is proposed to estimate whether the computed solution satisfies the necessary conditions of optimality in the sense of Pontryagin’s Minimum Principle. The efficiency of the proposed method is demonstrated through its application to the Williams-Otto semibatch reactor. Finally, elucidated relationships between three important tuning parameters (i.e., the accuracy tolerance ε, the minimum and maximum number of allowable LGR points in each subinterval Nmin and Nmax) and the performance can guide further applications of the proposed methodology to large-scale problems.

1. INTRODUCTION Multiproduct and multipurpose plants comprising batch and continuous processes with frequent product specification changes can be found everywhere in the chemical industry. Optimal operation of such a transient process requires the solution of a dynamic optimization problem to determine the control policy under specific constraints.1 Although extensive research on dynamic optimization has existed, large-scale nonlinear dynamic models development, and efficient solution methods for dynamic optimization problems are still two main barriers for the widespread industrial implementation of dynamic optimization.2,3 Numerical methods for solving dynamic optimization problems can be typically divided into two categories, namely, indirect methods and direct ones. Direct methods, which rely on the discretization of the original dynamic optimization problem and the application of nonlinear programming technique, seem well-suited for solving many practical problems.4 Direct methods are quite broad and cover several different approaches, such as the control vector parametrization (CVP or sequential) approach,5,6 the direct multiple shooting approach,7 and the full discretization (or simultaneous) approach.8 What distinguishes these approaches is the choice of which kind of variables to be discretized and how to approximate dynamic equations in the optimization framework. However, using direct methods to solve a dynamic optimization problem with discontinuous control profiles, which contains a series of so-called arcs and switching points between them, may lead to unsatisfactory results.1,6,9 For example, when the control with lower/upper bound appears linearly in dynamical equations and objective functions, the optimal solution is often characterized by bang−bang type control profiles. © 2014 American Chemical Society

Discontinuities in control profiles may lead to fundamental problems in discretization. In particular, a discretization which is too coarse on a uniform grid or which locally resolves the control profile in the wrong place on a nonuniform grid cannot meet prespecified accuracy requirements. On the other hand, a discretization which is too fine often causes inappropriately high computational cost and might result in robustness problems. In general, the discretization resolution is determined by visually inspecting the obtained solution, which is neither efficient nor precise. With the hope of achieving a trade-off between the solution accuracy and the computational cost,10 several coarse-to-fine policy search algorithms have been proposed. The main idea of these algorithms is to design a sequence of NLPs so that a solution to the final problem is sufficiently close to the solution of the continuous time dynamic optimization problem. For instance, an adaptive CVP algorithm, which can generate problem-adapted discretizations over time spans, has been presented.11 In this algorithm, discretization of the control profile is successively refined based on a wavelet-based analysis of the optimal solution which can be obtained in the preceding optimization step. Therefore, a solution with comparable accuracy can be obtained with less computational effort when compared with a traditional algorithm which relies on an equidistant discretization. In more recent work, the adaptive CVP algorithm has been extended into an automated detection procedure for determining the control structure including the number, type, and Received: Revised: Accepted: Published: 7066

December 7, 2013 April 2, 2014 April 5, 2014 April 5, 2014 dx.doi.org/10.1021/ie404148j | Ind. Eng. Chem. Res. 2014, 53, 7066−7078

Industrial & Engineering Chemistry Research

Article

sequence of arcs as well as the switching times between them.12,13 It has been illustrated that the combination of adaptation and a structure detection algorithm in an iterative manner can further enhance the robustness and efficiency of the solution methods.12 However, seeking the optimal number of switching points will become more complex when higher dimensional problems with respect to the number of control and path constraints are tackled.13 As a major class of direct methods, pseudospectral (also known as global orthogonal collocation) methods have attracted extensive attention over the past several decades.14,15 In a pseudospectral method, state and control variables are approximated using an appropriate set of global trial (basis) functions (e.g., Lagrange or Chebyshev polynomials), and the dynamics are orthogonally collocated (i.e., the collocation points are the roots of an orthogonal polynomial).14 The most well developed pseudospectral methods are the Gauss pseudospectral method (GPM),16,17 the Lobatto pseudospectral method (LPM),18 and the Radau pseudospectral method (RPM),19,20 where the collocation is performed at the Legendre−Gauss (LG) point, the Legendre−Gauss−Lobatto (LGL) point, and the Legendre−Gauss−Radau (LGR) point. All three sets of collocation points arise from a different form of Gaussian quadrature, and state variables are approximated globally by Lagrange polynomials. For smooth optimal control problems, these methods exhibit higher accuracy and efficiency than traditional direct ones,14,21 but their usage in solving nonsmooth problems may cause difficulties because of the prefixed distributions of global orthogonal polynomials and collocation points. Finding a proper distribution of collocation points for pseudospectral methods has been widely investigated. In the knotting method,22 this has been addressed by adding knots to divide the time interval into two or more segments, and therefore the distribution of collocation points would change accordingly. As collocation points are much denser at boundary points, knots should be placed where the solution is subject to sudden changes to achieve desired accuracy. The idea of knotting has been extended by Gong and co-workers23 to propose a refinement strategy which employs a differentiation matrix to identify discontinuities in the solution and uses Gaussian quadrature rules to generate a mesh that is dense near the end point of the time interval of interest. More recently, several adaptive methods have been presented to improve the convergence rate, the solution accuracy, and the applicability of pseudospectral methods, especially for nonsmooth problems. For example, Darby et al.24 have proposed an adaptive pseudospectral method that allows for changes in both the number of subintervals and the degree of the approximating polynomial within a subinterval to achieve specified accuracy. Rao et al.25 have also presented an hp adaptive method where the approximation error is estimated using the difference between an interpolated value of the state and a LGR quadrature approximation to the integral of the dynamics. Taking full advantage of the exponential convergence rate of the Gaussian quadrature collocation method, the estimated error remains computationally tractable when a high-accuracy solution is desired and reduces significantly the number of collocation points required to meet a specified accuracy tolerance compared with the methods of Darby et al.24 Despite their advantages, one potential limitation of these adaptive methods is the accumulation of unnecessary subintervals around discontinuous points since the number of subintervals can only increase in their iterative procedures.24,26

To solve dynamic optimization problems with discontinuous control profiles efficiently, an integrated framework that combines an adaptive pseudospectral method with the structure detection procedure is presented. In order to find the potential sequence and type of arcs in control profiles from numerical solution automatically, the adaptive LGR method is adopted first to iteratively determine the points where optimal control profiles are subject to sudden changes. The derived information is then exploited by the structure detection procedure to reformulate the original nonsmooth problem as a multistage problem (MSP). In the MSP, each stage corresponds to a particular arc obtained in the previous one. Further, adjacent arcs with the same input type will be merged into one to remove unnecessary LGR points for the final approximation. The MSP formulation, which contains the length of each stage as a decision variable, allows discontinuities in the control at the interface between two stages, thus optimal values of switching points can be obtained. The advantage of the pseudospectral method over other direct methods is that it not only provides the state and control profiles but also generates the costate and Hamiltonian trajectories by applying the costate mapping theorem.27,28 This feature provides a practical way to validate the optimality of the resulting solution by monitoring variations of the Hamiltonian function at LGR points for each subinterval. Therefore, the overall solution procedure will terminate as soon as the variation of Hamiltonian function meets the user-defined tolerance. The effectiveness of the proposed methodology is demonstrated by solving the dynamic optimization problem of the Williams−Otto batch reactor.29−31

2. PROBLEM FORMULATIONS Without a loss of generality, consider the following general dynamic optimization problem in Bolza form: min J = ϕ(x(tf )) +

u(t ), t f

∫t

tf

g (x(t ), u(t )) dt

0

(1-a)

s . t . ẋ(t ) = f(x(t ), u(t )), x(t0) = x 0

(1-b)

h(x(t), u(t)) ≤ 0

(1-c)

φ(x(tf )) ≤ 0

(1-d)

where x(t) ∈ R denotes the state variables with initial conditions x0. The process model is formulated as the smooth function f(·). The control variables u(t) ∈ Rnu and the final time tf are decision variables for the optimization problem on the time interval t ∈ [t0,tf]. The problem includes path and terminal constraints represented by h(·) ∈ Rnh and φ(·) ∈ Rnφ, respectively. For simplicity, we assume that path constraints are formulated separately for the state variables, hx(·), and the control variables, hu(·). Further, all constraints are assumed to be simple bound (or box) constraints. More complex constraints and combined state and control path constraints can be converted into this formulation by adding additional equations.12 Following Pontryagin’s Minimum Principle,32 dynamic optimization problem of eq 1 can be reformulated with the Hamiltonian function H(t) as nx

min H(t ) = g ( ·) + λT f( ·) + μT h( ·)

u(t ), t f

x ̇ (t ) = 7067

∂H , x(t0) = x 0 ∂λ

(2-a) (2-b)

dx.doi.org/10.1021/ie404148j | Ind. Eng. Chem. Res. 2014, 53, 7066−7078

Industrial & Engineering Chemistry Research T

λ ̇ (t ) = −

∂ϕ ∂φ ∂H T , λ (t f ) = + νT ∂x ∂xf ∂xf

0 = μT h( ·), 0 = v T φ(x(t f ))

Article

control profiles consisting of K arcs, then the original time interval t ∈ [t0,tf] is divided into K subinterval [tk−1,tk], k = 1, 2, ..., K, where t1, t2, ..., tk−1 are division points with t0 < t1, < t2 < ... 0, if hi( ·) = 0,

⎧ =0, if φ (x(tf )) < 0, j ⎪ νj(t )⎨ j = 1, 2, ..., nφ ⎪>0, if φj (x(tf )) = 0, ⎩

τ=

2t − (tk + tk − 1) , k ∈ [1, 2, ..., K ] tk − tk − 1

(4)

From eq 4, we obtain dτ 2 = ,1≤k≤K dt tk − tk − 1

(3-a)

(k)

(5)

(k)

Define x (τ) and u (τ) as the state and control in subinterval k, respectively. Based on eq 5, the problem of eq 1 can be rewritten as an MSP formulation

(3-b)

K

min J = ϕ(x(K )( +1)) + (k)

u (τ ), tk

k=1 +1

∫−1

(3-c)

Complementarity conditions by eqs 3-b and 3-c indicate that a particular Lagrange multiplier is positive if the corresponding constraint is active and zero otherwise. Optimal control profiles of eq 1 are typically discontinuous and consist of a distinct sequence of arcs; each arc is characterized by a different set of active constraints. These arcs are separated by so-called switching points, which encompass the time events where the input profile is discontinuous and may even jump. On the basis of the NCO, behaviors of a particular control profile ui(t) can be classified into four different types as follows:12 1. ui(t) = ui,min: the control is at its lower bound 2. ui(t) = ui,max: the control is at its upper bound 3. ui(t) = ui,path: the control tracks an active state path constraint 4. ui(t) = ui,sens: the control is sensitivity seeking (i.e., the control is not determined by any active constraints, but it is adjusted to minimize the objective function) Clearly, discontinuities of optimal control profiles typically occur for problems with active inequality path constraints. These problems are known to be difficult to solve using either an indirect or a direct method.34 Specifically, if we use an indirect method to solve an inequality path constrained dynamic optimization problem, the first-order optimality conditions from the calculus of variations are require modification in order to properly account for possible discontinuities in the optimal solution. On the other hand, when a direct method is carried out, the control is generally approximated by using a piecewise smooth function which may cause an inefficient and ill-conditioning discretized problem.

s. t .

∑

tk − tk − 1 2

(k)

g (x (τ ), u(k)(τ )) dτ

(6-a)

t − tk − 1 (k) dx(k)(τ ) = k f(x (τ ), u(k)(τ )), x(1)( −1) dτ 2 = x0

(6-b)

h(x(k)(τ ), u(k)(τ )) ≤ 0

(6-c)

φ(x(K )( +1)) ≤ 0

(6-d)

x(k)( +1) = x(k + 1)( −1), 1 ≤ k ≤ K − 1

(6-e)

Equation 6-e is known as linkage constraints, ensuring the continuity of state variables at the interface of subintervals. The MSP formulation of eq 6 consists of a collection of single stage problems, each of which corresponds to a segment of the complete profiles. The complete profiles are then obtained by minimizing the total cost function, subjecting to constraints within each stage and linkage constraints connecting adjacent stages. 3.2. Legendre−Gauss−Radau Discretization. Once the MSP shown as eq 6 is obtained, the next step is to transfer it into a discrete nonlinear programming problem (NLP), which can be solved numerically by well-developed algorithms. In this paper, the transformation of MSP into NLP is accomplished by the multiple-interval LGR method,20,29 where the state variable is approximated with Nk LGR points in subinterval k ∈ [1, 2, ..., K] as x(k)(τ ) ≈ X(k)(τ ) Nk + 1

3. LEGENDRE−GAUSS−RADAU PSEUDOSPECTRAL METHOD 3.1. Multistage Problem Formulation. As previously mentioned, computational problems, which are associated with using direct methods to solve problems with discontinuous control profiles, usually occur because of the inaccurate approximation of a nonsmooth function by a finite number of smooth functions. To deal with these issues, one simple yet efficient way is to reformulate the original problem of eq 1 as a multistage problem, as described by Schlegel and Marquardt in the context of the CVP method.12 Specifically, assuming

=

∑

X(jk)L(j k)(τ ), L(j k)(τ )

j=1 Nk + 1

=

∏ l=1 l≠j

τ − τl(k) τj(k) − τl(k) (7)

where L(k) j (τ), j polynomials; τ(k) 1 ,

= 1, 2, ..., Nk + 1 is a basis of Lagrange (k) τ(k) 2 , ..., τNk are the LGR collocation points in subinterval k defined on τ(k) ∈ [tk−1, tk); and τ(k) Nk+1 = +1 is a noncollocated point defining the end of the subinterval k. It

7068



Article

should be noticed that any interpolating function, whose value matches the approximation for the control at the Nk LGR collocation points in each subinterval, can be accepted. For simplicity, the control can also be approximated by using the Lagrange interpolating polynomial as U(k)(τ ) =

Nk

(k)

∑ U(i k)Lî

(k) (τ ), Lî (τ ) =

Nk

∏ l=1

i=1

Optimal solution of the resulting NLP fulfills the Karush− Kuhn−Tucker (KKT) conditions of optimality based on the Lagrangian function which is defined as ⎧ ⎫ ⎪ ϕ(X (NK )+ 1) − ⟨ṽ , φ(X (NK )+ 1)⟩ ⎪ K K ⎪ ⎪ K Nk ⎪ + ∑ ∑ ⎛⎜ tk − tk − 1 ⎞⎟w(k)g(k) − ⟨μ(̃ k) , h(k)⟩⎪ j ⎪ j ⎪ ⎝ ⎠ j j 2 k=1 j=1 ⎪ ⎪ ⎪ K − 1 Nk ⎪ ⎪ − ∑ ∑ Λ̃(k), D(k) X(k) + D(k) X(k + 1) ⎪ j j ,1: Nk 1: Nk j , Nk + 1 1 ⎪ ⎪ ⎬ Lf = ⎨ k = 1 j = 1 tk − tk − 1 (k) ⎪ ⎪ fj ⎪ − ⎪ 2 ⎪ ⎪ Nk ⎪ ⎪ (K ) (K ) (K ) (K ) ⎪ −∑ Λ̃ j , D(jK,1:) NK X1: NK + Dj , NK + 1 X NK + 1 ⎪ ⎪ ⎪ j=1 ⎪ ⎪ tK − tK − 1 (K ) ⎪ ⎪ − fj ⎩ ⎭ 2

τ − τl(k) τi(k) − τl(k)

l≠i

(8)

U(k) 1

where (i = 1, 2, ..., Nk) is the approximation of the control variable at the Nk LGR point in subinterval k ∈ [1, 2, ..., K]. The cost function of eq 6-a is then approximated using the LGR quadrature as K

J = ϕ(X (NKK)+ 1) +

Nk

⎛ tk − tk − 1 ⎞ (k) (k) (k) ⎜ ⎟w g (X i , U i ) ⎠ i 2

∑∑⎝ k=1 i=1

(9)

X(K) Nk+1

w(k) i

where is the approximation of x(tf), (i = 1, 2, ..., Nk) is the LGR quadrature weights37 in subinterval k ∈ [1, 2, ..., K]. Differentiating X(k)(τ) in eq 7 with respect to τ, we have dX(k)(τ ) ≡ Ẋ (k)(τ ) = dτ

Nk + 1

∑

(18) φ

(10)

j=1

Based on eq 10, the differential equations of eq 6-b with Nk LGR points can be now written as Nk + 1

∑

Dij(k) X(jk) −

j=1

tk − tk − 1 (k) (k) f(X i , U i ) = 0 2

(11)

Dij(k) = Lj̇(k)(τi(k)), i = 1, 2, ..., Nk ; j = 1, 2, ..., Nk + 1 (12) 37 where D(k) in ij is the Nk × (Nk +1) differentiation matrix subinterval k ∈ [1, 2, ..., K]. Path constraints of eq 6-c in subinterval k are enforced at Nk LGR points as

h(X(i k), U(i k)) ≤ 0, i = 1, 2, ..., Nk

(13)

(14)

The continuities across mesh points are maintained by the following condition:

X (Nkk)+ 1 = X1(k + 1)

(15)

(k+1) Since the same variables are used for both X(k) , Nk+1 and X1 constraints of eq 15 can be removed and constraints of eq 11 can be written as27 t − tk − 1 (k) (k) (k) (k + 1) D(jk,1:) Nk X1: − k fj Nk + Dj , Nk + 1 X1 2 = 0, 1 ≤ k ≤ K − 1 (16)

(K ) (K ) (K ) D(jK,1:) NK X1: NK + Dj , NK + 1 X NK + 1 −

tK − tK − 1 (K ) fj = 0 2

h

nz

4. A UNIFIED FRAMEWORK FOR SOLVING PROBLEMS WITH DISCONTINUOUS CONTROL PROFILES 4.1. Adaptive LGR Method. The construction of the MSP of eq 6 is based on the assumption that the number and sequence of arcs in control profiles are known in advance. Unfortunately, for most dynamic optimization problems, this assumption is not valid, since the solution cannot be known in advance. In addition, the number of LGR points used for control and state parametrizations should be adapted for the specific arc type to obtain a solution with fewer degrees of freedom but higher accuracy. Therefore, how to divide subintervals and how to adjust the degree of approximating polynomial within each subinterval are two crucial components for applying pseudospectral methods to solve nonsmooth problems. Accordingly, a mesh refinement algorithm referred to as the hp adaptive LGR method,26 which allows for changes in both the number of subintervals and the degree of the approximating polynomial within a subinterval, is adopted in this paper. This method relies on an estimated error to decide whether the degree of the approximating polynomial within an interval should be increased or the interval should be further divided into more subintervals. Specifically, in subinterval k, (k) (k) given a set of Mk = Nk + 1 LGR points (τ̂(k) 1 , τ̂2 , ..., τ̂Mk ) with (k) (k) τ̂(k) 1 = τ1 = tk−1 and τ̂Mk+1 = tk, values of the state approximation (k) (k) given in eq 7 are denoted as (X(τ̂(k) 1 ), X(τ̂2 ), ..., X(τ̂Mk )). If the (k) derivative of X̂ matches the dynamics at each of the LGR points τ̂(k) i , we have

The terminal constraints of eq 6-d are approximated as φ(X (NKK)+ 1) ≤ 0

̃ (k)

where ṽ ∈ R , μ̃ ∈ R , and Λ ∈ R are Lagrange multipliers associated with discretized terminal constraints of eq 14, discretized path constraints of eq 13, and discretized process constraints of eq 11 in subinterval k, respectively. The operator ⟨a,b⟩ is used to denote the standard inner product between vectors a ∈ Rn and b ∈ Rn. The KKT optimality conditions can be obtained by differentiating the Lagrangian function of eq 18 with respect to each variable and equating them to zero. It is important to note that values of discrete multipliers provide information about the status of a particular constraint at the optimal solution, which can be further used to extract the solution structure of optimal control profiles.12

X(jk)Lj̇(k)(τ )

(k)

(17)

Finally, the NLP that arises from LGR approximation is used to minimize the cost function of eq 9, subjecting to algebraic constraints of eq 13, eq 14, eq 16, and eq 17. As long as the NLP problem is formulated, it can be solved efficiently by exploiting its sparse structure.38 7069


Industrial & Engineering Chemistry Research (k) X̂ (τi(̂ k)) = X(k)(τk − 1) +

(X

(k)

(τl(̂ k)),

U

tk − tk − 1 2

Article

Mk

∑ Ijl(̂ k)f l=1

(k)

(τl(̂ k)))

(19)

̂ where I(k) jl , j,l = 1, 2, ..., Mk, is the Mk × Mk LGR integration (k) matrix corresponding to the LGR points defined at (τ̂(k) 1 , τ̂2 , ..., (k) (k) (k) (k) (k) τ̂Mk ). On the basis of X (τ̂l ) and X̂ (τ̂l ), the absolute and (k) (k) relative errors of the state at (τ̂(k) 1 , τ̂2 , ..., τ̂Mk+1) are then defined as

(k)

Ei(k)(τl(̂ k)) = |Yî (τl(̂ k)) − Y i(k)(τl(̂ k))|

ei(k)(τl(̂ k)) =

(20)

Remark 1: The preceding procedure only demonstrates the basic idea of a ph adaptive LGR method. A detailed derivation and implementation of this method can be found in Rao et al.26 4.2. Control Structure Detection with Optimality Validation. Even though significant advances in the hp adaptive pseudospectral method have been achieved, it still exhibits potential for further improvement, especially for solving problems with discontinuous control profiles. For example, a primary issue of this adaptive method is the termination with several unnecessary subintervals, which may result in an accumulation of LGR points around discontinuities and even oscillation phenomena because of the lack of mechanisms for deleting or merging redundant subintervals. In addition, currently, the decision to stop the iterative procedure is based on the local error analysis of state approximations, i.e., verifying if the relative error derived from eq 22 is below the specified tolerance at all LGR points. Alternatively, this can be accomplished by monitoring changes in the cost function: when improvements in the cost values are sufficiently small, the iteration is terminated.11 However, neither of them offers a way to quantitatively measure the distance between the solution to the discretized NLP and the true solution to original dynamic optimization problem.23,31 Hence, from these stopping criteria, we have no idea whether the optimal control policy is actually achieved. Obviously, if it was possible to remove these redundant subintervals reasonably, a solution with fewer degrees of freedom but higher accuracy might be achieved, which was particularly valuable for practical applications. To this end, a methodology combining the hp adaptive pseudospectral method with the structure detection algorithm is presented. The structure detection algorithm, as described in the context of CVP method,12 was built on the information about the status of control and path constraints. As for the LGR method, such information can be obtained from values of these Lagrange multipliers at LGR points. Each multiplier is linked to a specific constraint and will be enforced to nonzero if the corresponding constraint is active. Specifically, subintervals with controls at their upper or lower bounds, i.e., ui,max or ui,min, can be detected by nonzero multipliers of control path constraints of eq 13. Analogously, subintervals containing active constraints ui,path are detected by nonzero multipliers of state path constraints, and the remaining subintervals are denoted as sensitivity-seeking arcs ui,sens. A more complicated situation can be encountered, especially in the presence of state constraints and/or multiple inputs. Adjacent subintervals with the same input type are merged into one arc to remove unnecessary subintervals, and a new discretization scheme for the next refinement step is now

Ei(k)(τl(̂ k)) 1+

max

|Y i(k)(τj(̂ k))|

j ∈ [1,2, ···, Mk + 1]

(21)

with l = 1, 2, ..., Mk; i = 1, 2, ..., nx. The maximum relative error in subinterval k is then defined by comparing two approximations of the state as26 (k) emax =

max

i ∈ [1,2, ···, nx], l ∈ [1,2, ···, Mk + 1]

ei(k)(τl(̂ k))

(22)

Let ε be the desired relative error accuracy tolerance, the kth subinterval is considered to be within the accuracy tolerance if the maximum violation in the differential-algebraic equations on [tk−1,tk] is below ε. On the other hand, if e(k) max > ε, i.e., ε is not satisfied in subinterval k, according to the convergence theory of pseudospectral method, the number of LGR points required in this subinterval is26 ⎛ ⎛ e (k) ⎞⎞ Nk̃ = Nk + ceil⎜⎜log N ⎜⎜ max ⎟⎟⎟⎟ k ⎝ ε ⎠⎠ ⎝

(23)

where ceil is the operator that rounds to the next highest integer. Let Nmin and Nmax be the user-specified minimum and maximum bounds on the number of LGR points within subinterval k, respectively. If Ñ k ≤ Nmax, then Nk is increased to Ñ k to achieve the desired error tolerance. If, on the other hand, Ñ k > Nmax, then Ñ k exceeds the upper limit, and the subinterval k should be further divided into more subintervals. In this case, the number of new subintervals, denoted as nk, can be computed through eq 24. ⎛ ⎛ Ñ ⎞ ⎞ nk = max⎜⎜ceil⎜ k ⎟ , 2⎟⎟ ⎝ ⎝ Nmin ⎠ ⎠

(24)

Equation 24 ensures that the sum of LGR points in the newly created subintervals equals the number of LGR points that would be used in subinterval k if Ñ k was accepted. In addition, each newly created subinterval would contain the minimum allowable number of LGR points. The iterative reduction of the approximation error can be achieved through either increasing the degree of the polynomial in a subinterval or increasing the number of subintervals. An overview of this procedure can be found as follows: 7070



Article

Figure 1. Framework for solving a dynamic optimization problem with discontinuous control profiles.

generated and solved again, where the previous solution serves as an initial guess. The main algorithm behind the adaptive LGR method is to design a sequence of NLPs such that a solution to the final problem is sufficiently close to the exact solution of the dynamic optimization problem. To stop the refinement as soon as the resulting solution has the optimal number of subintervals, a stopping criterion by monitoring variations of the Hamiltonian function at LGR points is presented. Specifically, as a key property of the LGR method, KKT conditions of the NLP of eq 18 are exactly equivalent to the discretized firstorder necessary conditions of the continuous Bolza problem in eq 1.28 This feature provides an easy way to estimate costate variables by using KKT multipliers of the NLP as λj(k)

=

(k) Λ̃

w(j k)

,

μj(k)

Based on eq 25, the discretized Hamiltonian function in subinterval k at LGR point j is defined as H (j k) = g (X(jk), U(jk)) + ⟨λj(k) , f(X(jk), U(jk))⟩ + ⟨μj(k) , h(X(jk), U(jk))⟩, 1 ≤ j ≤ Nk

Once the NLP algorithm converges, the last two terms of eq 26 equal zero because of feasibility and complementarity conditions. From optimal control theory, the Hamiltonian should be constant.39,40 H (j k) = H(tf ) = constant

(27)

Equation 27 can be used to verify the optimality of a

μj(̃ k) 2 = , tk − tk − 1 w(j k)

1 ≤ k ≤ K , 1 ≤ j ≤ Nk

(26)

computed solution. If the Hamiltonian function is not constant over all subintervals, the corresponding solution is not optimal, (25)

and subintervals should be further refined. 7071



Article

Figure 2. Optimal control profiles obtained by the CVP method with 32 and 256 parameters.

σH =

1 K × Nk

1 = K × Nk

K

solution of increasingly refined NLP generates a problemadapted discretization, leading to a better approximation for control profiles with fewer degrees of freedom.

Nk

∑ ∑ (H (j k) − mH)2 , mH k=1 j=1 K

Nk

5. CASE STUDY The Williams−Otto semibatch reactor29 has been chosen for illustrating the effectiveness of the proposed method. Reactions A + B → C, C + B → P + E, and P + C → G take place in the reactor. Reactant A already exists in the reactor, and reactant B is fed continuously to accelerate reactions. Through these exothermic reactions, products P and E as well as side-product G are formed. The heat generated from exothermic reactions is removed by a cooling jacket, where the cooling water temperature is the manipulated variable. During the batch, path constraints on the inlet flow rate of reactant B, denoted as FB,in, the reactor temperature TR, the reactor volume V, and the scaled cooling water temperature TW must be observed. The economic objective is to maximize yields of main products at the end of the batch, with control variables FB,in(t) and TW(t). In summary, the dynamic optimization problem can be defined as follows:30,31

∑ ∑ H (j k) k=1 j=1

(28)

Here, defined as eq 28, the mean square error of the discretized Hamiltonian is utilized to estimate the distance of the computed solution to the true solution. This measure relies on a multilevel solution strategy,28 where the optimization problem is solved with increasingly fine discretization by the adaptive LGR method. Then, the error can be used to decide if the obtained solution is sufficiently close to the true solution, i.e., when to stop the refinement, by simply checking the following condition |σH| ≤ tol

(29)

where tol is a user-defined tolerance error. It is found that the Hamiltonian function is sensitive to changes in the solution structure and thus provides an intuitive way to find whether the optimal number of subintervals should be chosen greater or smaller. An integrated framework for solving dynamic optimization problems with discontinuous control profiles can now be established by linking adaptation, structure detection, and optimality validation along the lines presented in previous sections. The procedure of this framework is described with the aid of Figure 1. Given an optimization problem, the algorithm starts by transcribing it into a NLP formulation by applying LGR discretization for each subinterval. At the solution of the NLP, the error e(k) max for each subinterval is calculated according to eq 22. If the accuracy tolerance is met, i.e., e(k) max ≤ ε, values of the Hamiltonian function at LGR points for each subinterval are estimated based on eqs 25 and 26. Further, if the condition |σH| ≤ tol is met, the procedure terminates and offers the final solution of the optimization problem. Otherwise, current subintervals will be classified and merged by applying the structure detection procedure, resulting in a new NLP formulation for the next refinement step. The repetitive

min

FB,in(t ), TW(t )

Φ(tf )

s . t . xȦ = −

x Ḃ =

7072

1 xAFB,in − k1η1xAx B 1000 V

1 (1 − x B)FB,in − k1η1xAx B − k 2η2x BxC 1000 V

xĊ = −

1 xCFB,in + k 7η1xAx B − k 3η2x BxC − k6η3xCx P 1000 V

x Ṗ = −

1 x PFB,in + k 2η2x BxC − k4η3xCx P 1000 V

x Ė = −

1 x EFB,in + k 3η2x BxC 1000 V dx.doi.org/10.1021/ie404148j | Ind. Eng. Chem. Res. 2014, 53, 7066−7078


Article

Figure 3. (a) Optimal control profiles obtained by the adaptive LGR method with the accuracy tolerance ε = 10−3. (b) Optimal control profiles obtained by the adaptive LGR method with the accuracy tolerance ε = 10−6.

xĠ = − TṘ =

1 xGFB,in + k5η3xCx P 1000 V

improvement in objective function. The resulting control profiles of this problem by using a coarse grid (32 parameters for each control) and fine grid (256 parameters for each control), are shown in Figure 2. It can be observed from Figure 2 that the resulting control profiles contain a sequence of arcs and show discontinuous jumps at some points. Compared with the coarser solution, profiles with the fine grid are much smoother owing to the finer discretization. Objective function values for two such parametrization methods are −4767.56 and −4768.27 with computational times of 65.6 and 1887.7 s, respectively. With 256 control parameters, the resulting NLP becomes intractable computationally, which is quite common for numerical solution methods. Furthermore, it is not able to find true locations of discontinuities in control profiles, although it is possible to capture the true solution with this choice of parametrization. The dynamic optimization problem of eq 30 is then solved by the adaptive LGR method with the accuracy tolerances ε = 10−3 and ε = 10−6, respectively. The minimum and maximum number of allowable LGR points in each subinterval are set as Nmin = 4 and Nmax = 24, respectively. Figure 3a and b depict the obtained optimal control profiles under the two different accuracy tolerances. For the relatively large tolerance ε = 10−3, the adaptive LGR method terminated with a global polynomial function in interval [0, 1000] consisting of 21 LGR points. It

1 (Tin − TR )FB,in + k 8η1xAx B + k 9η2x BxC 1000 V + k10η3xCx P − l1TR + l 2TW

V̇ =

1 FB,in 1000

Φ̇ = −5554.1(k 2η2x BxC − k4η3xCx P)V − k11η2x BxCV

V (t f ) ≤ 5 60 ≤ TR (t ) ≤ 90, ∀ t ∈ [0, tf ]

0 ≤ FB,in(t ) ≤ 5.784, ∀ t ∈ [0, tf ] 0.02 ≤ TW(t ) ≤ 0.1, ∀ t ∈ [0, tf ]

with nx = 9 states, nu = 2 controls, and the batch time tf = 1000 s. The above dynamic optimization problem has been solved by using the CVP method,11 with an iterative refinement of the parametrization and a termination criterion based on the 7073



Article

Figure 4. (a) Optimal control profiles based on the proposed methodology with the accuracy tolerance ε = 10−3. (b) Optimal control profiles based on the proposed methodology with the accuracy tolerance ε = 10−6.

can be seen from Figure 3a that accuracy tolerance ε = 10−3 may be too large to achieve a desirable solution. On the other hand, when the adaptive LGR method works with the accuracy tolerance ε = 10−6, the interval [0, 1000] is divided into 15 subintervals and contains 218 LGR points. Clearly, lower tolerance ε requires a larger number of LGR points and subintervals to generate a solution that is in better agreement with the optimal solution. However, as shown in Figure 3b, an accumulation of LGR points around the control discontinuity would occur, and most of these points are not actually required for representing the true solution. This phenomenon can be explained in that the exact location of the discontinuity cannot be captured by the adaptive LGR method, which can only improve the approximation accuracy by increasing the number of LGR points and subintervals. When ε = 10−6 is used, the values of the objective function and computational time are −4768.31 and 14.6 s, respectively. Compared with CVP methods, both of these two indicators have been improved, especially for the computational burden. Finally, the proposed solution framework is applied to solve the optimization problem with accuracy tolerances ε = 10−3 and ε = 10−6, respectively. For each scenario, the minimum and maximum number of allowable LGR points in a subinterval are also set to Nmin = 4 and Nmax = 24, respectively. On the other hand, the tolerance for variations of the Hamiltonian function is set as tol = 10−3 to result in an appropriate termination of the method and to keep NLPs within a reasonable size. Figure 4a

and b demonstrate optimal control profiles with two different accuracy tolerances. Six switching times between them are [52.33, 148.25, 359.9, 518.7, 538.68, 858.81] s. Under both accuracy tolerances, all relevant arcs of the control switching structure can be detected by the proposed methodology. The proposed method obtains an objective function of −4768.313 with the computational time of 9.62 s when the accuracy tolerance is set as ε = 10−6. The corresponding optimal profiles of constrained state variables Tr and V are depicted in Figure 5a and b. In order to demonstrate the optimality of the resulting control profiles, values of the Hamiltonian function at each LGR point were calculated and illustrated in Figure 6a and b. Apparently, the Hamiltonian function remains reasonably constant over the time interval [0, 1000] under different accuracy tolerances, indicating all relevant arcs of the control switching structure have already been presented in the obtained solution. Tables 1 and 2 show the quantitative comparisons of the proposed methodology and CVP based method30,31 and adaptive LGR method, respectively. Table 1 clearly shows that the proposed method in this paper is orders of magnitude faster than the CVP based method while achieving the optimal solution with high quality. Recalling Figure 3 and comparing it with Figure 4, it can be seen that the accumulation of unnecessary LGR points around the switching points can be tackled by the proposed methodology as the structure detection 7074



Article

Figure 5. (a) Optimal profiles for the constrained reactor temperature and reactor volume with the accuracy tolerance ε = 10−3. (b) Optimal profiles for the constrained reactor temperature and reactor volume with the accuracy tolerance ε = 10−6.

Figure 6. (a) The values of the Hamiltonian function at LGR points with the accuracy tolerance ε = 10−3. (b) The values of the Hamiltonian function at LGR points with the accuracy tolerance ε = 10−6.

Table 1. Comparison of the Proposed Method with CVP Based Method switch points (s) proposed method CVP method

τ1

τ2

τ3

τ4

τ5

τ6

objective function

CPU times (s)

52.33 52.34

148.25 148.26

359.9 359.9

518.7 518.7

538.68 538.68

858.81 858.79

−4768.31 −4768.27

9.62 1887.7

7075



Article

method, with the help of the control structure detection procedure, adjacent subintervals with the same input type would be merged into one arc to remove unnecessary subintervals. Tables 4 and 5 display the performance of the proposed method using various values of Nmin and Nmax while keeping the accuracy tolerance constant as ε = 10−5. As can be seen from Table 4, for Nmin = 2, 4, and 6, the proposed method would terminate with the same solution with 11, 8, and 6 refinement iterations, respectively. Interestingly, a smaller value of Nmin would require more refinement iterations to classify and merge those unnecessary subintervals based on the structure detection procedure. Compared with the results associated with the smaller value of Nmin, the method with Nmin = 8 and 10 would terminate with more LGR points, since Nmin is greater than the number of LGR points to achieve the required accuracy tolerances in some subintervals. Therefore, Nmin represents a tuning parameter that can trade off the conflict between the number of refinement iterations and the number of LGR points. Table 5 outlines the result accuracy and computational speed of the proposed method under different Nmax’s. For Nmax = 24, 28, and 32, the proposed method would result in the same solution after the same number of refinement iterations. This means Nmax = 24 is large enough to meet the required accuracy tolerance in each subinterval. Compared with Nmax = 20, the method with Nmax = 16 would take two additional iterations to reach the final solution. From the above analysis, one may conclude that a larger Nmax would be always preferred. However, a large Nmax may reduce the number of subintervals, and thus certain arcs of the control switching structure might be missed.

Table 2. Comparison of the Proposed Method with Adaptive LGR Method

proposed method adaptive LGR method

number of LGR points

number of subintervals

objective function

CPU times (s)

69

6

−4768.31

9.62

218

15

−4768.31

14.6

procedure can classify and merge the subintervals. Table 2 gives detailed quantitative information about the number of subintervals and LGR points for both methods. In total, compared with the CVP and adaptive LGR methods, less computational time will be needed to obtain the optimal solution with comparable accuracy based on the proposed method.

6. DISCUSSION To clearly demonstrate the capability of the proposed methodology, it is important to investigate influences of several important tuning parameters (i.e., the accuracy tolerance ε, the minimum and maximum number of allowable LGR points in each subinterval Nmin and Nmax) on its performance. Based on the aforementioned studied case, this section is dedicated to revealing relationships between the above three key parameters and the performance. The base parameters are Nmin = 4, Nmax = 24, and ε = 10−5. If one of these parameters is analyzed, the remaining parameters are kept as the base values. Table 3 summarizes the main performance of the method under different accuracy tolerance, ε. As described in eq 23, ε controls the growth of the number of LGR points in a subinterval; the decrease of ε would therefore lead to an increase in the computation burden, the number of LGR points, along with the number of refinement iterations and the CPU time. On the other hand, the accuracy of the solution, in terms of e(k) max, σH, and the value of the objective function can be improved accordingly. Hence, an appropriate value of ε should make a trade-off between computational effort and accuracy. From Table 3, ε = 10−5 and ε = 10−6 are excellent candidates because they would result in high quality solutions while maintaining a relatively small computational time. However, it should be noted that a small ε might lead to a computationally intractable problem for higher dimensional problems (i.e., problems with more states and controls) because of the increased memory and computational requirements associated with the increased LGR points and subintervals. As described in eq 24, Nmin controls the number of newly created subintervals, while Nmax is responsible for deciding whether a subinterval should be further divided into more subintervals. It can be expected that the number of subintervals would increase as either Nmin or Nmax decreases for the adaptive LGR method, in other words, more subintervals than actually required would be generated. However, in the proposed

7. CONCLUSIONS In this paper, a numerical solution framework combining the adaptive LGR method with the structure detection procedure was presented, to cope with computational difficulties that arise from dynamic optimization problems with discontinuous control profiles. Initially, the adaptive LGR method was applied to obtain a first solution with the possible sequence and type of subintervals in control profiles automatically. Subsequently, the structure detection procedure was performed to classify such subintervals and to remove unnecessary LGR points for the approximation of the state and control variables. In this way, the drawback of accumulation of LGR points around the switching points can be resolved. Furthermore, since the length of each subinterval is defined as a decision variable in optimization, the solution not only yields optimal control profiles but also provides precise values for discontinuous points presented in control structure. The overall procedure would terminate as long as the variation of Hamiltonian function reached the desired error tolerance, resulting in a solution that could capture switching points accurately. Studies on the dynamic optimization of the Williams-Otto reactor

Table 3. Summary of Accuracy and Speed of the Proposed Method under Different Accuracy Tolerances ε ε ε ε ε ε ε

= = = = =

−3

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



number of refinement iteration

35 39 52 69 95

6 6 6 6 9

6 6 8 11 17 7076

e(k) max 4.92 5.87 8.21 3.64 6.89

× × × × ×

10−4 10−5 10−6 10−7 10−8

σH

objective function

CPU time (s)

0.0014 0.0011 2.26e-004 2.95e-005 2.93e-005

−4768.30 −4768.30 −4768.31 −4768.31 −4768.31

4.43 4.79 7.38 9.62 20.78



Article

Table 4. Summary of Accuracy and Speed of the Proposed Method Using Different Nmin Nmin




2 4 6 8 10

52 52 52 63 77

6 6 6 6 6

11 8 6 5 5

σH

objective function

CPU time (s)

2.26e-004 2.26e-004 2.26e-004 2.42e-005 3.95e-005

−4768.31 −4768.31 −4768.31 −4768.31 −4768.31

15.88 7.38 6.66 5.40 6.21

σH

objective function

CPU time (s)

6.26e-004 6.26e-004 2.26e-004 2.26e-004 2.26e-004

−4768.31 −4768.31 −4768.31 −4768.31 −4768.31

16.18 15.53 7.38 7.38 7.38

e(k) max 8.21 8.21 8.21 4.75 7.75

× × × × ×

10−6 10−6 10−6 10−6 10−6

Table 5. Summary of Accuracy and Speed of the Proposed Method under Different Nmax Nmax




16 20 24 28 32

55 55 52 52 52

7 7 6 6 6

13 11 8 8 8

7.67 7.67 8.21 8.21 8.21

AUTHOR INFORMATION

Corresponding Authors

*Tel.: +86 532 86981169. Fax: +86 532 86981718. E-mail: [email protected]. *Tel.: +1 412 268 2261. Fax: +1 412 268 7139. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Basic Research Program of China (973, Grant No. 2012CB215006), the Fundamental Research Funds for the Central Universities (Grant No. 27R1305010A), and China Postdoctoral Science Foundation (2013M541964); their supports are hereby acknowledged.

■

× × × × ×

10−6 10−6 10−6 10−6 10−6

(5) Vassiliadis, V.; Sargent, R.; Pantelides, C. Solution of a class of multistage dynamic optimization problems: 1. Problems without path constraints. Ind. Eng. Chem. Res. 1994, 33, 2111−2122. (6) Vassiliadis, V.; Sargent, R.; Pantelides, C. Solution of a class of multistage dynamic optimization problems: 2. Problems with path constraints. Ind. Eng. Chem. Res. 1994, 33, 2123−2133. (7) Bock, H. G.; Plitt, K. J. A Multiple Shooting Algorithm for Direct Solution of Optimal Control Problem. In Proceedings 9th IFAC World Congress; Pergamon Press: Elmsford, NY, 1984. (8) Biegler, L. T. Solution of dynamic optimization problems by successive quadratic programming and orthogonal collocation. Comput. Chem. Eng. 1984, 8, 243−247. (9) Barton, P. I.; Allgor, R. J.; Feehery, W. F.; Galan, S. Dynamic Optimization in a Discontinuous World. Ind. Eng. Chem. Res. 1998, 37, 966−981. (10) Betts, J.; Huffmann, W. Mesh refinement in direct transcription methods for optimal control. Optim. Control Appl. Methods 1998, 19, 1−21. (11) Schlegel, M.; Stockmann, K.; Binder, T.; Marquardt, W. Dynamic optimization using adaptive control vector parametrization. Comput. Chem. Eng. 2005, 29, 1731−1751. (12) Schlegel, M.; Marquardt, W. Detection and exploitation of the control switching structure in the solution of dynamic optimization problems. J. Process Control 2006, 16, 275−290. (13) Schlegel, M.; Marquardt, W. Adaptive switching structure detection for the solution of dynamic optimization problems. Ind. Eng. Chem. Res. 2006, 45, 8083−8094. (14) Huntington, G. T.; Rao, A. V. Comparison of local and global collocation methods for optimal control. J. Guid. Control Dyn. 2008, 31, 432−436. (15) Ross, I. M.; Karpenko, M. A review of pseudospectral optimal control: From theory to flight. Ann. Rev. Control 2012, 36, 182−197. (16) Benson, D. A. A Gauss pseudospectral transcription for optimal control. Ph.D. Thesis, Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA, 2004. (17) Benson, D. A.; Huntington, G. T.; Thorvaldsen, T. P.; Rao, A. V. Direct trajectory optimization and costate estimation via an orthogonal collocation method. J. Guid. Control Dyn. 2006, 29, 1435− 1440. (18) Elnagar, G.; Kazemi, M.; Razzaghi, M. The pseudospectral Legendre method for discretizing optimal control problems. IEEE Trans. Autom. Control 1995, 40, 1793−1796. (19) Kameswaran, S.; Biegler, L. T. Convergence rates for direct transcription of optimal control problems using collocation at Radau points. Comput. Optim. Appl. 2008, 41, 81−126. (20) Garg, D.; Patterson, M. A.; Hager, W. W.; Rao, A. V.; Benson, D. A.; Huntington, G. T. A unified framework for the numerical solution of optimal control problems using pseudospectral methods. Automatica 2010, 46, 1843−1851.

highlighted main features of the presented method. Generally speaking, the illustrative example demonstrated that the present method allowed one to detect the switching structure and to obtain an accurate control profile resolution with an acceptable computational effort simultaneously. On the other hand, based on the case study of the Williams-Otto reactor, relationships between three important tuning parameters (i.e., the accuracy tolerance ε, the minimum and maximum number of allowable LGR points in each subinterval Nmin and Nmax) and the performance were also elucidated. Of course, these relationships can guide further applications of the proposed methodology to large-scale problems.

■

e(k) max

REFERENCES

(1) Srinivasan, B.; Palanki, S.; Bonvin, D. Dynamic optimization of batch processes I. Characterization of the nominal solution. Comput. Chem. Eng. 2003, 27, 1−26. (2) Binder, T.; Blank, L.; Bock, H. G.; Bulirsch, R.; Dahmen, W.; Diehl, M.; Kronseder, T.; Marquardt, W.; Schlöder, J. P.; von Stryk, O. Introduction to model based optimization of chemical processes on moving horizons. In Online Optimization of Large Scale Systems; Grötschel, M., Krumke, S. O., Rambau, J., Eds.; Springer-Verlag: Berlin, 2001. (3) Engell, S. Feedback control for optimal process operation. J. Process Control. 2007, 17, 203−219. (4) Betts, J. Practical Methods for Optimal Control Using Nonlinear Programming; SIAM Series on Advances in Design and Control; SIAM: Philadelphia, PA, 2001. 7077



Article

(21) Fahroo, F.; Ross, I. M. Advances in pseudospectral methods for optimal control. In Proceedings AIAA guidance, navigation, and control conference; Pergamon Press: Honolulu, HI, 2008. (22) Ross, I. M.; Fahroo, F. Pseudospectral knotting methods for solving optimal control problems. J. Guid. Control Dyn. 2004, 27, 397− 405. (23) Gong, Q.; Fahroo, F.; Ross, I. M. Spectral algorithm for pseudospectral methods in optimal control. J. Guid. Control Dyn 2008, 31, 460−471. (24) Darby, C. L.; Hager, W. W.; Rao, A. V. An hp-adaptive pseudospectral method for solving optimal control problems. Optim. Control Appl. Methods 2011, 32, 476−502. (25) Rao, A. V.; Patterson, M. A.; Hager, W. W.; A ph-collocation scheme for optimal control. Automatica, 2013. Submitted for Publication, to be found under the URL: http://vdol.mae.ufl.edu/ SubmittedJournalPublications/PattersonHagerRaoPH.pdf. (26) Darby, C. L.; Hager, W. W.; Rao, A. V. Direct trajectory optimization using a variable low-order adaptive pseudospectral method. J. Spacecr. Rockets 2011, 48, 433−445. (27) Benson, D. A.; Huntington, G. T.; Thorvaldsen, T. P.; Rao, A. V. Direct trajectory optimization and costate estimation via an orthogonal collocation method. J. Guid. Control Dyn. 2006, 29, 1435− 1440. (28) Darby, C. L.; Garg, D.; Rao, A. V. Costate estimation using multiple-interval pseudospectral methods. J. Spacecr. Rockets 2011, 48, 856−866. (29) Forbes, J. F. Model structure and adjustable parameter selection for operations optimizations. Ph. D. Thesis, McMaster University, Hamilton, Canada, 1994. (30) Hannemann, R.; Marquardt, W. Continuous and discrete composite adjoints for the Hessian of the Lagrangian in shooting algorithms for dynamic optimization. Siam J. Sci. Comput. 2010, 31, 4675−4695. (31) Hannemann, R.; Marquardt, W. How to verify optimal controls computed by direct shooting methods? -A tutorial. J. Process Control. 2012, 22, 494−507. (32) Pontryagin, L. S.; Boltyanskiy, V. G.; Gamkrelidze, R. V.; Mishchenko, Y. F. The Mathematical Theory of Optimal Processes; Wiley-Interscience: New York, 1962. (33) Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Hemisphere Publ.: New York, 1975. (34) Kameswaran, S.; Biegler, L. T. Simultaneous dynamic optimization strategies: Recent advances and challenges. Comput. Chem. Eng. 2006, 30 (10−12), 1560−1575. (35) Biegler, L. T. An overview of simultaneous strategies for dynamic optimization. Chem. Eng. Process. 2007, 46, 1043−1053. (36) Abramowitz, M.; Stegun, I. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables; Dover Publications: New York, 1965. (37) Garg, D.; Patterson, M. A.; Darby, C. L.; Francolin, C.; Huntington, G. T.; Hager, W. W.; Rao, A. V. Direct Trajectory Optimization and Costate Estimation of Finite-Horizon and InfiniteHorizon Optimal Control Problems Using a Radau Pseudospectral Method. Comput. Optim. Appl. 2011, 49, 335−358. (38) Patterson, M. A.; Rao, A. V. Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems. J. Spacecr. Rockets 2012, 39, 364−377. (39) Tanartkit, P.; Biegler, L. T. A nested, simultaneous approach for dynamic optimization problems. II. The outer problem. Comput. Chem. Eng. 1997, 21, 1365−1388. (40) Biegler, L. T.; Cervantes, A. M.; Wächter, A. Advances in simultaneous strategies for dynamic process optimization. Chem. Eng. Sci. 2002, 57, 575−593.

7078


The Combination of Adaptive Pseudospectral Method and Structure

Recommend Documents