Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr

Jérémy Rouot
email: jeremy.rouot@grenoble-inp.org

Towards Geometric Time Minimal Control without Legendre Condition and with Multiple Singular Extremals for Chemical Networks

Keywords: Mass action chemical systems, Pontryagin Maximum Principle, Geometric optimal control, Time minimal synthesis, McKeithan type chemical scheme Mathematics Subject Classification: 49K15, 58K45, 92B05

This article deals with the problem of maximizing the production of a species for a chemical network by controlling the temperature. Under the socalled mass kinetics assumption the system can be modeled as a single-input control system using the Feinberg-Horn-Jackson graph associated to the reactions network. Thanks to Pontryagin's Maximum Principle, the candidates as minimizers can be found among extremal curves, solutions of a (non smooth) Hamiltonian dynamics and the problem can be stated as a time minimal control problem with a terminal target of codimension one. Using geometric control and singularity theory the time minimal syntheses (closed loop optimal control) can be classified near the terminal manifold under generic conditions. In this article, we focus to the case where the generalized Legendre-Clebsch condition is not satisfied, which paves the road to complicated syntheses with several singular arcs. In particular, it is related to the situation for a weakly reversible network like the McKeithan scheme of two reactions: T+M A B .

Introduction

The optimization of the production is an important problem in chemical and biological engineering and the control can be either the temperature (batch or closed reactor) or by feeding the reactor (semi-batch or open case). In this article we shall concentrate on the first case. Moreover we assume that the dynamics is modeled using the mass action kinetics assumptions and hence given at constant temperature T by a polynomial system based only on the Feinberg-Horn-Jackson graph associated to the chemical network. Also in the 70's those researchers obtained (under the so-called zero deficiency assumption) in a series of seminal articles [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF][START_REF] Horn | General mass action kinetics[END_REF] a complete description of the dynamics, at constant temperature. If this dynamics is well understood, in the optimal problem the temperature is not constant and the analysis becomes very intricate. Thanks to the Pontryagin Maximum Principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF] candidates as minimizers can be found among extremals solutions of a (non smooth) Hamiltonian dynamics and optimal solutions are concatenation of bang arcs, with minimum and maximum temperature, and the so-called singular arcs, defined as a solution of a smooth Hamiltonian constrained dynamics [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. Moreover maximizing the production of one species during the batch can be restated as producing a fixed amount of this species while minimizing the batch duration. In this frame, the problem is a time minimal control problem, with a terminal target manifold of codimension one.

For applications, the time optimal control has to be computed as a closed loop feedback and this leads to the problem of computing the time minimal synthesis, for a single-input control system. At the end of the 80's, geometric optimal control has produced an important literature to compute the time minimal syntheses for the fixed end point case, where the initial point is localized near the terminal point. This was done using the Lie algebraic structure of the control system, mainly for single-input (smooth) system, where the control appears linearly, some seminal references are [START_REF] Sussmann | Synthesis, presynthesis, sufficient conditions for optimality and subanalytic sets[END_REF][START_REF] Schättler | On the local structure of time-optimal bang-bang trajectories in R 3[END_REF][START_REF] Boscain | Optimal Syntheses for Control Systems on 2-D Manifolds[END_REF] either in general context or in view of application formed by a sequence of two irreversible reactions: A → B → C and in relation with an industrial project [START_REF] Bonnard | Time minimal control of batch reactors[END_REF].

Our aim is to extend this work to more complicated reaction schemes and to deal in particular with weakly reversible chemical schemes. More precisely we shall concentrate on a McKeithan type scheme of the form T+M A B assuming that the coefficient governing the dynamics are given by Arrhenius law. This scheme was already studied in the context of control theory using stabilization techniques with "feeding" types control [START_REF] Sontag | Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction[END_REF] and ad hoc observer design [START_REF] Chavez | Observer design for a class of nonlinear systems, with applications to chemical and biological networks[END_REF]. In our case, this network is a test bed case for our very general approach.

The key point for this extension is the analysis of singular trajectories and their role in the synthesis. This is connected with an important question and the need to extend the standard synthesis related to the turnpike phenomenon [START_REF] Sussmann | The Structure of Time-Optimal Trajectories for Single-Input Systems in the Plane: the C ∞ Nonsingular Case[END_REF] to deal with cases, where the strict Legendre-Clebsch condition is not satisfied, a situation encountered in a recent application in MRI [START_REF] Bonnard | Time minimal saturation of a pair of spins and application in magnetic resonance imaging[END_REF].

The organization of this article is the following. In Section 2, we recall briefly the Feinberg-Horn-Jackson theory to model the dynamics of chemical networks of constant temperature and the properties of the dynamics under the zero deficiency assumption [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF][START_REF] Horn | General mass action kinetics[END_REF], which can be applied to the McKeithan scheme. The stability properties are recalled and can be applied to control stabilization and observer design [START_REF] Sontag | Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction[END_REF][START_REF] Chavez | Observer design for a class of nonlinear systems, with applications to chemical and biological networks[END_REF]. In Section 3, we present the fundamental results of the time minimal control problem, which are relevant to our study: Pontryagin Maximum Principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF], regular and singular extremals [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF], [START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF]. The general turnpike theorem [START_REF] Sussmann | The Structure of Time-Optimal Trajectories for Single-Input Systems in the Plane: the C ∞ Nonsingular Case[END_REF] is recalled and extensions are presented in relation with the problem with terminal manifold of codimension one and when the strict Legendre-Clebsch condition is not satisfied. The concept of conjugate and focal points is introduced based on [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF]. In Section 4, we analyze the time minimal control problem for a two reactions McKeithan scheme. To compute the time minimal syntheses, we present techniques and results from [START_REF] Bonnard | Time minimal control of batch reactors[END_REF][START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF], which have to be extended to analyze the problem. The computational complexity of the problem is discussed and symbolic computations are presented to cope with this complexity.

Mathematical model and stability properties of the McKeithan network

2.1 Mass action kinetics networks and dynamics using the Feinberg-Horn-Jackson graph ( [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF], [START_REF] Horn | General mass action kinetics[END_REF])

We consider a set of m chemical species {X1, . . . , Xm} and the state of the dynamics is the vector c = (c1, . . . , cm) ∈ R m ≥0 representing the molar concentration. Let R be a set of reactions, each reaction being denoted by y → y and of the form

m i=1 αiXi -→ m i=1 βiXi,
where αi, βi are the stoichiometric coefficients and the vectors y = (α1, . . . , αm) and y = (β1, . . . , βm) are the vertices of the so-called Feinberg-Horn-Jackson oriented graph associated to the network, edges being oriented according to y → y . Each reaction is characterized by a reaction rate K(y → y ) and the system is said simple (or mass kinetics) if the rate of the reaction is of the form:

K(y → y ) = k(T ) c y , c y = c α 1 1 . . . c αm m (1.1) and k(T ) = A e -E/(RT )
is the Arrhenius law, A is the exponential factor, E is the activation energy, both depending on the reaction, R is the gas constant and T is the temperature. Note that different rate formulae can be used to deal in particular with biomedical systems (see for instance [START_REF] Sontag | Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction[END_REF]). The dynamics of the system, taking into account the whole network is:

ċ(t) = f (c(t), T ) = y→y K(y → y ) (y -y). (1.2)

More explicit representation of the dynamics

Definition 1. The stoichiometric subspace is S := span{y -y ; y → y ∈ R} and the sets (c(0) + S) ∩ R m ≥0 are called the (strictly if > 0) positive stoichiometric compatibility classes.

From [START_REF] Anderson | Global asymptotic stability for a class of nonlinear chemical equations[END_REF] we have. Lemma 1. Let c(t) be a solution of (1.2) with initial condition c(0) ∈ R m ≥0 . Then c(t) belongs for all t ≥ 0 to the strictly positive compatibility class (c(0) + S) ∩ R m ≥0 .

Definition 2. Having labeled the set of vertices by i = 1, . . . , n, with corresponding stoichiometric vector (y1, . . . , yn), the complex matrix is Y := (y1, . . . , yn). The incidence connectivity matrix A := (aij) contains the Arrhenius coefficients ki of the reactions using the rule: k1 = a21 indicates a reaction with kinetics constant k1 from the first node to the second, that is y1 →

k 1 y2.
With the mass kinetics assumption, the dynamics can be expressed as

ċ(t) = f (c(t), T ) = Y Ãc Y ,
where à is the Laplacian matrix in graph theory defined by

à = A -diag n i=1 ai1, . . . , n i=1 ain (1.3)
and we denote c Y = (c y 1 , . . . , c yn ) .

The McKeithan scheme ([23, 27])

It is given by the reaction scheme:

T + M C0 C1 . . . CN k1 kp,0 kp,1 kp,N k-1,0 k-1,1 k-1,N .
The matrix Y is given by

Y =          1 0 . . . . . . 0 1 0 • 0 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 0 0 . . . . . . 0 1          and the matrix A = (aij) is defined by a21 = k1, a1i = k-1,i-2, i = 2, . . . , m (m = N + 2), ai,i-1 = kp,i-3, i = 3, .
. . , m and all others aij are zero. The stoichiometric subspace is defined by:

{c : T + C0 + . . . + CN = M + C0 + . . . + CN = 0}
and we note: δ1 = T + C0 + . . . + CN and δ2 = M + C0 + . . . + CN the constants associated to first integrals of the dynamics.

Consider the special case N = 2 so that the reaction scheme is denoted by

T+M A B k 1 k 2 k 3 k 4
, and restricting to the stoichiometric class (δ1 and δ2 being fixed), one gets with: x := [A], y := [B], [•] denoting the respective concentrations, that the dynamics is described by the equations

ẋ = k1(δ1x -y)(δ2 -x -y) -(k2 + k3)x ẏ = k2x -k4y. (1.4)
Definition 3. The deficiency of the network is: δ = n -l -s, where n is the number of vertices, l is the number of connected components and s is the dimension of the stoichiometric subspace. The network is called strongly connected if for each pair (i, j) of vertices such that there exists an oriented path joining i to j there exists a path joining j to i.

Using [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF], refined by [START_REF] Anderson | Global asymptotic stability for a class of nonlinear chemical equations[END_REF][START_REF] Sontag | Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction[END_REF], one has the following result.

Theorem 1. The graph associated to the McKeithan scheme is strongly connected and with deficiency zero. In each strictly positive compatibility class there exists in this domain an unique globally asymptotically stable equilibrium.

Application to stabilization and observer design for the McKeithan scheme

This stability result has consequences to control and observation properties of the network, see [START_REF] Sontag | Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction[END_REF][START_REF] Chavez | State-estimators for chemical reaction networks of Feinberg-Horn-Jackson zero deficiency type[END_REF], that we recall briefly. The dynamics (1.4) can be converted into a control system of the form

ċ(t) = f (x(t)) + u(t) y(t) = h(c(t)), (1.5)
where u(•) is a feeding control and h is a polynomial observation function. Asymptotic stability of equilibrium in each strictly positive compatibility class will allow to get stabilization result for a single equilibrium. Moreover it leads to design under a mild assumption (detectability) a simple observer. We refer to [START_REF] Sontag | Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction[END_REF] and [START_REF] Chavez | Observer design for a class of nonlinear systems, with applications to chemical and biological networks[END_REF] for the detailed presentation of those results, the geometric construction being clear.

3 The optimal control problem and Pontryagin Maximum Principle

Statement and notation for the optimal control problem

The system is written as dc dt = f (c, T ) (see (1.4)) and controlling the temperature leads to T ∈ [Tm, TM ]. In the sequel, we shall use the terminology direct for the corresponding control problem. In practice, thermodynamics has to be used to model the heat exchanges, in relation with the heat produced by the reactions [START_REF] Tóth | Reaction kinetics: exercises, programs and theorems[END_REF] or the heat exchange device used in the experiments, depending upon the technical achievements. To avoid this part of the study and without losing any mathematically generality, we shall use v as the control variable setting v = u, where v := ki(T ) for some reaction i, where ki(T ) = Aie -E i /(RT ) (see (1.1)). This leads to deal with the so-called indirect control system: q(t) = F (q(t)) + u(t) G(q(t)),

where q = (c, v) ∈ R n is the extended state variable, F = f (c, v) ∂ ∂q , G = ∂ ∂v , u-≤ u ≤ u+
, where [u-, u+] can be normalized to [-1, +1]. Note that the bounds v ∈ [vm, vM ] will not be taken into account in our study. The map v → v is the standard Goh transformation in optimal control, see [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF].

The optimal control problem of physical interest is the problem of maximizing the production of one species and using a proper variable labeling, the optimal problem is therefore of the Mayer type:

q = F (q) + u G(q), max |u|≤1 q1(t f ),
where t f is the time duration of the batch and q1 is the desired product.

Note it will be show (thanks to the Maximum Principle) that a "dual" formulation is min

u(•) t f , c1(t f ) = d,
where d > 0 is the desired amount of the species X1 during the batch.

Maximum Principle [24]

Notations and concepts

Consider a general control system of the form q = X(q, u), q ∈ R n , where X is a real analytic (C ω ) and the control is u : [0, t f (u)] → [-1, 1]. The set of admissible controls U is the set of bounded measurable mappings. If q(0) = q0 (initial state), we denote by q(•, q0, u) (in short q(•)) the solution starting from q0. Fixing t f , the accessibility set in time t f is the set A(q0, t f ) = ∪ u(•)∈U q(t f , q0, u). The extremity mapping (in time t f ) is the map: E q 0 ,t f : u(•) → q(t f , q0, u) defined on a domain of U; the set U is endowed with the L ∞ -norm topology.

Pontryagin Maximum Principle (PMP)

Statement in the time minimal case with q(t f ) ∈ N := smooth terminal manifold.

Notation: H(q, p, u) = p • X(q, u) denotes the pseudo-Hamiltonian (Hamiltonian lift of the vector field X), p is the adjoint vector in R n \ {0} and • is the scalar product. We denote by M (q, p) = max |u|≤1 H(q, p, u).

Statement of the PMP: If (q * , u * ) is an optimal control-trajectory pair on [0, t * f ] then there exists p * ∈ R n \ {0} such that a.e.: q * (t) = ∂H ∂p (q * (t), p * (t), u * (t)), ṗ * (t) = -∂H ∂q (q * (t), p * (t), u * (t)), H(q * (t), p * (t), u * (t)) = M (q * (t), p * (t)).

(1.6) Moreover M (q * (t), p * (t)) is a positive constant and p * satisfies the transversality condition p * (t f ) ⊥ T q * (t f ) N.

(1.7)

Statement in the Mayer case: As before, except that the transversality condition is replaced by

p * (t f ) = ∂ϕ ∂q (x * (t f )),
where ϕ is the Mayer cost function to minimize.

Definition 4. An extremal (q, p, u) is a solution of (1.6) on [0, t f ]. It is called a BC- extremal is q(0) = q0

Computations of singular extremals

First case. Consider the case q = X(q, u), where H = p • X(q, u) and the condition ∂H ∂u = 0 is satisfied. Denote by z(•) = (q(•), p(•)) the reference extremal. From the maximization condition, the Legendre-Clebsch condition ∂ 2 H ∂u 2 ≤ 0 has to be fulfilled. If this inequality is strict, one can use the implicit function theorem to compute the singular control as the dynamic feedback: z → us(z) and plugging such us into H(q, p, u) leads to define the true (or maximized) Hamiltonian. Second case. Let q = F (q) + u G(q). One introduces the following notations. If X and Y are two real analytic vector fields, the Lie bracket is defined by [X, Y ](q) = ∂X ∂q (q)Y (q) -∂Y ∂q (q)X(q). The extremal lift of X is HX (z) = p • X(q), z = (q, p) and the Poisson bracket is defined by {HX , HY }

(z) = p • [X, Y ](q).

Computation of singular extremals:

The condition

∂ ∂u d dt 2 ∂H ∂u = {{HG, HF }, HG} ≥ 0 (resp. > 0) is called the (resp. strict) generalized Legendre-Clebsch condition.
Recall the following.

Proposition 1 ([19]

). The generalized Legendre-Clebsch condition is a necessary optimality condition for the time minimal control problem with fixed extremities.

To compute the singular extremal, we differentiate twice t → HG(z(t)) and we get

HG(z) = {HG, HF }(z) = 0, {{HG, HF }, HF }(z) + u {{HG, HF }, HG}(z) = 0. (1.8)
Assume (in relation with the Legendre-Clebsch condition) {{HG, HF }, HG} = 0. The corresponding extremal is called of order 2 and the singular control u is computed as us(z), using relation (1.8). Plugging such us(z) into H(q, p, u) leads to define the true singular Hamiltonian, denoted by Hs(z). One has: Proposition 2. Singular extremals of order 2 are the solutions of :

dq dt = ∂Hs ∂p (z), dp dt = - ∂Hs ∂q (z) (1.9)
with the constraints {HG, HF }(z) = HG(z) = 0.

Moreover, in order to be admissible, the singular control is given by

us(z) = - {{HG, HF }, HF }(z) {{HG, HF }, HG}(z) (1.10)
and has to satisfy the admissibility constraint |us(z)| ≤ 1.

Definition 5. Let z = (q, p) be a singular extremal of order 2 and M = HF = h the constant value of the Hamiltonian. The extremal is called exceptional if h = 0. If h > 0, the extremal is called hyperbolic (resp. elliptic) if {{HG, HF }, HG} > 0 (resp. < 0).

The case of chemical networks

Recall that for our network ċ

= f (c, v), v = ki and it is extended into q = F (q) + u G(q) with F = f (c, v) ∂ ∂c and G = ∂ ∂v . Denote H = pc • f (c, v) and H = p • (F + u G)
the respective Hamiltonian lifts, with p = (pc, pv). One has the following relation between the corresponding singular extremals.

Lemma 2. The pair (q, p) is solution of

q = ∂H ∂p , ṗ = - ∂H ∂q , ∂H ∂u = 0
if and only if pv = 0 and (c, pc, pv) is solution of

ċ = ∂ H ∂pc , ṗc = ∂ H ∂c , ∂ H ∂v = 0,
and moreover the following relations are satisfied

d dt ∂H ∂u = {HG, HF } = - ∂ H ∂v (1.11) ∂ ∂u d 2 dt 2 ∂H ∂u = {{HG, HF }, HG} = - ∂ 2 H ∂v 2 .
(1.12)

In particular this gives the correspondence between both singular Hamiltonian and the respective Legendre-Clebsch and generalized Legendre-Clebsch conditions.

The case n = 3

We have q = (x, y, v) and introduce the following determinants:

D = det(G, [G, F ], [[G, F ], G]), D = det(G, [G, F ], [[G, F ], F ]), D = det(G, [G, F ], F ). Using p • G = p • [G, F ] = p • ([[G, F ] , F ] + u [[G, F ] , G]) = 0
and eliminating p leads to the following. Proposition 3. If n = 3, the singular control is given by the feedback us(q) = -D (q)/D(q) and singular trajectories are defined by the vector field Xs(q) = F (q) -D (q) D(q) G(q). Singular trajectories are hyperbolic if DD > 0, elliptic if DD < 0 and exceptional if D = 0.

Optimality status: the case n = 3. We use [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF] to describe the optimality status of singular trajectories. The system is q = F + u G and we relax the bound |u| ≤ 1, assuming u ∈ R so that singular arcs are admissible. We assume the following:

(H0) F, G are linearly independent, (H1) G, [G, F ] are linearly independent.
One picks a smooth singular arc z(•) = (q(•), p(•)) defined on [0, t f ] so that p is unique up to a non zero multiplicative scalar and q(•) is a one-to-one immersion. We have: Proposition 4. There exists a C 0 -neighborhood U of the reference singular arc t → q(t), t ∈ [0, t f ] so that q(•) is time minimal (resp. maximal) if q(•) is hyperbolic (resp. elliptic) up to the first conjugate time t1c with respect to all trajectories with the same extremities contained in U. In the exceptional case, the reference singular arc is time minimal (and time maximal).

Algorithm to compute the first conjugate point in the hyperbolic-elliptic case. Let V (t), t ∈ [0, t f ] be the solution of the (variational) Jacobi equation:

δq(t) = ∂X ∂q (z(t)) δq(t)
with initial condition V (0) = G(q(0)), where t → q(t) is the reference singular arc. The first conjugate point t1c is the first t > 0 such that V (t) is collinear to G(q(t)).

See [5, p. 123] for a proof and the geometric interpretation.

Small time classification of regular extremals

In this section, we recall the seminal results coming from singularity theory due to [START_REF] Ekeland | Discontinuité des champs hamiltoniens et existence des solutions optimales en calcul des variations[END_REF][START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF] to analyze the small time extremal curves near the switching surface.

Definition 6. We denote by σ+ (resp. σ-) a bang arc with constant control u = 1 (resp. u = -1) and σs an admissible singular arc. We denote by σ1σ2 an arc σ1 followed by σ2. The surface Σ : HG(z) = 0 is called the switching surface and let Σ ⊂ Σ given by HG(z) = {HG, HF }(z) = 0. Let z(•) = (q(•), p(•)) be a reference curve on [0, t f ]. We note t → Φ(t) := HG(z(t)) the switching function, which codes the switching times.

Deriving twice with respect to time the switching function, one gets

Φ(t) = {HG, HF }(z(t)) , (1.13)
Φ(t) = {{HG, HF }, HF }(z(t)) + u(t) {{HG, HF }, HG}(z(t)).

(1.14)

From this, we derive.

Lemma 3. Assume that t is an ordinary switching time, that is Φ(t) = 0 and Φ(t) = 0. Then near z(t) every extremal projects onto σ+σ-if Φ(t) < 0 and σ-σ+ if Φ(t) > 0.

The situation is more complex for higher contact with Σ.

Definition 7. The fold case is characterized by Φ(t) = Φ(t) = 0 and Φ(t) = 0 (replacing u by ±1 in (1.14)). Hence z(t) ∈ Σ . Assume that locally Σ is a regular surface of codimension two. We have three cases:

• Parabolic case: Φ+(t) Φ-(t) > 0.
• Hyperbolic case: Φ+(t) > 0 and Φ-(t) < 0.

• Elliptic case: Φ+(t) < 0 and Φ-(t) > 0.

Denote by us(•) the singular control defined by (1.14) as Φ(t) = 0. In the hyperbolic case, through z(t), assuming {{HG, HF } , HG} (z(t)) = 0, there exists a singular arc, which is strictly admissible that is |us(t)| < 1. This arc is hyperbolic if {{HG, HF } , HG} (z(t)) > 0 (strict Legendre-Clebsch condition) and elliptic if this quantity is < 0. In the parabolic case, it can be absent or not admissible that is

|us(t)| > 1.
One has from [START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF],

σ+

Global case

It is based on [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF] and presented here for n = 3. It uses proposition 4. Under our assumptions (H0) -(H1), and if moreover the singular control is strictly admissible on [0, t f ], that is |us(t)| < 1, the reference singular arc can be immersed in a C 0domain up to the first conjugate time t1c and where the time minimal policy is of the form σ±σsσ±, see Fig. 1.2 for the interpretation of the first conjugate time t1c for the problem in the q-space with q = (c, v), v = u.

The curse of the non strict Legendre-Clebsch condition

More complicated and challenging situation is to analyze situations, where the strict Legendre-Clebsch condition is not satisfied that is {{HG, HF } , HG} (z(t)) vanishes at some times. Due to the complexity and in relation with our application, we shall assume that n = 3.

In this case, one has

D = 0 since {{HG, HF } , HG} = p • [[G, F ] , G] and recall that D = det(G, [G, F ], [[G, F ] , G]).
The singular flow is defined by the vector field Xs(q) = F (q) -D (q) D(q) G(q). Using a time reparameterization, it is related to the one-dimensional foliation Fs : D(q)F (q) -D (q)G(q). Complexity of the analysis is due to non isolated singularities contained in the set D (q) = D(q) = 0.

In relation with our optimal problem with terminal manifold N of codimension one given by x = d, our singularity resolution will be concerned by analyzing arcs initiating from the set S : n • [G, F ] (q) = 0, where n = (1, 0, 0) is the normal vector to N . Singularities can be roughly classified into two types: local in relation with singularities of S and propagated along the singular flow, and Lagrangian singularities in relation with the concept of conjugate-focal points. This will be developed in the next section.

Geometric techniques to analyze the 2d-McKeithan Network

Notations and definitions

The system is written q = F (q) + u G(q), |u| ≤ 1, with q = (c, v), G = ∂ ∂v . The terminal manifold N is given by c1 = d. The problem is to determine small time synthesis near a given point q0 ∈ N , which can be identified to 0. More precisely one wants to classify the syntheses under generic assumptions near the terminal manifold in relation with the Lie algebraic structure of {F, G} at q0. Our approach developed in our series of articles [START_REF] Bonnard | Time minimal synthesis for planar systems in the neighborhood of a terminal manifold of codimension one[END_REF][START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF][START_REF] Bonnard | Time minimal control of batch reactors[END_REF][START_REF] Launay | The generic local structure of time-optimal synthesis with a target of codimension one in dimension greater than two[END_REF] is to use the construction of semi-normal form for the action of the pseudo-group G of local diffeomorphisms and feedback transformation u → -u (so that σ+ and σ-can be exchanged). Additionally recall that the pseudo-group G f formed by local diffeomorphims and feedback actions of the form u = α(q) + β(q)v leaves the singular flow invariant [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF]. These groups act on the jet space of F at zero, G being identified to ∂ ∂v and the terminal manifold N to c1 = d. Note that the problem is flat that is G is tangent to N . For the action of the pseudo-group G on the jet spaces of F we refer to [START_REF] Martinet | Singularities of smooth functions and maps[END_REF] (we shall work in the C k category, where k ≥ 1 is not precised) and the semi-normal form is related to a semi-algebraic stratification on the jet spaces of (F, G, N ). One has N : c1 = d and the initial state is such that c1(0) < d. Denote in general N ⊥ = {(q, p); p • v = 0, ∀v ∈ TqN } and let n be the outward normal to N so that n = (1, 0, . . . , 0) and N ⊥ = (q, n(q)). Let z = (q, p) be a BC-extremal on [t f , 0], t f < 0 and z(0) ∈ N ⊥ (this convention is used since we integrate backwards from the terminal manifold). The set of minimizing switching points can be stratified into stratum of first kind if the optimal curves are tangent and second kind if they are transverse. The splitting locus L is the set of points, where the optimal control is not unique and the cut locus C is the closure of the set of points, where the optimal trajectories loses its optimality, see [START_REF] Boltyanskii | Sufficient conditions for optimality and the justification of the dynamic programming method[END_REF][START_REF] Brunovsky | Existence of Regular Syntheses for General Problems[END_REF] for the introduction of those concepts in the frame of semianalytic geometry.

We introduce the following triplet (F, G, N ) in the C ω -category, with N = {f (q) = 0}, where f is a (local) submersion. Fixing q0 ∈ R n , we denote by j k F (q0), j k G(q0), j k f (q0) the respective k-jets, that is the Taylor expansions at order k. We say that (F, G, f ) has at q0 a singularity of codimension i if (j k F (q0), j k G(q0), j k f (q0)) ∈ Σi, where Σi is a semi-algebraic submanifold of codimension i in the jets space.

Taking q0 ∈ N with a singularity of codimension i, an unfolding is a C 0 change of coordinates near q0 such that small time minimal synthesis is described by a system q = F (q, λ) + u G(q, λ), |u| ≤ 1, q ∈ R n-m and λ is a (vector) parameter.

Local syntheses: general tools and classification

We shall present the main steps to compute the time minimal syntheses, restricting our study to the 3-dimensional case, but it can be clearly extended to the n-space with the concept of codimension [START_REF] Launay | The generic local structure of time-optimal synthesis with a target of codimension one in dimension greater than two[END_REF].

The system is written: q = F (q) + u G(q), |u| ≤ 1 and let q = (x, y, z) be the coordinates, G being identified to ∂ ∂z and G being tangent to N , which can be identified to x = 0 and n = (1, 0, 0) is the outward normal to N . Recall that (generic) singular trajectories are given by q = F (q) + us(q)G(q), us(q) = -D (q)/D(q) and they can be classified into: hyperbolic, elliptic, exceptional, see section 3.

The first step is to stratify the terminal manifold into:

• S : singular locus defined by {q ∈ N, n • [G, F ](q) = 0}, • E : exceptional locus defined by {q ∈ N, n • F (q) = 0}.
Note that since the problem is flat, n • G(q) = 0 if q ∈ N so that N ⊥ ⊂ Σ, where Σ is the switching surface. Generic case. The first case is when F (q0) and [G, F ] (q0) are independent and not tangent to N and the synthesis follows from Lemma 3. It is given by

σ+ if n • [G, F ] (q0) < 0 and σ-if n • [G, F ](q0) > 0 since the final point is a virtual switching point.
Generic hyperbolic singular case and the concept of focal points. Now we present the basis of our analysis, the so-called hyperbolic situation. In this case, based on [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF], a semi-algebraic normal form is constructed to obtain the corresponding local syntheses and it is extended (in the jet space) along the reference singular arc in order to obtain the concept of focal point based on the notions of turnpike and conjugate points presented in Section 3. Semi-normal form: one has q0 = 0 and we make the following assumptions:

• The tangent space to N is G(0), [G, F ] (0).
• The set of points δ, where [G, F ] is tangent to N is a simple curve passing through 0 and transverse to G. • D(0) and D (0) are non zero.

With those assumptions, through 0, there exists a simple BC-singular extremal σs transverse to N . One can choose local coordinates so that G is identified to ∂ ∂z , δ to the axis (Oy) can be identified to t → (t, 0, 0), whose image is the (Ox)-axis. Note that N is identified to x = 0.

The semi-normal form is constructed in a tubular neighborhood of the reference singular curve σs and the vector field F is developed in the jet space with respect to (y, z) since σs : t → (t, 0, 0) is the x-coordinate. One has the following semi-normal form

ẋ = 1 + a(x)z 2 + 2b(x)yz + c(x)y 2 + R1 ẏ = d(x)y + e(0)z + R2 ż = (u -us |σ (x)) + f (x)y + g(0)z + R3,
where a(x) = 0, e(0) = 0 and R1 (resp. R2, R3) are terms of order ≥ 3 (resp. ≥ 2) in (y, z) and us is the singular control.

Furthermore we make the following assumption:

• The reference arc is hyperbolic on [t f , 0] so that a(x) < 0.

According to [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF] the time minimal synthesis near q0 is of the form σ+σsσ-, see Fig. Local syntheses and the semi-bridge phenomenon.

We use [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF] and the Mathematica's program given in Appendix 6.1. One fixes a point q0 ∈ S and we make the explicit computations of the switching and cut loci near q0 using the jet expansion of F at q0, G and N being normalized respectively to ∂ ∂v and x = 0. The computations are in the semi-algebraic category, easily implementable using symbolic algorithms and the cut locus and switching loci are stratifiable.

Switching locus. We denote by K the set of ordinary switching points for BCextremals, K+ being a switching σ-σ+ and K-being a switching σ+σ-, while W, W+, W-corresponds respectively to switching points of optimal extremals. More precisely, near N the stratification of W is W = W1 ∪ W2 where W1 is the first kind stratum composed by the hyperbolic singular arcs and W2 = Ws ∪ W+ ∪ W-is the second kind strata, where Wε, ε ∈ {-1, 1} is composed by the ordinary switching points of the policy σ-εσε and Ws is the first switching point of the Bang-Bang-Singular's policy.

From [START_REF] Agrachev | Symplectic methods for optimization and control[END_REF] we have two situations describe by Fig. ponds to a non optimal situation and has to be rejected to determine W . Using the theory and computations of [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF] only the stratification of W localized near q0 can be computed.

Singular locus. Near q0, one determines the singular locus Γs restricting to admissible singular trajectories, which are optimal. In other words, the test is : hyperbolicity and strict admissibility |us| < 1.

Cut locus. The cut locus C is the set of points where optimality is lost. One part of the cut locus is formed by the splitting locus L, where two minimizers intersects. It can be stratified into strata corresponding between intersections σ+, σ-intersections between σ+σ-and σ-. . . Again such intersections are described in [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF].

The local syntheses.

It is based on the classification of section 3.3.

Generic case. Since G is tangent to N , every final point is a virtual switching point. Denoting Φ(t) the switching function on [t f , 0], with

Φ(0) = n • G(q0) one has: if Φ(0) := n • [G, F ](q0) < 0 (resp. > 0) the terminating arc is σ-(resp. σ+).
Codimension one case. We have different cases corresponding to the fold case, since in the parabolic case one must distinguish between a parabolic point corresponding either to a non admissible hyperbolic or elliptic arc. The cases are represented on Fig. 1.5-1.6.

The syntheses can be represented by foliations by 2d-planes. Note that the role of σ+ and σ-can be interchanged since the semi-normal forms are computed using the transformation u → -u.

Focal points.

The synthesis in the hyperbolic case can be extended to a tubular neighborhood of σs for t ∈ [t f , 0] introducing the concept of focal point as follows. By assumptions, W := δ (0) belongs to span{G(0), [G, F ](0)}. Let λ1, λ2 be two scalars such that W = λ1G(0) + λ2[F, G](0) and by assumption λ1 = 0. Denoting by Xs(q) := F (q) -D (q) D(q) G(q), we introduce the variational equation δq(t) = ∂Xs ∂q (σs(t)) δq(t)

and let W (•) be a solution on [t f , 0] such that W (0) = W . Definition 8. Let t 1f be the first time in [t f , 0] such that: det(W (t), G(σs(t)), F (σs(t))) = 0. Then t 1f is called the first focal point along σs.
As for the fixed end point problem we have the following: Proposition 5. Provided t ∈]t 1f , 0], then in a tubular neighborhood of σs the set S = {exp(tXs(σ))} is a smooth surface and the synthesis is given by Fig. 1.3. Algorithm: Besides the definition, which leads to compute focal point using a singular value decomposition, an equivalent computation is to determine t = t 1f such that W (t) becomes collinear to G(σs(t)).

Two codimension-two cases.

The situation is more intricate. It is analyzed using the semi-normal form constructed in [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF]. A model is

ẋ = 1 + a z 2 + α1 xy 2 + α2 yz 2 + α3 xz 2 ẏ = b z ż = c z -us(0) -usx x -usy y, (1.15)
where us is the singular control, us(0) ∈ {1, 3}, usx = ∂us ∂x (0), usy = ∂us ∂y (0) > 0 and a, b, c = 0. The switching function Φ(t) := p3(t) is developed at order 3 and factorized as tP (t), where P is polynomial of order two with two roots t1, t2, which determine the switching points.

If a > 0, we are in the elliptic situation and if a < 0 in the hyperbolic situation and the reference singular arc is identified to σs : t → (t, 0, 0) being not admissible. One has |us(0)| > 1 and we can assume us(0) > 1.

Figures presented in this section are obtained with the Mathematica's program of Appendix 6.1, we use the following values:

α1 = α2 = α3 = usy = c = b = 1.
The case a > 0 and us(0) = 3. We have the following three situations The case a < 0 and us(0) = 1. There are two generic cases described by Fig. 1.9-1.10. The two cases are discriminated by the existence or not of the singular arc in the transition. In the second case, the switching locus Σ has two strata : Σ = W+ ∪ Ws, where W+ corresponds to optimal policies σ-σ+ and Ws to policies σ-σ+σs, Σ is not

C N σ- σ+ K - N σ - σ + σ - σ - W - C σ + σ + I W- N σ+ σ- σ+
C 1 . The stratum Ws ∪ Γs is C 1 but not C 2 . σ - N σ + σs N σ - σ - σ - σ + σ + σ + σs σ - W + K + N σ - σ - σ + σ + σs y < 0 y = 0 y > 0

Beyond the strict generalized Legendre-Clebsch condition

In our previous work in the 90's we restricted our study to the case, where the strict generalized Legendre-Clebsch condition is satisfied, see [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF]. Now our program is to extend this analysis when this condition is not satisfied.

W + N σs σ + σ + σ + σ - σ - σ - σ - σ - σ - Ws σ - σ - σ - N σ + W + σ - σ - σ - σ + W + N

The semi-bridge phenomenon as transition between two saturations

A bridge in optimal control was introduced in [START_REF] Bonnard | Time minimal saturation of a pair of spins and application in magnetic resonance imaging[END_REF] as a policy of the form Bang-Singular-Bang-Singular, where the second bang is related as a saturation due to the violation of the strict Legendre-Clebsch condition. In our context a semi-bridge is interpreted as a path to a bridge due to saturation between optimal singular arcs.

A tutorial model for the semi-bridge. Consider the case

ẋ = 1 + a y -3c yz + c z 3 ẏ = z ż = u, |u| ≤ 1, (1.16)
where a, c = 0 are parameters. Lie brackets computations. We have

F = (1 + ay -3cyz + cz 3 ) ∂ ∂x + z ∂ ∂y , G = ∂ ∂z , [G, F ] = 3c(y -z 2 ) ∂ ∂x -∂ ∂y , [[G, F ], G] = -6cz ∂ ∂x , [[G, F ], F ] = a ∂ ∂x , D(q) = det(G, [G, F ], [[G, F ], G]) = -6cz, D (q) = det(G, [G, F ], [[G, F ], F ]) = a, D (q) = det(G, [G, F ], F ) = ay -2cz 3 + 1.

Stratification of N .

We have the following properties Observe that E : n • F (q) = 0 is the defined by y(a -3cz) + cz 3 + 1 = 0.

• N : x = 0, G = ∂ ∂z • S : n • [G, F ] = 0 ∩ {x = 0} is the parabola: y = z 2 . • n • [[G, F ] , G] (0) = 0, n • [[G, F ] , F ] (0) = 0.
One will localize our study in a neighborhood V of 0 such that: E ∩ V = ∅. We represent such a situation on Fig. 1.11 One has us = a/(6cz) and the points usat are given by |us| = 1.

We have two saturating points, one in the hyperbolic domain and one in the elliptic domain. Near the fold, the two saturation phenomenons are glued together using the curvature of S, the normal to N being given by n = (1, 0, 0), while q0 follows S. Note that since the control is blowing up at the fold, we encounter the case a > 0, us(0) = 3 (see Fig. 1.7) and the bifurcation phenomenon of C, which splits into the case σ+, σ-and the case σ+σ-and σ-intersecting minimizers. The two strata will be denoted C12 and C1. The model detects the two strata.

The stratification of the terminal manifold near the semi-bridge and the local time minimal synthesis are described in Fig. 1.12. Local syntheses by symbolic computations. We use the Mathematica's program in Appendix 6.1 to derive the local syntheses for the semi-bridge model.

usat usat y v u = -1 u = +1 S (iv) (i) (ii) (vi) (v) (i) :
We compute the time minimal synthesis in a neighborhood U of q0 ∈ N where N : x = 0 is the target. Since the problem is flat, q0 is either an ordinary switching point, a fold point, or a codimension two case.

Denote n the unit normal of N at q0 such that n belongs to the half-space containing the set {X + u Y, |u| ≤ 1}, i.e. n = (1, 0, 0). Assume X(q0) and [Y, X](q0) are not tangent to N . Then q0 is an ordinary point and the optimal control is given by u

(q0) = -sign p • [Y, X](q0).
We take q0 = (0, s 2 0 , s0) and we look at the behavior of BC-extremals reaching N near q0. The switching surface W , the splitting locus C and the trajectories σ± are computed via symbolic computations, expanding the jets of F and G up to some order.

• Near a parabolic point on an hyperbolic singular arc.

BC-extremals σ-switch while BC-extremals σ+ don't. These computations are represented in Fig. 1.13 for q0 = (0, (-0.05) 2 , -0.05). • Near a parabolic point on an elliptic singular arc.

BC-extremals σ± reflect on the switching surface W±, hence a cut locus appears between the switching surfaces. We represent in Fig. 1.14 with q0 = (0, 0.04 2 , 0.04), 0.04 > zsat = a/(6c) this cut locus and the switching surface W-, together with the elliptic singular surface Γs. • Near a saturated point on an hyperbolic singular arc.

We have q0 = (0, z 2 sat , -zsat) where zsat = a 6c . The synthesis is represented in Fig. 1.15 and it corresponds to the synthesis described earlier by Fig. 1.9. Remark 1. We can investigate other singularities by the same method, for instance:

• Cusp singularity y = z 3 . • Singularity y 2 = z 2 . • Case D = D = 0 is q = {0}.

Global analysis.

Take u = ε ∈ {-1, 1}. The dynamics (1.16) is integrable and we have

p3(t) = 1 2 t -2ct 2 + (a -6cεs0)t -6cs 2 0 + 6cw0 .
Denoting by Γε(q0), ε ∈ {-1, 1} the surface composed by the trajectories σε passing through q0, then Γε(0, w0, s0) is parameterized by

x =t (a -3cs0) + cs 3 0 + 1 + 1 2 t 2 s0(a -3cs0) + 3cε(s 2 0 -w0) + o(t 3 ) y =w0 + s0t + εt 2 2 + o(t 3 ) z =s0 + εt + o(t 3 ).
Solving p3 = 0 with respect to w0, the switching surface Kε for BC-extremals of controle u = ε on N is parameterized by:

x = 1 6 t 2 - a 2 c + 6aεs0 + 6as0 -18cεs 2 0 -9cs 2 0 + 1 6 t 6as 2 0 -12cs 3 0 + 6 , y =t (ε + 1)s0 - a 6c + 1 6 (3ε + 2)t 2 + s 2 0 z =s0 + εt + o(t 3 ).
The switching condition of a trajectory on Wε is given by the sign of [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF]. For -zsat < z < zsat, the optimal policy is σ+σ-. The case z = -zsat was given in Fig. 1.15 and is also represented here with the surface Ws. Construction of a semi-normal form and the need of integrability assumption.

Notations and assumption: We normalize in a small neighborhood of N . We assume that F is transverse to N . Let q = (x, y, z), q0 = (0, 0, 0), G = ∂ ∂z and N := {x = 0}. We denote by σ0 the reference trajectory of F identified to t → (t, 0, 0) so that F |σ 0 = ∂ ∂x . Assuming x small enough one can write the system as:

ẋ = 1 + f (y, z) + R1 ẏ = g(y, z) + R2 ż = h(y, z) + R3,
where Ri(x, y, z) = o(|x|) and the corresponding model is

ẋ = 1 + f (y, z) ẏ = g(y, z) ż = h(y, z) + u.
Singular flow on the model: Furthermore, we assume ∂g ∂z = 0. Since the singular flow is feedback invariant, setting z = g(y, z), z → z and using a proper feedback the system takes the form

ẋ = 1 + f (y, z) ẏ = z ż = u.
Computing, we get:

[G, F ] = -fz ∂ ∂x + ∂ ∂y [[G, F ] , G] = -f z 2 ∂ ∂x , [[G, F ] , F ] = -fyz ∂ ∂x + fy ∂ ∂y .
(1.17)

Hence D = -f z 2 , D = -fyfz + fyz and us = - D D = fyfz -fyz f z 2 .
Therefore, we have dy dz = z us and assuming the separability condition

dy dz = h1(y)h2(z)
one can integrate the singular trajectories such that q(0) ∈ S as follows, parameterizing by z instead of t using the relation Assuming that y(0) can be expressed as y(0) = p1(z(0)) and moreover that the left-hand side is invertible, one gets y = p2(z, z(0)). This defines the singular leaf F through q(0) ∈ S using the equation

dx dz = 1 + f (y, z) z .
Using x(0) = 0, one gets:

x(z, z(0)) = z z(0)

(1 + f (p2(z, z(0)), z)) dz.

Remark 2. This requires an integrability condition to obtain the leaf F. If it is not satisfied, integration is obtained through a proper numeric integration.

Application.

Singular set. Back to the tutorial model (1.16), we fix a = c = 1, we are in the separable case and using our previous algorithm one gets Proposition 6. 1. Integrating the singular flow, we get:

x(t, z(0)) = t + 2ε(2a+3ck) (kt+z(0) 2 ) 5/2 -10c(kt+z(0) 2 ) 3 t 15k 2 y(t, z(0)) = 2ε 3k kt + z(0) 2 3/2 z(t, z(0)) = ε kt + z(0) 2 1/2 ,
where ε = ±1, k = 12c/a. 2. Parameterizing by z, the integral singular leaf F with q(0) ∈ S is given

y(z, z(0)) = 2c (z 3 -z(0) 3 )/a + z(0) 2 x(z, z(0)) = g(z) -g(z(0)),
where g(z) = 3cz 2 5a 2 2a(2c + 1) z 3 + 30cz(0) 2 (a -2cz(0)) z + 30c 2 z 4 + 5a(az(0) 2 -2cz(0) 3 + 1) .

We represent on Fig. 1.20 the stratification of the singular leaf F with a = c = 1 using hyperbolic and elliptic cases and the admissibility conditions |us| ≤ 1.

The 2d-McKeithan scheme

Recall that the network is T

+ M A B k 1 k 2 k 3 k 4
and the system q = F (q)+u G(q),

using the coordinates q = (x, y, v), x = [A], y = [B], v = k1 and reduced to the stoichiometric class T + M + A = δ1 and T + M + B = δ2 takes the form ẋ = -β2xv α 2 -β3xv α 3 -δ3v (x + y) + δ4v + v (x + y) 2 ẏ = β2xv α 2 -β4yv α 4 with 0 ≤ x ≤ δ1, 0 ≤ y ≤ δ2, δ3 = δ1 + δ2, δ4 = δ1 δ2, k2 = β2 v α 2 , k3 = β3 v α 3 , k4 = β4 v α 4 .

Lie brackets computations

We have:

• F = -β2xv α 2 -β3xv α 3 -δ3v(x+y)+δ4v+v (x+y) 2 ∂ ∂x + β2xv α 2 -β4yv α 4 ∂ ∂y , • G = ∂ ∂v , • [G, F ] = x(α2β2v α 2 -1 + α3β3v α 3 -1 + δ3) + δ3y -δ4 -x 2 -2xy -y 2 ∂ ∂x + α4β4yv α 4 -1 -α2β2xv α 2 -1 ∂ ∂y , • [[G, F ] , G] = x(α2β2(α2 -1)v α 2 -2 + α3β3(α3 -1)v α 3 -2 ) ∂ ∂x + yα4β4(α4 - 1)v α 4 -2 -xα2β2(α2 -1)v α 2 -2 ∂ ∂y , • [[G, F ] F ] = -xv α 2 β2δ3(α2 -1) -yv α 2 β2δ3(α2 -1) + v α 2 (α2β2δ4 + xy(2α2β2 - 2β2) -β2δ4) + x 2 (α2β2 -β2)v α 2 + y 2 (α2β2 -β2)v α 2 + yv α 3 (β3δ3 -α3β3δ3) + v α 3 (α3β3δ4 -β3δ4)+x 2 (β3 -α3β3)v α 3 +y 2 (α3β3 -β3)v α 3 +yv α 4 (α4β4δ3 -β4δ3)+ xy(2β4 -2α4β4)v α 4 +y 2 (2β4 -2α4β4)v α 4 ∂ ∂x + x(v α 2 +α 3 -1 (α2β2β3 -α3β2β3)+ v α 2 +α 4 -1 (α4β2β4 -α2β2β4)) + xv α 2 (α2β2δ3 -β2δ3) + yv α 2 (α2β2δ3 -β2δ3) + v α 2 (-α2β2δ4+xy(2β2-2α2β2)+β2δ4)+x 2 (β2-α2β2)v α 2 +y 2 (β2-α2β2)v α 2 ∂ ∂y .

Singular arcs

One has:

• D(q) = ((α4 -1)α4β4yv α 4 -3 -(α2 -1)α2β2xv α 2 -3 )(α2β2xv α 2 + α3β3xv α 3 + δ3v(x + y) -δ4v -vx 2 -2vxy -vy 2 ) + x(α 2 2 β2v α 2 -α2β2v α 2 + (α3 - 1)α3β3v α 3 )(α2β2xv α 2 -3 -α4β4yv α 4 -3 ), • D (q) = β2v α 2 -2 (α2β2xv α 2 + α3β3xv α 3 + δ3v(x + y) -δ4v -vx 2 -2vxy - vy 2 )(α2β3xv α 3 -α2β4xv α 4 + (α2 -1)δ3v(x + y) + δ4(v -α2v) -α2vx 2 - 2α2vxy -α2vy 2 -α3β3xv α 3 + α4β4xv α 4 + vx 2 + 2vxy + vy 2 ) + (α2β2xv α 2 -1 - α4β4yv α 4 -1 )((α2 -1)β2v α 2 (δ4 -(x + y)(δ3 -x -y)) + (α3 -1)β3v α 3 (y(y -δ3) + δ4 -x 2 ) + (α4 -1)β4yv α 4 (δ3 -2(x + y))), • D (q) = (β2xv α 2 -1 -β4yv α 4 -1 )(α2β2xv α 2 + α3β3xv α 3 + δ3v(x + y) -δ4v -vx 2 - 2vxy -vy 2 ) -(α2β2xv α 2 -1 -α4β4yv α 4 -1 )(β2xv α 2 + β3xv α 3 + δ3v(x + y) -δ4v - vx 2 -2vxy -vy 2 ),
and the singular control is given by: us = -D (q)/D(q).

Classification of local syntheses for the McKeithan network

We consider the case max [A], with x = [A]. We proceed as follows.

Stratification of the terminal manifold: x = d.

Singular locus

S : n • [G, F ](q) = 0 and x = d with n = (1, 0, 0). It is given by:

S : α2β2dv α 2 -1 + α3β3dv α 3 -1 + dδ3 -δ4 -d 2 + y(δ3 -2d) -y 2 = 0. (1.18)
Denoting by ∆ the discriminant of the polynomial function y → n • [G, F ](q) ∩ x = d, a singularity can occur for ∆ = 0. One has Lemma 4. Assume αi, βi, δi > 0, i=1,2 and d, v > 0. Then we have ∆ = (δ1 -δ2) 2 + 4d (α2β2v α 2 -1 + α3β3v α 3 -1 ) > 0 so that there is no ramification and S contains at most two real positive branches.

Definition 9. A semi-bridge occurs at a point q ∈ S if n • [[G, F ] , G] (q) = 0.
Computing, a semi-bridge occurs if

v α 3 -α 2 = - (α2 -1)α2β2 (α3 -1)α3β3 . (1.19)

Exceptional locus

It is given by E : n • F (q) = 0 and x = d. Computing, one gets:

E : -β2dv α 2 -β3dv α 3 + d 2 v -dδ3v + y (2dv -δ3v) + δ4v + vy 2 = 0.
(1.20)

The discriminant of the polynomial n • X(q) in y is ∆ = v(4d(β2v α 2 + β3 v α 3 ) + v(δ1 -δ2) 2 ) > 0 and E contains at most two real positive branches. Fig. 1.21 gives a picture of stratification of x = d for the McKeithan system with a focus where S is folded. We represent the sets S and E and the stratification of S in hyperbolic, elliptic and parabolic points. From this, we deduce the admissible points. We are in the flat case and a point of N is either an ordinary point or a fold point. For an ordinary point, the final optimal control is regular and is equal to -sign Φ(0). At a fold point on S, the final optimal control may be singular (see Definition 7).

The optimal syntheses near the fold can be described using the techniques of the tutorial model (1.16) and symbolic computations.

Conclusion

In this article, we have presented the general tools to analyze a Mayer problem to optimize the yield of chemical networks. Using the McKeithan network, we have completed the analysis of [START_REF] Bonnard | Time minimal control of batch reactors[END_REF] for the simple scheme A → B → C.

The main point is to extend a geometric techniques from [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF] to consider the case when the strict Legendre-Clebsch is not satisfied. The starting point is to analyze the semi-bridge phenomenon.

Note that with the same techniques we can analyze the situations up to codimension 2 describe in details in [START_REF] Launay | The generic local structure of time-optimal synthesis with a target of codimension one in dimension greater than two[END_REF] (in particular the exceptional locus E). Furthermore, it can be extended to an important problem of codimension 3 occurring for reversible chemical networks and the existence of equilibria, see Fig.

Fig. 1 . 1 .

 11 Fig. 1.1. Fold case in the elliptic 2D-situation and the antiturnpike phenomenon.

Fig. 1 . 2 .

 12 Fig. 1.2. (left) Time minimal policy in the q-space. (right) Conjugate point in the c-space.

4. 1

 1 Time minimal syntheses near the terminal manifold 4.1.

  Fig. 1.3.

1 . 4 .Fig. 1 . 4 .

 1414 Fig. 1.4. (left) Crossing. (right) Reflecting.

Fig. 1 . 6 .

 16 Fig. 1.5. Fold case.

Fig. 1 . 7 . 0 Fig. 1 . 8 .

 17018 Fig. 1.7. Bifurcation of the cut (a > 0, us(0) = 3). The cut locus C splits into C1 ∪ C12 where C1 : σ-∩ σ+ and C12 : σ+σ-∩ σ-.

Fig. 1 . 9 .

 19 Fig. 1.9. Saturating case 1: -1 = a < 0, 1 = usx > 0. Top figures are obtained using the Mathematica's program given in the Appendix 6.1 and these situations are sketched in the bottom figures.

Fig. 1 . 10 .

 110 Fig. 1.10. Saturating case 2: -1 = a < 0, 1 = usx < 0. Top figure is obtained using the Mathematica's program given in the Appendix 6.1 and the corresponding situations are sketched in the bottom figures.

Fig. 1 . 11 .

 111 Fig. 1.11. Stratification of N for the system (1.16) (a = c = 1). Green dotted line: elliptic, red line: hyperbolic, crosses: saturating values of the singular control.

Fig. 1 . 12 .

 112 Fig. 1.12. Stratification of the target x = 0 for a semi-bridge model and the associated local syntheses.

Fig. 1 . 13 .

 113 Fig. 1.13. Time minimal synthesis near an hyperbolic fold point for q0 = (0, (-0.05) 2 , -0.05). Only the BC-extremals with optimal control -1 switch for the tutorial model (1.16) with a = c = 1. The tangent space of the switching surface at q0 is represented.

Fig. 1 . 14 .

 114 Fig. 1.14. Time minimal synthesis near an elliptic fold point for q0 = (0, 0.04 2 , 0.04) for the tutorial model (1.16) with a = 1, c = 5.

Fig. 1 . 15 .

 115 Fig. 1.15. Time minimal synthesis near a saturating hyperbolic fold point q0 = (0, z 2 sat , -zsat) for the tutorial model (1.16) with a = c = 1. Trajectories σ-switch, while σ+ don't. The surface Ws is the first switching point of the sequence σ+σ-σs.

,

  and this determinant is equal to as 2 0 -2cs 3 0 + 1 (6cεs0 -a)/(6c). In the case a = c = 1, we get that σ+ (resp. σ-) switches on W+ (resp. W-) for z ∈]zsat, 1[ (resp. z ∈] -zsat, 1[) as illustrated in Fig.1.16 (resp. Fig.1.17).

Fig. 1 . 16 .

 116 Fig. 1.16. The trajectories σ+ starting from q0 = (0, w0, s0) such that w0 < s 2 0

Fig. 1 . 17 .

 117 Fig. 1.17. The trajectories σ-starting from q0 = (0, w0, s0) such that w0 > s 2 0

Fig. 1 .

 1 Fig. 1.18. (left) σ-reflect on W-for z > zsat. (right) At z = a 18c < zsat such that us(z) = 3, σ-cross W-and the cut locus disappears (see Fig.1.7). For -zsat < z < zsat, the optimal policy is σ+σ-. The case z = -zsat was given in Fig.1.15 and is also represented here with the surface Ws.

1

 1 Fig. 1.18. (left) σ-reflect on W-for z > zsat. (right) At z = a 18c < zsat such that us(z) = 3, σ-cross W-and the cut locus disappears (see Fig.1.7). For -zsat < z < zsat, the optimal policy is σ+σ-. The case z = -zsat was given in Fig.1.15 and is also represented here with the surface Ws.

  Fig. 1.18. (left) σ-reflect on W-for z > zsat. (right) At z = a 18c < zsat such that us(z) = 3, σ-cross W-and the cut locus disappears (see Fig.1.7). For -zsat < z < zsat, the optimal policy is σ+σ-. The case z = -zsat was given in Fig.1.15 and is also represented here with the surface Ws.

Fig. 1 . 19 .

 119 Fig. 1.19.For z > zsat σ+ reflect on W+, hence we have a cut locus. For z < -zsat, the optimal policy is Bang-Singular.

Fig. 1 . 20 .

 120 Fig. 1.20. Stratification of the singular set Γs for the system (1.16) with parameters a = c = 1.

Fig. 1 . 21 .

 121 Fig. 1.21. Stratification of the surface x = d for the McKeithan reaction. Dotted line: elliptic, red line: hyperbolic.

  .m --Tested with Mathematica v11.3 Linux *) (* This program computes the singular set, switching surface \ and splitting locus to derive the time minimal synthesis in a neighborhood of N. *) (* Lie derivative of w wrt v *) LieDerive[v_List,w_List,var_List]:=Transpose[D[w,#]& /@ var] . v /; \ Length[v] == Length[var] (* Lie bracket [v,w] *) LieBracket[v_List,w_List,var_List]:=\ LieDerive[w,v,var] -LieDerive[v,w,var] /; \ YXY = LieBracket[YX,Y,{x,y,v}]; YXX = LieBracket[YX,X,{x,y,v}]; DD = Det[{Y,YX,YXY}]; DDp = Det[{Y,YX,YXX}]; DDpp = Det[{Y,YX,X}]; aux = Solve[(YX[[1]])==0,{y}]; Sy[v_] = Expand[(y /. aux[[1]]) /. x->0] (*.......... Gammas ..........X,Y,us,q0p0,ord]; ps = {t0val,w0val,s0val}; gms = mTaylor[#,{t,w0,s0},ps,ordG]& /@ qpt[[1;;3]] //Expand gmsall = mTaylor[#,{t,w0,s0},ps,ordG]& /@ qpt //Expand (*.......... Gamma+,Gamma-..........*) u = ep; q0p0 = {0,w0,s0,1,0,0}; pt = {t0val,w0val,s0val}; qpt = Refine /@ IntSmallTime[X,Y,u,q0p0,ord]; qpt = mTaylor[#,{t,w0,s0},pt,ordG] & /@ qpt; gmall = Refine /@ qpt //Refine gmp = gmall[[1;;3]] /. {ep->1} gmm = gmall[[1;;3]] /. {ep->-1} (*.......... W+,W-..........*) ord = 3; qpt = Refine /@ IntSmallTime[X,Y,u,q0p0,ord]; eq = mTaylor[qpt[[6]],vars,{t0val,w0val,s0val},ord] // Simplify aux = Solve[eq==0,vw]; sols = Flatten[vw /. aux]; sols = Refine[Expand[mTaylor[#,{t,w0,s0},{0,w0val,s0val},ord] & /@ sols]] // Simplify solp = mTaylor[ # /. {ep->1},vars,{t0val,w0val,s0val},ordG] & /@ sols; solm = mTaylor[ # /.{ep->-1},vars,{t0val,w0val,s0val},ordG] & /@ sols qptp = (qpt /. ep->1 /. vw-># ) & /@ solp; qpwps = mTaylor[#,{t,w0,s0},pt,ordG] & /@ qptp // Flatten; qptm = (qpt /. ep->-1 /. vw-># ) & /@ solm; qpwms = mTaylor[#,{t,w0,s0},pt,ordG] & /@ qptm //Flatten (*.......... Ws ..........@ IntSmallTime[X,Y,us,q0p0,ord]; qpt = mTaylor[#,{t,w0,s0},pt,ord] & /@ qpt; q0p02 = qpt /. t->ss; qpt2 = Refine /@ IntSmallTime[X,Y,1,q0p02,ord]; qpt2 = mTaylor[#,{t,w0,s0,ss},{0,w0val,s0val,0},ord] & /@ qpt2; p3t2 = qpt2[[6]] aux = Solve[p3t2==0,ss]; sols = Flatten[ss /. aux]; sols = Refine[Expand[mTaylor[#,{t,w0,s0},{0,w0val,s0val},ord] & /@ sols]] qpws = qpt2 /. {ss-> #} & /@ sols; qpws = mTaylor[#,{t,w0,s0},pt,ordG] & /@ qpws //Flatten ws = Refine /@ qpws[[1;;3]] (*.......... C1 ..........*) ord = 3; q0p0 = {0,w0,s0,1,0,0}; pt = {0,w0val,s0val}; qpt = IntSmallTime[X,Y,u,q0p0,ord]; qpt = mTaylor[#,{t,w0,s0},pt,ord] & /@ qpt; (* u=-1 *) {xp,yp,zp} = qpt[[1;;3]] /.{ep->-1}; (* u=+1 *) {xm,ym,zm} = qpt[[1;;3]] /.{ep->1,s0->s0p,w0->w0p}; eqns = {xp==xm,yp==ym,zp==zm}; aux = Solve[eqns,{t,w0p,s0p}] // Expand; aux2 = {t,w0p,s0p} /. aux; c1sols = mTaylor[#,{t,w0,s0},pt,ord] & /@ aux2; (*.......... C12 ..........*) ord = 2; q0p0 = {0,w0,s0,1,0,0}; pt = {0,w0val,s0val}; qpt = IntSmallTime[X,Y,u,q0p0,ord]; qpt = mTaylor[#,{t,w0,s0},pt,ord] & /@ qpt; vars = {t,ss1,ss2,w0,s0,w0p,s0p}; pt2 = {0,0,0,w0val,s0val,w0val,s0val}; (* u=-1 *) {x1,y1,z1} = qpt[[1;;3]] /.{ep->-1}; (* u=-1 & switch *) ptfin = (qpt /.{ep->-1,w0->w0p,s0->s0p})/.t->ss1; qpt2 = IntSmallTime[X,Y,1,ptfin,ord] /. t->ss2; {x3,y3,z3} = mTaylor[#,vars,pt2,ord] & /@ qpt2[[1;;3]]; eqns = {x1==x3,y1==y3,z1==z3,ss2+ss1==t,ptfin[[6]]==0}; aux = Solve[eqns,{t,ss1,ss2,w0p,s0p}]; aux2 = {t,ss1,ss2,w0p,s0p} /. aux; c12sols = mTaylor[#,vars,pt2,ord] & /@ sols;

  and p(t f ) satisfies the transversality condition. An extremal control is called regular if |u(t)| ≤ 1 a.e. and singular if ∂H ∂u = 0 everywhere. An extremal is said exceptional if M = 0. A regular extremal control is called "bangbang" if u(•) is piecewise constant on [0, t f ] (i.e. the number of switches is finite).

Length[v] == Length[w] == Length[var]

(* Poisson bracket {f,g} *) PoissonBracket[f_, g_,q_List,p_List] /; Length[q] == Length [p] (*--------------------------------------------------------*) (*--------------------------main--------------------------*) (*--------------------------------------------------------*