Elena K Kostousova
email: kek@imm.uran.ru

State Estimation for Control Systems with a Multiplicative Uncertainty through Polyhedral Techniques

Keywords: Differential systems, reachable sets, set-valued state estimation, multiplicative uncertainty, polyhedral estimates, parallelepipeds, parallelotopes, interval analysis

The paper deals with polyhedral estimates for reachable tubes of differential systems with a multiplicative uncertainty, namely linear systems with set-valued uncertainties in initial states, additive inputs and coefficients of the system. We present nonlinear parametrized systems of ordinary differential equations (ODE) which describe the evolution of the parallelotope-valued estimates for reachable sets (time crosssections of the reachable tubes). The main results are obtained for internal estimates. In fact, a whole family of the internal estimates is introduced. The properties of the obtained ODE systems (such as existence and uniqueness of solutions, nondegeneracy of estimates) are investigated. Using some optimization procedure we also obtain a differential inclusion which provides nondegenerate internal estimates. Examples of numerically constructed external and internal estimates are presented.

Introduction

The problem of constructing trajectory tubes (in particular, reachable tubes) is an essential theme in control theory [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF]. Since practical construction of such tubes may be cumbersome, different numerical methods are devised for this cause, in particular, methods based on approximations of sets either by arbitrary polytopes with a large number of vertices or by unions of points [START_REF] Bushenkov | Multi-dimensional Images Given by Mappings: Construction and Visualization[END_REF], [START_REF] Baier | A Computational Method for Non-convex Reachable Sets Using Optimal Control[END_REF], [START_REF] Artstein | Feedback and Invariance under Uncertainty via Set-Iterates[END_REF] (here and below, we mention, as examples, only some references from numerous publications; see also references therein). Such methods, as well as the methods based on different schemes of discrete approximations of initial set-valued problems [START_REF] Baier | Approximation of Reachable Sets by Direct Solution Methods of Optimal Control Problems[END_REF], [START_REF] Veliov | Second Order Discrete Approximations to Strongly Convex Differential Inclusions[END_REF] and numerical methods of solving the Hamilton-Jacobi-Bellman equation [START_REF] Taras'yev | Approximation Schemas and Finite-Difference Operators for Constructing Generalized Solutions of Hamilton-Jacobi Equations[END_REF], are devised to obtain approximations as accurate as possible. But they may require much calculations, especially for large dimensional systems; also smaller step-sizes create a heavy computational load. It is appropriate to mention approximations by polytopes based on support functions or supporting points [START_REF] Baier | Computing Aumanns Integral[END_REF], [START_REF] Krastanov | Dynamic Interactive System for Analysis of Linear Differential Inclusions[END_REF].

Other techniques are based on estimates of sets by domains of some fixed shape such as ellipsoids, parallelepipeds, zonotopes [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF], [START_REF] Filippova | Trajectory Tubes of Nonlinear Differential Inclusions and State Estimation Problems[END_REF], [START_REF] Gusev | Estimates of Reachable Sets of Multidimensional Control Systems with Nonlinear Interconnections[END_REF], [START_REF] Kornoushenko | Interval Coordinatewise Estimates for the Set of Accessible States of a Linear Stationary System. I-IV[END_REF]- [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF], [START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF], [START_REF] Kuntsevich | Calculation and Control of Attainability Sets for Linear and Certain Classes of Nonlinear Discrete Systems[END_REF]- [START_REF] Polyak | Ellipsoidal Parameter or State Estimation under Model Uncertainty[END_REF]. Fair results in this area were obtained for linear systems with set-valued initial states and set-valued additive uncertain inputs. The main advantage of the mentioned techniques is that they enable to obtain approximate solutions using relatively simple tools (up to explicit formulas). Note that more accurate approximations and even exact representations of the solutions may be obtained by using the whole families of such simple estimates (as it was proposed by A.B. Kurzhanski) [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF]- [START_REF] Kurzhanski | On Ellipsoidal Techniques for Reachability Analysis. Part II: Internal Approximations. Box-valued Constraints[END_REF], [START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF], [START_REF] Kostousova | External and Internal Estimation of Attainability Domains by Means of Parallelotopes[END_REF]. The methods of interval analysis which use subpavings of interval vectors [START_REF] Jaulin | Applied Interval Analysis[END_REF] serve the same purpose, but such methods may require much computations and memory for large dimensional systems.

It is also important to study linear systems when system matrices are uncertain too. This leads to the multiplicative uncertainty and additional difficulties due to nonlinearity of the problem (in particular, reachable sets -cross-sections of reachable tubes -can be non-convex). There are some results for such systems with different types of bounds on uncertainties [START_REF] Barmish | The Propagation of Parametric Uncertainty via Polytopes[END_REF], [START_REF] Digailova | On the Joint Estimation of the Model and State of an Under-Determined System from the Results of Observations[END_REF], [START_REF] Filippova | On the Estimation of Trajectory Tubes of Differential Inclusions[END_REF], including constructing external ellipsoidal estimates [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF], [START_REF] Polyak | Ellipsoidal Parameter or State Estimation under Model Uncertainty[END_REF] and external interval (in other terms, coordinate-wise or box-valued) estimates [START_REF] Kornoushenko | Interval Coordinatewise Estimates for the Set of Accessible States of a Linear Stationary System. I-IV[END_REF], [START_REF] Kuntsevich | Calculation and Control of Attainability Sets for Linear and Certain Classes of Nonlinear Discrete Systems[END_REF], [START_REF] Nazin | Interval Parameter Estimation Under Model Uncertainty[END_REF].

We construct polyhedral (parallelepiped-valued and parallelotope-valued) estimates for reachable sets and reachable tubes of differential systems with parallelepiped-valued uncertainties in initial states and in additive uncertain inputs and with interval uncertainties in coefficients of the system. In contrast to interval analysis, faces of our estimates may be not parallel to coordinate planes. The main results are obtained for the internal estimates. Using constructions from [START_REF] Kostousova | On Polyhedral Estimates for Reachable Sets of Discrete-Time Systems with Bilinear Uncertainty[END_REF], [START_REF] Kostousova | On Polyhedral Estimates for Trajectory Tubes of Dynamical Discrete-Time Systems with Multiplicative Uncertainty[END_REF], we obtain nonlinear parametrized systems of ordinary differential equations (ODE) which describe the evolution of centers and matrices of the parallelotope-valued internal estimates for the reachable sets. So, in fact, the whole family of internal estimates is introduced (but, unfortunately, unlike the case of linear systems [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF]- [START_REF] Kurzhanski | On Ellipsoidal Techniques for Reachability Analysis. Part II: Internal Approximations. Box-valued Constraints[END_REF], [START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF], [START_REF] Kostousova | External and Internal Estimation of Attainability Domains by Means of Parallelotopes[END_REF], this family does not ensure exact representations of the reachable sets in general). The properties of the obtained ODE systems (such as existence and uniqueness of solutions for fixed values of parameters, nondegeneracy of estimates) are investigated. Using some optimization procedure we also obtain a differential inclusion which provides nondegenerate internal estimates. ODE for external estimates were obtained earlier [START_REF] Kostousova | Outer Polyhedral Estimates for Attainability Sets of Systems with Bilinear Uncertainty[END_REF]. Here we remind these results for completeness of the exposition. Results of numerical simulations are presented.

The following notation is used below: R n -the n-dimensional vector space; ⊤ -the transposition symbol;

x 1 = n i=1 |x i |, x 2 = (x ⊤ x) 1/2 , x ∞ = max 1≤i≤n |x i | -vector norms for x = (x 1 , x 2 , . . . , x n ) ⊤ ∈ R n ; e i =
(0, . . . , 0, 1, 0, . . . , 0) ⊤ -the unit vector oriented along the axis 0x i (the unit stands at iposition); e = (1, 1, . . . , 1) ⊤ ; R n×m -the space of real n×m-matrices A = {a j i } = {a j } (with columns a j ); I -the unit matrix; 0 -the zero matrix (vector); Abs A = {|a j i |} for A = {a j i }; diag π, diag {π i } -the diagonal matrix

A with a i i = π i (π i -the components of the vector π); det A -the determinant of A ∈ R n×n ; tr A = n i=1 a i i -the trace of A; A = max 1≤i≤n m j=1 |a j i |the matrix norm for A∈R n×m induced by the vector norm x ∞ ; the notation of the type k = 1, . . . , N is used instead of k = 1, 2, . . . , N for shortening.

Problem Formulation

Consider the following system (x ∈ R n is the state):

ẋ = A(t) x + w(t), t ∈ T = [0, θ]. (1)
Here the initial state x(0) = x 0 ∈ R n , the input (control/disturbance) w(t) ∈ R n (which is assumed to be a Lebesgue measurable function) and the measurable matrix function A(t) ∈ R n×n are unknown but subjected to given set-valued constraints

x 0 ∈ X 0 , w(t) ∈ R(t), a.e. t ∈ T, (2)
A(t) ∈ A(t) = {A ∈ R n×n | A(t) ≤ A ≤ A(t)}, a.e. t ∈ T, (3)
where X 0 , R(t) are given convex compact sets in R n , the set-valued map R(t) is continuous, the matrix functions A(t), A(t) are continuous. Matrix and vector inequalities (≤, <, ≥, >) here and below are understood componentwise. The interval constraints (3) can be rewritten in the form

A(t) ∈ A(t) = {A| Abs (A-Ã(t)) ≤ Â(t)}, Ã = (A+A)/2, Â = (A-A)/2. ( 4
)
Let X (t)=X (t, 0, X 0 ) be a reachable set of system (1)-( 3) at time t>0 that is the set of all points x∈R n , for each of which there exist x 0 , w(•), A(•) that satisfy (2)-( 3) and generate a solution x(•) of (1) that satisfies x(t)=x. The multivalued function X (t), t∈T , is known as a trajectory (or reachable) tube X (•).

We presume the sets X 0 , R(t) to be parallelotopes (then the sets X (t) are not obliged to be parallelotopes) and look for external and internal parallelepipedvalued or parallelotope-valued (shorter, polyhedral ) estimates P ± (t) for X (t).

By a parallelepiped P(p, P , π)⊂R n we mean a set such that P =P(p,

P , π) = {x ∈ R n | x=p+ n i=1 p i π i ξ i , ξ ∞ ≤1}, where p∈R n ; P ={p i }∈R n×n is such that det P =0, p i 2 =1 1 ; π ∈ R n , π ≥ 0.
It may be said that p determines the center of the parallelepiped, P -the orientation matrix, p i -the "directions" and π i -the values of its "semi-axes". We call a parallelepiped nondegenerate if π>0.

By a parallelotope P[p, P ] ⊂ R n we mean a set

P = P[p, P ] = {x ∈ R n | x = p+ P ζ, ζ ∞ ≤ 1}
, where p ∈ R n and the matrix P = {p i } ∈ R n×m , m ≤ n, may be singular. We call a parallelotope P nondegenerate, if m = n and det P = 0.

Each parallelepiped P(p, P , π) is a parallelotope P[p, P ] with P = P diag π; each nondegenerate parallelotope is a parallelepiped with P = P diag { pi -1 2 }, π i = pi 2 or, in a different way, with P = P , π = e, where e = (1, 1, . . . , 1) ⊤ .

We call P an external (internal) estimate for X ⊂ R n if P ⊇ X (P ⊆ X ).

Assumption 1. The set X 0 = P 0 = P[p 0 , P0 ] = P(p 0 , P 0 , π 0 ) is a parallelepiped, the sets R(t) = P[r(t), R(t)] are parallelotopes where R(t) ∈ R n×m , m ≤ n; r(•), R(•) and A(•), A(•) are continuous vector and matrix functions.

Problem 1. Find some external P + (t) and internal P -(t) polyhedral estimates 2 for reachable sets X (t):

P -(t) ⊆ X (t) ⊆ P + (t), t ∈ T .
3 Auxiliary Discrete Time Systems. Primary Estimates

Following arguments similar to [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF], [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF]Sec. 3.2] we obtain ODE for the estimates. The first step in this way is to construct estimates for reachable sets X [k] of auxiliary discrete time systems -the Euler approximations 3 of the initial system:

x[k] = A[k-1]x[k-1] + w[k-1], k=1, . . . , N ; x[0] ∈ P 0 ; w[k] ∈ R[k] = h N R(t k ); A[k] ∈ A[k] = {I + h N A | A ∈ A(t k )}, (5)
t k =kh N , h N =θN -1 . It is known that X [k] satisfy the relations X [k] = A[k-1]• X [k-1]+R[k-1], k=1, . . . , N , X [0] = P 0 , which involve two operations with sets -multiplying an interval matrix A = {A ∈ R n×n |A≤A≤A} on a set X ⊂ R n : A • X = {y ∈ R n | y = Ax, A ∈ A, x ∈ X } and the Minkowski sum [25, p.93].
In [START_REF] Kostousova | Outer Polyhedral Estimates for Attainability Sets of Systems with Bilinear Uncertainty[END_REF], [START_REF] Kostousova | On Polyhedral Estimates for Reachable Sets of Discrete-Time Systems with Bilinear Uncertainty[END_REF], [START_REF] Kostousova | On Polyhedral Estimates for Trajectory Tubes of Dynamical Discrete-Time Systems with Multiplicative Uncertainty[END_REF], the ways of constructing primary polyhedral estimates for A • P and P 1 + P 2 (where P, P 1 , P 2 are parallelepipeds or parallelotopes) are described; hence we have the corresponding recurrence relations for external and internal estimates P ± [k] for X [k]. Passing to the limit as N → ∞, we obtain the corresponding nonlinear ODE systems for parallelotopes/parallelepipeds P ± (t).

Internal Estimates

We come to the following ODE system for parallelotopes P -(t)=P[p -(t), P -(t)]:

dp - dt = Ã(t) p -+ r(t), p -(0) = p 0 ; (6)
d P - dt = Ã(t) P -+ diag ν(t, P -; J(t)) • B( P -) + R(t) Γ (t), P -(0) = P0 , ν i (t, P -; J) = âji i (t) • η ji (t, P -), i = 1, . . . , n, η(t, P -) = max{0, Abs p -(t) -(Abs P -)e}, B = diag β( P -) • P -, β i ( P -) = 1 / (e i ⊤ (Abs P -) e), i = 1, . . . , n, (7)
2 Our estimates will satisfy the upper and lower semigroup properties and the superreachability and subreachability properties similarly to [25, Lemmas 3.3.2, 3.3.4] and [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF]Remark 8.2] respectively; these properties are analogues to the semigroup property [25, p.9] for X (t). Also our set-valued estimates will be continuous. 3 In connection with the Euler approximations, the papers [START_REF] Dontchev | Error Estimates for Discretized Differential Inclusions[END_REF], [START_REF] Wolenski | The Exponential Formula for the Reachable Set of a Lipschitz Differential Inclusion[END_REF] may be mentioned which analyzed the numerical error of the set-valued method.

(the operation of maximum is understood componentwise). Here Γ (t) ∈ R m×n is an arbitrary Lebesgue measurable matrix function satisfying Γ (t) = max 1≤i≤m n j=1 |γ j i | ≤ 1, a.e. t ∈ T , and J = {j 1 , . . . , j n } is an arbitrary permutation of numbers {1, . . . , n} or even a measurable vector function J(•) with values J(t) being arbitrary permutations of numbers {1, . . . , n}. Let G and J be the sets of all such functions Γ (•) and J(•) respectively.

Later on it is useful to mark out the following case.

Assumption 2. Either R(t) are singletons (then the function w(•) ≡ r(•) may be assumed to be measurable) or

Γ (•) ∈ G is such that R(t)Γ (t) ≡ 0, t ∈ T .
Theorem 1. Let the above assumptions about the system (1), ( 2), (4) be satisfied, P 0 be a nondegenerate parallelotope (det P0 = 0) and J(•) ∈ J, Γ (•) ∈ G.

Then the system ( 6), ( 7) has a unique solution (p -(•), P -(•)) at least on some subinterval To prove P -(t) ⊆ X (t) we verify the subreachability property:

T 1 = [0, θ 1 ] ⊆ T ,
P -(t) ⊆ X (t, s, P -(s)), ∀s, t : 0 ≤ s ≤ t ≤ θ 1 . Fix t ∈ T 1 . If x * ∈ P -(t)
, then there exists ξ such that Abs ξ ≤ e and x * = p -(t) + P -(t) ξ. Consider x * = x * (t) as a function of t (when ξ is fixed). Evidently, x * (s) ∈ P -(s) also for arbitrary s ≤ t. It remains to check that it is possible to find functions A(τ ) = Ã(τ ) + ∆A(τ ) ∈ A(τ ) and w(τ ) ∈ R(τ ), τ ∈ [s, t], such that x * (τ ) will satisfy (1) for τ ∈ [s, t]. Differentiating x * (τ ) with the account ( 6), [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF] we have (the argument τ is omitted for short): ẋ * = Ãx * + w + diag ν B ξ = Ax * + w + q, where w = r + RΓ ξ ∈ R (because Γ ξ ∞ ≤ Γ ξ ∞ ≤ 1), A = Ã + ∆A, q = diag ν B ξ -∆A (p -+ Pξ). The desired equality q = 0 is achieved if we take ∆A in the form ∆A = diag α D, where D = {e j1 • • • e jn } ⊤ is a matrix corresponding to a permutation J={j 1 , . . . , j n }=J(τ ) of rows of the unit matrix, and components of the vector α = α(τ ) are calculated by formulas α i = 0, if |p - ji | ≤ e ji ⊤ (Abs P -)e (i.e. if ν i = 0), and α i = ν i e i ⊤ Bξ/(p - ji + e ji ⊤ Pξ) otherwise. Inequalities |α i | ≤ âji i , i = 1, . . . , n (which ensure A ∈ A) are obtained by obvious estimates.

⊓ ⊔

Note that Theorem 1 describes the whole family of estimates P -(•) where J(•) and Γ (•) are parameters.

Remark 1. Obviously, we have X (t) ⊇ X 0 (t) ≡ P 0-(t), t ∈ T , where X 0 (t) are reachable sets of the system (1) under assumptions x 0 ∈ P 0 , w(•) ≡ r(•) and A(•) ≡ Ã(•), and parallelotopes P 0-(t) are determined by ( 6), [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF] when ν ≡ 0, Γ ≡ 0. We call these parallelotopes P 0-(t) trivial internal estimates for X (t).

The following corollary compares internal estimates P -(t) for X (t) satisfying ( 6), [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF] with trivial internal estimates P 0-(t) for X (t) in the sense of volume. We like to remind that volume of a nondegenerate parallelotope P = P[p, P ] ⊂ R n is equal to vol P = 2 n | det P |.

Corollary 1. Under conditions of Theorem 1, we have vol P

-(t) = vol P 0-(t) • exp (ψ 1 (t) + ψ 2 (t)), t ∈ T 1 , where ψ 1 (t) = t 0 ν(τ, P -(τ ); J(τ )) ⊤ β( P -(τ ))dτ , ψ 2 (t) = t 0 tr (Ξ(τ, P -(τ ))Γ (τ ))dτ , Ξ(t, P -) = ( P -) -1 R(t).
Therefore under additional Assumption 2 we have: (i) vol P -(t) ≥ vol P 0-(t), and vol P -(t) > vol P 0-(t) iff ψ 1 (t) > 0;

(ii) if it is turned out that P -(t) ∋ 0 for all t ∈ T , then P -(t) ≡ P 0-(t), t ∈ T .

Proof. The expression for vol P -(t) follows from the equality det P -(t) = m 0 (t)• exp (ψ 1 (t) + ψ 2 (t)), where m 0 (t) = det P0 exp t 0 tr Ã(τ ) dτ , which, in turn, can be obtained similarly to [17, p. 293]. Namely, we use the change of variables P -(t) = Φ(t) P (t), where Φ satisfies Φ = Ã Φ, Φ(0) = I, and obtain the relation

d det P dt / det P = ν(t, P -(t); J(t)) ⊤ β( P -(t)) + tr (( P -) -1 R Γ ) (8)
on the base of the known relation d det P/dt = det P tr (P -1 Ṗ ) and ( 7); then the Liouville formula det Φ(t)= exp t 0 tr Ã(τ ) dτ is used. Assumption 2 yields ψ 2 (t) ≡ 0. Thus (i) is evident, (ii) is true because we have ψ 1 (t) ≡ 0 similarly to [20, Corollary 1].

⊓ ⊔

In general case (without Assumption 2) we can obtain some differential inclusions which determine internal estimates on the whole time interval T . Consider two ways to do that. Following the first way, fix J(•)∈J and Γ (•)∈G; if a denominator of some row of the matrix B vanishes, replace this row by a suitable set. The corresponding differential inclusion determines P -(t) on T (see Theorem 2 below), but there is no guarantee for P -(t) to be nondegenerate. The second way allows (using some considerations of "local" optimality of the estimate volume) to construct Γ (•) to ensure nondegenerate estimates on T (Theorem 3).

Following the first way, consider the matrix differential inclusion

d P - dt ∈ Ã(t) P -+ diag ν(t, P -; J(t)) • B( P -) + R(t) Γ (t), P -(0) = P0 , (9)
where ν(t, P -; J(t)) is defined in [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF], and B( P -) is the set of all matrices B( P -) such that each row e i ⊤ B (i=1, . . . , n) satisfy the following conditions e i ⊤ B = e i ⊤ P -/(e i ⊤ (Abs P -) e), if e i ⊤ (Abs P -) e = 0, arbitrary row such that e

i ⊤ B 1 ≤ 1, if e i ⊤ (Abs P -) e = 0. ( 10
)
Theorem 2. For arbitrary J(•) ∈ J, Γ (•) ∈ G, there exists a solution (p -(•), P -(•)) of the system (6), ( 9)-( 10) which is determined on the whole T , and all solutions of this system determine internal parallelotope-valued estimates P -(t) for X (t), t ∈ T .

Proof. Existence and extendability are obtained using [10, p. 66, Theorem 6].

The inclusions P -(t) ⊆ X (t) are proved similarly to Theorem 1.

⊓ ⊔

Following the second way under conditions of Theorem 1, assume, without loss of generality, that det P0 >0. The idea of local optimization arises from [START_REF] Digailova | On the Joint Estimation of the Model and State of an Under-Determined System from the Results of Observations[END_REF] and consists in finding the maximal possible velocity of increasing det P -(t) (therefore vol P -(t)) at time t, by the choice of the value Γ , when the value P -= P -(t) has already been found. Consider the set Γ (t, P -) of matrices Γ (t, P -), which are solutions to the following optimization problem: max{tr (Ξ(t, P -)Γ ) | Γ ∈ R m×n s.t. Γ ≤ 1}. This set Γ (t, P -) may be described in the following form:

Γ (t, P -) = {Γ (t, P -) = {γ i k (t, P -)} | γ i k (t, P -) = sign (ξ k i (t, P -)) l i k , k = 1, . . . , m, i = 1, . . . , n, L = {l i k } ∈ L}. (11)
Here

Ξ(t, P -) = {ξ k i (t, P -)} = ( P -) -1 R(t) ∈ R n×m , and L is a set of matrices L = {l i k } ∈ R m×n satisfying conditions l i k ≥ 0, k = 1, . . . , m, i = 1, . . . , n; l i k = 0 if i ∈ I k (t, P -), k=1, . . . , m, i=1, . . . , n; n i=1 l i k = 1, k=1, . . . , m; I k (t, P -) = Argmax {|ξ k i (t, P -)| | i=1, .
. . , n}, k=1, . . . , m, where sign z is equal to -1, 0, 1 for z < 0, z = 0, z > 0 respectively. Consider the matrix differential inclusion

d P - dt ∈ Ã(t) P -+diag ν(t, P -; J(t))•B( P -)+ R(t) Γ (t, P -), P -(0) = P0 , (12)
where ν(t, P -; J) and B( P -) are the same as in [START_REF] Chernousko | Ellipsoidal Bounds on Reachable Sets of Dynamical Systems with Matrices Subjected to Uncertain Perturbations[END_REF].

Theorem 3. Let the above conditions be satisfied and P 0 be a nondegenerate parallelepiped with det P0 > 0. Then, for each function J(•) ∈ J, there exists a solution (p -(•), P -(•)) of system (6), ( 11)-( 12), which is determined on the whole interval T , and all solutions of this system determine parallelotopes P -(t) which turn out to be internal nondegenerate parallelepiped-valued estimates for X (t), t ∈ T .

Proof. We use arguments similar to [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF]Theorem 5.2]. Existence and extendability are obtained using [10, p. 66, Theorem 6]. In particular, ( 8) is used to see that the function P = Φ -1 Pis such that det P (t) is a nondecreasing function and therefore P (t) can not leave the domain where det P (t) ≥ det P (0) > 0; consequently, solutions to [START_REF] Filippova | On the Estimation of Trajectory Tubes of Differential Inclusions[END_REF] are defined on the whole interval T and determine nondegenerate parallelepiped-valued estimates. The inclusions P -(t) ⊆ X (t) are proved similarly to Theorem 1. ⊓ ⊔ Remark 2. We can choose J(•) in ( 7), ( 9) and ( 12) in different ways, in particular as a constant. A simple way is also to apply a "local" optimization which arises from [START_REF] Digailova | On the Joint Estimation of the Model and State of an Under-Determined System from the Results of Observations[END_REF]. Fix a natural number N and introduce a grid T N of times τ k = kh N , k=0, . . . , N , h N = θN -1 . Let us, for each τ k ∈ T N , solve the optimization problem which is to maximize ν(τ k , P -; J) ⊤ β( P -) over all possible permutations J = {j 1 , . . . , j n } assuming that P -= P -(τ k ) has already been found. Then we can sequentially construct the piecewise constant function J(t) ≡ J(τ k ) ∈

Argmax J ν(τ k , P -(τ k ); J) ⊤ β( P -(τ k )), t ∈ [τ k , τ k+1
), k = 0, . . . , N -1, and find P -(•). Note that the described procedure is not obliged to give the estimates P -(t) with maximal volume even if N → ∞.

External Estimates

In [START_REF] Kostousova | Outer Polyhedral Estimates for Attainability Sets of Systems with Bilinear Uncertainty[END_REF], the ODE systems of two types were obtained for external estimates for X (t) in the form of parallelepipeds P + (t) = P(p + (t), P (t), π + (t)), where P (t) is a fixed matrix function. Let us restate here, for completeness of the exposition, the ODE system for the more accurate estimates of the type II:

dp + dt = Ṗ P -1 p + + P (Φ (+) -Φ (-) )/2 + r, p + (0) = p 0 ; ( 13
)
dπ + dt = (Φ (+) + Φ (-) )/2 + Abs (P -1 R) e, π + (0) = Abs (P (0) -1 P 0 ) π 0 , where Φ

(±) i = max ξ∈Ξ ± i
±P -1 ( Ã-Ṗ P -1 )x + Abs (P -1 ) Â Abs x i ,

x = p + + P diag π + ξ; Ξ ± i = {ξ | ξ∈E(P(0, I, e)), ξ i =±1}, i=1, . . . , n, (14)
the symbol E(P) denotes the set of all vertices of a parallelepiped P =P(p, P , π), namely the set of points of the form

x = p + n j=1 p j π j ζ j , ζ j ∈ {-1, 1}.
Theorem 4. Let Assumption 1 be satisfied and P (t) ∈ R n×n be an arbitrary continuously differentiable function such that det P (t) = 0, t ∈ T . Then the system ( 13),( 14) has a unique solution (p + (•), π + (•)) on T , and the parallelepipeds P + (t) = P(p + (t), P (t), π + (t)) are the external estimates for the reachable sets X (t) of the system (1), ( 2), (4): X (t) ⊆ P + (t), t ∈ T .

Proof. The existence, uniqueness and extendability of the solution follow from [10, pp. 7,8,10], the inclusions -from [18, Theorem 1]. ⊓ ⊔ Remark 3. In fact, Theorem 4 describes the whole family of estimates where P (•) is a parameter. Some heuristic ways of choosing P (•) were indicated in [START_REF] Kostousova | Outer Polyhedral Estimates for Attainability Sets of Systems with Bilinear Uncertainty[END_REF] (in particular, (i) find P (•) from relations Ṗ = Ã(t)P , P (0)=P 0 , or (ii) put P (t)≡I).

Examples

Consider some examples. The estimates were calculated using the Euler approximations [START_REF] Barmish | The Propagation of Parametric Uncertainty via Polytopes[END_REF] with N = 100 (in fact, the estimates for X [k] are presented in figures below). But it would be emphasized that different schemes of approximation can be used for solving the obtained differential systems and finding the estimates. by discrete analogous to Theorem 3 (similar to [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF]Example 6.1]). Fig. 1 (b) shows the initial set P 0 (dashed line), the external estimate for X [N ] (thin line) and four internal ones. Three of them correspond to "quasistationary" functions Γ (•) (similarly to [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF]Example 6.1]), the last-named (thick line) corresponds to Theorem 3. For comparison, the trivial internal estimate P 0-[N ] is shown too (dashed thick line); it is the "smallest" of the presented internal estimates. 2, where drawings are similar to Fig. 1. Since parallelotopes here are three-dimensional, we present their two-dimensional projections on coordinate plains. The reachable sets belong to the intersection of external estimates and contain the internal ones.

Example 1. Let à ≡ 0 1 -1.5 0 ,  ≡ 0 0 0.1 0 , R ≡ P( 0
It must be admitted that the proposed estimates may turn out to be rather conservative. But we can calculate them easily via integration of the ODE, and they can give useful information, while it is hard to calculate exact reachable sets. Improved external (possibly nonconvex) estimates in the form of the union of parallelepipeds can be constructed for systems with constant coefficients [START_REF] Kostousova | Outer Polyhedral Estimates for Attainability Sets of Systems with Bilinear Uncertainty[END_REF].

Fig. 1 .

 1 Fig. 1. External and internal estimates for X [•] (a) and X [N ] (b) in Example 1.

  1, 1) ⊤ , I, (0.2, 0.2, 0.2) ⊤ ), θ = 0.4. Such system may be interpreted as a simple ecological model of dynamics of a number of microorganisms which have 3 stages of development, provide division at the last stage and produce from 2 to 8 descendants[33, p. 112]. The additive control describes injecting a preparation to reduce the population. Estimates for X [•] and X [N ] are shown in Fig.

3 Fig. 2 .

 32 Fig. 2. Projections of external and internal estimates for X [•] and X [N ] in Example 2.

  where 0<θ 1 ≤θ, and we have det P -(t) =0, t ∈ T 1 . The corresponding nondegenerate parallelotopes P -(t)=P[p -(t), P -(t)], t ∈ T 1 , are internal estimates for the reachable sets X (t) of the system (1), (2), (4): P -(t)⊆X (t), t∈T 1 . Under Assumption 2, the subinterval T 1 coincide with T .Proof. Here and below we give mostly only sketches; more details can be found in[START_REF] Kostousova | On Polyhedral Estimates for Reachable Sets of Differential Systems with a Bilinear Uncertainty. (Russian; to apperar)[END_REF]. The existence, uniqueness and extendability of the solution are obtained using the known results[10, pp. 7,8,10]. In particular, under Assumption 2, we verify that P -(t) can not leave the domain where det P

-(t) ≥ δ for some δ > 0 and then use [10, p. 10, Theorem 4].

The normality condition p i

2=1 may be omitted to simplify formulas (it ensures the uniqueness of the representation of a nondegenerate parallelepiped).

Acknowledgments. The work was supported by the Program of the Presidium of the Russian Academy of Sciences No. 17 "Dynamic Systems and Control Theory" and by the Russian Foundation for Basic Research (grant 12-01-00043).


Florian Dietrich

Samuel Martin

Marc Jungers

Transient cluster formation in generalized Hegselmann-Krause opinion dynamics

We analyse the generalized Hegselmann-Krause model of opinion dynamics. The asymptotic state of such a system has been well studied in the literature, however the transient state is still poorly understood. Predicting which groups of agents will form clusters remains to be studied. We present sufficient conditions to detect cluster formation in the transient phase of the multi-agent system. We also give a procedure to know how much time a cluster stays consistent, i.e., before it merges with other agents in the system. Our criterion can be computed locally using variables obtained from the initial conditions. Finally our results are illustrated by a numerical example.

I. INTRODUCTION

The analysis and design of cooperative behaviours in networked dynamic systems has recently received a lot of attention. It finds application in technical fields such as cooperative robotics [START_REF] Cortes | Coverage control for mobile sensing networks[END_REF], mobile sensor networks [START_REF] Curtin | Progress toward autonomous ocean sampling networks[END_REF] or distributed algorithms [START_REF] Tsitsiklis | Problems in decentralized decision making and computation[END_REF] and takes inspiration from collective motion in nature [START_REF] Reynolds | Flocks, herds and schools: A distributed behavioral model[END_REF]- [START_REF] Cucker | On the mathematics of emergence[END_REF].

A more recent direction of research is to use multi-agent systems as a tool to model social networks. The consensus system has been widely applied to model opinion dynamics. Regarding the linear consensus system, several conditions for convergence towards global agreement have been developed [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF]- [START_REF] Blondel | Convergence in Multiagent Coordination, Consensus, and Flocking[END_REF].

By contrast, when the communication strength between agents depends on their current opinions, as in the Hegselmann-Krause bounded confidence model [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF], [START_REF] Deffuant | Mixing beliefs among interacting agents[END_REF], it has been shown that clusters may asymptotically form. Groups of agents reach local consensuses but disagree with the rest of the agents [START_REF] Blondel | On Krause's multi-agent consensus model with state-dependent connectivity[END_REF]- [START_REF] Bhattacharyya | On the convergence of the Hegselmann-Krause system[END_REF]. The study of cluster formation is of particular relevance in social systems since it is linked to the formation and preservation of local norms and cultures (see [START_REF] Mäs | Individualization as driving force of clustering phenomena in humans[END_REF] and references therein). Cluster patterns may also lead to conflict between distant clusters [START_REF] Durkheim | The division of labor in society[END_REF]. From a technical point of view, cluster formation in multi-agent systems allows to apply model reduction tools and simplify the analysis of the system [START_REF] Bıyık | Area aggregation and time-scale modeling for sparse nonlinear networks[END_REF].

In static networks of interactions, the detection of clusters or community patterns have been largely explored [START_REF] Morȃrescu | Opinion dynamics with decaying confidence: Application to community detection in graphs[END_REF], [START_REF] Newman | Modularity and community structure in networks[END_REF]. However, interactions and opinions evolve in time. Methods to predict the birth of opinion clusters remains to be developed. Moreover, results on the convergence to clusters of opinions do not provide a way to predict a priori which groups of agents will converge towards

The authors are with Université de Lorraine, CRAN, UMR 7039, 2 Avenue de la Forêt de Haye, Vandoeuvre-lès-Nancy, France and CNRS, CRAN, UMR 7039. {florian.dietrich, samuel.martin, marc.jungers} at univ-lorraine.fr a cluster. Finally, previous analyses of opinion dynamics always focus on the asymptotic behaviour of the systems, be it a global consensus or local clusters. However, simulations have shown that despite global consensus may occur asymptotically, meta-stable clusters may form in temporarily the transient phase [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF]. This is in particular the case for the generalized Hegselmann-Krause model where the communication range is not a bounded set [START_REF] Chatterjee | Generalized opinion dynamics from local optimization rules[END_REF] (see Section II-B for an illustration of this fact). It is of interest to focus on the transient state rather than the asymptotic state because in many real world applications, the time frame of interest (political election, war, etc.) is finite. Finite time analysis has analogously been reported in epidemiology in [START_REF] Barbillon | Network impact on persistence in a finite population dynamic diffusion model: application to an emergent seed exchange network[END_REF], where extinction is the stable asymptotic state, but the study rather focuses on quasi-equilibrium occurring between the emergence of a disease and its extinction.

In this paper, we provide a new criterion to detect the formation of opinion clusters in systems of opinion dynamics. The criterion depends only on the initial set of opinions and can be computed locally. The results apply to various forms of opinion-dependent communication strengths, possibly with unbounded support such as the generalized Hegselmann-Krause model. The criterion we offer allows to detect cluster formation in the transient dynamics. It also gives a duration on how much time we can differentiate a local cluster from the rest of the agents before global consensus occurs.

This paper is organized as follows. In Section II we formalise the multi-agent system we use in this paper by introducing notation and definitions. We also describe experimental observations that display transient cluster formation. In Section III, several assumptions are stated in order to obtain simple bounds on the dynamics of the system. These bounds are used in Section IV to give conditions to detect transient cluster formation. We then give a numerical illustration in Section V to show how cluster formation is detected using our main result. Finally, we conclude in Section VI.

Throughout the paper, the notation R will denote the set of real numbers, R + the set of the positive real (zero included) and N the set of integers. The notation |.| will denote the cardinality of a set of integers, and the absolute value for a real number.

II. PROBLEM FORMULATION

Let us consider a system of n agents, numbered from 1 to n and forming a discrete-time multi-agent system. The opinion or state of agent i at time t ∈ N is a scalar denoted x i (t) ∈ R. When there is no ambiguity, we drop the time notation and write x i = x i (t) and x + i = x i (t+1), the opinion at time t + 1, which is obtained according to the following consensus update :

x + i = x i + n j=1 α ij n k=1 α ik (x j -x i ), (1)
where α ij is the weight associated with the influence of agent j over i. In this paper we will use state dependent weights of the form

α ij = f (|x i -x j |), (2)
where the function f : R + → R + is called the influence function. Examples of suitable influence functions are provided in Section III-B.

A. Notation and Definitions

This subsection gives two definitions used to describe the multi-agent system. Namely, two key distances are defined for a given group of agents : the diameter of the group, and the distance separating the group from the other agents of the system. For the rest of the paper, let A ⊆ {1, . . . , n} be a group of agents. The group A will be the cluster candidate.

Definition 1 (Group diameter δ): The diameter of the group A is the maximum distance between two agents in

A δ = max (i,j)∈A 2 |x i -x j |. (3)
Definition 2: The distance ∆ is defined as the minimum distance between all agents in group A and all the other agents in the system, which can be formally written as

∆ = min i∈A, j / ∈A |x i -x j |. (4)
We just defined the main parameters of a group of agents, and we will then look at how these parameters evolve.

B. Experimental observations

As it has been emphasized in the introduction, the asymptotic behaviour of such a system is well known. In particular, if the influence function is bounded above and remains positive over R + , the system always asymptotically converges towards consensus [START_REF] Blondel | Convergence in Multiagent Coordination, Consensus, and Flocking[END_REF]. Nevertheless, to characterize the transient behaviour remains an open question. One typical observation is that before merging toward a unique opinion, agents converge into some temporary local agreement. Figure 1 displays an instance of opinion dynamics (1) when the influence function is the decreasing exponential f = x → exp(-x). It first shows that two clusters of similar opinions rapidly form, these clusters remain far appart over a relatively long period of time while they get closer. Finally the two groups merged into one to achieve the expected consensus.

We will now provide a definition to characterize this transient cluster formation. During the first phase when the agents form clusters, the cluster formation is faster than the dynamics which tends to bring the clusters together. Formally speaking and using the previous definition, we can foresee the evolution of the diameter δ of a group is decreasing faster than the distance ∆ between the two groups.

C. Problem statement

From the previous observations, the problem of the detection of transient clusters arises. Once a cluster is detected, it is also important to know for how long the cluster stays before merging with other agents.

The criterion we present to answer this question is based only on the evolution over time of the couple (δ, ∆) defined in equations ( 3) and ( 4). In accordance with the previous observation, we say that cluster is forming when a group of agents has a diameter decreasing faster than the distance separating the group from the other agents of the system. Definition 3 (Cluster formation): Given a group of agents A, this group is forming a cluster at time t if its ratio δ/∆ is decreasing from t to t + 1, i.e. if it verifies

δ(t + 1) ∆(t + 1) ≤ δ(t) ∆(t) . ( 5
)
We propose in the rest of this paper an upper bound of the evolution on the ratio δ/∆, by finding an upper bound on the evolution of δ (Section III-D), and then a lower bound on the evolution of the distance ∆ (Section III-E). As those bounds are obtained by considering worst case scenarios, we can then form a worst case system whose evolution bounds the evolution of the real system. By mean of this worst case system, which is easily computable locally, we can predict how long a group matches the definition of a cluster stated above (Theorem 1 in Section IV).

III. PRELIMINARIES

A. Assumptions

In this subsection, several assumptions are stated for use in later lemmas and propositions. We first introduce the function h defined as follows to lighten the equations :

h = R + → R + x → x f (x) , (6)
with f the influence function defined in [START_REF] Curtin | Progress toward autonomous ocean sampling networks[END_REF].

Assumption 1: Function f defined in (2) is non-negative, non-increasing, of class C 2 on R + and verifies f (0) > 0.
The non-increasing and non-negative properties of the influence function is common in terms of opinion dynamics because the further apart two opinions are, the less they tend to influence each others. The classical Hegselmann-Krause model from [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF] does not satisfy the previous assumption because of the discontinuity of its influence function, hence it is not C 2 . However we propose in this paper smoothed versions of the Hegselmann-Krause function to match the previous assumption.

The next assumption is more technical and allows us to get locally computable results, but it is not essential to guaranty the occurrence of transient cluster formations.

Assumption 2: Function h defined in [START_REF] Cucker | On the mathematics of emergence[END_REF] has the following properties :

-there exists x such that h is concave on [0, x] and convex on [ x, +∞), -lim x→+∞ h(x) is finite, -and for all x ∈ arg max x∈R + h(x), we have h ( x) < 0.

According to the previous assumption, x is an inflexion point of function h. Aiming at preserving the order of agents, we consider the following assumption. Assumption 3: Function f is log concave, with log the natural logarithm, meaning function log(f ) is concave. Assumption 3 ensures that the multi-agent system is order preserving, as shown in [29, Th 2]. It guarantees that agents are staying in the same order over time. More formally the order preservation property is stated in the following lemma.

Lemma 1: Suppose that Assumption 3 holds. The multi-agent system (1) is order preserving, i.e., for any agent i and j in the multi-agent system,

x i ≤ x j implies x + i ≤ x + j .
In the sequel of the paper, a renumbering allows to write without loss of generality

x 1 ≤ x 2 ≤ • • • ≤ x n .
It is also natural to consider for candidate clusters, groups of consecutive opinions, that is two integers m, M ∈ {1, . . . , n} such that m ≤ M and A = {m, • • • , M }. One immediate consequence is that under Assumption 3 we have for all time

δ = x M -x m . (7)
Examples of influence functions satisfying Assumptions 1-3 are detailed in the following subsection.

B. Examples of suitable influence functions

Here we give two examples of functions satisfying Assumptions 1-3. The class of decreasing exponentials satisfies Assumptions 1-3 and is a classical class of influence function [START_REF] Chatterjee | Generalized opinion dynamics from local optimization rules[END_REF].

Another class of functions satisfying Assumptions 1-3 is the one of smoothed generalized Hegselmann-Krause influence functions on infinite support. We choose to consider such functions using sigmoid according to the formula

f GHK (x) = 1 -sig(α(x -1)) 1 -sig(-α) , sig(x) = 1 1 + e -x
where α is a parameter to control the smoothness of the curve. Figure 2 shows that the larger the α the closer we are to an original Hegselmann-Krause function as in [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF], with a communication range of 1. Note that the Hegselmann-Krause model cannot be considered here due to its discontinuity, which does not verify Assumption 1.

C. Basic properties

This subsection gives properties of function h defined in [START_REF] Cucker | On the mathematics of emergence[END_REF] and will be used in the proofs of this paper.

Lemma 2 (Function h): Under Assumptions 1 and 2 the following properties hold for function h defined in (6) :

arg max x∈R + h(x) is a singleton and its unique element is denoted x, -h is non-decreasing on [0, x] and non-increasing on [ x, +∞).

Proof: Under Assumption 2, function h is convex on interval [ x, +∞). We then prove that h is non-increasing on the same interval. Assume the opposite is true. Then, there exist x, y ∈ [ x, +∞) with x < y such that h(x) < h(y). By the mean value theorem, there exits z ∈ [ x, +∞) such that h (z) > 0. Since h is convex over this interval, the graph of h stays above its tangents. As a consequence, function h would diverge if it were not non-increasing.

As h is non-increasing on [ x, +∞), it falls that all elements of arg max x∈R + h(x) are in [0, x]. Let us pick an element x ∈ arg max x∈R + h(x) and show that it is unique. By concavity of h over [0, x], h is non-increasing over the same interval. Since h ( x) = 0, h is non-negative over [0, x] and non-positive over [ x, x]. This implies that h is non-decreasing over [0, x] and non-increasing over [ x, x] and also over [ x, ∞).

Under Assumption 1, h is C 2 and its Taylor development up to the order 2 at x is

h(x) = h( x) + (x -x) 2 2 h ( x) + o(x -x) 2 ,
where the first order term is zero because, by definition,

x is an extremum of function h, so h ( x) = 0. Then by dividing the previous equation by the quantity (x -x) 2 , we get for

x = x h(x) -h( x) (x -x) 2 = h ( x) 2 + o(1)
.

Let α = -1 4 h ( x) > 0, by Assumption 2. From definition of o(1) it comes that it exists η > 0 such that for x ∈ [ x -η , x + η]\{ x} we have h(x) -h( x) (x -x) 2 ≤ h ( x) 2 -α = h ( x) 4 < 0,
according to Assumption 2, then as a consequence,

∀x ∈ [ x -η, x + η]\{ x}, h(x) < h( x). ( 8
)
Since h is non-decreasing over [0, x] and non-inscreasing over [ x, +∞), using equation ( 8), we obtain ∀x ∈ [0, x]\{ x}, h(x) < h( x).

D. Upper bound on group diameter evolution

In this subsection we give an upper bound on the evolution of the diameter δ defined in (3) of a selected group of agents A, over two iterations. This bound is defined as the minimum of two bounds obtained via different approaches, one being less conservative for small values of δ.

Proposition 1: Under Assumptions 1-3 and for ∆ ≥ x, we have

δ + ≤ µ(δ, ∆) δ , (9)
where

µ(δ, ∆) = 1 - |A| f (δ) s A (δ, ∆) + (n -|A|) s B (δ, ∆) µ B (δ, ∆), (10)
with

s A (δ, ∆) = f (δ) + (|A| -1)f (0) + (n -|A|) f (∆), (11)
s B = s B (δ, ∆) = f (0) + (|A| -1)f (δ) + f (∆), (12)
µ B (δ, ∆) = 1 δ h(∆)+min{M h(∆)-h(∆+δ), 0} , (13)
M = M (δ, ∆) =    ε 0 (s B -ε 0 ) if ε 0 ≤ s B 2 1 otherwise , (14)
ε 0 = ε 0 (δ, ∆) = δ |A| max [0,δ] |f | + (n -|A|) max [∆,+∞) |f | , (15)
with function h defined in [START_REF] Cucker | On the mathematics of emergence[END_REF].

Remark 1: The previous equations are easy to compute locally because the only variables encountered are cardinalities and distances δ and ∆.

This result can be interpreted as follows : if the computation of µ(δ, ∆) gives a value smaller than 1, it gives the information that the group of agents is contracting, meaning its diameter is decreasing.

Proof: First, let us detail the expression of δ + given in [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF], where we split the global sum into three, each corresponding to a part of the system : the selected group A, the agents below A put in a group named B 1 , and the agents above A grouped in B 2 . More formally we have

B 1 = {1, • • • , m -1} and B 2 = {M + 1, • • • , n}. Then we can write δ + = δ -d A + d B , (16)
where

d A = i∈A h(x M -x i ) σ M + h(x i -x m ) σ m , (17)
with the notation d B = d B1 + d B2 and

d B1 = i∈B1 h(x m -x i ) σ m - h(x M -x i ) σ M , (18)
d B2 = i∈B2 h(x i -x M ) σ M - h(x i -x m ) σ m , (19)
and

σ i = N j=1 f (|x i -x j |). (20)
To obtain an upper bound on δ + , we will seek separate bounds on previously defined quantities d A and d B .

It can be first verified that the quantities s A and s B respectively defined in [START_REF] Hendrickx | Convergence of type-symmetric and cut-balanced consensus seeking systems[END_REF] and ( 12) satisfy

s A (δ, ∆) ≥ max{σ m , σ M }, and s B (δ, ∆) ≤ min{σ m , σ M },
where σ i has been defined in [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: convergence conjectures and theorems[END_REF]. Such bounds deal with the denominators of d A and d B .

To obtain a lower bound on d A defined in [START_REF] Deffuant | Mixing beliefs among interacting agents[END_REF], we first remark that, for i ∈

A f (x M -x i ) ≥ f (δ) and f (x i -x m ) ≥ f (δ), leading to h(x M -x i ) ≥ (x M -x i )f (δ), and h(x i -x m ) ≥ (x i -x m )f (δ).
By summing the two previous equations it comes

h(x M -x i ) + h(x i -x m ) ≥ h(δ),
allowing us to find the following lower bound for d

A d A ≥ |A| h(δ) s A (δ, ∆) .
To cope with d B , we propose the following upper bound

d B ≤ (n -|A|) max u∈J h(u) s B (δ, ∆) -min u∈J h(u + δ) s A (δ, ∆) ,
with u an appropriate substitution and J = [∆, +∞). Thanks to Lemma 2, for ∆ ≥ x we have max u∈J h(u) = h(∆) and inf u∈J h(u + δ) = 0. This allows to get a part of the expression of µ(δ, ∆) in ( 10) by combining the two previous equations.

To obtain equation ( 14) defining M (δ, ∆), we rewrite the terms of the sum d B1 in [START_REF] Blondel | On Krause's multi-agent consensus model with state-dependent connectivity[END_REF], knowing the reasoning is the same for d B2 defined in [START_REF] Chatterjee | Generalized opinion dynamics from local optimization rules[END_REF]. We define, for an agent j ∈ B 1 , the quantity Q j as

Q j = h(x m -x j ) σ M + ε - h(x M -x j ) σ M ,
where ε = σ m -σ M , and the previous equation can be rewritten as

Q j = 1 σ M h(x m -x j ) -h(x M -x j ) + 1 σ M + ε - 1 σ M h(x m -x j ).
At this point and in order to bound d B , we want to obtain an upper bound on Q j for an agent j ∈ B 1 .

First, under the assumption ∆ ≥ x and according to Lemma 2, we have h(x m -x j ) ≤ h(∆) because with j ∈ B 1 , x m -x j ≥ ∆.

Next we suppose ∆ ≥ x, with x the inflexion point of function h defined in Assumption 2. Then for agent i ∈ B 1 we remark that x m -x j ≥ ∆ and x M -x j ≥ ∆ + δ. Then using Assumption 2, we also know that function h is convex on [ x, +∞), implying that its first derivative h is non-decreasing, allowing us to write

h (x m -x j ) -h (x M -x j ) ≤ 0.
This means that u → h (u) -h (u + δ) is non-increasing on [ x, +∞) and allows us to get the following upper bound

h(x m -x j ) -h(x M -x j ) ≤ h(∆) -h(∆ + δ).
Finally, we need to find an upper bound of

1 σ M +ε -1 σ M
and for that we need to bound the absolute value of ε, which is, by definition of σ m and σ M in [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: convergence conjectures and theorems[END_REF], given by

ε = n i=1 f (|x i -x m |) -f (|x i -x M |).
Then by using the triangle inequality and the mean value theorem, |ε| can be bounded by

|ε| ≤ δ |A| max [0,δ] |f | + (n -|A|) max [∆,+∞) |f | .
The right member of the previous equation is the definition of ε 0 (δ, ∆) given in [START_REF] Blondel | Convergence in Multiagent Coordination, Consensus, and Flocking[END_REF]. Finally, let us show that

1 σ M + ε - 1 σ M ≤ M (δ, ∆) s B (δ, ∆) ,
where M (δ, ∆) is defined in [START_REF] Morȃrescu | Coordination in networks of linear impulsive agents[END_REF]. It comes directly that

1 σ M + ε - 1 σ M ≤ 1 s B (δ, ∆)
.

Then we show that if ε 0 < s B we have

1 σ M + ε - 1 σ M ≤ ε 0 s B (s B -ε 0 )
.

For that we consider the case where ε < 0, the opposite case giving a null upper bound that does not add precision. So for ε < 0 we have

1 σ M + ε - 1 σ M = |ε| σ M (σ M -|ε|) ,
and then for

ε 0 < σ M such that ε ≤ |ε 0 | that |ε| σ M (σ M -|ε|) ≤ ε 0 σ M (σ M -ε 0 ) ≤ ε 0 s B (s B -ε 0 )
.

By observing that for s B ∈]ε 0 , 2ε 0 ] we have

1 s B ≤ ε 0 s B (s B -ε 0 ) ,
and we obtain the quantity M (δ, ∆) as in [START_REF] Morȃrescu | Coordination in networks of linear impulsive agents[END_REF] and with it the quantity µ(δ, ∆) in ( 10) by combination of the different results we just provided.

E. Lower bound on the evolution of inter-group distance

This subsection gives a lower bound on the evolution of the distance ∆ related to a selected group of agents, ∆ defined in (4) being the minimum distance separating the agents in the group from the others. Proposition 2: Under Assumptions 1-3 and for ∆ ≥ x we have

∆ + ≥ η(δ, ∆) ∆, (21)
with

η(δ, ∆) = 1 -|A| η A (δ, ∆) -(n -|A|) η B (δ, ∆), (22)
where

η A (δ, ∆) = h(∆) σ B (δ, ∆) , (23)
η B (δ, ∆) = 1 ∆ max h(∆) σ A (δ, ∆) , h(2∆ + δ) σ B (δ, ∆) , (24)
σ A (δ, ∆) = f (0) + (|A| -1) f (δ) + f (∆), (25)
σ B (δ, ∆) = f (0) + f (∆) + (|A| -1) f (δ + ∆). (26)
Remark 2: The previous equations are, just like the upper bound on δ + , easy to compute locally because only local variables (specific to the selected group of agents A) intervene. These variables are distances δ and ∆, and cardinalities.

Proof: We first introduce the following distance

∆ 1 = min i∈A,j∈B1 |x i -x j |,
and ∆ 2 for group B 2 . We will prove the result dealing with ∆ 1 , knowing the expression is similar for ∆ 2 . We begin by writing the detailed expression of

∆ + 1 ∆ + 1 = ∆ 1 + D A + D B where D A = j∈A h(x j -x m ) σ m - h(x j -x m-1 ) σ m-1 ,
and where

D B = D B1 + D B2 with D B1 = j∈B1 h(x m-1 -x j ) σ m-1 - h(x m -x j ) σ m , D B2 = j∈B2 h(x j -x m ) σ m - h(x j -x m-1 ) σ m-1 ,
with σ m and σ M defined in [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: convergence conjectures and theorems[END_REF].

It can be verified that quantities σ A (δ, ∆) and σ B (δ, ∆) defined in ( 25) and ( 26) satisfy

σ A (δ, ∆) ≤ σ m , and σ B (δ, ∆) ≤ σ m-1 .
Then by using Lemma 2, we can get the following lower bounds for ∆ ≥ x

D A ≥ - |A| h(∆) σ B (δ, ∆) = -|A| η A (δ, ∆), D B1 ≥ - |B 1 | h(∆) σ A (δ, ∆) , D B2 ≥ - |B 2 | h(2∆ + δ) σ B (δ, ∆) , giving D B ≥ -(n -|A|) max h(∆) σ A (δ, ∆) , h(2∆ + δ) σ B (δ, ∆) ,
leading to [START_REF] Durkheim | The division of labor in society[END_REF] and finally to [START_REF] Bhattacharyya | On the convergence of the Hegselmann-Krause system[END_REF].

F. Upper bound on evolution of ratio δ/∆

In this subsection, we combine the two previous bounds on δ and ∆ to bound the ratio of these two quantities.

Proposition 3: Under Assumptions 1-3 and for ∆ ≥ x we have

δ + ∆ + ≤ ρ(δ, ∆) δ ∆ , where ρ(δ, ∆) = µ(δ, ∆) η(δ, ∆) , (27)
with µ(δ, ∆) and η(δ, ∆) defined in Propositions 1 and 2, from which this proposition falls directly.

IV. MAIN RESULT Before introducing the main result allowing to detect cluster formation, we introduce a few notation in order to ease its formulation.

We first define the domain D of the plane (δ, ∆) where the quantity ρ(δ, ∆) defined in ( 27) is smaller than one, meaning that the group diameter is decreasing faster than the distance in between the group and the other agents. More formally, this domain is defined as

D = {(δ, ∆) ∈ (R + ) 2 / ρ(δ, ∆) ≤ 1}, (28)
where ρ(δ, ∆) is defined in [START_REF] Newman | Modularity and community structure in networks[END_REF]. We also introduce the worst case function W as

W = (R + ) 2 → R 2 (δ, ∆) → µ(δ, ∆) δ , η(δ, ∆) ∆ . (29)
As domain D does not provide enough information about the duration of a transient cluster formation, we now define a sequence of decreasing (in the sense of inclusion) domains D τ , for τ ∈ N, such that

D 0 = {(δ w , ∆ w ) ∈ D / [0, δ w ] × [∆ w , +∞) ⊆ D}, (30)
and

D τ +1 = {(δ w , ∆ w ) ∈ D τ ∩W -1 (D τ ) / [0, δ w ] × [∆ w , +∞) ⊆ D τ ∩W -1 (D τ )}, (31)
with W -1 (D τ ) the pre-image of D τ by function W . The definition of D τ allows to show that when W (δ, ∆) ∈ D τ then (δ + , ∆ + ) ∈ D τ , by Proposition 1 and 2. Iterating this property over time is the argument to prove the following theorem.

Theorem 1: Given a selected group of agents A satisfying (δ(0), ∆(0)) ∈ D T , for T ∈ N and D T defined in (31), we have ∀t ∈ {0, . . . , T },

δ(t + 1) ∆(t + 1) ≤ δ(t) ∆(t) , (32)
in other words, the group A is forming a cluster according to Definition 3 during at least T iterations.

This theorem allows us to detect when a transient cluster is forming and gives how many iterations it will be guaranteed to remain a cluster before other agents merge with it.

Remark 3: Notice that it makes sense to detect cluster formation for a limited period of time only since the decrease of ratio δ/∆ may occur for some time but may stop when ∆ becomes small enough. This is illustrated in Section V.

Proof:

Using Proposition 3 and equation ( 28), statement (32) is implied by ∀t ∈ {0, . . . , T }, (δ(t), ∆(t)) ∈ D.

Instead of proving the previous equation, we show that for t ∈ {0, . . . , T }, (δ(t), ∆(t)) ∈ D T -t . We will prove it by induction. The basis case is obvious because if (δ(0), ∆(0)) ∈ D T ⊆ D then the inequality (32) stands for t = 0.

We prove the inductive step by taking t ∈ {0, . . . , T -1} and assuming (δ(t), ∆(t)) ∈ D T -t , then we have

W (δ(t), ∆(t)) = δ w (t + 1), ∆ w (t + 1) ∈ D T -t-

V. NUMERICAL ILLUSTRATION

In the present section, we run numerical simulations to validate the main result of the paper. We choose an influence function f = x → exp(-|x|). The number of agents in the system is n = 30, and the number of agents in group A is set at |A| = 10. These parameters enable us to compute the domains D and D 0 to D 10 as defined in equations ( 28), ( 30) and (31), and drawn in Figure 3. We then consider a system with two different sets of initial conditions. In the first case a group of agent denoted A lies in set D while the other denoted A do not. Groups A and A are composed of 10 agents and are initially distributed over two extreme points respectively spaced by δ(0) = 1.5 and δ(0) = 2.5. All the remaining agents are located at one isolated point at distance ∆(0) = 5.7 and ∆(0) = 3.5 from A and A , respectively. These initial conditions are displayed for each group in Figure 3 with squares.

The time evolution of the two sets of initial conditions are depicted in Figure 4 and Figure 5 and also in Figure 3. In the configuration displayed in Figure 4, group A quickly contracts to form a transient cluster that will merge with the rest of the agents after a greater amount of iterations. By contrast, in Figure 5, the group A do not form a transient cluster despite having δ(0) < ∆(0). This illustrates the following fact. The property of having a diameter smaller than the distance to the rest of the agents is not a sufficient criterion to detect a cluster.

Initially, the couple (δ(0), ∆(0)) for group A is in D 3 . Using Theorem 1, we know that group A is forming a cluster during at least 3 iterations, which is conservative compared the actual evolution of group A. This is because Theorem 1 deals with the worst case scenario. At time t = 3, the distance ∆(t) has increased and (δ(t), ∆(t)) for group A is now in domain D 9 . As a consequence, knowing the position of the couple (δ(t), ∆(t)) at time t = 3, Theorem 1 ensures that the transient cluster is forming for at least 9 more iterations.

Unlike group A, the dynamics of group A do not present a transient cluster formation. This is consistent with the (δ(t), ∆(t)) trajectory of group A as displayed in Figure 3 : group A starts out of domain D and never reaches it. This underlines this importance of a formal criterion to detect whether a cluster is forming or not.

VI. CONCLUSIONS

In this paper we analysed transient cluster formation for general models of opinion dynamics including the generalized Hegselmann-Krause model. We proposed a criterion for detecting cluster formation. Our criterion can be computed locally using the initial conditions and only requires the knowledge of the diameter of the cluster candidate and its distance to the rest of the agents. Moreover, our main result provides a lower bound on how many iterations a cluster forms before merging with other agents. Finally our results are illustrated by a numerical example.

Fig. 1 .

 1 Fig.1. Simulation over 700 iterations of a multi-agent system with an influence function f (x) = exp(-|x|) and with 20 agents, with 10 agents initially above 9 and the others below 2, which eventually merge and converge to one final opinion, the consensus. One can see 3 phases : rapid convergence of agents into clusters, clusters remaining far apart and rapid convergence of clusters toward consensus. Note the log x-scale.

Fig. 2 .

 2 Fig. 2.Generalized Hegselmann-Krause influence functions satisfying Assumption 1-3 and Hegselmann-Krause influence function.

Fig. 3 .

 3 Fig.3. Trajectories, in the (δ, ∆) plane, of groups A and A composed of 10 agents in a system of 30 agents, whose time evolutions are shown in Figures4 and 5.The domain D in dashed red and domains Dτ in black are also drawn. The black dashed-dotted line stands for the first bisector.

Time

Fig. 4 .

 4 Fig.[START_REF] Reynolds | Flocks, herds and schools: A distributed behavioral model[END_REF]. Evolution of a system of N = 30 agents using influence function f (x) = exp(-|x|). Group A in blue is composed of 10 agents where 5 are initially at x = 0 and 5 at x = 1.5. The rest of the agents are in green and initially at x = 7.2. Note the log x-scale.

Fig. 5 .

 5 Fig.5. Evolution of a system of N = 30 agents using influence function f (x) = exp(-|x|). Group A in red is composed of 10 agents where 5 are initially at x = 0 and 5 at x = 2.5. The rest of the agents are in green and initially at x = 6. Note the log x-scale.


Andrzej Palczewski
email: a.palczewski@mimuw.edu.pl

Risk minimizing strategies for tracking a stochastic target *

Keywords: optimal portfolios, stochastic target, benchmark tracking

We consider a stochastic control problem of beating a stochastic benchmark. The problem is considered in an incomplete market setting with external economic factors. The investor preferences are modelled in terms of HARA-type utility functions and trading takes place in a finite time horizon. The objective of the investor is to minimize his expected loss from the outperformance of the benchmark compared to the portfolio terminal wealth, and to specify the optimal investment strategy. We prove that for considered loss functions the corresponding Bellman equation possesses a unique solution. This solution guaranties the existence of a well defined investment strategy. We prove also under which conditions the verification theorem for the obtained solution of the Bellman equation holds.

Introduction

We analyze the optimal portfolio and investment policy for an investor who is concerned about his wealth relative to the performance of a given benchmark. The benchmark evolves stochastically over time and the investor's objective is to minimize his loss with respect to this benchmark by investing in a portfolio of stochastically evolving financial instruments. Since the benchmark is not necessarily perfectly correlated with the investment opportunities, we are in the framework of an incomplete market, and there is no investment policy under which the investor can outperform the benchmark with certainty.

The portfolio problem where the objective is to exceed the performance of a selected target benchmark is sometimes called an active portfolio management. It is well known that many professional investors apply this benchmarking procedure. However, many small investors follow a benchmarking procedure as well, by trying to beat inflation, exchange rates, or other market indices.

The problem of an investment portfolio which outperforms a given benchmark has been studied for a long time. For objectives such as maximizing the probability that the investor's wealth achieves a certain performance goal relative to the benchmark, before falling below to a predetermined shortfall, or minimizing the expected time to reach the performance goal, the problem is studied by Browne [START_REF] Browne | Reaching Goals by a Deadline: Digital options and Continuous-Time Active Portfolio Management[END_REF], [START_REF] Browne | Beating a moving target: optimal portfolio strategies for outperforming a stochastic benchmark[END_REF]. For the special case where the benchmark is perfectly correlated with the investment opportunities, these problems over a finite-horizon are analyzed in [START_REF] Browne | Reaching Goals by a Deadline: Digital options and Continuous-Time Active Portfolio Management[END_REF], and for a more general model that the benchmark is not perfectly correlated with the investment opportunity in [START_REF] Browne | Beating a moving target: optimal portfolio strategies for outperforming a stochastic benchmark[END_REF]. The problem of finding the minimal initial data of a controlled process which guarantees to reach the benchmark with a given probability of success or, more generally, with a given level of expected loss was first introduced by Föllmer and Leukert [START_REF] Föllmer | -Quantile hedging[END_REF] in the context of quantile hedging. This approach has been then extended to the stochastic target problem studied by Soner and Touzi [START_REF] Soner | Stochastic target problems, dynamic programming and viscosity solutions[END_REF][START_REF] Soner | Dynamic programming for stochastic target problems and geometric flows[END_REF], and in a number of papers by Bouchard et. al. [START_REF] Bouchard | Stochastic target problems with controlled loss[END_REF][START_REF] Bouchard | Optimal control under stochastic target constraints[END_REF][START_REF] Bouchard | A stochastic target approach for P&L matching problems[END_REF].

In opposition to the majority of previously mentioned papers, in this paper, we study a loss minimization objective when the prices of financial instruments are functions of external economic factors. A similar problem but without taking into account economic factors is solved by Browne [START_REF] Browne | Beating a moving target: optimal portfolio strategies for outperforming a stochastic benchmark[END_REF]. The absence of economic factors makes the problem much simpler as the HJB equation is reduced in that case to an ODE. In the presence of external factors the HJB equation becomes a multidimensional nonlinear PDE for which the existence of solutions is a challenging problem. We solve this problem using the well developed theory of quasilinear parabolic equations. We also show that under suitable regularity assumptions the verification theorem holds. Hence, the obtained solution to the HJB equation is a solution to the optimization problem. The plan of the paper is as follows. In Section 2, we present the portfolio problem arising from the active portfolio management. In Section 3, we show that, under additional assumptions on the loss function, we can find a smooth solution to the HJB equation and construct effectively an optimal investment strategy. Section 4 is devoted to the formulation and proof of the verification theorem.

2 The portfolio problem.

We consider the portfolio problem in which the prices of securities are functions of external state variables (economic factors). Our goal is to construct a portfolio which can outperform a stochastic benchmark. We consider a general setting of the problem. In particular, the risk factors which define the dynamics of the benchmark can be different from the risk factors in the dynamics of securities. Hence, the problem is an incomplete market problem.

The setting of the market is as follows: we have a market defined on a probability space (Ω, F, P) with the filtration (F t ) t∈[0,T ] generated by d-dimensional standard Wiener process W (t) = (W 1 , . . . , W d ) (in what follows we treat W as a column vector). On that probability space we have N stochastic processes describ-ing the prices of securities with the dynamics

dS i (t) S i (t) = µ i (t, R)dt + d j=1 σ ij (t, R)dW j (t), i = 1, 2, . . . , N, (1)
where µ i and σ ij depend on an M -dimensional vector of economic factors R.

We assume that the dynamics of factors R follow the Markovian diffusion process

dR m (t) = µ r m (t, R)dt + d i=1 b mi (t, R)dW i (t), m = 1, 2, . . . , M. (2)
It is convenient to switch to vector notation and introduce the matrices σ = (σ ij ), B = (b ij ) and the column vectors µ = (µ 1 , . . . , µ N ) ′ , µ r = (µ r 1 , . . . , µ r M ) ′ . (Here and in what follows x ′ denotes the transpose of the matrix or vector x.) Assumption 2.1. For the model of security prices we assume that coefficients µ(t, r) and σ(t, r) are deterministic functions bounded and continuous for t ∈ [0, T ] and r ∈ R M . For the model of economic factors we make typical assumptions which guarantee the existence of strong solutions to equation (2), i.e., we assume that µ r (t, r) and B(t, r) are deterministic continuous functions of their arguments, which in addition fulfil the estimates

µ r (t, r 1 ) -µ r (t, r 2 ) + B(t, r 1 ) -B(t, r 2 ) ≤ c r 1 -r 2 , (3)
µ r (t, r) 2 + B(t, r) 2 ≤ c 2 (1 + r 2 ), (4)
for t ∈ [0, T ], r, r 1 , r 2 ∈ R M , where c is a positive constant.
We analyze the stochastic target problem in an incomplete market assuming that the dimension of risk factors is high, i.e. dimmension d of the Wiener process W is high, and the number of securities and economic factors much lower. This means that d ≫ N and d ≫ M . To guarantee well-posedness and solvability of the optimization problem, we have to make additional assumptions. Assumption 2.2. About the model of securities dynamics we assume a "partial invertibility" of the model, i.e., the matrix Σ = σσ ′ is nonsingular and Σ -1 (t, r) is bounded for t ∈ [0, T ] and r ∈ R M . This in fact means that securities are driven by N risk factors and limited to these N dimensions the security market is complete. About the model of dynamics of economic factors we assume that the matrix BB ′ is positive definite. Strictly speaking, we postulate that there exist positive constants

ν 1 , ν 2 such that for any x ∈ R M 0 < ν 1 x 2 ≤ x ′ BB ′ x ≤ ν 2 x 2 .
The stochastic benchmark is modelled as a general log-normal stochastic process H(t) which fulfils the equation

dH(t) := H(t) µ H dt + ξdW (t) ,
where ξ = (ξ 1 , . . . , ξ d ) ′ is a column vector, and coefficients µ H and ξ are deterministic functions of t ∈ [0, T ] and r ∈ R M .

We consider now a portfolio V (t) consisting of assets S i (t), i = 1, . . . , N . Denoting by π i the fraction of the total wealth V (t) invested in the security S i , we can write

dV π (t) = V π (t) µ V dt + θ V dW (t) , (5)
where after introducing the column vector π = (π 1 , . . . , π N ) ′ we have µ V = µ ′ π and θ V = σ ′ π.

We define now the new process

X(t) = H(t) V π (t) . (6)
This approach enables us to consider in the same framework losses and gains. Such an approach is not quite new in the financial literature. It is applied by Browne [START_REF] Browne | Reaching Goals by a Deadline: Digital options and Continuous-Time Active Portfolio Management[END_REF][START_REF] Browne | Beating a moving target: optimal portfolio strategies for outperforming a stochastic benchmark[END_REF]. A similar quotient is used by Dai Pra, Runggaldier and Tolotti [START_REF] Pra | Pathwise optimality for benchmark tracking[END_REF] who optimize the quadratic loss in the benchmark tracking problem.

For the process X(t) we obtain the following equation of evolution

dX X = µ X dt + θdW, (7)
where

θ = ξ -θ V , µ X = µ H -µ V -θ ′ θ V .
We optimize the process X(t) with respect to the vector of strategies π. Admissible strategies for our problem are defined as follows.

Definition 2.1. Let U be a complete, separable metric space and 0 < T < ∞. We define the set of admissible strategies Π(t, x, r) as fulfilling the conditions:

1. π : [t, T ] × Ω → U ⊆ R N is measurable, bounded and {F τ } τ ≥t -adapted, for each π ∈ Π(t, x, r), 2. X(t) = x (budget constraint), 3. R(t) = r.
We consider the optimization problem in the framework of utility theory. This means that we fix a utility function g and optimize the terminal value of process X measured by g. The optimization problem is of the form

min π∈Π(t,x,r) E g(X(T ))|X(t) = x, R(t) = r . (8)
In fact, it is better to call g a loss function as our goal is to minimize losses and not to maximize gains.

Under Assumptions 2.1, 2.2 and Definition 2.1, equation ( 7) admits the unique solution and the value function

u(t, x, r) := min π∈Π(t,x,r) E g(X(T ))|X(t) = x, R(t) = r (9)
is well defined.

With this value function we arrive at the following Hamilton-Jacobi-Bellman equation

∂ t u + inf π∈Π µ X x∂ x u + (µ r ) ′ ∇ r u + 1 2 θ ′ θx 2 ∂ xx u+ + 1 2 BB ′ (∇ r ⊗ ∇ r u) + xθ ′ B ′ ∇ r ∂ x u = 0, u(T, x, r) = g(x), (10)
where ∇ r denotes the gradient operator with respect to vector variable r (a column vector), ∂ t and ∂ x denote differential operator with respect to scalar variables t and x, respectively. ∇ r ⊗ ∇ r has the following meaning: when x and y are ndimensional column vectors then x⊗y denotes the n×n matrix xy ′ , i.e. ∇ r ⊗∇ r u is a matrix of all second order derivatives of u with respect to variables r i and r j , i, j = 1, . . . , M .

Smooth solutions of the HJB equation

To obtain smooth solutions to the HJB equation [START_REF] Soner | Dynamic programming for stochastic target problems and geometric flows[END_REF], we make the following assumption.

Assumption 3.1. The loss function g(x) is from the generalized HARA class and is given by the expression g(x) = cx α , for α > 1 and x ∈ [0, ∞).

Remark 3.1. In fact, from the technical point of view, we can assume only that g is such that α = -1. The assumption α > 1 is essential when we want to interpret g(x) as a loss function.

Under the above assumption, we postulate that the value function can be factorized in the form u(t, x, r) = g(x)q(t, r).

Substituting the above factorization into equation [START_REF] Soner | Dynamic programming for stochastic target problems and geometric flows[END_REF], we obtain the following PDE problem for q:

∂ t q + inf π∈Π αµ X q + (µ r ) ′ ∇ r q + α(α -1) 2 θ ′ θq+ + 1 2 BB ′ (∇ r ⊗ ∇ r q) + αθ ′ B ′ ∇ r q = 0, q(T, r) = 1. (12)
From equation ( 12) we can derive formally the optimal investment strategy

π * = π 0 + π 1 ∇ r q q , (13)
where

π 0 = 1 1 + α Σ -1 (µ + ασξ), π 1 = 1 1 + α Σ -1 σB ′ .
Substituting expression (13) into the HJB equation ( 12) we obtain

∂ t q + αµ * q + (µ r ) ′ ∇ r q + α(α -1) 2 (θ * ) ′ θ * q+ + 1 2 BB ′ (∇ r ⊗ ∇ r q) + α(θ * ) ′ B ′ ∇ r q = 0. (14)
In this equation, µ * denotes the value of µ X , and θ * the value of θ evaluated at the point of the optimal strategy π * . After rearrangements, the HJB equation ( 14) takes the form

∂ t q + A 0 (∇ r ⊗ ∇ r q) + A 1 ∇ r q ⊗ ∇ r q q + A 2 ∇ r q + A 3 q = 0, (15)
where

A 0 = 1 2 BB ′ , A 1 = 1 2 α(α + 1)(π 1 ) ′ Σπ 1 -αBσ ′ π 1 , A 2 =µ r + α(α + 1)(π 1 ) ′ Σπ 0 -α 2 (π 1 ) ′ σξ + α Bξ -(π 1 ) ′ µ -Bσ ′ π 0 , A 3 = 1 2 α 2 (π 0 ) ′ Σπ 0 -2(π 0 ) ′ σξ + ξ ′ ξ + α µ H -µ ′ π 0 - 1 2 ξ ′ ξ .
Equation (15) has to be solved in the strip 0 ≤ t ≤ T with the terminal condition q(T, r) = 1.

Equation ( 15) is a quasilinear parabolic equation which possesses a solution provided this solution is bounded away from zero. To find this solution we make the substitution z = ln q.

For the new function z we obtain the equation

∂ t z + A 0 (∇ r ⊗ ∇ r z) + (A 0 + A 1 )∇ r z ⊗ ∇ r z + A 2 ∇ r z + A 3 = 0, (17)
with the terminal condition z(T, r) = 0.

To solve equation (17) with condition (18), we use well known results in the theory of quasilinear parabolic equations.

Let us consider a boundary value problem for a n-dimensional quasilinear parabolic equation

∂ t w - n i,j=1 a ij (t, x)∂ x i x j w + a(t, x, w, ∂ x w) = 0, for (t, x) ∈ O T , w(t, x) = ψ(t, x), for (t, x) ∈ Γ T , (19) in a bounded domain O T = [0, T ] × O,
Γ T = ∂O × [0, T ] ∪ O × {t = 0}.
Theorem 7.4 in Chapter 6 of the book by Ladyzhenskaya, Solonnikov and Uraltseva [START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF] guarantees that, under suitable on coefficients a ij , a and boundary function ψ, there exists a unique solution to the boundary value problem ( 19) in H 1+β/2,2+β (O T ).

To apply the above result to equation ( 17), we have to make additional assumptions.

Assumption 3.2. Let functions µ, µ r , µ H , σ, B, ξ and Σ -1 be Hölder continuous functions of t with the Hölder exponent β/2, and Hölder continuous functions of r with the Hölder exponent β , for some β > 0. Now we can prove our main theorem.

z(t, r) = 0 for (t, r) ∈ ∂O × [0, T ]. (20)
The above defined augmented problem fulfils already the assumptions of Theorem 7.4 in Chapter 6 of [START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF]. Due to this theorem, there exists a unique solution of equation ( 17) with terminal condition (18) and boundary condition (20). This solution together with its derivatives can be estimated in O T with constants which depend only on constants present in the Assumptions, and not on the size of domain O. Hence, the solution which exists in any bounded domain O T belongs to H 1+β/2,2+β (O T ) and is uniformly bounded together with its derivatives independently of the size of the domain. Then we can consider a increased sequence of bounded smooth domains O n that fill in the whole R M and solutions z n to problem (17), ( 18), (20) with O replaced by O n . By the standard Arzela-Ascoli theorem, we can choose a subsequence of z n which converges to a function which is a solution to (17)-( 18) on [0, T ] × R M .

Remark 3.2. In many situations, economic factors should be restricted to nonnegative values only. In that cases, Theorem 3.1 is still applicable as we can construct a sequence of bounded smooth domains O n approximating the space R M + .

Corollary 3.1. Let us observe that due to Theorem 3.1 fuction z(t, r) and its derivatives are bounded. Returning back to the original function q, we conclude that q(t, r) is bounded away from zero and the quotient q rm /q is bounded. It follows than that the optimal investment strategy given by expression ( 13) is bounded and admissible in accordance with Definition 2.1.

Verification theorem

Theorem 3.1 guaranties a smooth solution to the terminal problem ( 17)-( 18). This solution is not necessarily a solution to the optimization problem. To prove the optimality, we need some additional results. First, we have to show that the solution to problem (17)-( 18) is a function of class C 1,2 on [0, T ] × R M .

To this end, we can use Theorem 8.1 from Chapter 6 of [START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF] which guarantees the existence of a unique solution in H 1+β/2,2+β (Q T ), where Q T = [0, T ] × R n , to the Cauchy problem for a quasilinear parabolic equation

∂ t w - n i,j=1 a ij (t, x)∂ x i x j w + a(t, x, w, ∂ x w) = 0, in (0, T ] × R n , w(0, x) = ψ(x), in R n , (21)
if in addition to the assumptions of Theorem 7.4 from [START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF] the coefficients of the equation can be uniformly estimated on every bounded set with the bound independent of the size of this set.

To apply the above mentioned theorem, we make the following assumption Assumption 4.1. Let functions µ, µ r , µ H , σ, B, ξ and Σ -1 , in addition to being Hölder continuous, be uniformly bounded for all (t, r), t ∈ [0, T ], r ∈ R M . Proof: The proof follows straightforwardly from Theorem 8.1 in Chapter 6 of [START_REF] Ladyzhenskaya | Linear and Quasilinear Equations of Parabolic Type[END_REF]. That theorem states that the solution to problem (21) is unique and belongs to

H 1+β/2,2+β ([0, T ] × R M ) for some β > 0. It is obvious that such a solution is a function of class C 1,2 .
To use the classical stochastic verification theorem (cf. Theorem 5.1 in Chapter 5 of the book by Yong and Zhou [START_REF] Yong | Stochastic Controls[END_REF]), we have to prove the following simple lemma.

Lemma 4.1. Let z(t, r) be the unique solution to the boundary value problem (17)-( 18), which exists due to Theorem 4.1 in [0, T ] × R M . Let q(t, r) = exp z(t, r) and π * be given by equation (13). Then

inf π∈Π αµ X q + (µ r ) ′ ∇ r q + α(α -1) 2 θ ′ θq + 1 2 BB ′ (∇ r ⊗ ∇ r q) + αθ ′ B ′ ∇ r q = = αµ * q + (µ r ) ′ ∇ r q + α(α -1) 2 (θ * ) ′ θ * q + 1 2 BB ′ (∇ r ⊗ ∇ r q) + α(θ * ) ′ B ′ ∇ r q,
where µ * denotes the value of µ X , and θ * the value of θ evaluated at the point of the optimal strategy π * .

Proof: The proof is straightforward as the left and right hand sides of the equation in the Lemma are left hand sides of equations ( 12) and ( 14), respectively. But equation ( 14) has been obtained from equation (12) upon substitution (13). The fact that π * is an admissible investment strategy has been already stated in Corollary 3.1.

From Theorem 4.1 and Lemma 4.1, we easily obtain the verification theorem.

Theorem 4.2. Under assumptions of Theorem 4.1, the function u(t, x, r) = g(x)q(t, r)

is the value function to the optimization problem [START_REF] Soner | Stochastic target problems, dynamic programming and viscosity solutions[END_REF], where g(x) fulfils the conditions of Assumption 3.1 and q(t, r) = exp z(t, r) with z(t, r) being the solution of the boundary value problem (17)-( 18), which exists due to Theorem 4.1.

Conclusions

In this paper, we have solved the stochastic optimization problem for the loss minimization with the state variable being the ratio of a stochastic benchmark to an investment portfolio. The control parameter of this problem is the portfolio investment strategy. The problem is solved in a market model of N securities being log-normal stochastic processes and depending on M external economic factors. The stochastic benchmark is also a log-normal process but the set of risk factors on which this benchmark depends can be larger than the set of risk factors of the securities making the whole problem an incomplete market problem. The stochastic optimization problem has been reduced to the HJB equation which, in this case, is a multidimensional quasilinear parabolic equation. Using the general theory of such equations, we have proved that under suitable regularity conditions the HJB equation possesses a unique solution which is sufficiently smooth to guarantee the fulfilment of the stochastic verification theorem. Hence, the solution to the HJB equation is a unique solution to the initial optimization problem. This is a natural extention of similar results obtained in a less general setting without the dependence of security prices on external economic factors.

Theorem 3 . 1 .

 31 Under Assumptions 2.1, 2.2, 3.1, 3.2 and the assumptions of Definition 2.1, there exists a unique solution z(t, r) to the terminal problem (17)-(18) and |z|, |∂ t z|, |∂ x i z|, |∂ x i x j z| are bounded in [0, T ] × R M . This solution belongs to H 1+β/2,2+β (O T ), where O T = [0, T ] × O and O is a bounded domain in R M . Proof: Let us consider the terminal problem (17)-(18) in O T = [0, T ] × O, where O is a fixed bounded domain in R M . To solve this problem, we supplement equation (17) and terminal condition (18) with the boundary condition

Theorem 4 . 1 .

 41 Under the assumptions of Theorem 3.1 and Assumption 4.1, there exists a unique solution to the terminal problem (17)-(18). This solution belongs to H 1+β,2+β/2 ([0, T ] × R M ) and the estimates of Theorem 3.1 hold for (t, r) ∈ [0, T ] × R M . In particular, the solution is a C 1,2 function on [0, T ] × R M .

  where O is a bounded domain in R n with the boundary of class H 2+β , and

grant no. NN-201-547838 is gratefully acknowledged.

A natural question which arises is the extension of the obtained results to a market with less restrictive assumpions. The most severe of these assumptions is the boundedness of the coefficients in the whole domain and the lifting of these restrictions will be the subject of future research.