T Bakir

B Bonnard

J Rouot

Geometric Optimal Control Techniques to Optimize the Production of Chemical Reactors using Temperature Control

Keywords: Mass action chemical systems, Zero deficiency theorem, Pontryagin maximum principle, Geometric optimal control

The dynamics of mass reaction kinetics chemical systems is modeled by the Feinberg-Horn-Jackson graph and under the "zero deficiency assumption", the behavior of the solutions is well known and splits into two cases: if the system is not weakly reversible there exists no equilibrium, nor periodic solution and if the network is weakly reversible in each stoichiometric subspace there exists only one equilibrium point and this point is asymptotically stable. By varying the temperature, one gets a single input control system and in this article we study the problem of maximizing the production of one species during the batch time. Our aim is to present the geometric techniques and results based on the Pontryagin maximum principle to compute the closed loop optimal solution. The complexity of the problem is illustrated by using two test bed examples: a sequence of two irreversible reactions and the McKeithan scheme.

Introduction

Important developments in optimal control come from the Pontryagin maximum principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF] which combined later (in the eighties) with the geometric techniques led to seminal theoretical results, see for instance [START_REF] Agrachev | Control Theory from the Geometric Viewpoint[END_REF] and [START_REF] Jurdjevic | Geometric control theory[END_REF]. They were applied to many industrial problems in space mechanics, quantum control [START_REF] Bonnard | Optimal Control with Applications in Space and Quantum Dynamics[END_REF] and an important modern application is the control of biological and chemical systems [START_REF] Schättler | Geometric optimal control. Theory, methods and examples[END_REF].

The objective of the article is to present the techniques of this area to analyze the problem of maximizing the production of one chemical species [X] of a chemical network where the sequence of reactions occurs in a batch chemical reactor and the control is the temperature. Thanks to the Maximum Principle it can be formulated as minimizing the batch time t f to produce a fixed amount [X] = d and this leads to a time minimal control problem for a single input control system: dc dt = f (c, T ), where c is the concentrations vector, T is the temperature, with c(0) = c 0 and c(t f ) ∈ N, N being the terminal manifold: {[X](t f ) = d}. The dynamics model under the mass action kinetics assumption produces, at constant temperature T , a polynomial system, which can be deduced from the Feinberg-Horn-Jackson graph [START_REF] Tóth | Reaction kinetics: exercises, programs and theorems[END_REF]. Thanks to the so-called seminal zero deficiency theorem obtained by the contributions of Feinberg, Horn, Jackson in the 70 s, the behaviors of solution can be split in many applications into two cases. The first case is the non weakly reversible case, and there exists no equilibrium point nor periodic trajectory, an example being a sequence of consecutive reactions A → B → C . . .. The second case is the weakly reversible case, where under a mild additional assumption, in each stoichiometric class there exists an unique equilibrium point c * which is asymptotically stable, see [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF][START_REF] Horn | General mass action kinetics[END_REF]. An example of this case being the so-called McKeithan scheme:

T + M A B .
If the dynamics is simple at fixed temperature, the optimal control by varying the temperature can be very complicated as illustrated by the following tutorial example. The rate of each reaction using Arrhenius law depends upon k(T ) = A e -E/(RT ) where A, E are parameters, the second E being the energy of activation. Consider the following scheme

A B C k 1 k 2
of a sequence of two first order irreversible reactions. Denoting by c 1 := [A], c 2 := [B] the respective concentrations, the dynamics takes the form:

ċ1 = -k 1 c 1 , ċ2 = k 1 c 1 -k 2 c 2 , k i = A i e -E i /(RT ) , i = 1, 2.
Since T → v(T ) k 1 (T ) is a bijection, v can be chosen as the control variable. The dynamics is irreversible and is very simple (note it is a consequence of the zero deficiency theorem) and integrating the dynamics at constant temperature leads to two situations only: k 1 k 2 or k 1 = k 2 (resonant case). But assume that we want by controlling the temperature to maximize the ratio: z = c 2 /c 1 . This leads to maximize z(t f ) where [0, t f ] is the batch duration. Using the dynamics:

ż = v -βv α z + vz, α = E 2 /E 1 , β = A 2 /A α 1 with z(0) = [B]/[A] |t=0
which can be taken as z(0) = 0 and v ∈ [v m , v M ] is associated to T ∈ [T m , T M ]. Denoting H = ż (which plays the role of Hamiltonian in the optimal problem) the problem leads to maximize H over v ∈ [v m , v M ] in order to maximize z(t f ). Clearly one gets three cases if α > 1 (and one case if α < 1): the maximum is depending upon z and can be either v = v m , v = v M or an intermediate value v = v s called singular in optimal control defined by: ∂H ∂v = 1 -αβv α-1 z + z = 0. In this case a simple graph analysis of H tells us that an optimal policy can be of the form σ M σ s σ m where σ M , σ s , σ m being arcs associated to v M , v s , v m respectively, σ 1 σ 2 representing an arc σ 1 followed by σ 2 . This corresponds on [0, t f ] to the concatenation of the controls v M , v s , v m .

A similar analysis in the case of α < 1 leads to the optimal control: v = v M over [0, t f ], associated to maximal temperature T M .

From this example we deduce two facts. The first one is that the optimal policy depends upon α = E 2 /E 1 v.s. k 1 /k 2 for the dynamics at constant T . Second, the discussion depends upon the competition between the Hamiltonians H associated to v m , v M , v s , to provide the maximum of H. This second point will give the geometric frame to perform the analysis in relation with singularity theory: competition between Hamiltonians dynamics to compute the optimal solutions, see the earliest seminal contributions [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]. In particular, the complexity of the problem is directly illustrated by the two test bed cases.

The organization of the article is the following. In section 2, the mathematical model is recalled based on the Feinberg-Horn-Jackson graph and the zero deficiency theorem is precisely stated [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF][START_REF] Horn | General mass action kinetics[END_REF] and we present the two test bed cases, in particular the McKeithan network [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] Mckeithan | Kinetic proofreading in T-cell receptor signal transduction[END_REF]. The section 3 is devoted to the Pontryagin maximum principle and the classification of the solutions into regular and singular solutions. The general properties are discussed, based on [START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF] in the regular case and [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 the singular case. This leads to compute time optimal 2 sequences useful in our study, based on basic concepts of singularity theory, in particular seminormal forms [START_REF] Martinet | Singularities of smooth functions and maps[END_REF]. The section 4 presents the main technical tool to handle the problem that is singularity theory and classification of the optimal syntheses based on [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]. This is applied in section 5 to discuss the optimal syntheses for the two test bed cases. In the final section we discuss in conclusion the mathematical and computational obstacles to complete our study.

Mathematical model and the zero deficiency theorem

In this section, we make a brief presentation of the model with main properties based on the seminal works of Feinberg-Horn-Jackson in the 70's. For more details, references are the earliest articles [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF][START_REF] Horn | General mass action kinetics[END_REF] and [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] Anderson | Global asymptotic stability for a class of nonlinear chemical equations[END_REF] for more recent contributions.

Mass action kinetics and the Feinberg-Horn-Jackson graph

We consider a set of m chemical species {X 1 , . . . , X m }, and the state of the dynamics is the vector c = (c 1 , . . . , c m ) ∈ R m ≥0 representing the molar concentrations. Let R be a set of reactions, each reaction being denoted y → y and is of the form:

m i=1 α i X i → m i=1 β i X i ,
where α i , β i are the stoichiometric coefficients and the vectors y = (α 1 , . . . , α m ) and y = (β 1 , . . . , β m ) are forming the vertices of the so-called Feinberg-Horn-Jackson oriented graph associated to the network, edges being oriented according to y → y . Each reaction y → y 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 , (1)
c y = c 1 α 1 . . . c m α m (2)
and

k(T ) = Ae -E/(RT ) (3)
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 another rate formulae can be used in particular to deal 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]). One can label the sequences of reaction R by i = 1, . . . , ñ and this is defining a set of (increasing) functions of the temperature denoted k i (T ), i = 1 . . . ñ, each defined by a set of parameters (A i , E i ) forming the set Λ.

The dynamics of the system, taking into account the whole network is:

ċ(t) = f (c(t), T ) = y→y ∈R K(y → y ) (y -y). (4)
Note that if the reactor is fed at some constant rate r, the more general dynamics takes the form ċ

(t) = f (c(t), T ) + r (c 0 1 , . . . , c 0 m ) -c(t) . (5)
where (c 0 1 , . . . , c 0 m ) are the fixed concentrations of the feeding reactor. This led to a more general system which can be set in an unique frame of chemical network (introducing additional reaction), but in this article, from the purpose of control theory, we shall consider only the case (4) which is called the batch (or closed) case.

More explicit representation of the dynamics and the zero deficiency theorem

First of all, we introduce the so-called stoichiometric subspace S = span{yy , y → y ∈ R} and the set (c(0) + S ) ∩ R m ≥0 are called the (strictly if > 0) positive stoichiometric compatibility classes. From (4) it is clear that S is invariant for the dynamics but more precisely one has.

Lemma 1 ([2]

). Let c(t) be a solution of (4) with initial condition c(0

) ∈ R m >0 . Then c(t) ∈ R m >0 for all t > 0.
Before to state our theorem the dynamics has to be rewritten using the following concepts.

• Having labeled the set of vertices by i = 1, • • • , n, whose corresponding stoichiometric vectors (y 1 , . . . , y n ), the complex matrix is Y = (y 1 , . . . , y n )

• The incidence connectivity matrix: A (a i j ) contains the Arrhenius coefficients k i of the reactions using the rule: k 1 = a 21 indicating a reaction with constant k 1 from the first node to the second, that is

y 1 y 2 k 1 .
With the mass kinetics assumption (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] for extension) the dynamics (4) can be expressed as

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

à = A -diag        n i=1 a i1 , . . . , n i=1 a in        (7)
and we denote:

c Y = (c y 1 , . . . , c y n ) . (8
            ċ1 . . . ċN             =                             -k 1 0 • • • • • • 0 k 1 -k 2 . . . 0 . . . 0 k 2 . . . . . . . . . . . . . . . . . . . . . 0 0 • • • 0 k N-1 -k N             . ( 9
)
One can introduce the normalized concentration c i /(c 1 (0) + . . . + c N+1 (0)), so that c i ∈ [0, 1]. Observe also that the c N+1 evolution can be deleted and this amounts to restrict to the stoichiometric class.

For the case A B C k 1 k 2 further normalizing coordinates are:

x = ln c 1 , y = c 2 /c 1 , v = k 1 ( 10
)
and the dynamics takes the form:

ẋ = -v, ẏ = v -βyv α + vy, v = u, (11)
where

α = E 2 /E 1 , β = A 2 /A α 1 .
Equations can be extended to a more general scheme:

n 1 A n 2 B n 3 C . . . k 1 k 2 k 3
of N consecutive irreversible reactions of order n i . We obtain the following matrices:

Y = diag(n 1 , . . . , n N+1 ), A =                       0 • • • • • • • • • 0 k 1 0 • • • • • • 0 0 k 2 0 • • • 0 . . . . . . . . . . . . . . . 0 0 0 k N 0                       . Introducing "reversibility" A B C k 1 k 2 k 3 leads to ċ1 = -vc 1 + βv α c 2 ċ2 = vc 1 -βv α c 2 -β v α c 2 (12)
where [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] Mckeithan | Kinetic proofreading in T-cell receptor signal transduction[END_REF]). It is given by:

α = E 3 /E 1 , α = E 2 /E 1 , β = A 3 /A α 1 , β = A 2 /A α 1 . Case 2: McKeithan scheme (
T + M C 0 C 1 . . . C N k 1 k p,0 k p,1 k p,N k -1,0 k -1,1 k -1,N
.

The matrix Y is defined by:

Y =                                     1
ẋ = k 1 (δ 1 -x -y)(δ 2 -x -y) -(k 2 + k 3 ) x ẏ = k 2 x -k 4 y. (13)
2.4. The zero deficiency theorem [START_REF] Feinberg | On chemical kinetics of a certain class[END_REF][START_REF] Horn | General mass action kinetics[END_REF] Definition and notation. Using graph theory concepts we introduce the following: the deficiency of the network is: δ = nls, with n = number of vertices, l = number of connected components and s = dimension of the stoichiometric subspaces. The network is called weakly reversible 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.

Theorem 1. Let any simple reaction network of deficiency zero.

(i) If the network is not weakly reversible, then for arbitrary kinetics, the differential equation cannot have a positive equilibrium nor a periodic trajectory that is contained in R m

>0

(ii) If the network is weakly reversible, there exists within each strictly positive stoichiometric compatibility class precisely one equilibrium c * , the equilibrium is locally asymptotically stable with (pseudo-Helmholtz) Lyapunov function V(c, c * ) = i c i (ln c iln c * i -1) + c * i . Moreover, there is no non trivial periodic orbit.

Remark 1. Note that additional properties in the McKeithan network led to global stability property (see [START_REF] Anderson | Global asymptotic stability for a class of nonlinear chemical equations[END_REF]) and the proof of the existence of equilibrium related to the Perron-Frobenius therorem [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].

Application.

• Case 1 A B C k 1 k 2 : δ = 3 -1 -2 =
0 and the network is not weakly reversible.

• Case 2 T+M A B k 1 k 2 k 3 k 4 : δ = 3 -1 -2 = 0.
The network is weakly reversible.

Equilibrium can be easily found from [START_REF] Boscain | Optimal Syntheses for Control Systems on 2-D Manifolds[END_REF] by solving ẋ = ẏ = 0 and we get.

Lemma 2. Let (δ 1 , δ 2 ) ∈ R * + 2 fixed.
There exists a unique equilibrium point X eq = (x e , y e ) ∈ R 2 + for the system (13) given by

x e = 2δ 1 δ 2 k 1 k 4 ∆ + δ p α + β , y e = 2δ 1 δ 2 k 1 k 2 ∆ + δ p α + β where α = k 1 (k 2 + k 4 ), β = k 4 (k 2 + k 3 ), δ m = δ 1 -δ 2 , δ p = δ 1 + δ 2 and ∆ = δ 2 m α 2 + 2δ p αβ + β 2 . Proposition 1.
Take the equilibrium point X eq . The linearized system at X eq is defined by the matrix

L = -k 2 (k 2 +k 3 )-∆ k 2 +k 4 k 4 (k 2 +k 3 )-∆ k 2 +k 4 k 3 -k 4 .
We have det L = ∆ > 0 and Tr L = -∆+k 2 (k 2 +k 3 )

k 2 +k 4 + k 4 < 0.
The point X eq is globally asymptotically stable and

• if (Tr L) 2 -4 det L < 0, then X eq is a stable focus (see Fig.1 (left)),
• if (Tr L) 2 -4 det L > 0, then X eq is a stable node (see Fig. 1 (right)).

Pontryagin maximum principle and extremal curves

Statement and notation for the optimal control problem

The system is written as dc dt = f (c, T ) and controlling the temperature leads to T ∈ [T m , T M ]. In practice, thermodynamics has to be used to model the heat exchanges, in relation with the heat produce by the reactions [START_REF] Tóth | Reaction kinetics: exercises, programs and theorems[END_REF] or the heat exchange model used in the experiments. To avoid this part of the study and without loss of generality for the mathematical point of view, we shall choose v as the control variable, where v is defined by v = k i (T ) for some reaction i. This leads to deal with the single-input affine control system:

q(t) = F(q(t)) + u(t) G(q(t)) (14)
where

q = (c, v) ∈ R n , F = F(c, v) ∂ ∂c , G = ∂ ∂v , u -< u < u + and the bounds v ∈ [v m , v M ]
will not be taken into account. Note that the map : v → v is the standard Goh-transformation in optimal control and plays an important role (see [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]). We choose the labeling such that the optimization of the production is given by: max

c 1 (t f ) ( 15
)
where t f is the time duration of the batch time. An equivalent formulation will be (thanks to the maximum principle)

min t f , c 1 (t f ) = d,
where d ≥ 0 is the desired amount of the species [X 1 ] during the batch duration t f .

Maximum principle [28]

3.2.1. Geometric preliminaries Consider a single-input control system: q(t) = f (q(t), u(t)), q ∈ R n , f real analytic and the set of input is the set U of bounded measurable mappings u : [0, t f (u)] → [-1, +1]. Fixing q(0) = q 0 we denote by q(•, q 0 , u) the solution associated to u(•) and starting from q(0). The accessible set in time t f is the set: A(q 0 , t f ) = ∪ u(•) q(t f , q 0 , u), image at time of the extremity mapping: E q 0 ,t f : U u(•) → q(t f , q 0 , u) ∈ R n . The set of inputs is endowed with the L ∞ -norm topology. The Maximum Principle is a parametrization of the boundary of accessibility set.

3.2.2. Needle or Weierstrass variations. First order Pontryagin cone K 1 (t f ) Take u(•) ∈ U such that the solution written shortly q(•) is defined on [0, t f ] and let t be a Lebesgue time on [0, t f ]. A needle variation of u(•) with data (t, v, ε) is the L 1 perturbation defined by: u ε = v on [t -ε, t] and equal to the control u(•) elsewhere with ε > 0 and v constant ∈ [-1, +1]. We shall denote in short by q ε (•) the response to u ε (starting at time 0 from q 0 ). Consider at final time t f , the curve ε → q ε (t f ). The first order Pontryagin cone K 1 (t f ) is the smallest convex closed cone containing tangent vectors d dε |ε=0 q ε (t f ) for all perturbations. The key point of the Maximum Principle is to use the Pontryagin cone as an approximation of the accessibility set, see [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF].

Statement of the maximum principle

One needs the following. The pseudo-Hamiltonian (or unmaximized Hamiltonian) is H(q, p, u) = p • f (q, u), where p ∈ R n \{0} is the adjoint vector. We denote M(q, p) = max |u|≤1 H(q, p, u). Consider first the time minimal control problem, with q(0) = q 0 fixed and q(t f ) ∈ N: terminal analytic manifold. One has.

Theorem 2. Assume (q * (•), u * (•)) is an optimal time minimal solution on [0, t * f ], then there exists p * (•) nonzero such that a.e. on [0, t * f ] q * (t) = ∂H ∂p ((q * (t), p * (t), u * (t)), (16a) ṗ *
(t) = - ∂H ∂q ((q * (t), p * (t), u * (t)) (16b)
H((q * (t), p * (t), u * (t)) = M((q * (t), p * (t)). ( 16c
)
Moreover t → M((q * (t), p * (t)) is constant and non negative and at the final time one has the transversality condition:

p * (t f ) ⊥ T q * (t f ) N. ( 17
)
Remark 2. If we replace the time minimal problem by a Mayer problem:

min u(•) ϕ(q(t f )), (18)
the same conditions (16a), (16b), (16c) hold, while (17) is replaced by:

p * (t f ) = - ∂ϕ ∂q (q * (t f )). ( 19
)
Definition 1. An extremal triplet (q, p, u) is a solution of (16a), (16b), (16c) and it is called a BC-extremal if the transversality condition is satisfied. An extremal control is called regular if |u(t)| = 1 a.e. and singular if ∂H ∂u (q, p, u) = 0 everywhere. It is called exceptional if M = 0. A regular extremal is called bang-bang if the number of switches is finite.

Computation of singular extremals and properties

Notations. Let X, Y be two real analytic vector fields of R n , the Lie bracket is defined by: [X, Y](q) = ∂X ∂q (q)Y(q) -∂Y ∂q (q)X(q). The Hamiltonian lift of X is

H X (z) = p • X(q), z = (q, p). The Poisson bracket is defined by {H X , H Y }(z) = p • [X, Y](q).
Computation of singular extremals.

Definition 2. Let H(z, u) = p • f (q, u). The condition ∂ 2 H ∂u 2 ≤ 0 (resp < 0) is called (resp. strict) Legendre-Clebsch condition. Assume f (q, u) = F(q) + u G(q). The condition ∂ ∂u d dt 2 ∂H ∂u = {{H G , H F }, H G } ≥ 0 (resp. > 0) is called the (resp. strict) generalized Legendre-Clebsch condi- tion.
Computation. If the strict Legendre-Clebsch condition is satisfied one uses ∂H ∂u = 0 to compute the singular control û(z) using the implicit function theorem. Plugging such û(z) in H(z, u) leads to the true Hamiltonian denoted Ĥ(z). In the affine case, one has ∂ 2 H ∂u 2 = 0 and we proceed as follows. The relation ∂H ∂u = 0 leads to H G (z) = p • G(q) = 0. Deriving twice with respect to time along the extremal solution leads to:

H G (z) = {H G , H F }(z) = 0, {{H G , H F }, H F }(z) + u {{H G , H F }, H G }(z) = 0. ( 20
)
Assume {{H G , H F }, H G }(z) 0, the singular extremal is called of order 2 and the singular control u is computed as û(z), using relation [START_REF] Horn | General mass action kinetics[END_REF]. Plugging such û(z) into H(z, u) leads to a true singular Hamiltonian denoted H s (z), one has: Lemma 3. Singular extremals of order 2 are the solutions of:

dq dt = ∂H s ∂p (z), dp dt = - ∂H s ∂q (z) (21)
with the constraints H G (z) = {H G , H F }(z) = 0.
Moreover in order to be admissible the singular control given by

u s (z) = - {{H G , H F }, H F }(z) {{H G , H F }, H G }(z) (22)
has to satisfy the admissibility constraints |u s (z)| ≤ 1.

Definition 3. Let (z, u) be a singular extremal of order 2 and M = H F = 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 {{H G , H F }, H G } > 0 (resp. < 0).
Time optimality status of singular extremals for the point to point problem. (N: single point).

One needs the following seminal results based on [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF] and [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], see also [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] for a tutorial presentation.

Proposition 2 (High-order maximum principle). The generalized Legendre-Clebsch condition is a necessary small time optimality condition.

Further results need the following.

Statement of the limit problem. We consider the time minimal problem with fixed end-point for the affine system:

dq dt = F(q) + u G(q). ( 23
)
The limit problem is the following:

• the constraints |u| ≤ 1 on the control are relaxed and we admit a specific class of impulse controls. More precisely, a trajectory is a finite concatenations of arcs corresponding to bounded measurable controls and finite jumps in the G direction. The set of admissible controls extending U is denoted U .

Assumptions. Consider the affine system q = F(q) + u G(q) (u ∈ R) and let (q, p, u) be an extremal of order 2 on [0, t f ]. Assume the following on U ⊂ R n .

(H 0 ) F, G are linearly independent and the reference extremal curve t → q(t) is one-to-one. Since the concept of singular-extremal is feedback invariant (see [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF]) one may set u ≡ 0. Let us denote by

K 1 (t) = span{ad k F • G(q(t)); k ∈ N} with ad F • G = [F, G].
Then it is known that K 1 (t) is the first-order Pontryagin cone along q(•) and the codimension is non zero.

Moreover let us then assume that:

(H 1 ) ∀t ∈ [0, t f ], K 1 (t)
is of codimension one and is spanned by the vectors:

ad k F •G(q(t)), k = 0, . . . , n -2. (H 2 ) If n ≥ 3, ∀t ∈ [0, t f ], F(q(t)) span{ad k F • G(q(t)); k = 0, . . . , n -3}. Remark 3. Note ∀t ∈ [0, t f ], [[G, F], G](q(t)) K 1 (t).
Theorem 3. Let (q(•), p(•), u(•)) be a singular extremal defined on [0, t f ] and satisfying assumptions (H 1 ), (H 2 ). Note that the adjoint vector p is then unique up to a non-zero factor, and for any time t in [0, t f ], p(t) is orthogonal to K 1 (t); so ((q(t), p(t)) is of order 2. Then there exists a C 0 -neighborhood U of the reference trajectory q(•) such that q(•) is a time-minimizing (resp. maximizing) trajectory with respect to all solutions of (23) contained in U and joining q(0) to q(t f ) (the set of admissible control being U ) if (q(•), p(•), u(•)) is exceptional or hyperbolic (resp. elliptic) and if t f < t 1c where t 1c is the first conjugate time to t f along q(•).

Proof. See [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].

Algorithm to compute conjugate times. They are described in [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] and practically implemented in the software HamPath [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF].

Relation between affine and non affine case using the Goh transformation for chemical networks.

Recall that for our network:

ċ = f (c, v), v = k i is extended into q = F(q) + u G(q), F = f (c, v) ∂ ∂c , G = ∂ ∂v with q = (c, v). Denote H = p c • f (c, v
) and H = p • (F(q) + u G(q)), p = (p c , p v ). One has the following relation between the corresponding singular extremals. Lemma 4. The pair (q, p) is solution of:

q = ∂H ∂p , ṗ = - ∂H ∂q , ∂H ∂u = 0 if and only if p v = 0 and (c, p c , v) is solution of: ċ = ∂ H ∂p c , ṗc = ∂ H ∂c , ∂ H ∂v = 0
and moreover the following relations are satisfied

d dt ∂H ∂u = {H G , H F } = - ∂ H ∂v (24)
∂ ∂u

d 2 dt 2 ∂H ∂u = {{H G , H F }, H G } = - ∂ 2 H ∂v 2 . ( 25
)
In particular [START_REF] Ledzewicz | Singular controls and chattering arcs in optimal control problems arising in biomedicine[END_REF] relates the Legendre-Clebsch and the generalized Legendre-Clebsch condition of both systems. 10

Computation in small dimensions.

• n = 2. Note q = (x, y). Introduce the following determinants

D = det(G, [[G, F], G]), D = det(G, [[G, F], F]) D = det(G, F), D s = det(G, [G, F]).
Eliminating p leads to:

-S: Singular arcs located in D s = 0.

-Singular control: u s (q) = -D (q)/D(q).

-Hyperbolic: DD > 0.

-Elliptic: DD < 0.

• n = 3. Note q = (x, y, z). 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).
Eliminating p leads to:

-Singular control: u s (q) = -D (q)/D(q).

-Singular vector fields: q = F(q) -D (q) D(q) G(q). -Hyperbolic: DD > 0.

-Elliptic: DD < 0.

-Exceptional: D = 0.

Small time classification of extremals. Construction of semi-normal forms

In this section, we present the basic results and techniques from singularity theory initiated in [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 small time extremals curves, which will be the basic tool to analyze the optimal control problem. It is the construction of semi-normal form using the action of the pseudo-group of diffeomorphisms or symplectomorphisms, combined with specific feedback transformations.

Classification of regular extremals. Notations: we denote by σ + (resp. σ -) a bang arc with constant control u = +1 (resp. -1) and σ s is the singular arc. We denote σ 1 σ 2 an arc σ 1 followed by σ 2 . The surface Σ : p • G(q) = 0 is called switching surface and let Σ ⊂ Σ given by: p

• G(q) = p • [G, F](q) = 0. If z(t) = (q(t), p(t))
is an extremal curve on [0, t f ], we note Φ(t) = p(t) • G(q(t)) the switching function (which codes the switching times). Differentiating twice with respect to time one gets:

Φ(t) = p(t) • [G, F](q(t)) , (26)
Φ(t) = p(t) • ([[G, F], F](q(t)) + u(t) [[G, F], G](q(t))). (27)
From this calculus we deduce.

Ordinary switching time. A time t ∈]0, t f [ is called an ordinary switching time and z(t) ∈ Σ an ordinary switching point if Φ(t) = 0 and Φ(t) 0. Clearly we have.

Lemma 5. In the ordinary case, near z(t) every extremal solution projects onto σ + σ -if Φ < 0 and σ -σ + if Φ > 0.

Fold case. If Φ(t) = Φ(t) = 0, then z(t) ∈ Σ . The situation is more intricate and is technically complicated to analyze, see [START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF] for the details. Let u = ε, ε = ±1 and Φε (t

) = p(t) • ([[G, F], F](q(t)) + ε[[G, F], G](q(t))
). If Φε (t) 0 for ε = ±1, the point z(t) is called a fold point and we have three cases assuming Σ is a regular surface of codimension two.

Case 1: parabolic case:

Φ+ (t) Φ-(t) > 0.
Case 2: hyperbolic case: Φ+ (t) > 0 and Φ-(t) < 0.

Case 3: elliptic case: Φ+ (t) < 0 and Φ-(t) > 0.

Denote by u s (t) the singular control given by

p(t) • ([[G, F], F](q(t)) + u s (t)[[G, F], G](q(t))) = 0.
From the above classification, we deduce. In the hyperbolic or elliptic case, through

z(t) if p(t) • [[G, F], G](q(t))
0 there exists a singular extremal which is strictly admissible that is |u s (t)| < 1. Moreover we are in the hyperbolic case (resp. elliptic) if this quantity is > 0 (resp. < 0). In the parabolic case, it can be absent (for instance in the linear case:

F(q) = A q, G(q) = b) or not admissible that is |u s (t)| > 1.
One has Theorem 4. In the neighborhood of z(t) every extremal projects onto:

• In the parabolic case: σ + σ -σ + or σ -σ + σ -.

• In the hyperbolic case: σ ± σ s σ ± .

• In the elliptic case: every extremal is bang-bang, i.e. of the form σ + σ -σ + σ -. . . but the number of switches is not uniformly bounded.

Note that the elliptic case opens the road to the so-called Fuller phenomenon [START_REF] Kupka | The ubiquity of Fuller's phenomenon[END_REF][START_REF] Fuller | Study of an optimum non linear system[END_REF] where a connection with a singular arc can be realized but with an infinite number of switches.

The general tool to analyze this situation is the concept of semi-normal form and we shall make a simplified presentation using the planar case. Also this will lead to extensions which will be crucial in our next section of classifications of extremals near a manifold of codimension one.

Introducing semi-normal forms. We consider a 2D-system: q = F(q) + u G(q), q = (x, y). We note C : det(F, G) = 0 the collinearity set and we restrict our system to U : R 2 \ C. Using the previous computations, singular trajectories are located on S : det(G, [G, F]) = 0. The direction of the adjoint vector p is obtained with p • G = 0 and using the previous notation the singular control is u s = -D /D. One picks a reference singular curve σ s (t). In a C 0 -neighborhood of this curve we can choose coordinates q = (x, y) such that: σ : t → (t, 0) and G = ∂ ∂y . This leads to the following semi-normal form: where the terms . . . are not relevant in our study. Note if we relax the bound |u| ≤ 1, we get the limit problem so that u s (x) can be set to 0 using a feedback transformation.

ẋ = 1 + εy 2 + . . . , ẏ = (u -u s (x)) + . . . , u = +1 u = -1 singul ar Σ p • G(q) < 0 p • G(q) > 0
Clearly the optimality status of the reference singular arc can be immediately deduced from this form. If ε < 0, we are in the hyperbolic case and the singular direction is time minimizing in a C 0 -neighborhood of the reference singular extremal and time minimizing in the elliptic case. The adjoint vector is oriented with the convention p • F > 0 and we deduce the standard turnpike phenomenon characteristic of the hyperbolic situation, that is, provided |u s | ≤ 1, the optimal policy is of the form σ ± σ s σ ± . (Note it is valid not only for small time). It is represented on Fig. 2 in the (q, p) space. Of course, this situation doesn't cover the various situations encountered in the applications. One important case being the saturation phenomenon where |u s | = 1, which will be discussed later.

Next we shall introduce the bridge phenomenon.

Connecting singular arcs: the bridge phenomenon. The aim of this paragraph is to analyze the connection between singular extremals, extending the previous classification. For the sake of the tutorial aspect we shall use a planar semi-normal form.

Definition 4. A bridge is a bang arc σ b connecting two singular arcs σ 1 s , σ 2 s so that σ 1 s σ b σ 2 s is an admissible extremal curve.

Birth of the model. We consider the (limit) model: ẋ = u, ẏ = 1x 2 y so that σ : (0, t) is a singular arc which is hyperbolic if y > 0 and elliptic if y < 0. Using the notation q = F + u G, singular trajectories are located on S : det(G, [G, F])(q) = 0 = {xy = 0}. Hence we have a ramification at (0, 0) between two singular lines: the vertical axis which can be followed with the control u s = 0 and the horizontal axis which can be followed by the "u = ∞" control.

To get our model we must bend this axis to make it tractable with a singular control u s → ∞ when going to 0. This leads to the following.

The (symmetric) bridge model. One takes

F = (1 -x 2 y) ∂ ∂y , G = -(y -1) ∂ ∂x + x ∂ ∂y . ( 28
)
The solutions of G are circles contained on (0, 1), note a parameter λ can be introduced using 1 → λ. Also one takes the bound |u| ≤ 1, but an additional parameter is |u| ≤ M. Lie brackets computations give us. Lemma 6. The singular set S : {q; det(G, [G, F])(q) = 0} is the union of the vertical axis (Oy) and the algebraic curve defined by: l(x, y) = x 2 (1 -2y) + 2y 3 -4y 2 + 2y + 1 = 0.

The equation l(x, y) = 0 is equivalent to x = ± √ 2y 3 -4y 2 +2y+1 2y-1

. The intersection with the vertical axis x = 0 is (O, y 0 ), y 0 -0.297157.

We represent on Fig. 3 the two singular trajectories: the set l = 0 denoted L and the vertical axis (Oy) denoted L . Further computations lead to. Proposition 3. The vertical singular line L is hyperbolic for y > y 0 and elliptic for y < y 0 . The singular line L is hyperbolic and the singular control along L ∩ {x ≤ 0} is given by:

D•D < 0 D•D > 0 y 0 L L y x d e t ( G , [ G , F ] ) = 0
u s (q) = - D (q) D(q) = A B
where A = -8y 7 + 30y 6 -44y 5 + 38y 4 -27y 3 + 15y 2 -6y + 2 and B = (2y -1)(2y 3 -4y 2 + 2y + 1) 8y 3 -15y 2 + 9y -3 .

We represent on Fig. 4 this singular control showing in particular the existence of two saturating points |u| = 1.

To conclude, this leads to construct a bridge using a simple numeric simulation, see Fig. 5, and an extremal solution which projects onto σ -σ L σ b + σ L , see Fig. 6. Note that time minimality can be analyzed using standard discussion, moreover mathematical estimates can be obtained near (O, y 0 ), making λ → +∞ and M → +∞, see Fig. 7 for the effect of λ → +∞ so that y 0 → 0 and L tends to the horizontal axis.

Remark 4 (Geometric remark). Interaction between the two singular curves is coded by the singularities of the flow q = F(q) -D (q) D(q) G(q) since D + u s D = 0. It can be analyzed using a time reparametrization using the dynamics: q = D(q)F(q) -D (q)G(q), with singular equilibrium in D = D = 0. This geometric remark is crucial for the extension to the non planar case.

Time minimal synthesis near the terminal manifold

The basic technical study applicable to chemical batch reactions was developed in the series of papers [START_REF] Bonnard | Towards a geometric theory in the time minimal control of chemical batch reactors[END_REF][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]. The problem is the following. Consider the system: q = F(q)+u G(q), |u| ≤ 1 with terminal manifold N = f -1 (0) where f : R n → R is a submersion. One consider the following local problem: take q 0 ∈ N, compute, in a small neighborhood U of q 0 (the size of U is not a priori known), the optimal closed control u * (q) to steer (staying in U) from q ∈ N to N ∩ U, in minimum time. Such a problem is called the local time minimal synthesis, with terminal manifold. Note that a lot of central work in this area was done in the 80 s, see [START_REF] Sussmann | Synthesis, presynthesis, sufficient conditions for optimality and subanalytic sets[END_REF][START_REF] Ledzewicz | Singular controls and chattering arcs in optimal control problems arising in biomedicine[END_REF][START_REF] Boscain | Optimal Syntheses for Control Systems on 2-D Manifolds[END_REF] when the terminal manifold is a point, to classify in small dimension the optimal synthesis in relation with the Lie algebraic structure of (F, G) at the given point q 0 . If all those contributions are valuable in our study, we shall make the computation explicit in the case where the terminal manifold N : f = 0 is of codimension one and using also different tools, that is the construction of semi-normal forms for the action of the pseudo-group G formed by local diffeomorphism and feedback transformation u → -u (so that σ + and σ -can be exchanged). This technique is coming from singularity theory, see [START_REF] Martinet | Singularities of smooth functions and maps[END_REF] for a useful introduction of this area. Note we shall work in the C k -category, where k ≥ 1 will not be precised and the semi-normal form is related to a semi-algebraic stratification on the jets spaces of (F, G, N). One needs to introduce some definitions and concepts.

Notations and definitions

Note that for chemical systems N : c 1 = d and the chemical species are such that c 1 (0) < d. Define N ⊥ = {(q, p); p • v = 0, ∀v ∈ T q N} and let n be the normal to N oriented as the outward normal to N (using the chemical model). Let z = (q, p) be a BC-extremal, t ∈ [t f , 0], t f < 0 (since we integrate backwards from the terminal manifold). A time s is called switching time if s belongs to the closure of the switching time on [t f , 0] (in the general setting, switches can accumulate due to the Fuller phenomenon [START_REF] Kupka | The ubiquity of Fuller's phenomenon[END_REF]). We shall denote by K the switching points and by W the set of switching points for minimizing curves. If W is stratified, a stratum is of the first kind if the optimal trajectories are tangent and of the 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 trajectory loses its optimality. The reader can refer to [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] and [START_REF] Sussmann | Synthesis, presynthesis, sufficient conditions for optimality and subanalytic sets[END_REF] for the details and results about the concepts in the frame of semi-analytic geometry.

Next one needs concepts from singularity theory.

Concepts and definitions. F, G and N are in the C w -category, N = f -1 (0), f is a submersion. The set of triplet (F, G, f ) is endowed with C ∞ -Whitney topology. We denote by j k F(q 0 ) (resp. j k G(q 0 ), j k f (q 0 )) the k-jet of F (resp. G, f ) that is the Taylor expansion at order k. We say that the system (G, F, f ) has at q 0 a singularity of codimension i if j k F(q 0 ), j k G(q 0 ), j k f (q 0 ) ∈ Σ i , where Σ i is a semi-algebraic submanifold of codimension i in the jet space. Taking a point q 0 with a singularity of codimension i, an unfolding is a C 0 -change of coordinates φ near q 0 such that (small) time minimal synthesis is described by a system ẋ =

F( x, λ) + u G( x, λ) , |u| ≤ 1, x ∈ R n-m
and λ is a parameter.

Local syntheses

We shall present the main step to compute the time minimal synthesis and thanks to the concept of unfoldings we shall restrict our study to the 3-dimensional case, which is also the situation corresponding to our test bed cases. The system is written q = F(q) + u G(q), |u| ≤ 1 and let q = (x, y, z) be the coordinates. The terminal manifold is N and we suppose that the problem is flat that is G being identified to ∂ ∂z and G is tangent everywhere to N. If n is normal to N, outwardly oriented, we introduce the following to stratify the terminal manifold.

• 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, that is N ⊥ ⊂ Σ : switching surface.

Generic case

Take q 0 ∈ N, q 0 S, q 0 E, then (n, q 0 ) is an ordinary switching point and according to section 3.2.5, near q 0 every BC-extremal is of the form σ

+ if n • [G, F](q 0 ) < 0 and of the form σ -if n • [G, F](q 0 ) > 0.
This gives the local synthesis.

Generic hyperbolic singular case

One can take q 0 = 0 and we make the following assumptions at 0.

• The tangent space to N at 0 is G(0), [G, F](0).

• The set of points L 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 single BC-extremal denoted σ transverse to N. Near 0 one can make the following normalization (that is using adapted coordinates): G being identified to ∂ ∂z , L to the axis (Oy) and σ can be identified to: t → (t, 0, 0) with image being the (Ox)-axis. Note that the target is identified to x = 0. The system is written as:

ẋ = 1 + i+ j>0 a i j (x)y i z j ẏ = i+ j>0 b i j (x)y i z j ż = (u -û(x)) + i+ j>0 c i j (x)y i z j
where û(x) is the singular control associated to σ. One further normalization leads to take

[G, F] | σ = [G, F](0).
Here we obtain a semi-normal form:

ẋ = 1 + a(x)z 2 + 2b(x)yz + c(x)y 2 + ε 1 ẏ = d(x)y + e(0)z + ε 2 ż = (u -û(x)) + f (x)y + g(0)z + ε 3 with b(0) = 0, e(0) 0 and ε 1 (resp. ε 2 , ε 3 ) are terms of order ≥ 3 (resp. ≥ 2) in y, z.
Computing, one has n • [[F, G], G](0) = 2a(0). We assume:

• a(0) < 0 (hence the singular arc is hyperbolic) and moreover û(0) ∈] -1, +1[ that is the singular arc is strictly admissible.

Hence (q 0 , n) is an hyperbolic fold point and using 3.2.5 an extremal near such a point is σ ± σ s σ ± . Moreover arcs σ + and σ -are defined on [t f , 0] such that σ + (0), σ -(0) near q 0 are σ + in the domain z < 0 and σ -in the domain z > 0 and are not intersecting since ż ∼ (uû(0)), hence: z(t) = (εû(0)) t + s, s = z(0) and ε = ±1. Using the normalization, one can evaluate the switches of BC-arc σ + and σ -for t < 0. One has:

p 1 (t) = 1 + o(1), p 2 (t) = o(1), ṗ3 (t) ∼ -2a(0)p 1 + (ε -û(0)) z + . . . . (29)
Hence, we have p 3 (t) ∼ -2a(0) t (s + (εû(0)) t/2 + o(s, t)). Using εs < 0, one gets that p 3 (t) is not vanishing for t < 0. Hence we deduce the following which gives the C 0 -unfolding.

Proposition 4. Under arc assumptions near 0, the optimal synthesis is described by the model:

ẋ = 1 + a(0)z 2 , ẏ = 0, ż = (u -û(0)
) and is represented in each invariant plane y = constant by Fig. 8.

Assume now that û(0) ] -1, +1[ then P = (q 0 , n) is a parabolic fold point and the synthesis can be again deduced from our normalizations. Each extremal is bang-bang with at most two switchings, and one is located at the terminal point on N.

Using the previous estimates one can compute the optimal switching locus W -near N and the optimal synthesis is represented on Fig. 9 in an invariant plane y = constant.

σ - W - σ + σ 0 + σ + 0 N ∩ y =constant

Saturate case.

One shall briefly present the transition between the hyperbolic and parabolic case due to the saturation of the control |û(0)| = 1, see [START_REF] Bonnard | Classification générique de synthèses temps minimales avec cible de codimension un et applications[END_REF] for a complete study. From our previous analysis the model is: ẋ = 1 + a(x)z 2 + 2b(x)yz + c(x)y 2 ẏ = d(x)y + c(0)z ż = (u -1)ûx xûy y with a(0) < 0, b 0 and ûx = ∂û ∂x (0), ûy = ∂û ∂y (0) and ûx , ûy are generically non zero and we may assume ûy > 0. Computations lead to the two cases represented in Figs. 10-11. Note that the existence of two syntheses is discriminated by vanishing (case 2) and the vanishing (case 1) of the singular curve.

- σ- σ- σ- σ- σ- σ- σ- σ- N σ+ W+ σ- σ- σ- σ+ W+ N y = constant < 0, v 0 < v sat y = 0, v 0 = v sat y = constant > 0, v 0 > v sat
E 0 σ - N E 0 σ - N σ - Case 1: v 0 < v - Case 2: v 0 > v - Figure 12: Local synthesis near E 0 = E ∩ (v = v 0 ) in the leaf v = v 0 when α > 1.

The exceptional case

The exceptional case corresponds either to bad controllability properties of the system, related to singular arcs or bad accessibility property of the target. We made a concise presentation related to our study for chemical systems, see [START_REF] Launay | The generic local structure of time-optimal synthesis with a target of codimension one in dimension greater than two[END_REF] for the details.

The generic case. It is a situation with (minimal) contact one of bang arcs with the target manifold. In this case, restricting to a planar case thanks to the concept of unfoldings we have two cases represented on Fig. 12 obtained by obvious topological analysis. We note E 0 the point where arcs σ + , σ -are tangent. In the case 2, the black domain represents the set of the domain with no optimal trajectory.

The codimension one case. The situation is more intricate, since we cannot use a planar invariant foliation. We have several cases corresponding to two situations. Situation 1. Bang arc with a contact of order two with the target at a point of E. We have two cases represented on Fig. 1314.

Situation 2. It concerns E ∩ S which is a singular arc tangent to N and according to our classification it concerns a singular exceptional arc. This complicated case is discussed 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].

c 1 v E σ - σ - σ - σ - σ - σ - E - A 1 v - O
σ + σ + σ + O c 1 v A 1 E + E

Application to the two test bed cases

The case

A k 1 -→ B k 2 -→ C
The system q = F(q) + u G(q), using the coordinates q = (x, c 2 , v), x = ln c 1 takes the form ẋ = -v ċ2 = ve xβv α c 2 v = u [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 in the coordinates q = (x, y, v) y = c 2 /c 1 one gets

ẋ = -v ẏ = v -βv α y + vy v = u. (31)
It is integrated on [t f , 0] together with the adjoint equation:

ṗ = -p ∂F ∂q (q) + u ∂G ∂q (q) .
Without losing any generality one may assume |u| ≤ 1.

Lie brackets computations

One uses the coordinates q = (x, y, v) and we get: Proof. One has: ẏ = v 1 -βv α-1 y + y and assume there exists τ , ẏ(t) > 0 for t < τ , ẏ(τ ) = 0 and ẏ(τ) < 0 for [τ , τ]. Differentiating and evaluating at τ one has

F(q) = -v ∂ ∂x + (v -βv α y + vy) ∂ ∂y , G(q) = ∂ ∂v , [G, F](q) = ∂ ∂x + (-1 + αβv α-1 y -y) ∂ ∂y , [[G, F] , F] (q) = (α -1)βv α ∂ ∂y , [[G, F] , G] (q) = α(α -1)βv α-
ÿ(τ ) = v -(α -1)βyv α-2 v + 1 with v = -1 along σ -. Hence since α > 1, ÿ(τ ) > 0. A contradiction.
Lemma 12. A sequence σ s σ -with non empty subarc cannot exist for a BC-extremal if α > 1.

Proof. Assume such a sequence exists with σ s on [t 1 , t 2 ] and σ -on [t 2 , t 3 ]. One has the relations:

Φ(t 2 ) = Φ(t 2 ) = 0, Φ(t 3 ) = 0.
Since Φ(t 2 ) = Φ(t 3 ) = 0, from Rolle's theorem, there exists t 2 < t 2 < t 3 such that Φ(t 2 ) = 0. Using Rolle's theorem again with Φ(t 2 ) = Φ(t 2 ), there exists t 2 < t 2 < t 2 such that Φ(t 2 ) = 0. Note that from Lemma 7 and Lemma 10 one has along σ -:

Φ(t) = p 2 (t)(α -1)βv α-2 αy (-1 -u s )
with u s = -v 2 /(αy).

Hence if Φ(t 2 ) = 0 one has u s (t 2 ) = -1 and σ -meets the saturating set S : u s = -1.

At the junction time t 2 between σ s and σ -, one has u s (t 2 ) ≥ -1, ċ2 (t 2 ) > 0 from Lemma 9 and along σ

-: ẏ = ċ2 /c 1 + c 2 /c 1 v ≥ 0 and v < 0. Moreover along σ -, sign Φ(t) = sign (-1 -u s ), since ẏ ≥ 0, v ≤ 0, -1 -u s is a decreasing function of time.
Since Φ(t 2 ) ≤ 0 it turns out that for t ∈]t 2 , t 3 [ one has Φ(t 2 ) < 0. This contradicts Φ(t 2 ) > 0 for t 2 > t 2 .

From Corollary 1 and Lemma 12 we deduce.

Theorem 5 (α > 1). Every optimal trajectory has at most two swtichings and is of the form σ + σ -σ s where each arc of the sequence can be empty.

5.1.4. Local optimal synthesis near N, the case α < 1 It will be easily follows from the classification of Section 4.1. We represent on Fig. 15 the stratification of the terminal manifold by S : n • [G, F](q) = 0 and E : n • F(q) = 0. The optimal policy near N is σ + . Theorem 6. In the case α < 1, the local optimal policy is σ + . 5.1.5. Local optimal synthesis near N, the case α > 1

E S u = +1 E + A 1 O v c 1
Again we use the classification of Section 4.1 but the situation is more intricate. It is represented on Fig. 16 and we refer to Section 4.1 for the corresponding figures.

• Singular hyperbolic arcs reaching the manifold are on S and we have a saturating point where u s = -1.

• The boundary of the accessibility set is associated to E, with E -accessible.

The global optimal policy follows from Theorem 5: an optimal trajectory being of the form σ + σ -σ s where each arc of the sequence can be empty. It can be numerically obtained, integrating backwards on [t f , 0] using the local synthesis resolution of our analysis.

• Near a point E 0 E -of E, the optimal synthesis is described by a C 0 invariant foliation F : (v = v 0 ) the leaves of which are given by Fig. 12; near E -there is no such foliation and the synthesis is given by Fig. 13.

• Near the point S sat of S, the optimal synthesis is given by Fig. 11.

The McKeithan scheme and the computational complexity

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 = k 1 and reduced to the stoichiometric class 0

c 1 A 1 v sat E - S sat S E v u = -1 parabolic hyperbolic u = +1 v -

T + M

+ A = δ 1 and T + M + B = δ 2 takes the form ẋ = -β 2 xv α 2 -β 3 xv α 3 -δ 3 v (x + y) + δ 4 v + v (x + y) 2 ẏ = β 2 xv α 2 -β 4 yv α 4 with 0 ≤ x ≤ δ 1 , 0 ≤ y ≤ δ 2 , δ 3 = δ 1 + δ 2 , δ 4 = δ 1 δ 2 , k 2 = β 2 v α 2 , k 3 = β 3 v α 3 , k 4 = β 4 v α 4 .

Lie brackets computations

One gets:

• F = -β 2 xv α 2 -β 3 xv α 3 -δ 3 v(x + y) + δ 4 v + v (x + y) 2 ∂ ∂x + β 2 xv α 2 -β 4 yv α 4 ∂ ∂y , • G = ∂ ∂v , • [G, F] = x(α 2 β 2 v α 2 -1 + α 3 β 3 v α 3 -1 + δ 3 ) + δ 3 y -δ 4 -x 2 -2xy -y 2 ∂ ∂x + α 4 β 4 yv α 4 -1 - α 2 β 2 xv α 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: and the singular control is given by: u s = -D (q)/D(q). Then we have ∆ = (δ 1 -δ 2 ) 2 + 4d (α 2 β 2 v α 2 -1 + α 3 β 3 v α 3 -1 ) > 0 so that there is no ramification and S contains at most two real positive branches. Definition 5. A semi-bridge occurs at a point q ∈ S if n • [[G, F] , G] (q) = 0.

• D(q) = ((α 4 -1)α 4 β 4 yv α 4 -3 -(α 2 -1)α 2 β 2 xv α 2 -3 )(α 2 β 2 xv α 2 + α 3 β 3 xv α 3 + δ 3 v(x + y) -δ 4 v - vx 2 -2vxy -vy 2 ) + x(α 2 2 β 2 v α 2 -α 2 β 2 v α 2 + (α 3 -1)α 3 β 3 v α 3 )(α 2 β 2 xv α 2 -3 -α 4 β 4 yv α 4 -3 ), • D (q) = β 2 v α 2 -2 (α 2 β 2 xv α 2 + α
Computing, a semi-bridge occurs if v α 3 -α 2 = -(α 2 -1)α 2 β 2 (α 3 -1)α 3 β 3 .

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

E : -β 2 dv α 2 -β 3 dv α 3 + d 2 v -dδ 3 v + y (2dv -δ 3 v) + δ 4 v + vy 2 = 0. ( 36
)
The discriminant of the polynomial n • X(q) in y is ∆ = v 4d(β 2 v α 2 + β 3 v α 3 ) + v(δ 1 -δ 2 ) 2 > 0 and E contains at most two real positive branches. Fig. 17 gives a picture of stratification of x = d for the McKeithan system with a point where S is folded.

Outline of the normalization. We choose the coordinates q = (x, y, z), the point q 0 is (0, 0, 0) so that the system q = F(q) + u G(q), |u| ≤ 1 is such that N : {x = 0}, G = ∂ ∂z and we have the following:

• the set S ∩ N : n • [G, F] (q) = 0, x = 0 is y = z 2 ,
• the point q 0 satisfies n • [[G, F] , G] (q 0 ) = 0 and n • [[G, F], F](q 0 ) 0,

• the trajectory associated to u = 0 defined on [t f , 0] with γ(0) = 0 is γ : t → (t, 0, 0). A brief recap of the dynamics is presented at constant temperature using the Feinberg-Horn-Jackson graph and is shown to be simple using the zero deficiency assumption.

The optimal control problem is analyzed in the frame of singularity theory and a first crucial step is to make the calculation of the closed loop optimal control for a time minimal control problem for which the terminal manifold is of codimension one and for a single input system. A classification of the local optimal syntheses is presented based on generic assumptions on Lie brackets relations. It is obtained from on explicit computations using adapted coordinates to simplify the calculations which amount to compute semi-normal form.

This approach is shown to be relevant to applications, since the classification can be used for different chemical networks. It is tested on two important schemes for applications: a sequence of two irreversible reactions and the McKeithan scheme. The classification presented in this article allows to completely solve the first case and has to be completed to make a complete analysis of the McKeithan scheme. Such completion is a good exercise for the reader to check if the tutorial objective of this article is fulfilled.

Additionally, for the first case a complete global solution of the problem is obtained based on a description of the (global) switching policy. For the McKeithan network, there is a serious difficulty to fulfill this scope which is revealed by the complexity of the singular flow, which contains in particular many singularities. They have to be analyzed in relation with the bridge phenomenon and the existence of equilibria. This is an important and challenging mathematical problem.

Note also that the graph structure of the network is not exploited directly in the optimal problem. It can be an important issue to handle the computational analysis of the Lie algebra structure associated to the system, which is crucial to understand the optimal control problem. This paves the road for further studies to deal in particular with nonzero deficiency network.

Figure 1 :

 1 Figure 1: Phase portrait for the McKeithan model. (left) Focus; (right) Node.

Figure 2 :

 2 Figure 2: Fold case in the hyperbolic case in dimension 2 and the turnpike phenomenon

Figure 3 :

 3 Figure 3: (Oy) and L are two singular trajectories for the system (28) represented by two continuous lines. The filled regions are the points where D • D > 0.

Figure 4 :

 4 Figure 4: Singular control u s for x < 0. There exists two saturating points.

Figure 5 :

 5 Figure 5: Time evolution of an extremal for the problem (28) of reaching the origin in minimal time starting from (x(0), y(0)) = (-1, -1).

Figure 6 :LD • D < 0 D • D < 0 D • D > 0 Figure 7 :

 60007 Figure 6: Trajectory associated to the extremal of Fig.5. σ b+ is a bridge connecting two switching points of the singular set and the singular control saturates at the point S sat .

Figure 8 :

 8 Figure 8: Hyperbolic case. |û(0)| < 1.

Figure 9 :

 9 Figure 9: Parabolic case: û(0) ] -1, 1[.

Figure 10 :

 10 Figure 10: Saturation case: Case 1: ûx (0) > 0.

Figure 11 :

 11 Figure 11: Saturation case: Case 2: ûx (0) < 0.

Figure 13 :

 13 Figure 13: Contact order 2. Non accessible points: grey region, E -: accessible. Local synthesis near E -.

Figure 14 :

 14 Figure 14: Contact order 2. Non accessible points: grey region, E + : non accessible. Local synthesis near E + .

Figure 15 :

 15 Figure 15: Stratification of the terminal manifold N for α < 1. Non-accessible points are in the filled region and on the dotted-line of E.

Figure 16 :

 16 Figure 16: Stratification of the manifold N for α > 1. Non-accessible points are in the filled region and on the dotted-line of E and points on E -are accessible.

3 β 3

 3 xv α 3 + δ 3 v(x + y) -δ 4 vvx 2 -2vxyvy 2 )(α 2 β 3 xv α 3α 2 β 4 xv α 4 + (α 2 -1)δ 3 v(x + y) + δ 4 (v -α 2 v) -α 2 vx 2 -2α 2 vxy -α 2 vy 2 -α 3 β 3 xv α 3 + α 4 β 4 xv α 4 + vx 2 + 2vxy + vy 2 ) + (α 2 β 2 xv α 2 -1 -α 4 β 4 yv α 4 -1 )((α 2 -1)β 2 v α 2 (δ 4 -(x + y)(δ 3xy)) + (α 3 -1)β 3 v α 3 (y(y -δ 3 ) + δ 4x 2 ) + (α 4 -1)β 4 yv α 4 (δ 3 -2(x + y))), • D (q) = (β 2 xv α 2 -1 -β 4 yv α 4 -1 )(α 2 β 2 xv α 2 + α 3 β 3 xv α 3 + δ 3 v(x + y) -δ 4 vvx 2 -2vxyvy 2 ) -(α 2 β 2 xv α 2 -1 -α 4 β 4 yv α 4 -1 )(β 2 xv α 2 + β 3 xv α 3 + δ 3 v(x + y) -δ 4 vvx 2 -2vxyvy 2 ),

5. 2 . 3 .Lemma 13 .

 2313 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 β 2 dv α 2 -1 + α 3 β 3 dv α 3 -1 + dδ 3 -δ 4 d 2 + y(δ 3 -2d)y 2 = 0. (35)Denoting by ∆ the discriminant of the polynomial function y → n • [G, F](d, y, v), a singularity can occur for ∆ = 0. One has Assume α i , β i , δ i > 0, i=1,2 and d, v > 0.

Figure 18 :

 18 Figure 18: Stratification of N for the system (37). Green dotted line: elliptic, red line: hyperbolic, crosses: saturating values of the singular control.

  If σ + is BC arc on [t 1 , t 2 ] and if t 2 is a switching time then t 1 cannot be a switching time if α > 1. Lemma 11. Let σ -be an extremal on [t 2 , t 3 ]. Assume ẏ(t 2 ) > 0 then ẏ(t) > 0 for t ≥ t 2 , assuming α > 1.

	Corollary 1. Proof. Assume t 1 , t 2 be switching times. Then one must have
	Φ(t 1 ) ≥ 0 and Φ(t 2 ) ≤ 0.
	A contradiction since Φ(t) > 0 from Lemma 10.
	2 y ∂ ∂y .
	21

Recall D(q) = det(G(q), [G, F](q), [[G, F] , G] (q)) , D (q) = det(G(q), [G, F](q), [[G, F] , F] (q)) , D (q) = det(G(q), [G, F](q), F(q)). [START_REF] Tóth | Reaction kinetics: exercises, programs and theorems[END_REF] and the singular control is defined by: D (q) + u s (q) D(q) = 0.

Hence, Lemma 7. We have D(q) = α(α -1)βv α-2 y, D (q) = (α -1)βv α , D (q) = (α -1)βv α y, so that,

1) The singular control is given by u s = -v 2 /(αy) and is negative. Moreover v decreases along a singular arc. The singular control saturates at S :

2) The singular trajectories are hyperbolic.

One needs also

Lemma 8. Let z(•) = (q(•), p(•)) be a BC-extremal on [t f , 0] in the coordinates q = (x, c 2 , 0) so that p(0) = (0, 1, 0). Then at any t < 0 one has p 1 (t) > 0 and p 2 (t) > 0.

Proof. The adjoint equation gives ṗ1 = -p 2 ve x , ṗ2 = p 2 βv α and the result follows.

Lemma 9. Along a BC-singular extremal one has with α > 1, ċ2 > 0 and c 2 increases.

Proof. From Lemma 8, in the coordinates q = (x, c 2 , 0), one has p 1 (t) > 0 and p 2 (t) > 0 for any t < 0 and eventually p 1 (0) = 0, p 2 (0) = 1. Using p • [G, F](q) = 0 we get the relation:

and for t < 0, p 1 > 0, p 2 > 0, hence we obtain

Since α > 1, (33) implies that the right member of (34) is positive. The result is also true at the final time if (eventually) c 2 belongs to

Properties of the switching function

In this section we use the coordinates q = (x, y, v) and we consider the switching function: We write ẋ = 1 + a y + b yz + c z 3 + ε 1 (x, y),

where

Moreover, along S we have D(q) • D (q) = -6cd 2 z and the singular control is given by u s (q) = a d 6 c 1 z .

Since the system (F, G, N) is normalized, we have

Also in a neighborhood of q 0 and along S, we can choose c > 0 so that the set S splits into hyperbolic for z < 0 and elliptic for z > 0. Fig. 18 gives the stratification of N near the q 0 for the system (37) where a = 1, c = 1, b = -3, d = 1 and |u| ≤ 5 and where the terms ε i , i = 1, 2, 3 are zeros.

Conclusion

In this article, we have presented the mathematical tools from geometric optimal control to maximize the production of chemical networks using temperature control with applications to two test-bed cases concerning two reactions schemes.


Alexander Okhotin
email: alexander.okhotin@spbu.ru

Elizaveta Sazhneva
email: sazhneva.eliza@yandex.ru

State complexity of GF(2)-concatenation and GF(2)-inverse on unary languages

The paper investigates the state complexity of two operations on regular languages, known as GF(2)-concatenation and GF(2)inverse (Bakinova et al., "Formal languages over GF(2)", LATA 2018), in the case of a one-symbol alphabet. The GF(2)-concatenation is a variant of the classical concatenation obtained by replacing Boolean logic in its definition with the GF(2) field; it is proved that GF(2)-concatenation of two unary languages recognized by an m-state and an n-state DFA is recognized by a DFA with 2mn states, and this number of states is necessary in the worst case, as long as m and n are relatively prime. This operation is known to have an inverse, and the state complexity of the GF(2)-inverse operation over a unary alphabet is proved to be exactly 2 n-1 + 1.

Introduction

Union and concatenation of formal languages are defined in terms of conjunction and disjunction: a string is in K ∪ L if it is in K or in L, and a string w is in K • L, if, for some partition w = uv, u ∈ K and v ∈ L-a disjunction of |w|+1 conjunctions. New variants of these two operations, obtained by replacing disjunctions with exclusive OR, have recently been proposed by Bakinova et al. [START_REF] Bakinova | Formal languages over GF(2)[END_REF]. Union (K ∪L) is thus replaced with symmetric difference (K L), whereas for concatenation (K • L), once the disjunction is replaced with exclusive OR, the condition of the existence of a partition turns into the condition that the number of partitions must be odd.

K • L = { w | # of partitions w = uv, with u ∈ K and v ∈ L, is non-zero } K L = { w | # of partitions w = uv, with u ∈ K and v ∈ L, is odd }
The latter operation is called GF(2)-concatenation, because it is actually a weighted concatenation with weights in the GF(2) field. For example, {ε, a} • {ε, a} = {ε, a, aa}, but {ε, a} {ε, a} = {ε, aa}, because two partitions of a cancel each other. Notably, GF(2)-concatenation is invertible: for every language L ⊆ Σ * with ε ∈ L, there exists a unique language L -1 ⊆ Σ * that Supported by Russian Science Foundation, project 18-11-00100.

Union Concatenation Star

Unambiguous ( ) mn -1 [START_REF] Jirásková | State complexity of unambiguous operations on deterministic finite automata[END_REF] (unamb•) m2 n-1 -2 n-2 [START_REF] Daley | Orthogonal concatenation: Language equations and state complexity[END_REF] (unamb * ) 3 8 2 n + 1 [6] Classical (∪) mn [START_REF] Maslov | Estimates of the number of states of finite automata[END_REF] (•) m2 n -2 n-1 [START_REF] Maslov | Estimates of the number of states of finite automata[END_REF] ( * )

3 4 2 n [9] GF(2) ( ) mn [2] ( ) m • 2 n [1]
( -1 ) 2 n + 1 [START_REF] Bakinova | Formal languages over GF(2)[END_REF] Table 1. State complexity of unambiguous, classical and GF(2)-variants of union, concatenation and star.

satisfies L L -1 = L -1 L = {ε}. For instance, {ε, a} -1 = a * , because {ε, a} -1 a * = {ε}: indeed, every non-empty string in the latter GF(2)concatenation has two partitions. Symmetric difference is the GF(2)-union.

Using GF(2)-operations instead of the classical operations gives rise to a new variant of formal language theory. For instance, GF(2)-grammars, defined by Bakinova et al. [START_REF] Bakinova | Formal languages over GF(2)[END_REF] and subsequently studied by Makarov and Okhotin [START_REF] Makarov | On the expressive power of GF(2)-grammars[END_REF], are incomparable in power to classical grammars with union and concatenation, but still have a parsing algorithm working in time O(n ω ), with ω < 3, and can be parsed by circuits of depth O((log n) 2 ).

The family of regular languages is closed under both the GF(2)-concatenation and the GF(2)-inversion operations. For a pair of languages recognized by an mstate and an n-state DFA, their GF(2)-concatenation is recognized by a DFA with m • 2 n states; this number of states is necessary in the worst case, witnessed by automata over a 2-symbol alphabet [START_REF] Bakinova | Formal languages over GF(2)[END_REF]. Similarly, the GF(2)-inverse of a language recognized by an n-state DFA is recognized by a DFA with 2 n + 1 states, and this bound is tight for alphabets containing at least 3 symbols [START_REF] Bakinova | Formal languages over GF(2)[END_REF].

To compare with the classical case, classical concatenation has state complexity m2 n -2 n-1 , and classical Kleene star, or the quasi-inverse, has state complexity 3 4 2 n [START_REF] Maslov | Estimates of the number of states of finite automata[END_REF]. Another point of comparison is with unambiguous concatenation and unambiguous star, defined by restricting the arguments, so that each string has a unique representation; classical operations and GF(2)-operations are two incomparable generalizations of the unambiguous operations. Unambiguous concatenation has state complexity m2 n-1 -2 n-2 [START_REF] Daley | Orthogonal concatenation: Language equations and state complexity[END_REF], whereas the state complexity of the unambiguous star is 3 8 2 n + 1 [START_REF] Jirásková | State complexity of unambiguous operations on deterministic finite automata[END_REF]. The state complexity of unambiguous, classical and GF(2)-variants of the three main operations on formal languages is compared in Table 1. All results refer to the case of DFA.

The goal of this paper is to investigate the state complexity of the GF(2)operations in the case of a unary alphabet [START_REF] Bakinova | Formal languages over GF(2)[END_REF]. In general, unary state complexity is substantially different from the case of multiple-symbol alphabets. The tradeoffs between different types of automata over a unary alphabet were studied by Chrobak [START_REF] Chrobak | Finite automata and unary languages[END_REF], Mereghetti and Pighizzini [START_REF] Mereghetti | Optimal simulations between unary automata[END_REF], Geffert et al. [START_REF] Geffert | Converting two-way nondeterministic unary automata into simpler automata[END_REF], Kunc and Okhotin [START_REF] Kunc | Describing periodicity in two-way deterministic finite automata using transformation semigroups[END_REF], Okhotin [START_REF] Okhotin | Unambiguous finite automata over a unary alphabet[END_REF], and others. The state complexity of operations on unary DFA was first investigated by Yu et al. [START_REF] Yu | The state complexity of some basic operations on regular languages[END_REF], who proved that concatenation is representable with mn states, and this bound is tight for relatively prime m and n; the state complexity of star on unary languages is (n -1) 2 + 1.

How do the GF(2)-operations stand in comparison? For the GF(2)concatenation of unary languages recognized by DFA with m and n states, the results generally resemble the classical case: it shall be established in Section 2 that 2mn states are sufficient, and for relatively prime m, n, this number of states is necessary. On the other hand, the case of GF(2)-inverse of a unary language is substantially different from the classical case; the state complexity turns out to be 2 n-1 + 1, which is established in Sections 3-6 by determining a connection between the states of a DFA recognizing a GF(2)-inverse and the coefficients of a certain associated sequence of polynomials over GF(2).

GF(2)-concatenation

As usual, a DFA is defined as a quintuple (Σ, Q, q 0 , δ, F ), where Σ is the input alphabet and Q is a finite set of states, with initial state q 0 ∈ Q, transition function δ : Q × Σ → Q and accepting states F ⊆ Q. This paper considers DFA over a unary alphabet Σ = {a}, where the transition function defines a sequence of states q 0 , q 1 , . . ., with q i+1 = δ(q i , a). Let j be the least number with q j = q i for some i < j. This is the point where the automaton starts to behave periodically; the states q 0 , . . . , q i-1 are called the tail, and the periodic part q i , . . . , q j-1 is called the cycle. If the tail is empty, the automaton is called cyclic.

The known construction for GF(2)-concatenation of two given DFA, A = (Σ, P, p 0 , η, E) and B = (Σ, Q, q 0 , δ, F ), works as follows [START_REF] Bakinova | Formal languages over GF(2)[END_REF]. The language L(A) L(B) is recognized by a DFA C with the set of states P × 2 Q . In a state (p, S), with p ∈ P and S ⊆ Q, the automaton simulates the computation of A in the first component p, while S is the set of all states reached an odd number of times in the ongoing simulated computations of B. The initial state of C is (p 0 , {q 0 }) if ε ∈ L(A) and (p 0 , ∅) otherwise. Its transition function, π : (P × 2 Q ) × Σ → P × 2 Q , is defined on a pair (p, S) as follows. Each currently simulated computation of B is continued, represented by the following set S . S = { q | the number of states q ∈ S, with q = δ(q, a), is odd } If the simulated automaton A passes through an accepting state, then the transition π((p, S), a) also adds a new computation to the set S .

π((p, S), a) = (η(p, a), S ), if η(p, a) / ∈ E (η(p, a), S {q 0 }), if η(p, a) ∈ E a state (p, S) is marked as accepting, if S contains an odd number of accepting states of B. F = (p, S) |S ∩ F | is odd
This completes the known construction, which is valid for every alphabet.

In the case of a unary alphabet, the state complexity of GF(2)-concatenation on two unary DFA is first investigated in the special case of both automata being cyclic, each with a unique accepting state. Under these restrictions, the state complexity depends only on the number of states in the given automata. Lemma 1. Let A and B be two minimal cyclic DFA over a unary alphabet, one with the cycle of length m and the other with the cycle of length n, and each with a single accepting state. Then, the GF(2)-concatenation L(A) L(B) is recognized by a cyclic DFA with period of length 2mn gcd(m,n) , and this is the minimal DFA for this language.

Proof. Let A = ({a}, P, 0, η, E), with P = {0, . . . , m -1}, η(i, a) = i + 1 mod m for all i ∈ P , and Assume that e = 0. After reading a e , the automaton reaches the state (e, {0}). From this point on, consider the states reached by the automaton after reading repetitive blocks a m . The state in the first component is rejecting until the last symbol in the block, and hence the states in both components are cyclically shifted by m, until the accepting state e reappears in the first component. At the last step, a new state 0 is added to the second component, while all pre-existent states in the second component have been shifted by m (mod n). Altogether, the following states are visited. π((0, ∅), a e ) = (e, {0}) π((e, {0}), a m ) = (e, {0,

E = {e}. Similarly, let B = ({a}, Q, 0, δ, F ), with Q = {0, . . . , n -1}, δ(i, a) = i + 1 mod n for all i ∈ Q, and F = {f }. Let C = (Σ, Q , q 0 , π, F )
k 1 }),
where

k 1 ≡ m (mod n) π((e, {0, k 1 }), a m ) = (e, {0, k 1 , k 2 })
where k 2 ≡ 2m (mod n), etc.

The subset continues to grow until some j-th step, with δ(k j , a m ) = 0. This means that (j + 1)m ≡ 0 (mod n). Since j is the least such number, it must be j = n gcd(m,n) -1. Therefore, π((e, ∅), a jm ) = (e, {0, k 1 , k 2 , . . . k j }).

From this point on, consider further computations upon reading repetitive blocks a m . The states in the second component keep cyclically shifting, and the states 0 added in the end of each block cancel out these states in the same order as they were added.

π((e, {0, k 1 , k 2 , . . . k j }), a m ) = (e, {k 1 , k 2 , . . . k j }) π((e, {k 1 , k 2 , . . . k j }), a m ) = (e, {k 2 , k 3 , . . . k j }),
etc.

In the end, π((e, {0, k 1 , k 2 , . . . k j }), a jm ) = (e, ∅). Finally, after reading a m-e , the automaton returns to its initial state.

π((e, ∅), a m-e ) = (0, ∅)

Overall, the automaton for the GF(2)-concatenation is a cycle on 2mn gcd(m,n)

states.

Claim. The period 2mn gcd(m,n) is minimal.

The idea of the argument is that C has a block of m + n consecutive rejecting states, and that this block occurs only once in the automaton.

Theorem 1. Let A and B be any two DFA over a unary alphabet, with m and with n states, respectively. Then, the GF(2)-concatenation L(A) L(B) is recognized by a DFA with 2mn states.

For relatively prime m and n, this number of states is necessary in the worst case.

Proof (a sketch). Let A have the set of states P = {0, 1, . . . , m -1}, with accepting states E ⊆ P . For each state i ∈ P , define a DFA A i by setting i in A as the only accepting state. Then, L(A) = i∈E L(A i ), and the union is disjoint. Let the language be periodic starting from k, with period m -k.

Similarly, let the set of states of B be Q = {0, 1, . . . , n -1}, with accepting states F ⊆ Q. Let B j be B with j as the only accepting state, so that L(B) = i∈F L(B i ). Let the periodic part begin at , with period n -.

Then, the desired GF(2)-concatenation can be represented as follows.

L(A) L(B) = i∈E L(A i ) j∈F L(B j ) = i∈E L(A i ) j∈F L(B j ) = = i∈E j∈F L(A i ) L(B j )
Each of these |E| • |F | languages is periodic beginning from k + -1, with period 2(m -k)(n -); the proof is omitted due to space constraints. These languages are then joined into a single automaton with at most 2mn states.

For the lower bound, Lemma 1 with relatively prime m, n provides the desired witness languages.

Automaton for GF(2)-inverse

With respect to GF(2)-concatenation, every language L containing the empty string is invertible, in the sense that there is a language L -1 satisfying L L -1 = L -1 L = {ε}. The GF(2)-inverse operation, f (L) = L -1 , preserves regularity, and its state complexity is 2 n + 1 [START_REF] Bakinova | Formal languages over GF(2)[END_REF].

Theorem A (Bakinova et al. [START_REF] Bakinova | Formal languages over GF(2)[END_REF]Thm. 2]) For every language L over an alphabet Σ, with ε ∈ L, a string w ∈ Σ * is in L -1 if and only if it has an odd number of representations of the form w = w 1 w 2 . . . w k , with k 0 and w 1 , . . . , w k ∈ L \ {ε}.

As proved by Bakinova et al. [START_REF] Bakinova | Formal languages over GF(2)[END_REF], for every n-state DFA A = ({a}, Q, q 0 , δ, F ), with ε ∈ L(A), the language L(A) -1 is recognized by a DFA C = ({a}, 2 Q ∪ {q 0 }, q 0 , δ , F ) defined as follows. The states of C are all subsets of Q and a new initial state q 0 , Its transition function is δ : (2 Q ∪ {q 0 }) × {a} → 2 Q ∪ {q 0 }. The transition in the state q 0 produces a singleton state corresponding to a single computation of A.

δ (q 0 , a) = {δ(q 0 , a)} In a state S ⊆ Q, first let S = { q | # of states p ∈ S with δ(p, a) = q is odd }.
Then the transition is defined as follows.

δ (S, a) = S , if |S ∩ F | is even S δ(q 0 , a), if |S ∩ F | is odd
The set of accepting states is In the general case, the DFA for the GF(2)-inverse L(A) -1 of an n-state DFA A has 2 n + 1 states, and it is known that this number is necessary in the worst case, for alphabets with at least three symbols [START_REF] Bakinova | Formal languages over GF(2)[END_REF]. It turns out that in the unary case it is always sufficient to use 2 n-1 + 1 states.

F = S |S ∩ F | is odd ∪ {q 0 }.
This upper bound is easy to establish for non-cyclic automata.

Lemma 2. Let A = ({a}, Q, 0, δ, F ) be an n-state non-cyclic DFA with 0 ∈ F . Then, the DFA for the GF(2)-inverse L(A) -1 constructed as above, has at most 2 n-1 + 1 reachable states.

Proof. Indeed, no subset containing the state 0 is ever reached, since this state is not reachable by any transitions.

For cyclic automata A, a deeper analysis of the automaton for its inverse is needed, since the set of unreachable states is harder to specify. The first result to be established is the following dependence between the membership of individual states in the subsets. be the DFA recognizing the GF(2)-inverse of L(A), defined as above. For each i 1, let S i ⊆ Q be the state of C reached upon reading the string a i .

Lemma 3. Let A = ({a}, Q, 0, δ, F ) be a cyclic DFA with Q = {0, . . . , n -1}, δ(i, a) = i + 1 mod n for all i, and 0 ∈ F . Let C = ({a}, 2 Q ∪ {q 0 }, q 0 , δ , F ) {0} ∅ {2} {3} {1,4} {0,1,2} {0,1} {1} {0,2} {0,3} {0,1,4} {1,2} q 0 ' {2,3} {0,2,3} {0,1,3,4} {2,4} {1,3} {0,1,2,4} {1,2,3} {0,1,2,3,4} {2
,3,4} {3,4} {4} {1,3,4} {0,2,4} {0,1,3} {1,2,4} {0,1,2,3} {1,2,3,4} {0,2,3,4} {0,3,4} {0,4}
Denote the membership of the j-th state in S i by a Boolean value S j i ∈ {0, 1}, with S j i = 1 if j ∈ S i , and S j i = 0 otherwise. Then, the membership of state 1 in each set S i is determined by the set S i-1 by the following formula.

S 1 i = f ∈F \{0} S f i-1 (for i 2)
Furthermore, the membership of 0 in S i depends on its membership in the previous n -1 states as follows.

S 0 i = f ∈F \{0} S 0 i-f (for i n + 1)
Proof (a sketch). The formula for S 1 i is directly inferred from the definition of automata. The second formula is inferred from this one using the following two observations: first, 0 is in S i if and only if 1 is in S i-(n-1) ; second, a state f is in S i-n if and only if 0 is in S i-f .

Polynomials for GF(2)-inverse

The subsets reached by the automaton for the GF(2)-inverse of a cyclic language L have a useful characterization in terms of certain polynomials over GF [START_REF] Brzozowski | Quotient complexity of regular languages[END_REF]. First, based on the automaton for L, a sequence of polynomials shall be constructed, and then the leading coefficients of these polynomials shall correspond to the membership of state 0 in the subsets of the automaton for L -1 .

A few definitions are due. Let f (x) = a n-1 x n-1 + . . . + a 1 x + a 0 be a polynomial of degree n -1 over GF [START_REF] Brzozowski | Quotient complexity of regular languages[END_REF], with a 0 , . . . , a n-1 ∈ {0, 1} and with a n-1 = 1.

For every i 0, let p i (x) be the polynomial obtained by taking x i modulo f (x); this is a polynomial of degree at most n -2. The polynomials p i (x) form a sequence, in which the first term is p 0 (x) = 1, and every succeeding term is obtained from the previous term as follows: if a polynomial p i (x) does not contain the monomial x n-2 , then the next term is p i+1 (x) = x • p i (x); and if there is a monomial x n-2 in p i (x), then the next term is p i+1 (x) = x • p i (x) + f (x); in other words, x n-1 is replaced with a n-2 x n-2 . . . + a 1 x + a 0 .

The first n terms in the sequence {x i mod f (x)} i=0 are 1, x, x 2 , . . . , x n-2 , x n-1 + f (x). The form of the subsequent terms non-trivially depends on f . A polynomial f (x) of degree n-1 is called primitive, if this sequence contains all 2 n-1 -1 non-zero polynomials. For each n 2, primitive polynomials are known to exist.

For all i 0 and j ∈ {0, . . . , n -2}, let b i,j ∈ {0, 1} be the coefficient at the monomial x j in p i (x). The coefficient at the term x n-2 in p i (x) depends on the coefficients at the same term in the preceding n -1 polynomials in the sequence.

i,n-2 = n-2 k=0 b i-(n-1-k), n-2 a k
Let A = ({a}, Q, 0, δ, F ) be a cyclic unary DFA with the set of states Q = {0, . . . , n -1} and with 0 ∈ F . The corresponding polynomial over GF(2) is defined as f (x) = j∈F x n-1-j . In the sequence of polynomials x i modulo f (x), let b i,n-2 be the coefficient at x n-2 in the i-th polynomial. Then, by Lemma 4, each coefficient is expressed through the preceding n -1 coefficients as follows.

b i,n-2 = n-2 k=0 b i-(n-1-k),n-2 • (n -1 -k ? ∈ F ) = j∈F \{0} b i-j,n-2
This is the same recurrent formula as in Lemma 3.

Example 2 (continued from Example 1). For the 5-state cyclic automaton with accepting states F = {0, 3, 4}, the corresponding polynomial is f (x) = x 4 + x + 1. Then, the sequence x i modulo f (x) begins with the following polynomials.

p 0 (x) = 1, p 1 (x) = x, p 2 (x) = x 2 , p 3 (x) = x 3 , p 4 (x) = x + 1, p 5 (x) = x 2 + x, etc.
The following table puts the subsets reachable by the automaton alongside the polynomials in this sequence. The acceptance status of each subset is provided for reference.

i S i+2 p i (x) a i+2 ? ∈ L -1 0 {2} 1 - 1 {3}
x + 2 {1, 4}

x 2 + 3 {0, 1, 2} x 3 + 4 {2, 3}

x + 1 + 5 {1, 3, 4} x 2 + x -6 {0, 2, 4} x 3 + x 2 -7 {0, 1, 3} x 3 + x + 1 -8 {1, 2, 4} x 2 + 1 + 9 {0, 1, 2, 3} x 3 + x -10 {1, 2, 3, 4} x 2 + x + 1 -11 {0, 2, 3, 4} x 3 + x 2 + x + 12 {0, 3, 4} x 3 + x 2 + x + 1 + 13 {0, 4}

x 3 + x 2 + 1 -14 {0, 1}

x 3 + 1 +

For every i, the state 0 is in the subset S i+2 if and only if the polynomial contains the term x 3 . This is not a coincidence, and this correspondence shall now be established in the general case.

Upper bound for the GF(2)-inverse

The following two binary sequences turn out to be identical. First, there is the sequence {S 0 i } representing the membership of the state 0 in the subsets reached by the automaton for the GF(2)-inverse. The other sequence is the sequence {b i,n-2 } of coefficients at x n-2 . Lemma 5. Let A = ({a}, Q, 0, δ, F ) be a cyclic DFA with Q = {0, . . . , n -1} and 0 ∈ F . For each i 1, let S i ⊆ Q be the state of the automaton for the inverse given in Section 3, reached upon reading the string a i . For each i 1, let p i (x) be x i taken modulo f (x) = j∈F x n-1-j . Then, for every i 0, the state 0 is in S i+2 if and only if the monomial x n-2 is in p i (x).

The proof is by induction on i: the base cases are i ∈ {0, 1, . . . , n -2}, on which the sequences coincide. The induction step follows by Lemmata 3 and 4, since both sequences are defined by the same formulae on the same data. Lemma 6. Assume that the sequence {S 0 i } ∞ i=0 has period p beginning at , in the sense that S 0 i = S 0 i+p for all i . Then, the sequence of states {S i } ∞ i=0 has period p beginning at .

  be the DFA recognizing the GF(2)-concatenation L(A) L(B), with the set of states Q = { (p, S) | p ∈ P, S ⊆ Q }, and with accepting states F = { (p, S) | p ∈ P, f ∈ S }. Claim. For all p ∈ P and k ∈ {1, . . . , n gcd(m,n) -1}, there are exactly two reachable subsets S 1 , S 2 ⊆ Q with |S 1 | = |S 2 | = k and π((p, S 1 ), a) = π((p, S 2 ), a). For k = 0 or k = n gcd(m,n) , there is a unique reachable subset (p, S) with (p, S) ∈ Q and |S| = k. Accordingly, Q contains 2mn gcd(m,n) reachable states.

Example 1 .

 1 Consider the following 5-state unary DFA. its inverse, constructed by the above method, is shown in Figure 1. It has a cycle of length 15 and a tail of length 2 (along with 16 unreachable states).

Fig. 1 .

 1 Fig. 1. DFA for the GF(2)-inverse of the language in Example 1.

Lemma 4 .

 4 For every i n -1, the coefficient b i,n-2 depends on the earlier coefficients as follows, with all arithmetic in GF(2).

  b

 Lemma 7. Let A = ({a}, Q, 0, δ, F ) be an n-state cyclic DFA with 0 ∈ F . Then, the DFA for the GF(2)-inverse L(A) -1 , constructed as in Section 3, has at most 2 n-1 + 1 reachable states.

Proof (a sketch). The sequence of polynomials x i modulo f (x) = j∈F x n-1-j contains at most 2 n-1 -1 distinct polynomials. By Lemma 5, this sequence coincides with the sequence of state 0, and then, by Lemma 6, the sequence of subsets has the same total length of the tail and the period. This gives the desired upper bound on the number of states.

Lower bound for the GF(2)-inverse

The lower bound shall be established using cyclic witness languages. The following property of GF(2)-inverses of cyclic unary languages comes useful. Lemma 8. Let L ⊆ a * , with ε ∈ L, be a unary language recognized by an n-state cyclic DFA. Then, its GF(2)-inverse L -1 contains a string a , with n + 1, if and only if the number of representations a = a -j a j , with a -j ∈ L -1 , a j ∈ L \ {ε} and j < n, is odd.

The proof is by establishing the equivalence with the condition in Theorem A. For any language L ⊆ a * , with ε ∈ L, let α i = 1 if a i ∈ L -1 , and α i = 0 otherwise. Then the condition in Lemma 8 can be written down as the following formula.

Lemma 9. Let L ⊆ a * be a language, with ε ∈ L, recognized by an n-state cyclic DFA. Then, α i = j∈F \{0} α i-j for i n + 1.

As in Section 4, let f (x) = a n-1 x n-1 +. . .+a 1 x+a 0 be a primitive polynomial over GF [START_REF] Brzozowski | Quotient complexity of regular languages[END_REF], with a n-1 = a 0 = 1 (primitive polynomials of any degree over GF(2) are known to exist). For every i 0, let

Since f is primitive, by definition, all polynomials p 0 (x), . . . , p 2 n-1 -2 are pairwise distinct, and then p 2 n-1 -1 = p 0 (x) = 1.

It turns out that the sequence of coefficients at x n-2 has the same period as the sequence of full polynomials.

Lemma 10. The minimal period of the sequence {b

Lemma 11. The sequence {b i,n-2 } ∞ i=0 contains all binary substrings of length n -1, except for (0, . . . , 0). A cyclic automaton A f corresponding to this primitive polynomial f (x) = a n-1 x n-1 +. . .+a 1 x+a 0 is defined as A f = ({a}, Q, 0, δ, F ), with Q = {0, . . . , n-1}, δ(i, a) = i + 1 mod n for all i, and F = { n -1 -i | a i = 1 } (by the same principle as in Lemma 5).

Let L = L(A f ), and consider the sequence {α i } ∞ i=0 defined as above. The goal is to prove that (α 2 , . . . , α n ) = (0, . . . , 0).

Sum Concatenation Star

Unambiguous ( )

Table 2. State complexity of unambiguous, classical and GF(2)-variants of sum, concatenation and star: the case of a unary alphabet.

Lemma 12. Let L be a unary language with ε ∈ L and L = {ε}, a * , which is recognized by a DFA with n states. Then, the inverse L -1 contains a string of length between 2 and n.

Since the sequence {α i } ∞ i=2 begins with something other than n-1 zeroes, by Lemma 10, the binary substring (α 2 , . . . , α n ) occurs somewhere in the sequence {b i,n-2 } ∞ i=0 . By Lemma 9, the rest of the terms of the sequence {α i } are defined by the same formula as the sequence {b i,n-2 }, which makes the binary strings α 2 , . . . , α 2 n-1 and b 0,n-2 , . . . , b 2 n-1 -2 identical up to a cyclic shift. In particular, the period of the sequence {α i } ∞ i=2 is 2 n-1 -1. It remains to determine the length of the tail. Since the construction in Section 3 produces 2 n-1 + 1 states, the length of the tail is at most 2. It turns out that it cannot be shortened, because the strings a and a 2 n-1 -1 have different membership status. Lemma 13. α 1 = α 2 n-1 +1 , and therefore the length of the tail is 2.

The following theorem has thus been established. Theorem 2. For every n 2, there exists a language L, with ε ∈ L, recognized by n-state unary cyclic DFA, for which the minimal DFA recognizing its GF(2)inverse L -1 has 2 n-1 + 1 states.

Future work

The results of this paper are summarized and compared to related results in Table 2.

A problem proposed for future research is to determine the number of states in NFA needed to represent these operations. There are two different modes of nondeterminism involved: existential nondeterminism in NFA and parity nondeterminism in both GF(2)-operations. Intuitively, one kind of nondeterminism cannot help implementing another kind, and the following straightforward construction might actually turn out to be the best possible: first, determinize the arguments, with a blow-up of the order e (1+o(1)) √ n ln n [START_REF] Chrobak | Finite automata and unary languages[END_REF]; and then, apply the constructions for deterministic automata presented in this paper. Could this construction be substantially improved upon?