SECOND ORDER SUFFICIENT CONDITIONS FOR ...

SIAM J. CONTROL OPTIM. Vol. 41, No. 2, pp. 380–403

c 2002 Society for Industrial and Applied Mathematics

SECOND ORDER SUFFICIENT CONDITIONS FOR OPTIMAL CONTROL PROBLEMS WITH FREE FINAL TIME: THE RICCATI APPROACH∗ HELMUT MAURER† AND HANS JOACHIM OBERLE‡ Abstract. Second order sufficient conditions (SSC) for control problems with control-state constraints and free final time are presented. Instead of deriving such SSC from first principles, we transform the control problem with free final time into an augmented control problem with fixed final time for which well-known SSC exist. SSC are then expressed as a condition on the positive definiteness of the second variation. A convenient numerical tool for verifying this condition is based on the Riccati approach, where one has to find a bounded solution of an associated Riccati equation satisfying specific boundary conditions. The augmented Riccati equations for the augmented control problem are derived, and their modifications on the boundary of the control-state constraint are discussed. Two numerical examples, (1) the classical Earth-Mars orbit transfer in minimal time and (2) the Rayleigh problem in electrical engineering, demonstrate that the Riccati equation approach provides a viable numerical test of SSC. Key words. optimal control, control-state constraints, free final time, second order sufficient conditions, Riccati equation, Earth-Mars orbit transfer, Rayleigh problem AMS subject classifications. 49K15, 49K40, 65L10, 70M20, 94C99 PII. S0363012900377419

1. Introduction. In the last three decades, one can find an extensive literature on second order sufficient conditions (SSC) for optimal control problems with control and state constraints; cf. [2, 4, 6, 7, 8, 9, 10, 20, 25, 26, 31, 32, 35] and further literature cited in these papers. SSC have shown to be of fundamental importance for stability and sensitivity analysis of parametric optimal control problems; cf., e.g., [2, 5, 9, 10, 18, 19, 21, 23, 24, 33, 34]. SSC are usually expressed in terms of the positiveness of a quadratic form on a certain critical cone which is obtained through linearization of equality and inequality constraints. In general, such conditions are far too abstract to lend themselves to numerical verification. A practical test for SSC can be devised on the basis of a matrixvalued Riccati equation [23, 25, 37]. The main ideas underlying this approach are already exposed in the book of Bryson and Ho [1] for unconstrained control problems. This test requires the construction of a bounded solution to a Riccati equation which has to satisfy additional boundary conditions. An inherent difficulty arises from the fact that the coefficients of the Riccati equation depend on the accurate solution for state, control, and adjoint variables. Most of the cited papers deal with control problems on a fixed time interval. Extensions of the results to problems with free final time or with nonfixed time intervals have been discussed in [1, 4, 12, 26, 31, 32]. The method in Bryson and Ho [1, Chapter 6] uses heuristic arguments and also suffers from the drawback that, e.g., ∗ Received by the editors August 31, 2000; accepted for publication (in revised form) December 4, 2001; published electronically June 18, 2002. This work was supported by Deutsche Forschungsgemeinschaft under grant MA 891-3. http://www.siam.org/journals/sicon/41-2/37741.html † Westf¨ alische Wilhelms-Universit¨ at M¨ unster, Institut f¨ ur Numerische Mathematik, Einsteinstrasse 62, D-48149 M¨ unster, Germany ([email protected]). ‡ Universit¨ at Hamburg, Fachbereich Mathematik, Bundesstr. 55, D-20146 Hamburg, Germany ([email protected]).

380

SSC FOR CONTROL PROBLEMS WITH FREE FINAL TIME

381

time-optimal control problems are not tractable via this approach. More precisely, the function α defined in (6.6.13) of [1] is identically zero, and hence the quantity in (6.6.16) is not defined. A rigorous proof of the SSC in [1] may be found in Chamberland and Zeidan [4], where extensions of the results to control problems with mixed control-state constraints are also given. Again, however, the time-optimal case is not covered by these conditions. A remedy for this deficiency has been proposed in Hull [12] for unconstrained control problems. These conditions have been tested in a numerically unchallenging situation. We emphasize that a general approach for SSC on nonfixed time intervals has been developed by Osmolovskii [26, 31, 32], but the author does not offer any practical device to test his conditions numerically. The aim of the present paper is to develop verifiable SCC for control problems with free final time and mixed control-state constraints. Our approach is rather straightforward in the sense that it uses the well-known idea (cf., e.g., [11]) of reducing the free final time case to the fixed final time case by treating the free final time as an augmented state variable. It is not surprising that this procedure will lead to an augmented set of Riccati equations and boundary conditions. For unconstrained control problems, this derivation has already been described in [22]. The organization of the paper is as follows. In section 2, we recall known SSC for control problems with fixed final time [25, 37]. Section 3 describes the effect of the time transformation on the augmented variables and functions of the problem. A straightforward calculation shows that the Riccati equation for the augmented problem splits into three separate parts. A salient feature of the approach is that it suffices to solve a reduced form of the Riccati equation on the boundary of the controlstate constraint. The boundary conditions for the Riccati equation are worked out in some cases of practical interest. In particular, we derive additional sign conditions of the Riccati solution at the initial and final time which turn out to be crucial in the numerical test. In sections 4 and 5, we apply the numerical methodology to two practical and challenging examples. A highly accurate numerical solution to both examples is obtained via the multiple shooting method [3, 29]. The classical problem of a planar Earth-Mars orbit transfer in minimal time [14, 16, 17] is treated in section 4. The augmented Riccati equation test succeeds in confirming the optimality of the numerical solution. Section 5 presents a modification of the Rayleigh problem, which has been solved in [13, 23, 36] on fixed time intervals. Surprisingly, when no control constraints are imposed, the free final time problem has several local minima and one local maximum. The augmented Riccati test is capable of proving optimality for both local minima. Then the Rayleigh problem, subject to control constraints, is studied. We derive the reduced Riccati equations on the boundary of the constraint and compute a bounded solution which satisfies the extra boundary condition. Let us mention two further applications and extensions. First, on the basis of SSC, it is rather straightforward to perform a computational sensitivity analysis for both examples. The numerical techniques in [21, 24, 23, 33, 34] indicate that sensitivity differentials of optimal solutions with respect to parameters can be obtained through the solution of an additional linear boundary value problem (BVP). The second extension concerns optimal control problems with pure state constraints to which the techniques of this paper apply as well. 2. Second order conditions for control problems with fixed final time. We consider the following autonomous control problem (CP) subject to mixed control-state constraints: for a given final time T > 0, determine a control

382

HELMUT MAURER AND HANS JOACHIM OBERLE

function u ∈ L∞ (0, T ; Rm ) and a state function x ∈ W 1,∞ (0, T ; Rn ) that minimize the functional T F (x, u) = g(x(0), x(T )) + (2.1) L(x(t), u(t))dt 0

subject to (2.2)

x(t) ˙ = f (x(t), u(t))

for a.e. t ∈ [0, T ] ,

(2.3) (2.4)

ϕ(x(0), x(T )) = 0 , C(x(t), u(t)) ≤ 0 for a.e. t ∈ [0, T ] .

It is assumed that the functions g : Rn × Rn → R, L : Rn × Rm → R, f : Rn × Rm → Rn , ϕ : Rn × Rn → Rr , 0 ≤ r ≤ 2n, and C : Rn × Rm → Rk are C 2 functions on appropriate open sets. In this section, the final time T is supposed to be specified. Further, we assume that there exists a feasible pair of functions (x0 , u0 ) ∈ W 1,∞ (0, T ; Rn ) × L∞ (0, T ; Rm ) satisfying the constraints (2.2)–(2.4). The first order necessary conditions for an optimal pair (x0 , u0 ) are well known in the literature [11, 28]. The unconstrained Hamiltonian function H 0 , respectively, the augmented Hamiltonian H, are defined as (2.5) H 0 (x, u, λ) = L(x, u) + λ∗ f (x, u) ,

H(x, u, λ, µ) = H 0 (x, u, λ) + µ∗ C(x, u),

where λ ∈ Rn denotes the adjoint variable and µ ∈ Rk is the multiplier associated with the control-state constraint (2.4); the asterisk denotes the transpose. Henceforth, partial derivatives will often be denoted by subscripts. In the following, we shall make the hypothesis that first order conditions are satisfied in normal form with a nonzero cost multiplier. Hence we assume that there exist Lagrange multipliers (for convenience, we shall drop the lower subscript zero) (λ, µ, ρ) ∈ W 1,∞ (0, T ; Rn ) × L∞ (0, T ; Rk ) × Rr such that the following first order necessary conditions hold for a.e. t ∈ [0, T ]: (2.6)

˙ λ(t) = −Hx (x0 (t), u0 (t), λ(t), µ(t))∗ ,

(2.7)

(−λ(0), λ(T )) = ∇(x(0),x(T )) (g + ρ∗ ϕ)(x0 (0), x0 (T )) ,

(2.8)

Hu (x0 (t), u0 (t), λ(t), µ(t)) = 0,

(2.9)

µ(t) ≥ 0

(2.10)

H 0 (x0 (t), u0 (t), λ(t)) ≡ const.

and µ(t)∗ C(x0 (t), u0 (t)) = 0,

We shall use the notation [t] to denote arguments of functions at the reference solution x0 (t), u0 (t), λ(t), µ(t). To introduce regularity assumptions, we consider for β ≥ 0 the set of β-active constraints Iβ (t) := {i ∈ {1, . . . , k} | C i [t] ≥ −β} , where C i denotes the ith component of the vector C . In particular, for β = 0, we obtain the set of active indices I0 (t) = {i ∈ {1, . . . , k} | C i [t] = 0} .

383


The following regularity assumption concerns the linear independence of gradients for active constraints; cf. [18, 19, 21, 25, 37]. (A1) For some β > 0, the gradients Cui [t] are uniformly linear independent for all i ∈ Iβ (t) a.e. on [0, T ]. Further, we consider a margin δ ≥ 0 and define the set of indices Jδ (t) := {i ∈ {1, . . . , k} | µi (t) > δ},

jδ (t) := card(Jδ (t)),

where µi denotes the ith component of the multiplier µ. It is obvious that Jδ (t) ⊂ I0 (t) holds for all δ ≥ 0 . In particular, the strict complementarity condition µi (t) > 0 is valid for all indices i ∈ Jδ (t) . It will be convenient to introduce the notation C δ [t] = (C i [t]) i∈Jδ (t) . We assume that the following modified strict Legendre–Clebsch condition [18, 19, 21, 25, 37] is satisfied, where | . | denotes the euclidean norm. (A2) For some δ > 0, there exists c > 0 such that, for all t ∈ [0, T ], the estimate v ∗ Huu [t]v ≥ c|v|2 holds for all v ∈ Rm satisfying Cuδ [t]v = 0. SSC can now be derived by studying the behavior of the second variation on the variational system associated with (2.2)–(2.4). In what follows, we shall use the abbreviation ϕ[ 0, T ] = ϕ(x0 (0), x0 (T )) and similar notation. The variational system of equations (2.2)–(2.4) is the set of functions (y, v) ∈ W 1,2 (0, T ; Rn ) × L2 (0, T ; Rm ) satisfying y(t) ˙ = fx [t]y(t) + fu [t]v(t),

(2.11) (2.12)

a.e. t ∈ [0, T ] , Dx(0) ϕ[0, T ]y(0) + Dx(T ) ϕ[0, T ]y(T ) = 0 ,

(2.13)

Cxi [t]y(t) + Cui [t]v(t) = 0 ∀ i ∈ Jδ (t) , a.e. t ∈ [0, T ] .

Moreover, we introduce the function G(x(0), x(T )) := g(x(0), x(T )) + ρ∗ ϕ(x(0), x(T )) and define the (2n, 2n)-matrix (2.14)

Γ[0, T ] :=

2 D(x(0),x(T ))

G(x0 (0), x0 (T )) =

G00 [0, T ] G0T [0, T ] GT 0 [0, T ] GT T [0, T ]

2 2 with obvious notation G00 [0, T ] = D(x(0),x(0)) G[0, T ], G0T [0, T ] = D((x(0),x(T )) G[0, T ], etc. Then the so-called second variation is given by the quadratic form 1 T y(t) Hxx [t] Hxu [t] J 2 (y, v) = dt (y(t)∗ , v(t)∗ ) Hux [t] Huu [t] v(t) 2 0 1 y(0) + (y(0)∗ , y(T )∗ )Γ[0, T ] (2.15) . y(T ) 2

The next theorem summarizes the SSC for a weak local minimum which are to be found in [21, 25, 35, 37]. Theorem 2.1 (SSC for control problems with fixed final time). Let (x0 , u0 ) be admissible for problem (CP). Suppose that there exist multipliers (λ, µ, ρ) ∈ W 1,∞ (0, T ; Rn ) × L∞ (0, T ; Rk ) × Rr such that the following conditions hold:

384


(1) the necessary conditions (2.6)–(2.10) are satisfied; (2) assumptions (A1) and (A2) hold; (3) there exist γ0 > 0 such that the quadratic form in (2.15) can be estimated from below as J 2 (y, v) ≥ γ0 ( ||y||21,2 + ||v||22 )

(2.16)

for all variations (y, v) ∈ W 1,2 (0, T ; Rn ) × L2 (0, T ; Rm ) satisfying the variational system (2.11)–(2.13); (4) if u0 is continuous, then one may choose β = 0 and δ = 0 in assumptions (A1) and (A2) and in condition (3). Then for all constants 0 < γ < γ0 with γ0 as in (2.16) there exists α > 0 such that F (x, u) ≥ F (x0 , u0 ) + γ ( ||x − x0 ||21,2 + ||u − u0 ||22 ) holds for all admissible (x, u) with ||x − x0 ||1,∞ + ||u − u0 ||∞ ≤ α . In particular, (x0 , u0 ) provides a strict weak local minimum for problem (CP). The SSC in the previous theorem usually are not suitable for a direct numerical verification in practical control problems. Let us mention that, for a discretized version of the control problem (CP), optimization techniques have been developed that allow us to check the positiveness condition by computing the reduced Hessian; cf. [2]. In order to obtain verifiable sufficient conditions for the control problem in function spaces, the SSC in Theorem 2.1 are strengthened in the following way. Consider a symmetric matrix function Q ∈ W 1,∞ (0, T ; Mn×n ) . For every vari˙ − ation y(t) satisfying the linearized state equation (2.11), we have y(t)∗ Q(t)(y(t) fx [t]y(t) − fu [t]v(t)) ≡ 0. Adding the last identity to the second variation J 2 (y, v) in (2.15) and performing a partial integration, we find that the definiteness condition (2.16) in Theorem 2.1 holds if the following two conditions (a) and (b) are satisfied. Condition (a). There exist a symmetric matrix Q ∈ W 1,∞ (0, T ; Mn×n ) and γ > 0 such that ˙ y Q(t) + Q(t)fx [t] + fx [t]∗ Q(t) + Hxx [t] Hxu [t] + Q(t)fu [t] ∗ ∗ (y , v ) v Huu [t] Hux [t] + fu [t]∗ Q(t) (2.17)

≥ γ |(y, v)|2

holds uniformly in t ∈ [0, T ] for all vectors (y, v) ∈ Rn × Rm with (2.18)

Cxi [t]y + Cui [t]v = 0 ∀ i ∈ Jδ (t) .

Condition (b). The boundary condition G00 [0, T ] + Q(0) G0T [0, T ] ξ0 (ξ0∗ , ξ1∗ ) (2.19) >0 GT 0 [0, T ] GT T [0, T ] − Q(T ) ξ1 is valid for all (ξ0 , ξ1 ) ∈ Rn × Rn \ {0} satisfying (2.20)

Dx(0) ϕ[0, T ]ξ0 + Dx(T ) ϕ[0, T ]ξ1 = 0 .

A first consequence is that the definiteness condition (2.16) holds if the matrix in (2.17) is positive definite on the whole space Rn × Rm and if conditions (2.19)


385

and (2.20) are satisfied. First, this leads to the requirement that the strict Legendre– Clebsch condition Huu [t] ≥ c · Im

∀ t ∈ [0, T ] ,

c > 0,

is valid on the whole interval [0, T ]. Second, by evaluating the Schur complement of this matrix and using the continuous dependence of ODEs on system data, the estimate (2.16) follows from the following assumption: there exists a solution of the Riccati equation (2.21) Q˙ = −Qfx [t]−fx [t]∗ Q−Hxx [t]+(Hxu [t]+Qfu [t])Huu [t]−1 (Hxu [t]+Qfu [t])∗ , for a.e. t ∈ [0, T ] such that the matrix function Q(t) is bounded on [0, T ] and satisfies the boundary conditions (2.19) and (2.20); cf. [25, Theorem 5.2]. However, in some applications, these conditions are too strong since the Riccati equation (2.21) may fail to have a bounded solution; cf. the Rayleigh problem in [23]. A weaker condition can be obtained by introducing the following modified or reduced Riccati equation. For δ ≥ 0, recall the definition of the vector C δ [t] = (C i [t])i∈Jδ (t) ,

jδ (t) = ca rd(Jδ (t)) .

Then the matrix Cuδ [t] of partial derivatives has dimension jδ (t) × m. For simplicity, the time argument will be omitted in what follows. The pseudoinverse of the matrix Cuδ is given by the (m × jδ )-matrix (Cuδ )+ := (Cuδ )∗ (Cuδ (Cuδ )∗ )−1 , which exists in view of the linear independence assumption (A1). Furthermore, let (Cuδ )⊥ denote an (m × (m − jδ ))-matrix whose column vectors form an orthogonal basis of the kernel Ker(Cuδ ). Consider then the following matrices (cf. [9, 10, 25, 37]): (2.22) (2.23) (2.24) (2.25)

Dδ := −(Cuδ )+ Cxδ , P δ := (Cuδ )⊥ , Aδ := fx + fu Dδ , δ := Hxx + Hxu Dδ + (Dδ )∗ Hux + (Dδ )∗ Huu Dδ , Hxx δ Hxu := Hxu + (Dδ )∗ Huu , δ (−1) (Huu ) := P δ ( (P δ )∗ Huu P δ )−1 (P δ )∗ .

δ (−1) ) in (2.25) is well defined by virtue of asNote that the (m × m)-matrix (Huu sumption (A2). It follows that the estimate (2.16) holds if there exists a symmetric matrix function Q(t) that solves the Riccati equation

(2.26)

δ δ δ (−1) δ Q˙ = −QAδ − (Aδ )∗ Q − Hxx + (Hxu + Qfu )(Huu ) (Hxu + Qfu )∗

for a.e. t ∈ [0, T ] such that Q(t) is bounded on [0, T ] and satisfies the boundary conditions (2.19), (2.20). In general, it is rather tedious to elaborate this Riccati equation explicitly. To facilitate the numerical treatment in practical applications, we discuss special cases in more detail. On interior arcs with C[t] < 0, we have jδ (t) = 0, and thus the Riccati equation (2.26) reduces to the one introduced in (2.21). Consider now a boundary arc with jδ (t) = m, where we have as many control components as active constraints. Due to assumption (A1), the pseudoinverse is given by (Cuδ )+ = (Cuδ )−1 , and hence

386


the matrix P δ = (Cuδ )⊥ = 0 in (2.22) vanishes. Then the matrices in (2.22)–(2.25) become (2.27) (2.28)

δ (−1) Dδ = −(Cuδ )−1 Cxδ , Aδ = fx − fu (Cuδ )−1 Cxδ , (Huu ) = 0, δ δ −1 δ δ −1 δ ∗ δ −1 δ Hxx = Hxx − Hxu (Cu ) Cx + [(Cu ) Cx ] [Huu (Cu ) Cx − Hux ] ,

and thus the Riccati equation (2.26) reduces to the linear ODE (2.29)

δ Q˙ = −QAδ − (Aδ )∗ Q − Hxx .

Another special case arises for pure control constraints where we have Cx ≡ 0 . Then (2.28) and (2.29) simplify to the linear ODE (2.30)

Q˙ = −Qfx − fx∗ Q − Hxx .

The Rayleigh problem in [23] provided an illustrative application of this approach. Let us also evaluate the boundary conditions (2.19) and (2.20) in a special case of practical interest. Suppose that the boundary conditions are separated and that some components for the initial and final state are fixed according to (2.31) xk (0) = ak for k ∈ K0 ⊂ {1, . . . , n},

xk (T ) = bk for k ∈ KT ⊂ {1, . . . , n},

whereas the other components are free. Denote the complements of the index sets by K0c = {1, . . . , n} \ K0 , KTc = {1, . . . , n} \ KT . Then it is easy to see that the boundary conditions (2.19), (2.20) are satisfied if the following submatrices are positive definite: (2.32)

[ Q(0) ](i,j)∈K0c ×K0c > 0 ,

[ −Q(T ) ](i,j)∈KTc ×KTc > 0 .

By virtue of the continuous dependence of solutions to ODEs on systems data, one of these definiteness conditions can be relaxed. For example, it suffices to require only positive semidefiniteness [ Q(0) ](i,j)∈K0c ×K0c ≥ 0 . This relaxation will be convenient for the numerical verification of SSC applied to the examples in sections 4 and 5. 3. SSC for control problems with free final time. We consider again the control problem (CP) in (2.1)–(2.4), but in this section the final time will not be specified. It will always be assumed that the final time T > 0 is positive. The first order necessary conditions for problem (CP) with free final time are well known [11, 28] and extend those given in the last section. Let H be the Hamiltonian defined in (2.5). Then it is assumed that there exist multipliers in normal form, (λ, µ, ρ) ∈ W 1,∞ (0, T ; Rn ) × L∞ (0, T ; Rk ) × Rr , which satisfy (2.6)–(2.10). In addition, the following transversality condition associated with the free final time T holds: (3.1)

H[T ] = 0 ,

i.e., H[t] ≡ 0

∀ t ∈ [0, T ] .

These conditions can be obtained by transforming problem (CP) with free final time with fixed final time T˜ = 1. The transformation proceeds by T into a problem (CP)


387

augmenting the state dimension and by introducing the free final time as an additional state variable. Indeed, it is this transformation that will allow us to develop SSC for the free final time case on the basis of the SSC in Theorem 2.1 for fixed final time. Now define the new time variable τ ∈ [0, 1] by (3.2)

t=τ ·T ,

0≤τ ≤1.

We shall use the same notation x(τ ) := x(τ ·T ) and u(τ ) := u(τ ·T ) for the state and the control variable with respect to the new time variable τ . The augmented state x x ˜ := ∈ Rn+1 , xn+1 := T , xn+1 satisfies the differential equations (3.3)

dx/dτ = T · f (x(τ ), u(τ )) ,

dxn+1 /dτ ≡ 0 .

As a result of this time transformation, we consider the following augmented control on the fixed time interval [0, 1]: minimize the functional problem (CP) 1 ˜ x(τ ), u(τ )) dτ F (˜ x, u) = F (x, T, u) = g˜(˜ x(0), x ˜(1)) + (3.4) L(˜ 0

subject to (3.5)

d˜ x/dτ = f˜(˜ x(τ ), u(τ )),

(3.6)

ϕ(˜ ˜ x(0), x ˜(1)) = 0 , ˜ C(˜ x(τ ), u(τ )) ≤ 0,

(3.7)

a.e. τ ∈ [0, 1] ,

a.e. τ ∈ [0, 1] .

The functions herein are given by (3.8) (3.9) (3.10)

g˜(˜ x(0), x ˜(1)) := g(x(0), x(1)) , T · f (x, u) f˜(˜ x, u) := , 0

˜ x, u) := T · L(x, u) , L(˜

ϕ(˜ ˜ x(0), x ˜(1)) := ϕ(x(0), x(1)) ,

˜ x, u) := C(x, u) . C(˜

on the fixed time interval [0, 1] falls into the category The transformed problem (CP) of control problems treated in the preceding section. Thus we are able to obtain by evaluating the SSC in Theorem 2.1 for the SSC for the transformed problem (CP) augmented state variable x ˜ = (x, T ) . to those of problem (CP) on the First we relate the multipliers for problem (CP) becomes time interval [0, T ]. The Hamiltonian for problem (CP) (3.11)

˜ µ ˜ ∗ f˜(˜ ˜ x, u, λ, ˜ x, u) + λ ˜ x, u) H(˜ ˜) = L(˜ x, u) + µ ˜∗ C(˜ ˜∗ C(x, u), = T · [ L(x, u) + λ∗ f (x, u) ] + µ

˜ ∗ = (λ∗ , λn+1 ) ∈ Rn+1 denotes the augmented adjoint variable and µ where λ ˜ ∈ Rk is the multiplier for the constraint (3.7). Introducing the scaled multiplier µ := µ ˜/T , ˜ is related to the Hamiltonian H in (2.5) as follows: the Hamiltonian H ˜ µ ˜ x, u, λ, (3.12) H(˜ ˜) = T · [ L(x, u) + λ∗ f (x, u) + µ∗ C(x, u) ] = T · H(x, u, λ, µ) .

388


According to (2.6)–(2.10), the multipliers for problem (CP), ˜ µ (λ, ˜, ρ) ∈ W 1,∞ (0, 1; Rn+1 ) × L∞ (0, 1; Rk ) × Rr ,

˜= λ

λ λn+1

,

satisfy the following necessary conditions for a.e. τ ∈ [0, 1]: ˜ T [τ ] = −H[τ ] , dλn+1 /dτ = −H

(3.13)

˜ x [τ ]∗ = −T · Hx [τ ]∗ , dλ/dτ = −H

(3.14)

(−λ(0), λ(1)) = ∇(x(0),x(1)) (g + ρ∗ ϕ)(x0 (0), x0 (1)) ,

(3.15)

λn+1 (0) = λn+1 (1) = 0 ,

(3.16)

˜ u [τ ] = 0 , H

(3.17)

µ ˜(τ ) ≥ 0

(3.18)

˜ ] ≡ const. H[τ

and µ ˜(τ )∗ C[τ ] = 0, ∀ τ ∈ [0, 1].

Relations (3.13), (3.15), and (3.18) immediately yield ˜ = T · H[1] , 0 = H[1]

(3.19)

which proves the transversality condition (3.1). To check the Legendre–Clebsch condition in assumption (A2), one has to observe the scaling ˜ uu [τ ] = T · Huu [τ ] . H

(3.20)

we have to evaluate In order to apply the SSC in Theorem 2.1 to problem (CP), all terms in relations (2.11)–(2.32) for the tilde quantities defined in (3.8)–(3.10). Recalling the scaled multiplier µ0 = µ ˜0 /T in (3.12), we obtain the transformed quantities T · fx f T · fu , f˜u = , C˜x˜ = (Cx , 0) , C˜u = Cu , (3.21) f˜x˜ = 0 0 0 (3.22)

˜ x˜x˜ = H

T · Hxx Hx0

(Hx0 )∗ 0

,

˜ xu = H

T · Hxu Hu0

,

˜ uu = T · Huu . H

˜ T = H 0 , which can be Observe that the last relations make use of the identity H 0 seen from (3.11) with H denoting the unconstrained Hamiltonian. The preceding with free final time. considerations lead to the following SSC for problem (CP) Theorem 3.1 (SSC for control problems with free final time). Let (x0 , T0 , u0 ) Suppose that there exist multipliers with T0 > 0 be admissible for problem (CP). 1,∞ n ∞ k (0, 1; R ) × L (0, 1; R ) × Rr such that the following conditions hold: (λ, µ, ρ) ∈ W (1) the necessary conditions (3.13)–(3.19) are satisfied; (2) assumptions (A1) and (A2) hold with respect to the time interval [0, 1]; (3) there exists γ0 > 0 such that the quadratic form (2.15) expressed in terms of the transformed functions (3.21), (3.22) can be estimated from below as J 2 (˜ y , v˜) ≥ γ0 ( ||˜ y ||21,2 + ||˜ v ||22 ) for all variations (˜ y , v˜) ∈ W 1,2 (0, 1; Rn+1 ) × L2 (0, 1; Rm ) which satisfy the variational system (2.11)–(2.13) on the time interval [0, 1].


389

(4) if u0 is continuous, then one may choose β = 0 and δ = 0 in assumptions (A1) and (A2) and in condition (3). Then for all constants 0 < γ < γ0 there exists α > 0 such that F (x, T, u) ≥ F (x0 , T0 , u0 ) + γ ( ||x − x0 ||21,2 + |T − T0 |2 + ||u − u0 ||22 ) holds for all admissible (x, T, u) with ||x − x0 ||1,∞ + |T − T0 | + ||u − u0 ||∞ ≤ α. In particular, (x0 , T0 , u0 ) provides a strict weak local minimum for problem (CP). Note that this theorem immediately yields SSC for the original problem (CP) since we have identified the pair of state and control functions (x(t), u(t)) = (x(τ ·T ), u(τ ·T )) on the interval [0, T ] with the pair (x(τ ), u(τ )) on the interval [0, 1]. It is apparent that these conditions are not very handy in practical applications. Again, we may resort to the Riccati equations and boundary conditions developed in (2.21)–(2.32). Now we consider the augmented (n + 1, n + 1)-matrix Q R ˜ (3.23) , Q= R ∗ qT where Q is a symmetric (n×n)-matrix, R is an n-vector, and qT is a scalar. Inserting the transformed quantities (3.21) and (3.22) into the Riccati equation (2.21) for the ˜ we obtain a Riccati equation for Q, a linear equation for R, and a direct matrix Q, integration for qT on the interval [0, 1]; the argument τ is omitted for simplicity: (3.24) dQ/dτ = T · [ −Qfx − fx∗ Q − Hxx + (Hxu + Qfu )(Huu )−1 (Hxu + Qfu )∗ ] , (3.25) dR/dτ = −Qf − T fx∗ R − (Hx0 )∗ + (Hxu + Qfu )(Huu )−1 (Hu0 + T fu∗ R)∗ , (3.26) dqT /dτ = −2R∗ f +

1 · (Hu0 + T fu∗ R)(Huu )−1 (Hu0 + T fu∗ R)∗ . T

Clearly, the Riccati equation (3.24) evaluated on the interval [0, 1] agrees with the Riccati equation (2.21) on the interval [0, T ]. Note again that Hu0 ≡ 0 holds on totally interior arcs with C[t] < 0. We wish to draw attention to the fact that (3.25) and (3.26) are not identical to corresponding equations in Bryson and Ho [1, sections 6.6, 6.7] and Chamberland and Zeidan [4, formulae (28)–(30)], or Hull [12]. The modified Riccati equation (2.26) can be worked out on the time interval [0, 1] in a similar way using the transformed quantities (3.21) and (3.22). However, since this procedure is quite cumbersome, we restrict the discussion to the special case m = jδ (t), which was considered already in (2.27)–(2.29). Upon computing the matrices in (2.27) and (2.28), Aδ = fx − fu (Cuδ )−1 Cxδ , δ Hxx = Hxx − Hxu (Cuδ )−1 Cxδ + [(Cuδ )−1 Cxδ ]∗ [Huu (Cuδ )−1 Cxδ − Hux ] , for the tilde quantities (3.21), (3.22), we recognize that the Riccati equation (2.29) splits into the following three equations: (3.27)

δ ], dQ/dτ = T · [ −QAδ − (Aδ )∗ Q − Hxx

(3.28)

dR/dτ = −Qf − T · (Aδ )∗ R − (Hx0 )∗ − [Hu0 (Cuδ )−1 Cxδ ] ∗ ,

(3.29)

dqT /dτ = −2R∗ f .

390


These formulas simplify considerably if Cx ≡ 0 holds, i.e., if the constraint is a pure control constraint. Then we get Aδ = fx , and the last equations yield (3.30)

dQ/dτ = T · [ −Qfx − fx∗ Q − Hxx ] , dR/dτ = −Qf − T · fx∗ R − (Hx0 )∗ , dqT /dτ = −2R∗ f .

These equations will provide a convenient test for SSC when applied to the Rayleigh problem in section 5. It is rather tedious to write out the boundary conditions (2.19) and (2.20) in the general case. We shall only discuss the important case in which the initial and final states are fixed; i.e., x(0) = x0 and x(1) = x1 hold with prescribed x0 , x1 ∈ Rn . In this situation, the positive definiteness condition (2.32) evaluated for the augmented ˜ reduces to the following boundary conditions: matrix Q qT (0) > 0

(3.31)

and qT (1) < 0 .

These conditions constitute extra conditions for the free final time case and will turn out to be crucial for the numerical examples discussed in the next two sections. Note that we may relax one of these conditions, e.g., the initial condition, to qT (0) ≥ 0 by virtue of the continuous dependence of ODE solutions on initial data. 4. Planar Earth-Mars transfer with minimal flight time. Rocket flights in an inverse square law field have been studied extensively in the literature; see, e.g., Kelley [14, 15], Kenneth and McGill [16], Lawden [17], Moyer and Pinkham [27], and Oberle and Taubert [30]. We consider the classical Earth-Mars orbit transfer with minimal transfer time. The state variables are r: distance of the vehicle to the sun; w: radial component of the velocity; v: horizontal component of the velocity; m: mass of the vehicle. The control variable is ψ: angle of the thrust vector with respect to local horizon. The thrust is always at its maximal value βmax since the final time is minimized. All variables are scaled according to the dynamic model treated in Oberle and Taubert [30]. The optimal control problem is to minimize the final time T (4.1) 1 dt F (x, T, u) = T = 0

subject to the equations of motion,

(4.2)

dr/dt

= w,

dw/dt

=

dv/dt dm/dt

1 v2 c − 2 + βmax sin ψ, r r m wv c = − + βmax cos ψ, r m = −βmax ,

and the boundary conditions for the initial and final state (4.3)

r(0) = 1.0,

w(0) = 0.0,

v(0) = 1.0,

r(T ) = 1.525,

w(T ) = 0.0,

v(T ) = 1.0/

m(0) = 1.0, r(T ).

The constants are given by c = 1.872 and βmax = 0.075. The underlying physical data of the vehicle are given as follows: the initial mass is m0 = 679.78 kg; the (maximal) thrust = 0.56493 N; and the (constant) equivalent exit velocity is vc = 55809 m/s.


391

4.1. The BVP. Recall now the time transformation t = τ · T introduced in (3.2). Since there is no control constraint in this problem, the Hamiltonian (3.12) is given by 2 1 v c ˜ = T · 1 + λr w + λ w (4.4) − 2 + βmax sin ψ H r r m wv c + λv − + βmax cos ψ − λm βmax . r m Let us evaluate the first order optimality conditions (3.13)–(3.16) and omit for convenience the subscript zero referring to the optimal solution. The optimal control ψ is derived from condition (3.16) and the assumed Legendre–Clebsch condition (A2) as (4.5)

sin ψ = −

λw , λ2w + λ2v

cos ψ = −

λv . λ2w + λ2v

The adjoint equations (3.13) on the normalized time interval [0, 1] are given by

2 v wv 2 − 3 − λv 2 , dλr /dτ = T · λw r2 r r

v , dλw /dτ = T · − λr + λv r

(4.6) 2v w dλv /dτ = T · − λw + λv , r r

βmax c 2 2 . dλm /dτ = T · − λ + λ w v m2 The transversality conditions (3.14) and (3.19) yield (4.7)

λm (1)

=

0,

H[1]

=

1 −

βmax c λw (1)2 + λv (1)2 = 0 . m(1)

After transforming the state equations (4.2) to the normalized time interval [0, 1] according to (3.3), we obtain a two-point BVP consisting of (4.2)–(4.8). This BVP can be further simplified by eliminating the variables m(τ ) and λm (τ ). The variable m(τ ) is substituted according to (4.8)

m(τ ) = 1.0 − βmax T · τ,

and the variable λm can be dropped since it does not enter into the first three equations in (4.6). The reduced BVP then comprises the six ODEs with respect to the variables r, w, v, λr , λw , and λv on the fixed time interval [0, 1] and the trivial equation dT /dτ ≡ 0 in view of (3.3). The corresponding boundary conditions are the six boundary conditions (4.3) with respect to r, w, and v and the Hamiltonian boundary condition in (4.7). Once a solution of this BVP has been determined, the adjoint variable λm can be obtained through an integration of the last equation in (4.6): 1 βmax c λm (τ ) = T · λw (τ )2 + λv (τ )2 dτ . 2 τ m(τ ) Alternatively, λm can be eliminated from the condition H ≡ 0.

392

HELMUT MAURER AND HANS JOACHIM OBERLE 1.6

0.35 0.3

1.5

0.25

radial velocity w

distance r

1.4

1.3

1.2

0.2 0.15 0.1

1.1

0.05

1 0

0.2

0.4

τ = t/T

0.6

0.8

0 0

1

1.15

0.2

0.4

τ = t/T

0.6

0.8

1

0.8

1

6

1.1

5

4

1

control ψ

tangent. velocity v

1.05

0.95

3

0.9

2

0.85 1

0.8 0.75 0

0.2

0.4

τ = t/T

0.6

0.8

1

0 0

0.2

0.4

τ = t/T

0.6

Fig. 4.1. Earth-Mars transfer: State variables r, w, v and control variable ψ.

The code BNDSCO in Oberle and Grimm [29] provides the following initial and final values for the adjoint variables and the final time:

(4.9)

λr (0)

= −0.52729 67236 × 101 ,

λr (1)

= −0.37511 95452 × 101 ,

λw (0)

= −0.26088 76037 × 101 ,

λw (1)

=

λv (0)

= −0.56884 53434 × 101 ,

λv (1)

= −0.35509 23985 × 101 ,

T

=

0.40004 02548 × 101 ,

0.33199 21219 × 101 .

Figure 4.1 displays the corresponding state variables r, w, v and the control variable ψ, while Figure 4.2 shows the adjoint variables. 4.2. SSC. Let us first check the strict Legendre–Clebsch condition in assumption ˜ ψ ψ [τ ] = T · Hψ ψ [τ ] in (3.20). We obtain (A2), taking into account the scaling H Hψ ψ =

βmax c 2 λw + λ2v m

and find ˜ ψ ψ [τ0 ] = 0.07390 1871 > 0, ˜ ψ ψ [τ ] | τ ∈ [ 0, 1] } = H min {H

τ0 = 0.50935 25818 ,

which verifies the strict Legendre–Clebsch condition (A2). To prove the SSC in Theorem 3.1, it remains to show that the Riccati equations (3.24)–(3.26) possess a bounded solution such that the sign conditions (3.31) hold.

393

SSC FOR CONTROL PROBLEMS WITH FREE FINAL TIME -

2.5

-

5 4

3

3 3.5

-

-

2

λw

λr

-

4

1 0

4.5

-1 -

-

5

5.5

-2 0

0.2

0.4

τ = t/T

0.6

0.8

1

-3

1

0

0.2

0.4

0.2

0.4

τ = t/T

0.6

0.8

1

0.6

0.8

1

2 1.8

0

1.6 -1

1.4 1.2

λ

λv

m

-2 -3

1

0.8 0.6

-4

0.4 -5 -6

0.2 0

0.2

0.4

τ = t/T

0.6

0.8

1

0 0

τ = t/T

Fig. 4.2. Earth-Mars transfer: Adjoint variables λr , λw , λv , λm .

The reader is reminded that we have eliminated the state variable m so that the remaining state variables r, w, v have fixed final values. The symmetric Riccati matrix ˜ in (3.23) is given in the form Q   q11 q12 q13 r1  q12 q22 q23 r2  Q R ˜=  = Q  q13 q23 q33 r3  . R∗ q r1 r2 r3 q T We refrain from writing down the Riccati equations (3.24)–(3.26) explicitly. The evaluation is rather tedious but can be simplified with the help of symbolic computations offered, for example, in the package MAPLE. It should be noted that the coefficients of the Riccati equation are functions of the nominal trajectory characterized by (4.9). ˜ We merely indicate how to find appropriate initial values for Q(0) such that the sign conditions (3.31) hold: qT (0) > 0,

qT (1) < 0.

We succeeded using a rather heuristic optimization technique. Starting with initial estimates for Q(0) which allowed the integration of (3.24)–(3.26) on [0, 1], we changed iteratively one component of Q(0) in order to minimize qT (1). Changing the indices of components in a cyclic way, we got the following initial values: q11 (0) q22 (0) r1 (0) qT (0)

= = = =

1.0, q12 (0) = 2.0, q13 (0) −50.0, q23 (0) = −10.0, q33 (0) 80.0, r2 (0) = 40.0, r3 (0) 10.0 > 0.

= 1.0, = −100.0, = 100.0,

For these data, a solution of the Riccati equations was found to exist on the whole interval [0, 1] with final value qT (1) = −18.090 44002 < 0 . Hence all assumptions for

394


5

35 30

0

25 -5

12

q11

20

q

-10

15 -15 10 -20

-25

5

0

0.2

0.4

τ = t/T

0.6

0.8

0 0

1

5

0.2

0.4

τ = t/T

0.6

0.8

1

0.6

0.8

1

0.8

1

0.6

0.8

1

0.6

0.8

1

0 -10

0 -20 -5

q22

q13

-30 -40

-10

-50 -15 -60 -20

0

0.2

0.4

τ = t/T

0.6

0.8

-70

1

0

-10

20

-20

15

-30

10

-40

33

25

5 0

0.4

τ = t/T

-50 -60

-5

-70

-10

-80

-15 -20

0.2

0

q

q23

30

-90 0

0.2

0.4

τ = t/T

0.6

0.8

1

-100

80

0

0.2

0.4

τ = t/T

0.6

60

75

40

70 20

65 60

r2

r

1

0

55

-20

50 45

-40

40 -60

35 30 0

0.2

0.4

τ = t/T

0.6

0.8

1

-80

0

5

60

0

qT

10

80

r3

100

40

20

0

0.4

0.2

0.4

τ = t/T

-5

-10

0

-20

0.2

-15

0.2

0.4

τ = t/T

0.6

0.8

1

-20

0

τ = t/T

Fig. 4.3. Solutions of the Riccati equations (3.24)–(3.26).

Theorem 3.1 are verified, and we draw the conclusion that the trajectory characterized by (4.9) provides a weak local minimum for problem (4.1)–(4.3). The component ˜ ] are shown in Figure 4.3. functions of Q[τ 5. Control of current in a tunnel-diode oscillator: Rayleigh problem with control constraints. The following Rayleigh problem has been treated in [13, 36] as a fixed final time control problem. SSC and sensitivity analysis for this model have been discussed in Maurer and Augustin [23]. In this section, we investigate a slightly modified problem with free final time. Figure 5.1 displays an electric circuit


395

Fig. 5.1. Tunnel-diode oscillator, x1 (t) = I(t).

(tunnel-diode oscillator), where L denotes inductivity, C denotes capacity, R denotes resistance, I denotes electric current, and D is a diode. The state variables are the electric current x1 (t) = I(t) at time t ∈ [ 0, T ] and x2 (t) := x˙ 1 (t) . The control u(t) is a suitable transformation of the voltage V0 at the generator. With an additional parameter c ≥ 0, the Rayleigh problem with free final time is defined as follows: minimize the functional T T Fc (x, T, u) = c · T + (5.1) ( u(t)2 + x1 (t)2 ) dt = ( c + u(t)2 + x1 (t)2 ) dt 0

0

subject to (5.2)

x˙ 1 (t) = x2 (t) , x˙ 2 (t) = −x1 (t) + x2 (t) ( 1.4 − 0.14 x2 (t)2 ) + 4 u(t) ,

(5.3)

x1 (0) = x2 (0) = −5 ,

(5.4)

| u(t) | ≤ 1

x1 (T ) = x2 (T ) = 0 ,

for t ∈ [0, T ] .

The solution of this problem with final time T = 4.5 specified and c = 0 may be found in [23]. In the following, we denote by Fc (T ) the optimal value of the control problem (5.1)–(5.4) for fixed final time T . The behavior of this function gives insight into the behavior of optimal solutions for free final time. 5.1. Unconstrained optimal solutions. We consider the unconstrained problem with control constraint (5.4) deleted. After applying the time transformation (3.2), the unconstrained Hamiltonian in (3.12) becomes ˜ = T · [ c + u2 + x2 + λ1 x2 + λ2 (−x1 + x2 (1.4 − 0.14x2 ) + 4u) ]. ˜ 0 (˜ x, u, λ) (5.5) H 1 2 Henceforth, we omit the lower index zero to denote the optimal solution. The control ˜ u0 [τ ] = 0, which yields is computed from the equation H u(τ ) = −2λ2 (τ ) . The transformed state equations (3.3) and adjoint equations (3.13) lead to the following ODEs in [0, 1] :

(5.6)

dx1 /dτ dx2 /dτ dλ1 /dτ dλ2 /dτ

= = = =

T T T T

· x2 , · (−x1 + 1.4x2 − 0.14x32 − 8λ2 ), · (−2x1 + λ2 ), · (−λ1 − 1.4λ2 + 0.42x22 λ2 ).

396


Fc (T )

F c (T )

30.2

30:0

.. ... .. .. ... ... .. ... ... ... .. ... ... ... ... ... .... ... .... .... .... ..... ..... ..... ..... ..... ...... ......... ........... .......................................... ..................................................................................................... ........................ ............... ..

... .. ........ .. ....... .. ....... .. ....... .. ....... . .. . . . . . .. ..... ... ..................... ...... ............ ...... ....... ....... ..... ...... .... ....... .... ....... .... ...... . . .... . . . . .... ..... ..... ....... ...... ........ ...... ....... ............ ........................ .........

30.1

30.0

T

29.9

2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0

29:9

T

2.0

2.1

2.2

2.3

2.4

2.5

Fig. 5.2. Graph of Fc (T ) for c = 1/16 and c = 0.

To compute the optimal solution of the control problem with fixed final time, we resolve the BVP which comprises (5.6) and the boundary conditions (5.7)

x1 (0) = x2 (0) = −5,

x1 (1) = x2 (1) = 0 .

The BVP is solvable only for final times T > T ∗ > 0 , where T ∗ is a suitable final time. The optimal value function Fc (T ) is shown in Figure 5.2 for T > T ∗ = 2 and two different values c = 1/16 and c = 0. Observe that, for c = 1/16, the graph of Fc (T ) in Figure 5.2 depicts two local minima and one local maximum with respect to the final time T , whereas no clear minimum can be discerned for c = 0 . For free final time T , the transversality condition (3.19) and the boundary con˜ 0 [1] = T · (c + u(1)2 + 4λ2 (1)u(1)) = T · (c − 4λ2 (1)2 ), from ditions (5.7) yield 0 = H which we get the boundary condition (5.8)

λ2 (1)2 = c/4 .

Now we solve the BVP (5.6)–(5.8) for the cases c = 1/16 and c = 0, using again the code BNDSCO in [29]. Case c = 1/16. We find three solutions: Solution 1:

T = 2.19460 79912, λ1 (0) = −9.01234 54748, λ2 (0) = −2.67606 29500,

Fc (T ) = 30.06097 62322, λ1 (1) = 0.97693 36044, λ2 (1) = 0.125.

Solution 2:

T = 2.46029 38602, λ1 (0) = −9.01228 20002, λ2 (0) = −2.67605 43511,

Fc (T ) = 30.07173 02593, λ1 (1) = 0.95904 74639, λ2 (1) = −0.125.

Solution 3:

T = 3.51535 36980, λ1 (0) = −9.01085 93855, λ2 (0) = −2.67586 16249,

Fc (T ) = 29.98534 49252, λ1 (1) = 0.15146 20116, λ2 (1) = −0.125.

It is obvious that these three solutions correspond to the two local minima and one local maximum shown in Figure 5.2. Note that λ2 (1) changes sign when passing from solution 1 to solutions 2 and 3. Now let us show that solution 3 indeed provides a local minimum. The respective optimal control, state, and adjoint variables are displayed in Figure 5.3. The

397

SSC FOR CONTROL PROBLEMS WITH FREE FINAL TIME u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 5 4 3 2 1 0

; ;

1 2

0.0

0.2

0.4

0.6

; ; ; ; ; ;

5 .................................................................... ................. .................................................. .......... ....... ..... ..... .... .... .... .... .... . . . .. ... .... ... .... ... ... .... ... .... ... . . .. ... ... ... ... ... ... .... .... .... . ... ..... .... .... ...................

1 2 3 4 5

4 3 2 1

0.2

0.4

0.6

0.8

; ; ; ; ; ;

0 1 2 3 4

6 0.0

5

2 4 6 8

. . .. . . . . . . . . . . . . . . . . . . . . . . . . . . .

0.0

1 0

....................................... .......... ..... ........ ... ...... ... ..... ... ..... .. ..... .. .... .. .... .. .... . . ... . ..... .. . ..... . . ..... .. . ..... . . ...... . .......... . . ........... . . ................. ........... .. ...................................................................... . . .

6

1.0

0.2

0.4

0.6

0.8

1.0

2

2

; ; ; ; ;

1.0

x2

x1 1 0

0.8

2 ....................................................................................................................... ............. .................................... ......... ......... ...... ..... .... ..... . . . . .... .... . . . ... .... .... .... .... ... ... ... .. ... . . .. ... ... .. ... .. ... ... ... ... ... . .. ... ...

10 0.0

0.2

0.4

0.6

0.8

1 0

; ; ; ;

1 2 3

1.0

........................... ......... ......... ......... ...... .......... ...... ............. ...... .................... ....... ............................................................................ ......... .......... ............. ........... . . . . . . . .. .... .... . .. ... . . .. .. . .. . . . .. . .. . . . . . . . . .. .. . . ....

4 0.0

0.2

0.4

0.6

0.8

1.0

Fig. 5.3. Optimal control, state, and adjoint variables for c = 1/16 and T = 3.5153536980.

˜ 0 [τ ] ≡ 2 · T > 0 . Next Legendre–Clebsch condition (A2) trivially holds in view of H uu we verify that the Riccati equations (3.24)–(3.26) have a bounded solution such that the boundary conditions (3.31) hold in the relaxed forms qT (0) ≥ 0 and qT (1) < 0 . The matrix (3.23) becomes   q 1 q 2 r1 Q R = (5.9) =:  q2 q4 r2  , Q R∗ q r1 r2 q T for which we evaluate the Riccati equations (3.24)–(3.26) as

(5.10)

dq1 /dτ dq2 /dτ dq4 /dτ dr1 /dτ dr2 /dτ dqT /dτ

T · [ 2q2 − 2 + 8q22 ], T · [ −q1 − (1.4 − 0.42x22 )q2 + q4 + 8q2 q4 ] , T · [ −2(q2 + (1.4 − 0.42x22 )q4 ) + 0.84x2 λ2 + 8q42 ] , −q1 x2 − q2 (−x1 + 1.4x2 − 0.14x32 − 8λ2 ) + T r2 −2x1 + λ2 + 8T q2 r2 , = −q2 x2 − q4 (−x1 + 1.4x2 − 0.14x32 − 8λ2 ) − T r1 −T · (1.4 − 0.42x22 )r2 − λ1 − λ2 (1.4 − 0.42x22 ) + 8T q4 r2 , = −2[ r1 x2 + r2 (−x1 + 1.4x2 − 0.14x32 − 8λ2 ) ] + 8T r22 . = = = =

398


It suffices to find a bounded solution of these Riccati equations satisfying qT (0) = 0. After several trials, we were successful with the initial values q1 (0) = 2.00684 76891, q2 (0) = 0.47018 97048, q4 (0) = −0.35197 44265, r1 (0) = 0, r2 (0) = 0, qT (0) = 0, for which we get the final values, q1 (1) = 0, q2 (1) = 0, q4 (1) = 0, r1 (1) = −0.125 = λ2 (1), r2 (1) = 0.02353 79884, qT (1) = −r2 (1) < 0, )||∞ ≤ 3 for all τ ∈ [0, 1]. Hence Theorem 3.1 asserts that and the bound ||Q(τ solution 3 is indeed a weak local minimum. We mention that the sufficient conditions in Hull [12] can also be checked numerically for solution 3. In a similar way, we can test the optimality of solution 1 with T = 2.1946 079912. We obtain a bounded solution of the Riccati equation for initial values q1 (0) = 1.59068 73787, q2 (0) = 0.33016 84322, q4 (0) = −0.39917 54076, r1 (0) = 0, r2 (0) = 0, qT (0) = 0 and final values q2 (1) = 0, q4 (1) = 0, q1 (1) = 0, r1 (1) = 0.125 = λ2 (1), r2 (1) = −1.15193 36046, qT (1) = r2 (1) < 0. )||∞ ≤ 5 for all τ ∈ [0, 1]. These values yield the bound ||Q(τ The situation is different for solution 2, which provides a local maximum with respect to the final time T . All initial values qT (0) ≥ 0 that we tested produce a solution of the Riccati equation with qT (1) ≥ 0 . Though this test does not exclude optimality of the solution, it is a rather strong indication of nonoptimality. Thus we may draw the conclusion that solution 2 behaves like a saddle point solution, which is a local minimum with respect to control for every fixed time but a local maximum with respect to final time. Case c = 0. Here the situation is more complicated since Figure 5.2 does not indicate a distinctive local minimum. The code BNDSCO of [29] provides, e.g., the following two solutions: Solution 1:

T = 2.29903 95815, λ1 (0) = −9.00409 78999, λ2 (0) = −2.67325 14651,

Fc (T ) = 29.9218 12616, λ1 (1) = 1.00813 43679, λ2 (1) = 0.

Solution 2:

T = 4.50237 87337, λ1 (0) = −9.00247 06599, λ2 (0) = −2.67303 08344,

Fc (T ) = 29.75107 51464, λ1 (1) = −0.0044561 06307, λ2 (1) = 0.

Solving the Riccati equation (5.10), e.g., for the final time T = 2.2990 395815, we find that the initial value qT (0) = 0 produces the final value qT (1) = 0 in all tested cases. Thus the sign conditions (3.29) cannot be verified, and the Riccati test is not able to detect whether solution 1 is a local minimum.

399


Fc (T )

47.5 47.0 46.5 46.0 45.5 45.0

Fc (T )

47.5

.. .. .. .. .. .. .. .. .. .. .. .. .. .. ... ... .. .. ... .. .. .. .. ... ... ... .... .... ..... ..... ..... ...... .......... .............. .............................. ................................................................................................................................................................................

47.0 46.5 46.0 45.5

T 3.5

4.0

4.5

5.0

45.0

5.5

.. .. .. .. .. .. .. .. .. .. .. .. .. .. ... .. .. .. .. .. ... .. .. .. ... ... ... ... .... ...... ..... ..... ..... .......... ........... ........................ ......................................................................................................................................................................................

3.5

4.0

4.5

5.0

T

5.5

Fig. 5.4. Graph of Fc (T ) for control constraint |u| ≤ 1: Cases c = 1/16 and c = 0.

5.2. Constrained optimal solutions. Now we consider solutions satisfying, in addition, the control constraint (5.4), −1 ≤ u(τ ) ≤ 1

∀ τ ∈ [0, 1] .

The augmented Hamiltonian (3.12) becomes, in view of (5.5), ˜ µ ˜ +µ ˜ x, u, λ, ˜ 0 (˜ ˜)= H (5.11) H(˜ x, u, λ) ˜1 (−u − 1) + µ ˜2 (u − 1) = T · [ c + u2 + x21 + λ1 x2 + λ2 (−x1 + x2 (1.4 − 0.14x22 ) + 4u) + µ1 (−u − 1) + µ2 (u − 1) ] , where µi := µ ˜/T, i = 1, 2, are the scaled multipliers. The state and adjoint equations agree with those given in (5.6). Again, we get the control law u(τ ) = −2λ2 (τ ) on interior arcs |u(τ )| < 1 . The unconstrained control u(τ ) depicted in Figure 5.3 suggests that the constrained control has one boundary arc with u(τ ) ≡ 1 and one boundary arc with u(τ ) ≡ −1. Thus we may assume the following solution structure of the optimal control:   1, 0 ≤ τ ≤ τ1       −2λ2 (τ ), τ1 ≤ τ ≤ τ2 u(τ ) = (5.12) . −1, τ2 ≤ τ ≤ τ3       −2λ2 (τ ), τ3 ≤ τ ≤ 1 The junction points τ1 , τ2 , τ3 are implicitly determined through the conditions that the control is continuous at these points. This leads to the junction conditions (5.13)

λ2 (τ1 ) = −0.5,

λ2 (τ2 ) = 0.5,

λ2 (τ3 ) = −0.5 .

Hence, on the interval [0, 1], we have to solve the multipoint BVP, which is composed by the state and adjoint equations (5.6) with control substituted from (5.12) as well as the boundary and junction conditions (5.7) and (5.13). The optimal value function F˜c (T ) for the constrained problem is depicted in Figure 5.4 for the values c = 1/16 and c = 0 . A distinctive minimum can only be detected in case c = 1/16 . Case c = 1/16. Again we use the code BNDSCO in [29] and obtain T = 4.54230 98018, 0.37589 55717, τ2 = λ1 (0) = −12.70813 77440, λ2 (0) = −4.59503 53190, 44.71797 06589. Fc (T ) =

τ1 = 0.22932 06694, τ3 = 0.63465 63122, λ1 (1) = 0.02860 41331, λ2 (1) = −0.125,

400

HELMUT MAURER AND HANS JOACHIM OBERLE u 1.0 0.5 0.0

; ;

0:5 1:0

. . .. . . .. . . .. . . .. . . .. . . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . . . .. . . . ... .. . . . . .. . . . . .. . ... . . . .. . . . . . .. . . . . .. . . . . . .. . . . . . .. . . . . .. . . . . .. . . . .. . . . . . .. . . . . .. . . .. . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . .

0.0

0.2

0.4

0.6

4.0

........................................................ ....... ................. ...... .................................................................. ..... ..... ..... ... . . . . ... .... . . . . .... ... ... ... ... ... ... .. ... . . .. ... ... ... ... .. ... ... ... ... ... . . .... .... .... ..... ....... ...................

; ; ;

2:0 4:0 6:0

2.0 0.0

0.2

0.4

0.6

0.8

; ;

2:0

0.0

4:0

1.0

0

; ; ; ; ; ;

2 4 6 8

10 12

0.2

0.4

0.6

0.8

1.0

2

...................................................... ......... .............. ....... ..................... ....................................................................... ..... ..... ..... ..... . . . .. ..... . . . ... ... . . ... ... .. ... ... ... .. .. .. . .. ... .. ... ... .. ... . . ... ... ... ... .... ..... ..... ....

0.0

...................... ...... .... ... .... .... ... .... .. .... ... .... ..... .. .... .. ..... .. ..... . ..... .. ..... . . .... .. .... . .... .. .... .. . .... ... .. ... .. ... . ..... .................................... .. . ......... ................... . ................................................. ... . . .. . . ... . .. .. .. ... .. .. .. .. . . .. ..

0.0

1 2

1.0

x2

x1 2.0 0.0

0.8

0.2

0.4

0.6

0.8

2 0

; ; ; ;

2 4 6

1.0

.......................... ...... ....... ..... ..... ..... .... ....... ..... .......... ..... ............... ........ ................................................................................................ ........... ............ . . . . . . .... ... . . . ... . . .. . .. . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. .. ... .....

8 0.0

0.2

0.4

0.6

0.8

1.0

Fig. 5.5. Optimal control, state, and adjoint variables for c = 1/16 and T = 4.54230 98018.

The corresponding optimal control, state, and adjoint variables are shown in Figure 5.5. In order to check the sufficient conditions in Theorem 3.1, we try to find a bounded solution of the Riccati equations (5.10) on the interior arcs [τ1 , τ2 ] and [τ3 , 1] and the modified Riccati equations (3.30) on the boundary arcs [0, τ1 ] and [τ2 , τ3 ] . The modified Riccati equations yield the following linear equations:

(5.14)

dq1 /dτ dq2 /dτ dq4 /dτ dr1 /dτ dr2 /dτ dqT /dτ

= = = = =

T · 2(q2 − 1) , T · [ −q1 − (1.4 − 0.42x22 )q2 + q4 ] , T · [ −2(q2 + (1.4 − 0.42x22 )q4 ) + 0.84x2 λ2 ] , −q1 x2 − q2 (−x1 + 1.4x2 − 0.14x32 + 4u) + T r2 − 2x1 + λ2 , −q2 x2 − q4 (−x1 + 1.4x2 − 0.14x32 + 4u) − T r1 −T · (1.4 − 0.42x22 )r2 − λ1 − λ2 (1.4 − 0.42x22 ) , = −2 [ r1 x2 + r2 (−x1 + 1.4x2 − 0.14x32 + 4u) ] .

We wish to find a solution satisfying the terminal values q1 (1) = q2 (1) = q4 (1) = 0 . A bounded solution of the Riccati equations then is obtained for the initial values q1 (0) = 2.39837121, q2 (0) = 0.89021498, q4 (0) = −1.26573031, r1 (0) = 0, r2 (0) = 0, qT (0) = 0,

401

SSC FOR CONTROL PROBLEMS WITH FREE FINAL TIME q1

.......................... .............. .... ... ... ... ... ... ... ... ... ........... . ... .................... . . ... . ........... .... ........ .. .... ..... .. . .... ... .. .... .... .... ............... ... ... .... .... .... .... .... .... ...... .......... ...................... ....... .... ... ... ... ... . ... .. ... .. .. .. ... .. .. .

2.5 2.0 1.5 1.0 0.5 0.0

0.0

0.2

0.4

0.6

0.8

q2 1.5 .... ... .... .. .... . .. .. .. . .. . .. . . .. . . . . . ........ . . . . ....... ........... .. .... .... . . ... ..... . . . ... ... . . .. ... . . ... ... .. . ... . .. . . .. .. . ... . . . .. ... .. . . . . ... .. . .. .. .. . . .. ........ ... ... .. ....... ......... ... .. .. . ... ...... . ... ..................................... .. ... ... ... ... . ..... . . ... ................... ... ... ... .... ..... ...........

1.0

0.5

0.0

1.0

0.0

;

0:1

;

0:2

;

0:3

2.0

................................ ...... ................................. ........... ......................................... ...... ..... .... .... .... ..... ..... .... ..... ..... .... .... .... .. .... ... ..... ... .... ... .... .. ... ... ... .. .. ... ... . .. . . .. . ... .. .. .. ... .. ... . .. ... .. ... .. ... .. ... .. .... ..... ...........

1.5 1.0 0.5

0.2

0.4

0.6

0.8

0.0

; ; ;

0:5 1:0

0.0

1:5

1.0

;

0:1

0.8

0.0 ......................... ..... ... ... .. .. .. .. . . .. . .. .. .. ... ... .. .................................................................................................. ... ..... . . .... . .... ... ..... ... .... .... .... .... ................................ .... .... .... .... ..... ..... .... ....... .......... ...........

;

0.2

0.4

0.6

0.8

1.0

;

0:1

;

0:2

..................................... ................. ............ .......... ........ ........ ........ ........ ...... ...... ...... ..... .... ..... .... ..... ..... .... .... .... ..... ..... ..... ...... ...... ..... ..... .... ..... .... .... ..... .... .... .... .... .... .... ..... .... . . . .................

0:2 0.0

1.0

qT

r2

0.0

0.6

...... . . . .. .. .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . ................ . . . .. .... ... . . ... .... .. . ... . . ... .. . ... ... . .. . ...................... . ... ... ............................ . .. . ... ......... . ... . . . ..... ..... . . ... ...... ........... . ... . ........ . ......................... . . ............................. ............................. . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..... .

0.0

0.2

0.1

0.4

q4

r1 0.0

0.2

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Fig. 5.6. Solutions q1 , q2 , q4 , r1 , r2 , qT of Riccati equations (5.10) and (5.14) for c = 1/16 and T = 4.54230 98018.

which produce the desired terminal values for q1 , q2 , q4 and r1 (1) = −0.125 = λ2 (1), r2 (1) = 0.146395866, qT (1) = −r2 (1) < 0 . )||∞ ≤ 4 for all τ ∈ [0, 1], which can be seen For these values, we get the bound ||Q(τ in Figure 5.6. It is interesting to note that it was not possible to obtain a bounded solution of the Riccati equation (5.10) on the whole interval. Thus the Riccati test developed in section 3 in its weaker form considerably facilitates the numerical check of SSC. Case c = 0. The code BNDSCO yields the solution T = 5.15173 31990, τ2 = 0.33186 51046, λ1 (0) = −12.70087 48310, λ2 (0) = −4.58996 79054, F˜c=0 (T ) = 44.70866 79043.

τ1 = 0.20192 87957, τ3 = 0.55797 03824, λ1 (1) = 0.09649 19382, λ2 (1) = 0,

Evaluating the Riccati equations (5.10) and (5.14) along this specific solution, we were not able to find a bounded solution satisfying the sign conditions qT (0) ≥ 0 and qT (1) < 0 . This fact confirms our impression gained from Figure 5.4 that no local minimum can be identified for c = 0.

402


Acknowledgments. We are indebted to D. Augustin, Ch. Beckmann, and J. Strade for numerical assistance with the examples. We would like to thank an anonymous reviewer for helpful remarks. REFERENCES [1] A. E. Bryson and Y. C. Ho, Applied Optimal Control, revised printing, Hemisphere, New York, 1975. ¨skens and H. Maurer, SQP-methods for solving optimal control problems with con[2] Ch. Bu trol and state constraints: Adjoint variables, sensitivity analysis and real-time control, J. Comput. Appl. Math., 120 (2000), pp. 85–108. [3] R. Bulirsch, Die Mehrzielmethode zur numerischen L¨ osung von nichtlinearen Randwertproblemen und Aufgaben der optimalen Steuerung, Report of the Carl-Cranz Gesellschaft, Oberpfaffenhofen, Germany, 1971. [4] M. Chamberland and V. Zeidan, Second order necessity and sufficiency theory for the free final time problem, in Proceedings of the 31st IEEE Conference on Decision and Control, Tucson, AZ, 1992, pp. 1518–1525. [5] A. L. Dontchev, W. W. Hager, A. B. Poore, and B. Yang, Optimality, stability and convergence in nonlinear control, Appl. Math. Optim., 31 (1995), pp. 297–326. [6] J. C. Dunn, Second-order optimality conditions in sets of L∞ functions with range in a polyhedron, SIAM J. Control Optim., 33 (1995), pp. 1603–1635. [7] J. C. Dunn, On L2 sufficient conditions and the gradient projection method for optimal control problems, SIAM J. Control Optim., 34 (1996), pp. 1270–1290. [8] J. C. Dunn, On second order sufficient conditions for structured nonlinear programs in infinitedimensional function spaces, in Mathematical Programming with Data Perturbations, A. V. Fiacco, ed., Lecture Notes in Pure and Appl. Math. 195, Marcel Dekker, New York, 1998, pp. 83–107. [9] U. Felgenhauer, Diskretisierung von Steuerungsproblemen unter stabilen Optimalit¨ atsbedingungen, Habilitationsschrift, Department of Mathematics, Technische Universit¨ at Cottbus, Cottbus, Germany, 1999. [10] U. Felgenhauer, On smoothness properties and approximability of optimal control functions, Ann. Oper. Res., 101 (2001), pp. 23–42. [11] M. Hestenes, Calculus of Variations and Optimal Control Theory, John Wiley, New York, 1966. [12] D. G. Hull, Sufficient conditions for a minimum of the free-final-time optimal control problem, J. Optim. Theory Appl., 68 (1991), pp. 275–287. [13] D. H. Jacobson and D. Q. Mayne, Differential Dynamic Programming, American Elsevier Publishing, New York, 1970. [14] H. J. Kelley, Gradient theory of optimal flight paths, ARS Journal, 30 (1960), pp. 947–954. [15] H. J. Kelley, Method of gradients, in Optimization Techniques, G. Leitmann, ed., Academic Press, New York, 1962, pp. 205–254. [16] P. Kenneth and R. McGill, Two-point boundary-value problem techniques, in Advances in Control Systems, Vol. 3, C. T. Leondes, ed., Academic Press, New York, 1966, pp. 69–109. [17] D. F. Lawden, Optimal Trajectories for Space Navigation, Academic Press, New York, 1967. [18] K. Malanowski, Two-norm approach in stability and sensitivity analysis of optimization and optimal control problems, Adv. Math. Sci. Appl., 2 (1993), pp. 397–443. [19] K. Malanowski, Stability and sensitivity of solutions to nonlinear optimal control problems, Appl. Math. Optim., 32 (1994), pp. 111–141. [20] K. Malanowski, Sufficient optimality conditions for optimal control subject to state constraints, SIAM J. Control Optim., 35 (1997), pp. 205–227. [21] K. Malanowski and H. Maurer, Sensitivity analysis for parametric control problems with control-state constraints, Comput. Optim. Appl., 5 (1996), pp. 253–283. [22] H. Maurer, Second order sufficient conditions for control problems with free final time, in Proceedings of 3rd European Control Conference, A. Isidori et al., eds., Rome, Italy, 1995, pp. 3602–3606. [23] H. Maurer and D. Augustin, Second order sufficient conditions and sensitivity analysis for the controlled Rayleigh problem, in Parametric Optimization and Related Topics IV, J. Guddat, H. Th. Jongen, F. Nozicka, G. Still, and F. Twilt, eds., Peter Lang Verlag, Frankfurt, Germany, 1996, pp. 245–259. [24] H. Maurer and H. J. Pesch, Solution differentiability for parametric nonlinear control problems with control-state constraints, J. Optim. Theory Appl., 86 (1995), pp. 285–309.


403

[25] H. Maurer and S. Pickenhain, Second order sufficient conditions for optimal control problems with mixed control-state constraints, J. Optim. Theory Appl., 86 (1995), pp. 649–667. [26] A. A. Milyutin and N. P. Osmolovskii, Calculus of Variations and Optimal Control, Transl. Math. Monogr. 180, AMS, Providence, RI, 1998. [27] H. G. Moyer and G. Pinkham, Several trajectory optimization techniques. II. Application, in Computing Methods in Optimization Problems, A. V. Balakrishnan and L. W. Neustadt, eds., Academic Press, New York, 1964, pp. 91–105. [28] L. W. Neustadt, Optimization: A Theory of Necessary Conditions, Princeton University Press, Princeton, NJ, 1976. [29] H. J. Oberle and W. Grimm, BNDSCO—A Program for the Numerical Solution of Optimal Control Problems, Internal Report 515–89/22, Institute for Flight Systems Dynamics, DLR, Oberpfaffenhofen, Germany, 1989. [30] H. J. Oberle and K. Taubert, Existence and multiple solutions of the minimum-fuel orbit transfer problem, J. Optim. Theory Appl., 95 (1997), pp. 241–262. [31] N. P. Osmolovskii, Quadratic conditions for nonsingular extremals in optimal control (a theoretical treatment), Russian J. Math. Phys., 2 (1995), pp. 487–516. [32] N. P. Osmolovskii, Second-order conditions for broken extremals, in Proceedings of the Conference on Calculus of Variations and Optimal Control, A. Ioffe, S. Reich, and I. Shafir, eds., Chapman & Hall/CRC, Boca Raton, FL, 2000, pp. 198–216. [33] H. J. Pesch, Real-time computation of feedback controls for constrained optimal control problems. I. Neighbouring extremals, Optimal Control Appl. Methods, 10 (1989), pp. 129–145. [34] H. J. Pesch, Real-time computation of feedback controls for constrained optimal control problems. II. A correction method based on multiple shooting, Optimal Control Appl. Methods, 10 (1989), pp. 147–171. [35] S. Pickenhain, Sufficiency conditions for weak local minima in multidimensional optimal control problems with mixed control-state restrictions, Z. Anal. Anwendungen, 11 (1992), pp. 559–568. [36] T. Tun and T. S. Dillon, Extensions of the differential dynamic programming method to include systems with state dependent control constraints and state variable inequality constraints, J. Appl. Sci. Engrg. A, 3 (1978), pp. 171–192. [37] V. Zeidan, The Riccati equation for optimal control problems with mixed state-control constraints: Necessity and sufficiency, SIAM J. Control Optim., 32 (1994), pp. 1297–1321.