Ind. Eng. Chem. Res. 2001, 40, 2089-2102
2089
Trajectory Optimization for Flat Dynamic Systems M. Guay,† S. Kansal,‡,§ and J. F. Forbes*,‡ Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6, and Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6
A number of important chemical engineering processes are operated in a transient manner (e.g., batch processes) and cannot be considered to reach a steady state. Optimizing the operations of such processes requires the solution of a dynamic optimization problem, producing time-based trajectories for process variables. A key characteristic of dynamic optimization problems is that the process model contains differential equations. Numerical solution techniques, which are currently in widespread use, are usually based on discretization schemes and can be computationally expensive. This paper proposes an alternative method for solving dynamic optimization problems in which the nonlinear process model is flat. The approach exploits, as appropriate, either the differential flatness or the orbital flatness of the process model to explicitly eliminate the differential equations from the optimization problem. The resulting optimization problem is solely algebraically constrained and can be solved using readily available optimization codes. The proposed approach is demonstrated on a range of benchmark problems taken from the literature. 1. Introduction The optimization of plant operations has been, and continues to be, of considerable interest to both academics and industrial practitioners alike. The majority of the developments have focused on determining the optimal steady-state operation of continuous processes. Unfortunately, there is a class of operations problems where such steady-state optimization approaches cannot be applied. Batch processes represent an important class of industrial operations that cannot be considered to ever reach a steady state. Batch processing is the basis for production in a variety of industries including specialty chemicals, pharmaceuticals, steel, and so forth. Furthermore, even in continuous operations, there are situations where the process is intentionally operated in a transient manner (e.g., grade transitions). Thus, there are a number of chemical engineering processes for which the operations optimization problem must be considered within a dynamic framework. In dynamic optimization problems, an objective function is optimized with respect to dynamic model equations (rather than the steady-state equations that are most often used for operations optimization problems) and other constraints. The solution of such dynamic optimization problems is usually some trajectory, in time, for the decision variables of the optimization problem. Such decision variables could include setpoints for the manipulated process variables (i.e., input variables) or measured process variables (i.e., output variables). It is then the task of the process controllers to track or implement these optimal trajectories. A general * Author to whom correspondence should be addressed. Phone: (780) 492-0873. Fax: (780) 492-2881. E-mail:
[email protected]. † Queen's University. ‡ University of Alberta. § Current address: Imperial Oil Ltd., Sarnia, ON, Canada.
form for dynamic process operations optimization problem is
min Φ(x(t),u(t)), t ∈ [t0, tf] x(t),u(t) subject to
f(x3 (t),x(t),u(t)) ) 0, x(t0) ) x0 c(x(t),u(t)) ) 0
(1)
g(x(t),u(t)) e 0 xL e x(t) e xU uL e u(t) e uU where Φ is the objective function, g and c are vectors containing algebraic functions that can include both initial-time (t0) and final-time (tf) constraints, x represents the process states, and u represents the process inputs. The subscripts L and U indicate lower and upper bounds, respectively. Further, the inequalities in problem 1 indicate a component-wise comparison. (This is a simplified version of the formulation given by Cuthrell and Biegler,6 but it contains all of the characteristics necessary for the developments and discussions of this paper.) Dynamic optimization problems, as given in problem 1, can be very difficult to solve analytically, and numerical solution techniques are typically used. Conventional nonlinear programming (NLP) techniques cannot be directly applied to problem 1 because of the presence of the differential equality constraints. Thus, special numerical techniques have to be employed in order to solve these problems. The most common solution techniques can be divided into two classes: (1) methods that rely on some form of approximation of the differential equations in order to pose the problem as a NLP and (2) methods that embed a differential equation solver into the optimization technique.
10.1021/ie0006312 CCC: $20.00 © 2001 American Chemical Society Published on Web 04/05/2001
2090
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
The transformation of problem 1 into a NLP form via discretization has been studied by a number of researchers (e.g., Biegler and co-workers,6,2,21 and Villadsen and Michelsen45). In these approaches to solving the dynamic optimization problem, the state and manipulated variables are parametrized in terms of some set of basis functions (e.g., simple polynomials45 or Lagrange interpolation polynomials21 in t), and the optimization horizon (i.e., t0, tf) is discretized in some fashion (e.g., via orthogonal collocation6). As a result of discretizing the optimization time horizon, the number of equality constraints is inflated by the fineness of this discretization. Thus, increasing the accuracy of the solution via a finer discretization of the optimization horizon can drastically increase the computational requirements of this approach. Methods that embed differential equation solvers within the optimization technique include control vector iteration (CVI),17 control vector parametrization (CVP),37 and iterative dynamic programming (IDP).7,23 In CVI, the input vector u(t) is assumed to be piecewiseconstant, and at each iteration of the optimization algorithm, the state equations are integrated forward in time; then, the results of the forward integration are used to integrate a set of adjoint equations backward in time to provide a correction to the current iterate for the input vector uk(t). CVP eliminates the need for the backward integration of the adjoint equations by expressing the input variable u(t) in terms of some predefined basis functions. IDP is a refinement of the general dynamic programming approach wherein at every iteration of the optimization algorithm (1) a set of discrete input variable profiles u(t) is used to integrate the state equations forward in time, (2) the grid is refined during a backward pass, (3) the grid is contracted around the best input profile identified during the backward pass, and (4) the procedure is repeated to convergence. In each of these approaches, at least one integration of the state equations is required at every iteration of the optimization algorithm, which can be computationally expensive. Furthermore, the accuracy of the solution depends on the discretization chosen for the optimization time horizon. Generally, numerical approaches can produce solutions that are arbitrarily close to the optimum of problem 1 via increasing the fineness of the discretization of the time grid used in the problem, increasing the order of the approximating polynomials, choosing appropriate basis functions, etc. However, increasing the solution accuracy usually drastically increases the computational requirements for solving the dynamic optimization problem. This tradeoff between solution accuracy and computational requirements is due to the differential equations contained in problem 1. If they could be transformed to algebraic equations, then conventional algebraic optimization techniques [i.e., linear programming (LP), quadratic programming (QP), and nonlinear programming (NLP)] could be used to efficiently and accurately solve the problem. Such a transformation is possible for flat nonlinear systems. In this paper, we consider an approach similar to that used by a number of researchers, including Agrawal et al.,1 Faiz et al.,8 Mahadevan et al.,24 and Oldenberg and Marquardt.30 The approach proposed in this paper extends beyond the class of differentially flat systems considered by other researchers to include the more general class of orbitally flat nonlinear systems and
involves a parametrization of the problem in the flat space that is computationally attractive. We show how problem 1 can be reduced to an algebraic optimization problem (i.e., LP, QP, or NLP) when no path constraints exist in the problem and to a semi-infinite optimization (SIO) problem in the presence of path and input constraints. A normalized form for the solution of a trajectory-generation problem for flat systems is developed. This form is general and applies to all differentially and orbitally flat systems. A formal approach for the solution of these problems is adopted and tested using examples reported in the dynamic optimization literature. 2. Mathematical Background Differential flatness, a notion introduced by Martin25 and studied further by Fliess et al.,9 Rathinam and Murray,33 and van Nieuwstadt,42 refers to the existence of so-called flat or linearizing outputs that summarize the dynamics of a nonlinear system. It is closely related to the general ability to linearize a nonlinear system by an appropriate choice of endogenous dynamic state feedback. Such feedback structures describe a special class of static or dynamic precompensators formed only of endogenous system variables, such as state variables and inputs and a finite number of their time derivatives. This special property has been used in the solution of a number of practical control problems involving chemical systems (e.g., Rothfuss35), mechanical systems (e.g., Rouchon36 and Murray28), and aircraft systems (e.g., Martin et al.26). A more general concept is the orbital flatness of nonlinear systems introduced by Fliess.9 This property refers to the linearizability of systems subject to an endogenous state feedback transformation and a state-dependent time-scaling transformation. It was shown by Fliess et al.11 and Guay15 that this concept can be viewed as a generalization of the differential flatness property to weakly locally accessible nonlinear systems. The orbital flatness of a simple chemical reactor due to Kravaris and Chung20 was demonstrated in Guay15. Differential flatness has been applied in the solution of tracking problems. Fliess and co-workers have provided a number of applications of differential flatness in the control of trailer systems (e.g., ref 9). The approach was also considered by Rothfuss35 in considering the assignment of trajectories for differentially flat chemical reactors. In this context, differential flatness becomes a very desirable property, as knowledge of the linearizing (or flat) outputs provides an algebraic parametrization of the dynamics of the system. Any trajectory in the flat output space can be mapped directly to corresponding state and input trajectories, without the explicit need to integrate the differential equations of the system. A number of researchers have considered the problem of trajectory assignment for flat systems. van Nieuwstadt and Murray41 first considered the problem of assigning trajectories in real time using flatness for various mechanical systems. They focused on a twodegree-of-freedom approach to controller design in which the flatness of the system was exploited for both assignment of the trajectories and tracking of the controller design. Optimal trajectory generation for realtime applications based on flatness was considered by Agrawal et al.1 In their work, linear controllable timeinvariant systems, the simplest form of flat systems,
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2091
were treated. They showed how trajectories in the flat output space can be assigned using polynomials parametrized by a relatively small number of constants. This approach greatly facilitates the computation of optimal trajectories in real time, as the dynamic optimization problem is transformed into an algebraic optimization problem. Faiz et al.8 extended this work to inequalityconstrained trajectory optimization problems for differentially flat nonlinear systems. The resulting SIO problem was solved via an off-line polytopic approximation/discretization scheme for the path constraints that solved the “locally reduced” problem using standard NLP techniques. Mahadevan et al.24 also restrict themselves to differentially flat systems and use a combination of discretization in time and polynomial parametrization of the flat outputs to transform the SIO to a NLP. Oldenberg and Marquardt30 employ a similar method in which the highest derivative for each of the flat outputs is parametrized by piecewise-constant functions with a discretized time horizon. As with the other approaches, the SIO is transformed into a NLP, which can be solved with standard optimization codes. In this work, the flat outputs are directly approximated with arbitrary polynomial functions in time, and the optimization problem is solved using either NLP or SIO techniques, depending on the presence of path constraints. No explicit discretization of the problem is performed to convert the problem to a conventional NLP form. In addition, the proposed approach extends beyond differential flatness to include the more general class of orbitally flat systems. This section continues by providing a brief introduction to both differential flatness and orbital flatness for nonlinear dynamic systems. In the next section, these concepts will be used to explicitly eliminate the differential equations from the dynamic optimization problem and transform the dynamic optimization problem into a normalized form. 2.1. Differential Flatness. Differential flatness is a term that applies to underdetermined systems of ordinary differential equations (ODEs). A general underdetermined ODE system of order k can be written as
f(t,x,x(1),...,x(k)) ) 0
(2)
where the elements of f ∈ R n-p are C∞ functions, x ∈ R n represents the dependent variables, t is the independent variable (usually time), x(k) stands for the kth time derivative of x with respect to t, and p (g1) is the number of equations by which the system is underdetermined. Differential flatness is then defined as follows:33 The system given by eq 2 is said to be differentially flat or simply flat if there exist variables y ) [y1, ..., yp]T given by an equation of the form
y ) h(t,x,x(1),...,x(m))
(3)
such that the original variables x can be recovered from y (locally) by an equation of the form
x ) g(t,y,y(1),...,y(l))
(4)
Further, the variables y ) [y1, ..., yp]T are referred to as the flat outputs. Differential flatness of a nonlinear system is very difficult to prove in general. There are no necessary and sufficient conditions for the existence of linearizing
outputs. As a result, one can only prove differential flatness when a carefully chosen set, or subset, of linearizing outputs can be found. For single-input systems, Shadwick40 showed that flatness is equivalent to static state-feedback linearizibility. In this special case, it is relatively easy to uncover the linearizing outputs by using the well-known GS algorithm.12 For systems where the number of states exceeds the number of inputs by one, Charlet et. al.4 demonstrated that dynamic feedback linearizability (i.e., flatness) is assured if the system is locally strongly accessible. In Martin and Rouchon,27 it was demonstrated that any driftless nonlinear systems with n states and n - 2 inputs is flat. In general, the multi-input case is difficult because one has to prove, at least partially, the dynamic feedback linearizability of the nonlinear system by endogenous feedback. To avoid this problem, Rathinam and Sluis32 showed that, for flat systems, the computation of flat outputs can be reduced to a one-dimensional problem. By a judicious choice of parametrization of the state space, multi-input integrability conditions can be reduced to a single integrability condition similar to the single-input situation. Application of the GS algorithm provides the remaining linearizing outputs. Unfortunately, this elegant result can be difficult to apply in situations where a suitable subset of flat outputs cannot be found. An alternative approach was developed by Guay et al.,14 who considered the problem of dynamic linearization by endogenous dynamic feedback. Using an exterior calculus approach, they developed a set of necessary and sufficient conditions for dynamic feedback linearizibility applicable to a large class of linearizable systems with arbitrary precompensators. Assuming a specific structure for the dynamic precompensator, these conditions yield an algorithm for identifying the linearizing outputs of control-affine systems. In general, the conditions for linearizability by static or dynamic state feedback can be quite restrictive. These restrictions can be relaxed somewhat by considering state-dependent time-scaling transformations along with state feedback and a local change of state coordinates. This problem is briefly introduced in the following section. 2.2.Time Scaling for Dynamic Systems. In the analysis of any continuous-time system, the time scale t is assumed to be specified and independent. For many analysis tasks, it is often cumbersome to restrict dynamics to evolve with respect to the actual time t. In a number of approaches, a new time scale τ is introduced to facilitate the analysis. This is the case in most singular-perturbation methods19 where the time scale of the slow dynamics are adjusted to facilitate stability analysis. Consider a smooth nonlinear system of the form
dx ) f(x,u) dt
(5)
where x ∈ M ⊂ R n and u ∈ U ⊂ R p. A general timescaling transformation τ can be defined by
dt t s(x,u), dτ
τ|t0 ) τ0
(6)
where t is the actual time. The function s(x,u) is assumed to be smooth such that 0 < s(x,u) < ∞ and is called a time-scaling function. Rewriting the original
2092
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
system given in eq 5 in the new time scale τ yields
dx t g(x,u) ) s(x,u) f(x,u) dτ
(7)
In eq 7, g(x,u) is well-defined because s(x,u) * 0. The restriction 0 < s(x,u) < ∞ is very important to ensure that new time τ is a strictly monotonically increasing function with respect to the actual time t. This guarantees that the stability and the state trajectory of the system are preserved by the transformation of time scale. Moreover, the time-scaling function s(x,u) is required to be C∞ with respect to x and u to preserve the smoothness of the vector fields. For single-input systems, necessary and sufficient conditions for the feedback linearization of control-affine nonlinear systems under state-dependent time-scaling functions were first derived by Sampei and Furuta.39 These conditions were later refined by Respondek34 and Guay.15 The algorithm presented by Guay15 identifies linearizing outputs for single-input control-affine nonlinear systems using an exterior calculus setting, which simplifies some aspects of the computations required. Recent work presented by Murray29 has considered the use of time-scale transformations to handle the problem of rate and magnitude saturations in trajectorygeneration and -tracking problems. The use of orbitalflatness transformations can serve as a valuable method of incorporating saturation problems and providing a more effective way of solving trajectory-generation problems. This is the topic of further research on the subject. In the following section, the flatness of nonlinear systems is exploited to transform the dynamic optimization problem into an algebraic optimization problem, which will be referred to as a normalized dynamic optimization (NDO) problem in this paper. 3. Normalized Form Optimization In this section, an algorithm for the solution of dynamic optimization problem 1 is presented that exploits the flatness of a system to transform the optimization problem into a problem formulation that is more easily solved.
Algorithm: Normalized Dynamic Optimization (NDO) 3.1. Identification of Flat Outputs. The flat outputs [y(t)] for the differential equation system x3 (t) ) f(x(t),u(t),z) are identified using the algorithm of Guay et al.14 or Rathinam and Murray.33 In case the system is not differentially flat in the original time scale t, the orbital flatness of the system is checked. For orbital flatness, a suitable time scaling is determined from the solution to eq 6, and if possible, the flat outputs [y(τ)] for the transformed system are identified. 3.2. System Variable Transformation. The original system variables x(t) and u(t) are written as functions of these flat outputs [y(t) or y(τ)] and a finite number of their derivatives
j (t)) x(t) ) r(y(t),y(1)(t),...y(k)(t)) t r(y j (t)) u(t) ) β(y(t),y(1)(t),...y(k)(t)) t β(y
(8)
for the case of differentially flat systems. Here, y(i)(t) stands for the ith derivative of y(t) with respect to t,
and k is the number of derivatives of y(t) required to represent the system in the form given in eq 8. The flat output y(t) and its derivatives up to order k have been collected in the vector y j (t). For the orbital flatness case that incorporates time scaling, the equivalent expressions take the form
x(t) ) r(y(τ),y(1)(τ),...y(k)(t)) t R(y j (τ)) u(t) ) β(y(τ),y(1)(τ),...y(k)(τ)) t β(y j (τ))
(9)
where the terminology is the same as in eq 8, except that y j is now a function of τ, y j (τ), instead of t. In this case, the new time scale τ is obtained as a function of t from the solution of eq 6. 3.3. Problem Transformation.1 The functions r and β are substituted for x(t) and u(t) in the original optimization problem. The normalized form of the optimization problem can then can be represented as
j (τ)),β(y j (τ))) min Φ(r(y y j subject to
g(r(y j (τ)),β(y j (τ))) e 0 c(r(y j (τ)),β(y j (τ))) ) 0 r(y j (τ0)) ) x0; τ ∈ [τ0, τf]
(10)
j (τ)) e xUB xLB e r(y uLB e β(y j (τ)) e uUB Note that problem 10 contains only algebraic constraints, yet is equivalent to problem 1. Problem 10 will be referred to as the normalized dynamic optimization (NDO) problem in the following discussions. 3.4. Parametrization. The flat outputs are parametrized by suitable functions of time (e.g., simple polynomials, Lagrange polynomials, etc.). In this paper, simple polynomials [i.e., y(τ) ) a + bτ + cτ2 + ...] are used to represent the flat outputs. Given an assumed structure for the flat outputs, the optimization problem is reduced to values for the parameter vector θ ) [aT, bT, cT, ...]T. Thus, problem 10 becomes
min Φ(θ, τ) θ subject to
fg(θ,τ) e 0 fc(θ,τ) ) 0 fx(τ0) ) x0; τ ∈ [τ0, τf]
(11)
xL e fx(θ,τ) e xU uL e fu(θ,τ) e uU The functions are not explicitly stated here for the purposes of brevity. In problem 11, fg(θ,τ) represents the transformed function g(x,u) in the inequality constraint of problem 1. The notation used for the other constraints is similar. The choice of polynomial order used to represent the flat outputs is dictated by the problem and should be made so that there is at least one degree of freedom for optimization. 3.5. Optimal Trajectory Determination. Determination of the optimal trajectories from problem 11 is a two-step procedure. In the first step, NDO problem 11 must be solved numerically. The second step consists of substituting the calculated values of the parameter
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2093
vector θ into the expression for the flat outputs y, which, in turn, are substituted into the expressions given in eqs 8 or 9 to determine the trajectories for the state and input variables. The optimization technique used in solving problem 11 depends on the presence of specific types of inequality constraints. Inequality constraints that represent restrictions on the path that trajectories might follow within the interval of optimization are often called path constraints. When there are no path constraints in problem 11, it is a conventional algebraic optimization problem that can be solved using standard LP, QP, or NLP codes, as appropriate. When path constraints are present in problem 11, it is a semi-infinite optimization problem. Good introductory treatments of semi-infinite optimization can be found in Polak31 and Hettich and Kortanek.16 Specifically, problem 11 is a one-parametric semi-infinite optimization problem,18 which is a subclass of semi-infinite optimization problems that is comparatively easy to solve. Example 3.1 provides a detailed illustration of the computations required in solving a dynamic optimization problem using the proposed method. Example 3.1. Parallel Reaction Problem. Consider the tubular chemical reactor problem that has been studied by several researchers2,7,21,22,37 in which the following reactions take place:
dx1 x1 1 ) - - x1 dτ u 2 dx2 x1 ) dτ u
(15)
A flat output y for this system is
y ) x1 + x2
(16)
Differentiating the flat output y is with respect to τ yields
dy dx1 dx2 1 ) + ) - x1 dτ dτ dτ 2 which can be rearranged to give
x1 ) - 2y˘
(17)
By substituting eq 17 into eq 16 and simplifying, we obtain
x2 ) y + 2y˘
(18)
Equation 17 or 18 can then be substituted into either expression in eq 15 and simplified to yield an expression for the input variable
A9 8B k 1
u)
8C A9 k
-2y˘ y˘ + 2y¨
(19)
2
Assuming that x1 t CA/CAf and x2 t CB/CAf, the system dynamics can be represented by
dx1 1 ) -ux1 - u2x1 dt 2 dx2 ) ux1 dt
(12)
where the control u ) k1L/v (L is the reactor length and v is the space velocity within the reactor). The dynamic optimization problem for this problem is
The new time scale τ is obtained as a solution to differential eq 14
[
x1(0) ) 1
max y(τf) + 2y˘ (τf) y(τ),y˘ (τ),y¨ (τ) y(0) ) 1 -7 y¨ (τ) g y˘ (τ), τ ∈[τ0, τf] 10
subject to
(21)
(13)
where τ is obtained as a function of t from the solution of eq 20 and the final time τf corresponds to t ) 1. For illustration purposes, suppose that y(τ) in problem 21 can be represented in quadratic form
y(τ) ) a + bτ + cτ2
x2(0) ) 0
(22)
Then, the expressions for its derivatives up to second order are
0 e u(t) e 5 This system is not differentially flat, but it can be shown to be orbitally flat with the time scaling
dt 1 ) s(x,u) ) 2, t|τ)0 ) 0 dτ u
(20)
y˘ (τ) e 0; τ ∈ [τ0, τf]
1 x˘ 1(t) ) -u(t)x1(t) - u2(t)x1(t) 2 x˘ 2(t) ) u(t) x1(t)
2
Substitution of the expressions for x1, x2, and u into problem 13 and simplification yields the NDO problem
max x2(tf), tf ) 1 u(t) subject to
]
y˘ (τ) + 2y¨ (τ) dt 1 ) ) dτ u2(t) - 2y˘ (τ)
(14)
The system, in the new time scale τ, can be represented as
y˘ (τ) ) b + 2cτ y¨ (τ) ) 2c
(23)
and the time scaling can be determined via the analytical solution of eq 20 to be
t)-
2c τ + - ln(b + 2cτ) 4 b + 2cτ
(24)
2094
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
Substituting eqs 22-24 into the NDO problem 21 yields 2 max a + (2 + τf)b + (4τf + τf )c a,b,c
subject to
a)1 7 7 b+21+ τcg0 10 10 b + 2τc e 0
(
)
(25)
2c - ln(b) ) 0 b -
τf 2c - ln(b + 2cτf) ) 1 + 4 b + 2cτ f
where the final two equations in problem 25 represent the initial and final conditions on the solution to the time-scaling problem given in eq 14. Problem 25 is a NLP problem in two decision variables, as the parameter a is required to be one, and can be solved using available nonlinear semi-infinite optimization codes. If the system had been differentially flat, then problem 25 would reduce to a simple linear semi-infinite optimization problem. The solution of problem 25 required solving a differential equation for the time scaling, which was the result of the system’s orbital flatness. Differentially flat systems will not require the solution of any such differential equations. In this example, the differential equation could be solved analytically; however, this will not be the general case. When an analytical solution is not possible, a numerical differential equation solver must be embedded within the optimization problem formulation. The additional computational load that results from embedding a numerical differential equation solver within the problem formulation is comparatively small, as the time-scaling differential equation is scalar. 4. Case Studies In this section, a number of benchmark case studies taken from the literature are investigated. These case studies are intended to provide a comparison of the proposed dynamic optimization method with the more established approaches. The comparison will be made in terms of the accuracy of the solution and the size of the optimization problem to be solved. 4.1. Single-Integrator System. This is the singleintegrator system studied in Goh and Teo,13 Luus,23 and Dadebo and McAuley.7 The optimization problem for this system can be written as
min x2(tf), tf ) 1 u(t) subject to
x˘ 1(t) ) u(t) x˘ 2(t) ) x12(t) + u2(t) x1(0) ) 1 x2(0) ) 0 x1(1) ) 1
(26)
Figure 1. Trajectory comparison for the single-integrator problem.
One choice of flat output for this system is y(t) t x1(t). Note that, because the initial condition for x2(t) is known, an expression for x2(t) can be obtained as an integral of a function of the flat output and its first derivative. Then, the system variables can be expressed in terms of the flat output y(t) and its derivatives as follows:
x2(t) )
x1(t) t y1(t)
(27)
∫0t[y2(t) + y˘ 2(t)] dt
(28)
u(t) ) y˘ (t)
(29)
Parametrizing the flat output y(t) as
y(t) ) a + bt + ct2 + dt3 + et4
(30)
and collecting the coefficients of the parametrization into the vector θ t [a, b, c, d, e]T enable the dynamic optimization problem 26 to be written as
subject to where
[
T min θ Hθ θ 0 0 0 1 1 1
[
1 1 /2 1 H ) /3 1 /4 1 /5
1 1
1
/2 4 /3 5 /4 6 /5 7 /6
1
/3 5 /4 23 /15 5 /3 61 /35
] []
0 1 θ) 1 1
1
/4 6 /5 5 /3 68 /35 17 /8
1
/5 /6 61 /35 17 /8 151 /63
7
]
(31)
Note that problem 26 does not contain any path constraints and is an equality-constrained QP, which can easily be solved either analytically or using a number of readily available optimization codes. The solution to problem 31 is θ ) [1, -0.4621, 0.4996, -0.0749, 0.0375]T, which yields an objective function value of x2(1) ) 0.9242343146. The analytical solution to this problem is identical to 10 significant figures. Figure 1 gives a comparison between the state and input
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2095
trajectories determined using the optimal parameter values in eqs 27-29 and the trajectories determined analytically. Note that there is no appreciable difference between the two solutions. The analytical solution can be found by reformulating problem 26 as a finite-time optimal control problem. Problem 26 is equivalent to
min u(t)
x1(0) ) 1
( x
(32)
max x2(tf), tf ) 1 u(t) 1 x˘ 1(t) ) -u(t) x1(t) - u2(t) x1(t) 2 x1(0) ) 1
(33)
where the manipulated variable is defined as u(t) t k1L/v (L is the reactor length and v is the space velocity), the first state variable is defined as x1 t CA/ CAf, and the second state variable is defined as x2 t CB/ CAf. The need to solve the time-scaling ODE, which was presented in example 3.1 for a discussion of the orbital flatness of this system, can be eliminated in this case by selecting the linearizing output to be y(t) t x1(t). (Note the different choice of flat output used here in comparison to that used in example 3.1.) The system variables can be expressed in terms of the linearizing output y(t) and its derivatives as follows. First, rewrite the expression
(34)
in terms of y(t) and y˘ (t)
Hence
(35)
(38)
and substituting into eq 37 yields
( x
)
b + 2ct 1-2 (a + bt + ct2) a + bt + ct2
x˘ 2(t) ) -1 +
Integrating using the initial condition x2(0) ) 0 gives
x2(t) )
( x
∫0t -1 +
1-2
)
b + 2cσ a + bσ + cσ2 (a + bσ + cσ2) dσ
which can be solved by quadrature. Collecting the coefficients into the vector θ t [a, b, c]T enables the dynamic optimization problem 33 to be written as
( x
∫01 -1 +
subject to
b + 2cσ 1-2 a + bσ + cσ2
) (a + bσ + cσ2) dσ
a)1
( x
0 e -1 +
0 e u(t) e 5
y˘ (t) 1 2 )0 u (t) + u(t) + 2 y(t)
(37)
Representing the linearizing output y(t) in quadratic form
max θ
x2(0) ) 0
1 x˘ 1(t) ) - u2(t) x1(t) - u(t) x1(t) 2
)
y(t) t a + bt + ct2
Problem 32 can be solved via least-squares optimal control theory.3 This classical approach requires the solution of a Riccati differential equation. Note that the proposed method converted an optimal control problem into an equality-constrained QP in only five decision variables (i.e., the coefficients in the polynomial representation for the flat output), which yielded very high solution accuracy despite the low-order approximation. Thus, this example has some interesting implications for the solution of linear optimal control problems. 4.2. Parallel Reaction Problem. A detailed description for this tubular reactor problem is given in example 3.1, in which the optimization problem is given as
x˘ 2(t) ) u(t) x1(t)
(36)
y˘ (t) 1-2 y(t) y(t)
x˘ (t) ) -1 +
x1(1) ) 1
subject to
x
y˘ (t) 1-2 y(t)
Substituting eq 36 into the ODE for the second state yields
∫01x12(t) + u2(t) dt x˘ 1(t) ) u(t)
subject to
u(t) ) - 1 +
)
(39)
b + 2ct 1-2 e5 a + bt + ct2
Note that problem 33 does contain input path constraints and is a semi-infinite optimization problem. The objective function value obtained for the optimal solution to this problem is x2(1) ) 0.57189, and the optimal polynomial coefficient values are θ ) [1, -1.03180, 0.081453]T. The optimal state and input trajectories determined using this technique are shown in Figure 2. The simple, quadratic approximation produced accurate results with very few computations. Also, the solution to the NDO problem 39 required only the integration of the objective functional and was computationally inexpensive. Increasing the order of the polynomial used to represent the flat output to tenth order produced an optimum solution with an objective function value of x2(1) ) 0.57304 and optimal polynomial coefficients of
θ ) [1, -1.06433, 0.15822, 0.022151, 0.15497, -0.097468, -0.16193, -0.093541, -0.19312, -0.10974, 0.41606]T which are slightly less than the reported literature values.
2096
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
A flat output for the transformed system can be defined as y(τ˜ ) ) x1(τ˜ ) + x2(τ˜ ). Then, representing the output y(t˜) in cubic form
y(τ˜ ) t a + bτ˜ + cτ˜ 2 + dτ˜ 3
(44)
allows the system variables to be expressed as
x1(τ˜ ) t 0.032258
(
)
1600cτ˜ + 2400dτ˜ 2 + 31τf + 31τfcτ˜ 2 + 31τfdτ3 τf x2(τ˜ ) ) -25.80645
Figure 2. Optimal trajectories for the parallel reaction problem.
4.3. Consecutive Reaction Problem. Consider the batch reactor problem given in Ray,37 which includes the consecutive reaction
A9 8B9 8C k k 1
2
The operating objective for this reactor is to obtain the reactor temperature progression that maximizes the intermediate product for a fixed batch time. The optimization problem for the system can be stated as
max x2(tf), tf ) 1 u(t) subject to
u(τ˜ ) ) -
[
]
τ˜ (2c + 3dτ˜ ) τf
x1(τ˜ )2 x˘ 1(τ˜ )
Note that b ) 0 ensures that x2(0) ) 0 and a ) 1 ensures that x1(0) ) 1. Collecting the remaining coefficients into the vector θ t [c, d]T enables the dynamic optimization problem 33 to be written as
min -25.8064 θ subject to
(
τf dt t 2 , t|τ˜)0 )0, t|τ˜)1 )1 dτ˜ u (τ˜ ) 0.9092 e -
x˘ 2(t) ) -0.03875u2(t) x2(t) + u(t) x12(t)
)
2c + 3d τf
x1(τ˜ )2
x˘ 1(t) ) -u(t) x12(t)
(45)
x˘ 1(τ˜ )
(46)
e 7.4831
0 e τ˜ e 1
(40)
x1(0) ) 1 x2(0) ) 0 0.9092 e u(t) e 7.4831
where the manipulated variable is defined as u t 4000e-2500/T and the state variables are defined as x1 t cA and x2 t cB, respectively. This system can be shown to be orbitally flat using the time-scaling transformation
dt 1 , t(0) ) 0 t dτ u2(τ)
(41)
To synchronize the two time scales, an additional linear change of time scale can be defined as
τ˜ )
τ τf
(42)
such that the new time scale τ˜ is defined on [0, 1]. Then, eq 41 can be re-expressed as
τf dt t 2 , t(0) ) 0 dτ˜ u (τ˜ )
(43)
Note that the transformed optimization problem remains a dynamic optimization problem; however, the transformed problem requires the integration of only one differential equation, which is needed for the synchronization of the two time scales. The computational advantage of the proposed approach is that it requires the integration of only a scalar ordinary differential equation, regardless of the underlying dimension of the optimization problem. Problem 46 can be solved via semi-infinite optimization; however, the final time in the new time scale (τf) becomes an additional decision variable in the optimization problem. The optimal objective function value for the solution is x2(1) ) 0.61019, and the optimal parameter values are θ ) [1, 0, -0.12919, 0.03883]T and τf ) 6.00053. The calculated system trajectories are shown in Figure 3. 4.4. Crane Container Problem. In this complex problem, containers are to be optimally transferred from a ship to a cargo truck using a crane. This problem has posed significant challenges in the optimal control literature and has been considered in Sakawa and Shindo,38 Goh and Teo,13 Luus,23 Teo et al.,43 and Dadebo and McAuley.7 The optimization problem for this system can be represented as
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2097
x5(t) )
y˘ 1(t) 9
x6(t))
y˘ 2(t) 9
u1(t)) -y2(t) y˘ 2(t) - 27.0756y1(t) u2(t)) x4(t) )
∫0t9u1(σ) + 9(17.25656)x3(σ) dσ x1(t) )
Figure 3. Optimal trajectories for the consecutive reaction problem.
∫0t9x4(σ) dσ
Note that this system is not flat, as defined in Martin and Rouchon.27 The system can be considered partially flat by concentrating on the subsystem
x˘ 2(t) ) 9x5(t)
∫01[x32(t) + x62(t)] dt
9 min u(t) 2
x˘ 3(t) ) 9x6(t)
(50)
x˘ 5(t) ) 9u2(t)
subject to x˘ 1(t) ) 9x4(t) x˘ 2(t) ) 9x5(t)
x˘ 6(t) )
x˘ 3(t) ) 9x6(t) x˘ 4(t) ) 9[u1(t) + 17.25656x3(t)] x˘ 5(t) ) 9u2(t) x˘ 6(t) )
y˘ 1(t) 9
2y˘ 1(t) y˘ 2(t) 81 (49)
-9[u1(t) + 27.0756x3(t)] x2(t)
-
18x5(t) x6(t) x2(t) x(0) ) [0 22 0 0 -1 0]T
18x5(t) x6(t) x2(t)
y1(t) t a0 + a1t + a2t2 + ... + a8t8
(47)
y2(t) t b0 + b1t + b2t2 + ... + b8t8
x(1) ) [10 14 0 2.5 0 0] -2.834 e u1(t) e 2.834 -0.809 e u2(t) e 0.713 |x4(t)| e 2.5 |x5(t)| e 1.0
where u1(t) and u2(t) are the torques of the hoist motor and the trolley drive motor, respectively. The objective function represents the swing of the container during the transfer and is to be minimized for safety reasons. This problem is normalized with the following choice of flat outputs
y2(t) t x3(t)
x2(t)
-
which is flat and independent of x1(t) and x4(t). The state variables x1(t) and x4(t) can be recovered by integrating appropriate expressions in the other state and input variables. In this case, it is still possible to express the optimization problem in an algebraic form by choosing an appropriate parametrization of the output trajectories. For this problem, the two chosen outputs are parametrized as high-order polynomials
T
y1(t) t x2(t)
-9[u1(t) + 27.0756x3(t)]
Then, the integrals for x1(t) and x4(t) can be solved analytically by symbolic computation. To solve this problem, the coefficients can be collected into the vector θ t [a0, a1,..., a8, b0, b1,..., b8]T. From the final and initial time constraints, a0 ) 22, a1 ) -9, b0 ) 0, and b1 ) 0 to ensure that x1(0) ) 22, x5(0) ) -1, x3(0) ) 0, and x6(0) ) 0. The dynamic optimization problem is transformed to a semi-infinite optimization problem of the general form T min θ Hθ θ
subject to
w(θ) ) 0 p(θ,t) e 0
(52)
t ∈ [0,1] (48)
The state and input variables expressed in terms of the flat outputs are
(51)
where
H)
[
0 0
0 H1
]
[
2098
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
H1 ) 0 0 0 263/270
0 5 /6
0 461
283
25
197
905
5
26
97
13
23
160
61
0 0 0 0 0 0
and
[
/6 461
/35
0 /630
/432
97
/144
79
/126
/21
0 /42
/40 331 /594 283 502 /432 13/21 53/90 /891 13/24 25 23 151 /42 /40 331/594 13/24 /286 197 160 191 673 /360 /297 /360 /1287 65/126 905 /1782 61/120 1307/2574 383/756 197/390 /630
/144
0
53
/90
[
0 /360 /297
/1782
/120 1307 /2574 673 /1287 383/756 65 /126 197/390 145 298 /288 /585 145 /288 2303/4590
191
/360
x1(y1(1),y2(1),y˘ 1(1),y˘ 2(1)) y1(1) y (1) w(θ) ) 2 x4(y1(1),y2(1),y˘ 1(1),y˘ 2(1)) y˘ 1(1)/9 y˘ 2(1)/9
]
]
x˙ 1 ) a1q(x1)(x1in - x1) - [kp(x3) + kttM(x3)]h(x1,x2,x3)x1 - b1kd(x3)x2 x˘ 2 ) a1q(x1)(u1 - x2) - kd(x3)x2
]
Solving this semi-infinite optimization problem yields a value for the objective function of 0.0058995, which is slightly larger than the one reported in the literature (e.g., Goh and Teo13). The optimal coefficient values are
[][ ] 0 0 -0.588425 1.31594 1.31569 0.439151 -0.0892923 1.07472 -0.836412
min tfinal u1(t),u2(t) subject to
p(θ,t) ) {-y2(t) y˘ 2(t) - 27.0756y1(t) - [2y˘ 1(t) y˘ 2(t)]/81} - 2.834 -{-y2(t) y˘ 2(t) - 27.0756y1(t) - 2y˘ 1(t) y˘ 2(t)/81} - 2.834 (y¨ 1(t)/9) - 0.713 -(y¨ 1(t)/9) - 0.809 x4(y1(t),y2(t),y˘ 1(t),y˘ 2(t)) - 2.5 -x4(y1(t),y2(t),y˘ 1(t),y˘ 2(t)) - 2.5 (y˘ 1(t)/9) - 1 -(y˘ 1(t)/9) - 1
22 -9 1.69126 -14.9895 [ai] ) 45.8733 , [bi] ) -34.2773 -54.0050 91.1489 -34.4415
where, for other methods, accurate initial estimates of the entire input trajectories were required to arrive at a satisfactory answer. 4.5. Polymerization Reactor. In this example, trajectory generation for a polymerization reactor is considered. Using the model given in ref 44, we attempt to compute a trajectory that would steer the styrene polymerization reactor between two equilibrium points in minimum time. The optimization problem for this system is given by
x˘ 3 ) a1q(x1)(c1x3in - x3) + d1q(x1)(u2 - x3) - e1kp(x3) h(x1,x2,x3)x1 h(x1,x2,x3) ) kdi(x3)fCI CI )
x2 q(x1)Mi
q(x1) )
x1 Ws 1 - x1 - Ws + + Fm Fs Fp
kd(x3) ) 1.58 × 1015e-15395.71807/x3 ki(x3) ) 0.002184e-13832.09045/x3 kp(x3) ) 10510.0e-3553.043060/x3
(54)
kttM(x3) ) 2310.0e-6374.789512/x3 kttS(x3) ) 558000.0e-8636.035603/x3 ktc(x3) ) 627500.0e-845.0805870/x3
(53)
The optimal state and input trajectories are shown in Figure 4. These results differ from those reported by Dadebo and McAuley,7 who used relaxations of the final state constraints to allow solution of the problem. Thus, the solutions reported in the literature do not strictly satisfy all of the constraints. The lower values of the objective functions reported by other researchers result from the infeasibility of their solutions. Using the coefficient values given in eq 53, the final states are x(1) ) [10, 14, 2.0 × 10-9, 2.5, - 1.0 × 10-7, 0]T. In addition, Figure 4 shows that the path constraints were respected. Thus, the proposed method produced a feasible solution to problem 47 that has an objective function value that is very close to the results reported in the literature. The resulting semi-infinite optimization problem was solved using initial values for the coefficients of zero. This result contrasts drastically with the literature,
ktd(x3) ) 0 kdi(x3) ) 50178.96259e-7275.318740/x3 x1in ) 0.9, x3in ) 298.15 K FM ) 830, Fs ) 790, Fp ) 1025 a1 ) 0.1, Mi ) 0.164, b1 ) 0.76098,
e1 ) -385.652
c1 ) 1.0663, d1 ) 0.32345, Ws ) 0.1, f ) 0.6 x1(0) ) 0.68, x2(0) ) 0.00123 x1(tfinal) ) 0.785, x2(0) ) 2.67 × 10-4 0.0 e u1(t) e 0.001 298.0 e u2(t) e 400.0 The dynamics that govern this process are differentially flat with flat outputs y1 ) x1 and y2 ) x2. Exploiting the flatness of the dynamics, the optimization problem 54 can be represented in normalized form by parametrizing the flat output paths as
[
a a a a a y ) a11 a12 a13 a14 a15 21 22 23 24 25 ) Θτ
[]
]
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2099
1 t t2 t3
where y ) [ y1 y2 ]T. The normalized optimization problem is
min tf Θ subject to
y1(t) ) Θτ y1(0) ) 0.68, y1(tfinal) ) 0.785
(55)
y2(0) ) 0.00123, y2(0) ) 2.67 × 10-4 0.0 e u1(t) e 0.001 298.0 e u2(t) e 400.0 Expressions for the variables x3, u1, and u2, in terms of the coefficients in the parametrization of the flat outputs, can be determined by substituting the expressions for the flat outputs and solving the following set of nonlinear algebraic equations
y˘ 1 ) a1q(x1)(x1in - x1) - [kp(x3) + kttM(x3)]h(x1,x2,x3)x1 - b1kd(x3)x2 y˘ 2 ) [a1q(x1)(u1 - x2) - kd(x3)x2] y¨ 1 ) γ1y˘ 1 + γ2y˘ 2 + γ3[a1q(x1)(c1x3in - x3) + d1q(x1)(u2 - x3) - e1kp(x3) h(x1,x2,x3)x1]
Figure 4. Optimal trajectories for the crane-container problem.
4.6. Bioreactor Problem. This bioreactor example is taken directly from Mahadevan et al.24 and modified slightly to
min tf u(t)
where
∂q (x - x1) - a1q(x1) - [kp(x3) + ∂x1 1in
γ1 ) a1
subject to
∂h x kttM(x3)]h(x1,x2,x3) - [kp(x3) + kttM(x3)] ∂x1 1 ∂h γ2 ) -[kp(x3) + kttM(x3)] x - b1kd(x3) ∂x2 1
(
γ3 ) -
(Km + x1)YX/S
x˘ 2(t) ) u
)
x1(0) ) 0 x2(0) ) 2
∂kd ∂h kttM(x3)) x - b1 x ∂x3 1 ∂x3 2 Solution of the semi-infinite optimization problem 55 yields coefficient values for the flat outputs of
[
µmaxx1ξ1
ξ1 ) C/x2 + YX/S(Sf - x1)
∂kp ∂kttM + h(x1,x2,x3)x1 - (kp(x3) + ∂x3 ∂x3
Θ) 0.68 4.70 × 10-10 3.48 × 10-9 -1.50 × 10-13 0.00123 -2.80 × 10-7 2.73 × 10-11 -1.09 × 10-15
x˘ 1(t) ) -u(t)(Sf - x1)/x2 -
]
with a transition time of tfinal ) 6465.5 s. The problem, solved using Matlab, required less than 60 s on 266 MHz Intel Pentium II workstation. The state trajectories are all smooth and monotonic, as shown in Figure 5. The input constraints were all strictly satisfied, as shown in Figure 6.
ξ1(0) ) 1 C ) x2(0){ξ1(0) - YX/S[Sf - x1(0)]} µmax ) 0.5, YX/S ) 0.4, Sf ) 15 Km ) 1.2 0 e u(t) e 10 x2(t) e 10 ξ1(tf) x2(tf) g 50 A flat output for this system is
y1 ) (Sf - x1)x2
(56)
2100
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
We parametrize the flat output as a polynomial of the form
[]
y1(t) ) a1 + b1t + c1t2 + d1t3 + e1t5 1 t 2 ) [a1 b1 c1 d1 e1] t t3 t5 ) θTτ
and restate the optimization problem in normalized form
min tf θ Figure 5. Optimal state trajectories for the styrene reactor problem.
subject to
y1(t) ) θTτ x1(t) ) f1(y1(t), y˘ 1(t)) x2(t) ) f2(y1(t), y˘ 1(t)) u1(t) ) f3(y1(t), y˘ 1(t), y¨ 1(t)) ξ1(t) ) C/x2(t) + YX/S[Sf - x1(t)] x1(0) ) 0 x2(0) ) 2
(57)
x1(0) ) 1 C ) x2(0){ξ1(0) - YX/S[Sf - x1(0)]} µmax ) 0.5, YX/S ) 0.4, Sf ) 15 Km ) 1.2 0 e u(t) e 10 x2(t) e 10 Figure 6. Optimal input trajectories for the styrene reactor problem.
Simple algebraic manipulations show that the states and input of this system can be written as
x1 )
y˘ 1YX/SKm (µmaxC + µmaxYX/Sy1 - y1YX/S)
) f1(y1,y˘ 1)
x2 ) (µmaxC + µmaxYX/Sy1 - y˘ 1YX/S) y1 ) (µmaxCSf + µmaxSfYX/Sy1 - y˘ 1YX/SKm - y˘ 1YX/SSf) f2(y1,y˘ 1) u1 ) (y˘ 1Sfµmax2C2 + 2Cy˘ 1YX/SSfy1µmax2 2Cµmaxy˘ 12YX/SSf + Cy1(3)µmaxKmy1YX/S Cµmaxy˘ 12YX/SKm + y˘ 1YX/S2y12Sfµmax2 + µmaxy¨ 1Kmy12YX/S2 - 2µmaxy˘ 12YX/S2y1Sf 2µmaxy˘ 12YX/S2Kmy1 + y˘ 13YX/S2Km + y˘ 13YX/S2Sf)/ ((µmaxCSf + µmaxSfYX/Sy1 - y˘ 1YX/SKm - y˘ 1YX/SSf)2) ) f3(y1, y˘ 1, y¨ 1)
ξ1(tf) x2(tf) g 50 Solving SIO problem 57 yields the solution
θT )
[30.0, 0.0, 0.666428, 0.0317631, -0.000182772]
tf ) 14.1548. The resulting input and state trajectories are shown in Figure 7. The problem, solved using Matlab, required approximately 30 s on a 266-MHz Intel Pentium II computer. Initial values were randomly generated from a normal distribution with a mean of 0 and a standard deviation of 0.0001. The fifth-order polynomial chosen provided accurate results that were robust to changes in the initial values. The input and state constraints were all strictly satisfied, as shown in Figure 7. 4.7. Summary. Table 1 compares the results reported in the literature with those obtained using the proposed NDO approach for the first four benchmark examples taken from the literature. All problems were solved within Matlab, using the appropriate routines from the optimization toolbox. Selection of the order of the polynomial used to parametrize the flat outputs depends primarily on two
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2101
Figure 7. Optimal trajectories for the fed-batch bioreactor problem. Table 1. Comparison of Objective Function Values for the Illustrative Examples proposed analytical CVI CVP IDP IDP collocation
integrator
parallel
series
crane
0.92423 0.92423 0.9251813 0.924287 0.9244123 -
0.57304 0.573492 0.569102 0.573537
0.610189 0.6107757
0.5735321
0.6107752
0.0058995 0.0054013 0.004237 0.0050223 -
factors: (1) the problem constraints and (2) the required accuracy of the solution. The approach adopted for selecting the parametrization was to first determine the lowest-order polynomial that provided degrees of freedom for optimization. Then, an iterative procedure was employed whereby the polynomial order was systematically increased, and the resulting NDO problem was solved until the desired solution accuracy was obtained. In general, the proposed method produced results that were similar to or better than those reported in the literature. Once the appropriate polynomial order was determined for parametrizing the flat outputs, the proposed approach showed itself to be computationally efficient. Solution times for the crane container problem were approximately 30 s using Matlab on a 266-MHz Intel Pentium II computer, whereas Dadebo et al.7 reported a solution time of approximately 14 000 s using compiled Fortran on a 100-MHz Intel 486 computer. Comprehensive benchmarking of the various approaches to dynamic optimization is beyond the scope of this work; however, it appears from the case studies that the proposed approach might prove useful in on-line trajectory optimization problems because of its computational efficiency. 5. Conclusions This paper presented a new approach for calculating optimal trajectories for flat and partially flat systems. The proposed approach produced results that were similar to or better than those obtained using currently available methods based on discretization schemes or dynamic programming for a range of benchmark problems. Although the benchmark problems solved in this
paper were all fixed-time problems, the method can be used to solve minimum- or variable-time problems. In the consecutive reaction problem of section 4.3, a fixedtime problem was converted into a variable-time problem in the new time scale and solved using the proposed method. The power of the proposed approach arises from its transformation of a dynamic optimization problem into a simpler and lower-dimensional form. It appears that the method can produce arbitrary precision for sufficiently smooth problems, with small additional computational cost. Note that the computational complexity for the NDO approach scales with the number of coefficients in the expressions used to represent the flat outputs and not with the underlying model size or the fineness of the discretization, whereas for other dynamic optimization methods, computational complexity scales with model size, discretization fineness, and so forth. Thus, the proposed approach might prove to be computationally tractable for larger-scale dynamic optimization of flat and partially flat systems. The computational advantages offered by the proposed approach arise because the method exploits the geometry of the system’s dynamics to provide a coordinate transformation that greatly simplifies the optimization calculations. The determination of this transformation requires a substantial analysis stage prior to optimization. Thus, the NDO approach requires more pre- and post-optimization analyses than the conventional approaches to gain the computational advantages that are promised by the method. Furthermore, this additional pre-optimization analysis is similar to that required for controller design in a differential geometric framework, and the results could be used to formulate the trajectory-tracking controller that would implement the optimization results. Thus, the NDO approach offers promise for real-time applications. However, the problem that must be addressed is ensuring that a trajectory-tracking controller can be synthesized to enforce the calculated optimal trajectories. This is true regardless of the method used to determine such trajectories. Acknowledgment The authors gratefully acknowledge the financial assistance of the Natural Sciences and Engineering Research Council of Canada. Literature Cited (1) Agrawal, S. K.; Faiz, N.; Murray, R. M. Feasible Trajectories of Linear Dynamics with Inequality Constraints Using HigherOrder Representations; Proceedings of the IFAC World Congress, Beijing, China, 1999; International Federation of Automatic Control: Kidlington, Oxford, U.K., 1999. (2) Biegler, L. T. Solution of dynamic optimization problems by successive quadratic programming and orthogonal collocation. Comput. Chem. Eng. 1984, 8 (3/4), 243-248. (3) Brockett, R. W. Finite Dimensional Linear Systems; Wiley: New York, 1969. (4) Charlet, B.; Le´vine, J.; Marino, R. Sufficient Conditions for Dynamic State Feedback Linearization. SIAM J. Control Optim. 1991, 29, 38-57. (5) Crandall, M. G.; Evans, L. C.; Lions, P. Some properties of viscosity solutions of Hamilton-Jacobi Equations. Trans. Am. Math. Soc. 1984, 282, 487-502. (6) Cuthrell, J. E.; Biegler, L. T. Simultaneous optimization and solution methods for batch reactor control profiles. Comput. Chem. Eng. 1989, 13 (1/2), 49-62.
2102
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
(7) Dadebo, S. A.; McAuley, K. B. Dynamic optimization of constrained chemical engineering problems using dynamic programming. Comput. Chem. Eng. 1995, 19 (5), 513-525. (8) Faiz, N.; Agrawal, S. K.; Murray, R. M. Differentially Flat Systems with Inequality Constraints: An Approach to Real-Time Feasible Trajectory Generation. J. Guid., Nav. Control 2000, in press. (9) Fliess, M.; Le´vine, J.; Martin, P.; Rouchon, P. Line´arisation par Bouclage Dynamique et Transformations de Lie-Ba¨cklund. C. R. Se´ ances Acad. Sci. 1993, 317 (10), 981-986. (10) Fliess, M.; Le´vine, J.; Martin, P.; Rouchon, P. Flatness and Defect of Nonlinear Systems: Introductory Theory and Examples. Int. J. Control 1995, 61, 1327-1361. (11) Fliess, M.; Le´vine, J.; Martin, P.; Ollivier, F.; Rouchon, P. A remark on nonlinear accessibility conditions and infinite prolongations. Syst. Control Lett. 1997, 31, 77-83. (12) Gardner, R. B.; Shadwick, W. F. The GS Algorithm for Exact Linearization to Brunovsky Normal Form. IEEE Trans. Autom. Control 1992, 37, 224-230. (13) Goh, C. J.; Teo, L. K. Control Parametrization: A Unified Approach to Optimal Control Problems with General Constraints. Automatica 1988, 24, 3-18. (14) Guay, M.; McLellan, P. J.; Bacon, D. W. Conditions for Dynamic Feedback Linearization of Control-Affine Nonlinear Systems. Int. J. Control 1997, 68, 87-106. (15) Guay, M. An Algorithm for Orbital Feedback Linearization of Single-Input Control-Affine Nonlinear Systems. Syst. Control Lett. 1999, 38, 271-281. (16) Hettich, A.; Kortanek, K. Semi-Infinite Programming: Theory, Method and Applications. SIAM Rev. 1993, 35, 380-429. (17) Hicks, G. A.; Ray, W. H. Approximation methods for optimal control synthesis. Can. J. Chem. Eng. 1971, 49, 522-528. (18) Jongen, H. Th.; Stein, O. On Generic One-Parametric Semi-Infinite Optimization. SIAM J. Optim. 1997, 7 (4), 11031137. (19) Kokotovic, P. V.; Khalil, H. K.; O’Reilly, J. Singular Perturbation Methods in Control: Analysis and Design; Academic Press: New York, 1986. (20) Kravaris, C.; Chung, C.-B. Nonlinear State-Feedback Synthesis by Global Input/Output Linearization. AIChE J. 1987, 33, 592-603. (21) Logsdon, J. S.; Biegler, L. T. Accurate solution of differential algebraic equations. Ind. Eng. Chem. Res. 1989, 28, 1628-1639. (22) Logsdon, J. S.; Biegler, L. T. Decomposition Strategies for Large-Scale Dynamic Optimization Problems. Chem. Eng. Sci. 1992, 47, 851-864. (23) Luus, R. Application of iterative dynamic programming to state constrained optimal control problems. Hung. J. Ind. Chem. 1991, 245-254. (24) Mahadevan, R.; Agrawal, S. K.; Doyle, F. J., III. A Flatness Based Approach to Optimization in Fed-Batch Reactors. Proceedings of ADCHEM 2000, Pisa, Italy, June 14-16, 2000; International Federation of Automatic Control: Kidlington, Oxford, U.K., 1999; pp 111-116. (25) Martin, P. Contributions a` l′e´tude des syste`mes diffe´rentiellements plats. Ph.D. Thesis, Ecole des mines de Paris, Paris, France, 1992. (26) Martin, P.; Devasia, S.; Paden, B. A different look at output tracking control of a vtol aircraft. Proceedings of the 33rd IEEE Conference on Decision and Control; Institute of Electrical and Electronics Engineers: New Brunswick, NJ, 1994; pp 2376-2381. (27) Martin, P.; Rouchon, P. Any (controllable) driftless system with m inputs and m + 2 states is flat. Proceedings of the IEEE Conference on Decision and Control, New Orleans, LA, 1995; Institute of Electrical and Electronics Engineers: New Brunswick, NJ, 1995.
(28) Murray, R. M.; Rathinam, M.; Sluis, W. M. Differential Flatness of Mechanical Control Systems. Proceedings of the ASME International Mechanical Engineering Congress and Exposition; American Society of Mechanical Engineers: New York, 1995. (29) Murray, R. M. Geometric Approaches to Control in the Presence of Magnitude and Rate Saturations; CDS Technical Report 99-001; California Institute of Technology: Pasadena, CA, 1999. (30) Oldenburg, J.; Marquardt, W. Dynamic Optimization Based on Higher Order Differential Model Representations. Proceedings of ADCHEM 2000, Pisa, Italy, June 14-16, 2000; International Federation of Automatic Control: Kidlington, Oxford, U.K., 1999; pp 809-814. (31) Polak, E. Optimization Algorithms and Consistent Approximations; Springer: New York, 1997. (32) Rathinam, M.; Sluis, W. M. A Test for Differential Flatness by Reduction to Single Input Systems; Technical Report CDS95018; California Institute of Technology: Pasadena, CA, 1995. (33) Rathinam, M.; Murray, R. M. Differential Flatness of Two One-Forms in Arbitrary Number of Variables. Syst. Control Lett. 1999, 36, 317-326. (34) Respondek, W. Orbital Feedback Linerization of SingleInput Nonlinear Control Systems. Proceedings of the IFAC NOLCOS (Nonlinear Control Systems Design Symposium), Enschede, The Netherlands, July 1-3, 1998; International Federation of Automatic Control: Kidlington, Oxford, U.K., 1998; pp 499-504. (35) Rothfuss, R.; Rudolph, R. J.; Zeitz, M. Flatness Based Control of a Nonlinear Chemical Reactor Model. Automatica 1996, 32, 1433-1439. (36) Rouchon, P.; Fliess, M.; Le´vine, J.; Martin, P. Flatness, Motion Planning and Trailer Systems. Proceedings of the 32nd IEEE Conference on Decision and Control, San Antonio, TX, 1993; Institute of Electrical and Electronics Engineers: New Brunswick, NJ, 1993. (37) Ray, W. H. Advanced Process Control; McGraw-Hill: New York, 1981. (38) Sakawa, Y.; Shindo, Y. Optimal Control of Container Cranes. Automatica 1982, 18, 257-266. (39) Sampei, M.; Futura, K. On Time Scaling for Nonlinear Systems: Application to Linearization. IEEE Trans. Autom. Control 1986, 31, 459-462. (40) Shadwick, W. F.; Absolute Equivalence and Dynamic Feedback Linearization. Syst. Control Lett. 1991, 15, 35-39. (41) van Nieuwstadt, M.; Murray, R. M. Approximate Trajectory Generation for Differentially Flat Systems with zero Dynamics. Proceedings of the IEEE Conference on Decision and Control; Institute of Electrical and Electronics Engineers: New Brunswick, NJ, 1995; pp 4224-4230. (42) van Nieuwstadt, M.; Rathinam, M.; Murray, R. M. Differential Flatness and Absolute Equivalence of Nonlinear Control Systems. SIAM J. Control Optim. 1998, 36 (4), 1225-1239. (43) Teo, K. L.; Goh, C. J.; Wong, K. H. A Unified Computational Approach to Optimal Control Problems; Wiley: New York, 1991. (44) Viel, F.; Busvelle, E.; Gauthier, J. P. Stabilty of Polymerization Reactors using I/O Linearization and a High-Gain Observer. Automatica 1995, 31, 971-984. (45) Villadsen, J.; Michelson, M. L. Solution of Differential Equations by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, NJ, 1978.
Received for review June 30, 2000 Accepted January 24, 2001 IE0006312
