- Email: [email protected]
Time-optimal aircraft trajectories in climbing phase and singular perturbations
Proceedings of the 20th World Congress Proceedings of 20th The International Federation of Congress Automatic Control Proceedings of the the 20th World World Congress Proceedings of the 20th World Congress The International Federation of Automatic Control Available online at www.sciencedirect.com Toulouse, France, July 9-14, 2017 The International Federation of The International Federation of Automatic Automatic Control Control Toulouse, France, July 9-14, 2017 Toulouse, France, July 9-14, 2017 Toulouse, France, July 9-14, 2017
ScienceDirect
IFAC PapersOnLine 50-1 (2017) 1625–1630
Time-optimal aircraft trajectories in Time-optimal aircraft trajectories in Time-optimal aircraft trajectories in Time-optimal aircraft trajectories in climbing phase and singular perturbations climbing phase and singular perturbations climbing phase and singular perturbations climbing phase and singular perturbations
Olivier Cots ∗∗ Joseph Gergaud ∗∗ Damien Goubinat ∗∗ ∗∗ Olivier Cots ∗∗ Joseph Gergaud ∗∗ Damien Goubinat ∗∗ Olivier Olivier Cots Cots Joseph Joseph Gergaud Gergaud Damien Damien Goubinat Goubinat ∗∗ ∗ ∗ Toulouse Univ., INP-ENSEEIHT-IRIT, UMR CNRS 5505, 2 rue ∗ Toulouse Univ., INP-ENSEEIHT-IRIT, UMR CNRS 5505, 2 rue ∗ Toulouse31071 Univ., INP-ENSEEIHT-IRIT, CNRS Camichel, France (e-mail: UMR [email protected] Toulouse31071 Univ.,Toulouse, INP-ENSEEIHT-IRIT, UMR CNRS 5505, 5505, 2 2 rue rue& Camichel, Toulouse, France (e-mail: [email protected] & Camichel, 31071 Toulouse, France (e-mail: [email protected] [email protected]). Camichel, 31071 Toulouse, France (e-mail: [email protected] & & [email protected]). ∗∗ [email protected]). SA, 105 av du General Eisenhower, B.P. 1147, ∗∗ Thales Avionics [email protected]). ∗∗ Thales Avionics SA, 105 av du General Eisenhower, B.P. 1147, ∗∗ Thales Avionics SA, 105 du Eisenhower, 31047 Toulouse Cedex, [email protected]) Thales Avionics SA,France 105 av av(e-mail: du General General Eisenhower, B.P. B.P. 1147, 1147, 31047 Toulouse Cedex, France (e-mail: [email protected]) 31047 Toulouse Cedex, France (e-mail: [email protected]) 31047 Toulouse Cedex, France (e-mail: [email protected]) Abstract: This article deals with the singularly perturbed aircraft minimum time-to-climb Abstract: This article deals with the singularly perturbed aircraft minimum time-to-climb Abstract: This article deals the perturbed aircraft minimum time-to-climb problem. First, introduce reduced-order problem with affine dynamics with respect to the Abstract: Thiswe article dealsaa with with the singularly singularly perturbed aircraft minimum time-to-climb problem. First, we introduce reduced-order problem with affine dynamics with respect towith the problem. First, we introduce a reduced-order problem with affine dynamics with respect the control and analyze it with the tools from geometric control: maximum principle combined problem. First, we introduce a reduced-order problem with affine dynamics withcombined respect to towith the control and analyze it with the tools from geometric control: maximum principle control order and analyze analyze it with with the tools toolsThen, from geometric geometric control: maximum principle combined combined with second optimality conditions. we compute candidates as minimizers using multiple control and it the from control: maximum principle with second order optimality conditions. Then, we compute candidates as minimizers using multiple second order Then, we compute as multiple shooting andoptimality homotopyconditions. methods for both that we briefly compareusing numerically. second order optimality conditions. Then, we problems, compute candidates candidates as minimizers minimizers using multiple shooting and homotopy methods for both problems, that we briefly compare numerically. shooting and homotopy methods for both problems, that we briefly compare numerically. Particular attention is paid to the singular extremals and we discuss about their local optimality. shooting and homotopy methods for both problems, that we briefly compare numerically. Particular attention is paid to the singular extremals and we discuss about their local optimality. Particular attention is to extremals and we about local optimality. Particular attention is paid paidFederation to the the singular singular extremals and Hosting we discuss discuss about their their local optimality. © 2017, IFAC (International of Automatic Control) by Elsevier Ltd. All rights reserved. Keywords: Aircraft trajectory, geometric optimal control, singular perturbation. Keywords: Aircraft trajectory, geometric optimal control, singular perturbation. Keywords: Keywords: Aircraft Aircraft trajectory, trajectory, geometric geometric optimal optimal control, control, singular singular perturbation. perturbation. g Å ã βR 1. INTRODUCTION g Å Θ(h) ã ã βR g 1. INTRODUCTION INTRODUCTION Å g , 1. Θ(h) := ã Å P (h) P βR 0 Θ(h) βR 1. INTRODUCTION := ,, P (h) P Θ(h) Θ 0 0 := P (h) P 0 Θ0 , P (h) := P0 Θ A flight is composed of several phases which are take-off, 0 := Θ(h) Θ − βh, A flight is composed of several phases which are take-off, Θ 0 0 A flight is composed several phases which are take-off, climb, cruise, descent,of approach and landing. := Θ(h) Θ − βh, A flight is composed of several phases which are take-off, 0 Θ(h) − climb, cruise, cruise, descent, approach approach and landing. landing. := Θ Θ(h) := ΘP00 (h) − βh, βh, climb, climb, cruise, descent, descent, approach and and landing. (h) . ρ(h) := P P (h) := . ρ(h) P (h) RΘ(h) := ρ(h) . RΘ(h) . ρ(h) := RΘ(h) RΘ(h) constants: g is the graviThe remaining data are positive The remaining data are wing positive constants: is the the graviThe remaining data are positive constants: ggg is gravitational constant, S the area, R the specific constant The remaining data are wing positive constants: is the gravitational constant, S the area, R the specific constant tational S the wing area, R specific of air, β constant, the thermical gradient and Pthe pressure 0, Θ 0 , the constant tational constant, S the wing area, R the specific constant of air, air, β temperature the thermical thermical gradient and P P0 , Θ0 ,, the the pressure of the and and theβ atgradient the sea level. of air, β temperature the thermical gradient and P00 ,, Θ Θ00 , the pressure pressure the at the sea level. In this article, we are interested in the time optimal control and and the temperature at the sea level. In this article, we are interested in the time optimal control and the temperature at the sea level. In this article, are in time control of an aircraft during its climbing This phase is The underlying optimal control problem consists of miniIn article, we we are interested interested in the thephase. time optimal optimal control The underlying optimal control problem consists of cruise miniof this an aircraft aircraft during itsdynamics climbing phase. This phase is The underlying optimal of miniof an during its climbing phase. This phase is mizing the transfer time,control i.e. theproblem time to consists reach the determined by its own governed by the fourunderlying optimal control problem consists of cruise miniof an aircraft during itsdynamics climbinggoverned phase. This phase is The mizing the transfer time, i.e. the time to reach the determined by its own by the fourmizing the transfer time, i.e. the time to reach the cruise determined by its own dynamics governed by the fourphase, with fixed initial and final states. It is well known dimensional system mizing the transfer time, i.e. the time to reach the cruise determined by its own dynamics governed by the fourphase, with fixed initial and final states. is well known dimensional system system phase, with fixed and final It is known dimensional that the dynamics (1)–(4) slowIt mass m) phase, with fixed initial initial and contains final states. states. It (the is well well known dimensional system dh that the dynamics (1)–(4) contains slow (the mass m) 1 that the dynamics (1)–(4) contains slow (the mass m) dh = v sin γ (1) and fast (the air slope γ) variables , see Ardema (1976) the (the dynamics (1)–(4) contains11 , slow (the mass m) dh = vv sin sin γ γ (1) that and fast air slope γ) variables see Ardema (1976) dh dt = (1) 1 and fast (the air slope γ) variables , see Ardema (1976) Nguyen (2006). Thisγ)time scale separation is normally dt = v sin γ (1) and fast (the air slope variables , see Ardema (1976) 2 dt and Nguyen (2006). This time scale separation is normally 1 ρ(h)Sv T (h) dv dt = T (h) − 1 ρ(h)Sv 22 C + C u2 − g sin γ (2) treated and This time normally by a(2006). singular perturbation analysis byis dv and Nguyen Nguyen This time scale scale separation separation is replacing normally 1 2 2 (h) dv by aa(2006). singular perturbation analysis by replacing − 121 ρ(h)Sv =T CD +C CD u22 − − gg sin sin γ γ (2) (2) treated ρ(h)Sv Tm (h) − dv dt = m 2 C treated by singular perturbation analysis by replacing D1 + D2 u the fast state equation with the following: treated by a singular perturbation analysis by replacing D1 + C D2 u 2 m dt m − = C − g sin γ (2) the equation with D1 D2 dm dt m the fast fast state state equation with the thegfollowing: following: ms (v) T22(h) m dt = −C m with the following: dγequation 1 ρ(h)Sv dm (3) the fast state dm dγ 1 ρ(h)Sv (4’) u − gg cos γ, ε > 0. ε = = −C (v) T (h) (3) dm dt = −Css (v) T (h) 1 dγ ρ(h)Sv (3) cos γ, γ, εε > > 0. 0. (4’) − εε dγ 12 ρ(h)Sv vg cos m u dt = dt = −C (3) s (v) T (h) g = u − (4’) dt 1 ρ(h)Sv dγ v m γ, ε > 0. (4’) ε dt dt 2 2 dt = 2 time-to-climb m u − vv cos g 1 ρ(h)Sv dγ cos γ, (4) u − = Let (P ) denote the problem with the addi2 m dt ε dγ 2 cos γ, γ, (4) Let u− − vgg cos = 11 ρ(h)Sv dγ 2 the time-to-climb problem with the addi2 ρ(h)Sv m u dt = Let (P ) denote ε (4) 2 ) denote the time-to-climb problem with the addi(P ε > 0. From the control theory, m u − vv cos γ, dt = 22 m (4) tional Let (Pεεartificial ) denote parameter the time-to-climb problem with the addidt artificial εε > 0. From the theory, where v x := (h, v, m, γ) is composed of tional m variable dt the2 state tional artificial parameter > 0. From the control theory, for a fixed valueparameter of ε > 0, the candidates as control minimizers are tional artificial parameter ε > 0. From the control theory, := where the state variable x (h, v, m, γ) is composed of for aa fixed value of εε > 0, the candidates as minimizers are := where the state (h, v, γ) is composed of the altitude, thevariable true airx speed, them, weight and the air for fixed value of > 0, the candidates as minimizers are := where the state variable x (h, v, m, γ) is composed of selected among a set of extremals, solution of a Hamiltofor a fixed value aofset ε >of0,extremals, the candidates as minimizers are the altitude, the true trueInair airthis speed, thethe weight and the the air air selected among solution of a Hamiltothe altitude, the speed, the weight and slope of the aircraft. model, lift coefficient is selected among a set of extremals, solution of a Hamiltothe altitude, the true air speed, the weight and the air nian system given by the Pontryagin Maximum Principle selected among a set of extremals, solution of a Hamiltoslope of the aircraft. In this model, the lift coefficient is system given by the Pontryagin Maximum Principle slope of In lift coefficient is taken as the the aircraft. control variable u. The the thrust (h) and the nian system given Pontryagin Maximum Principle slope of aircraft. In this this model, model, the liftT is nian (PMP), see Pontryagin et al. (1962). The application of the nian system given by by the the Pontryagin Maximum Principle taken as the the control variable u. The thrust Tcoefficient (h) and and the (PMP), see Pontryagin et al. (1962). The application of the taken as the control variable thrust T (h) the fuel flow Cs (v) are given by u. theThe Base of Aircraft DAta (PMP), see Pontryagin et al. (1962). The application of the taken as the control variable u. The thrust T (h) and the PMP leads to define a Boundary Value Problem (BV ε) (PMP), see Pontryagin et al. (1962). The application of P the fuel flow C (v) are given by the Base of Aircraft DAta s (v) are given by the Base of Aircraft DAta PMP leads to define a Boundary Value Problem (BV P fuel flow C ε) (BADA) model: s (v) are given by the Base of Aircraft DAta PMP leads to define a Boundary Value Problem (BV P fuel flow C which can be solved using shooting methods combined s PMP leads to define ausing Boundary Valuemethods Problemcombined (BV Pεε )) (BADA) model: which can be solved shooting (BADA) model: which can methods combined (BADA) model: with numerical tests tousing check shooting second-order optimality conh which can be be solved solved using shooting methods combined with numerical tests to check second-order optimality con+ h22 CT3 ), T (h) := CT1 (1 − h with numerical tests to check second-order optimality conh ditions. 2 with numerical tests to check second-order optimality con:= (1 − + h C ), T (h) C h CT2 + h2 CTT3 ), T ditions. T 3 ), := C ditions. (1 − −C + h C T (h) (h) := CTT111 (1 T T3 2 v C ditions. T When dealing with singular perturbation, one classical Cs (v) := Cs1 (1 + CvvT22 ), When dealing with one classical := C + C Cvs2 ), When dealing with singular perturbation, one s (v) s1 (1 approach consists in singular applying perturbation, the PMP to problem (Pε ) := (v) C (1 + ), C When dealing with singular perturbation, one classical classical s s ), Cs (v) := Cs11 (1 + C approach consists in applying the PMP to problem (P s ε) Cs22 approach consists in applying the PMP to problem (P ε) C approach consists in applying the PMP to problem (Pthe s2 CDi (drag coefficients) where the constants CTi , Csi , with 1 The altitude ε) h and the true air speed v are fast compare to , C , with C (drag coefficients) where the constants C 1 T s D i i i The altitude h and the true air speed v are fast compare to the , C , with C (drag coefficients) where the constants C 1 are specific to the aircraft. The International Standard T s D i i i mass but slow compare to the air slope. In this article, we consider The altitude altitude h h and and the the true true air air speed speed v v are are fast fast compare compare to to the the withInternational CDi (drag coefficients) where the constants CTi , Csi ,The 1 The are specific to the aircraft. to the air slope. In this article, we consider are specific the aircraft. The International Standard Atmospheric model of theStandard pressure mass h andbut v asslow slowcompare variables. mass but slow compare to the air slope. In this article, we consider are specific to to the provides aircraft. the Theexpression International Standard mass but slow compare to the air slope. In this article, we consider Atmospheric model provides the expression of the pressure h and v as slow variables. 2 Atmospheric model of pressure P (h), the temperature Θ(h) the andexpression the air density See section 4 for the exact definition of (Pε ), ε > 0. v variables. Atmospheric model provides provides the expression of the theρ(h): pressure 22hh and and v as as slow slow variables. P (h), the temperature Θ(h) and the air density ρ(h): See section 4 for the exact definition of (Pε ), ε > 0. P (h), the temperature Θ(h) and the air density ρ(h): See section 4 for 2 P (h), the temperature Θ(h) and the air density ρ(h): See section 4 for the the exact exact definition definition of of (P (Pεε ), ), εε > > 0. 0. Copyright 1661Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © © 2017 2017, IFAC IFAC (International Federation of Automatic Control) Copyright © 2017 IFAC 1661 Copyright © 2017 1661 Peer review responsibility of International Federation of Automatic Copyright © under 2017 IFAC IFAC 1661Control. 10.1016/j.ifacol.2017.08.327
Proceedings of the 20th IFAC World Congress 1626 Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
and then solve the associated singular boundary value problem with matching asymptotic expansion techniques, see Ardema (1976) and Mo¨ıss´eev (1985). In section 3, to approximate the solutions of problem (Pε ), ε > 0, we consider a reduced-order problem denoted 3 (P0 ), where ε is put to zero and where γ is taken as the control variable. This technique is also applied in Nguyen (2006) and is the first step of a second approach presented in Mo¨ıss´eev (1985). The problem (P0 ) has one state variable less than (Pε ), ε > 0, and can be tackled with the tools from geometric optimal control. This allow first to study the classification of the extremals, solution of the PMP, and then to check, on the singular arcs, second-order conditions of optimality.
with p0 ≤ 0 and (¯ p(·), p0 ) = (0, 0). Definition 1. We call extremal a triplet (x(·), p(·), u(·)) where (x(·), p(·)) is an integral curve of H satisfying (6). It is called a BC-extremal if it satisfies also (7) and (8).
We complete this analysis showing numerically in section 4 that problem (P0 ) from section 3 is a good approximation of problem (Pε ) when ε = 1. We use the HamPath software, Caillau et al. (2011), combining multiple shooting techniques with homotopy methods to solve the family of boundary value problems (BV Pε ) for ε ∈ [1 , 100], starting from the simpler problem when ε = 100.
We have the following characterization of singular controls which allows a practical computation. Proposition 3. The control u(·) with its associated trajectory x(·) are singular on [0 , T ] if and only if there exists a non zero adjoint p(·) such that (x(·), p(·), u(·)) is solution a.e. on [0 , T ] of (9) x˙ = ∂p H, p˙ = −∂x H, 0 = ∂u H.
2. GEOMETRIC OPTIMAL CONTROL
If u(·) is singular then its associated adjoint satisfies for each 0 < t ≤ T , p(t) ⊥ im dEx0 ,t (u(·)). Definition 4. A singular extremal is a triple (x(·), p(·), u(·)) solution of (9). It is called:
2.1 The time optimal control problem One considers a time optimal control problem given by the following data. A control domain U ⊂ Rm . A smooth control system: x(t) ˙ = F (x(t), u(t)), (5) with x(t) ∈ Rn , u(t) ∈ U , m ≤ n. An initial state x0 ∈ Rn and a final state xf ∈ Rn . The optimal control problem can be written tf | Ex0 ,tf (u(·)) = xf , (Ptmin ) min tf >0, u(·)∈Utf
where Utf := L∞ ([0 , tf ], U ) is the set of admissible controls 4 and where Ex0 ,tf is the end-point mapping defined by: Ex0 ,t : Ut → Rn , Ex0 ,t (u(·)) := x(t, x0 , u(·)), where t → x(t, x0 , u(·)) is the trajectory solution of (5) with control u(·) and such that x(0, x0 , u(·)) = x0 . 2.2 Pontryagin Maximum Principle Let define the pseudo-Hamiltonian: H : Rn × (Rn )* × U −→ R (x, p, u) −→ H(x, p, u) := p, F (x, u) . By virtue of the maximum principle, see Pontryagin et al. (1962), if (¯ u(·), t¯f ) is solution of (Ptmin ) then the associated trajectory x ¯(·) is the projection of an absolutely continuous integral curve (¯ x(·), p¯(·)) of H := (∂p H, −∂x H) such that the following maximization condition holds for almost every t ∈ [0 , t¯f ]: x(t), p¯(t), u). (6) H(¯ x(t), p¯(t), u ¯(t)) = max H(¯ u∈U
The following boundary conditions must be satisfied: (7) x ¯(t¯f ) = xf , 0 ¯ ¯ ¯ H(¯ x(tf ), p¯(tf ), u ¯(tf )) = −p , (8) 3
See section 3 for the exact definition of (P0 ). The set of admissible controls is the set of L∞ -mappings on [0 , tf ] taking their values in U such that the associated trajectory x(·) is globally defined on [0 , tf ]. 4
2.3 Preliminaries about singular trajectories A complete presentation of singular extremals – which play a crucial role in our analysis – can be found in Bonnard and Chyba (2003). In the following definitions and proposition we assume U = Rm , i.e. there is no control constraint. Definition 2. A control u(·) ∈ UT is called singular on [0 , T ] if the derivative of Ex0 ,T at u(·) is not surjective.
2 H is negative definite (strict Legendre• Regular if ∂uu Clebsch condition). • Strongly normal if im dEx(t1 ),t2 −t1 (u(·)|[t1 ,t2 ] ) is of corank one for each 0 < t1 < t2 ≤ T . • Exceptional if H = 0 and normal if H = 0.
2.4 Singular trajectories in the regular case In the regular case, using the implicit function theorem, we can solve locally the equation ∂u H = 0 and compute the smooth singular control as a function u ¯(z), z := (x, p). Plugging such u ¯ in H defines a true Hamiltonian h(z) := H(z, u ¯(z)). One can define the exponential mapping expx0 : (t, p0 ) → πx (exp(th)(x0 , p0 )), where h := (∂p h, −∂x h), x0 is fixed, πx is the x−projection (x, p) → x, and exp(th)(x0 , p0 ) is the extremal z(·) at time t solution of z(s) ˙ = h(z(s)), s ∈ [0 , t], z(0) = (x0 , p0 ). The domain of expx0 is R+ × P where P := p0 ∈ (Rn )* | h(x0 , p0 ) = −p0 . This leads to the following definition. Definition 5. Let z(·) := (x(·), p(·)) be a reference extremal solution of h. The time tc is said to be geometrically conjugate if (tc , p(0)) is a critical point of expx(0) . We have the following standard test: Proposition 6. The time tc is geometrically conjugate if and only if there exists a Jacobi field J(·) := (δx(·), δp(·)) solution of the variational equation δ z(t) ˙ = dh(z(t))·δz(t), vertical at 0 and tc , i.e. δx(0) = δx(tc ) = 0, and such that δx(·) ≡ 0 on [0 , tc ]. Let z¯(·) be a reference extremal. If the maximized Hamiltonian, z → maxu∈Rm H(z, u), is well defined and smooth in a neighbourhood of z¯(·), then one necessarily has h(z) = H(z, u ¯(z)) = maxu∈Rm H(z, u) on this neighbourhood under the Legendre-Clebsch condition.
1662
Proceedings of the 20th IFAC World Congress Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
A trajectory x ¯(·) is called strictly C 0 -locally optimal if it realizes a strict local minimum t¯f of the cost tf w.r.t. all trajectories of the system close to x ¯(·) in C 0 ([0 , t¯f ], Rn ) (endowed with the uniform topology) and having the same endpoints as x ¯(·). The following result from Agrachev and Sachkov (2004) is crucial in our optimality analysis. Theorem 7. For a normal regular extremal defined on [0 , t¯f ], in the neighbourhood of which the maximized Hamiltonian is smooth, the absence of conjugate time on (0 , t¯f ] is sufficient for strict C 0 -local optimality. Under the additional strong regularity assumption, the extremal is not locally optimal in L∞ topology on [0 , t], for every t > tc , with tc the first conjugate time.
1627
and by σs a singular arc. We denote by σ1 σ2 an arc σ1 followed by an arc σ2 . Ordinary switching time. It is a time t such that two ˙ bang arcs switch with Φ(t) = 0 and Φ(t) = H01 (z(t)) = 0. According to the maximum principle, near Σ, the extremal ˙ ˙ is of the form σ− σ+ if Φ(t) > 0 and σ+ σ− if Φ(t) < 0. Fold case. It is a time where a bang arc has a contact of ¨ ± (t) := H001 (z(t)) ± H101 (z(t)) order 2 with Σ. Denoting Φ the second derivative of Φ with u(·) ≡ ±1, we have three ¨ ± at the switching time: ¨ ± = 0) depending on Φ cases (if Φ
¨ − < 0. A connection ¨ + > 0 and Φ (1) Hyperbolic case: Φ with a singular extremal is possible at Σs and locally each extremal is of the form σ± σs σ± (by convention each arc of the sequence can be empty). ¨ − > 0. The singular extremal ¨ +Φ (2) Parabolic case: Φ at the switching point is not admissible and every extremal curve is locally bang-bang with at most two switchings, i.e. σ+ σ− σ+ or σ− σ+ σ− . ¨ − > 0. A connection ¨ + < 0 and Φ (3) Elliptic case: Φ with a singular arc is not possible and locally each extremal is bang-bang but with no uniform bound on the number of switchings.
2.5 Singular trajectories in the case of affine systems
For optimality analysis, one restricts our study to a single input affine system: x˙ = F0 (x) + u F1 (x), |u| ≤ 1. Relaxing the control bound, singular trajectories are parameterized by the constrained Hamiltonian system: (10) x˙ = ∂p H, p˙ = −∂x H, 0 = ∂u H = H1 , with H1 (x, p) := p, F1 (x) the Hamiltonian lift of F1 . The singular extremals are not regular and the constraint H1 = 0 has to be differentiated at least twice along an extremal to compute the control. Introducing the Lie brackets of 2.7 Conjugate points in the affine case for singular arcs two vector fields F0 and F1 , computed with the convention F01 (x) := [F0 , F1 ](x) := dF1 (x)·F0 (x)−dF0 (x)·F1 (x), and This section relies on the work of Bonnard and Kupka related to the Poisson bracket of two Hamiltonian lifts H0 (1993). Let z(·) := (x(·), p(·)) be a reference singular and H1 by the rule H01 := {H0 , H1 } = H[F0 ,F1 ] , one gets: extremal defined on [0 , T ]. In the affine case, a Jacobi field J(·) := (δx(·), δp(·)) is a solution of the variational H1 = H01 = H001 + u H101 = 0. equation δ z(t) ˙ = dhs (z(t)) · δz(t), t ∈ [0 , T ], satisfying A singular extremal along which H101 = 0 is called of also dH1 (z(t)) · J(t) = dH01 (z(t)) · J(t) = 0. The Jacobi minimal order and the corresponding control is given by: field is said semi-vertical at time t if δx(t) ∈ RF1 (x(t)). H001 (z) A time tc ∈ (0 , T ] is a conjugate time if and only if there . us (z) := − exists a non trivial Jacobi field semi-vertical at 0 and tc . H101 (z) Plugging such us in H defines a true Hamiltonian, denoted We introduce the following assumptions: hs , whose solutions initiating from H1 = H01 = 0 define the singular extremals of order 1. They are related to the (A1) H101 = 0 along z(·), F0 and F1 are linearly independent along x(·) ¶and x(·) is injective. regular case using the Goh transformation, see Bonnard © and Chyba (2003). The Legendre-Clebsch condition is (A2) K(t) := Span adk F0 · F1 (x(t)) | k = 0, · · · , n − 2 replaced and we have the following additional necessary has codimension one, where ad F0 · F1 = F01 . condition of optimality deduced from the high-order max- (A3) Non-exceptional case: Hs (z(t)) is nonzero. imum principle, see Krener (1977). If the control us (·) is singular and non saturating, i.e. |us | < 1, the following Then, we have the following result from Bonnard and generalized Legendre-Clebsch condition must hold along Kupka (1993): let tc be the first conjugate time. Under assumptions (A1)-(A2)-(A3), the reference singular trajecthe associated singular extremal: tory x(·) is C 0 -locally time minimizing in the hyperbolic 2 ∂ ∂ ∂H = H101 ≥ 0. (11) case and time maximizing in the elliptic case on [0 , tc ). ∂u ∂t2 ∂u Moreover, x(·) is not time optimal on [0 , t] in L∞ topology for every t > tc . 2.6 Generic classification of the bang-bang extremals We consider the affine system x˙ = F0 (x)+u F1 (x), |u| ≤ 1. An important issue is to apply the results from Kupka (1987) to classify the extremal curves near the switching surface. The switching surface is the set Σ : H1 = 0, while the switching function is t → Φ(t) := H1 (z(t)), where z(·) is an extremal curve. Let Σs : H1 = H01 = 0. The singular extremals are entirely contained in Σs . A bangbang extremal z(·) on [0 , T ] is an extremal curve with a finite number of switching times 0 ≤ t1 < · · · < tn ≤ T . We denote by σ+ , σ− the bang arcs for which u(·) ≡ ±1
3. APPLICATION TO THE REDUCED MODEL 3.1 The reduced model when ε = 0 Putting ε = 0 in equation (4’), we get g 1 ρ(h) S v 0= u − cos γ. 2 m v Solving (12) considering γ is small gives 2mg . u ¯= ρ(h) S v 2
1663
(12)
Proceedings of the 20th IFAC World Congress 1628 Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
Plugging u ¯ in equations (1)-(2)-(3) and taking γ as the new control variable gives the following system: x˙ = F0 (x) + u F1 (x), with x := (h, v, m), u := γ, ã Å 2 m g2 ∂ T (h) 1 ρ(h) S v 2 − C D1 − C F0 (x) = D2 m 2 m ρ(h) S v 2 ∂v ∂ , − (Cs (v) T (h)) ∂m ∂ ∂ −g . F1 (x) = v ∂h ∂v The optimal control problem is then: min tf , ˙ = F0 (x(t)) + u(t) F1 (x(t)), x(t) umin ≤ u(t) ≤ umax , t ∈ [0 , tf ] a.e., (P0 ) x(0) = x 0, x(tf ) = xf . We consider a realistic case where the initial state is x0 := (3480, 151.67, 69000) 5 , the final state is xf := (9144, 191, 68100). The bounds in the control are given by the maximal authorized air slope, umax := 0.262 radian, and we complete symmetrically with umin := −umax . See table 1 for the chosen values of the constant parameters. Table 1. Medium-haul aircraft parameters. Data S g C T1 C T2 C T3 C D1 CD2
Value
Unit
122.6 9.81 141040 14909.9 6.997e−10 0.0242 0.0469
m2 m.s−2 N m m−2
Data
Value
Unit
Cs 1 Cs 2 R Θ0 β P0
1.055e−5
kg.s−1 .N−1 m.s−1 J.kg−1 .K−1 K K.m−1 Pa
441.54 287.058 288.15 0.0065 101325
3.2 Singular extremals Introducing D0 (x) := det(F1 (x), F01 (x), F0 (x)), D001 (x) := det(F1 (x), F01 (x), F001 (x)), D101 (x) := det(F1 (x), F01 (x), F101 (x)), the singular control is given in feedback form by D001 (x) us (x) = − , (13) D101 (x) whenever D101 (x) = 0, since, along a singular extremal p, F1 (x) = p, F01 (x) = p, F001 (x) + u F101 (x) = 0 and p(·) never vanishes. Assuming D0 (x) = 0, we can write
F001 (x) − F101 (x) = α F0 (x) + α1 F1 (x) + α01 F01 (x), F001 (x) + F101 (x) = β F0 (x) + β 1 F1 (x) + β 01 F01 (x), which gives two relations at x: D001 − D101 = α D0 and D001 + D101 = β D0 . Assuming D001 − D101 = 0 and D001 + D101 = 0, the classification from section 2.6 is given by α and β. For instance, the hyperbolic case is given by α < 0 and β > 0. Besides, the generalized LegendreClebsch condition (11) becomes D0 D101 ≥ 0. 5 The altitute h is given in meter (m), the true air speed v in meter per second (m.s−1 ) and the mass m in kilogram (kg).
3.3 Shooting function and BC-extremal We use the Bocop software, Bonnans et al. (2012), based on direct methods to determine the structure and to get an initial guess for the shooting method we describe hereinafter. This procedure is justified since the indirect methods are very sensitive with respect to the initial guess. The application of the direct method gives a trajectory of the form σ− σs σ+ . We introduce the true Hamiltonians: h± (x, p) := H0 (x, p) ± H1 (x, p), hs (x, p) := H0 (x, p) + us (x) H1 (x, p). We define the shooting function S : R5n+3 → R5n+3 by: H1 (z1 ) H01 (z1 ) h+ (z(tf )) + p0 S(p0 , tf , t1 , t2 , z1 , z2 ) := (14) , πx (z(tf )) − xf z(t ) − z 1 1 z(t2 ) − z2 where: z(t1 ) := exp(t1 h− )(x0 , p0 ), z(t2 ) := exp((t2 − t1 )hs )(z1 ), z(tf ) := exp((tf − t2 )h+ )(z2 ). The two first equations are junction conditions with the singular extremal, then we have the boundary conditions and finally the matching conditions. The shooting method consists in solving S(y) = 0, with y := (p0 , tf , t1 , t2 , z1 , z2 ). An important point is that to any zero of S is associated a unique BC-extremal of problem (P0 ). We use the HamPath code, Caillau et al. (2011), based on indirect methods to solve the shooting equation (14) and we find a solution y¯ satisfying S(¯ y ) ≈ 6e−9 and given 6 −2 by p¯0 ≈ (2.673e , 0.448, −0.327), t¯f ≈ 644.2, t¯1 = 17.43 and t¯2 = 628.5. The control law u(·) = γ(·) associated to y¯ is given in Fig. 1. We check a posteriori that we are in the hyperbolic case and that the strict generalized LegendreClebsch condition is satisfied along the singular arc. Hence, the singular extremal is time minimizing up to the first conjugate time.
u
t¯1
t
¯ t¯2 tf
Fig. 1. The control law u(·) = γ(·) for the trajectory of the form σ− σs σ+ associated to y¯. 6
We do not give z¯1 nor z¯2 which can be computed by numerical integration.
1664
Proceedings of the 20th IFAC World Congress Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
3.4 Second order conditions of optimality Since D0 D101 > 0 along the singular extremal associated to y¯ and because the associated trajectory is injective (m(·) is strictly decreasing), then assumptions (A1), (A2) and (A3) from section 2.7 are satisfied. To compute the first conjugate time, we can take into account that the dimension of the state is 3. Let x(·) be a reference singular trajectory contained in D101 = 0 defined on [t1 , t2 ]. We define a Jacobi field J(·) along x(·) as a non trivial solution of the variational equation ∂ δ x(t) ˙ = (F0 + us F1 )(x(t)) · δx(t), ∂x satisfying the initial condition J(t1 ) ∈ RF1 (x(t1 )). The first conjugate time is the first time tc > t1 such that det (J(tc ), F0 (x(tc )), F1 (x(tc ))) = 0. Finally, the singular extremal associated to y¯ is time minimizing on [t¯1 , t¯2 ] according to Fig. 2.
1629
defined by ti = ti−1 + ∆t, i = 1, · · · , k − 1 with k ∈ N* , ∆t = (tf − t0 )/k and t0 := 0. Then, the multiple shooting function Sεk (p0 , tf , z1 , · · · , zk−1 ) is given by: π (z(t , t x f k−1 , zk−1 )) − xf hε (z(tf , tk−1 , zk−1 )) + p0 z(t1 , t0 , z0 ) − z1 , Sεk (p0 , tf , z1 , · · · , zk−1 ) := .. . z(tk−1 , tk−2 , zk−2 ) − zk−1 where z0 := (x0 , p0 ), z : (t1 , t0 , z0 ) → exp((t1 − t0 )hε )(z0 ).
We first use multiple shooting method to solve the equa5 (y) = 0, y := (p0 , tf , z1 , · · · , z4 ). Let denote by tion S100 5 y¯100 the solution we find. Then, we define the homotopic function ϕ(y, ε) := Sε5 (y). We use the differential path following method from HamPath code to solve the equation 5 ϕ(y, ε) = 0 starting from the initial point (¯ y100 , 100) to the target ε = 1. The solution from the homotopy at ε = 1 is used to initialize the multiple shooting method and solve the equation S15 (y) = 0. Let denote by y¯15 this solution. Then, we get S15 (¯ y15 ) ≈ 5e−9 .
We compare the state variable γ when ε ∈ {1, 10} in Fig. 3 to emphasize the impact of the singular perturbation.
We can see in Fig. 4 that the slow variable h is very well approximated by the one from the reduced-order problem (P0 ). In Fig. 5, we see that what we call the outer solution of the fast variable in the singular perturbation theory, is very well approximated by the singular optimal control of problem (P0 ). Besides, we get a very small relative gap, around 0.14%, between the final times associated to both solutions of problems (P0 ) and (P1 ).
Λ
t¯1
t
t¯2
Fig. 2. The determinant Λ := det (J, F0 , F1 ), with J(t¯1 ) = F1 (x(t¯1 )), along the singular arc for the trajectory of the form σ− σs σ+ associated to y¯. 4. MODEL WITH SINGULAR PERTURBATION We are interested in this section in the following optimal control problem: min tf , ˙ = Fε (x(t), u(t)), u(t) ∈ R, t ∈ [0 , tf ] a.e., (Pε ) x(t) x(0) = x0 , x(tf ) = xf , where x := (h, v, m, γ), u is the lift coefficient, the initial state is x0 := (3480, 151.67, 69000, 0), the final state is xf := (9144, 191, 68100, 0) and where the dynamics Fε is given by equations (1)-(2)-(3)-(4’), ε ≥ 1 fixed. The values of the constant parameters are given by table 1. 4.1 Multiple shooting and homotopy on ε The application of the PMP gives the maximizing control: pγ uε (x, p) := 2 ε pv v CD2 and plugging uε in H(x, p, u) := p, Fε (x, u) gives the true Hamiltonian hε (x, p) := p, Fε (x, uε (x, p)). We use multiple shooting to deal with numerical instability due to the singular perturbation. We note (ti , zi ) the intermediate discretized times and points. The times ti are fixed and
4.2 Second order conditions of optimality We check now the second order conditions of optimality along the path of zeros. For a fixed ε ∈ [1 , 100] we denote by pε0 the initial costate solution of the associated shooting equation. Then, we compute the three Jacobi fields Ji (·) := (δxi (·), δpi (·)) such that δxi (0) = 0 and where (δp1 (0), δp2 (0), δp3 (0)) is a basis satisfying dhε (x0 , pε0 ) · δpi (0) = 0. The first conjugate time tc > 0 is the first time such that det (δx1 , δx2 , δx3 , Fε ) vanishes. The determinant being very big when ε is close to 1, we prefer to use a SVD decomposition and compute the smallest singular value. ε We denote by σmin (t) the minimal singular value at time t ε and σmax (t) the maximal one. Then, tc > 0 is a conjugate ε time if and only if σmin (tc ) = 0. According to Fig. 6 the BC-extremal along the path of zeros are locally optimal for ε ∈ [1.5 , 100] 7 . When ε = 1, we cannot conclude about the local optimality because of significant round-off errors, see Fig. 7. 5. CONCLUSION Combining theoretical and numerical tools from geometric control, we get a σ− σs σ+ hyperbolic trajectory, solution of the PMP for the reduced-order problem (P0 ). Despite the singular perturbation, we also present BC-extremals 7
1665
ε (·) for ε ∈ {1.5, 2, 5, 10} for readibility. We only represent σmin
Proceedings of the 20th IFAC World Congress 1630 Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
satisfying second order conditions of optimality for problems (Pε ), ε > 0. In conclusion, the numerical comparison between problems (P0 ) and (P1 ) shows that the reducedorder problem is a very good approximation of the original time-to-climb problem.
ε = 1.5
ε=2 ε σmin
ε=5 ε = 10
γ
ε = 10
Normalized time ε Fig. 6. The singular value σmin (·) for ε ∈ {1.5, 2, 5, 10}. ε=1
Normalized time
Fig. 3. The air slope γ(·) for ε = 10 (dashed line) and ε = 1 (solid line).
1 σmin
1 σmax
Normalized time h
ε ε Fig. 7. The singular values σmin (·) and σmax (·) for ε = 1.
Normalized time
Fig. 4. The altitude h(·) for ε = 1 (blue dashed line) and ε = 0 (red solid line).
γ
Normalized time
Fig. 5. The air slope γ(·) for ε = 1 (blue dashed line) and ε = 0 (red solid line). When ε = 0, γ(·) is the control law of problem (P0 ). REFERENCES A. A. Agrachev & Y. L. Sachkov, Control theory from the geometric viewpoint, vol 87 of Encyclopaedia of Mathematical Sciences, Springer-Verlag, Berlin (2004), xiv+412. M. D. Ardema, Solution of the Minimum Time-to-Climb Problem by Matched Asymptotic Expansions, AIAA Journal of Guidance, Control, and Dynamics, 14 (1976), no. 7, 843–850.
F. J. Bonnans, P. Martinon & V. Gr´elard. Bocop - A collection of examples, Technical report, INRIA, 2012. RR-8053. B. Bonnard & M. Chyba, Singular trajectories and their role in control theory, vol 40 of Mathematics & Applications, Springer-Verlag, Berlin (2003), xvi+357. B. Bonnard & I. Kupka, Th´eorie des singularit´es de l’application entr´ee/sortie et optimalit´e des trajectoires singuli`eres dans le probl`eme du temps minimal, Forum Math., 5 (1993), no. 2, 111–159. J.-B. Caillau, O. Cots & J. Gergaud, Differential continuation for regular optimal control problems, Optimization Methods and Software, 27 (2011), no. 2, 177–196. A. J. Krener, The high order maximal principle and its application to singular extremals, SIAM J. Control Optim., 15 (1977), no. 2, 256–293. I. Kupka. Geometric theory of extremals in optimal control problems. i. the fold and maxwell case, Trans. Amer. Math. Soc., 299 (1987), no. 1, 225–243. N. Mo¨ıss´eev, Probl`emes math´ematiques d’analyse des syst`emes, Translated from the Russian D. Embarek, edited by Mir, Tomash Publishers, Moscow, 1985. N. Nguyen, Singular arc time-optimal climb trajectory of aircraft in a two-dimensional wind field, In AIAA Guidance, Navigation and Control Conference and Exhibit, Guidance, Navigation, and Control and Co-located Conferences, august 2006. L. S. Pontryagin, V. G. Boltyanski˘ı, R. V. Gamkrelidze & E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Translated from the Russian by K. N. Trirogoff, edited by L. W. Neustadt, Interscience Publishers John Wiley & Sons, Inc., New York-London, 1962.
1666
ScienceDirect
IFAC PapersOnLine 50-1 (2017) 1625–1630
Time-optimal aircraft trajectories in Time-optimal aircraft trajectories in Time-optimal aircraft trajectories in Time-optimal aircraft trajectories in climbing phase and singular perturbations climbing phase and singular perturbations climbing phase and singular perturbations climbing phase and singular perturbations
Olivier Cots ∗∗ Joseph Gergaud ∗∗ Damien Goubinat ∗∗ ∗∗ Olivier Cots ∗∗ Joseph Gergaud ∗∗ Damien Goubinat ∗∗ Olivier Olivier Cots Cots Joseph Joseph Gergaud Gergaud Damien Damien Goubinat Goubinat ∗∗ ∗ ∗ Toulouse Univ., INP-ENSEEIHT-IRIT, UMR CNRS 5505, 2 rue ∗ Toulouse Univ., INP-ENSEEIHT-IRIT, UMR CNRS 5505, 2 rue ∗ Toulouse31071 Univ., INP-ENSEEIHT-IRIT, CNRS Camichel, France (e-mail: UMR [email protected] Toulouse31071 Univ.,Toulouse, INP-ENSEEIHT-IRIT, UMR CNRS 5505, 5505, 2 2 rue rue& Camichel, Toulouse, France (e-mail: [email protected] & Camichel, 31071 Toulouse, France (e-mail: [email protected] [email protected]). Camichel, 31071 Toulouse, France (e-mail: [email protected] & & [email protected]). ∗∗ [email protected]). SA, 105 av du General Eisenhower, B.P. 1147, ∗∗ Thales Avionics [email protected]). ∗∗ Thales Avionics SA, 105 av du General Eisenhower, B.P. 1147, ∗∗ Thales Avionics SA, 105 du Eisenhower, 31047 Toulouse Cedex, [email protected]) Thales Avionics SA,France 105 av av(e-mail: du General General Eisenhower, B.P. B.P. 1147, 1147, 31047 Toulouse Cedex, France (e-mail: [email protected]) 31047 Toulouse Cedex, France (e-mail: [email protected]) 31047 Toulouse Cedex, France (e-mail: [email protected]) Abstract: This article deals with the singularly perturbed aircraft minimum time-to-climb Abstract: This article deals with the singularly perturbed aircraft minimum time-to-climb Abstract: This article deals the perturbed aircraft minimum time-to-climb problem. First, introduce reduced-order problem with affine dynamics with respect to the Abstract: Thiswe article dealsaa with with the singularly singularly perturbed aircraft minimum time-to-climb problem. First, we introduce reduced-order problem with affine dynamics with respect towith the problem. First, we introduce a reduced-order problem with affine dynamics with respect the control and analyze it with the tools from geometric control: maximum principle combined problem. First, we introduce a reduced-order problem with affine dynamics withcombined respect to towith the control and analyze it with the tools from geometric control: maximum principle control order and analyze analyze it with with the tools toolsThen, from geometric geometric control: maximum principle combined combined with second optimality conditions. we compute candidates as minimizers using multiple control and it the from control: maximum principle with second order optimality conditions. Then, we compute candidates as minimizers using multiple second order Then, we compute as multiple shooting andoptimality homotopyconditions. methods for both that we briefly compareusing numerically. second order optimality conditions. Then, we problems, compute candidates candidates as minimizers minimizers using multiple shooting and homotopy methods for both problems, that we briefly compare numerically. shooting and homotopy methods for both problems, that we briefly compare numerically. Particular attention is paid to the singular extremals and we discuss about their local optimality. shooting and homotopy methods for both problems, that we briefly compare numerically. Particular attention is paid to the singular extremals and we discuss about their local optimality. Particular attention is to extremals and we about local optimality. Particular attention is paid paidFederation to the the singular singular extremals and Hosting we discuss discuss about their their local optimality. © 2017, IFAC (International of Automatic Control) by Elsevier Ltd. All rights reserved. Keywords: Aircraft trajectory, geometric optimal control, singular perturbation. Keywords: Aircraft trajectory, geometric optimal control, singular perturbation. Keywords: Keywords: Aircraft Aircraft trajectory, trajectory, geometric geometric optimal optimal control, control, singular singular perturbation. perturbation. g Å ã βR 1. INTRODUCTION g Å Θ(h) ã ã βR g 1. INTRODUCTION INTRODUCTION Å g , 1. Θ(h) := ã Å P (h) P βR 0 Θ(h) βR 1. INTRODUCTION := ,, P (h) P Θ(h) Θ 0 0 := P (h) P 0 Θ0 , P (h) := P0 Θ A flight is composed of several phases which are take-off, 0 := Θ(h) Θ − βh, A flight is composed of several phases which are take-off, Θ 0 0 A flight is composed several phases which are take-off, climb, cruise, descent,of approach and landing. := Θ(h) Θ − βh, A flight is composed of several phases which are take-off, 0 Θ(h) − climb, cruise, cruise, descent, approach approach and landing. landing. := Θ Θ(h) := ΘP00 (h) − βh, βh, climb, climb, cruise, descent, descent, approach and and landing. (h) . ρ(h) := P P (h) := . ρ(h) P (h) RΘ(h) := ρ(h) . RΘ(h) . ρ(h) := RΘ(h) RΘ(h) constants: g is the graviThe remaining data are positive The remaining data are wing positive constants: is the the graviThe remaining data are positive constants: ggg is gravitational constant, S the area, R the specific constant The remaining data are wing positive constants: is the gravitational constant, S the area, R the specific constant tational S the wing area, R specific of air, β constant, the thermical gradient and Pthe pressure 0, Θ 0 , the constant tational constant, S the wing area, R the specific constant of air, air, β temperature the thermical thermical gradient and P P0 , Θ0 ,, the the pressure of the and and theβ atgradient the sea level. of air, β temperature the thermical gradient and P00 ,, Θ Θ00 , the pressure pressure the at the sea level. In this article, we are interested in the time optimal control and and the temperature at the sea level. In this article, we are interested in the time optimal control and the temperature at the sea level. In this article, are in time control of an aircraft during its climbing This phase is The underlying optimal control problem consists of miniIn article, we we are interested interested in the thephase. time optimal optimal control The underlying optimal control problem consists of cruise miniof this an aircraft aircraft during itsdynamics climbing phase. This phase is The underlying optimal of miniof an during its climbing phase. This phase is mizing the transfer time,control i.e. theproblem time to consists reach the determined by its own governed by the fourunderlying optimal control problem consists of cruise miniof an aircraft during itsdynamics climbinggoverned phase. This phase is The mizing the transfer time, i.e. the time to reach the determined by its own by the fourmizing the transfer time, i.e. the time to reach the cruise determined by its own dynamics governed by the fourphase, with fixed initial and final states. It is well known dimensional system mizing the transfer time, i.e. the time to reach the cruise determined by its own dynamics governed by the fourphase, with fixed initial and final states. is well known dimensional system system phase, with fixed and final It is known dimensional that the dynamics (1)–(4) slowIt mass m) phase, with fixed initial initial and contains final states. states. It (the is well well known dimensional system dh that the dynamics (1)–(4) contains slow (the mass m) 1 that the dynamics (1)–(4) contains slow (the mass m) dh = v sin γ (1) and fast (the air slope γ) variables , see Ardema (1976) the (the dynamics (1)–(4) contains11 , slow (the mass m) dh = vv sin sin γ γ (1) that and fast air slope γ) variables see Ardema (1976) dh dt = (1) 1 and fast (the air slope γ) variables , see Ardema (1976) Nguyen (2006). Thisγ)time scale separation is normally dt = v sin γ (1) and fast (the air slope variables , see Ardema (1976) 2 dt and Nguyen (2006). This time scale separation is normally 1 ρ(h)Sv T (h) dv dt = T (h) − 1 ρ(h)Sv 22 C + C u2 − g sin γ (2) treated and This time normally by a(2006). singular perturbation analysis byis dv and Nguyen Nguyen This time scale scale separation separation is replacing normally 1 2 2 (h) dv by aa(2006). singular perturbation analysis by replacing − 121 ρ(h)Sv =T CD +C CD u22 − − gg sin sin γ γ (2) (2) treated ρ(h)Sv Tm (h) − dv dt = m 2 C treated by singular perturbation analysis by replacing D1 + D2 u the fast state equation with the following: treated by a singular perturbation analysis by replacing D1 + C D2 u 2 m dt m − = C − g sin γ (2) the equation with D1 D2 dm dt m the fast fast state state equation with the thegfollowing: following: ms (v) T22(h) m dt = −C m with the following: dγequation 1 ρ(h)Sv dm (3) the fast state dm dγ 1 ρ(h)Sv (4’) u − gg cos γ, ε > 0. ε = = −C (v) T (h) (3) dm dt = −Css (v) T (h) 1 dγ ρ(h)Sv (3) cos γ, γ, εε > > 0. 0. (4’) − εε dγ 12 ρ(h)Sv vg cos m u dt = dt = −C (3) s (v) T (h) g = u − (4’) dt 1 ρ(h)Sv dγ v m γ, ε > 0. (4’) ε dt dt 2 2 dt = 2 time-to-climb m u − vv cos g 1 ρ(h)Sv dγ cos γ, (4) u − = Let (P ) denote the problem with the addi2 m dt ε dγ 2 cos γ, γ, (4) Let u− − vgg cos = 11 ρ(h)Sv dγ 2 the time-to-climb problem with the addi2 ρ(h)Sv m u dt = Let (P ) denote ε (4) 2 ) denote the time-to-climb problem with the addi(P ε > 0. From the control theory, m u − vv cos γ, dt = 22 m (4) tional Let (Pεεartificial ) denote parameter the time-to-climb problem with the addidt artificial εε > 0. From the theory, where v x := (h, v, m, γ) is composed of tional m variable dt the2 state tional artificial parameter > 0. From the control theory, for a fixed valueparameter of ε > 0, the candidates as control minimizers are tional artificial parameter ε > 0. From the control theory, := where the state variable x (h, v, m, γ) is composed of for aa fixed value of εε > 0, the candidates as minimizers are := where the state (h, v, γ) is composed of the altitude, thevariable true airx speed, them, weight and the air for fixed value of > 0, the candidates as minimizers are := where the state variable x (h, v, m, γ) is composed of selected among a set of extremals, solution of a Hamiltofor a fixed value aofset ε >of0,extremals, the candidates as minimizers are the altitude, the true trueInair airthis speed, thethe weight and the the air air selected among solution of a Hamiltothe altitude, the speed, the weight and slope of the aircraft. model, lift coefficient is selected among a set of extremals, solution of a Hamiltothe altitude, the true air speed, the weight and the air nian system given by the Pontryagin Maximum Principle selected among a set of extremals, solution of a Hamiltoslope of the aircraft. In this model, the lift coefficient is system given by the Pontryagin Maximum Principle slope of In lift coefficient is taken as the the aircraft. control variable u. The the thrust (h) and the nian system given Pontryagin Maximum Principle slope of aircraft. In this this model, model, the liftT is nian (PMP), see Pontryagin et al. (1962). The application of the nian system given by by the the Pontryagin Maximum Principle taken as the the control variable u. The thrust Tcoefficient (h) and and the (PMP), see Pontryagin et al. (1962). The application of the taken as the control variable thrust T (h) the fuel flow Cs (v) are given by u. theThe Base of Aircraft DAta (PMP), see Pontryagin et al. (1962). The application of the taken as the control variable u. The thrust T (h) and the PMP leads to define a Boundary Value Problem (BV ε) (PMP), see Pontryagin et al. (1962). The application of P the fuel flow C (v) are given by the Base of Aircraft DAta s (v) are given by the Base of Aircraft DAta PMP leads to define a Boundary Value Problem (BV P fuel flow C ε) (BADA) model: s (v) are given by the Base of Aircraft DAta PMP leads to define a Boundary Value Problem (BV P fuel flow C which can be solved using shooting methods combined s PMP leads to define ausing Boundary Valuemethods Problemcombined (BV Pεε )) (BADA) model: which can be solved shooting (BADA) model: which can methods combined (BADA) model: with numerical tests tousing check shooting second-order optimality conh which can be be solved solved using shooting methods combined with numerical tests to check second-order optimality con+ h22 CT3 ), T (h) := CT1 (1 − h with numerical tests to check second-order optimality conh ditions. 2 with numerical tests to check second-order optimality con:= (1 − + h C ), T (h) C h CT2 + h2 CTT3 ), T ditions. T 3 ), := C ditions. (1 − −C + h C T (h) (h) := CTT111 (1 T T3 2 v C ditions. T When dealing with singular perturbation, one classical Cs (v) := Cs1 (1 + CvvT22 ), When dealing with one classical := C + C Cvs2 ), When dealing with singular perturbation, one s (v) s1 (1 approach consists in singular applying perturbation, the PMP to problem (Pε ) := (v) C (1 + ), C When dealing with singular perturbation, one classical classical s s ), Cs (v) := Cs11 (1 + C approach consists in applying the PMP to problem (P s ε) Cs22 approach consists in applying the PMP to problem (P ε) C approach consists in applying the PMP to problem (Pthe s2 CDi (drag coefficients) where the constants CTi , Csi , with 1 The altitude ε) h and the true air speed v are fast compare to , C , with C (drag coefficients) where the constants C 1 T s D i i i The altitude h and the true air speed v are fast compare to the , C , with C (drag coefficients) where the constants C 1 are specific to the aircraft. The International Standard T s D i i i mass but slow compare to the air slope. In this article, we consider The altitude altitude h h and and the the true true air air speed speed v v are are fast fast compare compare to to the the withInternational CDi (drag coefficients) where the constants CTi , Csi ,The 1 The are specific to the aircraft. to the air slope. In this article, we consider are specific the aircraft. The International Standard Atmospheric model of theStandard pressure mass h andbut v asslow slowcompare variables. mass but slow compare to the air slope. In this article, we consider are specific to to the provides aircraft. the Theexpression International Standard mass but slow compare to the air slope. In this article, we consider Atmospheric model provides the expression of the pressure h and v as slow variables. 2 Atmospheric model of pressure P (h), the temperature Θ(h) the andexpression the air density See section 4 for the exact definition of (Pε ), ε > 0. v variables. Atmospheric model provides provides the expression of the theρ(h): pressure 22hh and and v as as slow slow variables. P (h), the temperature Θ(h) and the air density ρ(h): See section 4 for the exact definition of (Pε ), ε > 0. P (h), the temperature Θ(h) and the air density ρ(h): See section 4 for 2 P (h), the temperature Θ(h) and the air density ρ(h): See section 4 for the the exact exact definition definition of of (P (Pεε ), ), εε > > 0. 0. Copyright 1661Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © © 2017 2017, IFAC IFAC (International Federation of Automatic Control) Copyright © 2017 IFAC 1661 Copyright © 2017 1661 Peer review responsibility of International Federation of Automatic Copyright © under 2017 IFAC IFAC 1661Control. 10.1016/j.ifacol.2017.08.327
Proceedings of the 20th IFAC World Congress 1626 Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
and then solve the associated singular boundary value problem with matching asymptotic expansion techniques, see Ardema (1976) and Mo¨ıss´eev (1985). In section 3, to approximate the solutions of problem (Pε ), ε > 0, we consider a reduced-order problem denoted 3 (P0 ), where ε is put to zero and where γ is taken as the control variable. This technique is also applied in Nguyen (2006) and is the first step of a second approach presented in Mo¨ıss´eev (1985). The problem (P0 ) has one state variable less than (Pε ), ε > 0, and can be tackled with the tools from geometric optimal control. This allow first to study the classification of the extremals, solution of the PMP, and then to check, on the singular arcs, second-order conditions of optimality.
with p0 ≤ 0 and (¯ p(·), p0 ) = (0, 0). Definition 1. We call extremal a triplet (x(·), p(·), u(·)) where (x(·), p(·)) is an integral curve of H satisfying (6). It is called a BC-extremal if it satisfies also (7) and (8).
We complete this analysis showing numerically in section 4 that problem (P0 ) from section 3 is a good approximation of problem (Pε ) when ε = 1. We use the HamPath software, Caillau et al. (2011), combining multiple shooting techniques with homotopy methods to solve the family of boundary value problems (BV Pε ) for ε ∈ [1 , 100], starting from the simpler problem when ε = 100.
We have the following characterization of singular controls which allows a practical computation. Proposition 3. The control u(·) with its associated trajectory x(·) are singular on [0 , T ] if and only if there exists a non zero adjoint p(·) such that (x(·), p(·), u(·)) is solution a.e. on [0 , T ] of (9) x˙ = ∂p H, p˙ = −∂x H, 0 = ∂u H.
2. GEOMETRIC OPTIMAL CONTROL
If u(·) is singular then its associated adjoint satisfies for each 0 < t ≤ T , p(t) ⊥ im dEx0 ,t (u(·)). Definition 4. A singular extremal is a triple (x(·), p(·), u(·)) solution of (9). It is called:
2.1 The time optimal control problem One considers a time optimal control problem given by the following data. A control domain U ⊂ Rm . A smooth control system: x(t) ˙ = F (x(t), u(t)), (5) with x(t) ∈ Rn , u(t) ∈ U , m ≤ n. An initial state x0 ∈ Rn and a final state xf ∈ Rn . The optimal control problem can be written tf | Ex0 ,tf (u(·)) = xf , (Ptmin ) min tf >0, u(·)∈Utf
where Utf := L∞ ([0 , tf ], U ) is the set of admissible controls 4 and where Ex0 ,tf is the end-point mapping defined by: Ex0 ,t : Ut → Rn , Ex0 ,t (u(·)) := x(t, x0 , u(·)), where t → x(t, x0 , u(·)) is the trajectory solution of (5) with control u(·) and such that x(0, x0 , u(·)) = x0 . 2.2 Pontryagin Maximum Principle Let define the pseudo-Hamiltonian: H : Rn × (Rn )* × U −→ R (x, p, u) −→ H(x, p, u) := p, F (x, u) . By virtue of the maximum principle, see Pontryagin et al. (1962), if (¯ u(·), t¯f ) is solution of (Ptmin ) then the associated trajectory x ¯(·) is the projection of an absolutely continuous integral curve (¯ x(·), p¯(·)) of H := (∂p H, −∂x H) such that the following maximization condition holds for almost every t ∈ [0 , t¯f ]: x(t), p¯(t), u). (6) H(¯ x(t), p¯(t), u ¯(t)) = max H(¯ u∈U
The following boundary conditions must be satisfied: (7) x ¯(t¯f ) = xf , 0 ¯ ¯ ¯ H(¯ x(tf ), p¯(tf ), u ¯(tf )) = −p , (8) 3
See section 3 for the exact definition of (P0 ). The set of admissible controls is the set of L∞ -mappings on [0 , tf ] taking their values in U such that the associated trajectory x(·) is globally defined on [0 , tf ]. 4
2.3 Preliminaries about singular trajectories A complete presentation of singular extremals – which play a crucial role in our analysis – can be found in Bonnard and Chyba (2003). In the following definitions and proposition we assume U = Rm , i.e. there is no control constraint. Definition 2. A control u(·) ∈ UT is called singular on [0 , T ] if the derivative of Ex0 ,T at u(·) is not surjective.
2 H is negative definite (strict Legendre• Regular if ∂uu Clebsch condition). • Strongly normal if im dEx(t1 ),t2 −t1 (u(·)|[t1 ,t2 ] ) is of corank one for each 0 < t1 < t2 ≤ T . • Exceptional if H = 0 and normal if H = 0.
2.4 Singular trajectories in the regular case In the regular case, using the implicit function theorem, we can solve locally the equation ∂u H = 0 and compute the smooth singular control as a function u ¯(z), z := (x, p). Plugging such u ¯ in H defines a true Hamiltonian h(z) := H(z, u ¯(z)). One can define the exponential mapping expx0 : (t, p0 ) → πx (exp(th)(x0 , p0 )), where h := (∂p h, −∂x h), x0 is fixed, πx is the x−projection (x, p) → x, and exp(th)(x0 , p0 ) is the extremal z(·) at time t solution of z(s) ˙ = h(z(s)), s ∈ [0 , t], z(0) = (x0 , p0 ). The domain of expx0 is R+ × P where P := p0 ∈ (Rn )* | h(x0 , p0 ) = −p0 . This leads to the following definition. Definition 5. Let z(·) := (x(·), p(·)) be a reference extremal solution of h. The time tc is said to be geometrically conjugate if (tc , p(0)) is a critical point of expx(0) . We have the following standard test: Proposition 6. The time tc is geometrically conjugate if and only if there exists a Jacobi field J(·) := (δx(·), δp(·)) solution of the variational equation δ z(t) ˙ = dh(z(t))·δz(t), vertical at 0 and tc , i.e. δx(0) = δx(tc ) = 0, and such that δx(·) ≡ 0 on [0 , tc ]. Let z¯(·) be a reference extremal. If the maximized Hamiltonian, z → maxu∈Rm H(z, u), is well defined and smooth in a neighbourhood of z¯(·), then one necessarily has h(z) = H(z, u ¯(z)) = maxu∈Rm H(z, u) on this neighbourhood under the Legendre-Clebsch condition.
1662
Proceedings of the 20th IFAC World Congress Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
A trajectory x ¯(·) is called strictly C 0 -locally optimal if it realizes a strict local minimum t¯f of the cost tf w.r.t. all trajectories of the system close to x ¯(·) in C 0 ([0 , t¯f ], Rn ) (endowed with the uniform topology) and having the same endpoints as x ¯(·). The following result from Agrachev and Sachkov (2004) is crucial in our optimality analysis. Theorem 7. For a normal regular extremal defined on [0 , t¯f ], in the neighbourhood of which the maximized Hamiltonian is smooth, the absence of conjugate time on (0 , t¯f ] is sufficient for strict C 0 -local optimality. Under the additional strong regularity assumption, the extremal is not locally optimal in L∞ topology on [0 , t], for every t > tc , with tc the first conjugate time.
1627
and by σs a singular arc. We denote by σ1 σ2 an arc σ1 followed by an arc σ2 . Ordinary switching time. It is a time t such that two ˙ bang arcs switch with Φ(t) = 0 and Φ(t) = H01 (z(t)) = 0. According to the maximum principle, near Σ, the extremal ˙ ˙ is of the form σ− σ+ if Φ(t) > 0 and σ+ σ− if Φ(t) < 0. Fold case. It is a time where a bang arc has a contact of ¨ ± (t) := H001 (z(t)) ± H101 (z(t)) order 2 with Σ. Denoting Φ the second derivative of Φ with u(·) ≡ ±1, we have three ¨ ± at the switching time: ¨ ± = 0) depending on Φ cases (if Φ
¨ − < 0. A connection ¨ + > 0 and Φ (1) Hyperbolic case: Φ with a singular extremal is possible at Σs and locally each extremal is of the form σ± σs σ± (by convention each arc of the sequence can be empty). ¨ − > 0. The singular extremal ¨ +Φ (2) Parabolic case: Φ at the switching point is not admissible and every extremal curve is locally bang-bang with at most two switchings, i.e. σ+ σ− σ+ or σ− σ+ σ− . ¨ − > 0. A connection ¨ + < 0 and Φ (3) Elliptic case: Φ with a singular arc is not possible and locally each extremal is bang-bang but with no uniform bound on the number of switchings.
2.5 Singular trajectories in the case of affine systems
For optimality analysis, one restricts our study to a single input affine system: x˙ = F0 (x) + u F1 (x), |u| ≤ 1. Relaxing the control bound, singular trajectories are parameterized by the constrained Hamiltonian system: (10) x˙ = ∂p H, p˙ = −∂x H, 0 = ∂u H = H1 , with H1 (x, p) := p, F1 (x) the Hamiltonian lift of F1 . The singular extremals are not regular and the constraint H1 = 0 has to be differentiated at least twice along an extremal to compute the control. Introducing the Lie brackets of 2.7 Conjugate points in the affine case for singular arcs two vector fields F0 and F1 , computed with the convention F01 (x) := [F0 , F1 ](x) := dF1 (x)·F0 (x)−dF0 (x)·F1 (x), and This section relies on the work of Bonnard and Kupka related to the Poisson bracket of two Hamiltonian lifts H0 (1993). Let z(·) := (x(·), p(·)) be a reference singular and H1 by the rule H01 := {H0 , H1 } = H[F0 ,F1 ] , one gets: extremal defined on [0 , T ]. In the affine case, a Jacobi field J(·) := (δx(·), δp(·)) is a solution of the variational H1 = H01 = H001 + u H101 = 0. equation δ z(t) ˙ = dhs (z(t)) · δz(t), t ∈ [0 , T ], satisfying A singular extremal along which H101 = 0 is called of also dH1 (z(t)) · J(t) = dH01 (z(t)) · J(t) = 0. The Jacobi minimal order and the corresponding control is given by: field is said semi-vertical at time t if δx(t) ∈ RF1 (x(t)). H001 (z) A time tc ∈ (0 , T ] is a conjugate time if and only if there . us (z) := − exists a non trivial Jacobi field semi-vertical at 0 and tc . H101 (z) Plugging such us in H defines a true Hamiltonian, denoted We introduce the following assumptions: hs , whose solutions initiating from H1 = H01 = 0 define the singular extremals of order 1. They are related to the (A1) H101 = 0 along z(·), F0 and F1 are linearly independent along x(·) ¶and x(·) is injective. regular case using the Goh transformation, see Bonnard © and Chyba (2003). The Legendre-Clebsch condition is (A2) K(t) := Span adk F0 · F1 (x(t)) | k = 0, · · · , n − 2 replaced and we have the following additional necessary has codimension one, where ad F0 · F1 = F01 . condition of optimality deduced from the high-order max- (A3) Non-exceptional case: Hs (z(t)) is nonzero. imum principle, see Krener (1977). If the control us (·) is singular and non saturating, i.e. |us | < 1, the following Then, we have the following result from Bonnard and generalized Legendre-Clebsch condition must hold along Kupka (1993): let tc be the first conjugate time. Under assumptions (A1)-(A2)-(A3), the reference singular trajecthe associated singular extremal: tory x(·) is C 0 -locally time minimizing in the hyperbolic 2 ∂ ∂ ∂H = H101 ≥ 0. (11) case and time maximizing in the elliptic case on [0 , tc ). ∂u ∂t2 ∂u Moreover, x(·) is not time optimal on [0 , t] in L∞ topology for every t > tc . 2.6 Generic classification of the bang-bang extremals We consider the affine system x˙ = F0 (x)+u F1 (x), |u| ≤ 1. An important issue is to apply the results from Kupka (1987) to classify the extremal curves near the switching surface. The switching surface is the set Σ : H1 = 0, while the switching function is t → Φ(t) := H1 (z(t)), where z(·) is an extremal curve. Let Σs : H1 = H01 = 0. The singular extremals are entirely contained in Σs . A bangbang extremal z(·) on [0 , T ] is an extremal curve with a finite number of switching times 0 ≤ t1 < · · · < tn ≤ T . We denote by σ+ , σ− the bang arcs for which u(·) ≡ ±1
3. APPLICATION TO THE REDUCED MODEL 3.1 The reduced model when ε = 0 Putting ε = 0 in equation (4’), we get g 1 ρ(h) S v 0= u − cos γ. 2 m v Solving (12) considering γ is small gives 2mg . u ¯= ρ(h) S v 2
1663
(12)
Proceedings of the 20th IFAC World Congress 1628 Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
Plugging u ¯ in equations (1)-(2)-(3) and taking γ as the new control variable gives the following system: x˙ = F0 (x) + u F1 (x), with x := (h, v, m), u := γ, ã Å 2 m g2 ∂ T (h) 1 ρ(h) S v 2 − C D1 − C F0 (x) = D2 m 2 m ρ(h) S v 2 ∂v ∂ , − (Cs (v) T (h)) ∂m ∂ ∂ −g . F1 (x) = v ∂h ∂v The optimal control problem is then: min tf , ˙ = F0 (x(t)) + u(t) F1 (x(t)), x(t) umin ≤ u(t) ≤ umax , t ∈ [0 , tf ] a.e., (P0 ) x(0) = x 0, x(tf ) = xf . We consider a realistic case where the initial state is x0 := (3480, 151.67, 69000) 5 , the final state is xf := (9144, 191, 68100). The bounds in the control are given by the maximal authorized air slope, umax := 0.262 radian, and we complete symmetrically with umin := −umax . See table 1 for the chosen values of the constant parameters. Table 1. Medium-haul aircraft parameters. Data S g C T1 C T2 C T3 C D1 CD2
Value
Unit
122.6 9.81 141040 14909.9 6.997e−10 0.0242 0.0469
m2 m.s−2 N m m−2
Data
Value
Unit
Cs 1 Cs 2 R Θ0 β P0
1.055e−5
kg.s−1 .N−1 m.s−1 J.kg−1 .K−1 K K.m−1 Pa
441.54 287.058 288.15 0.0065 101325
3.2 Singular extremals Introducing D0 (x) := det(F1 (x), F01 (x), F0 (x)), D001 (x) := det(F1 (x), F01 (x), F001 (x)), D101 (x) := det(F1 (x), F01 (x), F101 (x)), the singular control is given in feedback form by D001 (x) us (x) = − , (13) D101 (x) whenever D101 (x) = 0, since, along a singular extremal p, F1 (x) = p, F01 (x) = p, F001 (x) + u F101 (x) = 0 and p(·) never vanishes. Assuming D0 (x) = 0, we can write
F001 (x) − F101 (x) = α F0 (x) + α1 F1 (x) + α01 F01 (x), F001 (x) + F101 (x) = β F0 (x) + β 1 F1 (x) + β 01 F01 (x), which gives two relations at x: D001 − D101 = α D0 and D001 + D101 = β D0 . Assuming D001 − D101 = 0 and D001 + D101 = 0, the classification from section 2.6 is given by α and β. For instance, the hyperbolic case is given by α < 0 and β > 0. Besides, the generalized LegendreClebsch condition (11) becomes D0 D101 ≥ 0. 5 The altitute h is given in meter (m), the true air speed v in meter per second (m.s−1 ) and the mass m in kilogram (kg).
3.3 Shooting function and BC-extremal We use the Bocop software, Bonnans et al. (2012), based on direct methods to determine the structure and to get an initial guess for the shooting method we describe hereinafter. This procedure is justified since the indirect methods are very sensitive with respect to the initial guess. The application of the direct method gives a trajectory of the form σ− σs σ+ . We introduce the true Hamiltonians: h± (x, p) := H0 (x, p) ± H1 (x, p), hs (x, p) := H0 (x, p) + us (x) H1 (x, p). We define the shooting function S : R5n+3 → R5n+3 by: H1 (z1 ) H01 (z1 ) h+ (z(tf )) + p0 S(p0 , tf , t1 , t2 , z1 , z2 ) := (14) , πx (z(tf )) − xf z(t ) − z 1 1 z(t2 ) − z2 where: z(t1 ) := exp(t1 h− )(x0 , p0 ), z(t2 ) := exp((t2 − t1 )hs )(z1 ), z(tf ) := exp((tf − t2 )h+ )(z2 ). The two first equations are junction conditions with the singular extremal, then we have the boundary conditions and finally the matching conditions. The shooting method consists in solving S(y) = 0, with y := (p0 , tf , t1 , t2 , z1 , z2 ). An important point is that to any zero of S is associated a unique BC-extremal of problem (P0 ). We use the HamPath code, Caillau et al. (2011), based on indirect methods to solve the shooting equation (14) and we find a solution y¯ satisfying S(¯ y ) ≈ 6e−9 and given 6 −2 by p¯0 ≈ (2.673e , 0.448, −0.327), t¯f ≈ 644.2, t¯1 = 17.43 and t¯2 = 628.5. The control law u(·) = γ(·) associated to y¯ is given in Fig. 1. We check a posteriori that we are in the hyperbolic case and that the strict generalized LegendreClebsch condition is satisfied along the singular arc. Hence, the singular extremal is time minimizing up to the first conjugate time.
u
t¯1
t
¯ t¯2 tf
Fig. 1. The control law u(·) = γ(·) for the trajectory of the form σ− σs σ+ associated to y¯. 6
We do not give z¯1 nor z¯2 which can be computed by numerical integration.
1664
Proceedings of the 20th IFAC World Congress Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
3.4 Second order conditions of optimality Since D0 D101 > 0 along the singular extremal associated to y¯ and because the associated trajectory is injective (m(·) is strictly decreasing), then assumptions (A1), (A2) and (A3) from section 2.7 are satisfied. To compute the first conjugate time, we can take into account that the dimension of the state is 3. Let x(·) be a reference singular trajectory contained in D101 = 0 defined on [t1 , t2 ]. We define a Jacobi field J(·) along x(·) as a non trivial solution of the variational equation ∂ δ x(t) ˙ = (F0 + us F1 )(x(t)) · δx(t), ∂x satisfying the initial condition J(t1 ) ∈ RF1 (x(t1 )). The first conjugate time is the first time tc > t1 such that det (J(tc ), F0 (x(tc )), F1 (x(tc ))) = 0. Finally, the singular extremal associated to y¯ is time minimizing on [t¯1 , t¯2 ] according to Fig. 2.
1629
defined by ti = ti−1 + ∆t, i = 1, · · · , k − 1 with k ∈ N* , ∆t = (tf − t0 )/k and t0 := 0. Then, the multiple shooting function Sεk (p0 , tf , z1 , · · · , zk−1 ) is given by: π (z(t , t x f k−1 , zk−1 )) − xf hε (z(tf , tk−1 , zk−1 )) + p0 z(t1 , t0 , z0 ) − z1 , Sεk (p0 , tf , z1 , · · · , zk−1 ) := .. . z(tk−1 , tk−2 , zk−2 ) − zk−1 where z0 := (x0 , p0 ), z : (t1 , t0 , z0 ) → exp((t1 − t0 )hε )(z0 ).
We first use multiple shooting method to solve the equa5 (y) = 0, y := (p0 , tf , z1 , · · · , z4 ). Let denote by tion S100 5 y¯100 the solution we find. Then, we define the homotopic function ϕ(y, ε) := Sε5 (y). We use the differential path following method from HamPath code to solve the equation 5 ϕ(y, ε) = 0 starting from the initial point (¯ y100 , 100) to the target ε = 1. The solution from the homotopy at ε = 1 is used to initialize the multiple shooting method and solve the equation S15 (y) = 0. Let denote by y¯15 this solution. Then, we get S15 (¯ y15 ) ≈ 5e−9 .
We compare the state variable γ when ε ∈ {1, 10} in Fig. 3 to emphasize the impact of the singular perturbation.
We can see in Fig. 4 that the slow variable h is very well approximated by the one from the reduced-order problem (P0 ). In Fig. 5, we see that what we call the outer solution of the fast variable in the singular perturbation theory, is very well approximated by the singular optimal control of problem (P0 ). Besides, we get a very small relative gap, around 0.14%, between the final times associated to both solutions of problems (P0 ) and (P1 ).
Λ
t¯1
t
t¯2
Fig. 2. The determinant Λ := det (J, F0 , F1 ), with J(t¯1 ) = F1 (x(t¯1 )), along the singular arc for the trajectory of the form σ− σs σ+ associated to y¯. 4. MODEL WITH SINGULAR PERTURBATION We are interested in this section in the following optimal control problem: min tf , ˙ = Fε (x(t), u(t)), u(t) ∈ R, t ∈ [0 , tf ] a.e., (Pε ) x(t) x(0) = x0 , x(tf ) = xf , where x := (h, v, m, γ), u is the lift coefficient, the initial state is x0 := (3480, 151.67, 69000, 0), the final state is xf := (9144, 191, 68100, 0) and where the dynamics Fε is given by equations (1)-(2)-(3)-(4’), ε ≥ 1 fixed. The values of the constant parameters are given by table 1. 4.1 Multiple shooting and homotopy on ε The application of the PMP gives the maximizing control: pγ uε (x, p) := 2 ε pv v CD2 and plugging uε in H(x, p, u) := p, Fε (x, u) gives the true Hamiltonian hε (x, p) := p, Fε (x, uε (x, p)). We use multiple shooting to deal with numerical instability due to the singular perturbation. We note (ti , zi ) the intermediate discretized times and points. The times ti are fixed and
4.2 Second order conditions of optimality We check now the second order conditions of optimality along the path of zeros. For a fixed ε ∈ [1 , 100] we denote by pε0 the initial costate solution of the associated shooting equation. Then, we compute the three Jacobi fields Ji (·) := (δxi (·), δpi (·)) such that δxi (0) = 0 and where (δp1 (0), δp2 (0), δp3 (0)) is a basis satisfying dhε (x0 , pε0 ) · δpi (0) = 0. The first conjugate time tc > 0 is the first time such that det (δx1 , δx2 , δx3 , Fε ) vanishes. The determinant being very big when ε is close to 1, we prefer to use a SVD decomposition and compute the smallest singular value. ε We denote by σmin (t) the minimal singular value at time t ε and σmax (t) the maximal one. Then, tc > 0 is a conjugate ε time if and only if σmin (tc ) = 0. According to Fig. 6 the BC-extremal along the path of zeros are locally optimal for ε ∈ [1.5 , 100] 7 . When ε = 1, we cannot conclude about the local optimality because of significant round-off errors, see Fig. 7. 5. CONCLUSION Combining theoretical and numerical tools from geometric control, we get a σ− σs σ+ hyperbolic trajectory, solution of the PMP for the reduced-order problem (P0 ). Despite the singular perturbation, we also present BC-extremals 7
1665
ε (·) for ε ∈ {1.5, 2, 5, 10} for readibility. We only represent σmin
Proceedings of the 20th IFAC World Congress 1630 Olivier Cots et al. / IFAC PapersOnLine 50-1 (2017) 1625–1630 Toulouse, France, July 9-14, 2017
satisfying second order conditions of optimality for problems (Pε ), ε > 0. In conclusion, the numerical comparison between problems (P0 ) and (P1 ) shows that the reducedorder problem is a very good approximation of the original time-to-climb problem.
ε = 1.5
ε=2 ε σmin
ε=5 ε = 10
γ
ε = 10
Normalized time ε Fig. 6. The singular value σmin (·) for ε ∈ {1.5, 2, 5, 10}. ε=1
Normalized time
Fig. 3. The air slope γ(·) for ε = 10 (dashed line) and ε = 1 (solid line).
1 σmin
1 σmax
Normalized time h
ε ε Fig. 7. The singular values σmin (·) and σmax (·) for ε = 1.
Normalized time
Fig. 4. The altitude h(·) for ε = 1 (blue dashed line) and ε = 0 (red solid line).
γ
Normalized time
Fig. 5. The air slope γ(·) for ε = 1 (blue dashed line) and ε = 0 (red solid line). When ε = 0, γ(·) is the control law of problem (P0 ). REFERENCES A. A. Agrachev & Y. L. Sachkov, Control theory from the geometric viewpoint, vol 87 of Encyclopaedia of Mathematical Sciences, Springer-Verlag, Berlin (2004), xiv+412. M. D. Ardema, Solution of the Minimum Time-to-Climb Problem by Matched Asymptotic Expansions, AIAA Journal of Guidance, Control, and Dynamics, 14 (1976), no. 7, 843–850.
F. J. Bonnans, P. Martinon & V. Gr´elard. Bocop - A collection of examples, Technical report, INRIA, 2012. RR-8053. B. Bonnard & M. Chyba, Singular trajectories and their role in control theory, vol 40 of Mathematics & Applications, Springer-Verlag, Berlin (2003), xvi+357. B. Bonnard & I. Kupka, Th´eorie des singularit´es de l’application entr´ee/sortie et optimalit´e des trajectoires singuli`eres dans le probl`eme du temps minimal, Forum Math., 5 (1993), no. 2, 111–159. J.-B. Caillau, O. Cots & J. Gergaud, Differential continuation for regular optimal control problems, Optimization Methods and Software, 27 (2011), no. 2, 177–196. A. J. Krener, The high order maximal principle and its application to singular extremals, SIAM J. Control Optim., 15 (1977), no. 2, 256–293. I. Kupka. Geometric theory of extremals in optimal control problems. i. the fold and maxwell case, Trans. Amer. Math. Soc., 299 (1987), no. 1, 225–243. N. Mo¨ıss´eev, Probl`emes math´ematiques d’analyse des syst`emes, Translated from the Russian D. Embarek, edited by Mir, Tomash Publishers, Moscow, 1985. N. Nguyen, Singular arc time-optimal climb trajectory of aircraft in a two-dimensional wind field, In AIAA Guidance, Navigation and Control Conference and Exhibit, Guidance, Navigation, and Control and Co-located Conferences, august 2006. L. S. Pontryagin, V. G. Boltyanski˘ı, R. V. Gamkrelidze & E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Translated from the Russian by K. N. Trirogoff, edited by L. W. Neustadt, Interscience Publishers John Wiley & Sons, Inc., New York-London, 1962.
1666