Ioannis K Dassios
email: jdasios@math.uoa.gr

Krzysztof J Szajowski
email: krzysztof.szajowski@pwr.edu.pl

A non-autonomous stochastic discrete time system with uniform disturbances

Keywords: Bayes control, optimal, singular system, disturbances, Pareto distribution

The main objective of this article is to present Bayesian optimal control over a class of non-autonomous linear stochastic discrete time systems with disturbances belonging to a family of the one parameter uniform distributions. It is proved that the Bayes control for the Pareto priors is the solution of a linear system of algebraic equations. For the case that this linear system is singular, we apply optimization techniques to gain the Bayesian optimal control. These results are extended to generalized linear stochastic systems of difference equations and provide the Bayesian optimal control for the case where the coefficients of these type of systems are non-square matrices. The paper extends the results of the authors developed for system with disturbances belonging to the exponential family.

Introduction

Linear stochastic discrete time systems (or linear matrix stochastic difference equations), are systems in which the variables take their value at instantaneous time points. The horizon of control depends on the problem. The state at instance n depends on random disturbance and the chosen controls. Discrete time systems differ from continuous time ones in that their signals are in the form of sampled data. With the development of the digital computer, the stochastic discrete time system theory plays an important role in the control theory. In real systems, the discrete time system often appears when it is the result of sampling the continuous-time system or when only discrete data are available for use. The investigation aims are, when such system is under consideration, determining the control goals, performance measures and the information available at moments of controls' specification. The small deviations of the parameters can be treated as disturbances. As the random disturbance is admitted the performance measure will be the mean value of the deviation of the states from the required behavior of the system. When all the parameters of the system are known and the distribution of disturbances is well defined then the optimal control can be determined at least for the finite horizon case. The extension of the model to the adaptive one means that the disturbances are not precisely described. Adaptive control is the control method used by a controller which must adapt to a controlled system with parameters which vary, or are initially uncertain (c.f. Black et al. [START_REF] Black | Adaptive systems: History, techniques, problems, and perspectives[END_REF] or Tesfatsion [START_REF] Tesfatsion | A dual approach to Bayesian inference and adaptive control[END_REF] for the history of the adaptive control). Under some unification the model of adaptive control of the linear system is formulated as a control of the discrete time Markov process (cf. [START_REF] Duncan | Adaptive control of a partially observed discrete time Markov process[END_REF]).

It is assumed that the disturbance has a fixed probabilistic description which is determined by the assumption. In this paper it is assumed that the distribution function is known to be an accuracy of parameters and the disturbances additionally change the state of the system. It resembles the statistical problem of estimation. It was seminal paper by Wald [START_REF] Wald | Contributions to the theory of statistical estimation and testing hypotheses[END_REF] where the background of the modern decision theory was established (cf. [START_REF] Wald | Statistical Decision Functions[END_REF]Chapt. 7]). The decision theory approach to the control problems were immediately applied (see books by Sworder [START_REF] Sworder | Optimal adaptive control systems[END_REF], Aoki [START_REF] Aoki | Optimization of stochastic systems[END_REF], Sage and Melsa [START_REF] Sage | Estimation theory with applications to communications and control[END_REF]). The new class of control systems under uncertainty was called adaptive (cf. [START_REF] Tesfatsion | A dual approach to Bayesian inference and adaptive control[END_REF], [START_REF] Black | Adaptive systems: History, techniques, problems, and perspectives[END_REF]). In these adaptive control problems the important role have Bayesian systems. In this class of control models it is assumed that the preliminary knowledge of the disturbances is given by a priori distributions of their parameters. The aim is to construct the controls in a close form. The construction of the Bayes control is also auxiliary for the construction of minimax controls (see Szajowski and Trybu la [START_REF] Szajowski | Minimax control of a stochastic system with the loss function dependent on parameter of disturbances[END_REF], Porosiński and Szajowski [START_REF] Porosiński | A minimax control of linear systems[END_REF], Grzybowski [START_REF] Grzybowski | Minimax control of a system with actuation errors[END_REF], González-Trejo et al. [START_REF] González-Trejo | Minimax control of discrete-time stochastic systems[END_REF]). It is observed the interested in various models of disturbance structure (cf. Duncan & Pasik-Duncan [START_REF] Duncan | Discrete time linear quadratic control with arbitrary correlated noise[END_REF]) and the disturbance distributions (cf. Walczak [START_REF] Walczak | Bayes and minimax control of discrete time linear dynamical systems[END_REF][START_REF] Walczak | Bayesian control of a discrete-time linear system with uniformly distributed disturbances[END_REF]). Stochastic discrete time systems have many applications which we have described in [START_REF] Dassios | Bayesian optimal control for a non-autonomous stochastic discrete time system[END_REF] where the Bayes control of the linear system with quadratic cost function and the disturbances having the distribution belonging to the exponential family with conjugate priors is solved.

The paper is organized as follows: the description of the stochastic discrete time systems is subject of the section 1.1 and some remarks on disturbances are given in the section 1.2. In the section 2 we determine the Bayes control for the conjugate prior distribution π of the parameter λ as the solution of a singular linear system and provide optimal Bayesian control. We close the paper by studying the Bayes control of a class of generalized linear stochastic discrete time systems.

Stochastic discrete time systems

Let x n ∈ R m be the state of the system, u n ∈ R m be the control. Assume that

υ n ∈ V ⊂ R m , with υ n = (υ 1 n , υ 2 n , . . . , υ k n , 0, . . . , 0) T , is the disturbance at time n and α n , b n , c n ∈ R m×m . Consider a stochastic discrete time system (cf.
Kushner [START_REF] Kushner | Introduction to stochastic control[END_REF])

x n+1 = α n x n + b n u n + c n υ n , ∀n = 0, 1, . . . , N -1.
(

The horizon N of the control, the time up to which the system is controlled, is a random variable, independent of the disturbances υ 0 , υ 1 , . . ., and has the following known distribution

P {N = k} = p k , ∀k = 0, 1, ..., M, M i=0 p k = 1, p M = 0. ( 2
)
In the authors paper [START_REF] Dassios | Bayesian optimal control for a non-autonomous stochastic discrete time system[END_REF] it was considered the family of the exponentially distributed disturbances. Let us assume here that the disturbances υ i n have the uniform distributions on [0, λ i ] with parameter λ i ∈ + , i = 1, 2, . . . , k and

X n = (x 0 , x 1 , ..., x n ), U n = (u 0 , u 1 , ..., u n ), λ = (λ 1 , λ 2 , ..., λ k , 0, ..., 0) T .
For convenience U M will be denoted by U and called a control policy. Definition 1. The control cost for a given policy U (the loss function) is the following

L(U, X N ) = N i=0 (y T i s i y i + u T i k i u i ), (3)
where

k i ∈ R m×m ≥ 0 m,m , are symmetric matrices, s i ∈ R 2m×2m ≥ 0 2m,2m and
y i =   x i • • • λ   ∈ R 2m , ∀i = 0, 1, ..., M .
With 0 i,j we will denote the zero matrix i × j.

Let the prior distribution π of the parameter λ be given. It is considered the Pareto priors (see [4, Ch. 9.7], [START_REF] Philbrick | A practical guide to the single parameter pareto distribution[END_REF]) with parameters r i > 0,

β i > 2 g(λ i |β i , r i ) = β i r βi i λ βi+1 i I [ri,∞) (λ i ). ( 4
)
Denote E N , E λ the expectations with respect to the distributions of N and random vectors υ 0 , υ 1 , ... (when λ is the parameter), E π and E are the expectations with respect to the distribution π and to the joint distribution υ n and λ, respectively.

Definition 2. (see [START_REF] Kushner | Introduction to stochastic control[END_REF], [START_REF] Porosiński | Bayes control for a multidimensional stochastic system[END_REF], [START_REF] Sawitzki | Exact filtering in exponential families: discrete time[END_REF], [START_REF] Walczak | Bayes and minimax control of discrete time linear dynamical systems[END_REF][START_REF] Walczak | Bayesian control of a discrete-time linear system with uniformly distributed disturbances[END_REF]) Let L(•, •) be the loss function given by (3).

(a) The risk connected with the control policy U , when the parameter λ is given, is defined as follows

R(λ, U ) = E N E λ [L(U, X N ) | X 0 ] = E N E λ [ N i=0 y T i s i y i + u T i k i u i | X 0 ] . (b)
The expected risk r, associated with π and the control policy U , is equal to

r(π, U ) = E π [R(λ, U )] = E N E[ N i=0 y T i s i y i + u T i k i u i | X 0 ] .
(c) The expected risk r, associated with π and the control policy U , is equal to

r n (π, U n ) = E N [E π [R(λ, U )] = E N E[ N i=n y T i s i y i + u T i k i u i | X 0 ] .
Let the initial state x 0 and the distribution π of the parameter λ be given.

Definition 3. A control policy U * is called the Bayes policy when r(π, U * ) = inf U ∈℘π r(π, U ), where ℘ π is the class of the control policies U for which exists r(π, U ).

Filtering

Let us assume that the random variables υ n have the density p(υ n , λ) with respect to a σ-finite measure µ on R. The consideration is focused on the special case when each coordinate has the uniform distribution, i.e. the density p(υ n , λ) has the following representation:

p(υ n , λ) = k i=1 p(υ i n , λ i ), (5)
where p(υ i n , λ i ) = 1 λi I [0,λi] (υ i n ), for all i = 1, 2, ..., k. V * i is the set of the random variables υ i n . We have:

E λi [υ i n ] = λ i 2 = q i λ i and, E λi [(υ i n ) 2 ] = λ 2 i 3 = q 1,i λ 2 i ,
where q i , q 1,i are constants. Let λ have the a priori distribution π with density

g(λ | β, r) = k i=1 g i (λ i ; β i , r i ), (6)
where

g i (λ i | β i , r i ) is given by (4) where β ∈ S β k ⊂ R m , r ∈ S r k ⊂ R m with β = (β 1 , β 2 , . . . , β k , 0, . . . , 0) T ,
and r = (r 1 , r 2 , . . . , r k , 0, . . . , 0) T . When such the a priori density is assigned to λ i and then the object of filtering, to determine the Bayes control, is to produce a posteriori density for λ i after any new observations of the state of the system. We change the control after obtaining the new data. Hence, to determine the Bayes control, a posteriori density for λ must be obtained after any new observations.

The Bayes Control

Suppose the initial state x 0 is given, the disturbances have the distribution with the density given by ( 5) and the prior distribution π of the parameter λ is given by [START_REF] Duncan | Discrete time linear quadratic control with arbitrary correlated noise[END_REF]. Let the distribution of the random horizon N be given by [START_REF] Black | Adaptive systems: History, techniques, problems, and perspectives[END_REF]. Consider the problem of the Bayes control for the system (1) with the starting point at the moment n, when X n , U n-1 are given. The expected risk is then given by (c.f Defintion 2 (c); see [START_REF] Porosiński | Bayes control for a multidimensional stochastic system[END_REF], [START_REF] Dassios | Bayesian optimal control for a non-autonomous stochastic discrete time system[END_REF])

r n (π, U n ) = E M i=n (y T i s i y i + u T i k i u i ) | X n , U n-1 . ( 14
)
Let us denote ϕ k = M i=k p i . We have

r n = E[ M i=n ϕ i ϕ n (y T i s i y i + u T i k i u i ) | X n , U n-1 ].
For the above truncated problem we provide the following definitions:

Definition 4. The Bayes risk is defined as

W n = inf U n r n (π, U n ), (15)
where r(π, U n ) is the expected risk defined in the defintion 2 (c) and the formulae (14).

Definition 5. If there exists U n * = (u * n , u * n+1 , .., u * N ) such that W n = r(π, U n * ), then U n * will be called the Bayes policy and u * i , i = n, n + 1, ..., N the Bayes controls for truncated control problem.

Obviously, r(π, U 0 ) = r(π, U ), W 0 = r(π, U * ). For the solution of the Bayes control problem we derive the Bayes controls u * n for n = N, N -1, ..., 1, 0 recursively. Then U 0 * is the solution of the problem. From the Bellman's dynamic programming optimality principle we obtain the following Lemma, see [START_REF] Porosiński | Bayes control for a multidimensional stochastic system[END_REF]. Lemma 2. Assume the stochastic discrete time system (1). Then the Bayes risk W n has the form

W n = x T n A n x n + 2r T n B n x n + 2r T n C n r n , (16)
where

A n , B n , C n ∈ R m×m , D n ∈ R m with A n = f 1 (s n ), B n = f 2 (Q n , Q n 2 , s n ), C n = f 3 (Q n , Q n 1 , Q n 3 , Q n 4 , s n ).
The functions f j , j = 1, 2, 3 are strictly monotonic, differentiable. The constants Q n,i , Q n,i j , j = 1, 2, 3, 4, n = 0, 1, ..., N are given by [START_REF] González-Trejo | Minimax control of discrete-time stochastic systems[END_REF], [START_REF] Grzybowski | Minimax control of a system with actuation errors[END_REF] and s n is defined in (3).

Bayesian optimal control for stochastic discrete time systems

We can now prove the following theorem Theorem 1. Assume the stochastic discrete time system [START_REF] Aoki | Optimization of stochastic systems[END_REF]. Then, the Bayes control u * n is given by the solution of the linear system

K n u * n = L n , (17)
where

K n = k n + b T n A n+1 b n ( 18
)
and

L n = -b T n [A n+1 α n x n + (A n+1 c n Q n + B n+1 Q n 2 )r n ]. ( 19
)
The matrices k n , A n , Q n , are defined in (3), ( 16), the lemma 1, respectively and e = n-1 j=0 Q j r j .

Proof. From ( 15), the Bayes risk is given by W n = inf U n r(π, U n ). It is, equivalently,

W n = min U n E M i=n (y T i s i y i + u T i k i u i ) | X n , U n-1 .
We have

W n = min un u T n k n u n + E y T n s n y n | X n , U n-1 + min U n+1 E E k i=n+1 (y T i s i y i + u T i k i u i ) | X n , U n-1 . It means W n = min un u T n k n u n + E y T n s n y n | X n , U n-1 + E [W n+1 | X n , U n-1 ]
. Hence, the Bayes control u * n satisfies the equation (∇ is the gradient):

∇ un u T n k n u n + E y T n s n y n | X n , U n-1 + E [W n+1 | X n , U n-1 ] un=u * n = 0 m,1 .
By using [START_REF] Sworder | Optimal adaptive control systems[END_REF] we get

k n u * n + b T n A n+1 (α n x n + b n u n + c n E(υ n | X n , U n-1 )) + E b T n B n+1 r n+1 un=u * n | X n , U n-1 = 0 m,1 .
By the properties of conjugate priors for the uniform distribution (see the lemma 1 we have

k n u * n +b T n A n+1 (α n x n +b n u n +c n Q n r n )+E b T n B n+1 r n+1 un=u * n | X n , U n-1 = 0 m,1 ,
and at the end (k

n +b T n A n+1 b n )u * n = -b T n [A n+1 α n x n +(A n+1 c n Q n +B n+1 Q n 2 )r n ]. The proof is completed. (c) L n /
∈ colspanK n and K n is rank deficient, a Bayesian optimal control is given by û

* n = (K T n K n + E T E) -1 K T n L n . ( 25
)
Where E is a matrix such that K

T n K + E T E is invertible and E 2 = θ, 0 < θ << 1. (d) m > r, K n is full rank, a Bayesian optimal control is given by û * n = K T n (K n K T n ) -1 L n . (26)
(e) L n ∈ colspanK n and K n is rank deficient, a Bayesian optimal control is given by (25).

The proof is based on ideas similar to those used in prove [START_REF] Dassios | Bayesian optimal control for a non-autonomous stochastic discrete time system[END_REF]Th. 3] and is omitted here.

Conclusions

In this article we focused on developing the Bayesian optimal control for a class of non-autonomous linear stochastic discrete time systems of type [START_REF] Aoki | Optimization of stochastic systems[END_REF]. Firstly, we proved that the Bayes control of these type of systems is the solution of a linear system of algebraic equations which can also be singular. For this case we used optimization techniques to derive the optimal Bayes control for [START_REF] Aoki | Optimization of stochastic systems[END_REF]. In addition, we used these methods to obtain the Bayesian optimal control of the non-autonomous linear stochastic discrete time system of type [START_REF] Black | Adaptive systems: History, techniques, problems, and perspectives[END_REF], where the coefficients of this system are non-square matrices.

The further extension of this paper is to study to Bayes control problem of stochastic fractional discrete time systems. The fractional nabla operator is a very interesting tool when applied to systems of difference equations and has many applications especially in macroeconomics, since it succeeds to provide information from a specific year in the past until the current year. For all these there is some research in progress.

Acknowledgments

I. Dassios is supported by Science Foundation Ireland (award 09/SRC/E1780).

This is possible if for n = 0, 1, ...N -1 and a given x 0 , we can derive υ n from (1), i.e. the equations

If for a value of n, the matrix c n is singular, we will have to compute the Moore-Penrose Pseudoinverse c † n and then use the following expression

The Moore-Penrose pseudo-inverse can be calculated via the singular value decomposition of c n (see [START_REF] Rugh | Linear system theory[END_REF]). In these cases a posteriori density f (λ | X n , U n-1 ) of the parameter λ, after having observed X n and chosen U n-1 , has the same form as [START_REF] Duncan | Discrete time linear quadratic control with arbitrary correlated noise[END_REF] i.e.

For known X n and U n-1 , the conditional distribution of υ n has the density

where

, for n = 0, 1, ..., M -1, i = 1, 2, ..., k. In addition (see [START_REF] Walczak | Bayes and minimax control of discrete time linear dynamical systems[END_REF][START_REF] Walczak | Bayesian control of a discrete-time linear system with uniformly distributed disturbances[END_REF]) by direct calculation we get Lemma 1. The following equations are fulfilled:

where

Similarly like for the system with the disturbances belonging to the exponential family (see [START_REF] Dassios | Bayesian optimal control for a non-autonomous stochastic discrete time system[END_REF]) we get Theorem 2. Consider the system (1) and the matrices K n , L n as defined in (18), ( 19) respectively. Then (a) ∀n such that K n is full rank, the Bayes control u * n , is given by

Where E is a matrix such that

Bayesian optimal control for generalized stochastic discrete time systems

In this subsection we will expand the results of the section 2.1 by studying Bayesian optimal control for a class of linear stochastic discrete time systems with non-square coefficients. We consider the following non-autonomous linear stochastic discrete time system

Where Definition 6. We will refer to system (22) as a generalized stochastic linear discrete time system.

In the above definition we use the term "generalized " because the coefficients in the system ( 22) can be either square or non-square matrices.

Theorem 3. Consider the system (22) for r = m and assume the matrices K n , L n as defined in (18), ( 19) respectively. Then, ∀n such that (a) m < r, rank(K n ) = m and L n ∈ colspanK n , the Bayes control u * n , is given by u

(b) m < r, rank(K n ) = m and L n / ∈ colspan(K n ), a Bayesian optimal control is given by û


H T Banks

J E Banks

Neha Murad

J A Rosenheim

K Tillman

Modelling Pesticide Treatment Effects on Lygus hesperus in Cotton Fields

Keywords: inverse problems, generalized least squares, model selection, information content, residual plots, piecewise linear splines, Hemiptera, herbivory, pesticide

We continue our efforts on modeling of the population dynamics of herbivorous insects in order to develop and implement effective pest control protocols. In the context of inverse problems, we explore the dynamic effects of pesticide treatments on Lygus hesperus, a common pest of cotton in the western United States. Fitting models to field data, we consider model selection for an appropriate mathematical model and corresponding statistical models, and use techniques to compare models. We address the question of whether data, as it is currently collected, can support time-dependent (as opposed to constant) parameter estimates.

Introduction

When addressing questions in fields ranging from conservation science to agricultural production, ecologists frequently collect time-series data in order to better understand how populations are affected when subjected to abiotic or biotic disturbance [START_REF] Brabec | Population dynamics of aphids on cereals: Digging in the time-series data to reveal population regulation caused by temperature[END_REF][START_REF] Brawn | The role of disturbance in the ecology and conservation of birds[END_REF][START_REF] Pickett | Natural disturbance and patch dynamics: An introduction[END_REF]. Fitting models to data, which generally requires a broad understanding of both statistics and mathematics, is an important component of understanding pattern and process in population studies. In agricultural ecology, pesticide disturbance may disrupt predator-prey interactions [START_REF] Stark | How risky is risk assessment? The role that life history strategies play in susceptibility of species to pesticides and other toxicants[END_REF][START_REF] Stark | Incorporating ecologically relevant measures of pesticide effect for estimating the compatibility of pesticides and biocontrol agents[END_REF] as well as impose both acute and chronic effects on arthropod populations. In the past several decades, the focus of many studies of pesticide effects on pests and their natural enemies has shifted away from static measures such as the LC50, as authors have emphasized population metrics/outcomes [17-19, 26, 29]. Simple mathematical models, parameterized with field data, are often used to then predict the consequences of increasing or decreasing pesticide exposure in the field. Accuracy in parameter estimation and quantification of uncertainty in fitting data to models, which has recently received increasing attention in ecological circles [START_REF] Hilborn | The Ecological Detective: Confronting Models with Data[END_REF][START_REF] Motulsky | Fitting Models to Biological Data Using Linear and Nonlinear Regression: A Practical Guide to Curve Fitting[END_REF], depends critically on the appropriate model selection. In most cases, this includes selection of both statistical and mathematical models in fits-to-data -something that is not always fully explicitly addressed in the ecological literature. We first addressed this gap in [START_REF] Banks | Model comparison tests to determine data information content[END_REF][START_REF] Banks | Modeling populations of Lygus Hesperus on cotton fields in the San Joaquin Valley of California: The importance of statistical and mathematical model choice[END_REF] using data from pest population counts of Lygus hesperus Knight (Hemiptera: Meridae) feeding on pesticide-treated cotton fields in the San Joaquin Valley of California [START_REF] Rosenheim | Estimating the impact of Lygus hesperus on cotton: The insect, plant and human observer as sources of variability[END_REF].

In particular, in [START_REF] Banks | Model comparison tests to determine data information content[END_REF][START_REF] Banks | Modeling populations of Lygus Hesperus on cotton fields in the San Joaquin Valley of California: The importance of statistical and mathematical model choice[END_REF] we investigated the effect of pesticide treatments on the growth dynamics of Lygus hesperus. This was done by constructing mathematical models and then fitting these models to field data so as to estimate growth rate parameters of Lygus hesperus both in the absence and in the presence of pesticide application. Overall, compelling evidence was found for the untreated fields, using model comparison tests, that it may be reasonable to ignore nymph mortality (i.e., just count total number of L. hesperus and not distinguish between nymphs and adults). This would greatly simplify the models, as well as the data collection process.

In the present effort we further examine the importance of model selection and demonstrate how optimal selection of both statistical and mathematical models is crucial for accuracy in parameter estimation and uncertainty quantification in fitting data to models. This report further investigates these issues by testing different data sets from the same database as in [START_REF] Banks | Model comparison tests to determine data information content[END_REF][START_REF] Banks | Modeling populations of Lygus Hesperus on cotton fields in the San Joaquin Valley of California: The importance of statistical and mathematical model choice[END_REF] but with a varied number of pesticide applications in treated fields.

Methods

The data used came from a database consisting of approximately 1500 replicates of L. hesperus density counts, using sweep counts, in over 500 Pima or Acala cotton fields in 1997-2008 in the San Joaquin Valley of California. This data is described more fully in [START_REF] Banks | Modeling populations of Lygus Hesperus on cotton fields in the San Joaquin Valley of California: The importance of statistical and mathematical model choice[END_REF]. We selected subsets to analyze using the following criteria:

-In each replicate (corresponding to data collected during one season at one field) we considered data that was collected by pest control advisors (PCAs) between June 1 and August 30. -We considered data which had the pesticide applications that targeted beet armyworms, aphids, mites as well as Lygus. -We used only replicates where adult and nymph counts were combined into a total insect count. -No counts were made on the days of pesticide applications.

-Superposition of pesticide applications has not been incorporated in the algorithm, so we chose samples with at least a week gap between consecutive pesticide applications.

We consider inverse or parameter estimation problems in the context of a parameterized (with vector parameter q ∈ Ω κq ⊂ R κq ) N -dimensional vector dynamical system or mathematical model given by

dx dt (t) = g(t, x(t), q), (1)
x(t 0 ) = x 0 , (2)
with scalar observation process

f (t; θ) = Cx(t; θ), (3)
where θ = (q T , xT

0 ) T ∈ Ω κ θ ⊂ R κq+ Ñ = R κ θ , Ñ ≤ N ,
and the observation operator C maps R N to R 1 . The sets Ω κq and Ω κ θ are assumed known restraint sets for the parameters.

We make some standard statistical assumptions (see [START_REF] Banks | Modeling and Inverse Problems in the Presence of Uncertainty[END_REF][START_REF] Banks | Mathematical and Experimental Modeling of Physical and Biological Processes[END_REF][START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF][START_REF] Seber | Nonlinear Regression[END_REF]) underlying our inverse problem formulations.

-A1) Assume E i are independent identically distributed i.i.d with E(E i ) = 0 and cov(E i , E i ) = σ 2 0 , where i = 1, ..., n and n is the number of observations or data points in the given data set taken from a time interval [0, T ].

-A2) Assume that there exists a true or nominal set of parameters

θ 0 ∈ Ω ≡ Ω κ θ . -A3) Ω is a compact subset of Euclidian space of R κ θ and f (t, θ) is continuous on [0, T ] × Ω.
Denote as θ the estimated parameter for θ 0 ∈ Ω. The inverse problem is based on statistical assumptions on the observation error in the data. If one assumes some type of generalized relative error data model, then the error is proportional in some sense to the measured observation. This can be represented by a statistical model with observations of the form

Y i = f (t i ; θ 0 ) + f (t i ; θ 0 ) γ E i , γ ∈ [0, 2], (4)
with corresponding realizations

y i = f (t i ; θ 0 ) + f (t i ; θ 0 ) γ i , γ ∈ [0, 2], (5)
where the i are realizations of the E i , i = 1, ..., n.

For relative error models one should use inverse problem formulations with Generalized Least Squares (GLS) cost functional

J n (Y ; θ) = n i=1 Y i -f (t i ; θ) f (t i ; θ) γ 2 . ( 6
)
The corresponding estimator and estimates are respectively defined by

Θ GLS = argmin θ∈Ω n i=1 Y i -f (t i ; θ) f (t i ; θ) γ 2 , γ ∈ [0, 2], (7)
with realizations

θGLS = argmin θ∈Ω n i=1 y i -f (t i ; θ) f (t i ; θ) γ 2 , γ ∈ [0, 2]. (8)
GLS estimates θn and weights {ω j } n j=1 are found using an iterative method as defined below (see [START_REF] Banks | Modeling and Inverse Problems in the Presence of Uncertainty[END_REF]). For the sake of notation, we will suppress the superscript n (i.e., θGLS := θn GLS ).

1. Estimate θGLS by θ(0) using the OLS method (( 8) with γ = 0 ). Set k = 0.

Compute weights ωj

= f -2γ (t j , θ (k)
).

3. Obtain the k + 1 estimate for θGLS by θ(k+1)

:= argmin n j=1 ωj [y j - f (t j , θ)] 2 .
4. Set k := k + 1 and return to step 2. Terminate when the two successive estimates for θGLS are sufficiently close.

Mathematical Models

Our focus here is on the comparison of two different models for insect (L. hesperus) population growth/mortality in pesticide-treated fields. The simplest model (which we denote as model B) is for constant reduced growth due to effects of pesticides versus an added time-varying mortality (denoted by model A) to reflect this decreased total population growth rate. Model B is given by

dx dt = ηx x(t 1 ) = x 0 , (9)
where x 0 is defined as initial population count at time t 1 of initial observation and η is the reduced population growth rate in the presence of pesticides. Model A is given by

dx dt = k(t)x x(t 1 ) = x 0 , (10)
where t 1 is again the time of the first data point, and k(t) is a time dependent growth rate

k(t) = η + p(t) t ∈ P j , j ∈ {1, 2, 3, 4} η otherwise.
Here p(t) is composed of piecewise linear splines as described below, and P j = [t pj , t pj + 1/4], j = 1, . . . j * with t pj as the time point of the j th pesticide application. Observe that these t pj are not the same as the observation or data points t j . Also note that |P j | = 1/4 which is approximately the length of time of one week when t is measured in months. This reflects the general assumption that pesticides are most active during the 7 days immediately following treatment.

Clearly, η is a reduced constant growth rate of the total population in the presence of pesticides. In addition, t=0 refers to June 1 (as no data is present before June 1 in our database). Piecewise linear splines [START_REF] Schultz | Spline Analysis[END_REF] were used to approximate p(t) as follows. Consider m linear splines

p(t) = m i=1 λ i l i (t),
where

l i (t) = 1/h      t -t i-1 t i-1 ≤ t < t i t i+1 -t t i ≤ t ≤ t i+1 0 otherwise,
where h is the step size, h = |Pj | (m+1) . Piecewise linear spline representations are simple, yet flexible in that they allow the modeler to avoid assuming a certain shape to the curve being approximated. Incorporating a time-dependent component such as p(t) is useful when modeling a system with discontinuous perturbations (such as the removal of a predator, or the application of an insecticide). The addition of more splines (m > 3) provides a finer approximation, but demands more terms in the parameter estimates. We assume here that it is likely that m = 3 is sufficient. (Our subsequent findings suggest that perhaps even m = 2 is sufficient!). In our analysis, we first estimated the initial condition x 0 using model B (as this data point precedes any pesticide applications and provides a good estimate for x 0 ), and then fixed this parameter in all subsequent parameter estimates. Therefore, the parameters to be estimated in model A are θ = q = {η, λ 1 , λ 2 , λ 3 } whereas for model B we must only estimate θ = q = {η, 0, 0, 0} since model A reduces to model B when applying the constraint p(t) ≡ 0, i.e. λ i = 0 for i = 1, 2, 3.

Parameter Estimation

Using the model information provided in [START_REF] Banks | Modeling populations of Lygus Hesperus on cotton fields in the San Joaquin Valley of California: The importance of statistical and mathematical model choice[END_REF] we try to estimate parameters for new data sets and determine whether the fit-to-data provided by model A does provide a statistically significantly better fit than the fit provided by model B. A big part of the parameter estimation process is the minimization of the respective cost functions for both model A and model B. The constrained nonlinear optimization solver in Matlab, fmincon was initially being used for minimization of GLS cost functionals in model A while fminsearch was being used for minimization of cost functionals for model B. We later switched to lsqnonlin which gave faster and better results. Since model A is stiff in nature, Matlab solver ode15s was used whereas for model B ode45 was used. Both fmincon and fminsearch require an initial guess of parameters. While for model B fminsearch was able to find a minimum fairly quickly, the initial guess of θ = {η, λ 1 , λ 2 , λ 3 } for model A involved a fairly detailed process given below:

1. Create a trial file of data selected based on a specific set of rules. 1.5]. (Note that best way to choose (see [START_REF] Banks | Modeling and Inverse Problems in the Presence of Uncertainty[END_REF][START_REF] Banks | Mathematical and Experimental Modeling of Physical and Biological Processes[END_REF]) gamma is to consider the plot of residuals using both residual vs time and residual vs model plots to ascertain whether the scatter of the error appears to violate the statistical assumption of being i.i.d.). 4. Choose an initial condition x 0 to the exponential model by considering the plot of the data and model and visually estimating x 0 . We observe that this did not produce acceptable results so we eventually solved a separate inverse problem to estimate intrinsic growth rate ignoring the effects of pesticides and initial condition.

Choose a parameter space

Ω = [ , K] × [-K, -] 3 . 3. Choose a constant γ ∈ [0,

Model Comparsion: Nested Restraint Sets

Here we summarize the use of statistically based model comparison tests. These residual sum of squares model comparison tests as developed in [START_REF] Banks | Statistical methods for model comparison in parameter estimation problems for distributed systems[END_REF], described in [START_REF] Banks | Modeling and Inverse Problems in the Presence of Uncertainty[END_REF][START_REF] Banks | Mathematical and Experimental Modeling of Physical and Biological Processes[END_REF] and extended in [START_REF] Banks | An extension of RSS-based model comparison tests for weighted least squares[END_REF] to GLS problems is used in the same manner as used in [START_REF] Banks | Model comparison tests to determine data information content[END_REF][START_REF] Banks | Modeling populations of Lygus Hesperus on cotton fields in the San Joaquin Valley of California: The importance of statistical and mathematical model choice[END_REF]. This test is used to determine which of several nested models is the best fit to the data; therefore, this test can be applied to the comparison of models A and B. In these examples below we are interested in questions related to whether the data will support a more detailed or sophisticated model to describe it. In the next section we recall the fundamental statistical tests to be employed here.

Statistical Comparison Tests

In general, assume we have an inverse problem for the model observations f (t, θ) and are given n observations. As in (6), we define

J n (Y ; θ) = n i=1 Y i -f (t i ; θ) f (t i ; θ) γ 2 ,
where our statistical model has the form (4). Here, as before, θ 0 , is the nominal value of θ which we assume to exist. We use Ω to represent the set of all the admissible parameters θ. We make some further assumptions.

-A4) Observations are taken at {t j } n j=1 in [0, T ]. There exists some finite measure µ on [0, T ] such that

1 n n j=1 h(t j ) -→ T 0 h(t)dµ(t)
as n → ∞, for all continuous functions h.

-A5) J 0 (θ) = T 0 (f (t; θ 0 ) -f (t; θ)) 2 dµ(t) = σ 2 has a unique minimizer in Ω at θ 0 . Let Θ n = Θ n
GLS (Y ) be the GLS estimator for J n as defined in ( 7) so that

Θ n GLS (Y ) = argmin θ∈Ω J n (Y ; θ) and θn GLS = argmin θ∈Ω J n (y; θ),
where as above y is a realization for Y . One can then establish a series of useful results (see [START_REF] Banks | Statistical methods for model comparison in parameter estimation problems for distributed systems[END_REF][START_REF] Banks | Modeling and Inverse Problems in the Presence of Uncertainty[END_REF][START_REF] Banks | An extension of RSS-based model comparison tests for weighted least squares[END_REF] for detailed proofs).

Result 1: Under A1) to A5),

1 n Θ n = 1 n Θ n GLS (Y ) -→ θ 0 as n → ∞ with probability 1.
We will need further assumptions to proceed (these will be denoted by A7)-A11) to facilitate reference to [START_REF] Banks | Statistical methods for model comparison in parameter estimation problems for distributed systems[END_REF][START_REF] Banks | Modeling and Inverse Problems in the Presence of Uncertainty[END_REF]). These include:

-A7) Ω is finite dimensional in R p and θ 0 ∈ intΩ. -A8) f : Ω → C[0, T ] is a C 2 function. -A10) J = ∂ 2 J0 ∂θ 2 (θ 0 ) is positive definite. -A11) Ω H = {θ ∈ Ω|Hθ = c}
where H is an r × p matrix of full rank, and c is a known constant.

In many instances, including the examples discussed here, one is interested in using data to question whether the "nominal" parameter θ 0 can be found in a subset Ω H ⊂ Ω which we assume for discussions here is defined by the constraints of assumption A11). Thus, we want to test the null hypothesis H 0 : θ 0 ∈ Ω H , i.e., that the constrained model provides an adequate fit to the data. Define then

Θ n H (Y ) = argmin θ∈Ω H J n (Y ; θ) and θn H = argmin θ∈Ω H J n (y; θ).
Observe that J n (y; θn H ) ≥ J n (y; θn ). We define the related non-negative test statistics and their realizations, respectively, by

T n (Y ) = J n (Y ; θ n H ) -J n (Y ; θ n ) and Tn = T n (y) = J n (y; θn H ) -J n (y; θn ).
One can establish asymptotic convergence results for the test statistics T n (Y )see [START_REF] Banks | Statistical methods for model comparison in parameter estimation problems for distributed systems[END_REF]. These results can, in turn, be used to establish a fundamental result about much more useful statistics for model comparison. We define these statistics by

U n (Y ) = nT n (Y ) J n (Y ; θ n ) , (11)
with corresponding realizations ûn = U n (y).

We then have the asymptotic result that is the basis of our analysis-ofvariance-type tests.

Results 2: Under the assumptions A1)-A5) and A7)-A11) above and assuming the null hypothesis H 0 is true, then U n (Y ) converges in distribution (as n → ∞) to a random variable U (r), i.e.,

U n D -→ U (r),
with U (r) having a chi-square distribution χ 2 (r) with r degrees of freedom.

In any graph of a χ 2 density there are two parameters (τ, α) of interest. For a given value τ , the value α is simply the probability that the random variable U will take on a value greater than τ . That is, Prob{U > τ } = α where in hypothesis testing, α is the significance level and τ is the threshold.

We then wish to use this distribution U n ∼ χ 2 (r) to test the null hypothesis, H 0 , that the restricted model provides an adequate fit to represent the data. If the test statistic, ûn > τ , then we reject H 0 as false with confidence level (1α)100%. Otherwise, we do not reject H 0 . For our examples below, we use a χ 2 (3) table, which can be found in any elementary statistics text, online or the partial summary below. Typical confidence levels of interest are 90%, 95%, 99%, 99.9%, with corresponding (α, τ ) values given in Table 1 below.

To test the null hypothesis H 0 , we choose a significance level α and use χ 2 tables to obtain the corresponding threshold τ = τ (α) so that Prob{χ 2 (r) > τ } = α. We next compute ûn = τ and compare it to τ . If ûn > τ , then we reject H 0 as false; otherwise, we do not reject the null hypothesis H 0 .

We use a χ 2 (3) for our comparison tests as summarized here.

We can then formulate the null and alternative hypotheses: H 0 : The fit provided by model A does not provide a statistically significantly better fit to the data than the fit provided by model B. H A : The fit provided by model A does provide a statistically significantly better We carried out multiple inverse problems with varying values of for these data sets (see [START_REF] Banks | Modelling pesticide treatment effects on Lygus hesperus in cotton fields[END_REF]). We visually examined the resulting residual plots (residual vs time and residual vs observed output) and determined whether the scatter of the error appear to be i.i.d.. On examining the plots for a wide range of we observed that the statistical i.i.d assumptions were approximately satisfied for values ranging around 0.7 to 0.8. We therefore used these values of in the results reported here and in [START_REF] Banks | Modelling pesticide treatment effects on Lygus hesperus in cotton fields[END_REF].

We note here that one could, as an alternative to use of residual plots, include the parameter as a parameter to be estimated along with as often done in statistical formulations [START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF] for the joint estimation of = (, ), or one could attempt to estimate the form of the statistical model (i.e., the value of directly from the data itself as suggested in [START_REF] Banks | Applications of difference-based variance estimation methods[END_REF]. Both of these methods offer some advantages but include more complex inverse problem analysis. We have therefore chosen here to use the simpler but less sophisticated analysis of residuals in our approach.

Results

We examine the importance of model selection and how optimal selection of both statistical and mathematical models is crucial for accuracy in parameter estimation and fitting data to models testing different data sets from the same database with varied number of pesticide applications. Tables 2, 3 and 4 contain summaries of the results for the investigated replicates.

Concluding Remarks

The above results strongly support the notion that time varying reduced growth/mortality rates as opposed to constant rates provide substantially better models at the population levels for the description of the effects of pesticides on the growth rates.

It is interesting to note that our findings hold consistently across the differing levels of pesticide applications, even in the case of only one pesticide application. Another interesting observation is that this is consistent even when one uses total pest counts as opposed to individual nymph and adult counts. This, of course, has significant implications for data collection procedures. The model comparison techniques employed here are just one of several tools that one can use to determine aspects of information content in support of model sophistication/complexity. Of note are the use of the Akiake Information Criterion (AIC) and its variations [START_REF] Banks | Model comparison tests to determine data information content[END_REF][START_REF] Banks | Information content in data sets for a nucleated-polymerization model[END_REF][START_REF] Bozdogan | Model selection and Akaike's information criterion (AIC): The general theory and its analytical extensions[END_REF][START_REF] Bozdogan | Akaike's information criterion and recent developments in information complexity[END_REF][START_REF] Burnham | Model Selection and Inference: A Practical Information-Theoretical Approach[END_REF][START_REF] Burnham | Multimodel inference: understanding AIC and BIC in model selection[END_REF] for model comparison.

Table 1 :

 1 Chi-Square Table: χ 2 (3)

	α .10 .05	.01	.001
	τ 6.251 7.815 11.345 16.266
	fit to the data than the fit provided by model B.
	We considered such comparison tests for a number of data sets with a varying
	no. of pesticides applications among the fields. These included
	-Replicate number 296 (1 pesticide application at t = 0.5 months)
	-Replicate number 350 (2 pesticide applications at times t = 1 and t = 1.7
	months)
	-Replicate number 277 (3 pesticide applications at t = .4, .7, 1.9 months)
	-Replicate number 178 (4 pesticide applications at t = .13, .77, 1.9, 2.27
	months)
	-Replicate number 174 (4 pesticide applications at t = .13, .77, 1.83, 2.27
	months)

Table 2 :

 2 Parameter estimates for models A and B

	Replicate No. No. of pest. apps	Par. Est. for model A (η, λ 1 , λ 2 , λ 3 )	Par. Est. for model B (η, λ 1 , λ 2 , λ 3 )
	296 (n=19)	1	(1.3551, -7.8782, -7.1632, -9.9904)	(0.3548, 0, 0, 0)
	350 (n=19)	2	(3.0150, -16.7393, -17.9356, -8.580)	(0.7683, 0, 0, 0)
	277 (n=19)	3	(3.1255, -3.5165, -6.2870, -7.1868)	(1.7284, 0, 0, 0)
	178 (n=22)	4	(1.8201, -6.6352, -2.5786, -4.1605)	(0.6335, 0, 0, 0)
	174 (n=20)	4	(2.8600, -7.9867, -13.2780, -4.6992)	(0.4110, 0, 0, 0)

Table 3 :

 3 Cost Functional values for models A and B

	Replicate No.	Cost functional Model A	Cost functional Model B	Range of γ values x 0
	296	3.9577	7.3663	0.5 to 1.0	0.39
	350	7.4100	18.223	0.5 to 0.8	0.28
	277	5.5290	8.075	0.7 to 0.9	0.05
	178	5.7596	8.371	0.7 to 1.0	0.37
	174	6.7480	9.390	0.0 to 0.8	0.34

Table 4 :

 4 Statistical comparison test results for degrees of freedom r = 3

	Replicate No.	Test statistic ûn	Threshold τ	Signif. Level α	Null hyp. (H 0 ) Accept or reject	Confid. Level
	296	16.36	16.26	0.001	Rejected	99.9 %
	350	27.73	16.26	0.001	Rejected	99.9 %
	277	8.7458	7.815	0.05	Rejected	95 %
	178	9.9715	7.815	0.05	Rejected	95 %
	174	7.8281	7.815	0.05	Rejected	95 %

Acknowledgements

This research was supported by the Air Force Office of Scientific Research under grant number AFOSR FA9550-12-1-0188 and grant number AFOSR FA9550-15-1-0298. J.E. Banks is grateful to J.A. Rosenheim for hosting him during his sabbatical, and also extends thanks to the University of Washington, Tacoma for providing support for him during his sabbatical leave.


Alfred Schmidt

Eberhard Bänsch
email: baensch@math.fau.de

Mischa Jahn

Andreas Luttmann

Carsten Niebuhr

Jost Vehmeyer
email: vehmeyer@math.uni-bremen.de

Optimization of Engineering Processes Including Heating in Time-Dependent Domains

Keywords: Optimization with PDEs, time-dependent domain, heat equation, thermoelasticity, free surface flow

We present two models for engineering processes, where thermal effects and time-dependent domains play an important role. Typically, the parabolic heat equation is coupled with other equations. Challenges for the optimization of such systems are presented. The first model describes a milling process, where material is removed and heat is produced by the cutting, leading to thermomechanical distortion. Goal is the minimization of these distortions. The second model describes the melting and solidification of metal, where the geometry is a result of free-surface flow of the liquid and the microstructure of the re-solidified material is important for the quality of the produced preform.

Introduction

The optimization of industrial engineering processes often lead to the treatment of coupled, nonlinear systems of PDEs. Here we want to investigate applications where an important part of the nonlinearities is created by time-dependent domains, which are not prescribed but who are part of the solution itself. The treatment of such models in the context of PDE constraints in optimal control problems typically generates additional challenges in both, the solution of the forward problems and the treatment and storage of adjoint solutions.

In the following, we present two models for engineering processes, where heating of metal workpieces and time-dependent domains play an important role. Thus, a parabolic heat equation is coupled with other equations. Challenges for the optimization of such coupled systems are presented.

The first model describes a milling process, where material is removed and heat is produced by the cutting. This leads to a thermomechanical distortion of corresponding author the workpiece during the cutting process and leads to an incorrect removal of material. An optimization of the cutting path and speed, varying chip thickness and thus heat production, etc., should give reduced distortion during the process and lead to a correct workpiece shape.

The second model describes the melting and solidification of metal heated by a laser beam. Due to free-surface flow, the shape of the liquid part depends on capillary boundary conditions, and heat transport on the flow field. The microstructure of the re-solidified material, which is important for subsequent process steps, depends on the temperature gradients near the moving liquid-solid interface and its velocity. Accelerating the process for mass production on the one hand and improving the microstructure on the other hand compete for an optimized process.

Both optimization problems can be formulated in an abstract setting as min u∈U ad J(u, y) under the constraint y = S(u)

where u denotes the control, U ad the set of admissible controls, y the state, J the error functional to be minimized, and S is the control-to-state operator, given by a nonlinear PDE. J is typically given by a deviation of y (or something derived from it) from a desired function y d plus some regularization by a norm of u, like

J(u, y) = d(y, y d ) + α u p p . (2)
For both applications, we will first state and describe the primal problem giving the solution operator S, and later cover some details of the associated optimization problem.

Thermomechanics of Milling Processes

Application

During a milling process, heat is produced by the cutting tool and transfered into the workpiece, and mechanical load is generated by cutting forces. Due to the resulting thermomechanical deformation of the workpiece, the final shape of the processed workpiece deviates from the desired shape, making a postprocessing finishing necessary. Deformations are relatively large especially when producing fine structures like thin walls for lightweight constructions. In order to reduce the shape deviation, an optimization of the tool path and other process parameters is desirable, taking into account the thermomechanical deformations.

Model

The mathematical model for the thermomechanics of the process includes thermoelasticity of the workpiece, energy and forces introduced by the process, and most importantly the cutting process itself, which leads to a time-dependent domain, whose shape influences the process and vice versa. As typical under the assumption of small deformations, the model is formulated in a reference configuration.

Let Ω(t) ⊂ R 3 denote the time-dependent domain in the reference configuration, Q T := {(x, t) : x ∈ Ω(t), t ∈ (0, T )} the space-time-domain, θ : Q T → R the temperature, and v : Q T → R 3 the deformation.

The change of geometry by material removal as well as energy and forces introduced by the cutting process are provided by a process model [START_REF] Denkena | Identification of the specific cutting force for complex cutting tools under varying cutting conditions[END_REF][START_REF] Denkena | Prediction of Temperature Induced Shape Deviations in dry Milling[END_REF], taking into account the tool path and velocity, chip thickness, temperature and deformation, and other global and local parameters and properties. Here, we rely on a macroscopic model where microscopic processes like chip formation are not directly considered, but their effects considered via the process model. Let us denote by Γ (t) ⊂ ∂Ω(t) the contact zone of the cutting tool at time t, Γ T := {(x, t) : x ∈ Γ (t), t ∈ (0, T )}, q Γ : Γ T → R the normal heat flux, and g Γ : Γ T → R 3 the forces introduced at the cutting surface.

In the notion of optimal control problems, the state y consists of the domain, temperature, and deformation, y = (Ω, θ, v), while the control u consists of the process parameters like tool path, feed rate, rotational velocity, etc. The material removal and thus the domain Ω(t) depend on the cutting process (control u) and the deformation v.

The coupled model includes the parabolic heat equation and quasistatic, elliptic thermoelasticity

θ -∇ • (κ∇θ) = 0, (3)
-∇ • σ = f v (θ) (4)
on Ω T with stress tensor σ = 2µDv + λtr(Dv)I and Dv = 1 2 (∇v + ∇v T ) the symmetric gradient or strain tensor. On the contact zone Γ T , we have boundary conditions for heat flux and mechanical forces given by the process model,

κ∇θ • n = q Γ (u, θ, v), (5)
σn = g Γ (u, θ, v). ( 6
)
The workpiece is clamped, which is reflected by Dirichlet conditions on a subset

Γ D ⊂ ∂Ω \ Γ T , v = 0 on Γ D , (7)
while cooling conditions and free deformation apply on the rest of the boundary,

κ∇θ • n = r(θ ext -θ), (8)
σn = 0. (9)
Initially, the temperature is typically constant at room temperature, thus θ = θ 0 on Ω(0).

Numerical Discretization of the Forward Problem

The system (3-9) of PDEs is discretized using a finite element method on an adaptively locally refined tetrahedral mesh [START_REF] Niebuhr | Finite Element Simulation of macroscopic Machining Processes -Implementation of time-dependent Domain and Boundary Conditions, 14[END_REF][START_REF] Schmidt | Design of adaptive finite element software: The finite element toolbox ALBERTA[END_REF], using piecewise polynomial functions for the temperature and the components of deformation. Time-discretization is based on a semi-implicit time stepping scheme.

The time-dependent domain Ω(t) is approximated by a subset Ω h (t) of the triangulation, where the completely cut off elements are ignored. The cutting process is simulated by a dexel method [START_REF] Niederwestberg | Contact Zone Analysis Based on Multidexel Workpiece Model and Detailed Tool Geometry[END_REF], which is able to compute the interaction of the tool with the (deformed) workpiece very efficiently, giving Ωh (t) and Γh (t). At the same time, chip thickness and other cutting parameters are computed and the process model returns approximations of the heat flux q Γ,h and forces g Γ,h at the cutting surface Γh (t). Based on that, the finite element method computes Ω h (t) and Γ h (t) and projects the boundary data onto Γ h (t). Finite element approximations of temperature and deformation are computed on Ω h (t). This is done in every time step of the finite element method. The overall method is described in [START_REF] Denkena | Modeling a Thermomechanical NC-Simulation[END_REF][START_REF] Denkena | Prediction of Temperature Induced Shape Deviations in dry Milling[END_REF].

Figure 1 shows the mesh, temperature, and deformation from the simulation of a milling process. The mesh is adaptively refined in order to approximate the geometry Ω(t) well by Ω h (t). The process removes layers of material to mill a pocket into a rectangular bar of metal. Especially the final thin backward wall is prone to deformations larger than the given tolerance, making it hard to produce the desired shape.

Optimization

The goal of a process optimization is a compensation and reduction of geometry errors while not slowing down the process too much. Control is given by the variation of process parameters like tool path and velocity, speed of tool rotation, etc. Such variations result in changes of the tool entry situation and thus, reflected by the process model, changes in heat source and cutting forces, and finally a change in the thermomechanic deformation. Admissible controls in U ad are defined by restrictions on the machining process.

Given a prescribed final geometry Ω d of the workpiece, the optimization functional should include the deviation of the process geometry from the desired one, as well as the process duration. Considering the geometry error, different criteria are possible, especially comparing geometries during the whole cutting process or only in the end. For the latter, this would nevertheless include geometry error terms over time, as material which was removed before cannot be added later on again.

Another approach to an error functional includes geometry deviations near the cutting zone Γ T during the whole process:

J(u, Ω, θ, v) = T 0 Ω(t) -Ω d (t) 2 Γ (t) + λ u 2 . ( 10
)
We show the effect of such an error functional in a model situation, where a L-shaped geometry is produced from a rectangular plate. This can be seen as a slice through the original workpiece, see the left of Figure 2. The non-optimized control does not consider the thermal extension and leads to increased material removal resulting in a recessed surface after the workpiece has cooled down. Figure 2 shows on the right the geometry error in the contact zone over the process time. The general optimal control problem requires to find the spatial tool po- sitions and the cutting parameters which minimize [START_REF] Niebuhr | Finite Element Simulation of macroscopic Machining Processes -Implementation of time-dependent Domain and Boundary Conditions, 14[END_REF]. In a first investigation of the model problem, simple raster milling is performed and control is given by traditional setting parameters, i.e. cutting depth, radial and tangential feed and cutting velocity. Figure 3 shows the resulting surface for the non-optimized process and for the process with optimal parameters. Adjoint problem. The optimization shown above was done with the standard MATLAB optimization toolbox which just calls the finite element package to solve the forward problem. For a more involved algorithm, the computation and storage of the adjoint solution would be used. The corresponding system of adjoint problems consists of the adjoint (backward) heat equation on the time-dependent domain (now given from the forward problem), coupled to the quasistatic adjoint elasticity equation. Due to the rather long process time, the three-dimensional time-dependent domain, and adaptive meshes for approximation of domain and solution, the handling of such adjoint solutions in the optimization procedure is a challenge.

3 Material Accumulation by Laser Heating

Application

For the production of micro components (like micro-valves, etc., with diameters smaller than 1mm) by cold forming, a necessary pre-forming step is to accumulate enough material for a subsequent cold forming step. This can be done by partial melting and solidification of a half-finished product like a thin wire [START_REF] Vollertsen | Energy balance in laser based free form heading[END_REF]. Due to the small scale, the dominant surface tension of the melted material leads to a nearly spherical form, leading to an accumulated solid sphere attached to the wire after solidification which is called preform. Due to the industrial background, very high process speeds are requested. However, for the subsequent forming step, the microstructure of the material is important. Thus, besides the speed of the process and an accurate size of the preform, its microstructure is part of the optimization goal. Formation of dendritic structures and their spacings, or other phases, during the solidification are strongly influenced by the liquid-solid interface velocity and local temperature gradients [START_REF] Kurz | Fundamentals of Solidification[END_REF]. Thus, temperature, phase transitions, and the geometry are important aspects of the corresponding optimization problem.

Model

We consider the time-dependent domain Ω(t) consisting of the solid and partially melted parts of the material (metal). Melting and solidification are typically modeled by the Stefan problem, including temperature θ and energy density e as variables in the space-time-domain

Q T := {(x, t) : x ∈ Ω(t), t ∈ (0, T )}: ė + v • ∇e -∇ • (κ∇θ) = 0, θ = β(e), (11)
where β(s) := c 1 min(s, 0) + c 2 max(0, s -L) and L denotes the latent heat of solid-liquid phase transitions. β is only a piecewise-smooth function with a constant part, making [START_REF] Niederwestberg | Contact Zone Analysis Based on Multidexel Workpiece Model and Detailed Tool Geometry[END_REF] a degenerate parabolic equation.

The liquid subdomain Ω l (t) is given by all points where the temperature is above the melting temperature (which is assumed to be 0 after some scaling),

Ω l (t) := {x ∈ Ω(t) : θ(x, t) > 0}. ( 12
)
The shape of the melted (and later on re-solidified) material accumulation is mainly influenced by the surface tension of the liquid, together with gravitational forces etc., which means free-surface flow. Parabolic Navier-Stokes equations with capillary boundary condition is the main model component for this, with solenoidal velocity field v and pressure p in

Ω l (t), v + v • ∇v -∇ • σ = f v (θ), ∇ • v = 0, (13)
with stress tensor σ = 1 Re Dv -pI. Here, f v denotes the forces introduced by gravity due to a temperature-dependent density and the Boussinesq approximation. The shape of the liquid subdomain is given through the capillary boundary condition on the free surface Γ (t) of the melted subdomain, where the surface tension (proportional to the mean curvature of the surface) balances the normal stress. In a differential geometric PDE formulation, the mean curvature vector H Γ is given by the Laplace-Beltrami-operator -∆ Γ applied to the coordinates of the surface (represented by the embedding id : Γ (t) → R 3 ), giving another (nonlinear) elliptic equation in the coupled system. Additionally, the normal component of the fluid velocity should be equal to the normal velocity V Γ of the capillary surface. Both lead to the following equations on Γ (t):

σn = 1 W e H Γ = - 1 W e ∆ Γ id, v • n = V Γ . (14)
On the solid-liquid interface and in the solid subdomain, the velocity vanishes,

v = 0 in Ω(t) \ Ω l (t). ( 15
)
The heating is done through a laser pointing at a spot on the boundary, modeled by a time-and space-dependent energy density q L , and cooling conditions apply on the whole boundary,

κ∇θ • n = q L + r(θ ext -θ) on ∂Ω(t). (16)
A more detailed description of the model can be found in [START_REF] Jahn | Finite element simulation of a material accumulation process including phase transitions and a capillary surface[END_REF].

Numerical Discretization of the Forward Problem

The forward problem with given boundary conditions is discretized by a finite element method which combines a Stefan problem solver with a free-surface Navier-Stokes solver. The latter is based on the Navier code [START_REF] Bänsch | Finite element discretization of the Navier-Stokes equations with a free capillary surface[END_REF], the combined approach is described in more detail in [START_REF] Bänsch | An ALE finite element method for a coupled Stefan problem and Navier-Stokes equations with free capillary surface[END_REF][START_REF] Jahn | Finite element methods for problems with solid-liquid-solid phase transitions and free melt surface[END_REF]. Locally refined (triangular or tetrahedral) meshes are needed in order to approximate the large variations in temperature near the heating zone and the surrounding of the solid-liquid interface sufficiently well, while keeping the overall numerical costs acceptable. Due to the big changes in geometry, starting from a thin wire and ending in a relatively large spherical accumulation, several remeshings are necessary during the simulation in order to avoid a degeneration of mesh elements.

Figure 4 shows a typical mesh, temperature field, and velocity field during the melting, with the laser pointing to the center bottom of the material. Due to rotational symmetry, a 2D FEM with triangular meshes could be used. As the wire is melted from below, the growing sphere is moving upwards and thus the velocity vectors are pointing upwards, too. In Figure 5, we show several stages during the solidification after the heating is switched off (from a simulation with shorter heating period). In contrast to our results for energy dissipation during the heating process [START_REF] Jahn | Energy dissipation in laserbased free form heading: a numerical approach[END_REF], the cooling after switching off the laser is mainly done by conducting heat upwards into the wire, which results in a downward movement of the interface. The additional cooling by the boundary conditions leads later on to an additional solidification at the boundary of the liquid region. During solidification, the shape does not change much anymore, so the velocities are typically quite small and not shown here.

Optimization

The goals for an optimization of the process are on the one hand a high speed in order to be able to produce thousands of micro preform parts in a short time, while generating material accumulations of the desired size, and on the other hand generating a microstructure which is well suited for the subsequent cold forming step and usage properties of the work piece. Models for microstructure evolution during solidification show a dependency of the microstructure quality on the speed V Γ I of the liquid-solid interface Γ I (t) and the temperature gradients near the interface [START_REF] Kurz | Fundamentals of Solidification[END_REF]. A higher speed and larger gradient typically result in a microstructure which gives better forming characteristics and useful properties. Thus, the error functional has to include parts for geometry approximation, for speed, and for the microstructure generation during solidification,

J(u, Ω, θ, v) = T 0 Ω(t) -Ω d (t) 2 (17) +λ 1 T 0 V Γ I -V d 2 Γ I (t) + ∇θ -G d 2 Γ I (t) + λ 2 u 2 .
The control parameters u are given by the time-dependent intensity and location of the laser spot, thus they enter the system of equations via the heating energy density q L (u). An additional control variable might be the intensity r of outer cooling, for example by adjusting the flow velocity of a cooling gas. U ad is given by control restrictions which follow from limits of the process, like an upper bound of the laser power for preventing an evaporation of the material.

Adjoint Problem. For an efficient implementation of the optimization procedure, the solution of the adjoint problem will be used. Here, the system of adjoint equations includes the adjoint (linearized) Navier-Stokes system on the prescribed time dependent domain (from the forward solution), together with the (linear) Laplace-Beltrami equation on the prescribed capillary boundary. For the formulation of the adjoint Stefan problem, a regularization could be used which was used in [START_REF] Nochetto | A posteriori error estimation and adaptivity for degenerate parabolic problems[END_REF] for the derivation of an aposteriori error estimate.

As in our first application, also here the efficient handling of the adjoint problem including adaptively refined meshes with remeshings for the time dependent domain and the corresponding system of solutions in the domain and on the capillary surface poses a challenge for the overall numerical optimization method.

Fig. 1 .

 1 Fig. 1. Simulation of a milling process: mesh and temperature (top) and deformation (bottom, amplified by a factor 100). The tool is at the moment cutting near the backward left lower corner.

Fig. 2 .

 2 Fig. 2. Model geometry (left), deformation v in the cutting zone at t = 10.5s (middle), and maximal deformation in the cutting zone Γ (t) over time (right).

Fig. 3 .

 3 Fig. 3. L-shaped workpiece near the end of the milling process, with initial control (left) and optimized parameters (right). Colors depict the modulus of deformation.

Fig. 4 .

 4 Fig. 4. Melting the end of a thin wire: liquid material accumulation with mesh, solidliquid interface, and liquid flow field.

Fig. 5 .

 5 Fig. 5. Material accumulation at the end of a thin wire: Solid-liquid interface and isothermal lines at 5 different times during solidification.

Acknowledgments. The authors gratefully acknowledge the financial support by the DFG (German Research Foundation) for the subproject A3 within the Collaborative Research Center SFB 747 "Micro cold forming" and the project MA1657/21-3 within the Priority Program 1480 "Modeling, Simulation and Compensation of thermal effects for complex machining processes".

Furthermore, we thank the Bremen Institute for Applied Beam Technology (BIAS) and the Institute of Production Engineering and Machine Tools Hanover (IFW) for cooperation.


Vyacheslav Maksimov

Game Control Problem for Systems of Distributed Equations

Keywords:

We consider a game problem of guaranteed positional control for a distributed system described by the phase field equations under incomplete information on system's phase states. This problem is investigated from the viewpoint of the first player (the partner). For this player, a procedure for forming feedback controls is specified. This procedure is stable with respect to informational noises and computational errors and is based on the method of extremal shift and the method of stable sets from the theory of guaranteed positional control. It uses the idea of stable dynamical inversion of controlled systems.

Introduction

The control theory for distributed systems has been intensively developed in recent time as a part of mathematical control theory. At present, there exists a number of monographs devoted to control problems for distributed systems [START_REF] Fattorini | Infinite Dimensional Optimization and Control Theory[END_REF][START_REF] Lasiecka | Control theory for partial differential equations: continuous and approximation theories[END_REF][START_REF] Tröltzsch | Optimal Control of Partial Differential Equations. Theory, Methods and Applications[END_REF]. As a rule, the emphasis is on program control problems in the case when all system's parameters are precisely specified. Along with this, the investigation of control problems for systems with uncontrollable disturbances (the problems of game control) is also reasonable. Similar problems have been poorly investigated. In the early 70es, N.N.Krasovskii suggested an effective approach to solving game (guaranteed) control problems, which is based on the formalism of positional strategies. The detailed description of this approach for dynamical systems described by ordinary differential equations is given in [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF]. The goal of the present work is to illustrate possibilities of this approach for investigating a game problem for systems described by the phase field equations.

We consider a system modeling the solidification process and governed by the phase field equations (introduced in [START_REF] Caginalp | An Analysis of a Phase Field Model of a Free Boundary[END_REF])

∂ ∂t ψ + l ∂ ∂t ϕ = k∆ L ψ + Bu -Cv in Q = Ω × (t 0 , ϑ], τ ∂ ∂t ϕ = ξ 2 ∆ L ϕ + g(ϕ) + ψ, ϑ = const < +∞, (1)
This work was supported by the Russian Science Foundation (project 14-01-00539).

with the boundary condition ∂ ∂n ψ = ∂ ∂n ϕ = 0 on ∂Ω × (t 0 , ϑ] and the initial condition ψ(t 0 ) = ψ 0 , ϕ(t 0 ) = ϕ 0 in Ω. Here, Ω ⊂ R n is a bounded domain with the sufficiently smooth boundary ∂Ω, ∆ L is the Laplace operator, ∂/∂n is the outward normal derivative, (U, | • | U ) and (V, | • | V ) are Banach spaces, B ∈ L(U ; H) and C ∈ L(V ; H) are linear continuous operators, H = L 2 (Ω), and the function g(z) is the derivative of a so-called potential G(z). We assume that g(z) = az + bz 2 -cz 3 .

Systems of form [START_REF] Fattorini | Infinite Dimensional Optimization and Control Theory[END_REF] have been investigated by many authors. In what follows, for the sake of simplicity, we assume that k = ξ = τ = c = 1. Further, we assume that the following conditions are fulfilled: (A1) the domain Ω ⊂ R n , n = 2, 3, has the boundary of C 2 -class; (A2) the coefficients a and b are elements of the space L ∞ (T × Ω), T = [t 0 , ϑ], and vrai sup c(t, η) > 0 for (t, η) ∈ [t 0 , ϑ] × Ω; (A3) the initial functions ψ 0 and ϕ 0 are such that {ψ | • | H are the scalar product and the norm in H, respectively. Let some initial state

0 , ϕ o } ∈ R = {ψ, ϕ ∈ W 2 ∞ (Ω) : ∂ ∂n ψ = ∂ ∂n ϕ = 0 on ∂Ω}. Introduce the notation: W 2,1 p (Q) = u | u, ∂u ∂ηi , ∂ 2 u ∂ηi∂ηj , ∂u ∂t ∈ L p (Q) for p ∈ [1, ∞) is the standard Sobolev space with the norm u W 2,1 p (Q) ; (•, •) H and
x 0 = {ψ 0 , ϕ 0 } and functions u(•) ∈ L ∞ (T ; U ) and v(•) ∈ L ∞ (T ; V ) be fixed. A solution of system (1), x(•; t 0 , x 0 , u(•), v(•)) = {ψ(•; t 0 , ψ 0 , u(•), v(•)), ϕ(•; t 0 , ϕ 0 , u(•), v(•))}, is a unique function x(•) = x(•; t 0 , x 0 , u(•), v(•)) ∈ V (1) T = V 1 × V 1 , V 1 = W 2,1
2 (Q) satisfying relations [START_REF] Fattorini | Infinite Dimensional Optimization and Control Theory[END_REF]. By virtue of the corresponding embedding theorem, without loss of generality, one can assume that the space V

(1) T is embedded into the space C(T ; H × H). Therefore, the element x(t) = {ψ(t), ϕ(t)} ∈ H × H is the phase state of system (1) at the time t. The following theorem takes place Theorem 1. [6, p. 25, Assertion 5] Let conditions (A1)-(A3) be fulfilled. Then for any u(•) ∈ L ∞ (T ; U ) and v(•) ∈ L ∞ (T ; V ) there exists a unique solution of system (1).

The paper is devoted to the investigation of the game control problem. Let us give the informal formulation of this problem. Let a uniform net ∆ = {τ i } m i=0 , τ i = τ i-1 + δ, τ 0 = t 0 , τ m = ϑ with a diameter δ = τ i -τ i-1 be fixed on a given time interval T . Let a solution of system (1) be unknown. At the times τ i ∈ ∆, a part of the phase states x(τ i ) (namely φ(τ i )) is inaccurately measured. The results of measurements

ξ h i ∈ H, i ∈ [1 : m -1], satisfy the inequalities |ξ h i -φ(τ i )| H ≤ h. (2)
Here, h ∈ (0, 1) is a level of informational noise. Let the following quality criterion be given:

I(x(•), u(•)) = σ(x(ϑ)) + ϑ t0 χ(t, x(t), u(t)) dt, where σ : H × H → R
and χ : T × H × H × U → R are given functions satisfying the local Lipschitz conditions. Let also a prescribed value of the criterion, number I * , be fixed.

The control problem under consideration consists in the following. There are two players-antagonists controlling system (1) by means of u and v, respectively.

One of them is called a partner; another, an opponent. Let P ⊂ U and E ⊂ V be given convex bounded and closed sets. The problem undertaken by the partner is as follows. It is necessary to construct a law (a strategy) for forming the control u (with values from P ) by the feedback principle (on the base of measurements of ϕ(τ i )) in order to provide the prescribed value of the quality criterion for any (unknown) realization v = v(•).

To form the control u providing the solution of the problem, along with the information on the "part" of coordinates of the solution of system (1) (namely, on the values ξ h i satisfying inequalities (2)), it is necessary to obtain some additional information on the coordinate ψ(•), which is missing. To get such a piece of information during the control process, it is reasonable, following the approach developed in [START_REF] Maksimov | Dynamical Inverse Problems of Distributed Systems[END_REF][START_REF] Maksimov | On an Inverse Problem for Phase Field Equations[END_REF][START_REF] Maksimov | Dynamical State and Control Reconstruction for a Phase Field Model[END_REF], to introduce an auxiliary controlled system. This system is described by a parabolic equation (the form is specified below). The equation has an output w * (t), t ∈ T , and an input p h (t), t ∈ T . The input p h (•) is some new auxiliary control; it should be formed by the feedback principle in such a way that p h (•) "approximates" the unknown coordinate ψ(•) in the mean square metric. Thus, along with the block of forming the control in the real system (it is called an controller), we need to incorporate into the control contour one more block (it is called an identifier) allowing to reconstruct the missing coordinate ψ(•) in real time. Note that, in essence, the identifier block solves a dynamical inverse problem, namely, the problem of (approximate) reconstruction of the unknown coordinate ψ(•). In he recent time, the theory of inverse problems for distributed systems has been intensively developed. Among the latest investigations, it is possible to mark out the research [START_REF] Liu | Inverse problem for a linearized Jordan-Moore-Gibson-Thompson equation[END_REF].

P T (•) = {u(•) ∈ L 2 (T ; U ) : u(t) ∈ P a.e. t ∈ T }, E T (•) = {v(•) ∈ L 2 (T ; V ) : v(t) ∈ E a.e. t ∈ T }. Elements of the product T × H are called positions, H = H × H × R × H × H × R.
Any function (perhaps, multifunction) U : T × H → P is said to be a positional strategy of the partner. The positional strategy corrects the controls at discrete times given by some partition of the interval T . Any function V : T × H × H → H is said to be a strategy of reconstruction. The strategy V is formed in order to reconstruct the unknown component ψ(•).

Consider the following ordinary differential equation

q(t) = χ(t, x(t), u(t)), q(t 0 ) = 0. (3)
Introducing this new variable q, we reduce the control problem of Bolza type to a control problem with a terminal quality criterion of the form I = σ(x(ϑ)) + q(ϑ). In this case, the controlled system consists of phase field equation ( 1) and ordinary differential equation (3). The scheme of an algorithm for solving the problem undertaken by the partner is as follows. In the beginning, auxiliary systems M 1 and M 2 (models) are introduced. The system M 1 has an input u * (•) and an output w(•); the system M 2 , an input p h (•) and an output w * (•), respectively. The model M 2 with its control law V forms the identifier, whereas the model M 1 and system (1) (with their control laws) form the controller. The process of synchronous feedback control of systems ( 1), (3), M 1 , and M 2 is organized on the interval T . This process is decomposed into m -1 identical steps. At the ith step carried out during the time interval [τ i , τ i+1 ), the following actions are fulfilled. First, at the time τ i , according to some chosen rules V and U, the elements

p h i ∈ V(τ i , ξ h i , w * (τ i )), u h i ∈ U(τ i , ξ h i , p h i , ψ h i , w(τ i )) are calculated.
Here, ψ h i is the result of measuring q(τ i ). Then (till the moment τ i+1 ), the control p h (t) = p h i , τ i ≤ t < τ i+1 , is fed onto the input of the system M 2 ; the control u h (t) = u h i , τ i ≤ t < τ i+1 , onto the input of system (1), [START_REF] Tröltzsch | Optimal Control of Partial Differential Equations. Theory, Methods and Applications[END_REF]. Under the action of these controls, as well as of the given control u * (t), τ i ≤ t < τ i+1 , and the unknown control of the opponent v(t), τ i ≤ t < τ i+1 , the states x(τ i+1 ), q(τ i ), w(τ i+1 ), and w * (τ i+1 ) are realized at the time τ i+1 . The procedure stops at the time ϑ.

Let models M 1 and M 2 with phase trajectories w(•) and w * (•) be fixed. A solution x(•) of system (1) starting from an initial state (t * , x * ) and corresponding to piecewise constant controls u h (•) and p h (•) (formed by the feedback principle) and to a control

v t * ,ϑ (•) ∈ E t * ,ϑ (•) is called an (h, ∆)-motion x h ∆,w (•) = x h ∆,w (•; t * , x * , U, V, v t * ,ϑ (•)
). This motion is generated by the positional strategies U and V. Thus, the motions x h ∆,w (•), q h ∆ (•),w(•), and w * (•) are formed simultaneously. So, for t ∈ [τ i , τ i+1 ), we define

x h ∆,w (t) = x(t; τ i , x h ∆,w (τ i ), u h τi,τi+1 (•), v τi,τi+1 (•)), q h ∆ (t) = q(t; τ i , q h ∆ (τ i ), u h τi,τi+1 (•)), w(t) = w(t; τ i , w(τ i ), u * τi,τi+1 (•)), w * (t) = w * (t; τ i , w * (τ i ), p h τi,τi+1 (•))
, where

u h (t) = u h i ∈ U(τ i , ξ h i , p h i , q h (τ i ), w(τ i )), p h (t) = p h i ∈ V(τ i , ξ h i , w * (τ i )) (4) for t ∈ [τ i , τ i+1 ), i ∈ [i(t * ) : m -1], |ξ h i -φ h ∆,w (τ i )| H ≤ h, |ψ h i -g h ∆ (τ i )| ≤ h, (5)
i(t * ) = min{i : τ i > t * }, u h (t) = u h * ∈ P, p h (t) = p h * ∈ H for t ∈ [t * , τ i(t * )
). The set of all (h, ∆)-motions is denoted by X h (t * , x * , U, V, ∆, w). Problem. It is necessary to construct a positional strategy U : T ×H → P of the partner and a positional strategy V : T × H × H → H of reconstruction with the following properties: whatever a value ε > 0 and a disturbance v T (•) ∈ E T (•), one can find (explicitly) numbers h * > 0 and δ * > 0 such that the inequalities

|I(x h ∆,w (•), u h T (•)) -I * | ≤ ε ∀x h ∆,w (•) ∈ X h (t 0 , x 0 , U, V, ∆, w) (6)
are fulfilled uniformly with respect to all measurements ξ h i and ψ h i with properties (5), if h ≤ h * and δ = δ(∆) ≤ δ * .

Algorithm for solving the Problem

To solve the Problem, we use ideas from [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF], namely, the method of a priori stable sets. In our case, this method consists in the following. Let a trajectory of model M 1 , w(•), possessing the property σ(w 1 (ϑ)) + w 2 (ϑ) = I * be known. Then, a feedback strategy U providing tracking the prescribed trajectory of M 1 by the trajectory of real system (1) is constructed. This means that the (h, ∆)motion x h ∆ (•) formed by the feedback principle (see ( 4)) remains at a "small" neighborhood of the trajectory w(•) during the whole interval T . This property of the (h, ∆)-motion allows us to conclude that the chosen strategy solves the considered control problem.

Let us pass to the realization of this scheme. Define

Φ(t, x, u, v) = {Bu - Cv, χ(t, x, u)} ∈ H × R, Φ u (t, x, v) = u∈P Φ(t, x, u, v), H * (t; x) = v∈E Φ u (t, x, v), H * (•; x) = {u(•) ∈ L 2 (T ; H × R) : u(t) ∈ H * (t; x) f or a. a. t ∈ T }.
As a model M 1 , we take the system including two subsystems, i.e., the phase field equation

∂ ∂t w (1) + l ∂ ∂t w (2) = ∆ L w (1) + u 1 in Ω × (t 0 , ϑ], (7)
∂ ∂t w (2) = ∆ L w (2) + g(w (2) ) + w (1) with the boundary condition ∂ ∂n w (1) = ∂ ∂n w (2) = 0 on ∂Ω × (t 0 , ϑ] and the initial condition w (1) (t 0 ) = ψ 0 , w (2) (t 0 ) = ϕ 0 in Ω , as well as the ordinary differential equation

ẇ(3) (t) = u 2 (t), w (3) ∈ R, w (3) (0) = 0. ( 8
)
By the symbol w(•), we denote the solution of system ( 7), [START_REF] Maksimov | On an Inverse Problem for Phase Field Equations[END_REF]. Then, the model M 1 has the control u(•) = {u 1 (•), u 2 (•)}. As a model M 2 , we use the equation

∂ w * ∂t = ∆ L w * + p h + g(w * ) in Ω × (t 0 , ϑ] (9)
with the boundary condition ∂w * ∂n = 0 on ∂Ω ×(t 0 , ϑ] and the initial condition w * (t 0 ) = ϕ 0 in Ω. Condition 1. There exists a control u(

•) = u * (•) = {u * 1 (•), u * 2 (•)} ∈ H(t; w (1) (t), w (2) (t)) for a.a. t ∈ T such that I * = σ(w (1) (ϑ)) + w (2) (ϑ).
The strategies U and V (see ( 4)) are defined in such a way:

U(t, ξ, p, ψ, w) = arg min{L(u, y) + (ψ -w (3) )χ(t, ξ, p, u) : u ∈ P }, (10)
V(t, ξ, w * ) = arg min{l(t, α, u, s) : u ∈ U d }, (11)
where w = {w (1) , w (2) , w

(3) }, L(u, y) = (-y, Bu) H , y = w (1) -p + l(w (2) -ξ), l(t, α, u, s) = exp(-2b * t)(s, u) H + α|u| 2 H , s = w * -ξ, b * = |a + 1/3b| L∞(Q) , α = α(h), U d = {p(•) ∈ L 2 (T ; H) : |p(t)| H ≤ d for a.a. t ∈ T }, d ≥ sup t∈T {|ψ(t; t 0 , x 0 , u(•), v(•))| H : u(•) ∈ P T (•), v(•) ∈ E T (•)}.

Condition 2. Let h, δ(h), and α(h) satisfy the conditions :

α(h) → 0, δ(h) → 0, (h + δ(h))α -1 (h) → 0 as h → 0.
Let us pass to the description of the algorithm for solving the Problem. Namely, we describe the procedure of forming the (h, ∆)-motion x h ∆,w (•) = {ψ h ∆,w (•), ϕ h ∆,w (•)} and trajectory g h ∆ (•) corresponding to some fixed partition ∆ and the strategies U and V, see [START_REF] Liu | Inverse problem for a linearized Jordan-Moore-Gibson-Thompson equation[END_REF] and (11). Before the algorithm starts, we fix a value h ∈ (0, 1), a function α = α(h): (0, 1) → (0, 1), and a partition

∆ h = {τ h,i } m h i=0 with diameter δ = δ(h) = τ i+1 -τ i , τ i = τ h,i , τ h,0 = t 0 , τ h,m h = ϑ.
The work of the algorithm is decomposed into m h -1 identical steps. We assume that

u h (t) = u h 0 ∈ U(t 0 , ξ h 0 , p h 0 , 0, w(t 0 )), p h (t) = p h 0 ∈ V(t 0 , ξ h 0 , ϕ 0 ) (|ξ h 0 -ϕ 0 | H ≤ h) on the interval [t 0 , τ 1 )
. Under the action of these piecewiseconstant controls as well as of an unknown disturbance v t0,τ1 (•), the (h, ∆)-

motion {x h ∆,w (•)} t0,τ1 = {ψ h ∆,w (•; t 0 , ψ 0 , u h t0,τ1 (•), v t0,τ1 (•)), ϕ h ∆,w (•, t 0 , ϕ 0 , u h t0,τ1 (•), v t0,τt (•))} t0,τ1 of system (1), the trajectory {q h ∆ (•)} t0,τ = {q(•; t 0 , q(t 0 )), u h t0,τ1 (•))} t0,τ1 of equation (3), the trajectory {w * (•)} t0,τ1 = {w * (•; t 0 , φ 0 , p h t0,τ1 (•))} t0,τ1
of the model M 2 , and the trajectory {w(•)} t0,τ1 = {w(•; t 0 , x 0 , u * t0,τ1 (•))} t0,τ1 of the model M 1 are realized. At the time t = τ 1 , we determine u h 1 and p h 1 from the conditions

u h 1 ∈ U(τ 1 , ξ h 1 , p h 1 , ψ h 1 , w(τ 1 )), p h 1 ∈ V(τ 1 , ξ h 1 , w * (τ 1 )) (|ξ h 1 -ϕ h ∆,w (τ 1 )| H ≤ h, |ψ h 1 -g h ∆ (τ 1 )| ≤ h), i.e., we assume that u h (t) = u h 1 and p h (t) = p h 1 for t ∈ [τ 1 , τ 2 )
. Then, we calculate the realization of the (h, ∆)-motion

{x h ∆,w (•)} τ1,τ2 = {ψ h ∆,w (•; τ 1 , ψ h ∆,w (τ 1 ), u h τ1,τ2 (•), v τ1,τ2 (•)), ϕ h ∆,w (•; τ 1 , ϕ h ∆,w (τ 1 ), u h τ1,τ2 (•), v τ1,τ2 (•))} τ1,τ2 ,
the trajectory {q h ∆ (•)} τ1,τ2 = {q(•; τ 1 , q h ∆ (τ 1 ), p h τ1,τ2 (•))} τ1,τ2 of equation ( 3), the trajectory {w * (•)} τ1,τ2 = {w * (•; τ 1 , w * (τ 1 ), p h τ1,τ2 (•))} τ1,τ2 of the model M 2 , and the trajectory {w(•)} τ1,τ2 = {w(•; τ 1 , w(τ 1 ), u * τ1,τ2 (•)} τ1,τ2 of the model M 1 . Let the (h, ∆)-motion x h ∆,w (•), the trajectory q h ∆ (•) of equation ( 3), the trajectory w * (•) of the model M 2 , and the trajectory w(•) of the model M 1 be defined on the interval [t 0 , τ i ]. At the time t = τ i , we assume that

u h i ∈ U(τ i , ξ h i , p h i , ψ h i , w(τ i )), p h i ∈ V(τ i , ξ h i , w * (τ i )) (12) (|ξ h i -ϕ h ∆,w (τ i )| H ≤ h, |ψ h i -g h ∆ (τ i )| ≤ h), i.e., we set u h (t) = u h i and p h (t) = p h i for t ∈ [τ i , τ i+1 ).
As a result of the action of these controls and an unknown disturbance v τi,τi+1 (•), the (h, ∆)-motion

{x h ∆,w (•)} τi,τi+1 = {ψ h ∆,w (•; τ i , ψ h ∆,w (τ i ), u h τi,τi+1 (•), v τi,τi+1 (•)), ϕ(•; τ i , ϕ h ∆,w (τ i ), u h τi,τi+1 (•), v τi,τi+1 (•))} τi,τi+1 , the trajectory {q h ∆ (•)} τi,τi+1 = {q(•; τ i , q h ∆ (τ i ), p h τi,τi+1 (•)} τi,τi+1 of equation (3), the trajectory {w * (•)} τi,τi+1 = {w * (•; τ i , w * (τ i ), p h τi,τ1+1 (•))} τi,τi+1
of the model M 2 , and the trajectory {w(•)} τi,τi+1 = {w(•; τ i , w(τ i ), u * τi,τi+1 (•))} τi,τi+1 of the model M 1 are realized on the interval [τ i , τ i+1 ]. This procedure stops at the time ϑ.

Theorem 2. Let Conditions 1 and 2 be fulfilled. Let also the models M 1 and M 2 be specified by relations [START_REF] Maksimov | Dynamical Inverse Problems of Distributed Systems[END_REF], [START_REF] Maksimov | On an Inverse Problem for Phase Field Equations[END_REF], and (9), respectively. Then, the strategies U and V of form [START_REF] Liu | Inverse problem for a linearized Jordan-Moore-Gibson-Thompson equation[END_REF] and (11) solve the Problem. Proof. To prove the theorem, we estimate the variation of the functional

Λ(t, x h ∆,w (•), q h ∆ (•), w(•)) = Λ 0 (x h ∆,w (t), q h ∆ (t), w(t)) + +0.5 t 0 Ω |∇π h ( , η)| 2 dη + l 2 Ω |∇ν h ( , η)| 2 dη d , where x h ∆,w (•) = {ψ h ∆,w (•), φ h ∆,w (•)}, w(•) = {w (1) (•), w (2) (•), w (3) (•)}, π h (t) = w (1) (t)-ψ h ∆,w (t), ν h (t) = w (2) (t)-ϕ h ∆,w (t), g h (t) = π h (t)+lν h (t), Λ 1 (x h ∆,w (t), w(t)) = 0.5|g h (t)| 2 H + 0.5l 2 |ν h (t)| 2 H , λ(q h ∆ (t), w (3) (t)) = 0.5|q h ∆ (t) -w (3) (t)| 2 , Λ 0 (x h ∆,w (t), q h ∆ (t), w(t)) = Λ 1 (x h ∆,w (t), w(t)) + λ(q h ∆ (t), w (3)
(t). It is easily seen that the functions π h (•) and ν h (•) are solutions of the system

∂π h (t, η) ∂t + l ∂ν h (t, η) ∂t = ∆ L π h (t, η) + u * 1 (t, η) -(Bu h )(t, η) + (Cv)(t, η), (13)
∂ν h (t, η) ∂t = ∆ L ν h (t, η) + R h (t, η)ν h (t, η) + π h (t, η) in Ω × (t 0 , ϑ],
with the initial condition π h (t 0 ) = ν h (t 0 ) = 0 in Ω and with the boundary condi-

tion ∂π h ∂n = ∂ν h ∂n = 0 on ∂Ω × (t 0 , ϑ]. Here, R h (t, η) = a(t, η)+b(t, η)(w (1) (t, η)+ ϕ h ∆,w (t, η)) -((w (1) (t, η)) 2 + w (1) (t, η)ϕ h ∆,w (t, η) + (ϕ h ∆,w ) 2 (t, η)).
Multiplying scalarly the first equation of (13) by g h (t), and the second one, by ν h (t), we obtain

(g h (t), g h t (t)) H + Ω {|∇π h (t, η)| 2 + l∇π h (t, η)∇ν h (t, η)} dη = = (g h (t), u * 1 (t) -Bu h (t) + Cv(t)) H , (14)
(ν h (t), ν h t (t)) H + Ω |∇ν h (t, η)| 2 dη ≤ (π h (t), ν h (t)) H + b|ν h (t)| 2 H for a.a. t ∈ T.
Here, we use the inequality vrai max

(t,η)∈T ×Ω {a(t, η)+b(t, η)(v 1 +v 2 )-(v 2 1 +v 1 v 2 +v 2 2 )} ≤ b, which is valid for any v 1 , v 2 ∈ R. It is evident that the inequality Ω l(∇π h (t, η), ∇ν h (t, η)) dη ≥ -0.5 Ω {|∇π h (t, η)| 2 + l 2 |∇ν h (t, η)| 2 }
≤ 2l 2 λ 2 Λ 0 (x h ∆,w (t), q h ∆ (t), w(t)) + γ (1) t + γ (2)
t , where γ

(1) t = (q h ∆ (t) -w (3) (t))(χ(t, x h ∆,w (t), u h (t)) -u * 2 (t)), γ (2)
t = (g h (t), u * 1 (t) - Bu h (t) + Cv(t)) H for a.a. t ∈ T . For a.a. t ∈ δ i = [τ i , τ i+1 ], it is easily seen that γ (1) t ≤ (q h ∆ (τ i ) -w (3) (τ i ))(χ(t, x h ∆,w (t), u h (t) -u * 2 (t)) + k 0 (t -τ i ) 1/2 , (19)
|χ(t, x h ∆,w (t), u h (t)) -χ(τ i , ξ h i , p h i , u h (t))| ≤ k 1 {h + (t -τ i ) 1/2 + |p h i -ψ h ∆,w (t)| 2 H }.
(20) From ( 19) and (20) we have for a.a.

t ∈ δ i γ (1) t ≤ (q h ∆ (τ i ) -w (3) (τ i ))(χ(τ i , ξ h i , ψ h i , u h (t)) -u * 2 (t))+ (21) +k 2 |q h ∆ (τ i ) -w (3) (τ i )|{h + (t -τ i ) 1/2 + |p h i -ψ h ∆,w (t)| H }.
Estimate the last term in the right-hand side of inequality (18). For a.a.

t ∈ [τ i , τ i+1 [ |g h (t) -y h i | H = |π h (t) + lν h (t) -y h i | H ≤ λ 1,i (t) + λ 2,i (t), (22)
where

y h i = w (1) (τ i )-p h i -l(w (2) (τ i )-ξ h i ), λ 1,i (t) = |w (1) (t)-ψ h ∆,w (t)-w (1) (τ i )+ p h i | H , λ 2,i (t) = l|w (2) (t) -ϕ h ∆,w (t) -w (2) (τ i ) + ξ h i | H . For a.a. t ∈ δ i = [τ i , τ i+1 ], it is easily seen that λ 1,i (t) ≤ |p h i -ψ h ∆,w (t)| H + t τi | ẇ(1) (τ )| H dτ , λ 2,i (t) ≤ lh + l t τi {| φh ∆,w (τ )| H + | ẇ(2) (τ )| H } dτ .
From this equation and ( 22), for a.a. t ∈ δ i , it follows that

|g h (t)-y h i | H ≤ lh+ t τi {l| φh ∆,w (τ )| H +| ẇ(1) (τ )| H +| ẇ(2) (τ )| H } dτ +|p h i -ψ h ∆,w (t)| H . (23)
By virtue of (11), taking into account results of [START_REF] Maksimov | On an Inverse Problem for Phase Field Equations[END_REF][START_REF] Maksimov | Dynamical State and Control Reconstruction for a Phase Field Model[END_REF], from estimate (23) we derive

m-1 i=0 τi+1 τi M (t; τ i ) dt ≤ k 3 (h + δ) + k 4 ϑ t0 |p h (τ ) -ψ h ∆,w (τ )| H dτ ≤ k 5 ν 1/2 (h), (24)
where M (t;

τ i ) = |g h (t) -y h i | H {|Bu h i | H + |Cv(t)| H + |u * 1 (t)| H } for a.a. t ∈ δ i , ν(h) = (h + δ(h) + α(h)) 1/2 + (h + δ(h))α -1 (h). Note that, for a.a. t ∈ δ i , γ (2) t ≤ (y h i , u * 1 (t) -Bu h i + Cv(t)) H + M (t; τ i ). ( 25
)
Let the symbol (•, •) H×R denote the scalar product in the space H × R. Let us define elements v e i from the conditions inf u∈P

(s i , Φ(τ i , p h i , ξ h i , u, v e i )) H×R ≥ sup v∈E inf u∈P (s i , Φ(τ i , ξ h i , u, v)) H×R -h, (26)
where

Φ(τ i , p h i , ξ h i , u h i , v(t)) = {Bu h i -Cv(t), χ(τ i , ξ h i , p h i , u h i )}, s i = {-y h i , ψ h i - w (3) (τ i )}. It is obvious (see Condition 1) that u * (t) ∈ H(t, w (1) (t), w (2) (t)) ⊂ u∈P Φ(t, w (1) (t), w (2) (t), u, v e i ) for a. a. t ∈ [τ i , τ i+1 ).
Then, for a.a. t ∈ δ i , there exists a control u (1) (t) ∈ P such that

Φ(t, w (1) (t), w (2) (t), u (1) (t), v e i ) = u * (t) for a. a. t ∈ [τ i , τ i+1 ]. (27)
Using the rule of definition of the strategy U, we deduce that The statement of the theorem follows from the last inequality.

(

2

 2 Problem statement Before passing to the problem formulation, we give some definitions. Furthermore, we denote by u a,b (•) the function u(t), t ∈ [a, b], considered as a whole. The symbol P a,b (•) stands for the restriction of the set P T (•) onto the segment [a, b] ⊂ T . Any strongly measurable functions u(•) : T → P and v(•) : T → E are called program controls of the partner and opponent, respectively. The sets of all program controls of the partner and opponent are denoted by the symbols P T (•) and E T (•) :

  dη (15) is fulfilled for a.a. t ∈ T . Let us multiply the first inequality of (14) by l 2 and add to the second one. Taking into account (15), we have for a.a.t ∈ T Bu h (t) + Cv(t)) H + l 2 (π h (t), ν h (t)) H + bl 2 |ν h (t)| 2 H . Note that π h (t) = g h (t) -lν h (t). In this case, for a.a. t ∈ T (π h (t), ν h (t)) H + b|ν h (t)| 2 H = (g h (t) -lν h (t), ν h (t)) H + b|ν h (t)| 2 (t, η)| 2 + l 2 |∇ν h (t, η)| 2 } dη ≤ (18)

	(g h (t), g h t (t)) H +l 2 (ν h (t), ν h t (t)) H +0.5	{|∇π h (t, η)| 2 +l 2 |∇ν h (t, η)| 2 } dη (16)
		Ω
		≤ (g h (t), u * 1 (t) -H =	(17)
	= (g h (t), ν h (t)) H + (b -l)|ν h (t)| 2 H ≤ 0.5(|g h (t)| 2 H + (0.5 + |b -l|)|ν h (t)| 2 H .
	Combining (16) and (17), we obtain
	d dt	Λ 0 (x h ∆,w (t), q h ∆ (t), w(t)) + 0.5

Ω

{|∇π h

  s i , Φ(τ i , p h i , ξ h i , u h i , v(t))) H×R ≤ inf , Φ(τ i , p h i , ξ h i , u, v e i )) H×R + 2h ≤ (31) ≤ (s i , Φ(t, p h i , ξ h i , u (1) (t), v e i )) H×R + 2h + L|ψ h i -w (3) (τ i )|(t -τ i ).Here, L is the Lipschitz constant of the function χ(•) in the corresponding domain. In this case, for a.a. t ∈ δ i , it follows from (27), (31) that(s i , Φ(τ i , p h i , ξ h i , u h i , v(t)) -u * (t)) H×R ≤ (32) ≤ 2h + L|ψ h i -w (3) (τ i )|{t -τ i + |ξ h i -w (1) (t)| H + |p h i -w (2) (t)| H }.By virtue of (32) and (21), for the interval [τ i , τ i+1 ], we derive (•),q h ∆ (•), w(•)) + k 7 {h 2 + t -τ i + |p h i -ψ h ∆,w (t)| 2 H } for a.a. t ∈ δ i . Using (34) and (24), by virtue of the Gronwall lemma, we obtain for t ∈ T Λ(t, x h ∆,w (•), q h ∆,w (•), w(•)) ≤ k 8

					γ	(1) t + γ	(2) t	≤	(33)
	≤ π * i (t) + k 6 {h 2 + t -τ i + Λ 0 (x h ∆,w (t), q h ∆ (t), w(t)) + |p h i -ψ h ∆,w (t)| 2 H },
	where
					π * i (t) = (s i , Φ(τ i , p h i , ξ h i , u h i , v(t)) -u * (t)) H×R .
	We deduce from inequalities (18) and (33) that
					dΛ(t, x h ∆,w (•), q h ∆ (•), w(•)) dt	≤ M (t; τ i )+	(34)
	+k 6 Λ(t, x h ∆,w m-1	τi+1
						M (t; τ
					i=0	τi
						u∈P	sup v∈E	(s i , Φ(τ i , p h i , ξ h i , u, v)) H×R + h. (28)
	Here (see (12)), u h i ∈ U(τ i , ξ h i , p h i , ψ h i , w(τ i )); v τi,τi+1 (•) is an unknown realization
	of the control of the opponent. In turn, from (26) we have
	sup v∈E	inf u∈P	(s i , Φ(τ i , p h i , ξ h i , u, v)) H×R ≤ inf u∈P	(s i , Φ(τ i , p h i , ξ h i , u, v e i )) H×R + h. (29)
	Moreover, it is evident that the equality
	inf u∈P	sup v∈E	(s i , Φ(τ i , p h i , ξ h i , u, v)) H×R = sup v∈E	inf u∈P	(s i , Φ(τ i , p h i , ξ h i , u, v)) H×R (30)
	is valid. From (28)-(30) we have
	(s i , Φ(τ i , p h i , ξ h i , u h i , v(t))) H×R ≤ inf

u∈P

(s i i ) dt + δ(h 2 + δ) ≤ k 9 ν 1/2 (h).