Giulia Cavagnari
email: giulia.cavagnari@unitn.it

Antonio Marigonda
email: antonio.marigonda@univr.it

Giandomenico Orlandi
email: giandomenico.orlandi@univr.it

Hamilton-Jacobi-Bellman Equation for a Time-optimal Control Problem in the Space of Probability Measures

Keywords: optimal transport, differential inclusions, time optimal control

In this paper we formulate a time-optimal control problem in the space of probability measures endowed with the Wasserstein metric as a natural generalization of the correspondent classical problem in R d where the controlled dynamics is given by a differential inclusion. The main motivation is to model situations in which we have only a probabilistic knowledge of the initial state. In particular we prove first a Dynamic Programming Principle and then we give an Hamilton-Jacobi-Bellman equation in the space of probability measures which is solved by a generalization of the minimum time function in a suitable viscosity sense.

Introduction

The controlled dynamics of a classical time-optimal control problem in finitedimension can be presented by mean of a differential inclusion as follows: ẋ(t) ∈ F (x(t)), for a.e. t > 0,

x(0) = x 0 ∈ R d , (1)
where F is a set-valued map from R d to R d . The problem in this setting is to The main motivation of this work is to model situations in which the knowledge of the starting position x 0 is only probabilistic (for example in the case of affection by noise) and this can happen even if the evolution of the system is deterministic.

We thus consider as state space the space of Borel probability measures with finite p-moment endowed with the p-Wasserstein metric W p (•, •), (P p (R d ), W p ). In [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] the reader can find a detailed treatment about Wasserstein distance.

Following this idea we choose to describe the initial state by a probability measure µ 0 ∈ P p (R d ) and for its evolution in time we take a time-depending probability measure on R d , µ := {µ t } t∈[0,T ] ⊆ P p (R d ), µ |t=0 = µ 0 . In order to preserve the total mass µ 0 (R d ) during the evolution, the process will be described by a (controlled) continuity equation

∂ t µ t + div(v t µ t ) = 0, 0 < t < T, µ |t=0 = µ 0 , (3)
where the time-depending Borel velocity field v t : R d → R d has to be chosen in the set of L 2 µt -selections of F in order to respect also the classical underlying control problem (1) which is the characteristic system of (3) in the smooth case.

It is well known that if v t (•) is sufficiently regular then the solution of the continuity equation is characterized by the push-forward of µ 0 through the unique solution of the characteristic system.

In Theorem 8.2.1 in [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] and Theorem 5.8 in [START_REF] Bernard | Young measures, superpositions and transport[END_REF], the so called Superposition Principle states that, if we conversely require much milder assumptions on v t , the solution µ t of the continuity equation can be characterized by the push-forward e t η, where e t : R d × Γ T → R d , (x, γ) → γ(t), Γ T := C 0 ([0, T ]; R d ) and η is a probability measure in the infinite-dimensional space R d × Γ T concentrated on those pairs (x, γ) ∈ R d × Γ T such that γ is an integral solution of the underlying characteristic system, i.e. of an ODE of the form γ(t) = v t (γ(t)), with γ(0) = x. We refer the reader to the survey [START_REF] Ambrosio | The flow associated to weakly differentiable vector fields: recent results and open problems, contribute in Nonlinear conservation laws and applications[END_REF] and the references therein for a deep analysis of this approach that is at the basis of the present work.

Pursuing the goal of facing control systems involving measures, we define a generalization of the target set S by duality. We consider an observer that is interested in measuring some quantities φ(•) ∈ Φ; the results of this measurements are the average of these quantities w.r.t. the state of the system. The elements of the generalized target set SΦ p are the states for which the results of all these measurements are below a fixed thershold.

Once defined the admissible trajectories in this framework, the definition of the generalized minimum time function follows in a straightforward way the classical one.

Since classical minimum time function can be characterized as unique viscosity solution of a Hamilton-Jacobi-Bellman equation, the problem to study a similar formulation for the generalized setting would be quite interesting. Several authors have treated a similar problem in the space of probability measures or in a general metric space, giving different definitions of sub-/super differentials and viscosity solutions (see e.g. [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF][START_REF] Ambrosio | On a class of first order Hamilton-Jacobi equations in metric spaces[END_REF][START_REF] Cardaliaguet | Deterministic differential games under probability knowledge of initial condition[END_REF][START_REF] Gangbo | Hamilton-Jacobi Equations in the Wasserstein Space[END_REF][START_REF] Gangbo | Optimal transport and large number of particles[END_REF]). For example, the theory presented in [START_REF] Gangbo | Optimal transport and large number of particles[END_REF] is quite complete: indeed there are proved also results on time-dependent problems, comparison principles granting uniqueness of the viscosity solutions under very reasonable assumptions.

However, when we consider as metric space the space P 2 (R d ), we notice that the class of equations that can be solved is quite small: the general structure of metric space of [START_REF] Gangbo | Optimal transport and large number of particles[END_REF] allows only to rely on the metric gradient, while P 2 (R d ) enjoys a much more richer structure in the tangent space (which is a subset of L 2 ).

Dealing with the definition of sub-/superdifferential given in [START_REF] Cardaliaguet | Deterministic differential games under probability knowledge of initial condition[END_REF], the major bond is that the "perturbed" measure is assumed to be of the form (Id R d + φ) µ in which a (rescaled) transport plan is used. It is well known that, by Brenier's Theorem, if µ L d in this way we can describe all the measures near to µ. However in general this is not true. Thus if the set of admissible trajectories contains curves whose points are not all a.c. w.r.t. Lebesgue measure (as in our case), the definition in [START_REF] Cardaliaguet | Deterministic differential games under probability knowledge of initial condition[END_REF] cannot be used.

In order to fully exploit the richer structure of the tangent space of P 2 (R d ), recalling that AC curves in P 2 (R d ) are characterized to be weak solutions of the continuity equation (Theorem 8.3.1 in [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]), we considered a different definition than the one presented in [START_REF] Cardaliaguet | Deterministic differential games under probability knowledge of initial condition[END_REF] using the Superposition Principle.

The paper is structured as follows: in Section 2 we give the definitions of the generalized objects together with the proof of a Dynamic Programming Principle in this setting. In Section 3 we focus on the main result of this work, namely we outline a Hamilton-Jacobi-Bellman equation in P 2 (R d ) and we solve it in a suitable viscosity sense by the generalized minimum time function, assuming some regularity on the velocity field. Finally, in Section 4 we illustrate future research lines on the subject.

Generalized Minimum Time Function

Definition 1 (Standing Assumptions). We will say that a set-valued function F : R d ⇒ R d satisfies the assumption (F j ), j = 0, 1, 2 if the following hold true

(F 0 ) F (x) = ∅ is compact and convex for every x ∈ R d , moreover F (•) is contin-
uous with respect to the Hausdorff metric, i.e. given x ∈ X, for every ε > 0 there exists δ > 0 such that |y -x| ≤ δ implies F (y) ⊆ F (x) + B(0, ε) and

F (x) ⊆ F (y) + B(0, ε). (F 1 ) F (•) has linear growth, i.e. there exist nonnegative constants L 1 and L 2 such that F (x) ⊆ B(0, L 1 |x| + L 2 ) for every x ∈ R d , (F 2 ) F (•) is bounded, i.e. there exist M > 0 such that y ≤ M for all x ∈ R d , y ∈ F (x). Definition 2 (Generalized target). Let p ≥ 1, Φ ⊆ C 0 (R d , R
) such that the following property holds (T E ) there exists x 0 ∈ R d with φ(x 0 ) ≤ 0 for all φ ∈ Φ.

We define the generalized target SΦ p as follows

SΦ p := µ ∈ P p (R d ) : R d φ(x) dµ(x) ≤ 0 for all φ ∈ Φ .
For an analysis of the properties of the generalized target see [START_REF] Cavagnari | Time optimal control problem in the space of probability measures[END_REF] or [START_REF] Cavagnari | Generalized control systems in the space of probability measures[END_REF] for deeper results.

Definition 3 (Admissible curves). Let

F : R d ⇒ R d be a set-valued function, I = [a, b] a compact interval of R, α, β ∈ P p (R d ).
We say that a Borel family of probability measures µ = {µ t } t∈I ⊆ P p (R d ) is an admissible trajectory (curve) defined in I for the system Σ F joining α and β, if there exists a family of Borel vector fields v = {v t (•)} t∈I such that 1. µ is a narrowly continuous solution in the distributional sense of the continuity equation ∂ t µ t + div(v t µ t ) = 0, with µ |t=a = α and µ |t=b = β.

2. J F (µ, v) < +∞, where J F (•) is defined as

JF (µ, v) :=        I R d 1 + I F (x) (vt(x)) dµt(x) dt, if vt L 1 µ t ∈ L 1 ([0, T ]), +∞, otherwise, (4)
where I F (x) is the indicator function of the set F (x), i.e., I F (x) (ξ) = 0 for all ξ ∈ F (x) and I F (x) (ξ) = +∞ for all ξ / ∈ F (x).

In this case, we will also shortly say that µ is driven by v.

When J F (•) is finite, this value expresses the time needed by the system to steer α to β along the trajectory µ with family of velocity vector fields v. Definition 4 (Generalized minimum time). Given p ≥ 1, let Φ ∈ C 0 (R d ; R) and SΦ p be the corresponding generalized target defined in Definition 2. In analogy with the classical case, we define the generalized minimum time function T Φ p :

P p (R d ) → [0, +∞] by setting T Φ p (µ 0 ) := inf {J F (µ, v) : µ is an admissible curve in [0, T ] , (5)
driven by v, with

µ |t=0 = µ 0 , µ |t=T ∈ SΦ p ,
where, by convention, inf ∅ = +∞. Given µ 0 ∈ P p (R d ), an admissible curve µ = {µ t } t∈[0, T Φ p (µ0)] ⊆ P p (R d ), driven by a time depending Borel vector-field v = {v t } t∈[0, T Φ p (µ0)] and satisfying µ |t=0 = µ 0 and µ |t= T Φ p (µ0) ∈ SΦ p is optimal for µ 0 if

T Φ p (µ 0 ) = J F (µ, v).
Some interesting results concerning the generalized minimum time function together with comparisons with the classical definition are proved in the proceedings [START_REF] Cavagnari | Time optimal control problem in the space of probability measures[END_REF] and in the forthcoming paper [START_REF] Cavagnari | Generalized control systems in the space of probability measures[END_REF].

Here we will focus our attention on the problem of finding an Hamilton-Jacobi-Bellman equation for our time-optimal control problem.

First of all we need to state and prove a Dynamic Programming Principle and, for this aim, the gluing result for solutions of the continuity equation stated in Lemma 4.4 in [START_REF] Dolbeault | A new class of transport distances between measures[END_REF] will be used.

Theorem 1 (Dynamic programming principle). Let p ≥ 1, 0 ≤ s ≤ τ , let F : R d ⇒ R d be a set-valued function, let µ = {µ t } t∈[0,τ ] be an admissible curve for Σ F . Then we have

T Φ p (µ 0 ) ≤ s + T Φ p (µ s ).
Moreover, if T Φ p (µ 0 ) < +∞, equality holds for all s ∈ [0, T Φ p (µ 0 )] if and only if µ is optimal for µ 0 = µ |t=0 .

Proof. The proof is based on the fact that, by Lemma 4.4 in [START_REF] Dolbeault | A new class of transport distances between measures[END_REF], the juxtaposition of admissible curves is an admissible curve. Thus, for every ε > 0 we consider the curve obtained by following µ up to time s, and then following an admissible curve steering µ s to the generalized target in time T Φ p (µ s ) + ε. We obtain an admissible curve steering µ 0 to the generalized target in time s + T Φ p (µ s ) + ε, and so, by letting ε → 0 + , we have T Φ p (µ 0 ) ≤ s + T Φ p (µ s ).

Assume now that T Φ p (µ 0 ) < +∞ and equality holds for all s ∈ [0, T Φ p (µ 0 )]. By taking s = T Φ p (µ 0 ) we get T Φ p (µ T Φ p (µ0) ) = 0, i.e., µ T Φ p (µ0) ∈ SΦ p . In particular, µ steers µ 0 to SΦ p in time T Φ p (µ 0 ), which is the infimum among all admissible trajectories steering µ 0 to the generalized target. So µ is optimal.

Finally, assume that µ is optimal for µ 0 and T Φ p (µ 0 ) < +∞. Starting from µ 0 , we follow µ up to time s. Since µ is still an admissible curve steering µ s to the generalized target in time T Φ p (µ 0 ) -s, we must have T Φ p (µ s ) ≤ T Φ p (µ 0 ) -s, and so T Φ p (µ s ) + s = T Φ p (µ 0 ), since the reverse inequality always holds true.

Hamilton-Jacobi-Bellman Equation

In this section we will prove that, under some assumptions, the generalized minimum time functional T Φ 2 is a viscosity solution, in a sense we will precise, of a suitable Hamilton-Jacobi-Bellman equation on P 2 (R d ). In this paper we assume the velocity field to be continuous for simplicity. In the forthcoming paper [START_REF] Cavagnari | Generalized control systems in the space of probability measures[END_REF] we prove a result of approximation of L 2 µ -selections of F with continuous and bounded ones in L 2 µ -norm that allows us to treat a more general case.

We recall that, given T ∈]0, +∞], the evaluation operator e t : R d × Γ T → R d is defined as e t (x, γ) = γ(t) for all 0 ≤ t < T . We set

T F (µ 0 ) := {η ∈ P(R d × Γ T ) : T > 0, η concentrated on trajectories of γ(t) = v(γ(t)), with v ∈ C 0 (R d ; R d ), v(x) ∈ F (x) ∀x ∈ R d and satisfies γ(0) η = µ 0 }, where µ 0 ∈ P 2 (R d ).
It is not hard to prove the following result.

Lemma 1 (Properties of the evaluation operator). Assume (F 0 ) and (F 1 ), and let L 1 , L 2 > 0 be the constants as in (F 1 ). For any

µ 0 ∈ P 2 (R d ), T ∈]0, 1], η ∈ T F (µ 0 ), we have: (i) |e t (x, γ)| ≤ (|e 0 (x, γ)| + L 2 ) e L1 for all t ∈ [0, T ] and η-a.e. (x, γ) ∈ R d × Γ T ; (ii) e t ∈ L 2 η (R d × Γ T ; R d ) for all t ∈ [0, T ];
(iii) there exists C > 0 depending only on L 1 , L 2 such that for all t ∈ [0, T ] we have e t -e 0 t 2

L 2 η ≤ C (m 2 (µ 0 ) + 1) .
In the case we are considering, where the trajectory t → e t η is driven by a sufficiently smooth velocity field, we recover as initial velocity what we expected.

Lemma 2 (Regular driving vector fields). Let µ = {µ t } t∈[0,T ] be an absolutely continuous solution of

     ∂ t µ t + div(vµ t ) = 0, t ∈]0, T [ µ |t=0 = µ 0 ∈ P 2 (R d ), where v ∈ C 0 b (R d ; R d ) satisfies v(x) ∈ F (x) for all x ∈ R d . Then if η ∈ T F (µ 0 ) satisfies µ t = e t η for all t ∈ [0, T ], we have that lim t→0 e t -e 0 t -v • e 0 L 2 η = 0 .
The proof is based on the boundedness of v and on the fact that, by hypothesis, γ ∈ C 1 , γ(t) = v(γ(t)). The conclusion comes applying Lebesgue's Dominated Convergence Theorem.

We give now the definitions of viscosity sub/-superdifferential and viscosity solutions that suit our problem. As presented in the Introduction, these concepts are different from the ones treated in [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF][START_REF] Ambrosio | On a class of first order Hamilton-Jacobi equations in metric spaces[END_REF][START_REF] Cardaliaguet | Deterministic differential games under probability knowledge of initial condition[END_REF][START_REF] Gangbo | Hamilton-Jacobi Equations in the Wasserstein Space[END_REF][START_REF] Gangbo | Optimal transport and large number of particles[END_REF], due mainly to the structure of P 2 (R d ).

we have A(s, p µ0 , η) ≤ B(s, p µ0 , η).

We recall that since by definition p µ0 ∈ L 2 µ0 , we have that p µ0 • e 0 ∈ L 2 η . Dividing by s > 0 the left hand side, we observe that we can use Lemma 2, indeed the velocity field v(•) associated to η ∈ T F (µ 0 ) satisfies all the hypothesis (the boundedness comes from (F 2 )) and so we have lim sup

s→0 + A(s, p µ0 , η) s = -1 - R d ×Γ T p µ0 • e 0 (x, γ), v • e 0 (x, γ) dη(x, γ).
Recalling that p µ0 ∈ D + δ T Φ 2 (µ 0 ) and using Lemma 1(iii), we have lim sup

s→0 + B(s, p µ0 , η) s = lim sup s→0 + B(s, p µ0 , η) e s -e 0 L 2 η • e s -e 0 s L 2 η ≤ Cδ,
where C > 0 is a suitable constant (we can take twice the upper bound on F given by (F 2 )).

We thus obtain for all η ∈ T F (µ 0 ) that

1 + R d ×Γ T p µ0 • e 0 (x, γ), v • e 0 (x, γ) dη(x, γ) ≥ -Cδ.
By passing to the infimum on η ∈ T F (µ 0 ) we have

-Cδ 1 + inf η∈T F (µ0) R d ×Γ T p µ0 • e 0 (x, γ), v • e 0 (x, γ) dη(x, γ) = 1 + inf η∈T F (µ0) R d p µ0 (x), v(x) dµ 0 (x) = -H F (µ 0 , p µ0 ), so T Φ 2 (•) is a subsolution, thus confirming Claim 1. Claim 2 : T Φ 2 (•) is a supersolution of H F (µ, D T Φ 2 (µ)) = 0. Proof of Claim 2. Given η ∈ T F (µ 0 ), let us define the admissible trajectory µ = {µ t } t∈[0,T ] = {e t η} t∈[0,T ] . Given q µ0 ∈ D - δ T Φ 2 (µ 0 ) ∩ C 0 b , we have R d ×Γ T q µ0 • e 0 (x, γ), e s (x, γ) -e 0 (x, γ) s dη(x, γ) ≤ 2δ e s -e 0 s L 2 η - T Φ 2 (µ 0 ) -T Φ 2 (µ s ) s .
Thus, using Lemma 2 and Lemma 1, we have

R d ×Γ T
q µ0 • e 0 (x, γ), v • e 0 (x, γ) dη(x, γ) ≤ 3Cδ -

T Φ 2 (µ 0 ) -T Φ 2 (µ s ) s ,
for all s > 0.

Then, by passing to the infimum on all admissible trajectories, we obtain -H F (µ 0 , q µ0 ) -1 = inf η∈T F (µ0) R d ×Γ T q µ0 • e 0 (x, γ), v • e 0 (x, γ) dη(x, γ)

≤ 3Cδ -sup η∈T F (µ0)
T Φ 2 (µ 0 ) -T Φ 2 (µ s ) s .

Thus

H F (µ 0 , q µ0 ) ≥ -3Cδ + sup η∈T F (µ0)

T Φ 2 (µ 0 ) -T Φ 2 (µ s ) s -1 .
By the Dynamic Programming Principle, recalling that

T Φ 2 (µ 0 ) -T Φ 2 (µ s ) s -1 ≤ 0
with equality holding if and only if µ is optimal, we obtain H F (µ 0 , q µ0 ) ≥ -C δ, which proves that T Φ 2 (•) is a supersolution, thus confirming Claim 2.

Conclusion

In this work we have studied a Hamilton-Jacobi-Bellman equation solved by a generalized minimum time function in a regular case. In the forthcoming paper [START_REF] Cavagnari | Generalized control systems in the space of probability measures[END_REF] an existence result is proved for optimal trajectories as well as attainability properties in the space of probability measures. Furthermore, a suitable approximation result allows to give a sense to a Hamilton-Jacobi-Bellman equation in a more general case.

We plan to study if it is possible to prove a comparison principle for an Hamilton-Jacobi equation solved by the generalized minimum time function, as well as to give a Pontryagin maximum principle for our problem.

2 )

 2 minimize the time needed to steer x 0 to a given closed target set S ⊆ R d , S = ∅, defining the minimum time function T : R d → [0, +∞] by T (x 0 ) := inf{T > 0 : ∃ x(•) solving (1) such that x(T ) ∈ S}. (Acknowledgements: The first two authors have been supported by INdAM -GNAMPA Project 2015: Set-valued Analysis and Optimal Transportation Theory Methods in Deterministic and Stochastics Models of Financial Markets with Transaction Costs.

Definition 5 (Sub-/Super-differential in P 2 (R d )). Let V : P 2 (R d ) → R be a function. Fix µ ∈ P 2 (R d ) and δ > 0. We say that p µ ∈ L 2 µ (R d ; R d ) belongs to the δ-superdifferential D + δ V (µ) at µ if for all T > 0 and η ∈ P(R d × Γ T ) such that t → e t η is an absolutely continuous curve in P 2 (R d ) defined in [0, T ] with e 0 η = µ we have

In the same way,

Definition 6 (Viscosity solutions). Let V : P 2 (R d ) → R be a function and

3. viscosity solution of H (µ, DV (µ)) = 0 if it is both a viscosity subsolution and a viscosity supersolution.

Definition 7 (Hamiltonian Function). Given µ ∈ P 2 (R d ), we define the map

Theorem 2 (Viscosity solution). Assume (F 0 ) and (F 2 ). Then

Proof. The proof is splitted in two claims.

Given η ∈ T F (µ 0 ) and set µ t = e t η for all t by the Dynamic Programming Principle we have T Φ 2 (µ 0 ) ≤ T Φ 2 (µ s ) + s for all 0 < s ≤ T Φ 2 (µ 0 ). Without loss of generality, we can assume 0 < s < 1. Given any


Konstantin Palagachev
email: konstantin.palagachev@unibw.de

Matthias Gerdts
email: matthias.gerdts@unibw.de

Exploitation of the value function in a bilevel optimal control problem

Keywords: bilevel optimal control, value function, pursuit-evasion Stackelberg game

The paper discusses a class of bilevel optimal control problems with optimal control problems at both levels. The problem will be transformed to an equivalent single level problem using the value function of the lower level optimal control problem. Although the computation of the value function is difficult in general, we present a pursuit-evasion Stackelberg game for which the value function of the lower level problem can be derived even analytically. A direct discretization method is then used to solve the transformed single level optimal control problem together with some smoothing of the value function.

Introduction

Bilevel optimization problems occur in various applications, e.g. in locomotion and biomechanics, see [START_REF] Hatz | Efficient numerical methods for hierarchical dynamic optimization with application to cerebral palsy gait modeling[END_REF][START_REF] Mombaur | Stability optimization of open-loop controlled walking robots[END_REF][START_REF] Albrecht | A bilevel optimization approach to obtain optimal cost functions for human arm movements[END_REF][START_REF] Albrecht | Modeling and numerical solution of inverse optimal control problems for the analysis of human motions[END_REF], in optimal control under safety constraints, see [START_REF] Knauer | Bilevel-Optimalsteuerung mittels hybrider Lösungsmethoden am Beispiel eines deckengeführten Regalbediengerätes in einem Hochregallager[END_REF][START_REF] Knauer | Fast and save container cranes as bilevel optimal control problems[END_REF][START_REF] Fisch | Development of a Framework for the Solution of High-Fidelity Trajectory Optimization Problems and Bilevel Optimal Control Problems[END_REF], or in Stackelberg dynamic games, compare [START_REF] Stackelberg | The Theory of Market Economy[END_REF][START_REF] Ehtamo | On applied nonlinear and bilevel programming for pursuit-evasion games[END_REF]. An abstract bilevel optimization problem (BOP) reads as follows:

Minimize F (x, y) with respect to (x, y) ∈ X × Y subject to the constraints G(x, y) ∈ K, H(x, y) = 0, y ∈ M (x),

where M (x) is the set of minimizers of the lower level optimization problem Minimize f (x, y) w.r.t. y ∈ Y s.t. g(x, y) ∈ C, h(x, y) = 0.

Herein, X, Y are (finite or infinite) Banach spaces, F, f :

X × Y → R, H : X × Y → V u , h : X × Y → V , G : X × Y → W u , g : X × Y → W are
sufficiently smooth functions into Banach spaces V u , V , W u , W , and K ⊂ W u , C ⊂ W are convex and closed cones. Bilevel optimization problems turn out to be very challenging with regard to both, the investigation of theoretical properties and numerical methods, compare [START_REF] Dempe | Foundations of bilevel programming[END_REF]. Necessary conditions have been investigated, e.g., in [START_REF] Ye | Necessary conditions for bilevel dynamic optimization problems[END_REF][START_REF] Dempe | KKT reformulation and necessary conditions for optimality in nonsmooth bilevel optimization[END_REF]. Typical solution approaches aim at reducing the bilevel structure into a single stage optimization problem. In the MPCC approach a single level optimization problem subject to complementarity constraints (MPCC) is obtained by replacing the lower level problem by its first order necessary conditions, compare [START_REF] Albrecht | Modeling and numerical solution of inverse optimal control problems for the analysis of human motions[END_REF]. However, if the lower level problem is non-convex, the MPCC is not equivalent in general to the original bilevel problem since non-optimal stationary points or non-global solutions may satisfy the necessary conditions as well. Still, the approach is often used owing to a well-established theory and the availability of numerical methods for MPCCs, especially for finite dimensional problems.

In this paper we focus on an equivalent transformation of the bilevel problem to a single level problem (see [START_REF] Dempe | A new equivalent single-level problem for bilevel problems[END_REF] for an alternative way). The equivalence can be guaranteed by exploitation of the value function V : X → R of the lower level problem, which is defined as

V (x) := inf y∈Y {f (x, y) | g(x, y) ∈ C, h(x, y) = 0}.
An equivalent reformulation of the bilevel optimization problem is then given by the following single level problem, compare [START_REF] Outrata | On the numerical solution of a class of Stackelberg problems[END_REF][START_REF] Ye | Necessary conditions for bilevel dynamic optimization problems[END_REF][START_REF] Ye | Optimal strategies for bilevel dynamic problems[END_REF]:

Minimize F (x, y) w.r.t. (x, y) ∈ X × Y subject to the constraints G(x, y) ∈ K, H(x, y) = 0, g(x, y) ∈ C, h(x, y) = 0, f (x, y) ≤ V (x).
The advantage of the value function approach is its equivalence with the bilevel problem. On the downside one has to be able to compute the value function, which in general might be intractable. Moreover, the value function is non-smooth in general (often Lipschitz continuous) and hence suitable methods from non-smooth optimization are required to solve the resulting single level problem. In Section 2 we discuss a class of bilevel optimal control problems that fit into the problem class BOP. In Section 3 we we are able to derive an analytical expression for the value function for an example and present numerical results. The new contribution of this paper is the discussion of a particular example, which combines the analytical expression of the value function of the lower level problem and a direct discretization method for the reformulated single level problem. This problem may serve as a test problem for theoretical and numerical investigations. The problem exhibits already most features of more challenging problems such as non-convexity, pure state constraints on the upper level problem as well as control constraints on both levels.

A class of bilevel optimal control problems

Let T > 0, be the fixed final time,

X := W 1,∞ ([0, T ], R nx ) × L ∞ ([0, T ], R nu ) × R np , n x , n u , n p ∈ N 0 , Y := W 1,∞ ([0, T ], R ny )×L ∞ ([0, T ], R nv )×R nq , n y , n v , n q ∈ N 0 , where L ∞ ([0, T ], R n )
denotes the Banach space of essentially bounded vectorvalued functions from [0, T ] into R n and W 1,∞ ([0, T ], R n ) is the Banach space of absolutely continuous vector-valued functions from [0, T ] into R n with essentially bounded first derivatives. Moreover, let the Banach spaces

V u := L ∞ ([0, T ], R nx ) × R n H , V := L ∞ ([0, T ], R ny ) × R n h , n H , n h ∈ N 0 , and the closed convex cones W u := {k ∈ L ∞ ([0, T ], R n G ) | k(t) ≤ 0 a.e. in [0, T ]}, W := {k ∈ L ∞ ([0, T ], R ng ) | k(t) ≤ 0 a.e. in [0, T ]}, n G , n g ∈ N 0 , be given. Let J, j : R nx × R ny × R np × R nq → R, F : R nx × R ny × R nu × R nv × R np × R nq → R nx , f : R ny × R nv × R np × R nq → R ny , Ψ : R nx × R ny × R nx × R ny × R np × R nq → R n H , ψ : R nx × R ny × R nx × R ny × R np × R nq → R n h , S : R nx × R ny × R nu × R nv × R np × R nq → R n G , s : R ny × R nv × R np × R nq → R ng .
be sufficiently smooth mappings. With these definitions the following class of bilevel optimal control problems (BOCP) subject to control-state constraints and boundary conditions fits into the general bilevel optimization problem BOP.

Minimize J(x(T ), y(T ), p, q) w.r.t. (x, u, p, y, v, q) ∈ X × Y subject to the constraints x (t) = F (x(t), y(t), u(t), v(t), p, q), (1) 0
= Ψ (x(0), y(0), x(T ), y(T ), p, q), (2)
0 ≥ S(x(t), y(t), u(t), v(t), p, q), (3)
(y, v, q) ∈ M (x(0), x(T ), p)

where M (x(0), x(T ), p) is the set of minimizers of the lower level problem OCP L (x(0), x(T ), p):

Minimize j(x(T ), y(T ), p, q) w.r.t. (y, v, q) ∈ Y subject to the constraints

y (t) = f (y(t), v(t), p, q), (4) 0
= ψ(x(0), y(0), x(T ), y(T ), p, q), (5)
0 ≥ s(y(t), v(t), p, q). (6)
Herein, (x, u, p) ∈ X are the state, the control, and the parameter vector of the upper level problem and (y, v, q) ∈ Y are the state, the control, and the parameter vector of the lower level problem. Please note that the lower level problem only depends on the initial and terminal states x(0), x(T ) and the parameter vector p of the upper level problem. The value function V is then a mapping from

R nx × R nx × R np into R defined by V (x 0 , x T , p) := inf (y,v,q)∈Y    j(x T , y(T ), p, q) y (t) = f (y(t), v(t), p, q), 0 = ψ(x 0 , y(0), x T , y(T ), p, q), 0 ≥ s(y(t), v(t), p, q)    .
Remark 1. In a formal way the problem class can be easily extended in such a way that the lower level dynamics f and the lower level control-state constraints s depend on x, u as well. However, in the latter case the value function of the lower level problem would then be a functional V : X → R, i.e. a functional defined on the Banach space X rather than a functional defined on the finite dimensional space R nx × R nx × R np . Computing the mapping V : X → R numerically would be computationally intractable in most cases.

Using the value function V we arrive at the following equivalent single level optimal control problem subject to control-state constraints, smooth boundary conditions, and an in general non-smooth boundary condition with the value function.

Minimize J(x(T ), y(T ), p, q) w.r.t. (x, u, p, y, v, q) ∈ X × Y subject to the constraints (1)-( 3), ( 4)-( 6), and j(x(T ), y(T ), p, q) ≤ V (x(0), x(T ), p).

(

) 7
It remains to compute the value function V and to solve the potentially nonsmooth single level optimal control problem. Both are challenging tasks owing to non-smoothness and non-convexity. The value function sometimes can be derived analytically as we shall demonstrate in Section 3. Otherwise, if Bellman's optimality principle applies, the value function satisfies a Hamilton-Jacobi-Bellman (HJB) equation, see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]. Various methods exist for its numerical solution, compare [START_REF] Osher | High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF][START_REF] Grüne | An adaptive grid scheme for the discrete hamilton-jacobi-bellman equation[END_REF][START_REF] Falcone | Convergence analysis for a class of high-order semi-Lagrangian advection schemes[END_REF][START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF][START_REF] Bokanowski | Reachability analysis and optimal control problems -Hamilton-Jacobi equations[END_REF]. The HJB approach is feasible if the state dimension n y does not exceed 5 or 6. If no analytical formula is available and if the HJB approach is not feasible, then a pointwise evaluation of V at (x(0), x(T ), p) can be realized by using suitable optimal control software, e.g. [START_REF] Gerdts | OCPID-DAE1 -optimal control and parameter identification with differential-algebraic equations of index 1[END_REF]. However, if the lower level problem is non-convex, then it is usually not possible to guarantee global optimality by such an approach. The single level problem can be approached by the non-smooth necessary conditions in [START_REF] De Pinho | Necessary Conditions for Optimal Control Problems Involving Nonlinear Differential Algebraic Equations[END_REF][START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF]. Alternatively, direct discretization methods may be applied. The non-smoothness in V in [START_REF] Dempe | A new equivalent single-level problem for bilevel problems[END_REF] has to be taken into account by, e.g., using bundle type methods, see [START_REF] Schramm | A version of the bundle idea for minimizing a nonsmooth function: Conceptual idea, convergence analysis, numerical results[END_REF], or by smoothing the value function and applying standard software. Finally, the HJB approach could also be applied to the single level problem again.

A follow-the-leader application

We consider a pursuit-evasion dynamic Stackelberg game of two vehicles moving in the plane. Throughout we assume that the evader knows the optimal strategy of the pursuer and can optimize its own's strategy accordingly. This gives rise to a bilevel optimal control problem. The lower level player (=pursuer P) aims to capture the upper level player (=evader E) in minimum time T . The evader aims to minimize a linear combination of the negative capture time -T and its control effort. The players have individual dynamics and constraints. The coupling occurs through capture conditions at the final time.

The bilevel optimal control problem

The evader E aims to solve the following optimal control problem, called the upper level problem (OCP U ):

Minimize -T + T 0 α 1 2 w(t) 2 + α 2 2 a(t) 2 dt (8)
subject to the constraints

x E (t) = v E (t) cos ψ(t), x E (0) = x E,0 , x E (T ) = x P (T ), (9)
y E (t) = v E (t) sin ψ(t), y E (0) = y E,0 , y E (T ) = y P (T ), (10)
ψ (t) = v E (t) tan δ(t), ψ(0) = ψ 0 , (11)
δ (t) = w(t), δ(0) = δ 0 , (12)
v E (t) = a(t), v E (0) = v E,0 , (13)
v E (t) ∈ [0, v E,max ], w(t) ∈ [-w max , w max ], a(t) ∈ [a min , a max ], (14)
(x P , y P , T ) ∈ M (x E (T ), y E (T )),
where M (x E (T ), y E (T )) denotes the set of minimizers of the lower level problem OCP L (x E (T ), y E (T )) below.

The equations of motion of E describe a simplified car model of length > 0 moving in the plane. The controls are the steering angle velocity w and the acceleration a with given bounds ±w max , a min , and a max , respectively. The velocity v E is bounded by the state constraint v E (t) ∈ [0, v E,max ] with a given bound v E,max > 0. The position of the car's rear axle is given by z E = (x E , y E ) and its velocity by v E . ψ denotes the yaw angle and α 1 , α 2 ≥ 0 are weights in the objective function. The initial state is fixed by the values x E,0 , y E,0 , ψ 0 , δ 0 , v E,0 . The final time T is determined by the lower level player P, who aims to solve the following optimal control problem, called the lower level problem OCP L (x E,T , y E,T ) with its set of minimizers denoted by M (x E,T , y E,T ):

Minimize T = T 0 1dt subject to the constraints z P (t) = v P (t), z P (0) = z P,0 , z P (T ) = (x E,T , y E,T ) , (15)
v P (t) = u P (t), v P (0) = v P (T ) = 0, (16)
u P,i (t) ∈ [-u max , u max ], i = 1, 2. (17)
Herein, z P = (x P , y P ) , v P = (v P,1 , v P,2 ) , and u P = (u P,1 , u P,2 ) denote the position vector, the velocity vector, and the acceleration vector, respectively, of P in the two-dimensional plane. z P,0 = (x P,0 , y P,0 ) ∈ R 2 is a given initial position. u max > 0 is a given control bound for the acceleration. The dynamics of the pursuer allow to move in x and y direction independently, which models, e.g., a robot with omnidirectional wheels.

The boundary conditions for vP,i ( T ) and ẑP,i ( T ) yield Ti = 2 ts,i and ts,i = ẑP,i ( T ) -ẑP,i (0) u max if ẑP,i ( T ) -ẑP,i (0) ≥ 0.

Case 3: The switching structure for control component

i ∈ {1, 2} is ûP,i (t) = -u max , if 0 ≤ t < ts,i , u max , otherwise.
This case can be handled analogously to Case 2 and we obtain Ti = 2 ts,i and ts,i = ẑP,i (0) -ẑP,i ( T )

u max if ẑP,i (0) -ẑP,i ( T ) ≥ 0.
The above analysis reveals the shortest times Ti , i ∈ {1, 2}, in which the i-th state can reach its terminal boundary condition. The minimum time T for a given terminal position is thus given by the value function V of OCP L (x E,T , y E,T ) (=minimum time function) with

V (x E,T , y E,T ) = max{ T1 , T2 } = 2 max    |x P,0 -x E,T | u max , |y P,0 -y E,T | u max    . (18)
That is, the final time is defined by the component i with the largest distance |ẑ P,i ( T ) -ẑP,i (0)|. For this component, the control is of bang-bang type with one switch at the midpoint of the time interval. The remaining control can be singular and it is not uniquely defined. The value function is locally Lipschitz continuous except at the point (x E,T , y E,T ) = (x P,0 , y P,0 ), compare Figure 1. This point, however, is of minor interest because interception takes place immediately. The equivalent single level problem (SL-OCP) reads as follows:

Minimize (8) subject to the constraints ( 9)-( 14), ( 15)-( 17) with (x E,T , y E,T ) = (x E (T ), y E (T )) and the non-smooth constraint

T ≤ V (x E (T ), y E (T )). ( 19
)

Numerical results

For the numerical solution of the single level problem SL-OCP we applied the direct shooting method OCPID-DAE1, [START_REF] Gerdts | OCPID-DAE1 -optimal control and parameter identification with differential-algebraic equations of index 1[END_REF]. The non-smooth constraint T ≤ V (x E (T ), y E (T )) with V from ( 18) was replaced by a continuously differentiable constraint which was obtained by smoothing the maximum function and the absolute value function in [START_REF] Knauer | Bilevel-Optimalsteuerung mittels hybrider Lösungsmethoden am Beispiel eines deckengeführten Regalbediengerätes in einem Hochregallager[END_REF]. Figure 2 shows a numerical solution of the pursuitevasion Stackelberg bilevel optimal control problem for the data v E,0 = 10, ψ E (0) = π/4, α 1 = 10, α 2 = 0, w max = 0.5, v E,max = 20, a min = -5, a max = 1, u max = 5, N = 50, T ≈ 18.01. Figure 3 shows several trajectories for the pursuer and the evader for different initial yaw angles covering the interval [0, 2π).

-1 -0.5 0 0.5 Remark 2. The constraint (19) may become infeasible under discretization. Instead, the value function V h of the discretized lower level optimal control problem should be used. However, since V h is hardly available for all kinds of discretizations, we use instead the relaxed constraint T ≤ V (x E (T ), y E (T )) + ε with some ε > 0.

Conclusions and Outlook

The paper discusses a specific bilevel optimal control problem and its reformulation as an equivalent single level problem using the value function of the lower level problem. For a sample problem it is possible to compute the value function analytically and to solve the overall bilevel problem numerically using a direct discretization method. This first numerical study leaves many issues open that have to be investigated in future research for the general problem setting. Amongst them are smoothness properties of the value function, representation of subdifferentials, the development of appropriate solution methods for nonsmooth problems, and the derivation of necessary (and sufficient) conditions of optimality for the class of bilevel optimal control problems.

  The work is supported by Munich Aerospace e.V.

Fig. 1 .

 1 Fig. 1. Value function of lower level problem with data xP,0 = yP,0 = 0, umax = 1.

Fig. 3 .

 3 Fig. 3. Left picture: Trajectories of pursuer (lines with '+') and evader (lines with boxes) for ψE(0) = π/4. Right picture: Trajectories of the pursuer (lines with '+') and the evader (lines with boxes) for different initial yaw angles of the evader.

The lower-level problem and its value function

The lower level problem admits an analytical solution. To this end, the Hamilton function (regular case only) reads as

The first order necessary optimality conditions for a minimum (ẑ P , vP , ûP , T ) are given by the minimum principle, compare [START_REF] Ioffe | Theory of extremal problems[END_REF]. There exist adjoint multipliers λ z , λ v with

and

In this case, the minimum principle provides no information on the singular control except feasibility. Notice furthermore that not all control components can be singular since this would lead to trivial multipliers in contradiction to the minimum principle. Hence, there is at least one index i for which the control component ûP,i is non-singular. In the non-singular case there can be at most one switch of each component ûP,i , i ∈ {1, 2}, in the time interval [0, T ], since λ v,i is linear in time. The switching time ts,i for the i-th control component computes to ts,i = c v,i /c z,i if c z,i = 0. We discuss several cases for non-singular controls.

Case 1: No switching occurs in ûP,i , i.e. ûP,i (t) ≡ ±u max for i ∈ {1, 2}. By integration we obtain vP,i (t) = ±u max t and thus vP,i ( T ) = 0 in contradiction to the boundary conditions. Consequently, each non-singular control component switches exactly once in [0, T ].

Case 2: The switching structure for control component

By integration and the boundary conditions we find vP,i (t) = u max t, if 0 ≤ t < ts,i u max (2 ts,i -t), otherwise ẑP,i (t) = ẑP,i (0) + 1 2 u max t 2 , if 0 ≤ t < ts,i ẑP,i (0) + u max t2 s,i -1 2 (2 ts,i -t) 2 , otherwise.


Mikhail Sumin
email: m.sumin@mail.ru

Stable Sequential Pontryagin Maximum Principle as a Tool for Solving Unstable Optimal Control and Inverse Problems for Distributed Systems

Keywords: Optimal boundary control, parabolic equation, minimizing sequence, dual regularization, stability, Pontryagin maximum principle

This article is devoted to studying dual regularization method as applied to parametric convex optimal control problem of controlled third boundary-value problem for parabolic equation with boundary control and with equality and inequality pointwise state constraints. These constraints are understood as ones in the Hilbert space L2. A major advantage of the constraints of the original problem which are understood as ones in L2 is that the resulting dual regularization algorithm is stable with respect to errors in the input data and leads to the construction of a minimizing approximate solution in the sense of J. Warga. Simultaneously, this dual algorithm yields the corresponding necessary and sufficient conditions for minimizing sequences, namely, the stable, with respect to perturbation of input data, sequential or, in other words, regularized Lagrange principle in nondifferential form and Pontryagin maximum principle for the original problem. Regardless of the fact that the stability or instability of the original optimal control problem, they stably generate a minimizing approximate solutions for it. For this reason, we can interpret these regularized Lagrange principle and Pontryagin maximum principle as tools for direct solving unstable optimal control problems and reducing to them unstable inverse problems.

Introduction

Pontryagin maximum principle is the central result of all optimal control theory, including optimal control for differential equations with partial derivatives. Its statement and proof assume, first of all, that the optimal control problem is considered in an ideal situation, when its input data are known exactly. However, in the vast number of important practical problems of optimal control, as well as numerous problems reducing to optimal control problems, the requirement of exact defining input data is very unnatural, and in many undoubtedly interest cases is simply impracticable. In similar problems, we can not, strictly speaking, to take as an approximation to the solution of the initial (unperturbed) problem with the exact input data, a control formally satisfying the maximum principle in the perturbed problem. The reason of such situation lies in the natural instability of optimization problems with respect to perturbation of its input data. As a typical property of optimization problems in general, including constrained ones, instability fully manifests itself in optimal control problems (see., e.g., [START_REF] Sumin | Stable Sequential Convex Programming in a Hilbert Space and its Application for Solving Unstable Problems[END_REF]). As a consequence, the mentioned above instability implies "instability" of the classical optimality conditions, including the conditions in the form of Pontryagin maximum principle. This instability manifests itself in selecting by them of arbitrarily distant "perturbed" optimal elements from their unperturbed counterparts in the case of an arbitrarily small perturbations of the input data. The above applies, in full measure, both to discussed below optimal control problem with pointwise state constraints for linear parabolic equation in divergent form, and to the classical optimality conditions in the form of the Lagrange principle and the Pontryagin maximum principle for this problem.

In this paper we discuss how to overcome the problem of instability of the classical optimality conditions in optimal control problems in the way of applying dual regularization method (see., e.g., [START_REF] Sumin | A Regularized Gradient Dual Method for the Inverse Problem of a Final Observation for a Parabolic Equation[END_REF][START_REF] Sumin | Duality-Based Regularization in a Linear Convex Mathematical Programming Problem[END_REF][START_REF] Sumin | Regularized Parametric Kuhn-Tucker Theorem in a Hilbert Space[END_REF]) and simultaneous transition to the concept of minimizing sequence of admissible elements as the main concept of optimization theory. In the role of the last, acts the concept of the minimizing approximate solution in the sense of J.Warga [START_REF] Warga | Optimal Control of Differential and Functional Equations[END_REF]. The main attention in the paper is given to the discussion of the so-called regularized or, in other words, stable, with respect to perturbation of input data, sequential Lagrange principle in the nondifferential form and Pontryagin maximum principle. Regardless of the stability or instability of the original optimal control problem, they stably generate minimizing approximate solutions for it. For this reason, we can interpret the regularized Lagrange principle and Pontryagin maximum principle that are obtained in the article as tools for direct solving unstable optimal control problems and reducing to them unstable inverse problems [START_REF] Sumin | Stable Sequential Convex Programming in a Hilbert Space and its Application for Solving Unstable Problems[END_REF][START_REF] Sumin | Dual Regularization and Pontryagin's Maximum Principle in a Problem of Optimal Boundary Control for a Parabolic Equation with Nondifferentiable Functionals[END_REF][START_REF] Sumin | On the Stable Sequential Kuhn-Tucker Theorem and its Applications[END_REF]. Thus, they contribute to a significant expansion of the range of applicability of the theory of optimal control in which a central role belongs to classic constructions of the Lagrange and Hamilton-Pontryagin functions. Finally, we note that discussed in this article regularized Lagrange principle in the nondifferential form and Pontryagin maximum principle may have another kind, more convenient for applications [START_REF] Sumin | On the Stable Sequential Kuhn-Tucker Theorem and its Applications[END_REF]. Justification of these alternative forms of the regularized Lagrange principle and Pontryagin maximum principle is based on the so-called method of iterative dual regularization [START_REF] Sumin | A Regularized Gradient Dual Method for the Inverse Problem of a Final Observation for a Parabolic Equation[END_REF][START_REF] Sumin | Duality-Based Regularization in a Linear Convex Mathematical Programming Problem[END_REF]. In this case, they take the form of iterative processes with the corresponding stopping rules when the error of input data is fixed and finite. Here these alternative forms are not considered.

Statement of Optimal Control Problem

We consider the fixed-time parametric optimal control problem (P δ p,r ) [START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF][START_REF] Casas | Pontryagin's Principle for Local Solutions of Control Problems with Mixed Control-State Constraints[END_REF] to the third boundary-value problem

g δ 0 (π) → min, π ≡ (u, w) ∈ D ⊂ L 2 (Q T ) × L 2 (S T ) , g δ 1 (π)(x, t) ≡ ϕ δ 1 (x, t)z δ [π](x, t) = h δ (x, t) + p(x,
: L 2 (Q T ) × L 2 (S T ) is a continuous convex functional, Q ⊂ Q ι,T is a compact set without isolated points with a nonempty interior, ι ∈ (0, T ), Q = cl intQ; and z δ [π] ∈ V 1,0 2 (Q T ) ∩ C(Q T ) is a weak solution
1 z t - ∂ ∂x i (a i,j (x, t)z xj ) + a δ (x, t)z + u(x, t) = 0 , (1)
z(x, 0) = v δ 0 (x), x ∈ Ω, ∂z ∂N + σ δ (x, t)z = w(x, t), (x, t) ∈ S T ,
corresponding to the pair π ≡ (u, w). The superscript δ in the input data of Problem (P δ p,r ) indicates that these data are exact (δ = 0) or perturbed (δ > 0), i.e., they are specified with an error, δ ∈ [0, δ 0 ], where δ 0 > 0 is a fixed number.

For definiteness, as target functional we take terminal one

g δ 0 (π) ≡ Ω G δ (x, z δ [π](x, T ))dx .
The input data for Problem (P 0 p,r ) are assumed to meet the following conditions:

a) It is true that a i,j ∈ L ∞ (Q T ), i, j = 1, . . . , n, a δ ∈ L ∞ (Q T ), σ δ ∈ L ∞ (S T ), v δ 0 ∈ C(Ω), ν|ξ| 2 ≤ a i,j (x, t)ξ i ξ j ≤ µ|ξ| 2 ∀(x, t) ∈ Q T , ν, µ > 0 , a δ (x, t) ≥ C 0 for a.e. (x, t) ∈ Q T , σ δ (x, t) ≥ C 0 for a.e. (x, t) ∈ S T ; b) It is true that ϕ δ 1 , h δ ∈ L ∞ (Q); ϕ δ 2 : Q × R 1 → R 1 is
Lebesgue measurable function that is continuous and convex with respect to z for a.e. (x,

t) ∈ Q, ϕ δ 2 (•, •, z(•, •)) ∈ L ∞ (Q) ∀z ∈ C(Q); G δ : Ω × R 1 → R 1 is
Lebesgue measurable function that is continuous and convex with respect to z for a.e.

x ∈ Ω, G δ (•, z(•, T )) ∈ L ∞ (Ω) ∀z(•, T ) ∈ C(Q); c) Ω ⊂ R n be a bounded domain with Lipschitz boundary S.
Assume that the following estimates hold:

|G δ (x, z) -G 0 (x, z)| ≤ C M δ ∀ (x, z) ∈ Ω × S 1 M , ϕ δ 1 -ϕ 0 1 ∞,Q ≤ Cδ , (2)
h δ -h 0 ∞,Q ≤ Cδ, |ϕ δ 2 (x, t, z) -ϕ 0 2 (x, t, z)| ≤ C M δ ∀ (x, t, z) ∈ Q × S 1 M , a δ -a 0 ∞,Q T ≤ Cδ, |v δ 0 -v 0 0 | (0) Ω ≤ Cδ, σ δ -σ 0 ∞,S T ≤ Cδ ,
where C, C M > 0 are independent of δ; S n M ≡ {x ∈ R n : |x| < M }. Let's note, that the conditions on the input data of Problem (P δ p,r ), and also the estimates of deviations of the perturbed input data from the exact ones can be weakened.

In this paper we use for discussing the main results, related to the stable sequential Lagrange principle and Pontryagin maximum principle in Problem (P 0 p,r ), a scheme of studying the similar optimization problems in the papers [START_REF] Sumin | Parametric Dual Regularization for an Optimal Control Problem with Pointwise State Constraints[END_REF][START_REF] Sumin | Stable Sequential Pontryagin Maximum Principle in Optimal Control Problem with State Constraints[END_REF] for a system of controlled ordinary differential equations. In these works, both spaces of admissible controls and spaces, where lie images of the operators that define the pointwise state constraints, represented as Hilbert spaces of squareintegrable functions. For this reason, we put the set D of admissible controls π into a Hilbert space also, i.e., assume that

D ⊂ Z ≡ L 2 (Q T ) × L 2 (S T ), π ≡ ( u 2 2,Q T + w 2 2,S T ) 1/2
. At the same time, we note that the conditions on the input data of Problem (P δ p,r ) allow formally to consider that the operators g δ 1 , g δ 2 , specifying the state constraints of the problem, act into space L p (Q) with any index p ∈ [1, +∞]. However, in this paper, taking into account the above remark, we will put images of these functional operators in the Hilbert space L 2 (Q) ≡ H.

Suppose that Problem (P 0 p,r ) has a solution (which is unique if g 0 0 is strictly (strongly) convex). Its solutions are denoted by π 0 p,r ≡ (u 0 p,r , w 0 p,r ), and the set of all such solutions is designated as U 0 p,r . Define the Lagrange functional, a set of its minimizers and the concave dual problem

L δ p,r (π, λ, µ) ≡ g δ 0 (π) + λ, g δ 1 (π) -h δ -p + µ, g δ 2 (π) -r , π ∈ D , U δ [λ, µ] ≡ Argmin {L δ p,r (π, λ, µ) : π ∈ D} ∀(λ, µ) ∈ H × H + , V δ p,r (λ, µ) → sup, (λ, µ) ∈ H × H + , V δ p,r (λ, µ) ≡ inf π∈D L δ p,r (π, λ, µ) .
Since the Lagrange functional is continuous and convex for any pair (λ, µ) ∈ H × H + and the set D is bounded, the dual functional V δ p,r , is obviously defined and finite for any (λ, µ) ∈ H × H + .

The concept of a minimizing approximate solution in the sense of J.Warga [START_REF] Warga | Optimal Control of Differential and Functional Equations[END_REF] is of great importance for the design of a dual regularizing algorithm for problem (P 0 p,r ). Recall that a minimizing approximate solution is a sequence

π i ≡ (u i , w i ), i = 1, 2, . . . such that g 0 0 (π i ) ≤ β(p, r) + δ i , π i ∈ D 0, i p,r
for some nonnegative number sequences δ i and i , i = 1, 2, . . ., that converge to zero. Here, β(p, r) is the generalized infimum, i.e., an S-function:

β(p, r) ≡ lim →+0 β (p, r), β (p, r) ≡ inf π∈D 0, p,r g 0 0 (π), β (p, r) ≡ +∞ if D 0, p,r = ∅ , D δ, p,r ≡ {π ∈ D : g δ 1 (π) -h δ -p 2,Q ≤ , min z∈H- g δ 2 (π) -r -z 2,Q ≤ }, ≥ 0 , D 00 p,r ≡ D 0 p,r , H -≡ {z ∈ L 2 (Q) : z(x, t) ≤ 0 for a.e. (x, t) ∈ Q}, H + ≡ -H -.
Obviously, in the general situation, β(p, r) ≤ β 0 (p, r), where β 0 (p, r) is the classical value of the problem. However, in the case of Problem (P 0 p,r ), we have β(p, r) = β 0 (p, r). Simultaneously, we may asset that β :

L 2 (Q) × L 2 (Q) → R 1 ∪
{+∞} is a convex and lower semicontinuous function. Note here that the existence of a minimizing approximate solution in Problem (P 0 p,r ) obviously implies its solvability.

From the conditions a) -c) and the theorem on the existence of a weak solution of the third boundary-value problem for a linear parabolic equation of the divergent type (see [8, ch. III, section 5] and also [START_REF] Plotnikov | Existence and Uniqueness Theorems and a Priori Properties of Weak Solutions[END_REF]), it follows that the direct boundary-value problem (1) and the corresponding adjoint problem are

uniquely solvable in V 1,0 2 (Q T ). Proposition 1. For any pair (u, w) ∈ L 2 (Q T ) × L 2 (S T ) and any T > 0 the direct boundary-value problem (1) is uniquely solvable in V 1,0 2 (Q T )
and the estimate

|z δ [π]| Q T + z δ [π] 2,S T ≤ C T ( u 2,Q T + v δ 0 2,Ω + w 2,S T ) ,
takes place, where the constant C T is independent of δ ≥ 0 and pair π ≡ (u,

w) ∈ L 2 (Q T ) × L 2 (S T ). Also the adjoint problem -η t - ∂ ∂x j a i,j (x, t)η xi + a δ (x, t)η = χ(x, t) , η(x, T ) = ψ(x), x ∈ Ω, ∂η ∂N + σ δ (x, t)η = ω(x, t), (x, t) ∈ S T is uniquely solvable in V 1,0 2 (Q T ) for any χ ∈ L 2 (Q T ), ψ ∈ L 2 (Ω), ω ∈ L 2 (S T
) and any T > 0. Its solution is denoted as η[χ, ψ, ω]. Simultaneously, the estimate

|η δ [χ, ψ, ω]| Q T + η δ [χ, ψ, ω] 2,S T ≤ C 1 T ( χ 2,Q T + ψ 2,Ω + ω 2,S T ) ,
is true, where the constant C 1 T is independent of δ ≥ 0 and a triple (χ, ψ, ω). Simultaneously, from conditions a) -c) and the theorems on the existence of a weak (generalized) solution of the third boundary-value problem for a linear parabolic equation of the divergent type (see, e.g., [START_REF] Casas | Pontryagin's Principle for Local Solutions of Control Problems with Mixed Control-State Constraints[END_REF]), it follows that the direct boundary-value problem is uniquely solvable in

V 1,0 2 (Q T ) ∩ C(Q T ). Proposition 2. Let us l > n + 1. For any pair (u, w) ∈ L l (Q T ) × L l (S T ) and any T > 0, δ ∈ [0, δ 0 ] the direct boundary-value problem (1) is uniquely solvable in V 1,0 2 (Q T ) ∩ C(Q T ) and the estimate |z δ [π]| (0) Q T ≤ C T ( u l,Q T + |v δ 0 | (0) Ω + w l,S T ) ,
takes place, where the constant C T is independent of pair π ≡ (u, w) and δ.

Further, the minimization problem for Lagrange functional

L δ p,r (π, λ, µ) → min, π ∈ D when (λ, µ) ∈ L 2 (Q) × L + 2 (Q) (3)
plays the central role in all subsequent constructions. It is usual problem without equality and inequality constraints. It is solvable as a minimization problem for weakly semicontinuous functional on the weak compact set D ⊂ L 2 (Q T ) × L 2 (S T ). Here, the weak semicontinuity is a consequence of the convexity and continuity with respect to π of the Lagrange functional. Minimizers π δ [λ, µ] ∈ U δ [λ, µ] for this optimal control problem satisfy the Pontryagin maximum principle under supplementary assumption of the existence of Lebesgue measurable with respect to (x, t) ∈ Q for all z ∈ R 1 and continuous with respect to z for a.e.

x, t gradients

∇ z ϕ δ 2 (x, t, z), ∇ z G δ (x, z) with the estimates |∇ z ϕ δ 2 (x, t, z)| ≤ C M , |∇ z G δ (x, z)| ≤ C M ∀z ∈ S 1
M where C M > 0 is independent of δ. Due to the estimates of the propositions 1, 2 and to the so called two-parameter variation [START_REF] Sumin | The First Variation and Pontryagin's Maximum Principle in Optimal Control for Partial Differential Equations[END_REF] of the pair π δ [λ, µ] that is needle-shaped with respect to control u and classical with respect to control w the following lemma is true.

Lemma 1. Let H(y, η) ≡ -ηy and the additional condition that specified above is fulfilled. Any pair

π δ [λ, µ] = (u δ [λ, µ], w δ [λ, µ]) ∈ U δ [λ, µ], (λ, µ) ∈ L 2 (Q) × L + 2 (Q) satisfies to (usual) Pontryagin maximum principle in the problem (3): for π = π δ [λ, µ] the following maximum relations H(u(x, t), η δ (x, t)) = max u∈U H(u, η δ (x, t)) for a.e. Q T , (4)
H(w(s, t), η δ (s, t)) = max w∈W H(w, η δ (s, t)) for a.e. S T hold, where

η δ (x, t), (x, t) ∈ Q T is a solution for π = π δ [λ, µ] of the adjoint problem -η t - ∂ ∂x j (a i,j (x, t)η xi ) + a δ (x, t)η = ϕ δ 1 (x, t)λ(x, t) + ∇ z ϕ δ 2 (x, t, z δ [π](x, t))µ(x, t), (x, t) ∈ Q T , η(x, T ) = ∇ z G δ (x, z δ [π](x, T )), x ∈ Ω, ∂η(x, t) ∂N + σ δ (x, t)η = 0, (x, t) ∈ S T .
Remark 1. Note that here and below, if the functions ϕ δ 1 , ∇ z ϕ δ 2 (•, •, z(•, •)), λ, µ ∈ L 2 (Q) are considered on the entire cylinder Q T , we set that the equalities ϕ δ 1 (x, t) =∇ z ϕ δ 2 (x, t, z(x, t)) = λ(x, t) = µ(x, t) = 0 take place for (x, t) ∈ Q T \Q; the same notation is preserved if these functions are taken on the entire cylinder.

In the next section we construct minimizing approximate solutions for Problem (P 0 p,r ) from the elements

π δ [λ, µ], (λ, µ) ∈ L 2 (Q) × L + 2 (Q).
As consequence, this construction leads us to various versions of the stable sequential Lagrange principle and Pontragin maximum principle. In the case of strong convexity and subdifferentiability of the target functional g 0 0 , these versions are statements about stable approximations of the solutions of Problem (P 0 p,r ) in the metric of

Z ≡ L 2 (Q T ) × L 2 (S T ) by the points π δ [λ, µ].
Due to the estimates (2) and the propositions 1, 2 we may assert that the estimates

|g δ 0 (π)-g 0 0 (π)| ≤ C 1 δ ∀π ∈ D, g δ 1 (π)-g 0 1 (π) 2,Q ≤ C 2 δ(1+ π ) ∀π ∈ Z , (5) h δ -h 0 2,Q ≤ Cδ, g δ 2 (π) -g 0 2 (π) 2,Q ≤ C 3 δ ∀π ∈ D , hold, in which the constants C 1 , C 2 , C 3 > 0 are independent of δ ∈ (0, δ 0 ], π.

Stable Sequential Pontryagin Maximum Principle

In this section we discuss the so-called regularized or, in other words, stable, with respect errors of input data, sequential Pontryagin maximum principle for Problem (P 0 p,r ) as necessary and sufficient condition for elements of minimizing approximate solutions. Simultaneously, this condition we may treat as one for existence of a minimizing approximate solutions in Problem (P 0 p,r ) with perturbed input data or as condition of stable construction of a minimizing sequence in this problem. The proof of the necessity of this condition is based on the dual regularization method [START_REF] Sumin | A Regularized Gradient Dual Method for the Inverse Problem of a Final Observation for a Parabolic Equation[END_REF][START_REF] Sumin | Duality-Based Regularization in a Linear Convex Mathematical Programming Problem[END_REF][START_REF] Sumin | Regularized Parametric Kuhn-Tucker Theorem in a Hilbert Space[END_REF] that is stable algorithm of constructing a minimizing approximate solutions in Problem (P 0 p,r ). Sketches of the proofs for the theorems in this section (Theorems 1, 2, 3) and some comments may be found in [START_REF] Sumin | Stable Sequential Pontryagin Maximum Principle in Optimal Control for Distributed Systems[END_REF][START_REF] Sumin | Subdifferentiability of Value Functions and Regularization of Pontryagin Maximum Principle in Optimal Control for Distributed Systems[END_REF].

Dual Regularization for Optimal Control Problem with Pointwise State Constraints

The estimates [START_REF] Warga | Optimal Control of Differential and Functional Equations[END_REF] give a possibility to organize for constructing a minimizing approximate solution in Problem (P 0 p,r ) the procedure of the dual regularization in accordance with a scheme of the paper [START_REF] Sumin | Stable Sequential Pontryagin Maximum Principle in Optimal Control Problem with State Constraints[END_REF]. In accordance with this scheme the dual regularization consists in the direct solving dual problem to Problem (P 0 p,r ) and its Tikhonov stabilization

R δ,α (δ) p,r (λ, µ
) ≡ V δ p,r (λ, µ) -α(δ) (λ, µ) 2 → max, (λ, µ) ∈ L 2 (Q) × L + 2 (Q) under consistency condition δ/α(δ) → 0, α(δ) → 0, δ → 0. This dual regular- ization leads to constructing minimizing approximate solution in Problem (P 0 p,r ) from the elements π δ [λ δ,α(δ) p,r , µ δ,α(δ) p,r ] ∈ Argmin {L δ p,r (π, λ, µ)
: π ∈ D}, where (λ δ,α p,r , µ δ,α p,r ) ≡ argmax{R δ,α p,r (λ, µ) : (λ, µ) ∈ L 2 (Q) × L + 2 (Q)} and δ → 0. We may assert that the following "convergence" theorem for the dual regularization method in Problem (P 0 p,r ) is valid. Theorem 1. Regardless of the properties of the solvability of the dual problem to Problem (P 0 p,r ) or, in other words, regardless of the properties of the subdifferential ∂β(p, r) (it is empty or not empty), it is true that exist elements

π δ ∈ U δ [λ δ,α(δ) p,r , µ δ,α(δ) p,r
] such that the relations

g 0 0 (π δ ) → g 0 0 (π 0 p,r ), g 0 1 (π δ )-h 0 -p → 0, g 0 2 (π δ )-r ≤ κ(δ), κ(δ) → 0, δ → 0 , (λ δ,α(δ) p,r , µ δ,α(δ) p,r
), (g δ 1 (π δ ) -h δ -p, g δ 2 (π δ ) -r) → 0, δ → 0 hold, in which the inequality g 0 2 (π δ ) -r ≤ κ(δ) is understood in the sense of ordering on a cone of nonpositive functions in L 2 (Q). Simultaneously, the equality

lim δ→+0 V 0 p,r (λ δ,α(δ) p,r , µ δ,α(δ) p,r ) = sup (λ,µ)∈H×H+ V 0 p,r (λ, µ)
is valid. If the dual of Problem (P 0 p,r ) is solvable, then the limit relation (λ

δ,α (δ) p,r , µ δ,α(δ) p,r
) → (λ 0 p,r , µ 0 p,r ), δ → 0 is valid also, where (λ 0 p,r , µ 0 p,r ) denotes minimumnorm solution of the dual problem. This theorem may be proved in exact accordance with a scheme of proving the similar theorem in [START_REF] Sumin | Stable Sequential Pontryagin Maximum Principle in Optimal Control Problem with State Constraints[END_REF]. We note only that, as in [START_REF] Sumin | Stable Sequential Pontryagin Maximum Principle in Optimal Control Problem with State Constraints[END_REF], this proving uses a weak continuity of the operators g δ 1 , g δ 2 that is consequence of the conditions on the input data of Problem (P 0 p,r ) and a regularity of the bounded solutions of the boundary-value problem (1) inside of the cylinder Q T [8, ch.III, theorem 10.1].

Stable Sequential Lagrange Principle for Optimal Control Problem with Pointwise State Constraints

We formulate in this subsection the necessary and sufficient condition for existence of a minimizing approximate solution in Problem (P 0 p,r ). Also, it can be called by stable sequential Lagrange principle in nondifferential form for this problem. Simultaneously, as we deal only with regular Lagrange function, the formulated theorem may be called by Kuhn-Tucker theorem in nondifferential form. Note that the necessity of the conditions of formulated below theorem follows from the theorem 1. At the same time, their sufficiency is a simple consequence of the convexity of Problem (P 0 p,r ) and the conditions on its input data. A verification of these propositions for similar situation of the convex programming problem in a Hilbert space may be found in [START_REF] Sumin | Stable Sequential Convex Programming in a Hilbert Space and its Application for Solving Unstable Problems[END_REF][START_REF] Sumin | On the Stable Sequential Kuhn-Tucker Theorem and its Applications[END_REF].

Theorem 2. Regardless of the properties of the subdifferential ∂β(p, r) (it is empty or not empty) or, in other words, regardless of the properties of the solvability of the dual problem to Problem (P 0 p,r ), necessary and sufficient conditions for Problem (P 0 p,r ) to have a minimizing approximate solution is that there is a sequence of dual variables (λ k , µ k ) ∈ H × H + , k = 1, 2, . . ., such that δ k (λ k , µ k ) → 0, k → ∞, and relations

π δ k [λ k , µ k ] ∈ D δ k , k p,r , k → 0, k → ∞ , (6)
(λ k , µ k ), (g δ k 1 (π δ k [λ k , µ k ]) -h δ k -p, g δ k 2 (π δ k [λ k , µ k ]) -r) → 0, k → ∞ (7)
hold for some elements

π δ k [λ k , µ k ] ∈ U δ k [λ k , µ k ]. The sequence π δ k [λ k , µ k ], k = 1, 2, .
. ., is the desired minimizing approximate solution and each of its weak limit points is a solution of Problem (P 0 p,r ). As (λ k , µ k ) ∈ H × H + , k = 1, 2, . . ., we can use the sequence of the points (λ

δ k ,α(δ k ) p,r , µ δ k ,α(δ k ) p,r
), k = 1, 2, . . ., generated by the dual regularization method of the theorem 1. If the dual of Problem (P 0 p,r ) is solvable, the sequence (λ k , µ k ) ∈ H × H + , k = 1, 2, . . ., should be assumed to be bounded. The limit relation

V 0 p,r (λ k , µ k ) → sup (λ,µ)∈H×H+ V 0 p,r (λ, µ) (8)
holds as a consequence of the relations (6), [START_REF] Sumin | On the Stable Sequential Kuhn-Tucker Theorem and its Applications[END_REF]. Furthermore, each weak limit point (if such points exist) of the sequence

(λ k , µ k ) ∈ H × H + , k = 1, 2, . . . is a solution of the dual problem V 0 p,r (λ, µ) → max, (λ, µ) ∈ H × H + .
Remark 2. If the functional g 0 0 is strongly convex and subdifferentiable on D then from the weak convergence of the unique in this case elements π δ k [λ k , µ k ] to unique element π 0 p,r as k → ∞, and numerical convergence g 0

0 (π δ k [λ k , µ k ]) → g 0 0 (π 0 p,r ), k → ∞ follows the strong convergence π δ k [λ k , µ k ] → π 0 p,r , k → ∞. Problem (P 0 p,r
) with the strongly convex g 0 0 for linear system of ordinary differential equations but with exact input data is studied in [START_REF] Sumin | Parametric Dual Regularization for an Optimal Control Problem with Pointwise State Constraints[END_REF].

Stable Sequential Pontryagin Maximum Principle for Optimal Control Problem with Pointwise State Constraints

Denote by U δ max [λ, µ] a set of the elements π ∈ D that satisfy all relations of the maximum principle (4) of the lemma 1. Under the supplementary condition of existence of continuous with respect to z gradients ∇ z ϕ δ 2 (x, t, z), ∇ z G δ (x, z) with corresponding estimates, it follows that the proposition of the theorem 2 may be rewritten in the form of the stable sequential Pontryagin maximum principle. It is obviously that the equality U δ max [λ, µ] = U δ [λ, µ] takes place under mentioned supplementary condition. Theorem 3. Regardless of the properties of the subdifferential ∂β(p, r) (it is empty or not empty) or, in other words, regardless of the properties of the solvability of the dual problem to Problem (P 0 p,r ), necessary and sufficient conditions for Problem (P 0 p,r ) to have a minimizing approximate solution is that there is a sequence of dual variables (λ k , µ k ) ∈ H × H + , k = 1, 2, . . ., such that δ k (λ k , µ k ) → 0, k → ∞, and relations (6), [START_REF] Sumin | On the Stable Sequential Kuhn-Tucker Theorem and its Applications[END_REF] hold for some elements

π δ k [λ k , µ k ] ∈ U δ k max [λ k , µ k ].
Moreover, the sequence π δ k [λ k , µ k ], k = 1, 2, . . ., is the desired minimizing approximate solution and each of its weak limit points is a solution of Problem (P 0 p,r ). As (λ k , µ k ) ∈ H × H + , k = 1, 2, . . ., we can use the sequence of the points (λ δ k ,α(δ k ) p,r , µ δ k ,α(δ k ) p,r

), k = 1, 2, . . ., generated by the dual regularization method of the theorem 1. If the dual of Problem (P 0 p,r ) is solvable, the sequence (λ k , µ k ) ∈ H × H + , k = 1, 2, . . ., should be assumed to be bounded. The limit relation (8) holds as a consequence of the relations (6), [START_REF] Sumin | On the Stable Sequential Kuhn-Tucker Theorem and its Applications[END_REF]. Remark 3. When the inequality constraint in Problem (P 0 p,r ) is absent, i.e., (P 0 p,r ) = (P 0 p ), and ϕ 2 (x, t) = r ≡ 0, ϕ 1 (x, t) ≡ 1, the target functional g 0 0 is taken, for example, in the form g 0 0 (π) ≡ π 2 ≡ u 2 + w 2 then Problem (P 0 p ) acquires the typical form of unstable inverse problem. In this case the stable sequential Pontryagin maximum principle of the theorem 3 becomes a tool for the direct solving such unstable inverse problem. Remark 4. In important partial case of Problem (P 0 p,r ) = (P 0 r ), when it has only the inequality constraint (ϕ δ 1 (x, t) = h δ (x, t) = p(x, t) = 0, (x, t) ∈ Q), "weak" passage to the limit in the relations of the theorem 3 leads to usual for similar optimal control problems Pontryagin maximum principle (see, e.g., [START_REF] Casas | Pontryagin's Principle for Local Solutions of Control Problems with Mixed Control-State Constraints[END_REF][START_REF] Raymond | Pontryagin's Principle for State-Constrained Control Problems Governed by Parabolic Equations with Unbounded Controls[END_REF]) with nonnegative Radon measures in the input data of the adjoint equation.

Here and below, we use the notations for the sets QT , ST , Qi,T and also for functional spaces and norms of their elements adopted in monograph[START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF].

Acknowledgments. This work was supported by the Russian Foundation for Basic Research (project no. 15-47-02294-r -povolzh'e -a) and by the Ministry of Education and Science of the Russian Federation within the framework of project part of state tasks in 2014-2016 (code no. 1727).


Vasile Dragan
email: vasile.dragan@imar.ro

Adrian-Mihail Stoica

Toader Morozan
email: teoader.morozan@imar.ro

A Criterion for Robust Stability with Respect to Parametric Uncertainties Modeled by Multiplicative White Noise with Unknown Intensity, with Applications to Stability of Neural Networks

Keywords: Stochastic systems, Markovian jumps, multiplicative white noise, robust stability, Lyapunov operators

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The stochastic systems subject both to Markovian jumps and to multiplicative white noise perturbations received a considerable attention over the last years. Relevant results include the stability of such systems, optimal control and filtering (see e.g. [START_REF] Huang | Infinite horizon H2/H∞ control for stochastic systems with Markovian jumps[END_REF], [START_REF] Mariton | Jump Linear Systems in Automatic Control[END_REF], [START_REF] Li | Indefinite stochastic linear quadratic control with Markovian jumps in indefinite time horizon[END_REF], [START_REF] Dragan | Mathematical Methods in Robust Control of Linear Stochastic Systems Springer[END_REF] and their references). In the present paper a robust stabilization problem of continuous-time linear dynamic systems with Markov jumps and corrupted with multiplicative (state-dependent) white noise perturbations is considered. The robustness analysis is performed with respect to the intensity of the white noise terms. It is proved that the robustness radius depends on the solution of an algebraic system of coupled Lyapunov matrix equations. The derived results are a generalization of the ones proved in [START_REF] Levit | Algebraic criterion for stochastic stability of linear systems with parametric action of white noise type[END_REF] for the case without Markovian jumps. The paper is organized as follows: in the next section the problem statement is presented. The third section includes some preliminary results concerning the Lyapunov operators associated to the considered class of stochastic systems. The main result is presented and proved in Section 4. In the last section the stability radius is determined for two relevant particular cases and a numerical example illustrates the theoretical developments.

The problem statement

Consider the system of stochastic linear differential equations:

dx(t) = A(η t )x(t)dt + r l=1 µ l b l (η t )c T l (η t )x(t)dw l (t) (1)
where {w l (t)} t≥0 , 1 ≤ l ≤ r, are one-dimensional independent standard Wiener processes defined on a given probability space (Ω, F, P); {η t } t≥0 is a homogeneous standard right continuous Markov process defined on the same probability space (Ω, F, P) and taking value in the finite set N = {1, 2, ..., N } and having the transition semigroup P (t) = e Qt , t ≥ 0, where Q ∈ R N ×N is a matrix whose elements q ij satisfy the condition

   q ij ≥ 0 if i = j, i, j ∈ N N j=1 q ij = 0, ∀i ∈ N. (2)
For more details we refer to [START_REF] Chung | Markov chains with stationary transition probabilities[END_REF], [START_REF] Doob | Stochastic Processes[END_REF], [START_REF] Friedman | Stochastic Differential Equations and Applications[END_REF], [START_REF] Oksendal | Stochastic Differential Equations[END_REF]. We also assume that {η t } t≥0 , {w l (t)} t≥0 , 1 ≤ l ≤ r, are independent stochastic processes. In (1) the matrices

A(i) ∈ R n×n , b l (i), c l (i) ∈ R n×1 , 1 ≤ l ≤ r, 1 ≤ i ≤ N are
known, while the scalars µ l ∈ R are unknown. The system (1) can be regarded as a perturbation of the so called nominal system

ẋ(t) = A(η t )x(t). (3)
The perturbed system (1) emphasizes the fact that the coefficients of the nominal system are affected by parametric uncertainties modeled by state multiplicative white noise perturbations with unknown intensity µ l . Often when we refer to the perturbed system (1) we shall say that it corresponds to the vector of parameters µ = (µ 1 , µ 2 , ..., µ r ). Assuming that the nominal system (3) is exponentially stable in mean square (ESMS) we want to find necessary and sufficient conditions which will be satisfied by the parameters µ l , 1 ≤ l ≤ r such that the perturbed system (1) to be also ESMS. In the special case N = {1} (no Markov jumps) the conditions derived in this note recover those derived in [START_REF] Levit | Algebraic criterion for stochastic stability of linear systems with parametric action of white noise type[END_REF]. The concept of exponential stability in mean square of the linear stochastic systems of type ( 1) and (3) may be found in [START_REF] Costa | Continuous-time Markov jump linear systems[END_REF] and [START_REF] Dragan | Mathematical Methods in Robust Control of Linear Stochastic Systems Springer[END_REF], respectively.

Some preliminaries

Let S N n = S n ⊗S n ⊗...⊗S n , where S n ⊂ R n×n is the subspace of symmetric matrices. Let S N n+ be the convex cone defined by

S N n+ = {X = (X(1), X(2), ..., X(N )) ∈ S N n |X(i) ≥ 0, 1 ≤ i ≤ N }.
Here, X(i) ≥ 0 means that X(i) is positive semidefinite, S N n+ is a closed convex cone with non empty interior. Its interior is

IntS N n+ = {X = (X(1), X(2), ..., X(N )) ∈ S N n+ |X(i) > 0, 1 ≤ i ≤ N }.
A T (i)X(i) + X(i)A(i) + N j=1 q ij X(j) + r l=1 µ 2 l c l (i)b T l (i)X(i)b l (i)c T l (i) +H(i) = 0; 1 ≤ i ≤ N ; (4)
A(i)Y (i) + Y (i)A T (i) + N j=1 q ji Y (j) + r l=1 µ 2 l b l (i)c T l (i)Y (i)c l (i)b T l (i) +H(i) = 0, 1 ≤ i ≤ N ;
(5) (iv) There exists X = (X(1), ..., X(N )) ∈ IntS N n+ satisfying the following system of LMIs

A T (i)X(i) + X(i)A(i) + N j=1 q ij X(j) + r l=1 µ 2 l c l (i)b T l (i)X(i)b l (i)c T l (i) < 0,
where

1 ≤ i ≤ N .
Then one associates the following Lyapunov operators to the nominal system (3) L :

S N n → S N n , L : S N n → S N n defined by L[X](i) = A T (i)X(i) + X(i)A(i) + N j=1 q ij X(j), 1 ≤ i ≤ N, (6)
and

L[X](i) = A(i)X(i) + X(i)A T (i) + N j=1 q ji X(j), 1 ≤ i ≤ N (7)
for all X = (X(1), ..., X(N )) ∈ S N n . Is is easy to check that L is the adjoint operator of L with respect to the usual inner product on S N n :

< X, Y >= N j=1 T r[X(i)Y (i)]. (8)
Invoking Proposition 3.20 and Theorem 3.21 from [START_REF] Costa | Continuous-time Markov jump linear systems[END_REF] (see also Theorem 3.2.2 and Theorem 3.2.4 from [START_REF] Dragan | Mathematical Methods in Robust Control of Linear Stochastic Systems Springer[END_REF]) one deduces that the nominal system (3) is ESMS if and only if the eigenvalues of the linear operator L are located in the half plane C -. This allows us to obtain the following result.

Corollary 3.2. If the nominal system (3) is ESMS then for each H ∈ S N n the equations

L[X] + H = 0 (9)
and

L[Y] + H = 0 ( 10
)
have unique solutions given by X

= -L -1 [H] = (-L -1 [H](1), ..., -L -1 [H](N )) ∈ S N n and Y = -L -1 [H] = (-L -1 [H](1), ..., -L[H](N )) ∈ S N n , respectively. If H ∈ S N n+ then X ∈ S N n+ , Y ∈ S N n+ . Moreover, if H ∈ IntS N
n+ then the unique solutions of ( 9) and [START_REF] Levit | Algebraic criterion for stochastic stability of linear systems with parametric action of white noise type[END_REF], respectively are in IntS N n+ i.e.

-L -1 [H](i) > 0 (11)
and

-L -1 [H](i) > 0 ( 12
)
for all 1 ≤ i ≤ N . Further, let us consider the ordered space (R d , R d + )
where the order relation is induced by the convex cone

R d + = {x = (x 1 , x 2 , ..., x d ) T ∈ R d |x i ≥ 0, 1 ≤ i ≤ d}. The interior IntR d + of the convex cone R d + consists of the set IntR d + = {x = (x 1 , x 2 , ..., x d ) T ∈ R d + |x i > 0, 1 ≤ i ≤ d}. If D = (d ij ) ij ∈ R d×d is the matrix of a linear operator D : R d → R d then DR d + ⊂ R d + if and only if d ij ≥ 0, 1 ≤ i, j ≤ d.
In this case D will be called positive matrix. Applying Theorems 2.6 and 2.7 from [START_REF] Dragan | Mathematical Methods in Robust Control of Discrete-time Linear Stochastic Systems[END_REF] in the special case of the ordered linear space (R d , R d + ) one obtains the following result.

Proposition 3.3. For a positive matrix D ∈ R d×d the following are equivalent:

(i) ρ(D) < 1, ρ(•) being the spectral radius;

(ii) There exists ψ ∈ IntR d + such that the equation

(I d -D)ζ = ψ (13) has a solution ζ ∈ IntR d + .

Main results

The equivalence (i) ↔ (iv) from Proposition 3.1 allows us to deduce that if the perturbed system (1) is ESMS for a value µ = (µ 1 , µ 2 , ..., µ r ) of the intensities of the white noises, then this system is ESMS for every value

µ ' = (µ ' 1 , ..., µ ' r ) which are satisfying |µ ' l | ≤ |µ l | for all 1 ≤ l ≤ r.
In this section we shall derive a set of necessary and sufficient conditions which guarantee the exponential stability in mean square of a perturbed system (1) corresponding to a set of unknown vector of intensities µ = (µ 1 , ..., µ r ). Using Proposition 3.1 one notices that the system (1) is ESMS if and only if for any H ∈ IntS N n+ the equation ( 4) has a solution X ∈ IntS N n+ . Using (6) we may rewrite (4) in the form

L[X] + H = 0 ( 14
)
where H = ( H(1), ..., H(N )),

H(i) = r l=1 µ 2 l b T l (i)X(i)b l (i)c l (i)c T l (i) + H(i), 1 ≤ i ≤ N. (15)
Further, we rewrite H in the form:

H = r l=1 N j=1 µ 2 l b T l (j)X(j)b l (j)Ξ lj + H ( 16
)
where

Ξ lj = (Ξ lj (1), ..., Ξ lj (N )) ∈ S N n+ with Ξ lj (i) = c l (j)c T l (j), if i = j 0, otherwise . (17)
Since the nominal system (3) is necessarily ESMS if the perturbed system ( 1) is ESMS, we deduce via Corollary 3.2 and ( 16) that the solution of the equation ( 14) satisfies

X = - r l=1 N j=1 b T l (j)X(j)b l (j)µ 2 l L -1 [Ξ lj ] -L -1 [H].
The i th component of this solution is

X(i) = - r l=1 N j=1 b T l (j)X(j)b l (j)µ 2 l L -1 [Ξ lj ](i) -L -1 [H](i), 1 ≤ i ≤ N. ( 18
)
Multiplying (18) on the left by b T k (i) and on the right by b k (i), one obtains

b T k (i)X(i)b k (i) = r l=1 N j=1 b T l (j)X(j)b l (j)m ki,lj + ν ki (19) 1 ≤ k ≤ r, 1 ≤ i ≤ N , where m ki,lj = -µ 2 l b T k (i)L -1 [Ξ lj ](i)b k (i) (20)
and

ν ki = -b T k (i)L -1 [H](i)b k (i). ( 21
)
One sees that (19) is a system of rN scalars equations with rN scalar unknowns.

Based on the fact that for each integer α ∈ {1, 2, ..., rN } there exists a unique pair of natural numbers (k, i) ∈ {1, 2, ..., r} × {1, 2, ..., N } such that α = (k -1)N + i we may write (19) as an equation of the form (13) on the space R rN . To this end we set

ζ = (ζ 1 , ζ 2 , ..., ζ rN ) T , Ψ = (Ψ 1 , Ψ 2 , ..., Ψ rN ) T , D = (d αβ ) 1≤α,β≤rN , ζ β = b T l (j)X(j)b l (j) if (l -1)N + j = β, (22)
Ψ α = ν ki = -b T k (i)L -1 [H](i)b k (i) if (k -1)N + i = α (23)
d αβ = m ki,lj = -µ 2 l b T k (i)L -1 [Ξ lj ](i)b k (i) (24)
if (k -1)N + i = α and (l -1)N + j = β. With these notations (19) may be written in a compact form:

(I -D)ζ = Ψ. ( 25
)
Since the matrix D defined by (24) depend upon the unknown parameters µ l , for each perturbed system of type (1) one may associate a matrix D = D(µ) as before. Now we are in a position to state and proof the following result.

Theorem 4.1 Assume b k (i) = 0, 1 ≤ k ≤ r, 1 ≤ i ≤ N .
Under this condition the following are equivalent:

(i) The perturbed system (1) corresponding to the set of parameters µ = (µ 1 , ..., µ r ) is ESMS;

(ii) The nominal system (3) is ESMS and the matrix D(µ) associated to the perturbed system (1) satisfies ρ(D(µ)) < 1.

Proof (i) ⇒ (ii) Let H ∈ IntS N n+ be arbitrary but fixed. If the perturbed system (1) is ESMS, then the nominal system (3) is also ESMS. Hence, the equation ( 14)-(15) has a solution X = (X(1), ..., X(N )) with X(i) > 0, 1 ≤ i ≤ N . Consider ζ, Ψ ∈ R rN defined via ( 22) and (23). Since b k (i) = 0, for all (k, i) ∈ {1, 2, ..., r} × {1, 2, ..., N } one deduces that ζ α > 0 and Ψ α > 0, 1 ≤ α ≤ rN . Further, one associates the matrix D whose elements are computed via (24). One may check that d αβ ≥ 0, for all 1 ≤ α, β ≤ rN . So, it follows that the equation (25) associated to the perturbed system (1) satisfies the conditions from Proposition 3.3 (ii). Hence, ρ(D) < 1, this shows that the assertion (ii) from the statement holds if (i) is satisfied. Now the implication (ii) ⇒ (i) will be proved. Let H ∈ IntS N n+ be arbitrary. It will be shown that the corresponding equation ( 4) has a solution X ∈ IntS N n+ . One notices that if the nominal system (3) is ESMS, then Ψ ∈ IntR rN and the matrix D ∈ R rN ×rN are well defined via ( 23) and ( 24

µ 2 l b T k (i)L -1 [Ξ lj ](i)b k (i) ζlj -b T k (i)L -1 [H](i)b k (i) (27) ∀(k, i) ∈ {1, 2, ..., r} × {1, 2, ..., N }. Defining X(i) = - r l=1 N j=1 µ 2 l ζlj L -1 [Ξ lj ](i) -L -1 [H](i), 1 ≤ i ≤ N. ( 28
)
and setting X = (X(1), ..., X(N )), it results that

X ∈ IntS N n+ since -L -1 [Ξ lj ](i) ≥ 0, -L -1 [H](i) > 0 and ζlj > 0 for all (l, j, i) ∈ {1, 2, ..., r} × {1, 2, ..., N } × {1, 2, ..., N }. From (28) it follows L[X] + r l=1 N j=1 µ 2 l ζlj Ξ lj + H = 0. ( 29
)
Based on (17) it results that Ξ lj (i) = 0 if i = j, obtaining thus the following componentwise version of ( 29)

L[X](i) + r l=1 µ 2 l ζli c l (i)c T l (i) + H(i) = 0, 1 ≤ i ≤ N. (30)
From ( 28) with (27

) it results that ζki = b T k (i)X(i)b k (i) for all 1 ≤ k ≤ r, 1 ≤ i ≤ N . Therefore one may rewrite (30) in the form L[X](i) + r l=1 µ 2 l c l (i)b T l (i)X(i)b l (i)c T l (i) + H(i) = 0 (31)
which is just (4). Thus the proof is complete. Remark 4.1 The result proved in Theorem 4.1 shows that, in order to decide if the perturbed system (1) corresponding to a vector of unknown parameters

µ l ∈ [-|μ l |, |μ l |], 1 ≤ l ≤ r
is ESMS, we have to check if the spectral radius of the matrix D = D(μ) (associated via (24) to the parameters μl ) is less then 1. One notices that (24) may be rewritten in the form

d αβ = μ2 l b T k (i)Z lj (i)b k (i) (32)
where (k, i), (l, j) ∈ {1, 2, ..., r} × {1, 2, ..., N } are such that (k -1)N + i = α and (l -1)N + j = β, Z lj = (Z lj (1), ..., Z lj (N )) being the unique solution of the equation

A T (i)Z lj (i) + Z lj (i)A(i) + N ι=1 q iι Z lj (ι) + Ξ lj (i) = 0 (33) 1 ≤ i ≤ N , with Ξ lj defined in (17).
According with [START_REF] Levit | Algebraic criterion for stochastic stability of linear systems with parametric action of white noise type[END_REF] one introduces the following definition. Definition 4.1 The vector of noise intensities µ 0 = (µ 0 1 , µ 0 2 , ..., µ 0 r ) is called critical for the system (1) if the system (1) corresponding to the noise intensities εµ 0 is ESMS if 0 < ε < 1 and it is not ESMS if ε ≥ 1.

Remark 4.2 Denote D(µ) the matrix D corresponding to the vector µ = (µ 1 , ..., µ r ) of the noise intensities. Based on (24) it follows that D(εµ 0 ) = ε 2 D(µ 0 ). Since the spectral radius of a positive matrix is an eigenvalue of that matrix one may infer that

ρ[D(εµ 0 )] = ε 2 ρ[D(µ 0 )].
If µ 0 is a critical vector of noise intensities then one obtains from Theorem 4.1 that ρ[D(εµ 0 )] < 1 for all 0 < ε < 1 obtaining thus that ρ[D(µ 0 )] < 1 ε 2 , 0 < ε < 1. So, we deduce that ρ[D(µ 0 )] ≤ 1. Since the perturbed system (1) corresponding to the vector µ 0 is not ESMS, one concludes via Theorem 4.1, that ρ[D(µ 0 )] = 1. Hence, the vector µ 0 = (µ 0 1 , ..., µ 0 r ) is a solution of the equation det(I rN -D(µ)) = 0.

In the space R r of the vector µ = (µ 1 , ..., µ r ) the critical vectors of noise intensities µ 0 are included in the boundary of the stability region.

Several special cases

The first special case analyzed here is r = 1 and N ≥ 2. In this case the system (1) becomes dx(t) = A(η t )x(t)dt + µb(η t )c T (η t )x(t)dw 1 (t).

(34)

The matrix D associated to the system (34) is D = µ 2 D 1 where

D 1 =     b T (1)Z 1 (1)b(1) b T (1)Z 2 (1)b(1) ... b T (1)Z N (1)b(1) b T (2)Z 1 (2)b(2) b T (2)Z 2 (2)b(2) ... b T (2)Z N (2)b(2) ... ... ... ... b T (N )Z 1 (N )b(N ) b T (N )Z 2 (N )b(N ) ... b T (N )Z N (N )b(N )     (35)
for each 1 ≤ j ≤ N , (Z j (1), ..., Z j (N )) is the unique solution of the following equation on S

N n A T (i)Z j (i) + Z j (i)A(i) + N ι=1 q iι Z j (ι) + Ξ j (i) = 0 (36) Ξ j (i) = 0 if i = j and Ξ j (i) = c(j)c T (j) if i = j.
In the special case of the system (34), the Theorem 4.1 yields to the following result.

Corollary 5.1 If b(i) = 0, ∀ 1 ≤ i ≤ N the following are equivalent:

(i) The perturbed system (34) is ESMS;

(ii) The nominal system (3) is ESMS and the parameter µ satisfies the condition µ 2 < 1 ρ[D1] . Remark 5.1 The previous Corollary shows that the exponential stability in mean square of the nominal system (3) is preserved for the perturbed system (34) if and only if the unknown parameter µ lies in the interval (-ρ

-1 2 [D 1 ], ρ -1 2 [D 1 ]
), which is the stability region in the case of perturbed system (34).

The second special case discussed here is r ≥ 1, N = 1. Now, the system (1) becomes

dx(t) = Ax(t)dt + r l=1 µ l b l c T l x(t)dw l (t). (37)
The matrix D associated to the system (37) via ( 24) is

D = Ddiag(µ 2 1 , ..., µ 2 r ), (38) where D
=     b T 1 Z 1 b 1 b T 1 Z 2 b 1 ... b T 1 Z r b 1 b T 2 Z 1 b 2 b T 2 Z 2 b 2 ... b T 2 Z r b 2 ... ... ... ... b T r Z 1 b r b T r Z 2 b r ... b T r Z r b r     (39)
for each 1 ≤ l ≤ r, Z l is the unique solution of the Lyapunov equation

A T Z l + Z l A + c l c T l = 0. (40)
The result proved in Theorem 4.1 yields the next result. Remark 5.2 a) The result stated in Corollary 5.2 is just the main result derived in [START_REF] Levit | Algebraic criterion for stochastic stability of linear systems with parametric action of white noise type[END_REF]. Its discrete-time version may be found in [START_REF] Morozan | Necessary and sufficient conditions of stability of stochastic discrete systems[END_REF].

b) Condition of the form b l = 0, 1 ≤ l ≤ r (as it is imposed in the general case in Theorem 4.1) is redundant in the case of system (37) because, if b l0 = 0 for some l 0 , it follows that the noise w l0 (t) does not affect the perturbed system.

In order to illustrate the above theoretical results one considers the dynamics of Hopfield neural network of form vi (t) = a i v i (t) + n j=1 b ij g j (v j (t)) + c i , i = 1, ..., N a i < 0 and the activation functions g i (•) are strictly increasing.Then its approximation around an equilibrium point v 0 is ẋ(t) = Ax(t) + Bf (x(t)) where x(t) = v(t) -v 0 , f (x) = g(x + v 0 ) -g(v 0 ) and where A = diag(a 1 , ..., a n ) the elements of f (•) being sector-type nonlinearities satisfying f k (x k ) (f k (x k ) -µ k x k ) ≤ 0, k = 1, ..., n. Then for the above system associate the linear approximation dx(t) = Ax(t)dt + µ

(

  iii) For any H ∈ IntS N n+ , there exists Y = (Y (1), y(2), ..., Y (N )) ∈ IntS N n+ solving the following equation on S N n

  ), respectively. If ρ(D) < 1 then, based on the implication (i) ⇒ (ii) from Proposition 3.3 it follows (25) has a unique solution ζ ∈ IntR rN + . Based on the components of the vectors Ψ , ζ and of the matrix D one may define Ψki = Ψ α , ζlj = ζ β , dki,lj = d αβ (26) if (k, i), (l, j) ∈ {1, 2, ..., r} × {1, 2, ..., N }, (k -1)N + i = α, (l -1)N + j = β. It is easy to see that Ψki = Ψ ki , dki,lj = m ki,lj . With these notations one obtains the following version of the equation (25)

Corollary 5 . 2

 52 The following are equivalent (i) The perturbed system (37) is ESMS for any value of the unknown parameters µ l ∈ (-|μ l |; |μ l |), 1 ≤ l ≤ r.(ii) A is a Hurwitz matrix and ρ[D] < 1, D being defined in (38)-(39) with µ l replaced by μl , 1 ≤ l ≤ r.

b

  l c T l x(t)dw l (t), n = 3, r = 2 A = diag (-0.5, -0.5, -0.5) b 1 = [ 0.5, 1, 1 ] T , c 1 = [ 1, 2, 1 ] T , b 2 = [ 1, 0.25, -1 ] T , c 2 = [ 0.25, -1 , 1 ] TApplying Corollary 5.2 it results that the above system is ESMS for all µ ∈ [-0.2857; 0.2857].

Acknowledgement

This paper has been partially supported by MEN-UEFISCDI, Program Partnerships, Projects PN-II-PT-PCCA-2013-4-1349 and PN-II-PT-PCCA-2011-3.1-1560.