# Well-posedness of the shooting algorithm for control-affine problems with a scalar state constraint

Maria Soledad Aronna, Joseph Frédéric Bonnans, Bean San Goh

► To cite this version:

Maria Soledad Aronna, Joseph Frédéric Bonnans, Bean San Goh. Well-posedness of the shooting algorithm for control-affine problems with a scalar state constraint. Journal of Computational and Applied Mathematics, 2023, 10.1007/s40314-023-02346-2. http://dx.doi.org/10.1007/s40314-023-02346-2

HAL Id: hal-04134841

https://inria.hal.science/hal-04134841

Submitted on 28 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# WELL-POSEDNESS OF THE SHOOTING ALGORITHM FOR CONTROL-AFFINE PROBLEMS WITH A SCALAR STATE CONSTRAINT

M. SOLEDAD ARONNA AND J. FRÉDÉRIC BONNANS AND BEAN SAN GOH

ABSTRACT. We deal with a control-affine problem with scalar control subject to bounds, a scalar state constraint and endpoint constraints of equality type. For the numerical solution of this problem, we propose a shooting algorithm and provide a sufficient condition for its local convergence. We exhibit an example that illustrates the theory.

## 1. INTRODUCTION

In this article we deal with an optimal control problem governed by the dynamics

$$ \dot{x}_t = f_0(x_t) + u_t f_1(x_t) \quad \text{for a.a. } t \in [0, T], $$

subject to the control bounds

$$ u_{\min} \le u_t \le u_{\max}, $$

with $u_{\min} < u_{\max}$, endpoint constraints like

$$ \Phi(x_0, x_T) = 0, $$

and a scalar state constraint of the form

$$ g(x_t) \le 0. $$

For this class of problems, we propose a shooting-like numerical scheme and we show a sufficient condition for its local quadratic convergence, that is also a second order sufficient condition for optimality (in a particular sense to be specified later on). Additionally, we solve an example of practical interest, for which we also prove optimality by applying second order sufficient optimality conditions obtained in [3].

This investigation is strongly motivated by applications since we deal with both control and state constraints, which naturally appear in realistic models. Many practical examples that are covered by our chosen framework can be found in the existing literature. A non exhaustive list includes the prey-predator model [25], the Goddard problem in presence of a dynamic pressure limit [39, 26], the optimal control of the atmospheric arc for re-entry of a space shuttle seen in [12], an optimal production and maintenance system studied in [33], and a recent optimization problem on running strategies [2]. We refer also to [16], [31], [38] and references therein.

<footnote>The first author was supported by FAPERJ (Brazil) through the Jovem Cientista do Nosso Estado Program and by CNPq (Brazil) through the Universal Program and the Productivity Scholarship. The second author was supported by the FiME Lab Research Initiative (Institut Europlace de Finance), and by the PGMO program.</footnote>
As it is commonly known, the application of the necessary conditions provided by Pontryagin's Maximum Principle leads to an associated two-point boundary-value problem (TPBVP) for the optimal trajectory and its associated multiplier [42]. A natural way for solving TPBVPs numerically is the application of shooting algorithms [34]. This type of algorithms has been used extensively for the resolution of optimal control problems (see e.g. [13, 8, 36] and references therein). In particular, shooting methods have been applied to control-affine problems both with and without state constraints. Some works in this direction are mentioned in the sequel. Maurer [30] proposed a shooting scheme for solving a problem with bang-singular solutions, which was generalized quite recently by Aronna, Bonnans and Martinon in [5], where they provided a sufficient condition for its local convergence. Both these articles [30] and [5] analyze the case with control bounds and no state constraints. Practical control-affine problems with state constraints were solved numerically in several articles, a non extensive list includes Maurer and Gillessen [32], Oberle [35], Fraser-Andrews [21] and the recent articles Cots [14] and Cots et al [15]. Up to our knowledge, there is no result in the existing literature concerning sufficient conditions for the convergence of shooting algorithms in the framework considered here.

The paper is organized as follows. In Sections 2 and 2.2 we introduce the problem and give the basic definitions. A shooting-like method and a sufficient condition for its local quadratic convergence are given in Sections 3 and 4, respectively. The algorithm is implemented in Section 5 where we solve numerically a variation of the regulator problem and we prove the optimality of the solution analytically.

**Notations.** Let $\mathbb{R}^k$ denote the $k$-dimensional real space, i.e. the space of column real vectors of dimension $k$, and by $\mathbb{R}^{k*}$ its corresponding dual space, which consists of $k$-dimensional row real vectors. With $\mathbb{R}_+^k$ and $\mathbb{R}_-^k$ we refer to the subsets of $\mathbb{R}^k$ consisting of vectors with nonnegative, respectively nonpositive, components. We write $h_t$ for the value of function $h$ at time $t$ if $h$ is a function that depends only on $t$, and by $h_{i,t}$ the $i$th component of $h$ evaluated at $t$. Let $h(t+)$ and $h(t-)$ be, respectively, the right and left limits of $h$ at $t$, if they exist. Partial derivatives of a function $h$ of $(t, x)$ are referred as $D_t h$ or $\dot{h}$ for the derivative in time, and $D_x h$, $h_x$ or $h'$ for the differentiation with respect to space variables. The same convention is extended to higher order derivatives. By $L^p(0, T)^k$ we mean the Lebesgue space with domain equal to the interval $[0, T] \subset \mathbb{R}$ and with values in $\mathbb{R}^k$. The notations $W^{q,s}(0, T)^k$ and $H^1(0, T)^k$ refer to the Sobolev spaces (see Adams [1] for further details on Sobolev spaces). We let $BV(0, T)$ be the set of functions with bounded total variation. In general, when there is no place for confusion, we omit the argument $(0, T)$ when referring to a space of functions. For instance, we write $L^\infty$ for $L^\infty(0, T)$, or $(W^{1,\infty})^{k*}$ for the space of $W^{1,\infty}$-functions from $[0, T]$ to $\mathbb{R}^{k*}$. We say that a function $h: \mathbb{R}^k \to \mathbb{R}^d$ is of class $C^\ell$ if it is $\ell$-times continuously differentiable in its domain.

## 2. THE PROBLEM

Let us consider $L^\infty(0, T)$ and $W^{1,\infty}(0, T; \mathbb{R}^n)$ as control and state spaces, respectively. We say that a control-state pair $(u, x) \in L^\infty(0, T) \times W^{1,\infty}(0, T; \mathbb{R}^n)$ is a trajectory if it satisfies both the state equation

$$ (2.1) \qquad \dot{x}_t = f_0(x_t) + u_t f_1(x_t) \quad \text{for a.a. } t \in [0, T], $$
and the finitely many *endpoint constraints* of equality type given by

$$ (2.2) \qquad \Phi(x_0, x_T) = 0. $$

Here $f_0$ and $f_1$ are assumed to be Lipschitz continuous and twice continuously differentiable vector fields over $\mathbb{R}^n$, $\Phi$ is of class $C^2$ from $\mathbb{R}^n \times \mathbb{R}^n$ to $\mathbb{R}^q$. Under these hypotheses, for any pair control-initial condition $(u, x_0)$ in $L^\infty(0, T) \times \mathbb{R}^n$, the state equation (2.1) has a unique solution. Additionally, we consider a *cost functional*

$$ \phi : \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}, $$

the *control bounds*

$$ (2.3) \qquad u_{\min} \le u_t \le u_{\max} \quad \text{for a.a. } t \in [0, T], $$

where $u_{\min} < u_{\max}$, and a *scalar state constraint*

$$ (2.4) \qquad g(x_t) \le 0 \quad \text{for all } t \in [0, T], $$

with the functions $\phi$ and $g: \mathbb{R}^n \to \mathbb{R}$ being of class $C^2$. A trajectory $(u, x)$ is said to be *feasible* if it satisfies (2.3)-(2.4).

**Remark 2.1** (On the control bounds). We allow $u_{\min}$ and $u_{\max}$ to be either finite real numbers, or to take the values $-\infty$ or $+\infty$, respectively, meaning that we also consider problems with control constraints of the form $u_t \le u_{\max}$ or $u_{\min} \le u_t$, as well as problems in the absence of control constraints.

Summarizing, this article deals with the optimal control problem in the Mayer form given by

$$ (P) \qquad \min \phi(x_0, x_T); \quad \text{subject to} \ (2.1)-(2.4). $$

**2.1. Types of minima.** Throughout this article, we make use of two notions of optimality that are *weak* and *Pontryagin minima* and are defined as follows.

**Definition 2.2** (Weak and Pontryagin minima). A *weak minimum* for (P) is a feasible trajectory $(u, x)$ for which there exists $\varepsilon > 0$ such that $\phi(x_0, x_T) \le \phi(\tilde{x}_0, \tilde{x}_T)$ for any feasible $(\tilde{u}, \tilde{x})$ verifying $\|(\tilde{u}, \tilde{x}) - (u, x)\|_{\infty} < \varepsilon$.

A feasible trajectory $(u, x)$ is called a *Pontryagin minimum* for (P) if for any $M > 0$, there exists $\varepsilon_M > 0$ such that $\phi(x_0, x_T) \le \phi(\tilde{x}_0, \tilde{x}_T)$ for any feasible $(\tilde{u}, \tilde{x})$ satisfying

$$ (2.5) \qquad \|\tilde{x} - x\|_{\infty} + \|\tilde{u} - u\|_{1} < \varepsilon_M, \qquad \|\tilde{u} - u\|_{\infty} < M. $$

Note that any Pontryagin minimum is also a weak minimum. Consequently, necessary conditions that hold for weak minima, also do it for Pontryagin one. This article provides a numerical scheme for approximating Pontryagin minima of (P). In order to achieve this, we make use of the auxiliary unconstrained transformed problem (TP) given in equations (3.3)-(3.11), which possesses neither control bounds nor state constraints and can be solved numerically in an efficient way. In Lemma 3.2 below we prove that transformed Pontryagin minima of (P) that verify certain structural hypotheses are weak minima of the unconstrained transformed problem (TP).
2.2. **Bang, constrained and singular arcs.** The contact set associated with the state constraint is defined as

$$ (2.6) \qquad C := \{t \in [0, T] : g(\hat{x}_t) = 0\}. $$

For $0 \le a < b \le T$, we say that $(a, b)$ is an active arc for the state constraint or, shortly, a C arc, if $(a, b)$ is a maximal open interval contained in C. A point $\tau \in (0, T)$ is a junction point of the state constraint if it is the extreme point of a C arc.

Similar definitions hold for the control constraint, with the difference that the control variable is only almost everywhere defined. The contact sets for the control bounds are given by

$$ \begin{array}{l} B_- := \{t \in [0, T] : \hat{u}_t = u_{\min}\}, \quad B_+ := \{t \in [0, T] : \hat{u}_t = u_{\max}\}, \\ B := B_- \cup B_+. \end{array} $$

Note that these sets are defined up to null measure sets. Additionally, observe that if $u_{\min} = -\infty$ then $B_- = \emptyset$ and, analogously, if $u_{\max} = +\infty$ then $B_+ = \emptyset$. We say that $(a, b)$ is a $B_-$ (resp. $B_+$) arc if $(a, b)$ is included, up to a null measure set, in $B_-$ (resp. in $B_+$), but no open interval strictly containing $(a, b)$ is. We say that $(a, b)$ is a $B_+$ arc if it is either a $B_-$ or a $B_+$ arc.

Finally, let $S$ denote the singular set given by

$$ (2.7) \qquad S := \{t \in [0, T] : u_{\min} < \hat{u}_t < u_{\max} \text{ and } g(\hat{x}_t) < 0\}. $$

We say that $(a, b)$ is an $S$ arc if $(a, b)$ is included, up to a null measure set, in $S$, but no open interval strictly containing $(a, b)$ is.

We call junction or switching times the points $\tau \in (0, T)$ at which the trajectory $(\hat{x}, \hat{u})$ switches from one type of arc $(B_-, B_+, C)$ or $S$ to another type. Junction/switching times are denominated by the type of arcs they separate. One can have, for instance, CS junction, $B_-B_+$ switching time, etc.

Throughout the remainder of the article, we assume that the state constraint is of first order, this is,

$$ (2.8) \qquad g'(\hat{x}_t)f_1(\hat{x}_t) \neq 0 \quad \text{on } C, $$

and we impose the following hypotheses on the control structure:

$$ (2.9) \qquad \left\{ \begin{array}{ll} (i) & \text{the interval } [0, T] \text{ is (up to a zero measure set) the disjoint union of finitely many arcs of type } B, C \text{ and } S, \text{ and the set } C \text{ does not contain isolated points,} \\ (ii) & \text{the control } \hat{u} \text{ is at uniformly positive distance of the bounds } u_{\min} \text{ and } u_{\max}, \text{ over } C \text{ and } S \text{ arcs,} \\ (iii) & \text{the control } \hat{u} \text{ is discontinuous at CS and SC junctions.} \end{array} \right. $$

*Remark 2.3.* Note that some problems (even if they are convex, like Fuller’s problem [22]) exhibit chattering phenomena with infinitely many switches, necessarily with some very short arcs. It is not clear how to deal with such problems with the method that we present in this article.

The example of the regulator problem studied in Section 5 fullfils the above hypothesis (2.9) (see as well the example given in [3, Remark 2]).

When a control satisfying hypothesis (2.9)(i) is, for instance, a concatenation of a bang and a singular arc, we call it a BS control. This denomination is extended for any finite sequence of arc types.
In order to formulate our shooting algorithm, we express the control as a function of the state on C arcs and we fix the control to its bounds on B arcs.

2.2.1. Expression of the control on constrained arcs. From $g(\hat{x}_t) = 0$ on C, we get

$$ (2.10) \qquad 0 = \frac{d}{dt} g(\hat{x}_t) = g'(\hat{x}_t)(f_0(\hat{x}_t) + \hat{u}_t f_1(\hat{x}_t)) \quad \text{on } C, $$

and, since (2.8) holds, we have that

$$ (2.11) \qquad \hat{u}_t = - \frac{g'(\hat{x}_t) f_0(\hat{x}_t)}{g'(\hat{x}_t) f_1(\hat{x}_t)} \quad \text{on } C. $$

### 3. SHOOTING FORMULATION

We now explain how to state a transformed problem with neither control bounds nor running state constraints that serves as an intermediate step to write a numerical scheme for problem (P). Afterwards the optimality system of the transformed problem is reduced to a nonlinear equation in a finite dimensional space.

The starting point is to estimate the arc structure of the control, i.e. the sequence of its different types of arcs and the approximate values of its junction times. This is done in practice by some direct method such as solving the nonlinear programming (NLP) associated to the discretization of the optimal control problem. Then we formulate a transformed problem in which the control is fixed to its bounds on B arcs, and is expressed as a function of the state on C arcs. So the optimization variables are now the control over singular arcs and the switching times. Subsequently, we express the optimality conditions of the transformed problem. Finally, by eliminating the control as a function of the state and costate, we reduce the optimality system to a finite dimensional equation.

So, let us assume for the remainder of the section that $(\hat{u}, \hat{x})$ is a Pontryagin minimum for (P). Additionally, without loss of generality and for the sake of simplicity of notation, we set $u_{\min} := 0$ and $u_{\max} := 1$. Recall further that $(\hat{u}, \hat{x})$ complies with the structural hypotheses (2.9) for the control $\hat{u}$, and that the state constraint is of first order, i.e. (2.8) holds true.

3.1. The transformed problem. We now state the transformed problem corresponding to (P), in the spirit of [5], and we prove that any Pontryagin minimum for the original problem (P) is transformed into a weak minimum of the unconstrained transformed problem.

For the Pontryagin minimum $(\hat{u}, \hat{x})$, let

$$ (3.1) \qquad 0 := \hat{\tau}_0 < \hat{\tau}_1 < \dots < \hat{\tau}_N := T $$

denote its associated switching times. Recall the definition of the sets C, $B_-$, $B_+$ and S given in Section 2.2 above. Set $\hat{I}_k := [\hat{\tau}_{k-1}, \hat{\tau}_k]$ for $k = 1, \dots, N$, and

$$ (3.2) \qquad \mathcal{I}(S) := \{k \in \{1, \dots, N\} : \hat{I}_k \text{ is a singular arc}\}. $$

Analogously, define $\mathcal{I}(C)$, $\mathcal{I}(B_-)$, and $\mathcal{I}(B_+)$. For each $k = 1, \dots, N$, consider a state variable $x^k \in W^{1,\infty}(0, T; \mathbb{R}^n)$, and for each singular arc $k \in \mathcal{I}(S)$, a control variable $u^k \in L^\infty(0, T)$. On the set B, we fix the control to the corresponding bound. Additionally, recall that from formula (2.11) we have that $\hat{u}_t = \Gamma(\hat{x}_t)$ on C, where $\Gamma$ is given by

$$ \Gamma(x) := - \frac{g'(x) f_0(x)}{g'(x) f_1(x)}. $$
After these considerations, we are ready to state the transformed problem. Define the optimal control problem (TP), on the time interval $[0, 1]$, by

(3.3) \qquad $\min \phi(x_0^1, x_1^N),$

$$
(3.4) \qquad \dot{x}^k = (\tau_k - \tau_{k-1})(f_0(x^k) + u^k f_1(x^k)) \quad \text{for } k \in \mathcal{I}(S),
$$

$$
(3.5) \qquad \dot{x}^k = (\tau_k - \tau_{k-1}) f_0(x^k) \quad \text{for } k \in \mathcal{I}(B_-),
$$

$$
(3.6) \qquad \dot{x}^k = (\tau_k - \tau_{k-1}) (f_0(x^k) + f_1(x^k)) \quad \text{for } k \in \mathcal{I}(B_+),
$$

$$
(3.7) \qquad \dot{x}^k = (\tau_k - \tau_{k-1}) (f_0(x^k) + \Gamma(x^k)f_1(x^k)) \quad \text{for } k \in \mathcal{I}(C),
$$

$$
(3.8) \qquad \dot{\tau}_k = 0 \quad \text{for } k = 1, \dots, N-1,
$$

$$
(3.9) \qquad \Phi(x_0^1, x_1^N) = 0,
$$

$$
(3.10) \qquad g(x_0^k) = 0 \quad \text{for } k \in \mathcal{I}(C),
$$

$$
(3.11) \qquad x_1^k = x_0^{k+1} \quad \text{for } k = 1, \dots, N-1.
$$

Remark 3.1. Since we use the expression (2.11) to eliminate the control from the expression (2.10) of the derivative of the state constraint equal to zero, we impose the entry conditions (3.10) in the formulation of (TP) in order to guarantee that the state constraint is active along $x^k$ for every $k \in \mathcal{I}(C)$.

Set

$$
(3.12) \qquad \begin{aligned}
\hat{x}_s^k &:= \hat{x}(\hat{\tau}_{k-1} + (\hat{\tau}_k - \hat{\tau}_{k-1})s) && \text{for } s \in [0,1] \text{ and } k=1,\dots,N, \\
\hat{u}_s^k &:= \hat{u}(\hat{\tau}_{k-1} + (\hat{\tau}_k - \hat{\tau}_{k-1})s) && \text{for } s \in [0,1] \text{ and } k \in \mathcal{I}(S).
\end{aligned}
$$

**Lemma 3.2.** Let $(\hat{u}, \hat{x})$ be a Pontryagin minimum of problem (P). Then the triple

$$
((\hat{u}^k)_{k \in \mathcal{I}(S)}, (\hat{x}^k)_{k=1}^N, (\hat{\tau}_k)_{k=1}^{N-1})
$$

is a weak solution of (TP).

*Proof.* Consider the feasible trajectories $((u^k), (x^k), (\tau_k))$ for (TP) satisfying

$$
(3.13) \qquad \|u^k - \hat{u}^k\|_{\infty} < \bar{\epsilon} \quad \text{and} \quad |\tau_k - \hat{\tau}_k| \le \bar{\delta} \quad \text{for all } k = 1, \dots, N,
$$

for some $\bar{\epsilon}, \bar{\delta} > 0$ to be determined later. Set $I_k := [\tau_{k-1}, \tau_k]$ and consider the functions $s_k: I_k \to [0, 1]$ given by $s_{k,t} := \frac{t - \tau_{k-1}}{\tau_k - \tau_{k-1}}$. Define $u: [0, T] \to \mathbb{R}$ by

$$
(3.14) \qquad u_t := \begin{cases} 0 & \text{if } t \in I_k, k \in \mathcal{I}(B_-), \\ 1 & \text{if } t \in I_k, k \in \mathcal{I}(B_+), \\ \Gamma(x^k(s_{k,t})) & \text{if } t \in I_k, k \in \mathcal{I}(C), \\ u^k(s_{k,t}) & \text{if } t \in I_k, k \in \mathcal{I}(S). \end{cases}
$$

Let $x: [0, T] \to \mathbb{R}^n$ be the state corresponding to the control $u$ and the initial condition $x(0) = x_0^1$.

We next show that if $\bar{\epsilon} > 0$ and $\bar{\delta} > 0$ are small enough, then $(u, x)$ is feasible for (P) and arbitrarily close to $(\hat{u}, \hat{x})$ in the sense of (2.5). Observe that $x(t) = x^k(s_{k,t})$ for all $k = 1, \dots, N$ and $t \in I_k$. Hence, $x$ satisfies the endpoint constraints. Furthermore, due to Gronwall's Lemma, $(u, x)$ verifies the estimate

$$
(3.15) \qquad \|u - \hat{u}\|_1 + \|x - \hat{x}\|_{\infty} < \mathcal{O}(\bar{\epsilon} + \bar{\delta}).
$$
Let us analyze the control constraints. Take $k = 1, \dots, N$. If $k \in \mathcal{I}(B_-) \cup \mathcal{I}(B_+)$, then $u_t \in \{0, 1\}$ for a.a. $t \in I_k$. On the other hand, by the hypothesis (2.9) on the control structure, there exists $\rho_1 > 0$ such that

$$ (3.16) \qquad \rho_1 < \hat{u}_t < 1 - \rho_1 \quad \text{over } C \text{ and } S \text{ arcs.} $$

Suppose now that $k \in \mathcal{I}(S)$. Then, in view of (3.13) and (3.16), we can see that the control constraints hold on $I_k$ provided that $\bar{\epsilon} \le \rho_1$. Finally, let $k$ be in $\mathcal{I}(C)$. Notice that in this case (3.16) is equivalent to

$$ (3.17) \qquad \rho_1 < \Gamma(\hat{x}_t) < 1 - \rho_1 \quad \text{on } \hat{I}_k. $$

Hence, by standard continuity arguments and for $\bar{\epsilon}, \bar{\delta}$ sufficiently small, we get that

$$ (3.18) \qquad 0 < \Gamma(x_t) < 1 \quad \text{on } I_k. $$

We therefore confirm that $(u, x)$ verifies the control constraints.

Let us now consider the state constraint. Take first $k \in \mathcal{I}(C)$. Then $g(x_{\tau_k}) = g(x_0^k) = 0$ and, by definition of $(u, x)$, we have that $\frac{d}{dt}g(x_t) = 0$ for all $t \in I_k$. Therefore, $x$ satisfies the state constraint on $I_k$ for $k \in \mathcal{I}(C)$. Next, observe that, due to (2.9), for any $t \in [0, T]$ sufficiently far from a $C$ arc, one has $g(\hat{x}_t) \le -\rho$ for some small $\rho > 0$. Thus, by (3.15) we get that $g(x_t) < 0$ for appropriate $\bar{\epsilon}, \bar{\delta}$. On the other hand, for $t \in [0, T]$ close to a $C$ arc, we reason as follows. Assume, without loss of generality, that $t$ is near an entry point $\tau_k$ of a $C$ arc. In view of hypothesis (2.9) and of the relation (2.10), we have that $\frac{d}{ds}|_{s=\tau_k}g(\hat{x}_s) > 0$, therefore, $\frac{d}{ds}|_{s=\tau_k}g(x_s) > 0$ as well, if $\bar{\epsilon}, \bar{\delta}$ are sufficiently small. Consequently, $g(x_t) < 0$. Hence, $x$ verifies the state constraint on $[0, T]$. With this, we conclude that $(u, x)$ is feasible for the original problem (P).

Finally, given $M > 0$ as in Definition 2.2, we can easily show that $\bar{\delta}$ and $\bar{\epsilon}$ can be taken in such a way that $(u, x)$ satisfies (2.5) for such $M$ and the corresponding $\epsilon_M$ provided by the Pontryagin optimality of $(\hat{u}, \hat{x})$. Consequently, $\phi(x_0, x_1) \ge \phi(\hat{x}_0, \hat{x}_1)$ or, equivalently,

$$ (3.19) \qquad \phi(x_0^1, x_1^N) \ge \phi(\hat{x}_0^1, \hat{x}_1^N), $$

which proves that $((\hat{u}^k)_{k \in \mathcal{I}(S)}, (\hat{x}^k)_{k=1}^N, (\hat{\tau}_k)_{k=1}^{N-1})$ is a weak solution of (TP), as desired. This concludes the proof. $\square$

3.2. **The shooting function.** We shall start by rewriting the problem (TP) in the following compact form, in order to ease the notation,

$$ (3.20) \qquad \min \tilde{\phi}(X_0, X_1), $$

$$ (3.21) \qquad \dot{X} = \tilde{f}_0(X) + \sum_{k \in \mathcal{I}(S)} U^k \tilde{f}_k(X), $$

$$ (3.22) \qquad \tilde{\Phi}(X_0, X_1) = 0, $$

where $X := ((x^k)_{k=1}^N, (\tau_k)_{k=1}^{N-1})$, $U := (u^k)_{k \in \mathcal{I}(S)}$, the vector field $\tilde{f}_0 : \mathbb{R}^{Nn+N-1} \to \mathbb{R}^{Nn+N-1}$ is defined as follows,

$$ (\tilde{f}_0(X))_{i=(k-1)n+1}^{kn} \quad \begin{cases} (\tau_k - \tau_{k-1})f_0(x^k) & \text{for } k \in \mathcal{I}(S) \cup \mathcal{I}(B_-), \\ (\tau_k - \tau_{k-1})(f_0(x^k) + f_1(x^k)) & \text{for } k \in \mathcal{I}(B_+), \\ (\tau_k - \tau_{k-1})(f_0(x^k) + \Gamma(x^k)f_1(x^k)) & \text{for } k \in \mathcal{I}(C), \end{cases} $$
and $(\tilde{f}_0(X))_{i=nN+1}^{Nn+N-1} := 0$. Additionally, for $k \in \mathcal{I}(S)$ the vector field $\tilde{f}_k : \mathbb{R}^{Nn+N-1} \to \mathbb{R}^{Nn+N-1}$ is given by

$$ (\tilde{f}_k(X))_{i=(k-1)n+1}^{kn} := (\tau_k - \tau_{k-1})f_1(x^k), $$

and $(\tilde{f}_k(X))_i := 0$ for the remaining index $i$, the new cost $\tilde{\phi} : \mathbb{R}^{2(Nn+N-1)} \to \mathbb{R}$ is

$$ \tilde{\phi}(X_0, X_1) := \phi(x_0^1, x_1^N), $$

and the function $\tilde{\Phi} : \mathbb{R}^{2(Nn+N-1)} \to \mathbb{R}^{d_{\Phi}}$ with $d_{\Phi} := q + |\mathcal{I}(C)| + n(N-1)$ is defined as

$$ \tilde{\Phi}(X_0, X_1) := \begin{pmatrix} \Phi(x_0^1, x_1^N) \\ (g(x_0^k))_{k \in \mathcal{I}(C)} \\ (x_1^k - x_0^{k+1})_{k=1}^{N-1} \end{pmatrix}. $$

The pre-Hamiltonian for problem (TP) is given by

$$ (3.23) \qquad \tilde{H} = P(\tilde{f}_0(X) + \sum_{k \in \mathcal{I}(S)} U^k \tilde{f}_k(X)) = \sum_{k=1}^{N} (\tau_k - \tau_{k-1}) H^k, $$

where $P$ denotes the costate associated to (TP),

$$ (3.24) \qquad H^k := p^k(f_0(x^k) + w^k f_1(x^k)), $$

with the notation $w^k$ defined as

$$ (3.25) \qquad w^k := \begin{cases} u^k & \text{if } k \in \mathcal{I}(S), \\ 0 & \text{if } k \in \mathcal{I}(B_-), \\ 1 & \text{if } k \in \mathcal{I}(B_+), \\ \Gamma(x^k) & \text{if } k \in \mathcal{I}(C), \end{cases} $$

and $p^k$ denotes the $n$-dimensional vector of components $P_{(k-1)n+1}, \dots, P_{kn}$. Note that $w^k$ is a variable only for $k \in \mathcal{I}(S)$, in which case it represents the control $u^k$.

### 3.3. Constraint qualification and first order optimality condition for (TP).

Since problem (TP) has only endpoint equality constraints, and the Hamiltonian is an affine function of the control, it is known that Pontryagin's Maximum Principle is equivalent to the first-order optimality conditions. So, the qualification condition is that the derivative of the constraint is onto at the nominal trajectory $(\hat{U}, \hat{X})$, see e.g. [11, Ch. 3]. This means that

$$ (3.26) \qquad \tilde{\Phi} : \mathbb{R}^{Nn+N-1} \times (L^\infty)^{|\mathcal{I}(S)|} \to \mathbb{R}^{d_{\tilde{\Phi}}}, \quad (X_0, U) \mapsto \tilde{\Phi}(X_0, X_1), $$

where $X_t$ is the solution of (3.21) associated to $(X_0, U)$ is such that

$$ (3.27) \qquad D\tilde{\Phi}(\hat{X}_0, \hat{U}) \text{ is surjective.} $$

Under this hypothesis, the first-order optimality condition in normal form is as follows, defining the endpoint Lagrangian associated to (TP) by:

$$ (3.28) \qquad \tilde{\ell}^\Psi := \phi(x_0^1, x_1^N) + \sum_{j=1}^{q} \Psi_j \Phi_j(x_0^1, x_1^N) + \sum_{k \in \mathcal{I}(C)} \gamma_k g(x_0^k) + \sum_{k=1}^{N-1} \theta_k(x_1^k - x_0^{k+1}). $$
**Theorem 3.3.** Let $(\hat{U}, \hat{X})$ be a weak solution for (TP) satisfying the qualification condition (3.27). Then there exists a unique $\tilde{\lambda} := (\tilde{\Psi}, P) \in \mathbb{R}^{d_{\Phi}^*} \times (W^{1,\infty})^{Nn+N-1}$ such that P is solution of

$$ (3.29) \qquad - \dot{P}_t = D_X \tilde{H}(\hat{U}_t, \hat{X}_t, P_t) \quad \text{a.e. on } [0, T], $$

with transversality conditions

$$ (3.30) \qquad \begin{aligned} P_0 &= -D_{X_0} \tilde{\ell}^{\tilde{\Psi}}(\hat{X}_0, \hat{X}_1), \\ P_1 &= D_{X_T} \tilde{\ell}^{\tilde{\Psi}}(\hat{X}_0, \hat{X}_1), \end{aligned} $$

and with

$$ (3.31) \qquad \tilde{H}_U(\hat{U}_t, \hat{X}_t, P_t) = 0. $$

*Proof.* This is a variant of [9, Theorem 1.174], where the cost function was supposed to be convex: one easily checks that the proof is essentially the same if the cost is differentiable. See also [6]. □

Since there is a unique associated multiplier, we omit from now on the dependence on $\tilde{\lambda}$ for the sake of simplicity of the presentation. Moreover, in some occasions, we omit the dependence on the nominal solution $(\hat{U}, \hat{X})$.

**3.4. Expression of the singular controls in problem (TP).** It is known that in this control-affine case, the control variable does not appear explicitly neither in the expression of $\tilde{H}_U$ nor in its time derivative $\dot{\tilde{H}}_U$ (see e.g. [37, 5]). The strengthened generalized Legendre-Clebsch condition [37] for (TP) reads

$$ (3.32) \qquad - \frac{\partial}{\partial U} \ddot{\tilde{H}}_U \succ 0. $$

Here $A \succ B$, where A and B are symmetric matrices of same size, means that $A-B$ is positive semidefinite. At this point, recall the definitions of $H^k$ and $p^k$ given in (3.24) and in the first line after (3.25), respectively. Simple calculations show that the l.h.s. of (3.32) is a $|\mathcal{I}(S)| \times |\mathcal{I}(S)|$-diagonal matrix with positive entries equal to

$$ (3.33) \qquad - (\tau_k - \tau_{k-1}) \frac{\partial}{\partial u^k} \ddot{H}_{u^k}^k \quad \text{for } k \in \mathcal{I}(S). $$

Then condition (3.32) becomes

$$ (3.34) \qquad \frac{\partial}{\partial u^k} \ddot{H}_{u^k}^k < 0 \quad \text{for } k \in \mathcal{I}(S). $$

Hence, thanks to (3.34), for each $k \in \mathcal{I}(S)$ one can compute the control $u^k$ from the identity

$$ (3.35) \qquad \ddot{H}_{u^k}^k = 0. $$

Apart from the previous equation (3.35), in order to ensure the stationarity $H_{u^k}^k = 0$, we add the following endpoint conditions:

$$ (3.36) \qquad 0 = H_{u^k}^k(0) = p_0^k f_1(x_0^k), \quad 0 = \dot{H}_{u^k}^k(0) = p_0^k[f_1, f_0](x_0^k) \quad \text{for } k \in \mathcal{I}(S). $$
3.5. **Lagrangians and costate equation.** The costate equation for $p^k$ is

$$ (3.37) \qquad \dot{p}^k = -(\tau_k - \tau_{k+1})D_{x^k}H^k, $$

with endpoint conditions

$$ (3.38) \qquad p_0^1 = -D_{x_0^1} \ell^\Psi = -D_{x_0^1} \phi - \sum_{j=1}^q \Psi_j D_{x_0^1} \Phi_j - \chi_{\mathcal{I}(C)}(1) \gamma_1 g'(x_0^1), $$

$$ (3.39) \qquad p_1^k = \theta^k \quad \text{for } k = 1, \dots, N-1, $$

$$ (3.40) \qquad p_0^k = \theta^{k-1} - \chi_{\mathcal{I}(C)}(k) \gamma_k g'(x_0^k) \quad \text{for } k = 2, \dots, N, $$

$$ (3.41) \qquad p_1^N = D_{x_1^N} \phi + \sum_{j=1}^q \Psi_j D_{x_1^N} \Phi_j, $$

where $\chi_{\mathcal{I}(C)}$ denotes the characteristic function associated to the set $\mathcal{I}(C)$. For the costate $p^{\tau_k}$ we have the dynamics

$$ (3.42) \qquad \dot{p}^{\tau_k} = -H^k + H^{k+1}, \quad p_0^{\tau_k} = 0, \quad p_1^{\tau_k} = 0 \quad \text{for } k = 1, \dots, N-1. $$

It is known that the pre-Hamiltonian of autonomous problems has a constant value along an optimal solution (see e.g. [42]). By similar arguments it is easily seen that each $H^k$ is a constant function of time along an optimal solution. Consequently, from (3.42) we get that $p^{\tau_k}$ vanishes identically and that

$$ (3.43) \qquad H_1^k = H_0^{k+1} \quad \text{for } k = 1, \dots, N-1. $$

3.6. **The shooting function and method.** The shooting function associated with (TP) that we propose here is

$$ (3.44) \qquad S: \mathbb{R}^{Nn+N-1} \times \mathbb{R}^{Nn+q+|\mathcal{I}(C)|,*} \to \mathbb{R}^{(N-1)n+N-1+q+|\mathcal{I}(C)|+2|\mathcal{I}(S)|} \times \mathbb{R}^{(N+1)n,*}, $$

$$ S: \mathbb{R}^{Nn+N-1} \times \mathbb{R}^{Nn+q+|\mathcal{I}(C)|,*} \to \mathbb{R}^{(N-1)n+N-1+q+|\mathcal{I}(C)|+2|\mathcal{I}(S)|} \times \mathbb{R}^{(N+1)n,*}, $$

$$ ((x_0^k), (\tau_k), (p_0^k), \Psi, \gamma) \mapsto \left( \begin{array}{l} \Phi(x_0^1, x_1^N) \\ (g(x_0^k))_{k \in \mathcal{I}(C)} \\ (x_1^k - x_0^{k+1})_{k=1, \dots, N-1} \\ p_0^1 + D_{x_0^1} \ell^\Psi \\ p_1^k - p_0^{k+1} - \chi_{\mathcal{I}(C)}(k) \gamma_k g'(x_0^k) \\ p^q - D_{x_0^N} \ell^\Psi \\ (H_1^k - H_0^{k+1})_{k=1, \dots, N-1} \\ (p_0^k f_1(x_0^k))_{k \in \mathcal{I}(S)} \\ (p_0^k[f_1, f_0](x_0^k))_{k \in \mathcal{I}(S)} \end{array} \right), $$

where $((x^k), (p^k))$ is the solution of the state and costate equations (3.4)-(3.7), (3.37) with initial values $(x_0^k)$, $(p_0^k)$, and control $(u^k)_{k \in \mathcal{I}(S)}$, given by the stationarity condition (3.35). Note that we removed the variable $\theta$ by combining equations (3.39) and (3.40).

The key feature of this procedure is that $\omega := ((x_0^k), (\tau_k), (p_0^k), \Psi, \gamma)$ satisfies

$$ (3.45) \qquad \mathcal{S}(\omega) = 0, $$

if and only if the associated solution $((x^k), (p^k), (u^k))$ verifies the Pontryagin's Maximum Principle for (TP). Briefly speaking, in order to find the candidate solutions of (TP), we shall solve (3.45).
Let us observe that the system (3.45) has $2Nn + N - 1 + q + |\mathcal{I}(C)|$ unknowns and $2Nn + N - 1 + q + |\mathcal{I}(C)| + 2|\mathcal{I}(S)|$ equations. Hence, as soon as a singular arc occurs, (3.45) has more equations than unknowns, i.e. it is overdetermined. We then follow [5], where the authors suggested to solve the shooting equations by the Gauss-Newton method. We recall the following convergence result for Gauss-Newton, see e.g. Fletcher [19], or alternatively Bonnans [10]. If $F$ is a $C^1$ mapping from $\mathbb{R}^n$ to $\mathbb{R}^p$ with $p > n$, the Gauss-Newton method computes a sequence $(y^j)$ in $\mathbb{R}^n$ satisfying $F(y^j) + DF(y^j)(y^{j+1} - y^j) = 0$. When $F$ has a zero at $\bar{y}$ and $DF(\bar{y})$ is onto, the sequence $(y^j)$ is well-defined provided that the starting point $y^0$ is close enough to $\bar{y}$ and in that case, $(y^j)$ converges superlinearly to $\bar{y}$ (quadratically if $DF$ is Lipschitz near $\bar{y}$). In view of the regularity hypotheses done in Section 2, we know that $S'$ is Lipschitz continuous.

## 4. SUFFICIENT CONDITION FOR THE CONVERGENCE OF THE SHOOTING ALGORITHM

The main result of this article is Theorem 4.5 of current section. It gives a sufficient condition for the local convergence of the shooting algorithm, that is also a sufficient condition for weak optimality of problem (TP), as stated in Theorem 4.4 below.

### 4.1. Second order optimality conditions for (TP).
We now recall some second order optimality conditions for (TP). Let us consider the quadratic mapping on the space $(L^\infty)^{|\mathcal{I}(S)|} \times (W^{1,\infty})^{Nn+N-1}$, defined as

$$ (4.1) \qquad \tilde{Q}(V, Z) := \frac{1}{2} D^2 \tilde{\ell}(Z_0, Z_1)^2 + \frac{1}{2} \int_0^1 [Z^\top \tilde{H}_{XX} Z + 2V \tilde{H}_{UX} Z] dt. $$

We next introduce the critical cone associated to (TP). Since the problem has only qualified equality constraints, this critical cone coincides with the tangent space to the constraints. Consider first the linearized state equation

$$ (4.2) \qquad \dot{Z} = \tilde{A}Z + \tilde{B}V \quad \text{a.e. on } [0, 1], $$

where $F(U, X) := \tilde{f}_0(X) + \sum_{k \in \mathcal{I}(S)} U^k \tilde{f}_k(X)$, $\tilde{A} := F_X$, $\tilde{B} := F_U$; and let the linearization of the endpoint constraints be given by

$$ (4.3) \qquad D\tilde{\Phi}(Z_0, Z_1) = 0. $$

The critical cone for (TP) is defined as

$$ (4.4) \qquad \tilde{C} := \{ (V, Z) \in (L^\infty)^{|\mathcal{I}(S)|} \times (W^{1,\infty})^{Nn+N-1} : (4.2)-(4.3) \text{ hold} \}. $$

The following result follows (see e.g. [29, 4] for a proof).

**Theorem 4.1 (Second order necessary condition).** If $(\hat{U}, \hat{X})$ is a weak minimum for (TP) that verifies (3.27), then

$$ (4.5) \qquad \tilde{Q}(V, Z) \ge 0 \quad \text{for all } (V, Z) \in \tilde{C}. $$

In the sequel we present some optimality conditions for (TP). The first one is a necessary condition due to Goh [23] and the second one, a sufficient condition from Dmitruk [17, 18]. The idea behind these results lies on the following observation. Note that $\tilde{H}_{UU}$ vanishes and, therefore, the quadratic mapping $\tilde{Q}$ does not contain a quadratic term on the control variation $V$. Consequently, the Legendre-Clebsch
necessary optimality condition on the positive semidefiniteness of $\tilde{H}_{UU}$ holds trivially and a second order sufficient condition cannot be obtained by strengthening inequality (4.5). In order to overcome this issue and derive necessary conditions for this singular case, Goh introduced a change of variables in [24] and applied it to derive necessary conditions in [23]. Some years later, Dmitruk [17] showed a second order sufficient condition in terms of the coercivity of the transformation $\tilde{\Omega}$ of $\tilde{Q}$ introduced below.

The Goh transformation for the linear system (4.2) is given by

$$ (4.6) \qquad Y_t := \int_0^t V_s ds, \quad \Xi_t := Z_t - \tilde{B}_t Y_t. $$

Notice that if $(V, Z) \in \tilde{C}$, then $(Y, \Xi)$ defined by the above transformation (4.6) satisfies (removing time indexes):

$$ (4.7) \qquad \dot{\Xi} = AZ + \tilde{B}U - \tilde{B}U - \dot{\tilde{B}}Y = A(\Xi + \tilde{B}Y) - \dot{\tilde{B}}Y, $$

and, therefore $\Xi$ is solution of the transformed linearized equation

$$ (4.8) \qquad \dot{\Xi} = \tilde{A}\Xi + \tilde{E}Y, $$

where $\tilde{E} := \tilde{A}\tilde{B} - \frac{d}{dt}\tilde{B}$, and $\Xi$ satisfies the transformed linearized endpoint constraints

$$ (4.9) \qquad D\tilde{\Phi}(\Xi_0, \Xi_1 + \tilde{B}_1h) = 0, $$

where we set $h := Y_1$.

Consider the function

$$ (4.10) \qquad \rho(\zeta_0, \zeta_1, h) := D^2\tilde{\ell}(\zeta_0, \zeta_1 + \tilde{B}_1h)^2 + h\tilde{H}_{UX,1}(2\zeta_1 + \tilde{B}_1h), $$

and the quadratic mapping

$$ (4.11) \qquad \tilde{\Omega}(Y, \tilde{h}, \Xi) := \frac{1}{2}\rho(\Xi_0, \Xi_1, \tilde{h}) + \frac{1}{2} \int_0^1 \left( \Xi^\top \tilde{H}_{XX} \Xi + 2Y \tilde{M} \Xi + Y \tilde{R} Y \right) dt, $$

for $(Y, \tilde{h}, \Xi) \in (L^2)^{|\mathcal{I}(S)|} \times \mathbb{R} \times (H^1)^{Nn+N-1}$ and

$$ (4.12) \qquad \tilde{M} := \tilde{f}_1^\top \tilde{H}_{XX} - \frac{d}{dt} \tilde{H}_{UX} - \tilde{H}_{UX} \tilde{A}, \quad \tilde{R} := \tilde{f}_1^\top \tilde{H}_{XX} \tilde{f}_1 - 2 \tilde{H}_{UX} \tilde{E} - \frac{d}{dt} (\tilde{H}_{UX} \tilde{f}_1). $$

Let us recall that the second order necessary condition for optimality stated by Goh [23] (and nowadays known as Goh condition) implies that if $(\hat{U}, \hat{X})$ is a weak minimum for (TP) verifying (3.27), then

$$ (4.13) \qquad \tilde{H}_{UX} \tilde{B} \text{ is symmetric,} $$

or, equivalently, $P \cdot D\tilde{f}_i \tilde{f}_j = P \cdot D\tilde{f}_j \tilde{f}_i$ for all $i, j = 1, \dots, N$. In the recent literature, this condition can be encountered as $P \cdot [\tilde{f}_i, \tilde{f}_j] = 0$. We shall mention that this necessary condition was first stated by Goh in [23] for the case with neither control nor state constraints, and extended in [4, 20] for problems containing control constraints.

Notice that when the control variable of (TP) is scalar (i.e. $|\mathcal{I}(S)| = 1$), then (4.13) is trivially verified since $\tilde{H}_{UX} \tilde{B}$ is also a scalar.

We claim that $(Df_i)f_j = 0$, for all $i, j \in \{1, \dots, N\}$ with $i \neq j$ and, consequently, the matrix in (68) is diagonal and the Goh condition holds trivially for (TP) even when $\mathcal{I}(S) > 1$. Indeed, let $k, \ell$ be disjoint elements of $\mathcal{I}(S)$. By the definition, the
nonzero components of $\tilde{f}_k(X)$ have coordinates in the set $I_k := \{(k-1)n+1, kn\}$. So, only the rows of $D\tilde{f}_k(X)$ in $I_k$ can be nonzero. However, since $I_k \cap I_\ell$ is empty, $\tilde{f}_\ell$ has only zero components in $I_k$. Our claim follows.

From (4.13) and [4, Theorem 4.4] we get the following result:

**Proposition 4.2.** For all $(V, Z) \in (L^\infty)^{|\mathcal{I}(S)|} \times (W^{1,\infty})^{Nn+N-1}$ solution of (4.2) and $(Y, \Xi)$ given by the Goh transform (4.6), it holds

$$ \tilde{Q}(V, Z) = \tilde{\Omega}(Y, Y_T, \Xi). $$

Define, for $(\Xi_0, Y, \tilde{h}) \in \mathbb{R}^n \times (L^2)^{|\mathcal{I}(S)|} \times \mathbb{R}$, the order function

$$ (4.14) \qquad \gamma(\Xi_0, Y, \tilde{h}) := |\Xi_0|^2 + \int_0^1 |Y_t|^2 dt + |\tilde{h}|^2. $$

**Definition 4.3** $(\gamma$-growth). A feasible trajectory $(\hat{U}, \hat{X})$ of (TP) satisfies the $\gamma$-growth condition in the weak sense if there exists a positive constant $c$ such that, for every sequence of feasible variations $\{(\delta X_0^k, V^k)\}_k$ converging to 0 in $\mathbb{R}^{Nn+N-1} \times (L^\infty)^{|\mathcal{I}(S)|}$, one has that

$$ (4.15) \qquad \tilde{\phi}(X_0^k, X_1^k) - \tilde{\phi}(\hat{X}_0, \hat{X}_1) \geq c\gamma(\Xi_0^k, Y^k, Y_T^k), $$

for $k$ large enough, where $(Y^k, \Xi^k)$ are given by Goh transform (4.6) and $X^k$ is the solution of the state equation (3.21) associated to $(\hat{X}_0 + \delta X_0^k, \hat{U} + V^k)$.

Consider the transformed critical cone

$$ (4.16) \qquad \tilde{P}_S^2 := \{ (Y, \tilde{h}, \Xi) \in (L^2)^{|\mathcal{I}(S)|} \times \mathbb{R} \times (H^1)^{Nn+N-1} : (4.8)-(4.9) \text{ hold} \}. $$

The following characterization of $\gamma$-growth holds (see [17] or [18, Theorem 3.1] for a proof).

**Theorem 4.4.** Let $(\hat{U}, \hat{X})$ be such that the qualification condition (3.27) holds. Then $(\hat{U}, \hat{X})$ is a weak minimum of (TP) that satisfies $\gamma$-growth in the weak sense if and only if (4.13) holds and there exists $c > 0$ such that

$$ (4.17) \qquad \tilde{\Omega}(Y, \tilde{h}, \Xi) \geq c\gamma(\Xi_0, Y, \tilde{h}) \quad \text{on } \tilde{P}_S^2. $$

We are now ready to state the following convergence result for the shooting algorithm.

**Theorem 4.5.** If $(\hat{U}, \hat{X})$ is a weak minimum of problem (TP) satisfying the constraint qualification (3.27) and the uniform positivity condition (4.17), then the shooting algorithm is locally quadratically convergent.

*Proof.* This is a consequence of the convergence result in [5, Theorem 5.4]. □

Note that in the proof of the above Theorem, it is established that the hypotheses imply that the derivative of the shooting function is injective.

## 5. APPLICATION TO A REGULATOR PROBLEM

Consider the following regulator problem, where $T = 5$:

$$ (5.1) \qquad \begin{array}{l} \min \frac{1}{2} \int_0^5 (x_{1,t}^2 + x_{2,t}^2) dt + \frac{1}{2} x_{1,5}^2, \\ \dot{x}_{1,t} = x_{2,t}, \quad \dot{x}_{2,t} = u_t \in [-1, 1], \end{array} $$
subject to the state constraint and initial conditions

$$ (5.2) \qquad x_{2,t} \ge -0.2, \quad x_{1,0} = 0, \quad x_{2,0} = 1. $$

To write the problem in the Mayer form, we introduce an auxiliary state variable given by the dynamics

$$ (5.3) \qquad \dot{x}_{3,t} = \frac{1}{2}(x_{1,t}^2 + x_{2,t}^2), \quad x_{3,0} = 0. $$

The resulting problem is then

$$ (5.4) \qquad \begin{array}{l} \min x_{3,5} + \frac{1}{2}x_{1,5}^2, \\ \dot{x}_1 = x_2, \quad \dot{x}_2 = u, \quad \dot{x}_3 = \frac{1}{2}(x_1^2 + x_2^2), \\ x_{1,0} = 0, \quad x_{2,0} = 1, \quad x_{3,0} = 0, \\ -1 \le u \le 1, \quad x_2 \ge -0.2. \end{array} $$

Using the optimal control solver BOCOP [40] we estimated that the optimal control $\hat{u}$ is a concatenation of a bang arc in the lower bound, followed by a constrained arc and ended with a singular one. Briefly, we can say that the optimal control has a B\_CS structure.

### 5.1. Checking local optimality.
In this subsection we compute analytically the optimal solution of (5.4) and check that it verifies the second order sufficient condition for state-constrained control-affine problems proved in [3, Theorem 5]. While the problem is convex, and hence, satisfying the first order optimality conditions is enough for proving optimality, the quoted second order conditions are of interest since they imply the quadratic growth, see [3, Definition 3].

To problem (5.4) we associate the functions $g: \mathbb{R}^3 \to \mathbb{R}$, $f_0, f_1: \mathbb{R}^3 \to \mathbb{R}^3$ given by

$$ g(x) := -x_2 - 0.2, \quad f_0(x) := \begin{pmatrix} x_2 \\ 0 \\ \frac{1}{2}(x_1^2 + x_2^2) \end{pmatrix}, \quad f_1(x) = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}. $$

So that the optimal control $\hat{u}$ is equal to -1 on the B\_ arc, and to $-\frac{g'(\hat{x})f_0(\hat{x})}{g'(\hat{x})f_1(\hat{x})} = 0$ on C according to formula (2.11), where $\hat{x}$ is the associated optimal state.

The pre-Hamiltonian of (5.4) reads

$$ H(u, x, p) := p_1x_2 + p_2u + \frac{1}{2}p_3(x_1^2 + x_2^2), $$

where $p := (p_1, p_2, p_3)$ is the costate.

Over the singular arc S, the inequality $u_{\min} < \hat{u}_t < u_{\max}$ and the minimum condition of Pontryagin's Maximum Principle (see e.g. [3, equation (2.12)]) imply that

$$ (5.5) \qquad H_u = 0. $$

Differentiating in time, (see e.g. [30, 5]), we obtain

$$ (5.6) \qquad H_u = pf_1, \quad \dot{H}_u = p[f_1, f_0], \quad \ddot{H}_u = p[[f_1, f_0], f_0] + \hat{u}p[[f_1, f_0], f_1], $$

where $[X, Y] := X'Y - Y'X$ denotes the Lie bracket associated with a pair of vector fields $X, Y: \mathbb{R}^n \to \mathbb{R}^n$. In this control-affine case, the control variable does not appear explicitly neither in the expression of $H_u$ nor in its time derivative $\dot{H}_u$.
So, if $p[[f_1, f_0], f_1]]$ takes only nonzero values along singular arcs, we obtain an expression for the control on singular arcs, namely

$$ (5.7) \qquad \hat{u} = - \frac{p[[f_1, f_0], f_0]](\hat{x})}{p[[f_1, f_0], f_1]](\hat{x})}. $$

The involved Lie brackets for this examples are

$$ (5.8) \qquad [f_1, f_0] = \begin{pmatrix} -1 \\ 0 \\ -x_2 \end{pmatrix}, \quad [[f_1, f_0], f_0] = \begin{pmatrix} 0 \\ 0 \\ x_1 \end{pmatrix}, \quad [[f_1, f_0], f_1] = \begin{pmatrix} 0 \\ 0 \\ -1 \end{pmatrix}. $$

On the other hand, the costate equation on the singular arc $S$ gives

$$ \begin{array}{l} \dot{p}_1 = -p_3 \hat{x}_1, \quad p_{1,5} = \hat{x}_{1,5}, \\ (5.9) \qquad dp_2 = -(p_1 + p_3 \hat{x}_2)dt + \nu, \quad p_{2,5} = 0, \\ \qquad \qquad \qquad \dot{p}_3 = 0, \quad p_{3,5} = 1, \end{array} $$

where $\mu$ is the multiplier associated with the state constraint, and $\nu$ is the density of $\mu$. Thus,

$$ (5.10) \qquad p_3 \equiv 1, $$

and from (5.7) and (5.8) we get

$$ (5.11) \qquad \hat{u} = \hat{x}_1 \quad \text{on } S. $$

Moreover, the first order optimality conditions imply that

$$ (5.12) \qquad 0 = H_u = p_2 \quad \text{on } C \cup S. $$

Let us write $\hat{\tau}_1, \hat{\tau}_2$ for the switching times, so that

$$ B_- = [0, \hat{\tau}_1], \quad C = [\hat{\tau}_1, \hat{\tau}_2], \quad S = [\hat{\tau}_2, 5]. $$

Since the control $\hat{u}$ is constantly equal to $-1$ on $B_-$, then

$$ (5.13) \qquad \hat{x}_{2,t} = 1 - t \quad \text{on } B_-, $$

until it saturates the state constraint at $\hat{\tau}_1$. Hence $1 - \hat{\tau}_1 = -0.2$, so it follows that

$$ (5.14) \qquad \hat{\tau}_1 = 1.2 $$

and $B_- = [0, 1.2]$. Consequently,

$$ (5.15) \qquad \hat{x}_{1,t} = t - \frac{t^2}{2} \quad \text{on } B_- $$

On $C = [1.2, \hat{\tau}_2]$, necessarily $\hat{x}_2 \equiv -0.2$, thus $u \equiv 0$ and

$$ (5.16) \qquad \hat{x}_{1,t} = 0.48 - 0.2(t - 1.2) = 0.72 - t/5 \quad \text{on } C. $$

On the singular arc $S = [\hat{\tau}_2, 5]$ we get, from the expression of $\hat{u}$ on $S$ (5.11), that $\ddot{\hat{x}}_1 = \dot{\hat{x}}_2 = \hat{u} = \hat{x}_1$. Thus

$$ (5.17) \qquad x_{1,t} = c_1 e^{5-t} + c_2 e^{t-5} \quad \text{on } S, $$

for some real constant values $c_1, c_2$. Therefore,

$$ (5.18) \qquad \hat{x}_{2,t} = \dot{\hat{x}}_{1,t} = -c_1 e^{5-t} + c_2 e^{t-5} \quad \text{on } S. $$

The stationarity condition $0 = H_u = p_2$ on $S$ yields $0 = \dot{p}_2 = -p_1 - \hat{x}_2$, where the second equality of latter equation follows from (5.9) and since $\nu \equiv 0$ on $S$. Thus $p_1 = -\hat{x}_2$ on $S$, so

$$ (5.19) \qquad p_{1,t} = c_1 e^{5-t} - c_2 e^{t-5} \quad \text{on } S. $$
The transversality condition for $p_1$ (see (5.9)) implies $c_1+c_2 = c_1-c_2$. Thus, $c_2 = 0$ and

$$ (5.20) \qquad \hat{x}_1 = -\hat{x}_2 = p_1 = c_1 e^{5-t} \quad \text{on } S. $$

Since $-0.2 = x_{2,\hat{\tau}_2}$, then $x_{1,\hat{\tau}_2} = 0.2$. Hence, from the expression of $\hat{x}_1$ on $C$ given in (5.16), we obtain $0.72 - \hat{\tau}_2/5 = 0.2$ and, consequently

$$ (5.21) \qquad \hat{\tau}_2 = 2.6. $$

Additionally, $-0.2 = x_{2,\hat{\tau}_2} = -c_1 e^{2.4}$ so that $c_1 = 0.2e^{-2.4}$. At time $\hat{\tau}_2$ we have that $p_1 = c_1 e^{2.4} = 0.2$ and, from the costate equation (5.9) and from (5.16), we have

$$ (5.22) \qquad \dot{p}_{1,t} = -\hat{x}_{1,t} = t/5 - 0.72 \le 2.6/5 - 0.72 < 0 \quad \text{on } C = [1.2, 2.6]. $$

So $p_1$ is decreasing on $C$ and, recalling (5.12), this implies that

$$ (5.23) \qquad \nu_t = p_{1,t} + \hat{x}_{2,t} > p_{1,\hat{\tau}_2} - 0.2 = 0 \quad \text{for } t \in [\hat{\tau}_1, \hat{\tau}_2). $$

Thus the complementarity condition for the state constraint [3, equation(3.4)(ii)] holds true.

In order to check that the strict complementarity hypothesis (i) of [3, Theorem 5] is satisfied, we will prove that the corresponding strict complementarity condition for the control constraint of [3, Definition 4] holds. In view of (5.22), we get

$$ (5.24) \qquad p_{1,t} = 0.2 + (t - \hat{\tau}_2)^2/10 - 0.72(t - \hat{\tau}_2). $$

Therefore $p_{1,\hat{\tau}_1} = 0.2 + (1.4)^2/10 + 0.72 \times 1.4 = 1.404$. On the $B_-$ arc, we have seen that $\hat{x}_1 > 0$ due to (5.15) and $x_{1,0} = 0$, so that $\dot{p}_1 = -\hat{x}_1 < 0$, and since $p_{1,\hat{\tau}_1} > 0$, it follows that $p_1$ has values greater than $p_{1,\hat{\tau}_1} = 1.404$. Therefore, since $\hat{x}_2 > -0.2$ over $B_-$,

$$ (5.25) \qquad \dot{p}_2 = -p_1 - \hat{x}_2 < 0 \quad \text{on } B_- $$

Since $p_{2,\hat{\tau}_1} = 0$, we get $p_2 > 0$ on $[0, \hat{\tau}_1)$ or, equivalently, $H_u > 0$ on $[0, \hat{\tau}_1)$. We conclude that the strict complementarity condition for the control constraint given in [3, Definition 4] holds. This completes the verification of [3, Theorem 5 - (i)].

Let us now verify the uniform positivity [3, equation (4.6)]. The dynamics for the linearized state is

$$ (5.26) \qquad \begin{aligned} \dot{z}_1 &= z_2, & \dot{z}_2 &= v, & \dot{z}_3 &= \hat{x}_1 z_1 + \hat{x}_2 z_2, \\ z_{1,0} &= z_{2,0} = z_{3,0} = 0. \end{aligned} $$

Let $C_S$ and $P_*^2$ denote the strict critical cone and the extended cone (invoking the notation used in [3]) at the optimal trajectory $(\hat{u}, \hat{x})$, respectively. Since $\hat{u}$ is $B_-CS$ then, for any $(v, z) \in C_S$, $v = 0$ on the initial interval $B_-$. Consequently,

$$ (5.27) \qquad y = 0 \quad \text{on } B_-, \text{ for all } (y, \xi, h) \in P_*^2. $$

On the other hand, the dynamics for the transformed linearized state is

$$ (5.28) \qquad \dot{\xi}_1 = \xi_2 + y, \quad \dot{\xi}_2 = 0, \quad \dot{\xi}_3 = \hat{x}_1\xi_1 + \hat{x}_2\xi_2 + \hat{x}_2y, $$

$$ (5.29) \qquad \xi_{1,0} = \xi_{2,0} = \xi_{3,0} = 0. $$

Take $(y, \xi, h) \in P_*^2$. Then

$$ (5.30) \qquad \xi_2 = 0 \quad \text{on } [0, 5]. $$

In view of [3, equation (3.15)(i)], we have that $-\xi_2 - y = 0$ on $C$. Thus, due to (5.30), we get that

$$ (5.31) \qquad y = 0 \quad \text{on } B_- \cup C. $$
Thus, from (5.28)-(5.31), we get $\xi_1 = \xi_3 = 0$ on $B_- \cup C$. Regarding the last component $h$ of the considered critical direction $(y, \xi, h) \in \mathcal{P}_*^2$, in view of the linearized cost equation [3, equation (3.16)] and due to the fact that there are no final constraints, we get that

$$ (5.32) \qquad \hat{x}_{1,5} \xi_{1,5} + \xi_{3,5} = 0. $$

Then, we deduce that there is no restriction on $h$. We obtain

$$ (5.33) \quad \mathcal{P}_*^2 = \left\{ \begin{array}{l} (y, \xi, h) \in L^2 \times (H^1)^n \times \mathbb{R} : y = \xi_1 = \xi_2 = 0 \text{ on } B_- \cup C, \\ \xi_2 = 0 \text{ on } [0, 5], \quad \dot{\xi}_1 = y \text{ and } \dot{\xi}_3 = \hat{x}_1 \xi_1 + \hat{x}_2 y \text{ on } S, \quad (5.32) \text{ holds} \end{array} \right\}. $$

The quadratic forms $Q$ and $\Omega$ are given by

$$ (5.34) \quad Q(v, z) := \int_0^T (z_1^2 + z_2^2)dt + z_{1,5}^2; \quad \Omega(y, h, \xi) := \int_0^T (\xi_1^2 + (\xi_2 + y)^2)dt + h^2. $$

Thus, $\Omega$ is a Legendre form on $\{(y, h, \xi) \in L^2 \times \mathbb{R} \times (H^1)^3 : (5.28) \text{ holds}\}$ and is coercive on $\mathcal{P}_*^2$. Hence Theorem 5 in [3] holds. Consequently, $(\hat{u}, \hat{x})$ is a Pontryagin minimum of problem (5.4) satisfying the $\gamma$-growth condition.

## 5.2. Transformed problem.
Here we transform the problem (5.4) to obtain a problem with neither control nor state constraints, as done in Section 3.1.

The optimal control associated with (5.4) has a $B_{-CS}$ structure, as said above. Then we triplicate the number of state variables obtaining the new variables $X_1, \dots, X_9$ and we consider two switching times that we write $X_{10}, X_{11}$. The new problem has only one control variable that corresponds to the singular arc and which we call $V$. The reformulation of (5.4) is as follows

$$ \begin{array}{l} \min X_{9,1} + \frac{1}{2}X_{7,1}^2, \\ \dot{X}_1 = X_{10}X_2, \\ \dot{X}_2 = -X_{10}, \\ \dot{X}_3 = \frac{1}{2}X_{10}(X_1^2 + X_2^2), \\ \dot{X}_4 = (X_{11} - X_{10})X_5, \\ \dot{X}_5 = 0, \\ \dot{X}_6 = \frac{1}{2}(X_{11} - X_{10})(X_4^2 + X_5^2), \\ \dot{X}_7 = (5 - X_{11})X_8, \\ \dot{X}_8 = (5 - X_{11})V, \\ \dot{X}_9 = \frac{1}{2}(5 - X_{11})(X_7^2 + X_8^2), \\ \dot{X}_{10} = 0, \\ \dot{X}_{11} = 0, \\ X_{1,0} = 0, \quad X_{2,0} = 1, \quad X_{3,0} = 0, \\ X_{1,1} = X_{4,0}, \quad X_{2,1} = X_{5,0}, \quad X_{3,1} = X_{6,0}, \\ X_{4,1} = X_{7,0}, \quad X_{5,1} = X_{8,0}, \quad X_{6,1} = X_{9,0}, \\ X_{2,1} = 0.2, \quad (\text{or } g(x_2(\tau_2)) = 0). \end{array} $$

From (5.11) we deduce that $V = X_7$.
We solved problem (5.35) by the shooting algorithm proposed in Section 3. The graphics of the optimal control and states are shown in Figure 1 in the original variables *u* and *x*, and the corresponding costate variables $p_1, p_2, p_3$ are displayed in Figure 2.

In our numerical tests we take 1000 time steps. The optimal switching times obtained numerically are $\hat{\tau}_1 = 1.2$ and $\hat{\tau}_2 = 2.6036023$, so these values agree with the ones founded analytically in (5.14) and (5.21), while the optimal cost is 0.3934884 and the shooting function evaluated at the optimal trajectory is $5. \times 10^{-6}$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Optimal control and state variables</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Costate variables</figcaption></figure>

*Remark 5.1.* In order to compare with existing methods, we have solved the same problem using GEKKO Python [7]. We have tested different variable definitions for the control input and we found the best result by setting the control as a *manipulated variables* [7]. The results are shown in Figure 3. On the right of this figure, we exhibit a zoom around the second switching point $\hat{\tau}_2$ of the optimal control. We can see that the method shows a discontinuity of the control at $\hat{\tau}_2$, but that the approximation of $\hat{\tau}_2$ is between 2.57 and 2.58, while we have calculated analytically that $\hat{\tau}_2 = 2.6$ and our shooting algorithm finds the approximate value 2.6036023. GEKKO was tested with up to 1000 time steps, while our shooting algorithm was run with 150 time steps. This indicates that our method is more accurate when it comes to finding switching times and approximating bang-singular solutions. This
feature of shooting methods has been already observed in the literature (see e.g. [41].

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Numerical solution using GEKKO. On the left: optimal control and state variables. On the right: zoom of the optimal control around the second switching time.</figcaption></figure>

## 6. CONCLUSION

We have shown that, for problems that are affine w.r.t. the control, and for which the time interval can be partitioned in a finite set of arcs, the shooting algorithm converges locally quadratically if a second order sufficient condition holds, and thus the method may be an efficient way to get a highly accurate solution.

The essential tool was to formulate a transformed problem. We think that this approach could be extended to more general problems with several state constraints and control variables, the latter possibly entering nonlinearly in the problem. Problems with vector control are important in complex models, see e.g. [28, 27] where they studied a control-affine problem with vector control, state constraints, but no singular arcs, generically, in the state-constrained solution. When considering several controls and state constraints, it becomes important to take into account the possibility of having high-order state constraints.

## ACKNOWLEDGEMENTS

The first author was supported by FAPERJ (Brazil) through the *Jovem Cientista do Nosso Estado* Program and by CNPq (Brazil) through the *Universal* Program and the Productivity in Research Scholarship. The second author was supported by the FiME Lab Research Initiative (Institut Europlace de Finance), and by the PGMO program.

We acknowledge the anonymous reviewers for their careful reading and comments that helped us improve this manuscript.

The first author thanks Gabriel de Lima Monteiro for his suggestions on the numerical implementations.

## STATEMENTS AND DECLARATIONS

The authors have no conflicts of interest to declare that are relevant to the content of this article.
REFERENCES

[1] R.A. Adams. Sobolev spaces. Academic Press, New York, 1975.

[2] A. Aftalion and J. Bonnans. Optimization of running strategies based on anaerobic energy and variations of velocity. *SIAM J. Applied Math.*, 74(5):1615–1636, 2014.

[3] M.S. Aronna, J.F. Bonnans, and Goh B.S. Second order analysis of control-affine problems with scalar state constraint. *Math. Program.*, 160(1):115–147, 2016.

[4] M.S. Aronna, J.F. Bonnans, A.V. Dmitruk, and P.A. Lotito. Quadratic order conditions for bang-singular extremals. *Numer. Algebra, Control Optim., AIMS Journal, special issue dedicated to Professor Helmut Maurer on the occasion of his 65th birthday*, 2(3):511–546, 2012.

[5] M.S. Aronna, J.F. Bonnans, and P. Martinon. A shooting algorithm for optimal control problems with singular arcs. *J. Optim. Theory Appl.*, 158(2):419–459, 2013.

[6] Aram Arutyunov and Dmitry Karamzin. A survey on regularity conditions for state-constrained optimal control problems and the non-degenerate maximum principle. *Journal of Optimization Theory and Applications*, 184:697–723, 2020.

[7] Logan Beal, Daniel Hill, R Martin, and John Hedengren. Gekko optimization suite. *Processes*, 6(8):106, 2018.

[8] H.G. Bock and K.J. Plitt. A multiple shooting algorithm for direct solution of optimal control problems*. *IFAC Proceedings Volumes*, 17(2):1603–1608, 1984. 9th IFAC World Congress: A Bridge Between Control Science and Technology, Budapest, Hungary, 2-6 July 1984.

[9] J.F. Bonnans. Convex and Stochastic Optimization. Springer-Verlag, Berlin, 2019.

[10] J.F. Bonnans, J.C. Gilbert, C. Lemaréchal, and C.A. Sagastizábal. *Numerical optimization: theoretical and practical aspects*. Springer Science & Business Media, Berlin, 2006.

[11] J.F. Bonnans and A. Shapiro. *Perturbation analysis of optimization problems*. Springer-Verlag, New York, 2000.

[12] B. Bonnard, L. Faubourg, G. Launay, and E. Trélat. Optimal control with state constraints and the space shuttle re-entry problem. *Journal of Dynamical and Control Systems*, 9(2):155–199, 2003.

[13] R. Bulirsch. Die Mehrzielmethode zur numerischen Lösung von nichtlinearen Randwertproblemen und Aufgaben der optimalen Steuerung. *Report der Carl-Cranz Gesellschaft*, 1971.

[14] O. Cots. Geometric and numerical methods for a state constrained minimum time control problem of an electric vehicle. *ESAIM: Control, Optimisation and Calculus of Variations*, 23(4):1715–1749, 2017.

[15] Olivier Cots, Joseph Gergaud, Damien Goubinat, and Boris Wembe. Singular versus boundary arcs for aircraft trajectory optimization in climbing phase. *ESAIM: M2AN (Online first)*, 2022.

[16] M.R. de Pinho, M.M. Ferreira, U. Ledzewicz, and H. Schaettler. A model for cancer chemotherapy with state-space constraints. *Nonlinear Analysis: Theory, Methods & Applications*, 63(5):e2591–e2602, 2005.

[17] A.V. Dmitruk. Quadratic conditions for a weak minimum for singular regimes in optimal control problems. *Soviet Math. Doklady*, 18(2), 1977.

[18] A.V. Dmitruk. Quadratic order conditions for a Pontryagin minimum in an optimal control problem linear in the control. *Math. USSR Izvestiya*, 28:275–303, 1987.

[19] R. Fletcher. *Practical methods of optimization*. John Wiley & Sons, New Jersey, 2013.

[20] H. Frankowska and D. Tonon. Pointwise second-order necessary optimality conditions for the Mayer problem with control constraints. *SIAM J. Control Optim.*, 51(5):3814–3843, 2013.

[21] G. Fraser-Andrews. Numerical methods for singular optimal control. *J. Optim. Theory Appl.*, 61(3):377–401, 1989.

[22] A. T. Fuller. Study of an optimum non-linear control system. *J. Electronics Control (1)*, 15:63–71, 1963.

[23] B.S. Goh. Necessary conditions for singular extremals involving multiple control variables. *J. SIAM Control*, 4:716–731, 1966.

[24] B.S. Goh. The second variation for the singular Bolza problem. *J. SIAM Control*, 4(2):309–325, 1966.

[25] B.S. Goh, G. Leitmann, and T.L. Vincent. Optimal control of a prey-predator system. *Math. Biosci.*, 19:263–286, 1974.
[26] K. Graichen and N. Petit. Solving the Goddard problem with thrust and dynamic pres-
sure constraints using saturation functions. In *17th World Congress of The International
Federation of Automatic Control*, volume Proc. of the 2008 IFAC World Congress, pages
14301–14306, Seoul, 2008. IFAC.

[27] Clara Leparoux, Bruno Hérissé, and Frédéric Jean. Optimal planetary landing with pointing and glide-slope constraints. In *2022 IEEE 61st Conference on Decision and Control (CDC)*, pages 4357–4362. IEEE, 2022.

[28] Leparoux, Clara, Hérissé, Bruno, and Jean, Frédéric. Structure of optimal control for planetary landing with control and state constraints. *ESAIM: COCV*, 28:67, 2022.

[29] E. S. Levitin, A. A. Milyutin, and N. P. Osmolovskiǐ. Higher order conditions for local minima in problems with constraints. *Uspekhi Mat. Nauk*, 33(6(204)):85–148, 272, 1978.

[30] H. Maurer. Numerical solution of singular control problems using multiple shooting tech-
niques. *J. Optim. Theory Appl.*, 18(2):235–257, 1976.

[31] H. Maurer and M.R. De Pinho. Optimal control of epidemiological seir models with 11-objectives and control-state constraints. *Pacific Journal of Optimization*, 12(2):415–436, 2016.

[32] H. Maurer and W. Gillessen. Application of multiple shooting to the numerical solution of optimal control problems with bounded state variables. *Computing*, 15(2):105–126, 1975.

[33] H. Maurer, J.-H.R. Kim, and G. Vossen. On a state-constrained control problem in optimal production and maintenance. In C. Deissenberg and R.F. Hartl, editors, *Optimal Control and Dynamic Games*, volume 7 of *Advances in Computational Management Science*, pages 289–308. Springer, US, 2005.

[34] D.D. Morrison, J.D. Riley, and J.F. Zancanaro. Multiple shooting method for two-point boundary value problems. *Comm. ACM*, 5:613–614, 1962.

[35] H.J. Oberle. Numerical computation of singular control problems with application to optimal heating and cooling by solar energy. *Appl. Math. Optim.*, 5(4):297–314, 1979.

[36] H.J. Pesch. A practical guide to the solution of real-life optimal control problems. *Control Cybernet.*, 23(1-2):7–60, 1994. Parametric optimization.

[37] H.M. Robbins. A generalized Legendre-Clebsch condition for the singular case of optimal control. *IBM J. of Research and Development*, 11:361–372, 1967.

[38] H. Schättler. Local fields of extremals for optimal control problems with state constraints of relative degree 1. *J. Dyn. Control Syst.*, 12(4):563–599, 2006.

[39] H. Seywald and E.M. Cliff. Goddard problem in presence of a dynamic pressure limit. *J. Guidance Control Dynam.*, 16(4):776–781, 1993.

[40] Inria Saclay Team Commands. Bocop: an open source toolbox for optimal control, 2017.

[41] Emmanuel Trélat. Optimal control and applications to aerospace: some results and challenges. *Journal of Optimization Theory and Applications*, 154:713–758, 2012.

[42] R. Vinter. *Optimal control. Systems & Control: Foundations & Applications*. Birkhäuser Boston, Inc., Boston, MA, 2000.

ESCOLA DE MATEMÁTICA APLICADA FGV EMAP, 22250-900 RIO DE JANEIRO,
BRAZIL

E-mail address: SOLEDAD.ARONNA@FGV.BR

INRIA SACLAY, L2S, CENTRALESUPELEC, 91190 GIF-SUR-YVETTE, FRANCE

E-mail address: FREDERIC.BONNANS@INRIA.FR

SCHOOL OF ELECTRICAL ENGINEERING, COMPUTING AND MATHEMATICAL SCIENCES CURTIN UNIVERSITY, PERTH, AUSTRALIA

E-mail address: GOH2OPTIMUM@GMAIL.COM