Andaluzia Matei
email: andaluziamatei@inf.ucv.ro

On the weak solvability and the optimal control of a frictional contact problem with normal compliance

Keywords: contact model, friction, normal compliance, weak solutions, optimal control, optimality condition, convergences

In the present work we consider a frictional contact model with normal compliance. Firstly, we discuss the weak solvability of the model by means of two variational approaches. In a first approach the weak solution is a solution of a quasivariational inequality. In a second approach the weak solution is a solution of a mixed variational problem with solution-dependent set of Lagrange multipliers. Nextly, the paper focuses on the boundary optimal control of the model. Existence results, an optimality condition and some convergence results are presented.

Introduction

The present paper focuses on the weak solvability and the boundary optimal control of the following contact model. Problem 1. Find a displacement field u : Ω → R 3 and a stress field σ : Ω → S 3 such that Div σ + f 0 = 0 in Ω,

σ = Fε(u) in Ω, u = 0 on Γ 1 , σν = f 2 on Γ 2 , -σ ν = p ν (u ν -g a ) on Γ 3 , σ τ ≤ p τ (u ν -g a ), σ τ = -p τ (u ν -g a ) uτ uτ if u τ = 0 on Γ 3 .
Notice that Ω is a bounded domain of R 3 with smooth enough boundary, partitioned in three measurable part Γ 1 , Γ 2 and Γ 3 ; u ν = u • ν, u τ = u -u ν ν, σ ν = (σν) • ν, σ τ = σν -σ ν ν, " • " denotes the inner product of two vectors,

• denotes the Euclidean norm, ν is the unit outward normal vector. Problem 1 is a contact problem with the normal compliance condition, associated to the Coulomb's law of dry friction. A normal compliance condition was firstly proposed in [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF]. Then, such a contact condition was used in many models, see e.g. the papers [START_REF] Andersson | A quasistatic frictional problem with normal compliance[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF][START_REF] Klarbring | A global existence result for the quasistatic frictional contact problem with normal compliance[END_REF][START_REF] Rochdi | Quasistatic viscoelastic contact with normal compliance and friction[END_REF].

In the normal compliance contact condition

-σ ν = p ν (u ν -g a ) on Γ 3 ,
p ν is a nonnegative prescribed function which vanishes for negative argument and g a > 0 denotes the gap (the distance between the body and the obstacle on the normal direction). When u ν < g a there is no contact and the normal pressure vanishes. When there is contact then u ν -g a is positive and represents a measure of the interpenetration of the asperities. Then, the normal compliance condition shows that the foundation exerts a pressure on the body which depends on the penetration. For details on the physical significance of the model we refer to [START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF]. The rest of the paper has the following structure. Section 2 is devoted to the weak solvability of the model by means of two variational approaches. In Section 3 we discuss an optimal control problem which consists of leading the stress tensor as close as possible to a given target, by acting with a control on a part of the boundary.

There are several works concerning the optimal control of variational inequalities, see for instance [START_REF] Barbu | Optimal Control of Variational Inequalities[END_REF][START_REF] Bonnans | Pontryagin's principle in the control of semiliniar elliptic variational inequalities[END_REF][START_REF] Friedman | Optimal Control for Variational Inequalities[END_REF][START_REF] Lions | Contrôle optimale des systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Mignot | Contrôle dans les inéquations variationnelles elliptiques[END_REF][START_REF] Mignot | Optimal control in some variational inequalities[END_REF][START_REF] Neittaanmaki | Optimization of Elliptic Systems: Theory and Applications[END_REF][START_REF] Sokolowski | Introduction to Shape Optimization[END_REF]. Nevertheless, only few works are devoted to the optimal control of contact problems, see [START_REF] Amassad | Optimal control of an elastic contact problem involving Tresca friction law[END_REF][START_REF] Capatina | Boundary optimal control for quasistatic bilateral frictional contact problems[END_REF][START_REF] Matei | Boundary optimal control for nonlinear antiplane problems[END_REF]. The present paper adds a new contribution.

On the weak solvability of the model

In this section we shall indicate two variational approaches in the study of Problem 1. Let us make the following assumptions. r) is measurable on Γ 3 , for any r ∈ R and p e (x, r) = 0 for all r ≤ 0, a.e. x ∈ Γ 3 .

Assumption 1 F : Ω×S 3 → S 3 , F(x, ε) = (F ijkl (x)ε jk ) for all ε = (ε ij ) ∈ S 3 , a.e. x ∈ Ω. F ijkl = F jikl = F klij ∈ L ∞ (Ω), 1 ≤ i, j, k, l ≤ 3. There exists m F > 0 such that F(x, τ ) : τ ≥ m F τ 2 for all τ ∈ S 3 , a.e. x in Ω. Assumption 2 f 0 ∈ L 2 (Ω) 3 , f 2 ∈ L 2 (Γ 2 ) 3 . Assumption 3 p e : Γ 3 × R → R + (e ∈ {ν, τ }). There exists L e > 0 such that |p e (x, r 1 ) -p e (x, r 2 )| ≤ L e |r 1 -r 2 | for all r 1 , r 2 ∈ R a.e. x ∈ Γ 3 . The mapping x → p e (x,
Assumption 4 g a ∈ L 2 (Γ 3 ), g a (x) ≥ 0, a.e. x ∈ Γ 3 . Assumption 5 m F > c 2 0 (L ν + L τ ).
Assumption 5 is a smallness assumption which was introduced mainly for mathematical reasons. However, for some materials and frictional contact conditions we have appropriate constants m F , L ν and L τ which fulfill Assumption 5. Notice that " : " denotes the inner product of two tensors and c 0 = c 0 (Ω, Γ 1 , Γ 3 )> 0 is a "trace constant" such that:

v L 2 (Γ3) 3 ≤ c 0 v V for all v ∈ V, (1)
where

V = { v ∈ H 1 (Ω) 3 | v = 0 a.e. on Γ 1 }.
In a first approach the weak solution is a solution of a quasivariational inequality having as unknown the displacement field.

Problem 2. Find a displacement field u ∈ V such that (Au, v -u) V + j(u, v) -j(u, u) ≥ (f , v -u) V for all v ∈ V. (2)
Herein,

A : V → V (Au, v) V = (Fε(u), ε(v)) L 2 (Ω) 3×3 s , j : V × V → R j(u, v) = Γ3 p ν (u ν -g a )|v ν | dΓ + Γ3 p τ (u ν -g a ) v τ dΓ, (f , v) V = Ω f 0 • v dx + Γ2 f 2 • v dΓ. (3)
We have the following existence and uniqueness result.

Theorem 1. Under Assumptions 1-5, Problem 2 has a unique weak solution.

For the proof we refer to Theorem 5.30 in [START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF].

The second approach is a mixed variational approach. The mixed variational formulations are related to modern numerical techniques in order to approximate the weak solutions of contact models. Referring to numerical techniques for approximating weak solutions of contact problems via saddle point technique, we send the reader to, e.g., [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF][START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF][START_REF] Wohlmuth | Discretization Methods and Iterative Solvers Based on Domain Decomposition[END_REF]. The functional frame is the following one.

V = {v ∈ H 1 (Ω) 3 | v = 0 a.e. on Γ 1 }; S = {w| Γ3 | w ∈ V }; D = S .
Notice that w| Γ3 denotes the restriction of the trace of the element w ∈ V to Γ 3 . Thus, S ⊂ H 1/2 (Γ 3 ; R 3 ) where H 1/2 (Γ 3 ; R 3 ) is the space of the restrictions on Γ 3 of traces on Γ of functions of H 1 (Ω) 3 . We use the Sobolev-Slobodeckii norm

ζ S = Γ3 Γ3 ζ(x) -ζ(y) 2
xy 3 ds x ds y

1/2 . For each ζ ∈ S, ζ ν = ζ • ν and ζ τ = ζ -ζ ν ν a.e. on Γ 3 . Let us consider f ∈ V , see (3)
a(•, •) : V × V → R, a(u, v) = Ω Fε(u(x)) : ε(v(x)) dx; b(•, •) : V × D → R b(v, µ) = µ, v| Γ3 .
Also, we define a variable set Λ = Λ(ϕ),

Λ(ϕ) = {µ ∈ D | µ, v| Γ3 ≤ Γ3 (p ν (x, ϕ ν (x) -g a )|v ν (x)| + p τ (x, ϕ ν (x) -g a ) v τ (x) )dΓ v ∈ V }.
Notice that •, • denotes the duality pairing between D and S. The second variational formulation of Problem 1 is the following one.

Problem 3. Find u ∈ V and λ ∈ Λ(u) ⊂ D such that a(u, v) + b(v, λ) = (f , v) V for all v ∈ V, b(u, µ -λ) ≤ 0 for all µ ∈ Λ(u).
In this second approach, a weak solution is a pair consisting of the displacement field and a Lagrange multiplier related to the friction force.

Theorem 2. Under Assumptions 1-4, Problem 3 has at least one solution.

The proof of Theorem 2, based on the abstract results we have got in [START_REF] Matei | On the solvability of mixed variational problems with solution-dependent sets of Lagrange multipliers[END_REF], can be found in the very recent paper [START_REF] Matei | Weak solutions via Lagrange multipliers for contact models with normal compliance[END_REF].

Remark 1. Treating the model in the first approach we can prove the existence and the uniqueness of the weak solution. But, the approximation of the weak solution is based on a regularization/penalization technique. Treating the model in the second approach we are led to a generalized saddle point problem. Recall that, for weak formulations in Contact Mechanics via saddle point problems, efficient algorithms can be written in order to approximate the weak solution (see primal-dual active set strategies). But, there are a few open questions here:

the study of the uniqueness of the weak solution of the mixed variational formulation Problem 3; a priori error estimates; algorithms.

Boundary optimal control

Let us discuss in this section a boundary optimal control problem related to our contact problem.

For a fixed function f 0 ∈ L 2 (Ω) 3 , we consider the following state problem.

(PS1) Let f 2 ∈ L 2 (Γ 2 ) 3 (called control) be given. Find u ∈ V such that (Au, v -u) V + j(u, v) -j(u, u) ≥ Ω f 0 (x) • (v(x) -u(x)) dx (4)
+ Γ2 f 2 (x) • (v(x) -u(x)) dΓ for all v ∈ V.
According to Theorem 1, for every control f 2 ∈ L 2 (Γ 2 ) 3 , the state problem (PS1) has a unique solution u ∈ V, u = u(f 2 ). In addition, the following estimation takes place:

u V ≤ 1 m F ( f 0 L 2 (Ω) 3 + c 0 f 2 L 2 (Γ2) 3 ),
where m F is the constant in Assumption 1 and c 0 appears in [START_REF] Amassad | Optimal control of an elastic contact problem involving Tresca friction law[END_REF]. Now, we would like to act a control on Γ 2 such that the resulting stress σ be as close as possible to a given target

σ d = Fε(u d )
where u d is a given function.

Let Q ∞ be the real Banach space

Q ∞ = {F = F ijkl | F ijkl = F jikl = F klij ∈ L ∞ (Ω), 1 ≤ i, j, k, l ≤ 3}
endowed with the norm F ∞ = max 1≤i,j,k,l≤3 F ijkl L ∞ (Ω) . According to [START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF], page 97,

Fτ L 2 (Ω) 3×3 s ≤ 3 F ∞ τ L 2 (Ω) 3×3 s for all τ ∈ L 2 (Ω) 3×3 s . Therefore, σ -σ d L 2 (Ω) 3×3 s ≤ 3 max 1≤i,j,k,l≤ 3 F ijkl L ∞ (Ω) u -u d V .
Thus, σ and σ d will be close to one another if the difference between the functions u and u d is small in the sense of V -norm.

To give an example of a target of interest, u d , we can consider u d = 0. In this situation, by acting a control on Γ 2 , the tension σ is small in the sense of L 2 -norm, even if the volume forces f 0 does not vanish in Ω.

Let α, β > 0 be two positive constants and let us define the following functional

L : L 2 (Γ 2 ) 3 × V → R, L(f 2 , u) = α 2 u -u d 2 V + β 2 f 2 2 L 2 (Γ2) 3 .
Furthermore, we denote

V ad = {[u, f 2 ] | [u, f 2 ] ∈ V × L 2 (Γ 2 ), such that (4) is verified}. (POC1) Find [u * , f * 2 ] ∈ V ad such that L(f * 2 , u * ) = min [u,f 2 ]∈V ad L(f 2 , u) .
A solution of (POC1) is called an optimal pair. The second component of the optimal pair is called an optimal control. Theorem 3. Problem (P OC1) has at least one solution (u * , f * 2 ).

Let us fix ρ > 0 and f 0 ∈ L 2 (Ω) 3 .

We introduce the following regularized state problem.

(PS2) Let f 2 ∈ L 2 (Γ 2 ) 3 (called regularized control) be given. Find u ∈ V such that (Au, v -u) V + j ρ (u, v) -j ρ (u, u) ≥ (f 0 , v -u) L 2 (Ω) 3 (5) +(f 2 , v -u) L 2 (Γ2) 3 for all v ∈ V.
Herein, j ρ : V × V → R is defined as follows,

j ρ (u, v) = Γ3 p ρ ν (x, u ν (x) -g a (x))( (v ν (x)) 2 + ρ 2 -ρ)dΓ + Γ3 p ρ τ (x, u ν (x) -g a (x))( v τ (x) 2 + ρ 2 -ρ)dΓ for all u, v ∈ V,
where p ρ e , e ∈ {ν, τ }, satisfies the following assumptions. Notice that the functional j ρ (•, •) has the following properties:

• for all u, v ∈ V , j ρ (u, v) ≥ 0; j ρ (u, 0 V ) = 0; • for all u ∈ V, j ρ (u, •) : V → R is a convex and Gâteaux differentiable functional; • j ρ (η 1 , v 2 )-j ρ (η 1 , v 1 )+j ρ (η 2 , v 1 )-j ρ (η 2 , v 2 ) ≤ c 2 0 (L ν +L τ ) η 1 -η 2 V v 1 - v 2 V for all η 1 , η 2 , v 1 , v 2 ∈ V.
• for all u, v ∈ V , there exists

∇ 2 j ρ (u, v) ∈ V such that lim h→0 j ρ (u, v + hw) -j ρ (u, v) h = (∇ 2 j ρ (u, v), w) V for all w ∈ V. (∇ 2 j ρ (u, v), w) V = Γ3 p ρ ν (x, u ν (x) -g a (x)) v ν (x)w ν (x) (v ν (x)) 2 + ρ 2 dΓ + Γ3 p ρ τ (x, u ν (x) -g a (x)) v τ (x) • w τ (x) v τ (x) 2 + ρ 2 dΓ.
The regularized state problem has a unique solution u ρ ∈ V that depends Lipschitz continuously on f . This is a straightforward consequence of an abstract result in the theory of the quasivariational inequalities, see e.g. Theorem 3.7, in [START_REF] Sofonea | Variational Inequalities with Applications. A Study of Antiplane Frictional Contact Problems[END_REF].

For every

f 2 ∈ L 2 (Γ 2 ) 3 , the problem (PS2) has a unique solution u ∈ V, u = u(f 2 ). In addition, u V ≤ 1 m F ( f 0 L 2 (Ω) 3 + c 0 f 2 L 2 (Γ2) 3 ).
There exists an unique z ∈ V such that

(z, v) V = Ω f 0 • v dx for all v ∈ V.
Furthermore, there exists an unique y(f 2 ) ∈ V such that

(y(f 2 ), v) V = Γ2 f 2 • v dΓ for all v ∈ V.
Let u ∈ V be the unique solution of (PS2). Let us define

∂ 2 j ρ (u, u) = {ζ ∈ V | j ρ (u, v) -j ρ (u, u) ≥ (ζ, v -u) V for all v ∈ V }.
Therefore,

z + y(f 2 ) -Au ∈ ∂ 2 j ρ (u, u).
Since j ρ (•, •) is convex and Gâteaux differentiable in the second argument, we can write ∂ 2 j ρ (u, u) = {∇ 2 j ρ (u, u)}.

Thus, we are led to the following operatorial equation

Au + ∇ 2 j ρ (u, u) = z + y(f 2 ).
Let us define the admissible set,

V ρ ad = {[u, f 2 ] | [u, f 2 ] ∈ V × L 2 (Γ 2 ) 3 , such that (5) is verified}.
Using the functional L, we introduce the regularized optimal control problem,

(POC2) Find [ū, f 2 ] ∈ V ρ ad such that L( f 2 , ū) = min [u,f 2 ]∈V ρ ad L(f 2 , u) .
Theorem 4. The problem (POC2) has at least one solution (ū, f 2 ).

A solution of (POC2) is called a regularized optimal pair and the second component f 2 is called a regularized optimal control.

The following result hold true.

Theorem 5. (An optimality condition) Any regularized optimal control f 2 verifies

f 2 = - 1 β γ(p( f 2 )),
where γ is the trace operator and p( f 2 ) is the unique solution of the variational equation

α(u( f 2 ) -u d , w) V = (p( f 2 ), Aw + D 2 2 j ρ (u( f 2 ), u( f 2 ))w) V for all w ∈ V, u( f 2 ) being the solution of (PS2) with f 2 = f 2 .
Herein, for all v ∈ V, writing u instead of u( f 2 ),

(D 2 2 j ρ (u, u)v, w) V = Γ3 ∂ 2 p ρ ν (x, u ν (x) -g a (x)) u ν (x)v ν (x)w ν (x) u ν (x) 2 + ρ 2 dΓ + Γ3 ∂ 2 p ρ τ (x, u ν (x) -g a (x)) u τ (x) • w τ (x)v ν (x) u τ (x) 2 + ρ 2 dΓ + Γ3 +p ρ ν (x, u ν (x) -g a (x)) v ν (x)w ν (x)ρ 2 (u ν (x) 2 + ρ 2 ) 3/2 dΓ + Γ3 p ρ τ (x, u ν (x) -g a (x)) v τ (x) • w τ (x)( u τ 2 + ρ 2 ) -(u τ • w τ )(u τ • v τ ) ( u τ (x) 2 + ρ 2 ) 3/2 dΓ.
The main tool in the proof of Theorem 5 is a Lions's Theorem, which we recall here for the convenience of the reader. Theorem 6. Let B be a Banach space, X and Y two reflexive Banach spaces. Let also be given two C 1 functions F : B × X → Y, L : B × X → R. We suppose that, for all β ∈ B, i) There exists a unique ũ(β) such that For the proof of Theorem 6 we refer to, e.g., [START_REF] Amassad | Optimal control of an elastic contact problem involving Tresca friction law[END_REF].

F (β, ũ(β)) = 0, ii) ∂ 2 F (β, ũ(β)) is an isomorphism from X onto Y. Then, J(β) = L(β, ũ (
Let us indicate in the last part of this section two convergence results. The first one involves the unique solution of the regularized state problem (PS2) and the unique solution of the state problem (PS1).

Theorem 7. Let ρ > 0, f 0 ∈ L 2 (Ω) 3 and f 2 ∈ L 2 (Γ 2 ) 3 be given. If u ρ , u ∈ V are the solutions of problems (P S2) and (P S1), respectively, then,

u ρ → u in V as ρ → 0.
Next, we have a convergence result involving the solutions of the problems (POC2) and (POC1). Theorems 3-8 are new results; their proofs will be published in [START_REF] Matei | Boundary optimal control for a frictional contact problem with normal compliance[END_REF].

Let us mention here some open questions:

• f 2 ρ → f * 2 in L 2 (Γ 2 )
3 as ρ → 0; • an optimality condition for (PS1); • to study the boundary optimal control of the model by means of the mixed variational formulation.

  , and let us define two bilinear forms a(•, •) and b(•, •) as follows:

Assumption 6 p ρ e : Γ 3 ×Assumption 8

 38 R → R + . The mapping x → p ρ e (x, r) is measurable on Γ 3 for any r ∈ R, and p ρ e (x, r) = 0 for all r ≤ 0, a.e. x ∈ Γ 3 . Assumption 7 p ρ e (x, •) ∈ C 1 (R) a.e. on x ∈ Γ 3 . There exists M e > 0 such that |p ρ e (x, r)| ≤ M e for all r ∈ R, a.e. x ∈ Γ 3 . In addition |∂ 2 p ρ e (x, r)| ≤ L e for all r ∈ R, a.e. x ∈ Γ 3 . There exists G e : R + → R + (e ∈ {ν, τ )} such that |p ρ e (x, r)p e (x, r)| ≤ G e (ρ) for all r ∈ R, a.e. x ∈ Γ 3 and lim ρ→0 G e (ρ) = 0.

∂ 2 F

 2 β)) is differentiable and, for every ζ ∈ B, dJ dβ (β)ζ = ∂ 1 L(β, ũ(β))ζ -p(β), ∂ 1 F (β, ũ(β))ζ Y ,Y ,where p(β) ∈ Y is the adjoint state, unique solution of (β, ũ(β)) * p(β) = ∂ 2 L(β, ũ(β)) in X .

Theorem 8 .

 8 Let [ū ρ , f 2 ρ ] be a solution of the problem (POC2). Then, there exists a solution of the problem (POC1),[u * , f * 2 ] , such that ūρ → u * in V as ρ → 0, f 2 ρ f * 2 in L 2 (Γ 2 ) 3 as ρ → 0.

This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, project number PN-II-ID-PCE-2011-3-0257, and by LEA Math Mode CNRS-IMAR.


Nina N Subbotina

Timofey B Tokmantsev

Eugenii A Krupennikov

Dynamic Programming to Reconstruction Problems for a Macroeconomic Model

Keywords: inverse problem, positive and negative discrepancy, optimal control problem, dynamic programming

Perturbed inverse reconstruction problems for controlled dynamic systems are under consideration. A sample history of the actual trajectory is known. This trajectory is generated by a control, which isn't known. Moreover, the deviation of the samples from the actual trajectory satisfies the known estimate of the sample error. The inverse problem with perturbed (inaccurate) sample of trajectory consists of reconstructing trajectories which are close to the actual trajectory in C. Controls generating the trajectories should be close in L2 to the normal control generating the actual trajectory and have the least norm in L2. A numerical method for solving this problem is suggested. The application of the suggested method is illustrated by the graphics.

Introduction

A model of macroeconomics is considered due to works by E.G. Al'brekht [START_REF] Al | Methods of construction and identification of mathematical models for macroeconomic processes // electronic journal[END_REF]. The model has the form of two nonlinear ordinary differential equations. The right-hand sides of the equations depend on control parameters. The rate of taxation, the refunding rate and the currency exchange course are included in control parameters because they determine economical conditions for production activity. A sample history of the actual trajectory of the model is known. A numerical method is suggested and verified to reconstruct the actual trajectory and the control generating it. It is based on the method of the dynamic programming. Results of numerical calculations of the solutions of the inverse problem are represented for statistic data obtained from a reports of companies sent to local statistic offices in Russia.

Macroeconomic Model

Consider a model of a macroeconomic system, where the symbol x 1 denotes the gross product, x 2 denotes production costs, G denotes profits.

Let dynamics of x 1 (t), x 2 (t) be of the form dx1 dt = u 1 (t) ∂G(x1,x2) ∂x1 , dx2 dt = u 2 (t) ∂G(x1,x2) ∂x2 [START_REF] Al | Methods of construction and identification of mathematical models for macroeconomic processes // electronic journal[END_REF] on a time interval t ∈ [0, T ]. Here u 1 (t), u 2 (t) are control parameters, satisfying the geometric restrictions

|u 1 | ≤ U 1 , |u 2 | ≤ U 2 , (2)
where U 1 > 0, U 2 > 0 are constants.

Known Data

We have got the following statistic data in the form of a table of parameters

x * 1 , x * 2 , G * measured at given instants t i , t i = 0, 1, . . . , N , t 0 = 0, t N = T, x * 1 (t 0 ), x * 1 (t 1 ), . . . , x * 1 (t N ), x * 2 (t 0 ), x * 2 (t 1 ), . . . , x * 2 (t N ), G * (t 0 ), G * (t 1 ), . . . , G * (t N ), where x * 1 (t i ) , x * 2 (t i ) are measurements of the actual trajectory x 1 * (•), x 2 * (•) of the system (1) on the interval [0, T ].

Hypothesis

Following the Albrekht's works, we assume that the mathematical model of the measured dynamics meets the following assertions

• the structure of the function G(x 1 , x 2 ) has the form of the polynomial

G(x 1 , x 2 ) = x 1 x 2 (a 0 + a 1 x 1 + a 2 x 2 ). (3)
• the given statistic data are measurements of the actual trajectory x * (•) = (x 1 * (•), x 2 * (•)) and profit function G(x * (•)) with errors, while estimate δ on admissible errors is known.

|x 1 * (t i ) -x * 1 (t i )| ≤ δ, |x 2 * (t i ) -x * 2 (t i )| ≤ δ, |G(x 1 * (t i ), x 2 * (t i )) -G * (t i )| ≤ δ, i = 0, 1, . . . , N. (4)
• such smooth continuous interpolations

y(•) = (y 1 (•), y 2 (•)) of the data x * (t i ) = (x * 1 (t i ), x * 2 (t i )), i = 0, 1, . . . , N are defined, that d 2 y i (t) dt 2 ≤ K, K > 0, t ∈ [0, T ], i = 1, 2, . (5)
y(•) -x * (•) c → 0, as δ → 0. (6)

Reconstruction Problems

The inverse problems are identification problem and reconstruction problem for the model, which supposes reconstructing such trajectories x δ (•) of system (1) generated by measurable controls u δ (•), satisfying (2), that

x δ (•) -x * (•) C = max t∈[0,T ]
x δ (t) -x * (t) → 0, as δ → 0;

u δ (•) -u * (•) 2 L2 = T 0 u δ (t) -u * (t) 2 dt → 0, as δ → 0; where x * (•) = (x 1 * (•), x 2 * (•)) is the actual trajectory on [0, T ] generated by "normal" control u * (•) = (u 1 * (•), u 2 * (•)), which has the minimal norm in L 2 ([0, T ], R 2 ).
The method suggested below is based on the dynamic programming [START_REF] Bellman | Dynamic Programming[END_REF] for auxiliary optimal control problems. It can be interpreted as a modification of Tikhonov method [START_REF] Tihonov | Solution of incorrectly formulated problems and the regularization method //Soviet Math[END_REF].

The other approach to solutions of the inverse problems with the help of optimal feedbacks [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF] in auxiliary optimal control problems was suggested in works by Osipov, Kryazhimskii [START_REF] Yu | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF].

6 Identification Problem for the Function G(x 1 , x 2 )

At first we consider the identification problem for parameters a 0 , a 1 , a 2 of the polynomial

G(x 1 , x 2 ) = x 1 x 2 (a 0 + a 1 x 1 + a 2 x 2 )
to obtain the best correspondence with the given statistic materials.

In order to do this, we apply the least square method to the statistic data

N i=0 [G * (t i ) -G(x * 1 (t i ), x * 2 (t i ))] 2 -→ min (a0,a1 ,a2)
.

Auxiliary Optimal Control Problems (AOCPs)

We introduce the following AOCPs to solve the reconstruction problem. Consider dynamics of the form

dx1 dt = u 1 ∂G(x1,x2) ∂x1 , dx2 dt = u 2 ∂G(x2,x2) ∂x2 , t ∈ [0, T ], u = (u 1 , u 2 ) ∈ P, (7)
P = {|u 1 | ≤ U 1 , |u 2 | ≤ U 2 , }. (8)
The set of admissible controls is defined as

U [0,T ] = {∀u(•) : [0, T ] → P -measurable}.
We introduce the α-regularized positive discrepancy functional

I + 0,x 0 1 ,x 0 2 (u(•)) = T 0 [(y 1 (t) -x 1 (t)) 2 + (y 2 (t) -x 2 (t)) 2 ] 2 + (9) α 2 (u 2 1 (t) + u 2 2 (t)) 2 dt,
where α is a small parameter. The functions y 1 (•), y 2 (•) are interpolations of statistic data.

We also consider the α-regularized negative discrepancy functional

I - 0,x 0 1 ,x 0 2 (u(•)) = T 0 - [(y 1 (t) -x 1 (t)) 2 + (y 2 (t) -x 2 (t)) 2 ] 2 + ( 10
)
α 2 (u 2 1 (t) + u 2 2 (t)) 2
dt.

Optimal Results in AOCPs

Let small parameters α > 0, δ > 0 be fixed and interpolations y 1 (•), y 2 (•) of the statistic data be known. The aim of the AOCPs at an initial state t = 0,

x 1 (0) = x 0 1 , x 2 (0) = x 0 2
is to minimize the cost functionals ( 10), ( 9) under the condition

x 1 (T ) = y 1 (T ), x 2 (T ) = y 2 (T ). (11)
The optimal results in the class U [0,T ] are equal to

V ± (0, x 0 1 , x 0 2 ) = inf u(•)∈U [0,T ] I ± t0,x 0 1 ,x 0 2 u(•) . (12)

Hamiltonian

Let's consider the AOCP for the negative discrepancy functional [START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF]. Let us denote

ω 1 (x) = ω 1 (x 1 , x 2 ) = ∂G(x 1 (t), x 2 (t)) ∂x 1 = a 0 x 2 + 2a 1 x 1 x 2 + a 2 x 2 2 , ω 2 (x) = ω 2 (x 1 , x 2 ) = ∂G(x 1 (t), x 2 (t)) ∂x 2 = a 0 x 1 + a 1 x 2 1 + 2a 2 x 1 x 2 . H α (t, x 1 , x 2 , s 1 , s 2 ) = min u∈P s 1 u 1 ω 1 (x 1 , x 2 ) + s 2 u 2 ω 2 (x 1 , x 2 )+ + α 2 (u 2 1 + u 2 2 ) 2 - (x 1 -y 1 (t)) 2 + (x 2 -y 2 (t)) 2 2 = = s 1 u 0 1 ω 1 (x 1 , x 2 ) + s 2 u 0 2 ω 2 (x 1 , x 2 ) + α 2 (u 0 1 2 + u 0 2 2 ) 2 (13) - (x 1 -y 1 (t)) 2 + (x 2 -y 2 (t)) 2 2 .
where for i = 1, 2,

u 0 i (x, s) =      -U i , if -siωi(x(t)) α 2 ≤ -U i , -s1ωi(x(t)) α 2 , if -s1ωi(x(t)) α 2 ∈ [-U i , U i ], U i , if -s1ωi(x(t)) α 2 ≥ U i .
So, for the simple case

u 0 i (x, s) ∈ [-U i , U i ], i = 1, 2, (14)
we get Hamiltonian of the form

H α (t, x 1 , x 2 , s 1 , s 2 ) = - 1 2α 2 (s 2 1 + s 2 2 ) - (x 1 -y 1 (t)) 2 + (x 2 -y 2 (t)) 2 2 .

Characteristics

Necessary optimality conditions for the AOCPs has the following form [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF][START_REF] Clarke | Optimization and nonsmooth analysis[END_REF]: the characteristic system

dx i dt = ∂H α (t, x, s) ∂s i , ds i dt = - ∂H α (t, x, s) ∂x i , i = 1, 2, t ∈ [0, T ], (15)
and the boundary conditions

x i (T ) = y i (T ), s i (T ) = ξ i , ω 2 i (x(T ))ξ i α 2 -ẏi (T ) ≤ δ, i = 1, 2. (16)

Characteristics for the Simple Case

Restrictions U 1 , U 2 for admissible controls are usually unknown. To simplify the explanations we assume that U 1 , U 2 are large enough to let interpolations y(t) provide the simple case(14) with boundary conditions (16).

The characteristic system for the simple case has the form:

dx 1 (t) dt = - s 1 (t) α 2 ω 2 1 (x 1 (t), x 2 (t)), (17)
dx 2 (t) dt = - s 2 (t) α 2 ω 2 2 (x 1 (t), x 2 (t)), ds 1 (t) dt = x 1 (t) -y 1 (t)+ + s 2 1 (t) α 2 F 1 (x 1 (t), x 2 (t)) + s 2 2 (t) α 2 F 2 (x 1 (t), x 2 (t)), ds 2 (t) dt = x 2 (t) -y 2 (t)+ + s 2 1 (t) α 2 F 3 (x 1 (t), x 2 (t)) + s 2 2 (t) α 2 F 4 (x 1 (t), x 2 (t)),
where

F 1 (x 1 , x 2 ) = 2a 1 x 2 (a 0 x 2 + 2a 1 x 1 x 2 + a 2 x 2 2 ), F 2 (x 1 , x 2 ) = (a 0 + 2a 1 x 1 + 2a 2 x 2 )(a 0 x 1 + a 1 x 2 1 + 2a 2 x 1 x 2 ), F 3 (x 1 , x 2 ) = (a 0 + 2a 1 x 1 + 2a 2 x 2 )(a 0 x 2 + 2a 1 x 1 x 2 + a 2 x 2 2 ), F 4 (x 1 , x 2 ) = 2a 2 x 1 (a 0 x 1 + a 1 x 2 1 + 2a 2 x 1 x 2 ),
boundary conditions

x 1 (T ) = y 1 (T ), x 2 (T ) = y 2 (T ), (18)
ξ - 1 ≤ s 1 (T ) = ξ 1 ≤ ξ + 1 , ξ - 2 ≤ s 2 (T ) = ξ 2 ≤ ξ + 2 ,
where

ξ - i = - ẏi (T )α 2 ω i (y(T )) 2 - δα 2 ω i (y(T )) 2 , ξ + i = - ẏi (T )α 2 ω i (y(T )) 2 + δα 2 ω i (y(T )) 2 , i = 1, 2. (19)

Solutions of Inverse Problems

Let us pick such characteristics (15)-(19) x α δ (•) and the realizations of extremal feedbacks u α δ [t] = u α (t, x α δ (t), s α δ (t)), generating them, which satisfy the relations:

x(0, ξ) -y(0) ≤ α + δ, I ± 0,x α δ (0) (u α δ [•]) = min
x(0,ξ)-y(0) ≤α+δ

I ± 0,x(0) (u α (•)) = V ± (0, x α δ (0)), (20)
u α (t) = u α (t, x(t, ξ), s(t, ξ)), t ∈ [0, T ].
We have got that these characteristics x α δ (•, ξ) and controls u α δ [•], generating them, provide solutions to the inverse problems [START_REF] Subbotina | The method of characteristics for Hamilton-Jacobi-Bellman equation[END_REF][START_REF] Krupennikov | Validation of a solution method for the problem of reconstructing the dynamics of a mecroeconomic system // Proceedings of the Institute of Mathematics and Mechanics[END_REF][START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF].

Assumptions

A1 Such constants α 0 > 0, δ 0 > 0 exist that state characteristics x 1 (t, ξ) and x 2 (t, ξ) of the form (15)-( 19) for all t ∈ [0, T ] belong to the compact set Φ:

Φ ⊃ Φ(δ, α) ∀ δ, α : 0 < δ ≤ δ 0 , 0 < α ≤ α 0 , (21)
Φ(δ, α) = (t, x) : t ∈ [0, T ], x = x(t, ξ), x(T, ξ) = y(T ), ω 2 i (x(T ))ξ i α 2 -ẏi (T ) ≤ δ, i = 1, 2 . A2 For (x 1 , x 2 ) ∈ Φ such constants ω i > 0, ωi > 0, i = 1, 2 exist, that 0 < ω 2 1 ≤ ω 2 1 (x 1 (t), x 2 (t)) ≤ ω2 1 , 0 < ω 2 2 ≤ ω 2 2 (x 1 (t), x 2 (t)) ≤ ω2 2 , t ∈ [0, T ].

Note

In the example below, one can choose such α 0 > 0, δ 0 > 0, r > 0, that

Φ = Φ r = {(t, x) : t ∈ [0, T ], x -y(t) ≤ r}, min 0≤t≤T y i (t) > 3r > 0, i = 1, 2,
and assumptions A1-A2 are true.

Main Result

Let us consider AOCPs for the system ( 7), ( 8) at initial states

x(0) ∈ {x : x -y(0) ≤ δ + α}
where the aim is to reach the target set {T, x = y(T )} and minimize the functional [START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF].

The following assertions are proven [START_REF] Krupennikov | Validation of a solution method for the problem of reconstructing the dynamics of a mecroeconomic system // Proceedings of the Institute of Mathematics and Mechanics[END_REF][START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF].

Lemma 1. Let x α δ (t) be a solution of the AOCP (7), ( 8), [START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF]. Let u α δ (t) be a control generating x α δ (t). If conditions A1-A2 are true in the problem, then such constant c > 0 exists that the following estimate takes place:

I 0,x α δ (0) (u α δ (•)) ≤ I 0,x * (0) (u * (•)) + ζ(α, δ), ζ(α, δ) = cδ(δ 2 + α 2 U 2 * ),
where

U * = max{U 1 , U 2 }.
We introduce the functions

φ(α, δ, h) = T M h T M h 2 + 2δ + α + ζ(α, δ) , ρ(h) = nU * T (K + M )h,
where K, M are constant parameters. Let us denote numerical approximations of the solution x α δ (•), u α δ (•) of AOCP ( 7), ( 8), [START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF]

as x h (•), u h (•).
Theorem 1. Let conditions A1 -A2 be true in AOCP [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], ( 8), [START_REF] Subbotina | A Study of the Stability of Solutions to Inverse Problems of Dynamics of Control Systems under Perturbations of Initial Data[END_REF]. Then there exists such constants M > 0, K > 0 and parameters h = h(δ) > 0, α = α(δ) > 0, δ > 0, satisfying the conditions lim

δ→0 h(δ) = 0, lim δ→0 α(δ) = 0, lim δ→0 2 α 2 φ(α, δ, h) + ρ(h) + T 2 (M h + α + 2δ + ζ(α, δ)) 2 = 0, (22)
that the following relations are true

lim δ→0 x h(δ) (•) -x * (•) C = 0, lim δ→0 u h(δ) (•) -u * (•) L2 = 0.

Numerical experiments

Results of application of the suggested numerical method via AOCP with the functional I - (0,x(0)) (u(•)) are exposed on the figures 1-6 below. Note that the results obtained via AOCP with the functional I + (0,x(0)) (u(•)) are not so satisfying (see figures [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF][START_REF] Subbotina | The method of characteristics for Hamilton-Jacobi-Bellman equation[END_REF]. This is because of the properties of characteristics in the considered AOCPs. We used the data on the industry of the Ural Region in Russia for the period 1970-1985 (10000 Rubles = 1) due to paper [

Perspectives

The suggested numerical method can be applied in the following directions.

-Identification and reconstruction of dynamic models of production activity for single firms, various branches of industry or industry and economics of a region. -Investigating properties of the examined object.

-A short-term and long-term prediction and analysis of scenarios of the process development in the future. -Analysis of the production plan and construction of feedback controls realizing the plan.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 12345678 Fig. 1. Trajectory bundle obtained with α 2 = 10 -4 , t ∈ [1, 1.5]

  1]:

	t	Year Gross regional 1 product x *	Costs x * 2 Profit G *
	0	1970 37.88	21.69	6.17
	1	1971 40.63	23.70	6.31
	2	1972 43.25	25.45	6.68
	3	1973 46.00	27.30	6.98
	4	1974 49.33	29.44	7.04
	5	1975 53.04	32.16	7.27
	6	1976 57.03	35.01	7.62
	7	1977 59.85	36.92	8.00
	8	1978 62.72	38.69	8.27
	9	1979 63.45	38.76	8.42
	10 1980 65.74	39.96	8.61
	11 1981 65.90	39.75	8.21
	12 1982 69.22	41.31	9.65
	13 1983 64.52	37.86	9.28
	14 1984 71.03	42.04	10.26
	15 1985 74.69	45.05	10.76

Acknowledgements. The work is supported partially by RFBR (projects No. 14-01-00168, No. 14-01-00486).


Vilmos Komornik
email: komornik@math.unistra.fr

Paola Loreti
email: paola.loreti@sbai.uniroma1.it

Michel Mehrenberger

Observability of a ring shaped membrane via Fourier series

We study the inverse Ingham type inequality for a wave equation in a ring. This leads to a conjecture on the zeros of Bessel cross product functions. We motivate the validity of the conjecture through numerical results. We do a complete analysis in the particular case of radial initial data, where an improved time of observability is available.

Introduction

We study the solution u = u(t, r, θ) of the wave equation in an annulus Ω of small radius a and big radius b:    u = ∆u, 0 < t < T, x ∈ Ω u = 0, 0 < t < T, x ∈ Γ = ∂Ω u(0, x) = u 0 (x), ∂ t u(0, x) = u 1 (x), x ∈ Ω.

(1)

a b Ω
We have the following classical observability estimate (see [START_REF] Komornik | Fourier Series in Control Theory[END_REF] for example): Proposition 1. There exists T 0 > 0 such that for T > T 0 the system (1) is observable: there exists a constant c > 0 such that we have

u 0 2 H 1 0 (Ω) + u 1 2 L 2 (Ω) ≤ c T 0 Γ |∂ ν u(t, x)| 2 dΓ dt, (2)
for all (u 0 , u 1 ) ∈ H 1 0 (Ω) × L 2 (Ω).

This proposition is a special case of a general result, which can be proved by micro-local analysis [START_REF] Bardos | Sharp suffcient conditions for the observation, control and stabilization of waves from the boundary[END_REF] with the critical time T 0 = 2 √ b 2 -a 2 , according to the geometrical ray condition.

2 √ b 2 -a 2
Our aim here is to study this problem using a Fourier series approach [START_REF] Ingham | Some trigonometrical inequalities with applications in the theory of series[END_REF][START_REF] Haraux | Séries lacunaires et contrôle semi-interne des vibrations d'une plaque rectangulaire[END_REF] Note that such an approach has been tackled for the whole disk (see [START_REF] Komornik | Fourier Series in Control Theory[END_REF]); the case of an annulus leads to new difficulties. In particular, we are not able to treat the general case, as it relies on very precise estimates of the location of the zeros of Bessel cross product functions, that we did not find in the literature. Instead, we state a conjecture on these zeros and give some numerical results to support the conjecture. In the case of radial functions or more generally functions with a limited number of modes in the angle direction, we are able to get the observability estimates with this method, even for smaller times T > 2(b -a), using an asymptotic result of MacMahon (see [START_REF]Handbook of mathematical functions[END_REF], p. 374). In Section 2, we give the expression of the solution. In Section 3, we formulate a theorem for the special case of radial functions. Section 4 is devoted to the statement of the conjecture and its numerical illustration. Finally, we prove in Section 5 the theorem of Section 3.

Expression of the solution

Let J ν (resp. Y ν ) be the Bessel functions of first (resp. second) kind of order ν. We recall the following proposition (see [START_REF] Komornik | Fourier Series in Control Theory[END_REF][START_REF] Lebedev | Uflyand Worked problems in applied mathematics[END_REF])

Proposition 2. (i) Let ν ∈ R * and 0 < α < 1. The nonnegative zeros γ ν,k , k ∈ N * of the Bessel cross product function H ν,α (x) = Y ν (x)J ν (αx) -J ν (x)Y ν (αx),
form a strictly increasing sequence:

0 < γ ν,α,1 < γ ν,α,2 < . . . γ ν,α,k < . . . , k ∈ N * .
(ii) The eigenfunctions of the Laplacian corresponding to (1) are R k,m (r)e imθ for k ∈ N * and m ∈ N, where

R k,m (r) = Y m (γ m,α,k )J m ( r b γ m,α,k ) -J m (γ m,α,k )Y m ( r b γ m,α,k ). (3)
(iii) For a dense set of initial data the solution of (1) is given by the formula

u(t, r, θ) = ∞ m=0 ∞ k=1 R k,m (r) c + k,m e imθ e i γ m,α,k b t + c - k,m e -imθ e -i γ m,α,k b t (4)
with complex coefficients c ± k,m , all but finitely of which vanish.

Observability estimates for radial solutions

We recall from [START_REF]Handbook of mathematical functions[END_REF] the following estimate of MacMahon (1894):

γ ν,α,k = πk 1 -α + 4ν 2 -1 8πα 1 -α k + O 1 -α k 3 . (5)
Thanks to this estimation, we can obtain the following theorem.

Theorem 1. Let T > 2(b -a).
For each positive integer M there exists a constant c M > 0 such that

u 0 2 H 1 0 (Ω) + u 1 2 L 2 (Ω) ≤ c M T 0 Γ |∂ ν u(t, x)| 2 dΓ dt, (6)
for all solutions of (1) of the form (4)

with c ± k,m = 0 whenever m ≥ M . Remark 1.
-This theorem covers the case of radial initial data, corresponding to the case

M = 1. -If 2 √ b 2 -a 2 > T > 2(b -a), then the constant c M tends to zero as M → ∞ because T 0 = 2 √ b 2 -a 2
is the critical observability time for general initial data.

-The theorem remains true if the integral in ( 6) is taken only over the outer boundary (the circle of radius b): see the estimate (16)) below. -We may obtain similar results by changing the boundary conditions. For example, we may take homogeneous Neumann condition on the inner boundary, and observe the solution only on the outer boundary (see [START_REF] Gottlieb | Eigenvalues of the Laplacian with Neumann boundary conditions[END_REF] for the corresponding asymptotic gap estimate that is needed).

Conjecture and numerical illustration

We then state the following conjecture.

Conjecture 1 Let 0 < α < 1 and ν > 1/2.
-There exists a positive integer

k ν (α) such that γ ν,α,k+1 -γ ν,α,k is decreasing for k ≤ k ν (α) and increasing for k ≥ k ν (α). -We have γ ν,α,kν (α)+1 -γ ν,α,kν (α) ≥ π 1 -α 2 -1/2 .
Using a Bessel Zeros Computer [3] we can evaluate the zeros γ ν,α,k for several parameters ν, α and k. We plot on Figure 1 (top and bottom left) the values k ν (α) versus ν for different values of α. We observe that k ν (α) increases with ν for a fixed α. The dependance seems to be almost linear, except for small values of α where the dependence seems to be quadratic (see Figure 1, bottom left). On Figure 1 (bottom right), we see the relative difference with the gap π √ 1-α 2 (corresponding to the optimal value of Proposition 1). We observe that γ ν,α,kν (α)+1 -γ ν,α,kν (α) decreases and approaches to this gap as ν increases.

Proof of Theorem 1

We first express the norm in terms of the Fourier coefficients. We have

u(t, r, θ) = ∞ m=0 ∞ k=1 R k,m (r) c + k,m e imθ e i γ m,α,k b t + c - k,m e -imθ e -i γ m,α,k b t with R k,m (r) = Y m (γ m,α,k )J m ( r b γ m,α,k ) -J m (γ m,α,k )Y m ( r b γ m,α,k ).
This leads to the equalities

u 0 (r, θ) = u(0, r, θ) = ∞ m=0 ∞ k=1 R k,m (r) c + k,m e imθ + c - k,m e -imθ .
and

u 1 (r, θ) = ∂ t u(0, r, θ) = ∞ m=0 ∞ k=1 iγ m,α,k b R k,m (r)e imθ c + k,m e imθ -c - k,m e -imθ .
Now, using the orthogonality of the eigenvectors of the Laplacian operator we obtain

u 1 2 L 2 (Ω) = b r=a 2π θ=0 |u 1 (r, θ)| 2 rdrdθ = 2π ∞ m=1 ∞ k=1 γ m,α,k b 2 c + k,m 2 + c - k,m 2 b a |R k,m (r)| 2 rdr + 2π ∞ k=1 γ 0,α,k b 2 c + k,0 -c - k,0 2 b a |R k,0 (r)| 2 rdr,
and a similar computation gives

Ω |u 0 | 2 dΩ = 2π ∞ m=1 ∞ k=1 c + k,m 2 + c - k,m 2 b a |R k,m (r)| 2 rdr + 2π ∞ k=1 c + k,0 + c - k,0 2 b a |R k,0 (r)| 2 rdr. Since u 0 2 H 1 0 (Ω) = Ω |∇u 0 | 2 dΩ + Ω |u 0 | 2 dΩ,
we have to compute also the the first integral on the right side. We have

Ω |∇u 0 | 2 dΩ = b r=a 2π θ=0 |∂ r u 0 | 2 rdrdθ + b r=a 2π θ=0 1 r |∂ θ u 0 | 2 drdθ. Since ∂ r u 0 (r, θ) = ∂ r u(0, r, θ) = ∞ m=0 ∞ k=1 R k,m (r) c + k,m e imθ + c - k,m e -imθ
and

∂ θ u 0 (r, θ) = ∂ θ u(0, r, θ) = ∞ m=0 ∞ k=1 imR k,m (r) c + k,m e imθ -c - k,m e -imθ ,
using the orthogonality of the eigenfunctions, we get

b r=a 2π θ=0 |∂ r u 0 | 2 rdrdθ = 2π ∞ m=1 ∞ k=1 c + k,m 2 + c - k,m 2 b a R k,m (r) 2 rdr, + 2π ∞ k=1 c + k,0 + c - k,0 2 b a R k,0 (r) 2 rdr and b r=a 2π θ=0 1 r |∂ θ u 0 | 2 drdθ = 2π ∞ m=1 ∞ k=1 c + k,m 2 + c - k,m 2 m 2 b a 1 r |R k,m (r)| 2 dr
Using all these results we obtain the equalities

E 0 := u 0 2 H 1 0 (Ω) + u 1 2 L 2 (Ω) = 2π ∞ m=1 ∞ k=1 c + k,m 2 + c - k,m 2 b a 1 + γ m,α,k b 2 + m 2 r 2 |R k,m (r)| 2 rdr + b a R k,m (r) 2 rdr + 2π ∞ k=1 ( b a c + k,0 + c - k,0 2 + c + k,0 -c - k,0 2 γ 0,α,k b 2 |R k,0 (r)| 2 rdr + c + k,0 + c - k,0 2 b a R k,0 (r) 2 rdr) (7)
and

Γ |∂ ν u(t, x)| 2 dΓ = ∈{a,b} 2π 0 |∂ r u(t, , θ)| 2 dθ = 2π ∈{a,b} ∞ m=1   ∞ k=1 c + k,m R k,m ( )e i γ m,α,k b t 2 + ∞ k=1 c - k,m R k,m ( )e -i γ m,α,k b t 2   + 2π ∈{a,b}   ∞ k=1 c + k,0 R k,0 ( )e i γ 0,α,k b t + ∞ k=1 c - k,0 R k,0 ( )e -i γ 0,α,k b t 2   Since asymptotic gap is 1 b π 1-α = π
b-a by the formula (5), we may apply Ingham's theorem (more precisely its version due to Haraux, see, e.g., [START_REF] Komornik | Fourier Series in Control Theory[END_REF]) to deduce from the last equality the existence of C 1,M > 0 such that

Γ |∂ ν u(t, x)| 2 dΓ ≥ C 1,M ∞ m=0 ∞ k=1 c + k,m 2 + c - k,m 2 ∈{a,b} R k,m ( ) 2 ,
for all complex sequences (c ± k,m ) with c ± k,m = 0 whenever m ≥ M . Now it remains to prove that there exists another constant C 2,M > 0 such that

∈{a,b} R k,m ( ) 2 ≥ C 2,M b a 1 + γ m,α,k b 2 + m 2 r 2 |R k,m (r)| 2 rdr + b a R k,m (r) 2 rdr , (8)
for m ≥ 1, and a similar inequality for m = 0.

We adapt an argument used in [START_REF] Komornik | Fourier Series in Control Theory[END_REF], p. 107. Let y satisfy the Bessel equation

x 2 y + xy + (x 2 -m 2 )y = 0. ( 9
)
Let c > 0. Multiplying the equation by 2y and integrating over (ca, cb) (instead of (0, c) as in [START_REF] Komornik | Fourier Series in Control Theory[END_REF]), we get

Recall that

R k,m (r) = Y m (γ m,α,k )J m ( r b γ m,α,k ) -J m (γ m,α,k )Y m ( r b γ m,α,k ). Now, we define y k,m (x) = Y m (γ m,α,k )J m (x) -J m (γ m,α,k )Y m (x)
which satisfies (9), and thus we have (10) with y = y k,m and c = γ m,α,k b . We have

y(cr) = R k,m (r), y(ca) = R k,m (a) = 0, y(cb) = R k,m (b) = 0, cy (cr) = R k,m (r) and thus 2c 2 b a rR k,m (r) 2 dr = b 2 R k,m (b) 2 -a 2 R k,m (a) 2 .
(11)

From (3) we have the relation

r 2 R k,m (r) + rR k,m (r) -m 2 R k,m (r) = -c 2 r 2 R k,m (r). (12)
Integrating by parts and using the relations

R k,m (b) = R k,m (a) = 0 we obtain that b a R k,m (r) 2 rdr = - b a rR k,m (r) R k,m (r) + 1 r R k,m (r) dr.
Using the relation (12) hence we conclude that

b a R k,m (r) 2 rdr = - b a rR k,m (r) 2 -c 2 + m 2 r 2 dr. (13)
Setting 5) and the inequalities γ k,α,m > 0, we conclude that

A = b a 1 + c 2 + m 2 r 2 |R k,m (r)
A ≤ C 3,M R k,m (b) 2 ≤ C 3,M R k,m (b) 2 + R k,m (a) 2 (16)
for a suitable constant C 3,M > 0. This proves (8) for m ≥ 1.

For m = 0, using (13) we deduce from (7) that

Fig. 2 - 1

 21 Fig. Index kν (α) where minimal gap occurs versus ν ≥ 1/2 (top left, logarithmic plot for α ∈ {0.001, 0.1, 0.2, 0.3, 0.5, 0.8}; top right, standard plot for α = 0.1 and α = 0.2; bottom left, α = 0.001 and comparison with function 0.25ν 2 ). Difference γ ν,α,kν (α)+1 -γ ν,α,kν (α) π/ √ 1-α 2

2 2 b a ry 2

 222 y y + 2x(y ) 2 + 2(x 2 -m 2 )yy )dx = 0.Integrating by parts the first and third terms, we obtaincb ca (-2x(y ) 2 + 2x(y ) 2 + 2xy 2 )dx = [x 2 y 2 + (x 2 -m 2 )y 2 ] cb ca , dx = [x 2 y 2 + (x 2 -m 2 )y 2 ] cb ca .The change of variable x = cr transforms this into 2c (cr)dr = [x 2 y 2 + (x 2 -m 2 )y 2 ] cb ca .

b 2

 2 R k,m (b) 2 -a 2 R k,m (a) 2 .

  and we conclude as before.


I R Rȃdulescu

D Cândea

A Halanay

A complex mathematical model with competition in leukemia with immune responsean optimal control approach

Keywords: leukemia, asymmetric division, competition, optimal control, treatment

This paper investigates an optimal control problem associated with a complex nonlinear system of multiple delay differential equations modeling the development of healthy and leukemic cell populations incorporating the immune system. The model takes into account space competition between normal cells and leukemic cells at two phases of the development of hematopoietic cells. The control problem consists in optimizing the treatment effect while minimizing the side effects. The Pontryagin minimum principle is applied and important conclusions about the character of the optimal therapy strategy are drawn.

Introduction

Leukemia is a cancer of the blood and bone marrow, characterized by large and uncontrolled growth of white blood cells. The most studied type of leukemia, Chronic myelogenous leukemia (CML), involves granular leukocyte precursors, namely the myelocyte line. The trigger of CML is a chromosomal abnormality, called the Philadelphia chromosome, that occurs in all cell lineages in about 90% of cases. The product of this chromosome is the formation of the BCR-ABl fusion protein which is thought to be responsible for the dysfunctional regulation of myelocyte proliferation and other features of CML. The standard treatment of CML in recent years is Imatinib, a molecular targeted drug ( [START_REF] Marley | Chronic myeloid leukemia:stem cell derived but progenitor cell driven[END_REF]), that has the effect that almost all patients attain hematological remission ( [START_REF] Lee | Chronic myelogenous leukemia[END_REF]) and 75% attain cytogenetic remission.

Nowadays, it is well known that both innate and adaptive immunity are implicated in the defense mechanisms against cancer and recent progress in cancer immunology suggest that the immune system plays a fundamental role in tumor progression [START_REF] Topalian | Cancer Immunotherapy Comes of Age[END_REF]. In CML, the biological literature reveals that T cells may play an important role in stemming the expansion of leukemic cells. This happens because leukemic cells express antigens that are immunogenic and can be recognized by cytotoxic T cells (CD8+ T cells or CTLs). In the paper [START_REF] Riether | Regulation of hematopoietic and leukemic stem cells by the immune system[END_REF], the authors found that leukemia-specific effectors CTLs were able to eliminate leukemic stem cells (LSCs) in vitro and in vivo in a setting with minimal leukemia load. The role of CD4+ T cells in leukemia is less clear, although in [START_REF] Chen | Development and dynamics of robust T-cell responses to CML under imatinib treatment[END_REF] the authors ascertained that some CML patients under imatinib-induced remission develop an anti-leukemia immune response involving both CD4+ and CD8+ T cells. Hence, our goal in this paper is to capture in a mathematical model the underlying dynamics of this disease by considering the evolution of healthy and leukemic cell populations along with one of the most important component of the cellular immune response to CML, namely T cell response.

Even if a variety of mathematical papers have applied a range of modeling approaches to study tumor-immune interactions in general (see, for example the recent review [START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF]), only a few described the specific leukemia-immune interaction. Leukemia-immune models have been formulated using mostly ordinary differential equations (ODE) ( [START_REF] Moore | A mathematical model for chronic myelogenous leukemia (CML) and T cell interaction[END_REF], [START_REF] Nanda | B cell chronic lymphocytic leukemia -a model with immune response[END_REF]) or delay differential equations (DDE) ( [START_REF] Berezansky | A mathematical model with time-varying delays in the combined treatment of chronic myeloid leukemia[END_REF], [START_REF] Deconde | Post-transplantation dynamics of the immune response to chronic myelogenous leukemia[END_REF], [START_REF] Kim | Dynamics and potential impact of the immune response to chronic myelogenous leukemia[END_REF], [START_REF] Peet | New computational tools for modeling chronic myelogenous leukemia[END_REF]). Some models that specifically study the immune response to CML are [START_REF] Kim | Dynamics and potential impact of the immune response to chronic myelogenous leukemia[END_REF], [START_REF] Berezansky | A mathematical model with time-varying delays in the combined treatment of chronic myeloid leukemia[END_REF], [START_REF] Neiman | A mathematical model of chronic myelogenous leukaemia[END_REF] and [START_REF] Moore | A mathematical model for chronic myelogenous leukemia (CML) and T cell interaction[END_REF]. However, none of the above papers have considered competition between healthy and leukemic cell populations, which is an important factor in CML dynamics.

Assumptions on the model

In this paper, we use a five-dimensional system of DDEs. The first four equations describing the healthy and leukemic cell populations are based on the Mackey and collaborators models of hematopoiesis ( [START_REF] Mackey | Periodic oscillations of blood cell population in chronic myelogenous leukemia[END_REF]). In the present model, we consider two types of hematopoietic cell populations: stem-like cell populations consisting of stem cells and progenitors with self-renew ability and mature cell populations formed by differentiated cells without self-renew ability. We include two types of cell populations, healthy and leukemic, each of them with two subpopulations of cells: stem-like and mature. In this way, we introduce space competition between the normal cell population and the CML one. We underline that in the case of leukemic cell populations mature cells are mostly unable to perform their functions.

The main difference from the Mackey model is considering the competition between healthy and CML cell populations and the fact that three types of division of a stem-like cell are considered: asymmetric division, symmetric self renewal and differentiation. In this paper we use the notation α = h, l with h for healthy and l for leukemia. Consequently, we assume that a fraction η 1α , α = h, l, of stem-like cell population is susceptible to asymmetric division: one daughter cell proceeds to differentiate and the other re-enters the stem cell compartment. A fraction η 2α , α = h, l, is susceptible to differentiate symmetrically with both cells that result following a phase of maturation and the fraction 1 -η 1α -η 2α , α = h, l, is susceptible to self-renewal so both cells that results after mitosis are stem-like cells. This four-dimensional competition model was introduced by Radulescu et al. in [START_REF] Radulescu | A study on stability and medical implications for a complex delay model for CML with cell competition and treatment[END_REF].

The fifth equation of the system models the anti-leukemia T cell immune response. The T cell population considered in this paper consists only of activated anti-leukemia T cells (involving both CD4+ and CD8+ T cells), which actively interact with CML cells. We do not consider other parts of the immune response, as the population of antigen presenting cells or the levels of citokines production. We assume that after encountering a leukemic cell, a T cell has two possibilities: either it inhibits the leukemic cell and activates a feedback function to stimulate the production of new T cells, or it is inhibited itself by the leukemic cell.

Description of the model

The state variables of the model are the healthy cell populations x 1 -stemlike and x 2 -mature, the CML cell populations x 3 -stem-like and x 4 -more differentiated and x 5 -the population of anti-leukemia T cells. The delays for healthy and leukemia stem-like cells are τ 1 and τ 3 for the duration of the cell cycle, independent of the type of division, and τ 2 and τ 4 for the time necessary for differentiation into mature leukocytes for healthy and, respectively, leukemia cells. τ is the duration of the cell cycle for T-cells and τ 5 = nτ with n the number of antigen depending divisions. We denote X τi = X(t -τ i ), where X = (x 1 , x 2 , x 3 , x 4 , x 5 ).

The optimal control model is

ẋ1 = f 1 (x 1 , x 2 , x 3 , x 4 , x 1τ1 , x 2τ1 , x 3τ1 , x 4τ1 ) ẋ2 = f 2 (x 2 , x 1τ2 , x 2τ2 , x 4τ2 ) ẋ3 = f 3 (x 1 , x 2 , x 3 , x 4 , x 5 , x 1τ3 , x 2τ3 , x 3τ3 , x 4τ3 , u 1 , k 1l , u 1τ3 ) ẋ4 = f 4 (x 3 , x 4 , x 5 , x 2τ4 , x 3τ4 , x 4τ4 , k 2l , u 1τ4 ) ẋ5 = f 5 (x 4 , x 5 , x 4τ5 , x 5τ5 , u 2 ) (1)
where

f 1 = -γ 1h x 1 -(η 1h + η 2h )k h (x 2 + x 4 )x 1 -(1 -η 1h -η 2h )β h (x 1 + x 3 )x 1 + +2e -γ 1h τ1 (1 -η 1h -η 2h )β h (x 1τ1 + x 3τ1 )x 1τ1 + η 1h e -γ 1h τ1 k h (x 2τ1 + x 4τ1 )x 1τ1 f 2 = -γ 2h x 2 + A h (2η 2h + η 1h )k h (x 2τ2 + x 4τ2 )x 1τ2 f 3 = -(γ 1l + f 1a )x 3 -[(η 1l + η 2l )k l ((x 2 + x 4 ) f u1 ) + (1 -η 1l -η 2l )β l ((x 1 + x 3 ) f u1 )] x 3 + +[2e -(γ 1l + f1a)τ3 (1 -η 1l -η 2l )β l ((x 1τ3 + x 3τ3 ) f u1τ 3 )+ +η 1l e -(γ 1l + f1a)τ3 k l ((x 2τ3 + x 4τ3 ) f u1τ 3 )]x 3τ3 -b 1 x 3 x 5 l 1 (x 3 + x 4 ) f 4 = -(γ 2l + f 2a ) x 4 + A l (2η 2l + η 1l )k l ((x 2τ4 + x 4τ4 ) f uτ 4 )x 3τ4 -b 2 x 4 x 5 l 1 (x 3 + x 4 ) f 5 = a 1 -a 2 x 5 -a 3 f u2 x 5 l 2 (x 4 ) + 2 n1 a 4 x 5τ5 l 2 (x 4τ5 )
subject to minimization of the cost functional min J(u),

where

J(u) = g(x(T )) + L(t, u(t), x 5 (t)) with g(x(T )) = A 1 x 3 (T ) + A 2 x 4 (T ) + E 1 x 3 (T )/x 1 (T ) + E 2 x 4 (T )/x 2 (T )-being
the weighted sum of the final tumor population and the ratio the ratio of leukemia cells and the healthy ones and

L(t, u(t), x 5 (t)) = T 0 [B 1 u 1 (t) + B 2 u 2 (t) + C 1 k 1 (t) + C 2 k 2 (t) -Dx 5 (t)] dt -
the cumulative drug toxicity and T cell amount.

The history of the state variables is given by

X(θ) = ϕ(θ), θ ∈ [-τ max , 0], τ max = max(τ 1 , τ 2 , τ 3 , τ 4 , τ 5 ).
The healthy and leukemic blood cell populations are seen in competition for resources and this is reflected in the fact that both feedback laws for selfrenewal and differentiation depend on the sum of healthy and leukemia cells.

Consequently, the rate of self-renewal is

β α (x 1 + x 3 ) = β 0α θ mα 1α θ mα 1α + (x 1 + x 3
) mα , with β 0α the maximal rate of self-renewal and θ 1α half of the maximal value and the rate of differentiation is

k α (x 2 + x 4 ) = k 0α θ nα 2α θ nα 2α + (x 2 + x 4
) nα ,with k 0α the maximal rate of differentiation and θ 2α is half of the maximal value. The rest of the parameters for healthy and CML cell populations are: γ 1α -the natural apoptosis, A α -an amplification factor and m α , n α parameters that control the sensitivity of β α respectively k α to changes in the size of stem-like and respectively mature populations. Table 1 contains a complete description of parameters of the model.

To model the influence of T cells on CML cells, we consider the feedback function l 1 (y) = 1 b 5 + y . Consequently, the last terms of the third and fourth equations represent the inhibition of CML cells by anti-leukemia T cells. We assumed that the inhibition of CML cell population by T cells increases with the number of leukemic cells up to a certain level and then reaches a maximal value of inhibition. A further increase in CML population will not modify this value.

As concerns the fifth equation, the first term a 1 is the natural supply of naive T cells, while the second term -a 2 x 5 indicates that T cells exit the population through death at the rate a 2 . Leukemia cells suppress anti-leukemia immune response. The precise mechanism is unknown, but it is assumed that the level of down regulation depends on the current leukemia population, so we consider that the immune system is regulated by the feedback function l 2 (y) = y b 5 + y 2 . This function ensures that T cells are stimulated by CML cells only if leukemia cell population has values in a certain range, called "the optimal load zone" (see [START_REF] Kim | Dynamics and potential impact of the immune response to chronic myelogenous leukemia[END_REF]). We take the rate of antigen stimulation as a feedback function depending on the level of the mature leukemic population, l 2 (x 4 ) and the third and the fourth terms, -a 3 u 2 x 5 l 2 (x 4 ) and 2 n1 a 4 x 5τ6 l 2 (x 4τ5 ) gives the rate at which naive T cells leave and re-enter the effector state after finishing the minimal developmental program of n 1 cell divisions (due to antigen stimulation). The time delay τ 5 = n 1 τ is the duration of this program. These terms represent the loss and respectively the production of T cells due to the competition with leukemic cells.

The treatment targeting the BCR-ABL gene is supposed to affect the apoptosis and the proliferation rates of leukemia cells ( [START_REF] Gottschalk | Imatinib (STI571)-Mediated Changes in Glucose Metabolism in Human Leukemia BCR-ABL-Positive Cells[END_REF]). In view of this fact, we consider the treatment functions

f u1 = 1 1 -u 1 , f 1a = (γ 1h -γ 1l ) k 1l and f 2a = cγ 2h k 2l , with u 1 , k 1l , k 2l : [0, T ] → [0, 1],
where u 1 (t), k 1l (t), k 2l (t) are the treatment effects.

The action of treatment on the proliferation rate will be considered through f u1 in the function of self-renew β l and in the function of differentiation or asymmetric division k l . Note that, in this way, both β l and k l became decreasing functions of u 1 . For more details, see [START_REF] Radulescu | A control delay differential equations model of evolution of normal and leukemic cell populations under treatment[END_REF]. The treatment acts on the apoptosis through the function f 1a on the stem cells and through the function f 2a on the mature ones. Also, from the law of the mass, we have f1a = t t-τ1 k 1l (s)ds.

Moreover, it seems that, in vivo, Imatinib is the trigger of complex mechanisms, some of them able to promote T cell expansion (see [START_REF] Riether | Regulation of hematopoietic and leukemic stem cells by the immune system[END_REF]). Imatinib's effect on T cell population is introduced in the form of a treatment function f u2 = 1 -u 2 , with the stimulatory effect u 2 : [0, T ] → [0, 1]. If no drug is given, then f u2 = 1 and a maximal effect takes place for u 2 (t) = 1, when T cell population is no longer inhibited by CML cell population.

The existence of an optimal control follows since one can transform the given problem into an optimal control problem for a system of ODEs whose solutions will be bounded together with their derivatives on compact intervals (see [START_REF] Benosman | Control of the dynamics of chronic myeloid leukemia by imatinib[END_REF])

Discretization of the Optimal Problem

In this section, we apply the numerical procedure from Gollmann et al. [START_REF] Gollmann | Theory and applications of optimal control problems with multiple time-delays[END_REF], in order to solve the delay optimal control problem (1)+(2) (see also [START_REF] Gollmann | Optimal control problems with delays in state and control variables subject to mixed control state constraints[END_REF], [START_REF] Gollmann | Optimal control problems with control and state delays and applications to growth processes[END_REF]). For that matter, we write the cost functional in the Mayer form

J(u, x) = h(x(T )), x = (x 1 , x 2 , x 3 , x 4 , x 5 ) ∈ R 5 .
We introduce the additional state variable z through the equation

ż(t) = B 1 u 1 (t) + B 2 u 2 (t) + C 1 k 1 (t) + C 2 k 2 (t) -Dx 5 (t), z(0) = 0.
Then, the cost functional (2) is rewritten as

J(u, x, z) = g(x(T )) + z(T ).
In the following, let τ > 0 be such that τ 1 = j 1 τ, τ 2 = j 2 τ, τ 3 = j 3 τ, τ 4 = j 4 τ, τ 5 = j 5 τ, j i ∈ N * , i = 1, 5, T = N τ and use the Euler integration method with a uniform step size τ > 0. Of course, τ can be refined in order to obtain an appropriate smaller step-size. Using the grid points t i = iτ , i = 0, N and the approximations x 1 (t i )

x 1i ∈ R, x 2 (t i ) x 2i ∈ R, x 3 (t i ) x 3i ∈ R, x 4 (t i ) x 4i ∈ R, x 5 (t i ) x 5i ∈ R, u 1 (t i ) u 1i , u 2 (t i ) u 2i and k 1 (t i ) k 1i , k 2 (t i ) k 2i the treatment function f 1a becomes j3 j=1
k 1li-j τ and the delay control problem (1)+( 2) is transformed into the nonlinear programming problem (NLP)

M inimize J = g(x N ) + z N (3) subject to                        x 1i -x 1i+1 + τ f 1 (x 1i , x 2i , x 3i , x 4i , x 1i-j1 , x 2i-j1 , x 3i-j1 , x 4i-k1 ) = 0 x 2i -x 2i+1 + τ f 2 (x 2i , x 1i-j2 , x 2i-j2 , x 4i-j2 ) = 0 x 3i -x 3i+1 + τ f 3 (x 1i , x 2i , x 3i , x 4i , x 5i , x 1i-j3 , x 2i-j3 , x 3i-j3 , x 4i-j3 , u 1i , k 1i , j3 j=1 k 1i-j τ, u 1i-j3 ) = 0 x 4i -x 4i+1 + τ f 4 (x 4i , x 5i , x 2i-j4 , x 3i-j4 , x 4i-j4 , k 2li , u 1i-j4 ) = 0 x 5i -x 5i+1 + τ f 5 (x 4i , x 5i , x 4i-j5 , x 5i-j5 , u 2i ) = 0 z i -z i+1 + τ (B 1 u 1i + B 2 u 2i + C 1 k 1i + C 2 k 2i -Dx 5i ) = 0 (4) -u 1i ≤ 0, u 1i -1 ≤ 0, -u 2i ≤ 0, u 2i -1 ≤ 0, (5)
-k 1i ≤ 0, k 1i -1 ≤ 0, -k 2i ≤ 0, k 2i -1 ≤ 0, i = 0, N -1..
Herein, the initial value profiles ϕ 1 ,ϕ 2 , ϕ 3 , ϕ 4 and ϕ 5 give the values x 1-i := ϕ 1 (-iτ ), i = 0, l 1 ,x 2-i := ϕ 2 (-iτ ), i = 0, l 2 , x 3-j := ϕ 3 (-iτ ), i = 0, l 3

x 4-j := ϕ 4 (-iτ ), i = 0, l 4 , x 5-j := ϕ 5 (-iτ ), i = 0, l 5 .

The variable to be optimized is represented by the vector w = (u 10 , u 20 , k 10 , k 20 , x 11 , .., , x 51 , z 1 , .., u

1 N -1 , u 2 N -1 , k 1 N -1 , k 2 N -1 , x 1 N , .., x 5 N , z N ) ∈ R 10N .

Numerical results

In the following figures, we plotted the trajectories of the healthy, respectively CML cell populations for the competition system, showing a comparison between the dynamics of a system without treatment and the dynamics of a system subject to optimal treatment. In the following simulations two aspects are combined, resulting four distinct manifestations of the disease:

-starting treatment in two different stages of the disease: a less severe stage when the population of leukemic cells and the healthy cells coexist (the leukemia cell population is still small) (S1) and a stage where healthy cell population disappeared and the number of leukemic cells is already very high (S2);

-two configurations of parameters describing two different forms of the disease for patients, configuration 1 and configuration 2 (see table 1).The configuration 2 corresponds to a more serious disease.

Considering multiple effects of treatment, on the apoptosis of leukemic stem cells, leukemic mature cells, proliferation rate and immune system, simulations show the impact of the disease on various optimal control solutions for four hypothetical patients (see figures 1, 2,3 and 4 for comparison between the dynamics of a system without treatment and with with optimal treatment). all situations (see figures). To solve the problem of optimal control the Matlab solver for NLP problems fmincon was used, selecting the ' interior-point' solver.

Conclusions

In this paper, an optimal control model for CML with the influence of the immune system and treatment was investigated. Based on clinical evidences and assumptions, the effects of Imatinib, the current first line treatment in CML, was considered. These effects include the decrease of leukemic proliferation and differentiation, the increase of leukemic apoptosis and some influences on the anti-leukemia immune response.

From figures 1, 2, 3 and 4 one can see one can see the decline of leukemic cells (i.e. molecular remission) and the increase of the number of healthy cells after approximately three months of treatment. Depending on the level of cell populations at diagnosis (S1 or S2) and on the leukemia severity (i.e. configuration 1 or configuration 2) the evolution of healthy cell population to a normal amount is more or less rapid.

The plots of optimal controls (5, 6,7 and 8) exhibit an optimal control effect different for four hypothetical patients. One can observe that the drug influence is slightly different for various manifestations of the disease. Consequently, for an optimal effect of treatment, the prescribed dose should be adapted considering the parameter's disease of a certain patient and the leukemic burden at diagnosis. Although the identification of most of the values parameters of the disease is a daunting task, there are some which can be computed by means of current methods and they might provide an important indication concerning dose adjustment and therapy management.

Fig. 1 :

 1 Fig.1: Simulations start from S1 for configuration 1 of parameters. Dashed line represents the dynamics of a system without treatment and continuous line represents the dynamics of a system with optimal treatment.

Fig. 2 :Fig. 3 :Fig. 4 : 2 Fig. 5 : 2 Fig. 6 : 2 Fig. 7 : 2 Fig. 8 :

 23425262728 Fig. 2: Simulations start from S2 for configuration 1 of parameters. Dashed line represents the dynamics of a system without treatment and continuous line represents the dynamics of a system with optimal treatment.

Table 1 :

 1 Description of parameters Hill coefficient of the β h feedback function 4 4 m l Hill coefficient of the β l feedback function 4 4 n h Hill coefficient of the k h feedback function 3 3 n l Hill coefficient of the k l feedback function 3 3 θ 1h Parameter for the β h feedback function(10 6 cells/kg) 1.6 1.6 θ 2h Parameter for the k h feedback function(10 6 cells/kg) 12 12 θ 1l Parameter for the β l feedback function(10 6 cells/kg) 0.5 0.5 θ 2l Parameter for the k l feedback function(10 6 cells/kg) 36 36 γ 1h Loss of stem cells due to mortality for healthy cells (days -1 ) 0.1 0.1 γ 1l Loss of stem cells due to mortality for leukemic cells (days -1 ) 0.04 0.01 η 1h Rate of asymmetric division for healthy cells

	Par. Description

Acknowledgments. This work was supported by CNCS-ROMANIA Grant ID-PCE-2011-3-0198.


Azmy S Ackleh
email: ackleh@louisiana.edu

Xinyu Li

Baoling Ma
email: baoling.ma@millersville.edu

Parameter Estimation in a Size-Structured Population Model with Distributed States-at-Birth

A least-squares method is developed for estimating parameters in a size-structured population model with distributed states-atbirth from field data. First and second order finite difference schemes for approximating the nonlinear-nonlocal partial differential equation model are utilized in the least-squares problem. Convergence results for the computed parameters are established. Numerical results demonstrating the efficiency of the technique are provided.

Introduction

It is often the case that direct observations of vital rates of individual organisms are not accessible and our knowledge of the vital rates is incomplete. Therefore, the inverse problem approach often plays an important role in deducing such information at the individual level from observation at the population level. In recent years substantial attention has been given to inverse problems governed by age/stage/size structured population models [START_REF] Ackleh | Parameter estimation in the nonlinear size-structured population model[END_REF][START_REF] Ackleh | Parameter identification in size-structured population models with nonlinear indivudual rates[END_REF][START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF][START_REF] Adams | Population dynamics models in plant-insect herbivore-pesticide interactions[END_REF][START_REF] Banks | Modeling and estimation in size structured population models[END_REF][START_REF] Banks | Estimation of growth and survival in size-structured cohort data[END_REF][START_REF] Banks | Estimation of growth rate distributions in size structured population models[END_REF][START_REF] Davidian | Nonlinear models for repeated measurement data: An overview and update[END_REF][START_REF] Fitzpatrick | Parameter estimation in conservation laws[END_REF]. Methodologies applied to solve such inverse problems include the least-squares approach [START_REF] Ackleh | Parameter estimation in the nonlinear size-structured population model[END_REF][START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF][START_REF] Banks | Estimation of growth rate distributions in size structured population models[END_REF] and the fixed point iterative technique [START_REF] Rundell | Determing the death rate for an age-structured population from census data[END_REF]. The least-squares approach has been often used in inverse problems governed by size-structured models. For example, in [START_REF] Banks | Estimation of growth and survival in size-structured cohort data[END_REF][START_REF] Banks | Estimation of growth rate distributions in size structured population models[END_REF] the authors used least-squares method to estimate the growth rate distribution in a linear size-structured population model. A similar technique was applied to a semi-linear size-structured model in [START_REF] Huyer | A size stuctured population model with dispersion[END_REF] where the mortality rate depends on the total population due to competition between individuals. Furthermore, such least-squares methodology has been applied for estimating parameters in general conservation laws [START_REF] Fitzpatrick | Parameter estimation in conservation laws[END_REF]. Therein the author utilizes monotone finite-difference schemes to numerically solve the conservation law and present numerical results for estimation the flux function from numerically generated data. And in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF], the authors solved an inverse problem governed by structured juvenile-adult model. Therein, the least-squares approach was used to estimate growth, mortality and reproduction rates in the adult stage from field data on an urban green tree frog population. The estimated parameters were then utilized to understand the long-term dynamics of this green tree frog population.

In this paper we consider the following nonlinear Gurtin-MacCamy type model with distributed states-at-birth: ∂ ∂t p(x, t; θ) + ∂ ∂x (g(x, t, Q(t; θ))p(x, t; θ)) = -µ(x, t, Q(t; θ))p(x, t; θ)

+ xmax xmin
β(x, y, t, Q(t; θ))p(y, t; θ)dy, (x, t) ∈ (x min , x max ) × (0, T ),

g(x min , t, Q(t; θ))p(x min , t; θ) = 0, t ∈ [0, T ], p(x, 0; θ) = p 0 (x), x ∈ [x min , x max ].
(1) Here, θ = (g, µ, β) is the vector of parameters to be identified. The function p(x, t; θ) is the parameter-dependent density of individuals of size x at time t.

Therefore, Q(t; θ) = xmax xmin p(x, t; θ)dx provides the total population at time t which depends on the vector of parameters θ = (g, µ, β). The functions g and µ represent the individual growth and mortality rates, respectively. It is assumed that individuals may be recruited into the population at different sizes with β(x, y, t, Q) representing the rate at which an individual of size y gives birth to an individual of size x. Henceforth, we will call the model (1) Distributed Size Structured Model and abbreviate it as DSSM.

The main goal of this paper is to develop a least-squares approach for estimating the parameter θ from population data and to provide convergence results for the parameter estimates. The paper is organized as follows. In Section 2, we set up a least-squares problem and present finite difference schemes for computing an approximate solution to this least-squares problem. In Section 3 we provide convergence results for the computed parameters. In Section 4, numerical examples showing the performance of the least-squares technique and an application to a set of field data on green tree frogs are presented.

The Least-Squares Problem and Approximation Schemes

Let

D 1 = [x min , x max ] × [0, T ] × [0, ∞) and D 2 = [x min , x max ] × [x min , x max ] × [0, T ] × [0, ∞) throughout the discussion. Let B = C 1 b (D 1 ) × C b (D 1 ) × C b (D 2 )
, where C b (Ω) denotes the Banach space of bounded continuous functions on Ω endowed with the usual supremum norm and C 1 b (Ω) is the Banach space of bounded continuous functions with bounded continuous derivatives on Ω and endowed with the usual supremum norm. Clearly, B is a Banach space when endowed with the natural product topology. Let c be a sufficiently large positive constant and assume that the admissible parameter space Θ is any compact subset of B (endowed with the topology of B) such that every θ = (g, µ, β) ∈ Θ satisfies (H1)-(H4) below.

(H1) g ∈ C 1 b (D 1 ) with g x (x, t, Q) and g Q (x, t, Q) being Lipschitz continuous in x with Lipschitz constant c, uniformly in t and Q. Moreover, 0 < g(x, t, Q) ≤ c for x ∈ [x min , x max ) and g(x max , t, Q

) = 0. (H2) µ ∈ C b (D 1 ) is Lipschitz continuous in x and Q with Lipschitz constant c, uniformly in t. Also, 0 ≤ µ(x, t, Q) ≤ c. (H3) β ∈ C b (D 2 ) is Lipschitz continuous in Q with Lipschitz constant c, uniformly
in x, y and t. Also, 0 ≤ β(x, y, t, Q) ≤ c and for every partition {x i } N i=1 of [x min , x max ], we have sup

(y,t,Q)∈[xmin,xmax]×[0,T ]×[0,∞) N i=1 |β(x i , y, t, Q) -β(x i-1 , y, t, Q)| ≤ c. (H4) p 0 ∈ BV ([x min , x max ]
), the space of functions with bounded total variation on [x min , x max ], and p 0 (x) ≥ 0.

We now define a weak solution to the model (1).

p 0 (x)φ(x, 0)dx = t 0 xmax xmin p(x, τ ; θ)[φ τ (x, τ ) + g(x, τ, Q(τ ; θ))φ x (x, τ ) -µ(x, τ, Q(τ ; θ))φ(x, τ )]dxdτ + t 0 xmax xmin xmax xmin β(x, y, τ, Q(τ ; θ))p(y, τ ; θ)φ(x, τ )dydxdτ.
(2) for every test function φ ∈ C 1 ((x min , x max ) × (0, T )) and t ∈ [0, T ].

We are interested in the following least-squares problem: given data X s which corresponds to the number of individuals at time t s , s = 1, 2, • • • , S, find a parameter θ = (g, µ, β) ∈ Θ such that the weighted least-squares cost functional

F (θ) = S s=1 |W (Q(t s ; θ)) -W (X s )|
2 is minimized over the admissible parameter space Θ, i.e., find θ * such that

θ * = arg min θ∈Θ F (θ) = arg min θ∈Θ S s=1 |W (Q(t s ; θ)) -W (X s )| 2 , (3)
where

W ∈ C([0, ∞)) is a weight function.
In order to numerically approximate the solution to the minimization problem (3), we first need to approximate the solution of model [START_REF] Ackleh | Parameter estimation in the nonlinear size-structured population model[END_REF]. To this end, we utilize similar finite-difference approximation schemes as those developed in [START_REF] Ackleh | A second order finite difference scheme for a size-structured population model with distributed states-at-birth[END_REF]. Suppose that the intervals [x min , x max ] and [0, T ] are divided into N and L subintervals, respectively. The following notations will be used throughout the pater: ∆x = (x max -x min )/N and ∆t = T /L. The discrete mesh points are given by

x i = x min + i∆x, t k = k∆t for i = 0, 1, • • • , N , k = 0, 1, • • • , L.
For ease of notation, we take a uniform mesh with constant sizes ∆x and ∆t. More general nonuniform meshes can be similarly considered. We shall denote by p k i (θ) and Q k (θ) the finite difference approximation of p(x i , t k ; θ) and Q(t k ; θ), respectively. For convenience we will also use the notation p k i and Q k without explicitly stating their dependence on θ. We also let

g k i = g(x i , t k , Q k ), µ k i = µ(x i , t k , Q k ), β k i,j = β(x i , y j , t k , Q k ). Here, Q k = N i=1 p k i ∆x.
We define the 1 , ∞ norms and TV (total variation) seminorm of the grid functions p k by

p k 1 = N i=1 |p k i |∆x, p k ∞ = max 0≤i≤N |p k i |, T V (p k ) = N -1 i=0 |p k i+1 -p k i |,
and the finite difference operators by

∆ + p k i = p k i+1 -p k i , 0 i N -1, ∆ -p k i = p k i -p k i-1 , 1 i N.
Throughout the discussion, we impose the following CFL condition concerning ∆x and ∆t:

(H5) ∆t ∆x + ∆t 1 c .
We discretize model (1) using the following first order explicit upwind finite difference scheme (FOEU):

p k+1 i -p k i ∆t + g k i p k i -g k i-1 p k i-1 ∆x = -µ k i p k i + N j=1 β k i,j p k j ∆x, 1 ≤ i ≤ N, 0 ≤ k ≤ M -1, g k 0 p k 0 = 0, 0 ≤ k ≤ M, p 0 i = p 0 (x i ), 0 ≤ i ≤ N. (4)
We could write the first equation in (4) equivalently as

p k+1 i = ∆t ∆x g k i-1 p k i-1 + 1 - ∆t ∆x g k i -µ k i ∆t p k i +   N j=1 β k i,j p k j ∆x   ∆t, 1 ≤ i ≤ N, 0 ≤ k ≤ M -1.
(5) It is easy to check that under assumptions (H1)-(H4) FOEU scheme converges to a unique weak solution of system (1) as proved in [START_REF] Ackleh | A second order finite difference scheme for a size-structured population model with distributed states-at-birth[END_REF]. The above approximation can be extended to a family of functions {p ∆x,∆t (x, t; θ)} defined by p ∆x,∆t (x, t; θ

) = p k i (θ) for (x, t) ∈ [x i-1 , x i ) × [t k-1 , t k ), i = 1, 2, • • • , N, k = 1, 2, • • • , M .
Since the parameter set is infinite dimensional, a finite-dimensional approximation of the parameter space is necessary for computing minimizers. Thus, we consider the following finite-dimensional approximations of (3): Let Q ∆x,∆t (t; θ) =

|W (Q ∆x,∆t (t s ; θ)) -W (X s )| 2 . (6)
Here, Θ m ⊆ Θ is a sequence of compact finite-dimensional subsets that approximate the parameter space Θ, i.e., for each θ ∈ Θ, there exist a sequence of θ m ∈ Θ m such that θ m → θ in the topology of B as m → ∞.

Remark 22 If the compact parameter space Θ is chosen to be finite dimensional, then the approximation space sequence can be taken to be Θ m = Θ.

Since the FOEU (4) scheme is first order it would require a large number of grid points to achieve high accuracy. Thus, we next propose a second order minmod finite difference scheme based on MUSCL schemes [START_REF] Ackleh | A second order finite difference scheme for a size-structured population model with distributed states-at-birth[END_REF][START_REF] Leveque | Numerical methods for conservation laws, birkhauser verlag, basel[END_REF][START_REF] Shen | High resolution schemes for a hierarchical sizestructured model[END_REF] to approximate the solutions of the DSSM model ( 1) in the least-squares problem. We begin by using the following second order approximations for the integrals:

Q k = N i=0 p k i ∆x = 1 2 p k 0 ∆x + N -1 i=1 p k i ∆x + 1 2 p k N ∆x and N j=0 β k i,j p k j ∆x = 1 2 β k i,0 p k 0 ∆x + N -1 j=1 β k i,j p k j ∆x + 1 2 β k i,N p k N ∆x.
Then we approximate the model (1) by

p k+1 i -p k i ∆t + f k i+ 1 2 -f k i-1 2 ∆x = -µ k i p k i + N j=0 β k i,j p k j ∆x, i = 1, 2, • • • , N, k = 0, 1, • • • , L -1, g k 0 p k 0 = 0, k = 0, 1, • • • , L, p 0 i = p 0 (x i ), i = 0, 1, • • • , N. (7)
Here, the limiter is defined as

f k i+ 1 2 = g k i p k i + 1 2 (g k i+1 -g k i )p k i + 1 2 g k i mm(∆ + p k i , ∆ -p k i ), i = 2, • • • , N -2, g k i p k i , i = 0, 1, N -1, N, (8)
where mm(a, b) = sign(a) + sign(b) 2 min(|a|, |b|).

Convergence Theory for the Parameter Estimation Problem Using FOEU

The results in this section pertain to the case when (1) is approximated by the FOEU scheme (4). Our future efforts will focus on extending these theoretical results to case when the model ( 1) is approximated by the higher order SOEM scheme [START_REF] Banks | Modeling and estimation in size structured population models[END_REF]. We establish the convergence results for the parameter estimation problem using an approach based on the abstract theory in [START_REF] Banks | Estimation of growth and survival in size-structured cohort data[END_REF]. To this end, we have the following theorem:

Theorem 31 Let θ r = (g r , µ r , β r ) ∈ Θ. Suppose that θ r → θ in Θ and ∆x r , ∆t r → 0 as r → ∞. Let p ∆xr,∆tr (x, t; θ r ) denote the solution of the finite difference scheme with parameter θ r and initial condition p 0 , and let p(x, t; θ) be the unique weak solution of the problem with initial condition p 0 (x) and parameter θ. Then

p ∆xr,∆tr (•, t; θ r ) → p(•, t; θ) in L 1 (x min , x max ), uniformly in t ∈ [0, T ]. Proof. Define p k,r i = p k i (θ r ).
From the fact that Θ is compact and using similar arguments as in [START_REF] Ackleh | A second order finite difference scheme for a size-structured population model with distributed states-at-birth[END_REF], there exist positive constants c 1 , c 2 , c 3 and c 4 such that

p k,r 1 ≤ c 1 , p k,r ∞ ≤ c 2 , T V (p k,r ) ≤ c 3 and N i=1 p m,r i -p n,r i ∆t r ∆x r ≤ c 4 (m -n),
where m > n. Thus, there exist p ∈ BV ([x min , x max ]) such that p ∆xr,∆tr (•, t; θ r ) → p(•, t) in L 1 (x min , x max ) uniformly in t. Hence, from the uniqueness of bounded variation weak solutions which can be established using similar techniques as in [START_REF] Ackleh | A second order finite difference scheme for a size-structured population model with distributed states-at-birth[END_REF], we just need to show that p(x, t) is the weak solution corresponding to the parameter θ.

In order to prove this, let φ ∈ C 1 ([x min , x max ] × [0, T ]) and denote the value of φ(x i , t k ) by φ k i . Multiplying equation ( 5) by φ k+1 i and rearranging some terms we have

p k+1,r i φ k+1 i -p k,r i φ k i = p k,r i (φ k+1 i -φ k i ) + ∆t ∆x [g k,r i-1 p k,r i-1 (φ k+1 i -φ k+1 i-1 ) -(g k,r i-1 p k,r i-1 φ k+1 i-1 -g k,r i p k,r i φ k+1 i )] -µ k,r i p k,r i φ k+1 i ∆t + N j=1
β k,r i,j p k,r j φ k+1 i ∆x∆t.

(9) Multiplying the above equation by ∆x,

summing over i = 1, 2, • • • , N , k = 0, 1, • • • , M -1,
and applying p k 0 = 0 and g k N = 0 we obtain,

N i=1 p L,r i φ L i -p 0,r i φ 0 i ∆x = M -1 k=0 N i=1 p k,r i φ k+1 i -φ k i ∆t ∆x∆t + M -1 k=0 N -1 i=0 g k,r i-1 p k,r i-1 φ k+1 i -φ k+1 i-1 ∆x ∆x∆t - M -1 k=0 N i=1 µ k,r i p k,r i φ k+1 i ∆x∆t + M -1 k=1 N i=1 N j=1 β k,r i,j p k,r j φ k+1 i ∆x∆t∆x.
(10) Using the fact that θ r → θ as r → ∞ in Θ, passing to the limit in [START_REF] Banks | Estimation Techniques for Distributed Parameter Systems[END_REF] we find that p(x, t) is the weak solution corresponding to the parameter θ.

Since W is a continuous on [0, ∞), as an immediate consequence of Theorem 31, we obtain the following.

Corollary 32 Let p ∆xr,∆tr (x, t; θ r ) denote the numerical solution of the finite difference scheme with parameter θ r → θ in Θ and ∆x r , ∆t r → 0 as r → ∞. Then F ∆xr,∆tr (θ r ) → F (θ), as r → ∞.

In the next theorem, we establish the continuity of the approximate cost functional in the parameter θ ∈ Θ (a compact set), so that the computational problem of finding an approximate minimizer has a solution.

Theorem 33 Let ∆x and ∆t be fixed. For each θ ∈ Θ, let p ∆x,∆t (x, t; θ) denote the solution of the finite difference scheme and θ r → θ as r → ∞ in Θ; then

p ∆x,∆t (•, t; θ r ) → p ∆x,∆t (•, t; θ) as r → ∞ in L 1 (x min , x max ) uniformly in t ∈ [0, T ].
Proof. Fix ∆x and ∆t. Define p k,θ r i and p k,θ i to be the solution of the finite difference scheme with parameter θ r and θ, respectively. Let

v k,θ i = p k,θ r i -p k,θ i . Then v k,θ i satisfy the following v k+1,θ i = ∆t ∆x g k,θ r i-1 p k,θ r i-1 -g k,θ i-1 p k,θ i-1 + (p k,θ r i -p k,θ i ) - ∆t ∆x g k,θ r i p k,θ r i -g k,θ i p k,θ i -∆t(µ k,θ r i p k,θ r i -µ k,θ i p k,θ i ) + N j=1 β k,θ r i,j p k,θ r j -β k,θ i,j p k,θ j ∆x∆t, 1 ≤ i ≤ N, 0 ≤ k ≤ M -1, v k+1,θ 0 = p k+1,θ r 0 -p k+1,θ 0 = 0, 0 ≤ k ≤ M -1. (11) Here, Q k,θ r = N i=1 p k,θ r i ∆x, g k,θ r i = g θ r (x i , t k , Q k,θ r
) and similar notations are used for µ k,θ r i and β k,θ r i,j . Using the first equation in [START_REF] Calsina | A model of physiologically structured population dynamics with a nonlinear individual growth rate[END_REF] and assumption (H5) we obtain

N i=1 |v k+1,θ i |∆x ≤ N i=1   1 -∆tµ k,θ r i +   N j=1 β k,θ r i,j ∆x   ∆t   |v k,θ i |∆x -∆t N i=1 g k,θ r i |v k,θ i | -g k,θ r i-1 |v k,θ i-1 | + ∆t N i=1 | g k,θ r i-1 -g k,θ i-1 p k,θ i-1 -g k,θ r i -g k,θ i p k,θ i | +∆t N i=1 |µ k,θ r i -µ k,θ i |p k,θ i ∆x + ∆t N i=1 N j=1 |β k,θ r i,j -β k,θ i,j |p k,θ j ∆x∆x.
(12) By assumption (H1), and ( 11)

N i=1 g k,θ r i |v k,θ i | -g k,θ r i-1 |v k,θ i-1 | = g k,θ r N |v k,θ N | -g k,θ r 0 |v k,θ 0 | = 0. ( 13
)
By assumptions (H2) and (H3) we have

N i=1   1 -µ k,θ r i ∆t +   N j=1 β k,θ r i,j ∆x   ∆t   |u k,θ i |∆x ≤ (1 + c(x max -x min )∆t) v k,θ 1 .
(14) By assumption (H1),

N i=1 | g k,θ r i-1 -g k,θ i-1 p k,θ i-1 -g k,θ r i -g k,θ i p k,θ i | ≤ sup i |g k,θ r i-1 -g k,θ i-1 | N i=1 |p k,θ i -p k,θ i-1 | + N i=1 g k,θ r i -g k,θ r i-1 -g k,θ i -g k,θ i-1 ∆x p k,θ i ∆x. ( 15
)
g k,θ r i -g k,θ r i-1 -g k,θ i -g k,θ i-1 ∆x ≤ 1 0 g θ r x (τ x i-1 + (1 -τ )x i , t k , Q k,θ r ) -g θ r x (τ x i-1 + (1 -τ )x i , t k , Q k,θ ) + 1 0 g θ r x (τ x i-1 + (1 -τ )x i , t k , Q k,θ ) -g θ x (τ x i-1 + (1 -τ )x i , t k , Q k,θ (16
) Assumption (H1), ( 15) and ( 16) yield

N i=1 | g k,θ r i-1 -g k,θ i-1 p k,θ i-1 -g k,θ r i -g k,θ i p k,θ i | ≤ sup i |g k,θ r i-1 -g k,θ i-1 | N i=1 |p k,θ i -p k,θ i-1 | + N i=1 1 0 g θ r x (τ x i-1 + (1 -τ )x i , t k , Q k,θ r ) -g θ r x (τ x i-1 + (1 -τ )x i , t k , Q k,θ ) dτ + 1 0 g θ r x (τ x i-1 + (1 -τ )x i , t k , Q k,θ ) -g θ x (τ x i-1 + (1 -τ )x i , t k , Q k,θ ) dτ p k,θ i ∆x. (17) Note that |Q k,θ r -Q k,θ | = N i=1 (p k,θ r j -p k,θ j )∆x ≤ N i=1 |v k,θ i |∆x = v k,θ 1 . (18)
By the assumption (H1) and the equation above, [START_REF] Rundell | Determing the death rate for an age-structured population from census data[END_REF] yields

N i=1 | g k,θ r i-1 -g k,θ i-1 p k,θ i-1 -g k,θ r i -g k,θ i p k,θ i | ≤ sup i |g k,θ r i-1 -g k,θ i-1 |T V (p k,θ ) + c v k,θ 1 + sup i 1 0 g θ r x ( xi , t k , Q k,θ ) -g θ x ( xi , t k , Q k,θ ) dx p k,θ 1 , (19) where xi = τ x i-1 + (1 -τ )x i . By assumption (H2) N i=1 |µ k,θ r i -µ k,θ i |p k,θ i ∆x ≤ sup i |µ k,θ r i -µ k,θ i | p k,θ 1 . (20)
And from assumption (H3) we obtain

N i=1 N j=1 |β k,θ r i,j -β k,θ i,j |p k,θ j ∆x∆x ≤ sup i,j |β k,θ r i,j -β k,θ i,j | p k,θ 1 . (21)
Using ( 13)-( 21) we arrive at

v k+1,θ 1 ≤ (1 + c(x max -x min )∆t) v k,θ 1 + ∆t sup i |µ k,θ r i -µ k,θ i | p k,θ 1 + sup i,j |β k,θ r i,j -β k,θ i,j p k,θ 1 + sup i |g k,θ r i-1 -g k,θ i-1 |T V (p k,θ ) + c v k,θ 1 + sup i 1 0 g θ r x ( xi , t k , Q k,θ ) -g θ x ( xi , t k , Q k,θ ) dx p k,θ 1 . (22) Note that |µ k,θ r i -µ k,θ i | ≤ |µ θ r (x i , t k , Q k,θ r ) -µ θ r (x i , t k , Q k,θ )| + |µ θ r (x i , t k , Q k,θ ) -µ θ (x i , t k , Q k,θ )|, |g k,θ r i-1 -g k,θ i-1 | ≤ |g θ r (x i-1 , t k , Q k,θ r ) -g θ r (x i-1 , t k , Q k,θ )| + |g θ r (x i-1 , t k , Q k,θ ) -g θ (x i-1 , t k , Q k,θ )|, |β k,θ r i,j -β k,θ i,j | ≤ |β θ r (x i , y j , t k , Q k,θ r ) -β θ r (x i , y j , t k , Q k,θ )| + |β θ r (x i , y j , t k , Q k,θ ) -β θ (x i , y j , t k , Q k,θ )|.
(23) Thus, by assumptions (H1)-(H4) and the equations (23) [START_REF] Shen | High resolution schemes for a hierarchical sizestructured model[END_REF], we have

sup i |µ k,θ r i -µ k,θ i | ≤ c v k,θ 1 + sup i |µ θ r (x i , t k , Q k,θ ) -µ θ (x i , t k , Q k,θ )|, sup i |g k,θ r i-1 -g k,θ i-1 | ≤ c v k,θ 1 + sup i |g θ r (x i-1 , t k , Q k,θ ) -g θ (x i-1 , t k , Q k,θ )|, sup i,j |β k,θ r i,j -β k,θ i,j | ≤ c v k,θ 1 + sup i,j |β θ r (x i , y j , t k , Q k,θ ) -β θ (x i , y j , t k , Q k,θ )|. Set δ k = 3c p k,θ 1 + cT V (p k,θ
) and

ρ r k = p k,θ 1 sup i |µ θ r (x i , t k , Q k,θ ) -µ θ (x i , t k , Q k,θ )| + sup i,j |β θ r (x i , y j , t k , Q k,θ ) -β θ (x i , y j , t k , Q k,θ )| + sup i 1 0 g θ r x ( xi , t k , Q k,θ ) -g θ x ( xi , t k , Q k,θ ) dx + sup i |g θ r (x i-1 , t k , Q k,θ ) -g θ (x i-1 , t k , Q k,θ )|T V (p k
Then we have

v k+1,θ 1 ≤ (1 + c(x max -x min )∆t) v k,θ 1 + ∆t(ρ r k + δ k ) v k,θ 1 . (24)
Since for each k, ρ r k → 0 as r → ∞, the result follows from (24). Next, we establish subsequential convergence of minimizers of the finite dimensional problem [START_REF] Adams | Population dynamics models in plant-insect herbivore-pesticide interactions[END_REF] to a minimizer of the infinite dimensional problem (3).

Theorem 34 Suppose that Θ m is a sequence of compact subsets of Θ. Moreover, assume that for each θ ∈ Θ, there exist a sequence of θ m ∈ Θ m such that θ m → θ as m → ∞. Then the function F ∆x,∆t has a minimizer over Θ m . Furthermore, if θ r m denotes a minimizer of F ∆xr,∆tr over Θ m and ∆x r , ∆t r → 0, then any subsequence of θ r m has a further subsequence which convergence to a minimizer of F . the agreement of best fit model solutions obtained using FOEU and SOEM schemes in solving DSSM with the corresponding data sets with no noise as well as with different noise levels. A comparison of the two finite difference methods performing in the least-squares parameter estimation process was provided in Table 1. It can be seen that using the SOEM scheme for approximating DSSM in the least-squares problem one could obtain a similar estimated value b (where the least-squares errors are at the same scale) with less than 10% of the CPU time compared to that using the FOEU scheme.

Fitting DSSM model to green tree frog population estimates from field data

In this example, we fit DSSM to a set of green tree frog population estimates obtained from capture-mark-recapture (CMR) field data during years 2004-2009 (as shown in Figure 2 (Left) [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF]. The purpose is to estimate individual level vital rates for adult green tree frogs, and then use these parameter estimates to gain understanding of the dynamics of this population. In [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF][START_REF] Ackleh | A nonautonomous juvenile-adult model: well-posedness and long-time behavior via a comparison principle[END_REF] the authors developed a model (referred to as JA model hereafter) to describe the dynamics of green tree frog population by dividing individuals into two stages, juveniles and adults, and set a least-squares problem to obtain the best fitted parameters to the CMR field data. With those estimated parameter values, they obtained the adult frog population curve (see Figure 2 (Left)). Due to the relative short duration for the tadpole stage and the lack of other information regarding tadpoles, the authors in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF] simply assumed constant vital rates for tadpole stage. To circumvent this issue, we consider the short duration of juvenile stage as part of the reproduction process and adopt the DSSM to describe the adult frog dynamics. Here, β(x, y, t, Q) in DSSM represents the rate at which an adult frog of size y gives birth to tadpoles that survive to metamorphose into frogs of size x.

As in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF], we take t = 0 to be the first week in January 2004. We chose ∆t = 1/52. Since there are 52 weeks every year ∆t represents one week. Let X s , s = 1, 2, • • • , 136, denote the observed number of frogs which was estimated statistically from CMR experiment data for 136 weeks during the breeding seasons in this six year experiment. The growth rate and mortality rates are assumed to It shows clearly that DSSM model output agrees with the population estimates resulting from field observations better than the JA model. Specifically, DSSM is more accurate in capturing the population dynamics when adult frog numbers are relatively low. Also, the DSSM fitting bears a smaller least-squares error of 33.5 compared to the JA model fitting which yielded an error of 37.8. Furthermore, the γ distribution that provided the best fit as presented in Figure 2 (Right) indicates that the newly metamorphosed frogs have body length between 1.5 cm and 2 cm and approximately 99.6% of adult frogs give birth to tadpoles that eventually metamorphose into frogs of size between 1.5 cm and 2.0 cm. We also applied a statistically based method to compute the variance in the estimated model parameters θ = (α 1 , • • • , α 8 ) similar to the work in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF] using standard regression formulations [START_REF] Davidian | Nonlinear models for repeated measurement data: An overview and update[END_REF]. Table 2 provides the standard deviation for α i s estimated above.

  xmax xmin p ∆x,∆t (x, t; θ)dx denote the finite difference approximation of the total population and consider the finite dimensional minimization problem arg min θ∈Θm F ∆x,∆t (θ) = arg min θ∈Θm S s=1

Fig. 1 .

 1 Fig. 1. (a) Comparison of data sets without noise and the corresponding best fit model solution using FOEU and SOEM schemes in solving DSSM. (b) (c) (d) Comparison of data at different noise levels (σ = 0.05, 0.10, 0.15) and the corresponding best fit model solutions using FOEU and SOEM schemes in solving DSSM.

Table 1 .

 1 Performance comparisons of FOEU and SOEM schemes in least-squares parameter estimates process

	Without noise	With noise
	σ = 0.05	σ = 0.10	σ = 0.15
	FOEU SOEM FOEU SOEM FOEU SOEM FOEU SOEM
	Estimated value of b 0.2507 0.2503 0.2674 0.2768 0.2853 0.2837 0.3033 0.3028
	Least-squares error 6.33-04 3.85e-04 0.0838 0.1871 0.3112 0.2857 0.6583 0.6473
	CPU time (in minutes) 285.48 20.67 264.59 16.45 367.49 22.47 368.71 15.33

Table 2 .

 2 Parameter estimation values and corresponding standard deviation Estimated value 0.486 2.973 0.0085 2.376 0.000 47.061 6.740 1.849 Standard deviation 0.0183 0.7044 0.1164 1.6373 0.125 2.8295 0.8000 1.4763

	α1	α2	α3	α4	α5	α6	α7	α8

Acknowledgments

The work of A.S. Ackleh, X. Li, and B. Ma are partially supported by the National Science Foundation under grant # DMS-1312963.

Proof. The proof here is a direct application of the abstract theory in [START_REF] Banks | Estimation Techniques for Distributed Parameter Systems[END_REF], base on the convergence of F ∆xr,∆tr (θ r ) → F (θ).

Numerical Results

In this section we present several numerical simulations to demonstrate the performance of the parameter estimation methodology. Although the theory presented here applies for the case of infinite dimensional parameter space Θ, for simplicity we restrict the unknown parameter space to finite-dimensional in the examples below.

Convergence of parameter estimates computed by FOEU and SOEM in the least-squares problem

In this example, we test the performance of the parameter-estimation technique using both FOEU and SOEM approximation schemes. As a first step in generating data, we choose ∆x = 0.0100, ∆t = 0.0025, x min = 0,

and the initial condition

0.8, 0.25 ≤ x ≤ 0.45, 2.5, 0.45 < x ≤ 0.65, 0.7, 0.65 < x < 0.8, 0, else.

We then solve system (1) with this choice of parameters in MATLAB using SOEM discretization and collect the resulting total population

Assume all parameters are known except for µ = b exp(Q) with b being an unknown parameter to be estimated. In our parameter estimation simulations, we fixed ∆x = 0.005, ∆t = 0.0025 for FOEU scheme. As for the SOEM scheme, the mesh sizes were chosen to be four times larger, that is, ∆x = 0.020 and ∆t = 0.010. We began with the above-mentioned data without noise in the least-squares problem described in Section 2. We then modified the data by adding normally distributed noise with mean zero and standard deviation σ = 0.05, 0.10, and 0.15, respectively, to the data. To solve the least-squares minimization problem we set θ = b and use the goal function

i.e., W = 1. For each data set the least-squares minimization process was performed to estimate b using both numerical schemes. Our simulation results corroborate the convergence results of computed parameters. Figure 1 demonstrates take the same forms as in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF]:

Here, the mortality rate was assumed to depend linearly on density as well as time since frogs hibernate during winter time. By monitoring program [START_REF] Pham | Population estimates of hyla cinerea (Schneider) (green treefrog) in an urban environment[END_REF] the breeding season begins around the middle of April and ends in early August. Thus, similar to [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF] the birth rate function was assumed to be

Here, the gamma distribution density function with the shape parameter α 7 and scale parameter α 8 , γ(x, α 7 , α 8 ), was chosen to model the size distribution of newly metamorphosed frogs. The constant is a positive small number that allows β to be extended to (x, y, t) ∈ [1.5, 6] × [1.5, 6] × [0, 1] and to satisfy the smoothness properties in (H4). We then also extend β periodically over one year intervals [t, t + 1], t = 1, 2, • • • , 5.

We have α 1 , • • • , α 8 as unknown constants to be estimated (i.e., θ = (α 1 , • • • , α 8 )). We chose the initial condition in DSSM to be p(x, 0) = 615.96 exp(-0.75x) which implies Q(0) = 257.2 (cf. [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF]). To solve the least-squares minimization problem, similar to [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF], we set the goal function to be

To guarantee that the estimated parameter values are biologically relevant, we set appropriate upper and lower bounds for each α i . That is, 2, we simulated DSSM and compared the resulting adult frog population approximations to the data as well as the population estimates from the JA model in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF]. The comparison results are demonstrated in Figure 2 (Left). 1) and the JA model in [START_REF] Ackleh | Fitting a structured juvenile-adult model for green tree frogs to population estimate from capture-markrecapture field data[END_REF]. Right: The probability density function that an adult frog gives birth to tadpoles that eventually metamorphose into frogs of size x.