Towards Geometric Time Minimal Control without Legendre Condition and with Multiple Singular Extremals for Chemical Networks
11Towards Geometric Time Minimal Controlwithout Legendre Condition and with MultipleSingular Extremals for Chemical Networks
Bernard Bonnard and J´er´emy Rouot Inria Sophia Antipolis and Institut de Math´ematiques de Bourgogne, 9 avenueSavary, Dijon, 21078 France, [email protected] [email protected] Summary.
This article deals with the problem of maximizing the production ofa species for a chemical network by controlling the temperature. Under the so–called mass kinetics assumption the system can be modeled as a single–input controlsystem using the Feinberg–Horn–Jackson graph associated to the reactions network.Thanks to Pontryagin’s Maximum Principle, the candidates as minimizers can befound among extremal curves, solutions of a (non smooth) Hamiltonian dynamicsand the problem can be stated as a time minimal control problem with a terminaltarget of codimension one. Using geometric control and singularity theory the timeminimal syntheses (closed loop optimal control) can be classified near the terminalmanifold under generic conditions. In this article, we focus to the case where thegeneralized Legendre-Clebsch condition is not satisfied, which paves the road tocomplicated syntheses with several singular arcs. In particular, it is related to thesituation for a weakly reversible network like the McKeithan scheme of two reactions:T+M A B .
Keywords:
Mass action chemical systems · Pontryagin Maximum Principle · Geometric optimalcontrol · Time minimal synthesis · McKeithan type chemical scheme
Mathematics Subject Classification: · · The optimization of the production is an important problem in chemical and bio-logical engineering and the control can be either the temperature (batch or closed a r X i v : . [ m a t h . O C ] J a n Bernard Bonnard and J´er´emy Rouotreactor) or by feeding the reactor (semi-batch or open case). In this article we shallconcentrate on the first case. Moreover we assume that the dynamics is modeledusing the mass action kinetics assumptions and hence given at constant tempera-ture T by a polynomial system based only on the Feinberg–Horn–Jackson graphassociated to the chemical network. Also in the 70’s those researchers obtained (un-der the so–called zero deficiency assumption) in a series of seminal articles [17, 18]a complete description of the dynamics, at constant temperature. If this dynamicsis well understood, in the optimal problem the temperature is not constant and theanalysis becomes very intricate. Thanks to the Pontryagin Maximum Principle [24]candidates as minimizers can be found among extremals solutions of a ( non smooth )Hamiltonian dynamics and optimal solutions are concatenation of bang arcs, withminimum and maximum temperature, and the so–called singular arcs, defined asa solution of a smooth Hamiltonian constrained dynamics [5]. Moreover maximi-zing the production of one species during the batch can be restated as producinga 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 loopfeedback and this leads to the problem of computing the time minimal synthesis , fora single–input control system. At the end of the 80’s, geometric optimal control hasproduced an important literature to compute the time minimal syntheses for the fixedend point case , where the initial point is localized near the terminal point. This wasdone using the
Lie algebraic structure of the control system, mainly for single–input(smooth) system, where the control appears linearly, some seminal references are[28, 26, 12] either in general context or in view of application formed by a sequenceof two irreversible reactions: A → B → C and in relation with an industrial project[8]. Our aim is to extend this work to more complicated reaction schemes and to dealin particular with weakly reversible chemical schemes. More precisely we shall con-centrate on a McKeithan type scheme of the form T+M A B assumingthat the coefficient governing the dynamics are given by Arrhenius law. This schemewas already studied in the context of control theory using stabilization techniqueswith “ feeding ” types control [27] and ad hoc observer design [14]. In our case, thisnetwork is a test bed case for our very general approach.The key point for this extension is the analysis of singular trajectories and theirrole in the synthesis. This is connected with an important question and the needto extend the standard synthesis related to the turnpike phenomenon [29] to dealwith cases, where the strict Legendre–Clebsch condition is not satisfied, a situationencountered in a recent application in MRI [6].The organization of this article is the following. In Section 2, we recall brieflythe Feinberg–Horn–Jackson theory to model the dynamics of chemical networks ofconstant temperature and the properties of the dynamics under the zero deficiencyassumption [17, 18], which can be applied to the McKeithan scheme. The stabilityproperties are recalled and can be applied to control stabilization and observerdesign [27, 14]. In Section 3, we present the fundamental results of the time minimalcontrol problem, which are relevant to our study: Pontryagin Maximum Principle[24], regular and singular extremals [5],[20]. The general turnpike theorem [29] is Geometric Techniques to Optimize Chemical Reactions 3recalled and extensions are presented in relation with the problem with terminalmanifold of codimension one and when the strict Legendre–Clebsch condition isnot satisfied. The concept of conjugate and focal points is introduced based on[7]. In Section 4, we analyze the time minimal control problem for a two reactionsMcKeithan scheme. To compute the time minimal syntheses, we present techniquesand results from [8, 9], which have to be extended to analyze the problem. Thecomputational complexity of the problem is discussed and symbolic computationsare presented to cope with this complexity. We consider a set of m chemical species { X , . . . , X m } and the state of the dynamicsis the vector c = ( c , . . . , c m ) (cid:124) ∈ R m ≥ representing the molar concentration. Let R be a set of reactions, each reaction being denoted by y → y (cid:48) and of the form m (cid:88) i =1 α i X i −→ m (cid:88) i =1 β i X i , where α i , β i are the stoichiometric coefficients and the vectors y = ( α , . . . , α m ) (cid:124) and y (cid:48) = ( β , . . . , β m ) (cid:124) are the vertices of the so–called Feinberg–Horn–Jacksonoriented graph associated to the network, edges being oriented according to y → y (cid:48) .Each reaction is characterized by a reaction rate K ( y → y (cid:48) ) and the system is said simple (or mass kinetics) if the rate of the reaction is of the form: K ( y → y (cid:48) ) = k ( T ) c y ,c y = c α . . . 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, bothdepending on the reaction, R is the gas constant and T is the temperature. Note thatdifferent rate formulae can be used to deal in particular with biomedical systems(see for instance [27]). The dynamics of the system, taking into account the wholenetwork is: ˙ c ( t ) = f ( c ( t ) , T ) = (cid:88) y → y (cid:48) K ( y → y (cid:48) ) ( y (cid:48) − y ) . (1.2) Definition 1.
The stoichiometric subspace is S := span { y − y (cid:48) ; y → y (cid:48) ∈ R} and the sets ( c (0) + S ) ∩ R m ≥ are called the (strictly if > ) positive stoichiometriccompatibility classes. From [2] we have. Bernard Bonnard and J´er´emy Rouot
Lemma 1.
Let c ( t ) be a solution of (1.2) with initial condition c (0) ∈ R m ≥ . Then c ( t ) belongs for all t ≥ to the strictly positive compatibility class ( c (0) + S ) ∩ R m ≥ . Definition 2.
Having labeled the set of vertices by i = 1 , . . . , n , with correspondingstoichiometric vector ( y , . . . , y n ) , the complex matrix is Y := ( y , . . . , y n ) . The in-cidence connectivity matrix A := ( a ij ) contains the Arrhenius coefficients k i of thereactions using the rule: k = a indicates a reaction with kinetics constant k fromthe first node to the second, that is y → k y . With the mass kinetics assumption, the dynamics can be expressed as˙ c ( t ) = f ( c ( t ) , T ) = Y ˜ Ac Y , where ˜ A is the Laplacian matrix in graph theory defined by˜ A = A − diag (cid:32) n (cid:88) i =1 a i , . . . , n (cid:88) i =1 a in (cid:33) (1.3)and we denote c Y = ( c y , . . . , c y n ) (cid:124) . It is given by the reaction scheme: T + M C C . . . C N k k p, k p, k p,N k − , k − , k − ,N . The matrix Y is given by Y = . . . . . .
01 0 · . . . . . . and the matrix A = ( a ij ) is defined by a = k , a i = k − ,i − , i = 2 , . . . , m ( m = N + 2), a i,i − = k p,i − , i = 3 , . . . , m and all others a ij are zero.The stoichiometric subspace is defined by: { c : T + C + . . . + C N = M + C + . . . + C N = 0 } and we note: δ = T + C + . . . + C N and δ = M + C + . . . + C N the constantsassociated to first integrals of the dynamics. Geometric Techniques to Optimize Chemical Reactions 5Consider the special case N = 2 so that the reaction scheme is denoted byT+M A B k k k k , and restricting to the stoichiometric class ( δ and δ beingfixed), one gets with: x := [ A ] , y := [ B ], [ · ] denoting the respective concentrations,that the dynamics is described by the equations˙ x = k ( δ x − y )( δ − x − y ) − ( k + k ) x ˙ y = k x − k y. (1.4) Definition 3.
The deficiency of the network is: δ = n − l − s , where n is the numberof vertices, l is the number of connected components and s is the dimension of thestoichiometric 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 apath joining j to i . Using [17], refined by [2, 27], one has the following result.
Theorem 1.
The graph associated to the McKeithan scheme is strongly connectedand with deficiency zero. In each strictly positive compatibility class there exists inthis domain an unique globally asymptotically stable equilibrium.
This stability result has consequences to control and observation properties of thenetwork, see [27, 15], that we recall briefly. The dynamics (1.4) can be convertedinto a control system of the form˙ c ( 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 classwill allow to get stabilization result for a single equilibrium. Moreover it leads todesign under a mild assumption (detectability) a simple observer. We refer to [27]and [14] for the detailed presentation of those results, the geometric constructionbeing clear. The system is written as d c d t = f ( c, T ) (see (1.4)) and controlling the temperatureleads to T ∈ [ T m , T M ]. In the sequel, we shall use the terminology direct for thecorresponding control problem. In practice, thermodynamics has to be used to model Bernard Bonnard and J´er´emy Rouotthe heat exchanges, in relation with the heat produced by the reactions [31] orthe heat exchange device used in the experiments, depending upon the technicalachievements. To avoid this part of the study and without losing any mathematicallygenerality, we shall use ˙ v as the control variable setting ˙ v = u , where v := k i ( T ) forsome reaction i , where k i ( T ) = A i e − E i / ( RT ) (see (1.1)). This leads to deal with theso–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]. Note that the bounds v ∈ [ v m , v M ] will not be taken into account in our study. The map v (cid:55)→ ˙ v is thestandard Goh transformation in optimal control, see [5].The optimal control problem of physical interest is the problem of maximizingthe production of one species and using a proper variable labeling, the optimalproblem is therefore of the
Mayer type :˙ q = F ( q ) + u G ( q ) , max | u |≤ q ( t f ) , where t f is the time duration of the batch and q is the desired product.Note it will be show (thanks to the Maximum Principle) that a “dual” formula-tion is min u ( · ) t f , c ( t f ) = d, where d > X during the batch. 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 )] (cid:55)→ [ − , U is the set of bounded measurable mappings. If q (0) = q (initial state), we denote by q ( · , q , u ) (in short q ( · )) the solution starting from q .Fixing t f , the accessibility set in time t f is the set A ( q , t f ) = ∪ u ( · ) ∈U q ( t f , q , u ). The extremity mapping (in time t f ) is the map: E q ,t f : u ( · ) (cid:55)→ q ( t f , q , u ) defined on adomain of U ; the set U is endowed with the L ∞ -norm topology. Statement in the time minimal case with q ( t f ) ∈ N := smooth terminalmanifold. Notation: H ( q, p, u ) = p · X ( q, u ) denotes the pseudo–Hamiltonian (Hamiltonian liftof the vector field X ), p is the adjoint vector in R n \ { } and · is the scalar product.We denote by M ( q, p ) = max | u |≤ H ( q, p, u ). Geometric Techniques to Optimize Chemical Reactions 7 Statement of the PMP:
If ( q ∗ , u ∗ ) is an optimal control–trajectory pair on [0 , t ∗ f ]then there exists p ∗ ∈ R n \ { } such that a.e.:˙ q ∗ ( t ) = ∂H∂p ( q ∗ ( t ) , p ∗ ( t ) , u ∗ ( t )) , ˙ p ∗ ( 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 transversalitycondition p ∗ ( t f ) ⊥ T q ∗ ( t f ) N. (1.7) Statement in the Mayer case:
As before, except that the transversality condition isreplaced 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) = q and p ( t f ) satisfies the transversality condition. An extremalcontrol is called regular if | u ( t ) | ≤ a.e. and singular if ∂H∂u = 0 everywhere. Anextremal is said exceptional if M = 0 . A regular extremal control is called ”bang–bang” if u ( · ) is piecewise constant on [0 , t f ] (i.e. the number of switches is finite). 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 themaximization condition, the Legendre–Clebsch condition ∂ H∂u ≤ z (cid:55)→ u s ( z ) and plugging such u s 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 H X ( z ) = p · X ( q ) , z = ( q, p ) andthe Poisson bracket is defined by { H X , H Y } ( z ) = p · [ X, Y ]( q ). Computation of singular extremals:
The condition ∂∂u dd t ∂H∂u = {{ H G , H F } , H G } ≥ >
0) is called the (resp.strict) generalized Legendre–Clebsch condition .Recall the following.
Proposition 1 ([19]).
The generalized Legendre–Clebsch condition is a necessaryoptimality condition for the time minimal control problem with fixed extremities.
To compute the singular extremal, we differentiate twice t (cid:55)→ H G ( z ( t )) and we get 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 . (1.8)Assume (in relation with the Legendre–Clebsch condition) {{ H G , H F } , H G } (cid:54) = 0.The corresponding extremal is called of order 2 and the singular control u is com-puted as u s ( z ), using relation (1.8). Plugging such u s ( z ) into H ( q, p, u ) leads todefine the true singular Hamiltonian, denoted by H s ( z ). One has: Bernard Bonnard and J´er´emy Rouot Proposition 2.
Singular extremals of order are the solutions of : d q d t = ∂H s ∂p ( z ) , d p d t = − ∂H s ∂q ( z ) (1.9)with the constraints { H G , H F } ( z ) = H G ( z ) = 0 . Moreover, in order to be admissible, the singular control is given by u s ( z ) = − {{ H G , H F } , H F } ( z ) {{ H G , H F } , H G } ( z ) (1.10)and has to satisfy the admissibility constraint | u s ( z ) | ≤ Definition 5.
Let z = ( q, p ) be a singular extremal of order and M = H F = h theconstant value of the Hamiltonian. The extremal is called exceptional if h = 0 . If h > , the extremal is called hyperbolic (resp. elliptic) if {{ H G , H F } , H G } > (resp. < ). Recall that for our network ˙ c = f ( c, v ) , v = k i and it is extended into˙ q = F ( q ) + u G ( q ) with F = f ( c, v ) ∂∂c and G = ∂∂v . Denote ˜ H = p c · f ( c, v ) and H = p · ( F + u G ) the respective Hamiltonian lifts, with p = ( p c , p v ). One has thefollowing relation between the corresponding singular extremals. Lemma 2.
The pair ( q, p ) is solution of ˙ q = ∂H∂p , ˙ p = − ∂H∂q , ∂H∂u = 0 if and only if p v = 0 and ( c, p c , p v ) is solution of ˙ c = ∂ ˜ H∂p c , ˙ p c = ∂ ˜ H∂ ˜ c , ∂ ˜ H∂v = 0 , and moreover the following relations are satisfied dd t ∂H∂u = { H G , H F } = − ∂ ˜ H∂v (1.11) ∂∂u d d t ∂H∂u = {{ H G , H F } , H G } = − ∂ ˜ H∂v . (1.12)In particular this gives the correspondence between both singular Hamiltonianand the respective Legendre–Clebsch and generalized Legendre–Clebsch conditions. Geometric Techniques to Optimize Chemical Reactions 9 n = 3 We have q = ( x, y, v ) and introduce the following determinants: D = det( G, [ G, F ] , [[ G, F ] , G ]) ,D (cid:48) = det( G, [ G, F ] , [[ G, F ] , F ]) ,D (cid:48)(cid:48) = det( G, [ G, F ] , F ) . Using p · G = p · [ G, F ] = p · ([[ G, F ] , F ] + u [[ G, F ] , G ]) = 0 and eliminating p leadsto the following. Proposition 3. If n = 3 , the singular control is given by the feedback u s ( q ) = − D (cid:48) ( q ) /D ( q ) and singular trajectories are defined by the vector field X s ( q ) = F ( q ) − D (cid:48) ( q ) D ( q ) G ( q ) . Singular trajectories are hyperbolic if DD (cid:48)(cid:48) > , ellipticif DD (cid:48)(cid:48) < and exceptional if D (cid:48)(cid:48) = 0 . Optimality status: the case n = 3. We use [7] to describe the optimality statusof singular trajectories. The system is ˙ q = F + u G and we relax the bound | u | ≤ u ∈ R so that singular arcs are admissible. We assume the following:( H ) F, G are linearly independent,( H ) G, [ G, F ] are linearly independent.One picks a smooth singular arc z ( · ) = ( q ( · ) , p ( · )) defined on [0 , t f ] so that p isunique up to a non zero multiplicative scalar and q ( · ) is a one-to-one immersion. Wehave: Proposition 4.
There exists a C –neighborhood U of the reference singular arc t (cid:55)→ 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 t c with respect to all trajectories withthe same extremities contained in U . In the exceptional case, the reference singulararc is time minimal (and time maximal). Algorithm to compute the first conjugate point in the hyperbolic–ellipticcase.
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 t c is the first t > V ( t ) is collinear to G ( q ( t )).See [5, p. 123] for a proof and the geometric interpretation. In this section, we recall the seminal results coming from singularity theory due to[16, 20] 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 = − ) and σ s an admissible singular arc. We denote by σ σ an arc σ followed by σ . The surface Σ : H G ( z ) = 0 is called the switching surface and let Σ (cid:48) ⊂ Σ given by H G ( z ) = { H G , H F } ( z ) = 0 . Let z ( · ) = ( q ( · ) , p ( · )) be a referencecurve on [0 , t f ] . We note t (cid:55)→ Φ ( t ) := H G ( z ( t )) the switching function, which codesthe switching times. Φ ( t ) = { H G , H F } ( z ( t )) , (1.13)¨ Φ ( t ) = {{ H G , H F } , H F } ( z ( t )) + u ( t ) {{ H G , H F } , H G } ( z ( t )) . (1.14)From this, we derive. Lemma 3.
Assume that t is an ordinary switching time, that is Φ ( t ) = 0 and ˙ Φ ( t ) (cid:54) = 0 . Then near z ( t ) every extremal projects onto σ + σ − if ˙ Φ ( t ) < and σ − σ + if ˙ Φ ( t ) > . The situation is more complex for higher contact with Σ . Definition 7.
The fold case is characterized by Φ ( t ) = ˙ Φ ( t ) = 0 and ¨ Φ ( t ) (cid:54) = 0 (replacing u by ± in (1.14) ). Hence z ( t ) ∈ Σ (cid:48) . Assume that locally Σ (cid:48) is a regularsurface of codimension two. We have three cases: • Parabolic case: ¨ Φ + ( t ) ¨ Φ − ( t ) > . • Hyperbolic case: ¨ Φ + ( t ) > and ¨ Φ − ( t ) < . • Elliptic case: ¨ Φ + ( t ) < and ¨ Φ − ( t ) > . Denote by u s ( · ) the singular control defined by (1.14) as ¨ Φ ( t ) = 0. In the hy-perbolic case, through z ( t ), assuming {{ H G , H F } , H G } ( z ( t )) (cid:54) = 0, there exists asingular arc, which is strictly admissible that is | u s ( t ) | <
1. This arc is hyperbolic if {{ H G , H F } , H G } ( z ( t )) > <
0. In the parabolic case, it can be absent or not admissible that is | u s ( t ) | > σ + σ − σ s Fig. 1.1.
Fold case in the elliptic 2D–situation and the antiturnpike phenomenon.
It is based on [7] and presented here for n = 3. It uses proposition 4. Under ourassumptions ( H ) − ( H ), and if moreover the singular control is strictly admissible on [0 , t f ], that is | u s ( t ) | <
1, the reference singular arc can be immersed in a C –domain up to the first conjugate time t c and where the time minimal policy is ofthe form σ ± σ s σ ± , see Fig. 1.2 for the interpretation of the first conjugate time t c for the problem in the q –space with q = ( c, v ), ˙ v = u . Geometric Techniques to Optimize Chemical Reactions 11 q (0) σ s σ + σ − s i n g u l a r a r c s σ s t c c (0) Fig. 1.2. (left)
Time minimal policy in the q –space. (right) Conjugate point in the c –space. More complicated and challenging situation is to analyze situations, where the strictLegendre–Clebsch condition is not satisfied that is {{ H G , H F } , H G } ( z ( t )) vanishesat some times. Due to the complexity and in relation with our application, we shallassume that n = 3.In this case, one has D = 0 since {{ H G , H F } , H G } = p · [[ G, F ] , G ] and recallthat D = det( G, [ G, F ] , [[ G, F ] , G ]). The singular flow is defined by the vector field X s ( q ) = F ( q ) − D (cid:48) ( q ) D ( q ) G ( q ). Using a time reparameterization, it is related to theone–dimensional foliation F s : D ( q ) F ( q ) − D (cid:48) ( q ) G ( q ). Complexity of the analysis isdue to non isolated singularities contained in the set D (cid:48) ( q ) = D ( q ) = 0.In relation with our optimal problem with terminal manifold N of codimensionone given by x = d , our singularity resolution will be concerned by analyzing arcsinitiating from the set S : n · [ G, F ] ( q ) = 0, where n = (1 , ,
0) is the normalvector to N . Singularities can be roughly classified into two types: local in relationwith singularities of S and propagated along the singular flow, and Lagrangiansingularities in relation with the concept of conjugate–focal points. This will bedeveloped in the next section. The system is written ˙ q = F ( q ) + u G ( q ), | u | ≤
1, with q = ( c, v ), G = ∂∂v . Theterminal manifold N is given by c = d . The problem is to determine small time synthesis near a given point q ∈ N , which can be identified to 0. More precisely onewants to classify the syntheses under generic assumptions near the terminal manifoldin relation with the Lie algebraic structure of { F, G } at q . Our approach developedin our series of articles [10, 9, 8, 21] is to use the construction of semi–normalform for the action of the pseudo–group G of local diffeomorphisms and feedbacktransformation u → − u (so that σ + and σ − can be exchanged). Additionally recallthat the pseudo–group G f formed by local diffeomorphims and feedback actions ofthe form u = α ( q ) + β ( q ) v leaves the singular flow invariant [4]. These groups act2 Bernard Bonnard and J´er´emy Rouoton the jet space of F at zero, G being identified to ∂∂v and the terminal manifold N to c = d . Note that the problem is flat that is G is tangent to N . For the actionof the pseudo–group G on the jet spaces of F we refer to [22] (we shall work in the C k category, where k ≥ F, G, N ). One has N : c = d and the initial state is such that c (0) < d . Denote in general N ⊥ = { ( q, p ); p · v =0 , ∀ v ∈ T q N } and let n be the outward normal to N so that n = (1 , , . . . ,
0) and N ⊥ = ( q, n ( q )). Let z = ( q, p ) be a BC–extremal on [ t f , t f < 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 ifthe optimal curves are tangent and second kind if they are transverse. The splittinglocus L is the set of points, where the optimal control is not unique and the cutlocus C is the closure of the set of points, where the optimal trajectories loses itsoptimality, see [3, 13] for the introduction of those concepts in the frame of semi–analytic 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 q ∈ R n , we denoteby j k F ( q ) , j k G ( q ) , j k f ( q ) the respective k –jets, that is the Taylor expansionsat order k . We say that ( F, G, f ) has at q a singularity of codimension i if( j k F ( q ) , j k G ( q ) , j k f ( q )) ∈ Σ i , where Σ i is a semi–algebraic submanifold of codi-mension i in the jets space.Taking q ∈ N with a singularity of codimension i , an unfolding is a C change ofcoordinates near q 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. We shall present the main steps to compute the time minimal syntheses, restrictingour study to the 3-dimensional case, but it can be clearly extended to the n –spacewith the concept of codimension [21].The system is written: ˙ q = F ( q ) + u G ( q ), | u | ≤ q = ( x, y, z ) bethe coordinates, G being identified to ∂∂z and G being tangent to N , which can beidentified to x = 0 and n = (1 , ,
0) is the outward normal to N . Recall that (generic)singular trajectories are given by ˙ q = F ( q ) + u s ( q ) G ( q ), u s ( q ) = − D (cid:48) ( q ) /D ( q ) andthey 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 ( q ) and [ G, F ] ( q ) are independent andnot tangent to N and the synthesis follows from Lemma 3. It is given by σ + if n · [ G, F ] ( q ) < σ − if n · [ G, F ]( q ) > virtual switching point. Generic hyperbolic singular case and the concept of focal points.
Nowwe present the basis of our analysis, the so–called hyperbolic situation. In this case,based on [9], a semi–algebraic normal form is constructed to obtain the correspondinglocal syntheses and it is extended (in the jet space) along the reference singular arc
Geometric Techniques to Optimize Chemical Reactions 13in order to obtain the concept of focal point based on the notions of turnpike andconjugate points presented in Section 3.
Semi–normal form: one has q = 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 passingthrough 0 and transverse to G . • D (0) and D (cid:48)(cid:48) (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 , δ tothe axis ( Oy ) can be identified to t (cid:55)→ ( t, , Ox )–axis. Notethat N is identified to x = 0.The semi–normal form is constructed in a tubular neighborhood of the referencesingular curve σ s and the vector field F is developed in the jet space with respect to( y, z ) since σ s : t (cid:55)→ ( t, ,
0) is the x–coordinate. One has the following semi-normalform ˙ x = 1 + a ( x ) z + 2 b ( x ) yz + c ( x ) y + R ˙ y = d ( x ) y + e (0) z + R ˙ z = ( u − u s | σ ( x )) + f ( x ) y + g (0) z + R , where a ( x ) (cid:54) = 0 , e (0) (cid:54) = 0 and R (resp. R , R ) are terms of order ≥ ≥ y, z ) and u s is the singular control.Furthermore we make the following assumption: • The reference arc is hyperbolic on [ t f ,
0] so that a ( x ) < q is of the form σ + σ s σ − , seeFig.1.3 and is described by the unfolding˙ x = 1 − a (0) z ˙ z = u − u s (0)˙ y = 0in a small neighborhood of q = 0. 0 σ + σ − σ s x = 0 x < Fig. 1.3.
Local syntheses and the semi-bridge phenomenon.
We use [9] and the
Mathematica ’s program given in Appendix 6.1. One fixes a point q ∈ S and we make the explicit computations of the switching and cut loci near q using the jet expansion of F at q , G and N being normalized respectively to ∂∂v and x = 0. The computations are in the semi-algebraic category, easily implementableusing symbolic algorithms and the cut locus and switching loci are stratifiable. Switching locus.
We denote by K the set of ordinary switching points for BC-extremals, K + being a switching σ − σ + and K − being a switching σ + σ − , while W, W + , W − corresponds respectively to switching points of optimal extremals. Moreprecisely, near N the stratification of W is W = W ∪ W where W is the first kindstratum composed by the hyperbolic singular arcs and W = W s ∪ W + ∪ W − is thesecond kind strata, where W ε , ε ∈ {− , } is composed by the ordinary switchingpoints of the policy σ − ε σ ε and W s is the first switching point of the Bang-Bang-Singular’s policy.From [1] we have two situations describe by Fig.1.4. The reflecting case corres- Kσ + σ − K σ − σ + Fig. 1.4. (left)
Crossing. (right)
Reflecting.ponds to a non optimal situation and has to be rejected to determine W .Using the theory and computations of [9] only the stratification of W localizednear q can be computed. Singular locus.
Near q , one determines the singular locus Γ s restricting to ad-missible singular trajectories, which are optimal. In other words, the test is : hyper-bolicity and strict admissibility | u s | < Cut locus.
The cut locus C is the set of points where optimality is lost. One partof the cut locus is formed by the splitting locus L , where two minimizers intersects. Itcan be stratified into strata corresponding between intersections σ + , σ − intersectionsbetween σ + σ − and σ − . . . Again such intersections are described in [9]. 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 switchingpoint. Denoting Φ ( t ) the switching function on [ t f , Φ (0) = n · G ( q ) one has:if ˙ Φ (0) := n · [ G, F ]( q ) < >
0) the terminating arc is σ − (resp. σ + ). Codimension one case.
We have different cases corresponding to the foldcase, since in the parabolic case one must distinguish between a parabolic pointcorresponding either to a non admissible hyperbolic or elliptic arc. The cases arerepresented on Fig.1.5-1.6.The syntheses can be represented by foliations by 2d-planes. Note that the roleof σ + and σ − can be interchanged since the semi-normal forms are computed usingthe transformation u (cid:55)→ − u . Geometric Techniques to Optimize Chemical Reactions 15 σ − σ − σ − σs N ∩ y = constant σ + σ + σ + Hyperbolic
C Nσ − σ + Elliptic
Fig. 1.5.
Fold case.
C Nσ − σ + W − Nσ + σ − σ + Non admissible elliptic case W + Nσ − σ + Non admissible hyperbolic case
Fig. 1.6.
Parabolic point with a non admissible singular arc.The C -unfoldings are: • Hyperbolic and Parabolic:˙ x = 1 + a y ˙ y = u − u s (0) , | u s (0) | < , and a < a > • Parabolic: ˙ x = 1 + a y ˙ y = u − u s (0) , | u s (0) | > , and a (cid:54) = 0. 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 := δ (cid:48) (0) belongs to span { G (0) , [ G, F ](0) } . Let λ , λ be two scalars such that W = λ G (0) + λ [ F, G ](0) and by assumption λ (cid:54) = 0. Denoting by X s ( q ) := F ( q ) − D (cid:48) ( q ) D ( q ) G ( q ), we introduce the variational equation˙ δq ( t ) = ∂X s ∂q ( σ s ( t )) δq ( t )and let W ( · ) be a solution on [ t f ,
0] such that W (0) = W . Definition 8.
Let t f be the first time in [ t f , such that: det( W ( t ) , G ( σ s ( t )) , F ( σ s ( t ))) = 0 . Then t f is called the first focal point along σ s . As for the fixed end point problem we have the following:6 Bernard Bonnard and J´er´emy Rouot
Proposition 5.
Provided t ∈ ] t f , , then in a tubular neighborhood of σ s the set S = { exp( tX s ( σ )) } 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 singularvalue decomposition, an equivalent computation is to determine t = t f 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 cons-tructed in [9]. A model is˙ x = 1 + a z + α xy + α yz + α xz ˙ y = b z ˙ z = c z − u s (0) − u sx x − u sy y, (1.15)where u s is the singular control, u s (0) ∈ { , } , u sx = ∂u s ∂x (0), u sy = ∂u s ∂y (0) > a, b, c (cid:54) = 0. The switching function Φ ( t ) := p ( t ) is developed at order 3 andfactorized as tP ( t ), where P is polynomial of order two with two roots t , t , whichdetermine the switching points.If a >
0, we are in the elliptic situation and if a < σ s : t (cid:55)→ ( t, ,
0) being not admissible.One has | u s (0) | > u s (0) > Mathematica ’s pro-gram of Appendix 6.1, we use the following values: α = α = α = u sy = c = b = 1. The case a > u s (0) = 3. We have the following three situations
C Nσ − σ + K − N σ − σ + σ − σ − W − C σ + σ + I W − Nσ + σ − σ + Fig. 1.7.
Bifurcation of the cut ( a > , u s (0) = 3). The cut locus C splits into C ∪ C where C : σ − ∩ σ + and C : σ + σ − ∩ σ − .Note that the cut appears when the trajectories reflect on the switching surface,see Fig.1.8. Geometric Techniques to Optimize Chemical Reactions 17 y < Fig. 1.8.
A situation in the elliptic parabolic case. The trajectories of Γ − reflect on W − leading to a cut locus between W − and W + . The case a < u s (0) = 1. There are two generic cases described byFig.1.9-1.10. The two cases are discriminated by the existence or not of the singulararc in the transition. In the second case, the switching locus Σ has two strata : Σ = W + ∪ W s , where W + corresponds to optimal policies σ − σ + and W s to policies σ − σ + σ s , Σ is not C . The stratum W s ∪ Γ s is C but not C . σ − N σ + σs N σ − σ − σ − σ + σ + σ + σs σ − W + K + N σ − σ − σ + σ + σs y < y = 0 y > Fig. 1.9.
Saturating case 1: − a <
0, 1 = u sx >
0. Top figures are obtainedusing the
Mathematica ’s program given in the Appendix 6.1 and these situationsare sketched in the bottom figures.
In our previous work in the 90’s we restricted our study to the case, where the strictgeneralized Legendre–Clebsch condition is satisfied, see [9]. Now our program is toextend this analysis when this condition is not satisfied.8 Bernard Bonnard and J´er´emy Rouot W + N σs σ + σ + σ + σ − σ − σ − σ − σ − σ − W s σ − σ − σ − N σ + W + σ − σ − σ − σ + W + N Fig. 1.10.
Saturating case 2: − a <
0, 1 = u sx <
0. Top figure is obtainedusing the
Mathematica ’s program given in the Appendix 6.1 and the correspondingsituations are sketched in the bottom figures.
The semi–bridge phenomenon as transition between two saturations
A bridge in optimal control was introduced in [6] as a policy of the form Bang-Singular-Bang-Singular, where the second bang is related as a saturation due to theviolation of the strict Legendre-Clebsch condition. In our context a semi–bridge isinterpreted as a path to a bridge due to saturation between optimal singular arcs.
A tutorial model for the semi–bridge.
Consider the case˙ x = 1 + a y − c yz + c z ˙ y = z ˙ z = u, | u | ≤ , (1.16)where a, c (cid:54) = 0 are parameters. Lie brackets computations.
We have F = (1 + ay − cyz + cz ) ∂∂x + z ∂∂y , G = ∂∂z , [ G, F ] = 3 c ( y − z ) ∂∂x − ∂∂y , [[ G, F ] , G ] = − cz ∂∂x , [[ G, F ] , F ] = a ∂∂x , D ( q ) = det( G, [ G, F ] , [[ G, F ] , G ]) = − cz,D (cid:48) ( q ) = det( G, [ G, F ] , [[ G, F ] , F ]) = a, D (cid:48)(cid:48) ( q ) = det( G, [ G, F ] , F ) = ay − cz + 1 . Stratification of N . We have the following properties • N : x = 0, G = ∂∂z • S : n · [ G, F ] = 0 ∩ { x = 0 } is the parabola: y = z . • n · [[ G, F ] , G ] (0) = 0, n · [[ G, F ] , F ] (0) (cid:54) = 0. Geometric Techniques to Optimize Chemical Reactions 19 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.Observe that E : n · F ( q ) = 0 is the defined by y ( a − cz ) + cz + 1 = 0.One will localize our study in a neighborhood V of 0 such that: E ∩ V = ∅ . Werepresent such a situation on Fig.1.11 One has u s = a/ (6 cz ) and the points u sat aregiven by | u s | = 1.We have two saturating points, one in the hyperbolic domain and one in theelliptic domain. Near the fold, the two saturation phenomenons are glued togetherusing the curvature of S , the normal to N being given by n = (1 , , q follows S . Note that since the control is blowing up at the fold, we encounter thecase a > , u s (0) = 3 (see Fig.1.7) and the bifurcation phenomenon of C , whichsplits into the case σ + , σ − and the case σ + σ − and σ − intersecting minimizers. Thetwo strata will be denoted C and C . The model detects the two strata.The stratification of the terminal manifold near the semi–bridge and the localtime minimal synthesis are described in Fig.1.12. u sat u sat y vu = − u = +1 S ( iv )( i ) ( ii )( vi )( v ) ( i ) : Fig.1 . a < , | u s (0) | = 1( ii ) : Fig.1 . a < , | u s (0) | > iii ) ( iii ) : Fig.1 . a > , | u s (0) | > iv ) : Fig.1 . center ) : a > , u s (0) = 3( v ) : Fig.1 . left ) : a < , | u s (0) | < vi ) : Fig.1 . right ) : a > , | u s (0) | < Fig. 1.12.
Stratification of the target x = 0 for a semi–bridge model and theassociated local syntheses.0 Bernard Bonnard and J´er´emy Rouot Local syntheses by symbolic computations.
We use the
Mathematica ’s program in Appendix 6.1 to derive the local synthesesfor the semi–bridge model.We compute the time minimal synthesis in a neighborhood U of q ∈ N where N : x = 0 is the target. Since the problem is flat, q is either an ordinary switchingpoint, a fold point, or a codimension two case.Denote n the unit normal of N at q such that n belongs to the half-spacecontaining the set { X + u Y, | u | ≤ } , i.e. n = (1 , , X ( q ) and [ Y, X ]( q ) are not tangent to N . Then q is an ordinary pointand the optimal control is given by u ( q ) = − sign p · [ Y, X ]( q ).We take q = (0 , s , s ) and we look at the behavior of BC-extremals reaching N near q . 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 someorder. • Near a parabolic point on an hyperbolic singular arc.BC-extremals σ − switch while BC-extremals σ + don’t. These computations arerepresented in Fig. 1.13 for q = (0 , ( − . , − . • Near a parabolic point on an elliptic singular arc.BC-extremals σ ± reflect on the switching surface W ± , hence a cut locus ap-pears between the switching surfaces. We represent in Fig. 1.14 with q =(0 , . , . , . > z sat = a/ (6 c ) 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 q = (0 , z sat , − z sat ) where z sat = a c . The synthesis is represented inFig.1.15 and it corresponds to the synthesis described earlier by Fig.1.9. Fig. 1.13.
Time minimal synthesis near an hyperbolic fold point for q =(0 , ( − . , − . − a = c = 1. The tangent space of the switching surface at q is represented. Geometric Techniques to Optimize Chemical Reactions 21 Fig. 1.14.
Time minimal synthesis near an elliptic fold point for q = (0 , . , . a = 1 , c = 5. Fig. 1.15.
Time minimal synthesis near a saturating hyperbolic fold point q =(0 , z sat , − z sat ) for the tutorial model (1.16) with a = c = 1. Trajectories σ − switch,while σ + don’t. The surface W s is the first switching point of the sequence σ + σ − σ s . These computations confirm the expected results and validate the symbolic pro-gram to investigate the McKeithan model.Remark 1.
We can investigate other singularities by the same method, for instance: • Cusp singularity y = z . • Singularity y = z . • Case D = D (cid:48) = 0 is q = { } . Global analysis. u = ε ∈ {− , } . The dynamics (1.16) is integrable and we have p ( t ) = 12 t (cid:0) − ct + ( a − cεs ) t − cs + 6 cw (cid:1) . Denoting by Γ ε ( q ) , ε ∈ {− , } the surface composed by the trajectories σ ε passing through q , then Γ ε (0 , w , s ) is parameterized by x = t (cid:0) ( a − cs ) + cs + 1 (cid:1) + 12 t (cid:0) s ( a − cs ) + 3 cε ( s − w ) (cid:1) + o ( t ) y = w + s t + εt o ( t ) z = s + εt + o ( t ) . Solving p = 0 with respect to w , the switching surface K ε for BC-extremals ofcontrole u = ε on N is parameterized by: x = 16 t (cid:18) − a c + 6 aεs + 6 as − cεs − cs (cid:19) + 16 t (cid:0) as − cs + 6 (cid:1) ,y = t (cid:16) ( ε + 1) s − a c (cid:17) + 16 (3 ε + 2) t + s z = s + εt + o ( t ) . The switching condition of a trajectory on W ε is given by the sign ofdet (cid:18) ∂K ε ∂t | t =0 , ∂K ε ∂t | t =0 , ∂Γ ε (0 , s , s ) ∂t | t =0 (cid:19) , and this determinant is equal to (cid:0) as − cs + 1 (cid:1) (6 cεs − a ) / (6 c ). In the case a = c = 1, we get that σ + (resp. σ − ) switches on W + (resp. W − ) for z ∈ ] z sat ,
1[ (resp. z ∈ ] − z sat , Fig. 1.16.
The trajectories σ + starting from q = (0 , w , s ) such that w < s switch on W + for z ∈ ] z sat ,
1[ ( a = c = 1).We deduce the time minimal synthesis from Fig.1.18-1.19 for all point q along S .More precisely, if z < − z sat the optimal policy is Bang-Singular-Bang; if z = − z sat ,the synthesis is given by Fig.1.15. If − z sat < z < z sat , we have first a σ + σ − ’s policy,then we encounter the case u s ( z ) = 3 (see Fig.1.7), then we have a cut for z > z sat . Geometric Techniques to Optimize Chemical Reactions 23 Fig. 1.17.
The trajectories σ − starting from q = (0 , w , s ) such that w > s switch on W − for z ∈ ] − z sat ,
1[ ( a = c = 1).Cut locus for z > z sat . u s ( z ) = 3. Fig. 1.18. (left) σ − reflect on W − for z > z sat . (right) At z = a c < z sat such that u s ( z ) = 3, σ − cross W − and the cut locus disappears (see Fig.1.7). For − z sat < z For z > z sat σ + reflect on W + , hence we have a cut locus. For z < − z sat ,the optimal policy is Bang-Singular.4 Bernard Bonnard and J´er´emy Rouot Construction of a semi–normal form and the need of integrabilityassumption. Notations and assumption: We normalize in a small neighborhood of N . We assumethat F is transverse to N . Let q = ( x, y, z ), q = (0 , , G = ∂∂z and N := { x = 0 } .We denote by σ the reference trajectory of F identified to t → ( t, , 0) so that F | σ = ∂∂x .Assuming x small enough one can write the system as:˙ x = 1 + f ( y, z ) + R ˙ y = g ( y, z ) + R ˙ z = h ( y, z ) + R , where R i ( x, y, z ) = o ( | x | ) and the corresponding model is˙ x = 1 + f ( y, z )˙ y = g ( y, z )˙ z = h ( y, z ) + u. Singular flow on the model: Furthermore, we assume ∂g∂z (cid:54) = 0. Since the singular flowis feedback invariant, setting z = g ( y, z ), z → z and using a proper feedback thesystem takes the form ˙ x = 1 + f ( y, z )˙ y = z ˙ z = u. Computing, we get:[ G, F ] = − f z ∂∂x + ∂∂y [[ G, F ] , G ] = − f z ∂∂x , [[ G, F ] , F ] = − f yz ∂∂x + f y ∂∂y . (1.17)Hence D = − f z , D (cid:48) = − f y f z + f yz and u s = − D (cid:48) D = f y f z − f yz f z . Therefore, we have d y d z = zu s and assuming the separability condition d y d z = h ( y ) h ( z )one can integrate the singular trajectories such that q (0) ∈ S as follows, paramete-rizing by z instead of t using the relation (cid:90) yy (0) h ( y ) d y = (cid:90) zz (0) h ( z ) d z. Geometric Techniques to Optimize Chemical Reactions 25Assuming that y (0) can be expressed as y (0) = p ( z (0)) and moreover that theleft–hand side is invertible, one gets y = p ( z, z (0)). This defines the singular leaf F through q (0) ∈ S using the equationd x d z = 1 + f ( y, z ) z . Using x (0) = 0, one gets: x ( z, z (0)) = (cid:90) zz (0) (1 + f ( p ( z, z (0)) , z )) d z. Remark 2. This requires an integrability condition to obtain the leaf F . If it is notsatisfied, 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 theseparable case and using our previous algorithm one gets Proposition 6. 1. Integrating the singular flow, we get: x ( t, z (0)) = t + ε (2 a +3 ck ) ( kt + z (0) ) / − c ( kt + z (0) ) t k y ( t, z (0)) = ε k (cid:0) kt + z (0) (cid:1) / z ( t, z (0)) = ε (cid:0) kt + z (0) (cid:1) / , where ε = ± , k = 12 c/a .2. Parameterizing by z , the integral singular leaf F with q (0) ∈ S is given y ( z, z (0)) = 2 c ( z − z (0) ) /a + z (0) x ( z, z (0)) = g ( z ) − g ( z (0)) , where g ( z ) = 3 cz a (cid:16) a (2 c +1) z +30 cz (0) ( a − cz (0)) z +30 c z +5 a ( az (0) − cz (0) + 1) (cid:17) . We represent on Fig.1.20 the stratification of the singular leaf F with a = c = 1using hyperbolic and elliptic cases and the admissibility conditions | u s | ≤ Recall that the network is T + M A B k k k k and the system ˙ q = F ( q )+ u G ( q ),using the coordinates q = ( x, y, v ), x = [ A ], y = [ B ], v = k and reduced to thestoichiometric class T + M + A = δ and T + M + B = δ takes the form˙ x = − β xv α − β xv α − δ v ( x + y ) + δ v + v ( x + y ) ˙ y = β xv α − β yv α with 0 ≤ x ≤ δ , 0 ≤ y ≤ δ , δ = δ + δ , δ = δ δ , k = β v α , k = β v α , k = β v α .6 Bernard Bonnard and J´er´emy Rouot Fig. 1.20. Stratification of the singular set Γ s for the system (1.16) with parameters a = c = 1. We have: • F = (cid:0) − β xv α − β xv α − δ v ( x + y )+ δ v + v ( x + y ) (cid:1) ∂∂x + (cid:0) β xv α − β yv α (cid:1) ∂∂y , • G = ∂∂v , • [ G, F ] = (cid:0) x ( α β v α − + α β v α − + δ ) + δ y − δ − x − xy − y (cid:1) ∂∂x + (cid:0) α β yv α − − α β xv α − (cid:1) ∂∂y , • [[ G, F ] , G ] = (cid:0) x ( α β ( α − v α − + α β ( α − v α − ) (cid:1) ∂∂x + (cid:0) yα β ( α − v α − − xα β ( α − v α − (cid:1) ∂∂y , • [[ G, F ] F ] = (cid:0) − xv α β δ ( α − − yv α β δ ( α − 1) + v α ( α β δ + xy (2 α β − β ) − β δ ) + x ( α β − β ) v α + y ( α β − β ) v α + yv α ( β δ − α β δ ) + v α ( α β δ − β δ )+ x ( β − α β ) v α + y ( α β − β ) v α + yv α ( α β δ − β δ )+ xy (2 β − α β ) v α + y (2 β − α β ) v α (cid:1) ∂∂x + (cid:0) x ( v α + α − ( α β β − α β β )+ v α + α − ( α β β − α β β )) + xv α ( α β δ − β δ ) + yv α ( α β δ − β δ ) + v α ( − α β δ + xy (2 β − α β )+ β δ )+ x ( β − α β ) v α + y ( β − α β ) v α (cid:1) ∂∂y . One has: • D ( q ) = (( α − α β yv α − − ( α − α β xv α − )( α β xv α + α β xv α + δ v ( x + y ) − δ v − vx − vxy − vy ) + x ( α β v α − α β v α + ( α − α β v α )( α β xv α − − α β yv α − ), • D (cid:48) ( q ) = β v α − ( α β xv α + α β xv α + δ v ( x + y ) − δ v − vx − vxy − vy )( α β xv α − α β xv α + ( α − δ v ( x + y ) + δ ( v − α v ) − α vx − α vxy − α vy − α β xv α + α β xv α + vx + 2 vxy + vy ) + ( α β xv α − − α β yv α − )(( α − β v α ( δ − ( x + y )( δ − x − y )) + ( α − β v α ( y ( y − δ ) + δ − x ) + ( α − β yv α ( δ − x + y ))), Geometric Techniques to Optimize Chemical Reactions 27 • D (cid:48)(cid:48) ( q ) = ( β xv α − − β yv α − )( α β xv α + α β xv α + δ v ( x + y ) − δ v − vx − vxy − vy ) − ( α β xv α − − α β yv α − )( β xv α + β xv α + δ v ( x + y ) − δ v − vx − vxy − vy ),and the singular control is given by: u s = − D (cid:48) ( q ) /D ( q ). 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 , , S : α β dv α − + α β dv α − + dδ − δ − d + y ( δ − d ) − y = 0 . (1.18)Denoting by ∆ the discriminant of the polynomial function y (cid:55)→ n · [ G, F ]( q ) ∩ x = d ,a singularity can occur for ∆ = 0.One has Lemma 4. Assume α i , β i , δ i > , i=1,2 and d, v > . Then we have ∆ = ( δ − δ ) +4 d ( α β v α − + α β v α − ) > so that there is no ramification and S contains atmost 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 α − α = − ( α − α β ( α − α β . (1.19) Exceptional locus It is given by E : n · F ( q ) = 0 and x = d . Computing, one gets: E : − β dv α − β dv α + d v − dδ v + y (2 dv − δ v ) + δ v + vy = 0 . (1.20)The discriminant of the polynomial n · X ( q ) in y is ∆ = v (4 d ( β v α + β v α ) + v ( δ − δ ) ) > E contains at most two real positive branches.Fig.1.21 gives a picture of stratification of x = d for the McKeithan system witha 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 admissiblepoints. We are in the flat case and a point of N is either an ordinary point or afold point. For an ordinary point, the final optimal control is regular and is equalto − sign Φ (0). At a fold point on S , the final optimal control may be singular (seeDefinition 7).The optimal syntheses near the fold can be described using the techniques ofthe tutorial model (1.16) and symbolic computations.8 Bernard Bonnard and J´er´emy Rouot E S D · D (cid:48)(cid:48) < D · D (cid:48)(cid:48) > u sat u sat equilibrium n · F > u sat yv Fig. 1.21. Stratification of the surface x = d for the McKeithan reaction. Dottedline: elliptic, red line: hyperbolic. In this article, we have presented the general tools to analyze a Mayer problem tooptimize the yield of chemical networks. Using the McKeithan network, we havecompleted the analysis of [8] for the simple scheme A → B → C .The main point is to extend a geometric techniques from [9] to consider the casewhen the strict Legendre-Clebsch is not satisfied. The starting point is to analyzethe semi-bridge phenomenon.Note that with the same techniques we can analyze the situations up to codi-mension 2 describe in details in [21] (in particular near the exceptional locus E ).Furthermore, it can be extended to an important problem of codimension 3 occur-ring for reversible chemical networks and the existence of equilibria, see Fig. 1.21. (* Strata.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 aneighborhood of N.*)(* Lie derivative of w wrt v *)LieDerive[v_List,w_List,var_List]:=Transpose[D[w, Geometric Techniques to Optimize Chemical Reactions 29 Length[v] == Length[w] == Length[var](* Poisson bracket {f,g} *)PoissonBracket[f_, g_,q_List,p_List] /; Length[q] == Length[p] := \Total @ Flatten @ Fold[List,0,MapThread[D[f, 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 ..........*)ord = 3;us = - DDp / DD;q0p0 = {0,w0,0,1,0,0};qpt = IntSmallTime[X,Y,us,q0p0,ord];ps = {t0val,w0val,s0val};gms = mTaylor[ Geometric Techniques to Optimize Chemical Reactions 31 qpt = mTaylor[ c12sols = mTaylor[ References 1. A. Agrachev , R.V. Gamkrelidze, Symplectic methods for optimization and con-trol , in Geometry of Feedback and Optimal Control, Monogr. Textbooks PureAppl. Math. , B. Jakubczyk and W. Respondek, eds., Marcel Dekker, NewYork, (1998) pp. 19–772. D. Anderson, Global asymptotic stability for a class of nonlinear chemical equa-tions , SIAM J. Appl. Math. no.68 (2008) 1464–1476.3. V.G. Boltyanskii, Sufficient conditions for optimality and the justification of thedynamic programming method , SIAM J. Control Optim., (1966) 326–361.4. B. Bonnard, Feedback equivalence for nonlinear systems and the time optimalcontrol problem , SIAM J. Control Optim. , no.6, (1991) 1300–1321.5. B. Bonnard, M. Chyba, Singular trajectories and their role in control the-ory , Math´ematiques & Applications Springer-Verlag Berlin Heidelberg, 2003xvi+357.6. B. Bonnard, O. Cots, J. Rouot, T. Verron, Time minimal saturation of a pair ofspins and application in magnetic resonance imaging , Math. Control Relat. Fields,AIMS (2019).7. B. Bonnard, I. Kupka, Th´eorie des singularit´es de l’application entr´ee/sortie etoptimalit´e des trajectoires singuli`eres dans le probl`eme du temps minimal , ForumMath. , no.2 111–159 (1993).8. B. Bonnard, G. Launay, Time minimal control of batch reactors , ESAIM: Control,Optimisation and Calculus of Variations (1998) 407–467.9. B. Bonnard, G. Launay, M. Pelletier, Classification g´en´erique de synth`eses tempsminimales avec cible de codimension un et applications , Annales de l’I.H.P. Anal-yse non lin´eaire no.1 (1997) 55–102.10. B. Bonnard, M. Pelletier, Time minimal synthesis for planar systems in theneighborhood of a terminal manifold of codimension one , J. of Mathematical Sys-tems, Estimation and Control no.3 (1995) 379–381.11. B. Bonnard, D. Sugny, Optimal Control with Applications in Space and Quan-tum Dynamics , AIMS Series on Applied Mathematics , 2012 xvii+298.12. U. Boscain, B. Piccoli, Optimal Syntheses for Control Systems on 2-D Manifolds ,Math´ematiques et Applications, Springer–Verlag Berlin Heidelberg, 2004 xiv–262.13. P. Brunovsky, Existence of Regular Syntheses for General Problems , J. Diff.Equations (1980) 317–343.14. M. Chavez, Observer design for a class of nonlinear systems, with applica-tions to chemical and biological networks , Ph.D. thesis, Rutgers University, NewBrunswick, NJ, 2003.15. M. Chavez, E. Sontag, State-estimators for chemical reaction networks ofFeinberg–Horn–Jackson zero deficiency type , Euro. J. Cont. (2002) 343–359.16. I. Ekeland, Discontinuit´e des champs hamiltoniens et existence des solutionsoptimales en calcul des variations , Inst. Hautes ´Etudes Sci. Publ. Math. , (1977)5–32. Geometric Techniques to Optimize Chemical Reactions 3317. M. Feinberg, On chemical kinetics of a certain class , Rational Mech. Anal., no.1, (1972).18. F. Horn, R. Jackson, General mass action kinetics , Arch. Ration. Mech. Anal., no.1, (1972) 81–116.19. A. J. Krener, The high order maximal principle and its application to singularextremals . SIAM Journal on Control and Optimization, no.2, (1977) 256–293.20. I. Kupka, Geometric theory of extremals in optimal control problems. I. Thefold and Maxwell case , Trans. Amer. Math. Soc. no.1, (1987) 225–243.21. G. Launay, M. Pelletier, The generic local structure of time-optimal synthesiswith a target of codimension one in dimension greater than two. J. Dynam. ControlSystems no.2 (1997) 165–203.22. J. Martinet, Singularities of smooth functions and maps , London MathematicalSociety Lecture Note Series, , Cambridge University Press, Cambridge-NewYork, 1982, xiv+256.23. T.W. McKeithan Kinetic proofreading in T-cell receptor signal transduction Proc. Natl. Acad. Sci. USA, (1995) 5042–5046.24. L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, E.F. Mishchenko, Themathematical theory of optimal processes , Oxford, Pergamon Press, 1964.25. H. Sch¨attler, U. Ledzewicz, Geometric optimal control. Theory, methods andexamples , Interdisciplinary Applied Mathematics, Springer-Verlag New York,2012 xx+640.26. H. Sch¨attler, On the local structure of time-optimal bang-bang trajectories in R , SIAM J.Control Optim., (1988) 186–204.27. E.D. Sontag Structure and stability of certain chemical networks and applica-tions to the kinetic proofreading model of T-cell receptor signal transduction IEEETrans. Automat. Contr. no.7 (2001) 1028–1047.28. H.J. Sussmann, Synthesis, presynthesis, sufficient conditions for optimality andsubanalytic sets , Nonlinear Controllability and Optimal Control, H.J. SussmannEd., Marcel Dekker, New York, (1990) 1–20.29. H.J. Sussmann, The Structure of Time–Optimal Trajectories for Single–InputSystems in the Plane: the C ∞ Nonsingular Case , SIAM J. Control and Opt. (1987) 433–465.30. H.J. Sussmann, The Structure of Time–Optimal Trajectories for Single–InputSystems in the Plane: the General Real Analytic case , SIAM J. Control and Opt. (1987) 868–904.31. J. T´oth, A. L. Nagy, D. Papp, Reaction kinetics: exercises, programs andtheorems , Springer-Verlag New York, 2018 xxiv+469.32. E. Tr´elat, E. Zuazua, The turnpike property in finite-dimensional nonlinearoptimal control , J. Differential Equations,258