Analytical Initialization of a Continuation-Based Indirect Method for Optimal Control of Endo-Atmospheric Launch Vehicle Systems

TL;DR

This paper introduces an analytical initialization strategy for indirect optimal control methods applied to endo-atmospheric launch vehicles, combining simplified analytical guidance with continuation techniques to improve convergence.

Contribution

It presents a novel approach that integrates analytical solutions and continuation algorithms to enhance indirect optimal control of launch vehicles.

Findings

Successful convergence to optimal solutions demonstrated.

Improved initial guess reduces computation time.

Method applicable to complex launch vehicle dynamics.

Abstract

In this paper, we propose a strategy to solve endo-atmospheric launch vehicle optimal control problems using indirect methods. More specifically, we combine shooting methods with an adequate continuation algorithm, taking advantage of the knowledge of an analytical solution of a simpler problem. This procedure is resumed in two main steps. We first simplify the physical dynamics to obtain an analytical guidance law which is used as initial guess for a shooting method. Then, a continuation procedure makes the problem converge to the complete dynamics leading to the optimal solution of the original system. Numerical results are presented.

Equations

$$\dot{r}$$

$$\dot{v}$$

$$\dot{\gamma }$$

$$\dot{\chi }$$

$$\dot{m}$$

$$\displaystyle\displaystyle\min\$$

$$\dot{\mathbf{x}}(t)$$

$$\mathbf{x}(0)$$

$$f^{0}=\displaystyle\sigma u^{2}-\bigg{(}\frac{f_{T}}{m}\cos(\alpha)-(d+\eta c_{m}u^{2})v^{2}-g\sin(\gamma)\bigg{)}$$

$$\dot{\mathbf{x}}=\frac{\partial H}{\partial \mathbf{p}},\dot{\mathbf{p}}=-\frac{\partial H}{\partial \mathbf{x}},H(t,\mathbf{x},\mathbf{p},\mathbf{u})=\underset{\mathbf{v}\in U}{\max }H(t,\mathbf{x},\mathbf{p},\mathbf{v})\text{vskip}6.0ptplus2.0ptminus2.0pt$$

$$\underset{\mathbf{v}\in U}{\max }H(t_f,\mathbf{x}(t_f),\mathbf{p}(t_f),\mathbf{v})=0$$

$$\mathbf{p}(0)\perp T_{{\mathbf{x}}_0}{M}_0,\mathbf{p}(t_f)\perp T_{{\mathbf{x}}_{t_f}}{M}_f$$

$$\dot{\mathbf{x}}(t)$$

$$\dot{\mathbf{p}}(t)$$

$$\min$$

$$\begin{array}[]{c}\displaystyle H=p_{r}\sin(\gamma)+p_{L}\frac{\cos(\gamma)\cos(\chi)}{r}+p_{l}\frac{\cos(\gamma)\sin(\chi)}{r\cos(L)}\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \displaystyle+p_{\gamma}c_{m}u_{1}+p_{\chi}\frac{c_{m}u_{2}}{\cos(\gamma)}+p^{0}(d+\eta c_{m}(u_{1}^{2}+u_{2}^{2}))\end{array}$$

$$p_{\gamma }=-2\eta p^0{u}_1,p_{\chi }=-2\eta p^0\cos (\gamma )u_2$$

$$p_r^{\prime }$$

$$u_1^{\prime \prime }=\frac{p_{\gamma }^{\prime \prime }}{2\eta }=\frac{d}{2\eta }\left(c_m{u}_1+\frac{\cos (\gamma )}{h_r}\right)-\frac{1}{2\eta r^2}\left(p_L\cos (\chi )+p_l\frac{\sin (\chi )}{\cos (L)}\right)\text{vskip}6.0ptplus2.0ptminus2.0pt$$

$$+P_1(u_1,u_2,u_1^{\prime },u_2^{\prime })$$

$$\{\begin{array}{l}{u}_1^{\prime \prime }=\frac{c_{m}d}{2\eta }\left(u_1+\frac{\cos (\gamma )}{c_m{h}_r}\right)+P_1(u_1,u_2,u_1^{\prime },u_2^{\prime })\text{vskip}6.0ptplus2.0ptminus2.0pt\\ u_2^{\prime \prime }=\frac{c_{m}d}{2\eta }{u}_2+P_2(u_1,u_2,u_1^{\prime },u_2^{\prime })\end{array}$$

$$W^{1,\infty }(0,t_f;{\mathbb{R}}^p):=\{\mathbf{z}\in L^{\infty }(0,t_f;{\mathbb{R}}^p):{\mathbf{z}}^{\prime }\in L^{\infty }(0,t_f;{\mathbb{R}}^p)\}$$

$$\mathbf{f}=\mathbf{f}(\mathbf{y},\mathbf{z})={\mathbf{f}}_0(\mathbf{y})+\sum _{i=1}^m{z}_i{\mathbf{f}}_i(\mathbf{y})$$

$${\mathbf{y}}^{\prime }(s)=A(s)\mathbf{y}(s)+B(s)\mathbf{z}(s)$$

$$\tiny\left(\begin{array}[]{ccccc}0&0&0&\cos(\gamma_{0})&0\\ -\frac{\cos(\gamma_{0})\cos(\chi_{0})}{r^{2}}&0&0&-\frac{\sin(\gamma_{0})\cos(\chi_{0})}{r}&-\frac{\cos(\gamma_{0})\sin(\chi_{0})}{r}\\ -\frac{\cos(\gamma_{0})\sin(\chi_{0})}{r^{2}\cos(L)}&\frac{\cos(\gamma_{0})\sin(\chi_{0})\tan(L)}{r\cos(L)}&0&-\frac{\sin(\gamma_{0})\sin(\chi_{0})}{r\cos(L)}&\frac{\cos(\gamma_{0})\cos(\chi_{0})}{r\cos(L)}\\ 0&0&0&0&0\\ 0&0&0&0&0\end{array}\right)$$

$$\scriptsize B(s)=\footnotesize\left(\begin{array}[]{ccccc}0&0&0&c_{m}(r)&0\\ 0&0&0&0&\frac{c_{m}(r)}{\cos(\gamma_{0})}\\ \end{array}\right)^{\top}$$

$$\text{Span}\{B_i(s)\mathbf{w}:\mathbf{w}\in {\mathbb{R}}^p,i\in \mathbb{N}\}={\mathbb{R}}^n$$

$$\textnormal{det}\big{(}B_{0}(s)\mathbf{w}_{1},B_{0}(s)\mathbf{w}_{2},B_{1}(s)\mathbf{w}_{1},B_{1}(s)\mathbf{w}_{2},B_{2}(s)\mathbf{w}_{1}\big{)}=\frac{c_{m}^{5}}{r^{3}\cos(L)}$$

$$\mathbf{i}$$

$$\footnotesize\left(\begin{array}[]{ccc}\cos(\gamma-\lambda_{1})&-\sin(\gamma-\lambda_{1})&-\sin(\chi-\lambda_{2})\cos(\gamma)\\ \sin(\gamma-\lambda_{1})&\cos(\gamma-\lambda_{1})&\sin(\chi-\lambda_{2})\sin(\gamma)\\ \sin(\chi-\lambda_{2})\cos(\gamma)&-\sin(\chi-\lambda_{2})\sin(\gamma)&\cos(\chi-\lambda_{2})\end{array}\right)$$

$$(\mathbf{v})_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}=\left(\footnotesize\begin{array}[]{c}v\cos(\gamma-\lambda_{1})\\ v\sin(\gamma-\lambda_{1})\\ v\sin(\chi-\lambda_{2})\cos(\gamma)\\ \end{array}\right)$$

$$\left(\frac{d}{dt}\mathbf{n}\right)_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}=\left(\footnotesize\begin{array}[]{c}0\\ \displaystyle-\frac{v}{R}\sin(\gamma-\lambda_{1})\\ \displaystyle-\frac{v}{R}\sin(\chi-\lambda_{2})\cos(\gamma)\\ \end{array}\right)$$

$$\frac{d}{dt}\mathbf{n}$$

$$\left(\frac{d}{dt}\mathbf{n}\right)_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}=\left[\frac{d}{dt}(M_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})})\cdot M_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}^{\top}\right]^{\top}\left(\footnotesize\begin{array}[]{c}1\\ 0\\ 0\\ \end{array}\right)$$

$$\displaystyle\left[\frac{d}{dt}(M_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})})\right]^{\top}\left(\scriptsize\begin{array}[]{c}1\\ 0\\ 0\\ \end{array}\right)=\frac{d}{dt}(\mathbf{n})_{(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})}\vskip 6.0pt plus 2.0pt minus 2.0pt$$

Full text

1 Introduction

Guidance of autonomous launch vehicles towards rendezvous points is a complex task often considered in missile applications, mainly for interception of targets. It represents an optimal control problem, whose aim consists in finding a control law enabling the vehicle to join a final point of the 3D space considering prescribed constraints as well as performance criteria. The rendezvous point may be a static point, for example an intermediate point at the middle of the flight mission, as well as a moving point if, for example, the task consists in intercepting a maneuvering target. This requires not only high numerical precision of concerned algorithms but also a real-time processing of optimal trajectories.

The most common approach to solve this kind of task resides on analytical guidance laws. They correct errors coming from perturbations and misreading of the system. For interceptor missile applications, in [Lin and Tsai, 1987], [Lin, 1991], [Indig et al., 2013] guidance laws that maximize the final velocity and so the lethality of the interceptor are designed. Varying velocity is also considered for the design of an energy optimal guidance law in [Song et al., 1999]. However, the trajectories induced by guidance laws are not optimal because of some approximations made.

Ensuring the optimality of trajectories can be achieved exploiting direct methods. These techniques consist in discretizing each component of the optimal control problem (the state, the control, etc.) to reduce the whole mathematical representation to a nonlinear constrained optimization problem. Since they are quite robust, they are widely used [Hargraves and Paris, 1987], [Ross and Fahroo, 2002], [Ross et al., 2003]. However, these methods are computationally demanding and can often be used offline uniquely.

In order to manage efficiently real-time processing of optimal control sequences for launch vehicle systems one may consider indirect methods. These use a mathematical study of the system (exploiting the Pontryagin Maximum Principle) to determine some necessary conditions of optimality. Indirect methods converge much faster than direct methods with a better precision. Since the problem is equivalent to the research of zeros of a function, the main difficulty remains their initialization. In [Lu et al., 2003], [Pan and Lu, 2010] the initialization problem is bypassed using finite differences algorithms and multiple-shooting methods respectively. By contrast, in [Pontani and Cecchetti, 2013] the authors use second-order conditions and conjugate point theory for a multistage launch vehicle problem; the initialization of the numerical method is based on the particle swarm algorithm. However, most of these approaches remain computationally demanding and not easily applicable in view of real-time processing.

In this paper, we propose to solve endo-atmospheric launch vehicle optimal control problems using indirect methods managing the issue coming from the initialization by both an analytical guidance law and a continuation method. Continuation procedures have shown to be reliable and robust for problems like coplanar orbit transfer and atmospheric reentry [Cerf et al., 2012], [Graichen and Petit, 2008]. This combination allows to preserve precision and fast numerical computations. The proposed approach is resumed in two main steps. We first simplify the physical dynamics to obtain an analytical guidance law which is used as initial guess for a shooting method. Then, a continuation procedure makes the problem converge to the complete dynamics leading to the optimal solution of the original system.

The paper is structured as follows. Section 2 introduces the dynamics of a general endo-atmospheric launch vehicle system, the optimal control formulation and the related numerical approach. Section 3 is devoted to the construction of a simplified dynamics able to initialize successfully a shooting method and to make the final algorithm converge efficiently to the complete dynamics. In Section 4 the continuation method to recover the original dynamics is presented with its indirect formulation and numerical tests. Finally, Section 5 contains conclusions and perspectives.

2 Dynamics, Optimal Control Problem and Numerical Method

2.1 Physical Model

Let $(\mathbf{I},\mathbf{J},\mathbf{K})$ be an inertial frame centered at the center of the planet $O$, $(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})$ be the NED frame and $(\mathbf{i},\mathbf{j},\mathbf{k})$ be the velocity frame. The endo-atmospheric launch vehicle system is modeled as an axisymmetric thrust propelled rigid body of mass $m$ [Pucci et al., 2015]. The coordinates $(r,l,L,v,\gamma,\chi)~{}\in~{}\mathbb{R}^{6}$ ($L$ is the latitude, $l$ is the longitude, $\gamma$ is the path angle and $\chi$ is the azimuth angle) are used to represent the position $\bm{\xi}=(r\cos(L)\cos(l),r\cos(L)\sin(l),r\sin(L))$ of the center of mass $G$ of the vehicle and its velocity $\mathbf{v}~{}=~{}v\cos(\gamma)\cos(\chi)\mathbf{e}_{L}+v\cos(\gamma)\sin(\chi)\mathbf{e}_{l}-v\sin(\gamma)\mathbf{e}_{r}$.

Neglecting the wind velocity, the Coriolis and the centripetal force (this is legitimate because of the short length of the considered trajectories), the dynamics takes the following form

[TABLE]

where $g$ is the modulus of the gravity, $\eta$ is an aerodynamic efficiency factor, $\alpha=\alpha_{\max}u$ is the angle of attack while $\beta$ is the angle of bank, $u$ stands for the normalized lift coefficient while $q=q(t)$ is the mass flow and $f_{T}=f_{T}(t)$ represents the modulus of the propulsion which depends directly on $q(t)$.

Based on a standard model of flight dynamics [Pucci et al., 2015], [Pepy and Hérissé, 2014], coefficients $d$ and $c_{m}$ are approximated by $d~{}=~{}d(r)~{}~{}=~{}\frac{1}{2m}\rho SC_{D_{0}}$ and $c_{m}=c_{m}(r)=\frac{1}{2m}\rho SC_{L_{\max}}$ where $S$ is the reference area, $C_{L_{\max}}$ is the maximal value of the lift coefficient and $C_{D_{0}}$ is the drag coefficient for $\alpha=0$, which is considered constant; finally, $\rho$ is the air density for which an exponential model $\rho~{}=~{}\rho(r)~{}=~{}\rho_{0}\exp(-(r~{}-~{}r_{T})/h_{r})$ is considered, where $h_{r}$ is a reference altitude. Since $q(t)$ is a predefined function of time, in this paper controls are only $u$ and $\beta$.

2.2 Optimal Control Problem and Maximum Principle

Consider now the Optimal Control Problem (OCP)

[TABLE]

where $\mathbf{x}(t)=(r,L,l,v,\gamma,\chi)(t)\in\mathbb{R}^{6}$, $\mathbf{u}(t)=(u,\beta)(t)\in\mathbb{R}^{2}$, $\mathbf{f}$ is the mapping defined by dynamics (2.1), $M_{0}$, $M_{f}$ are smooth submanifolds in $\mathbb{R}^{6}$ and the transfer time $t_{f}$ is not fixed. Finally,

[TABLE]

where $\sigma\geq 0$ is constant. By definition, $u^{2}$ takes its values in $[0,1]$. However, we do not consider any boundaries on $u$ preferring to penalize it, using $\sigma$, $\eta c_{m}v^{2}$ within the cost.

The Pontryagin Maximum Principle (PMP) [Boltyanskiy et al., 1962] states that, if $\mathbf{u}$ is optimal with response defined on $[0,t_{f}]$, and shortly denoted $\mathbf{x}(t)$, then there exists $p^{0}\leq 0$ and $\mathbf{p}~{}\in~{}AC([0,t_{f}],\mathbb{R}^{6})$ such that a.e. on $[0,t_{f}]$

[TABLE]

where $H=\left\langle\mathbf{p},\mathbf{f}(t,\mathbf{x},\mathbf{u})\right\rangle+p^{0}f^{0}(t,\mathbf{x},\mathbf{u})$ is the Hamiltonian and $\mathbf{p}$ satisfies the transversality conditions

[TABLE]

Treating (OCP) by indirect methods consists in solving

[TABLE]

with an appropriate value for $\mathbf{p}_{0}$.

2.3 Shooting and Continuation Method

It is known [Trélat, 2008] that finding a solution of (6) can be reduced to solve $G(\mathbf{p}_{0},t_{f})=0\;,\ G:\mathbb{R}^{n+1}\longrightarrow\mathbb{R}^{n+1}$ ($G$ is called shooting function) using Newton-type methods. This is the content of the well known shooting method in optimal control [Trélat, 2008], [Trélat, 2012], [Stoer and Bulirsch, 2013]. Its advantage is its extremely good numerical accuracy, relevant for aerospace applications [Trélat, 2012]. Since it relies on the Newton method, it inherits of the very quick convergence properties of the Newton method. The main drawback of the shooting method is that it may be difficult to initialize.

To overcome this difficulty, one can entrust with the robustness of the continuation method. It consists in deforming the problem into a simpler one that we are able to solve and then in solving a series of shooting problems, step by step by parameter deformation, to recover the original problem [Allgower and Georg, 2003]. This approach increases the efficiency of the shooting method because it allows to relax its initialization. The continuation parameter $\lambda$ may be a physical parameter (or several) of the problem, or an artificial one. The path consists of a convex combination of the simpler problem and of the original one, with $\lambda\in[0,1]$.

The main algorithm consists then in finding a solution of some simplification of (OCP) first and, from this, solving by continuation the original formulation (OCP).

3 New Guidance Law as Good Estimate for Continuation Method

Continuation methods allow us to solve iteratively (OCP) once the solution of some (usually) simpler optimal control problem is known. Here, we introduce and treat an efficient simpler problem coming from a modification of (2.1).

This modified version of (OCP) is designed with the hope that, on one hand, the shooting method can be easily initialized and, on the other hand, the continuation is feasible. It is at this stage that intuition of aerospace engineer is of primary importance, in order to design simplified problems and homotopy parameters resulting in a meaningful continuation procedure.

3.1 Simplified Problem and Analytical Smooth Controls

If one ignores the contributions of the curvature of the Earth, of the gravity and of the propulsion within (2.1), introducing the curvilinear abscissa $s(t):=\int_{0}^{t}v(t^{\prime})\;dt^{\prime}$ and a new variable $w:=\ln(v)$, the following simplification of (OCP) is obtained (the equation of $w^{\prime}$ can be neglected because it does not influence anymore the dynamics)

[TABLE]

where $u_{1}:=u\cos(\beta)$, $u_{2}:=u\sin(\beta)$ and with initial, final conditions $(r_{0},L_{0},l_{0},v_{0},\gamma_{0},\chi_{0})$, $(r_{f},L_{f},l_{f},\gamma_{f},\chi_{f})$.

The next step consists in finding an analytical solution of (3.1). This is achieved exploiting the PMP formulation, as showed in the following.

The associated Hamiltonian function is

[TABLE]

It can be remarked that, since the transfer time is not fixed and formulation (3.1) is autonomous, (8) takes zero as value for all times $t\in[0,t_{f}]$ [Boltyanskiy et al., 1962]. Since no constraints are considered on $u_{1}$, $u_{2}$, from applying the weaker version of PMP [Boltyanskiy et al., 1962], it follows explicitly that

[TABLE]

Simple calculations show that

[TABLE]

A further analysis streamlines this formulation letting one considering only normal trajectories.

Proposition 1

Assume that $\cos(\gamma)\neq 0$ a.e. in $[0,t_{f}]$. Then, no abnormal trajectories arise in (3.1), i.e. $p^{0}\neq 0$.

{pf}

If $p^{0}=0$, then from (9) it follows $p_{\gamma}=p_{\chi}\equiv 0$ in $[0,t_{f}]$. Using $H=0$, $p_{\gamma}^{\prime}=0$ and $p_{\chi}^{\prime}=0$ raises a homogeneous linear system of equations in $(p_{r},p_{L},p_{l})$ whose associated matrix has $-\frac{\cos(\gamma)}{r^{2}\cos(L)}$ as determinant which is not zero almost everywhere in $[0,t_{f}]$ by assumption. Then $(p_{r}(t),p_{L}(t),p_{l}(t),p_{\gamma}(t),p_{\chi}(t),p^{0})\equiv 0$ on $[0,t_{f}]$ which raises a contradiction with the PMP. □

The assumption $\cos(\gamma)\neq 0$ a.e. in $[0,t_{f}]$ is relevant because, if the trajectory takes $\cos(\gamma)=0$ as optimal value, a change of coordinates like $\tilde{\gamma}=\gamma+\frac{\pi}{2}$ can be done on the dynamics during the numerical integration to avoid this singularity. Moreover, other changes of coordinates that optimize the computational speed can be taken into account [Zhu et al., 2016]. Thus, by Proposition 1 and the PMP, we can substitute $p^{0}=-1$ within the previous simplified formulation.

Proceeding formally under smoothness assumption for optimal controls of this simplified problem, analyzing the second derivative of $u_{1}$ and $u_{2}$ allows to seek some relations. Indeed, using (9) and (3.1) it can be easily shown that

[TABLE]

where $P_{1}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$ is a polynomial of degree greater than one of $u_{1}$, $u_{2}$ and their first derivatives. Since $\frac{1}{r^{2}}\ll 1$ the previous expression can be approximated. Iterating the same procedure on $u_{2}$ one is led to solve ($P_{2}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$ identifies again a polynomial of degree greater than one of $u_{1}$, $u_{2}$ and their first derivatives)

[TABLE]

It is clear that, if additional assumptions allow to neglect the contribution of $P_{1}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$, $P_{2}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$, equations of system (11) can be solved independently using the assumption that parameters $c_{m}$, $d$ and $\cos(\gamma)$ are almost constant along the trajectory, justifying this approximation because scenarios with limited altitude are considered to solve the simplified problem (3.1). This is not limiting, because, once the shooting method on simplified problem has converged, a continuation procedure on the final point can be started changing the initial scenario.

3.2 Local Controllability of the Simplified System

An evident sufficient condition that allows to neglect the contribution of $P_{1}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$, $P_{2}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$ within system (11) is to deal with small (in the sense of $L^{\infty}$ norm) controls $\mathbf{z}:=(u_{1},u_{2})$ that have small derivatives. This is achieved studying the local controllability around control $\bar{\mathbf{z}}=(0,0)$ in

[TABLE]

where $\mathbf{z}^{\prime}$ is the first distributional derivative of $\mathbf{z}$. The set of admissible controls on $[0,t_{f}]$ is denoted by $\mathcal{U}_{\mathbf{y}_{0},t_{f},\mathbb{R}^{p}}$ where $\mathbf{y}_{0}$ denotes the initial point of the state variable $\mathbf{y}~{}:=~{}(r,L,l,\gamma,\chi)$.

The attention is focused on a special class of control system to which the dynamics of (3.1) belongs, the affine systems. A control system (2.2) is affine if

[TABLE]

where the $\mathbf{f}_{i}$’s are (at least) $C^{1}$ vector fields. This standard result holds (whose proof will be reported in [Bonalli et al., Possibly 2017])

Theorem 1

Let $\mathbf{y}_{0}\in\mathbb{R}^{n}$ and $\bar{\mathbf{z}}\in\mathcal{U}_{\mathbf{y}_{0},t_{f},\mathbb{R}^{p}}$. If the linearized system of (12) along $(\mathbf{y}_{\bar{\mathbf{z}}}(\cdot),\bar{\mathbf{z}}(\cdot))$ is controllable from $\mathbf{y}_{0}$ in time $t_{f}$, then (12) is locally controllable at $\mathbf{y}_{\bar{\mathbf{z}}}(t_{f})$ from $\mathbf{y}_{0}$ in time $t_{f}$ using controls $\mathbf{z}$ with $\|\mathbf{z}-\bar{\mathbf{z}}\|_{W^{1,\infty}}$ small enough.

All these results can be applied to the dynamics of (3.1) allowing to give conditions to solve easily (11). The idea is to compute the trajectory solution of (3.1) given by control $\bar{\mathbf{z}}~{}=~{}(0,0)\in~{}W^{1,\infty}(0,t_{f};\mathbb{R}^{2})$ and to consider scenarios with final point close to $\mathbf{y}_{\bar{\mathbf{z}}}(t_{f})$. If the linearized system along $(\mathbf{y}_{\bar{\mathbf{z}}}(\cdot),\bar{\mathbf{z}}(\cdot))$ is controllable from $\mathbf{y}_{0}$ in time $t_{f}$, then, thanks to Theorem 1, the dynamics of (3.1) is locally controllable at $\mathbf{y}_{\bar{\mathbf{z}}}(t_{f})$ from $\mathbf{y}_{0}$ in time $t_{f}$ with small controls in $W^{1,\infty}(0,t_{f};\mathbb{R}^{2})$. In other words, if the considered scenarios have final point close to $\mathbf{y}_{\bar{\mathbf{z}}}(t_{f})$ we can neglect the contribution of $P_{1}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$, $P_{2}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})$ within system (11) to solve it analytically. It must be noted that it is not restrictive considering scenarios with final point close to $\mathbf{y}_{\bar{\mathbf{z}}}(t_{f})$ to solve (3.1), because, as said before, a continuation procedure on the final point can be started changing considerably the initial scenario.

The linearized system of (3.1) has the form

[TABLE]

where $A(s)$ is the matrix

[TABLE]

and

[TABLE]

This system is controllable in any $s_{f}$. To show this, we consider the well known fact [Trélat, 2008] that, for a system with no constraint $\mathbf{y}^{\prime}(s)=A(s)\mathbf{y}(s)+B(s)\mathbf{z}(s)$ where $s~{}\mapsto~{}A(s)$ and $s\mapsto B(s)$ are of class $C^{\infty}$, if, defined the sequence $B_{0}(s)~{}:=~{}B(s)$, $B_{i+1}(s):=A(s)B_{i}(s)-\frac{dB_{i}}{ds}(s)$ for $i\in\mathbb{N}$, there exists $s\in[0,s_{f}]$ such that

[TABLE]

then the system is controllable in time $s_{f}$. Denoting $\mathbf{w}_{1}~{}=~{}(1,0)^{\top}$, $\mathbf{w}_{2}=(0,1)^{\top}$, it is easy to show that

[TABLE]

which is different from zero. Then, the claim follows.

This result validates the previous approach. Initial conditions for controls and their first derivatives must be sought to solve (11) (where we impose then $P_{1}=P_{2}=0$). Two conditions are given using the two last equations of (3.1). The second ones can be investigated exploiting first an analysis of the Line Of Sight (LOS).

3.3 LOS Analysis

The line of sight is the vector joining the current position $\bm{\xi}$ to the desired final point $\bm{\xi}_{f}$. It results to be useful defining its modulus $R~{}:=\|\bm{\xi}_{f}-\bm{\xi}\|$ and the associated unitary vector $\mathbf{n}~{}:=~{}\frac{\bm{\xi}_{f}-\bm{\xi}}{R}$. We denote $(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})$ the orthonormal frame where $\mathbf{i}_{n}=\mathbf{n}$, $\mathbf{j}_{n}$ is the unitary vector in the plane $(\mathbf{i}_{n},\mathbf{e}_{r})$ perpendicular to $\mathbf{i}_{n}$ and oriented by $\mathbf{j}_{n}\cdot\mathbf{e}_{r}<0$ and $\mathbf{k}_{n}=\mathbf{i}_{n}\wedge~{}\mathbf{j}_{n}$. Frames $(\mathbf{i},\mathbf{j},\mathbf{k})$, $(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})$ are functions of $(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})$ by definition [CNES, 1993], locally according to

[TABLE]

where $\lambda_{1}$, $\lambda_{2}$ are Euler angles. We denote $M_{(\mathbf{i},\mathbf{j},\mathbf{k})}$, $M_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}$ the associated transition matrices from frame $(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})$ to frame $(\mathbf{i},\mathbf{j},\mathbf{k})$ and from frame $(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})$ to frame $(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})$ respectively. The coordinates of a vector within a frame $(\mathbf{a},\mathbf{b},\mathbf{c})$ are denoted as $(\cdot)_{(\mathbf{a},\mathbf{b},\mathbf{c})}$. The following analysis seeks some approximate relations between the derivatives of $\lambda_{1}$, $\lambda_{2}$ and angles $\gamma$, $\chi$ under appropriate assumptions.

First-Order Approximation 1** **(FOA)

Along the considered trajectories, only weak variations from the LOS direction are allowed, i.e. $\gamma\cong\lambda_{1}$ , $\chi\cong\lambda_{2}$. Moreover, the contribution of the rotation of frame $(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})$ around frame $(\mathbf{I},\mathbf{J},\mathbf{K})$ is neglected, i.e. $\frac{d}{dt}(\mathbf{n})_{(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})}\cong\left(\frac{d}{dt}\mathbf{n}\right)_{(\mathbf{e}_{L},\mathbf{e}_{l},\mathbf{e}_{r})}$.

It must be remarked that this assumption is not a strong one because the initial scenario involved lives in a neighborhood of the null control and the considered trajectories are short enough to make this choice effective. It allows to manipulate transition maps obtaining easily the desired results.

Proposition 2

Under (FOA), the transition matrix $M$ from frame $(\mathbf{i},\mathbf{j},\mathbf{k})$ to frame $(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})$ is

[TABLE]

Moreover, we have

[TABLE]

{pf}

Matrix $M$ is obtained easily by approximating $\cos(\chi-\lambda_{2})\cong 1$ within the changing variable rotation $M_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}\circ M_{(\mathbf{i},\mathbf{j},\mathbf{k})}^{\top}$. The second relation follows by multiplying matrix $M$ by $(v,0,0)^{\top}$. □

The next step consists of making explicit the coordinates of the derivative of the unitary vector $\mathbf{n}$ within frame $(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})$. This requires two lemmas.

Lemma 1

It holds $\frac{d}{dt}\mathbf{n}=-\frac{1}{R}\left(I-\mathbf{n}\ \mathbf{n}\;^{\top}\right)\mathbf{v}$. Moreover, under (FOA), one has

[TABLE]

{pf}

Differentiating the definition of $\mathbf{n}$ one has

[TABLE]

The second expression follows by substituting $(\mathbf{v})_{(\mathbf{i}_{n},\mathbf{j}_{n},\mathbf{k}_{n})}$ within the previous equation. □

Lemma 2

Under (FOA), we have

[TABLE]

{pf}

One has

[TABLE]

Some calculations lead to

[TABLE]

Gathering together all the previous results, it is straightforward to show the following final result

[TABLE]

Relations (14) turn out to be useful to obtain initial condition for analytical controls as shown in the followings.

3.4 New Approximated Analytical Guidance Law

The previous results can be exploited to obtain an approximate solution of (11) and consequently new approximate analytical optimal controls $u_{1}$, $u_{2}$ for (3.1). We impose $c_{m}$, $d$ and $\cos(\gamma)$ to be constant along the trajectory (thanks to considerations of Section 3.1). We start considering control $u_{1}$. From (11), we have

[TABLE]

We set $b:=\sqrt{\frac{c_{m}d}{2\eta}}$. Assumption (FOA) leads to $R~{}\cong~{}(s_{f}~{}-~{}s)$. Plugging (15) into equation $\gamma^{\prime}=c_{m}u_{1}$ and integrating, one obtains

[TABLE]

It is straightforward to see that, thanks to (FOA) and Proposition 2, one has $\dot{R}~{}=~{}-v\cos(\gamma-\lambda_{1})$. Now, (14) is used to differentiate the quantity $R\sin(\gamma-\lambda_{1})$. Indeed

[TABLE]

Then, since under usual assumption $\dot{R}\cong-v$, it follows $(R\sin(\gamma-\lambda_{1}))^{\prime}=Rc_{m}u_{1}$. Integrating this last equation, using (16) and noticing that $R(s_{f})=0$, one has

[TABLE]

where $k_{1}(R)$ and $k_{2}(R)$ are gain parameters defined as

[TABLE]

Since $\sin(\gamma-\lambda_{1})\cong\gamma-\lambda_{1}$, the following guidance law for control $u_{1}$ is deduced

[TABLE]

where $k_{3}(R)=2+k_{1}(R)-k_{2}(R)$.

With the same argumentation, exploiting again (14), the following guidance law for control $u_{2}$ can be easily derived

[TABLE]

Relations (17), (18) provide an analytical guess for $(s_{f},p_{r}(0),p_{L}(0),p_{l}(0),p_{\gamma}(0),p_{\chi}(0))$ to initialize the shooting method applied to problem (3.1). Indeed, they may be chosen as $s_{f}=R(0)$, $p_{\gamma}(0)=2\eta u_{1}(0)$, $p_{\chi}(0)=2\eta\cos(\gamma_{0})u_{2}(0)$ and the guess values of $p_{r}(0)$, $p_{L}(0)$, $p_{l}(0)$ are obtained from (8) and (3.1).

It is interesting to note that (17), (18) generalize the guidance law proposed in [Lin, 1991].

4 Solution of the Original System by a Continuation Method

Since the procedure developed in Section 3 gives us a good guess to initialize the shooting method on a simplified version of (OCP), we proceed now to the analysis of the continuation method that allows to recover an optimal solution of (OCP) starting an homotopy procedure on the simplified problem (3.1).

Clearly, several types of different continuations can be implemented. However, we verified that passing from (3.1) to (OCP) through a continuation parameter $\lambda_{1}$ and, after that, changing the final point of the initial scenario directly on (OCP) with another parameter $\lambda_{2}$ results to be enough to obtain a fast convergence to the optimal solution.

4.1 Solution of the Original Problem

Without loss of generality, we focus on a particular instance of (OCP), whose cost (see (3)) takes the form

[TABLE]

[TABLE]

The system, which depends on $\lambda_{1}$, is introduced exploiting the change of variable $w:=\ln(v)$ that leaves the formulation consistent with the simplified model (3.1) as it is shown hereafter. If we denote $\tilde{\mathbf{x}}_{f}$ the final point of the initial test scenario (see Section 3.1) and $\mathbf{x}_{f}$ the final point of the desired scenario, we have

[TABLE]

where the final point takes the form $\mathbf{x}(t_{f})=\tilde{\mathbf{x}}_{f}+\lambda_{2}(\mathbf{x}_{f}-\tilde{\mathbf{x}}_{f})$, $(\lambda_{1},\lambda_{2})\in[0,1]^{2}$ and the evolution of the mass is given by $m(t)=m_{0}-\lambda_{1}\int^{t}_{0}q(s)\;ds$. The continuation algorithm implemented makes $(\lambda_{1},\lambda_{2})$ converge from $(0,0)$ to $(1,1)$.

As usual, the Hamiltonian of (4.1) takes the form $H=\left\langle\mathbf{p},\mathbf{f}\right\rangle+f^{0}$ where we assume that $p^{0}\neq 0$ (which is coherent with the simplified formulation (3.1)). A priori, the Hamiltonian is not uniformly equal to zero because $f_{T}=f_{T}(t)$, $q=q(t)$ depend explicitly on $t$. From transversality conditions we obtain

[TABLE]

Finally, since no constraints are considered on controls $\beta$, $u$, relations $\frac{\partial H}{\partial\beta}~{}=~{}0\;,\ \frac{\partial H}{\partial u}=0$ hold. From the first one

[TABLE]

An analytical formulation of $\beta$ is immediately obtained if a first-order approximation $\cos(\alpha)\cong 1$, $\sin(\alpha)\cong\alpha$ is considered (this remains proper if small values of $\alpha_{\max}$ are considered). Indeed, from (21) it follows $p_{\chi}\frac{\cos(\beta)}{\cos(\gamma)}~{}-~{}p_{\gamma}\sin(\beta)=0$. The same reasoning leads to an explicit relation for $u$.

Variable $w$ is introduced to obtain a guess of $p_{w}(0)$ for the simplified problem that arises from condition $(\lambda_{1},\lambda_{2})~{}=~{}(0,0)$. In this case, the Hamiltonian is uniformly equal to zero since it does not depend on time. Then, from the dual equations it follows that $\dot{p}_{w}=-H=0$ and from (20) one obtains $p_{w}(t)~{}=~{}0$ in $[0,t_{f}]$.Then, the procedure achieved in the previous section provides an analytical guess to initialize the sequence of shooting methods within the continuation algorithm that solves the original optimal control problem.

The continuation algorithm is summed up as follows.

Set $(\lambda_{1},\lambda_{2})=(0,0)$, $\tilde{\mathbf{x}}_{f}=(\tilde{r}_{f},\tilde{L}_{f},\tilde{l}_{f},\tilde{\gamma}_{f},\tilde{\chi}_{f})$ the final point reached by dynamics (4.1) (equivalent to (3.1) at this stage) with control $u=0$ at $t_{f}$ and compute the analytical guess for problem (4.1) with initial and final point $\mathbf{x}_{0}$, $\tilde{\mathbf{x}}_{f}$ respectively; 2. 2.

With $\lambda_{2}=0$ fixed, start a continuation procedure with $\lambda_{1}~{}\in~{}[0,1]$ on the adjoint formulation of (4.1), with initial and final point $\mathbf{x}_{0}$, $\tilde{\mathbf{x}}_{f}$ respectively, until $\lambda_{1}~{}=~{}1$; 3. 3.

With $\lambda_{1}=1$ fixed, start a continuation procedure with $\lambda_{2}~{}\in~{}[0,1]$ on the adjoint formulation of (4.1), with initial and final point $\mathbf{x}_{0}$, $\tilde{\mathbf{x}}_{f}+\lambda_{2}(\mathbf{x}_{f}-\tilde{\mathbf{x}}_{f})$ respectively, until $\lambda_{2}=1$.

At point 1 of the previous algorithm $t_{f}$ must be guessed: it should be chosen such that the assumptions given in Section 3.1 for the integration of (11) are satisfied.

4.2 Test Cases and Numerical Results

The previous approach is tested on three test cases. The final scenarios are adopted trying to consider real mission conditions. The following numerical values are chosen: the maximal curvature at altitude zero is chosen as $c_{m}(0)~{}=~{}0.00075\;\textnormal{m}^{-1}$, the drag coefficient at altitude zero is $d(0)~{}=~{}0.00005\;\textnormal{m}^{-1}$, the aerodynamic efficiency factor is $\eta=0.442$ and the reference altitude is $h_{r}=7500\;\textnormal{m}$. Moreover, the following functions are chosen to represent respectively the mass flow and the propulsion

[TABLE]

where $q_{0}=10\;\textnormal{kg}\cdot\textnormal{s}^{-1}$, $t_{sw}=20\;\textnormal{s}$ and $v_{e}=1500\;\textnormal{m}\cdot\textnormal{s}^{-1}$ is the constant fuel injection velocity. Finally, the maximal angle of attack is $\alpha_{\max}~{}=~{}\pi/6$. These conditions are representative of a realistic interceptor scenario. The numerical calculations are computed on a machine Intel(R) Xeon(R) CPU E5-1607 v2 @ 3.00GHz, with 8.00 Gb of RAM.

We impose an initial velocity of $v_{0}=1000\;\textnormal{m/s}$. If we represent the state by $\mathbf{y}=(r,L,l,\gamma,\chi)$ (the initial velocity is fixed), the three scenarios are the following (in standard units, where $r_{T}$ is the Earth radius)

• $\bm{S_{1}}:$

$\mathbf{y}_{0}=(r_{T}+3000,5454661/r_{t},46086/r_{t},-\pi/6,0)$ ,

$\mathbf{y}_{f}=(r_{T}+12000,5475000/r_{t},42000/r_{t},0,\pi/8)$ ;

• $\bm{S_{2}}:$

$\mathbf{y}_{0}=(r_{T}+3000,5454661/r_{t},46086/r_{t},\pi/4,0)$ ,

$\mathbf{y}_{f}=(r_{T}+12000,5485000/r_{t},36178/r_{t},-\pi/4,-\pi/2)$ ;

• $\bm{S_{3}}:$

$\mathbf{y}_{0}=(r_{T}+3000,5454661/r_{t},46086/r_{t},0,0)$ ,

$\mathbf{y}_{f}=(r_{T}+3000,5485000/r_{t},46086/r_{t},0,0)$ .

The associated optimal values recovered by the proposed approach are (in standard units) $(v(t_{f}),t_{f})_{\bm{S_{1}}}=(986.7,24.5)$, $(v(t_{f}),t_{f})_{\bm{S_{2}}}=(851.6,36.6)$ and $(v(t_{f}),t_{f})_{\bm{S_{3}}}$ $=(688.8,31.5)$. Using around 400 time steps within the explicit second-order Runge-Kutta solver for the integration of the ODEs inside the shooting function, the simulations take on average $0.27\;\textnormal{s}$ for $\bm{S_{1}}$, $0.73\;\textnormal{s}$ for $\bm{S_{2}}$ and $2.1\;\textnormal{s}$ for $\bm{S_{3}}$. These computing times remain in the order of $1\;\textnormal{Hz}$ giving the possibility to consider the whole procedure as a good approach for a real-time processing.

The same scenarios were treated also with a non-initialized direct method (AMPL/IPOPT with 200 time steps) [Fourer et al., 1993]. Modifying the boundaries and the initial guess of IPOPT, the three scenarios are solved by the direct method (obtaining the same optimal solutions provided previously) with computational times that vary respectively between $1\;\textnormal{s}$ and $2\;\textnormal{s}$ for $\bm{S_{1}}$, $2\;\textnormal{s}$ and $5\;\textnormal{s}$ for $\bm{S_{2}}$, $5\;\textnormal{s}$ and $10\;\textnormal{s}$ for $\bm{S_{3}}$. Our indirect approach reveals itself to be at least as efficient as a direct method, and sometimes more successful.

It is interesting to note that scenarios $\bm{S_{1}}$, $\bm{S_{2}}$ do not require step 3 of the previous algorithm to be solved: imposing $\tilde{\mathbf{x}}_{f}~{}:=~{}\mathbf{x}_{f}$ within step 1 , the initial guess is good enough to find the solution of (OCP) with a continuation on $\lambda_{1}$ only. From this, 3 iterations of the continuation method on $\lambda_{1}$ are required to solve $\bm{S_{1}}$, 10 iterations on $\lambda_{1}$ for $\bm{S_{2}}$ while 18 iterations on $\lambda_{1}$ and 6 iterations on $\lambda_{2}$ are needed to solve $\bm{S_{3}}$.

The trajectories and the optimal controls $u_{1}~{}:=~{}u\cos(\beta)$, $u_{2}~{}:=~{}u\sin(\beta)$ are shown in Figure 1. It is interesting to note that, even if no explicit constraint are considered on both the controls, they turn out to be bounded in $[-1,1]$. They are forced to be small due to the presence of $u^{2}$ inside the cost function.

5 Conclusions and Perspectives

In this paper we solve an optimal control problem for endo-atmospheric launch vehicle systems using indirect methods. The problem is solved mixing an analytical initialization with a continuation method on the dynamics. First, a simplified version of the original problem is solved. Then, a continuation method which recovers the complete optimal solution starts: every shooting method within the iterative procedure is initialized by the approximate solution found at the previous stage. This approach provides numerical precision and high computational. Moreover, numerical simulations show that this approach may compete with usual direct methods.

Clearly, some improvements can be brought. A first remark concerns the model. Indeed, the velocity of the wind and stabilizing components within the cost such as $\sigma u^{2}$ have not been considered in the previous analysis. However, such terms can be added to the formulation exploiting a new stage of continuation on them. Secondly, in this contribution, neither control constraints nor state constraints were taken into account. It is not hard to consider control constraints (a detailed geometric study [Bonalli et al., Possibly 2017] shows that neither abnormal nor singular arcs exist for this problem) while, as it is known, state constraints result to be more challenging (one could, for example, penalize these constraints within the cost function [Zhu et al., 2016]). Finally, as in previous sections, dynamics (2.1) reveals some singularities due to eulerian coordinates. In order to deal with this flaw, an inline change of coordinates can be considered during the numerical procedure. This strategy does not increase the computational effort and it allows to solve many more scenarios [Bonalli et al., Possibly 2017].

All these precious improvements will be reported in [Bonalli et al., Possibly 2017].