Nikolay Timofeev

Galina Timofeeva

Estimation of Loan Portfolio Risk on the Basis of Markov Chain Model

Keywords: Loan portfolio, Markov chain, incomplete information

A change of shares of credits portfolio is described by Markov chain with discrete time. A credit state is determined on as an accessory to some group of credits depending on presence of indebtedness and its terms. We use a model with discrete time and fix the system state through identical time intervals -once a month. It is obvious that the matrix of transitive probabilities is known incompletely. Various approaches to the matrix estimation are studied and methods of forecast the portfolio risk are proposed. The portfolio risk is set as a share of problematic loans. We propose a method to calculate necessary reserves on the base of the considered model.

Introduction

Markov chain models [START_REF] Markov | Rasprostranenie zakona bol'shih chisel na velichiny, zavisyaschie drug ot druga[END_REF] are widely used to explain the dynamics of state changes for different systems. Often they are used as a mathematical model for some random physical process.

Markov chains are used in finance and economics to model a variety of different phenomena, including asset prices, market crashes and credit portfolio dynamics [START_REF] Jones | Estimating Markov Transition Matrices Using Proportions Data: An Application to Credit Risk[END_REF][START_REF] Thyagarajan | Retail Banking Loan Portfolio Equilibrium Mix: A Markov Chain Model Analysis[END_REF].

If the transition probability matrix of Markov chain is known then dynamics of the system states probabilities is completely described by a system of difference equations. As a rule the transition probabilities are unknown and estimated during the system evolution.

The most important indicator of a bank loan portfolio quality is a probability of default which is closely connected with a share of the problematic loans [START_REF] Hanson | Confidence intervals for probabilities of default[END_REF][START_REF] Andrade | Structural models in consumer credit[END_REF]. A value of necessary reserves depends on quality and structure of the portfolio. On the one hand reserves should provide low probability of default, on the other hand they impact on profitability of the portfolio.

Let's assume that a change of shares of credits portfolio is described by Markov chain with discrete time. In this case the credit state is determined on as an accessory to some group of credits depending on presence of indebtedness and its terms. We will use a model with discrete time and fix the system state through identical time intervals -once a month. It is proposed that the transition probabilities vary a little. Thus, we consider a stable economic situation when the transition probabilities are constant.

It is obvious that the matrix of transitive probabilities is known incompletely. Its values are estimated using a data on changes the quality of loans (a migration analysis of the portfolio). Criterion of a choice is accuracy of the forecast of a share of problematic loans.

Mathematical Model

Dynamics of System States

We consider a system with k states, the probability that the system is in i-th state at moment t denote by x i (t), i = 1, . . . , k. Thus the following conditions hold:

0

≤ x i (t) ≤ 1, x 1 (t) + . . . + x k (t) = 1. ( 1
)
The dynamics of the system states probabilities is described by the discrete Markov chain model:

x j (t + 1) = k i=1 p ij x i (t), t = 0, 1, . . . , T, (2)
where p ij is the probability of transition from state i to state j in one step. The first-order stationary Markov model for credit transitions is somewhat restrictive as a credit quality responds to changes in economics. Using a higherorder Markov process or a nonstationary transition probability matrix may be more appropriate, but in such models one should estimates too many parameters and the requirements to statistical data increase quite substantially. Thus the simple Markov chain is usually used in stable economic situation and with not longer time horizon [START_REF] Jones | Estimating Markov Transition Matrices Using Proportions Data: An Application to Credit Risk[END_REF][START_REF] Thyagarajan | Retail Banking Loan Portfolio Equilibrium Mix: A Markov Chain Model Analysis[END_REF].

Let's denote by x(t) a vector of states probabilities x(t) = {x 1 (t), . . . , x k (t)} ⊤ , by P a matrix of the transition probabilities P = {p ij } and rewrite equation (2) in the vector form:

x(t + 1) = P ⊤ x(t), t = 0, 1, . . . , T.

(

When the transition probability matrix P is known incompletely there is a problem to estimate x(T ). It is assumed that we have information about the number of transitions from i-th state to j-th on t step, t = 1, . . . , m.

Ways to Select Groups

We consider two ways of describing a credit portfolio dynamics: a regular Markov chain in which we do not take into account repaid loans and renovation of the portfolio and a scheme which included "new loan" and "repaid loan" as possible states of a loan. In the first way we investigate a steady-state behavior of the loan portfolio shares, in the second way we study the profitability of loans from delivery to its repayment. Let us consider a fist way of describing the states of loan. For beginning we consider a simplified scheme with three groups of loans:

1. Loans without delay, including the new ones (S 1 ); 2. Overdue loans with 1 -65 days delay (S 2 ); 3. Non-performing (problematic) loans (S 3 ). A graph of the system state is in Fig. 1. An investigation of this scheme allows to define the most important features of the model. As a rule risk-managers use more detailed schemes which takes into account a number of overdue days on a loan. For example let consider the expanded scheme with 5 groups of loans [START_REF] Timofeev | Mathematical Model of the Vintage Analysis of Banking Credit Portfolio[END_REF]: loans without delay, including a new one; overdue loans with 1 to 35 days delay; overdue loans with 35 to 65 days delay; more than 65 days overdue loans (problematic loans); reconstructed loans.

In these schemes the renovation and repayment of credits are considered without allocating a separate state. New credits are included into the first group together with credits without delay.

A loan is "reconstructed" in case it was problematic in the previous period and a borrower has made partial payments under the credit.

In schemes with the repayments(amortization) there is the category "repaid credits". The scheme on the basis of the expanded scheme of loans has 6 groups of loans: loans without delay, including a new one; overdue loans with 1 to 35 days delay; overdue loans with 35 to 65 days delay; more than 65 days overdue loans (problematic loans); reconstructed loans; repaid loans.

Estimation of Transition Probabilities

For the estimation of the probability p ij one usually use the statistical data about transitions from one state to another. Let't denote by n i (t) the number of individuals who are in state i in period t, and by n ij (t) the number of individuals who were in state i in period t -1 and are in state j in period t. We can estimate the probability p ij of an individual being in state j in period t given that they were in state i in period t -1.

The probability of transition p ij (t) from any given state i is approximated by a proportion of individuals that started in state i and ended in state j as a proportion of all individuals in that started in state i:

w ij (t) = n ij (t) n i (t -1) . ( 4
)
If the Markov chain is stationary (i.e. p ij (t) ≡ p ij ) then one can use another estimate:

wij = T t=1 n ij (t) T -1 t=0 n i (t) (5)
Using the methods described above, it is possible to estimate a transition matrix using count data.

Anderson and Goodman [START_REF] Anderson | Statistical Inference About Markov Chains[END_REF] showed that the estimator w ij given by equation ( 5) is a maximum-likelihood estimator and calculated statistical moments of random values

ξ ij (t) = n ij -p ij n i (t -1).
On the t-th step n i (t -1) n i , i = 1, . . . , k are known. The statistical moment of w ij are following [START_REF] Timofeeva | Predicting the Components of Credit Portfolio Based on a Markov Chain Model[END_REF]:

E(w ij (t)) = p ij , (6)
V ar(w ij (t)) = pij (1-pij ) ni , Cov(w ij (t), w il (t)) = - pij p il ni , j = l, Cov(w ij (t), w cl (t)) = 0, i = c. (7)
Here E(ξ) is a mathematical expectation of a random value ξ, V ar(ξ) is its variance, Cov(ξ, ζ) is a covariance between ξ and ζ.

From relations [START_REF] Timofeev | Mathematical Model of the Vintage Analysis of Banking Credit Portfolio[END_REF] follows that w ij and w cl are noncorrelated if i = c. Suppose that instead of observing the actual count of transitions from the different states, we only observe the aggregate proportions y i (t), which represent the proportion of observations with the state i. The aggregate proportions y i (t) estimate the system state probabilities:

y i (t) ≈ N x i (t), i = 1, . . . , N , t = 1, . . . , T ,
where N is a number of all individuals.

If the time series of observations T are sufficiently long, it is possible to estimate a transition matrix P from aggregate data using the least square method [START_REF] Kalbfleisch | Least-Squares Estimation of Transition Probabilities From Aggregate Data[END_REF] and its generalizations [START_REF] Macrae | Estimation of Time-Varying Markov Processes with Aggregate Data[END_REF].

The maximum-likelihood estimates in asset prices model is used to estimate transition matrices in credit risk modeling with a decades-old methodology that uses aggregate proportions data [START_REF] Jones | Estimating Markov Transition Matrices Using Proportions Data: An Application to Credit Risk[END_REF].

To specificate an order and a number of states in Markov chain model one may use criteria [START_REF] Anderson | Statistical Inference About Markov Chains[END_REF], but they are based on the detailed analysis of the data such as n ijm (t), where n ijm (t) is a number of individuals in state i at t -2, in j at t -1 and in m at t.

Approaches to Estimation the Share of Problematic Loans

We considered two ways to forecast the share of problem loans taking into account the uncertainty of the transition probabilities matrix. There are a confidence estimation method and a simulation method.

Confidence Estimation

The confidence estimation method consists of two stages: construction a confidence set Z α for the transition probabilities matrix based on statistical data and relations ( 6)-( 7) and the analysis of all possible trajectories of the system taking into account that the probabilities are constant but uncertain. Let's estimate the elements of the transition probability matrix P . Denote the confidence region for

{p ij = z s , i = 1, . . . , k, j = 1, . . . , k, j = i} on m-th step by Z α ⊂ R K , s = 1, . . . , K, K = k(k -1).
Thus p ij are the transition probabilities then Z α ⊂ Z + ⊂ R K , where Z + is the set of all possible values of transition probabilities {p ij , j = i}

Z + = {z s : 0 ≤ z s ≤ 1, K s=1 z s ≤ 1} ⊂ R K .
Estimation for p ii follows from the equalities

p i1 + . . . + p ik = 1, i = 1, . . . , k. (8)
In the considered model p ij = z s are distributed approximately normal [START_REF] Anderson | Statistical Inference About Markov Chains[END_REF] with mean values equal to w ij zs and a covariance matrix G defined by relations [START_REF] Timofeev | Mathematical Model of the Vintage Analysis of Banking Credit Portfolio[END_REF] with substitution p ij by w ij .

Therefore we may use the confidence set Z α defined by joint restrictions:

Z α = {z ∈ Z + : (z -z) ⊤ G(z -z) ≤ b (α) K }, (9)
where b

(α)
K is the α-quantile of χ 2 distribution with K degrees of freedom. The next step is to solve the state estimation problem of a multistage deterministic system with uncertain matrix P :

x(t + 1) = P ⊤ x(t), t = m, . . . , T, x(m) = x * , P ⊂ Z. ( 10
)
and to find an information set

X(t, Z) = {x ∈ R k + : x = (P ⊤ ) T -m x * , P ∈ Z}. ( 11
)
We may construct information sets for system (10) using approaches proposed by Kurzanski and Tanaka [11]. This method is very time-consuming because the number of estimated elements of the matrix is large. We can perform calculations only for the scheme with 3 groups of loans, in which only 4 probabilities should be estimated. For schemes 2 and 3 the number of estimated probabilities are 9 and 17 respectively.

Example 1. Let us consider the estimation problem for the scheme with 3 groups of loans (Fig. 1). For this scheme the matrix of transition probabilities has a form

P =   1 -p 12 p 12 0 p 21 1 -p 21 -p 23 p 23 0 p 32 1 -p 32   , (12)
and only 4 transition probabilities p ij should be estimated. Denote them as z s , s = 1, . . . , 4:

p 12 = z 1 , p 23 = z 2 , p 32 = z 3 , p 21 = z 4 . ( 13
)
From a statistical data we estimate the mean values zs = w ij , where w ij defined by [START_REF] Anderson | Statistical Inference About Markov Chains[END_REF]. A covariance matrix is calculated using relations [START_REF] Timofeev | Mathematical Model of the Vintage Analysis of Banking Credit Portfolio[END_REF], where estimates w ij substituted instead of values p ij :

G =      z1(1-z1) n1 0 0 0 0 z2(1-z2) n2 0 -z2 z4 n2 0 0 z3(1-z3) n3 0 0 -z2 z4 n2 0 z4(1-z4) n2      . ( 14
)
Then find a confidence set for transition probabilities p ij in ellipsoidal form [START_REF] Kalbfleisch | Least-Squares Estimation of Transition Probabilities From Aggregate Data[END_REF]. For the considered scheme with 3 groups of loans we obtain

Z α = {z s ∈ R 4 + : (z -z) ⊤ G(z -z) ≤ b (α) 4 }. (15)
Then we find an information set for system (5) for a given T = 12 using ellipsoidal calculus [START_REF] Kurzhanskiy | Ellipsoidal Techniques for Reachability Analysis of Discrete-Time Linear Systems[END_REF]. Thus we get the confidence set for the portfolio shares after 12 months x( 12) and for the share of problematic loans x 3 [START_REF] Kurzhanskiy | Ellipsoidal Techniques for Reachability Analysis of Discrete-Time Linear Systems[END_REF] in particularly. For our data we obtain x 3 (12) ∈ [0.03; 0.152] with probability α = 0.95.

Simulation Method

This method has 4 stages:

1. We determine statistical moments (4) and ( 5) for the transition probabilities based on statistical data on transitions between loan states;

2. Generate the random vector of unknown probabilities as the Gaussian vector with the given statistical moments; 3. Simulate the dynamics of system (2) with the large number of runs; 4. On the basis of the calculations determine the mean value and the confidence interval for the share of problem loans.

A histogram of the predicted share of problem loans after 12 months one can see in Fig. 2. On the base of modeling we find the expected mean of value of the share of problem loans after 12 month: x(t + 12) = 0.09 and quantile q 0.95 = 0.135, i.e. x(t + 12) < 0.135 with the probability α = 0.95. This result is more precise than the confidence interval obtained by the confidence estimation, because the confidence interval for the share of problematic loans, obtained in section 4.1, depends on the form of confidence set for unknown parameters p ij and the ellipsoidal form is not optimal in this case.

Estimation of a reserve

According recommendations of the International Committee [START_REF]Basel Committee on Banking Supervision, International Convergence of Capital Measurement and Capital Standards: A Revised Framework[END_REF] bank portfolio managers should estimate a risk of the portfolio and form reserves in a proportion of the its value.

There are different approaches to estimate the necessary reserves. Some of them based on probability of "nonreconstruction" of the problematic credits. Let's denote by p

[t] ij a probability of transition from the i-th state to j-th state by t steps. It is an element of the t power of the transition probability matrix

P [t] = P t . Denote by D [t] = p [t]
mm , where m-th state is "problematic loans". The value 1 -D [T ] is called a probability of reconstruction during a period T .

The value of risk for j-th group of loans is defined as

rj = p [τ1] jm D [τ2] , if j = m, rm = D [τ1] .
The reserve W is equals to: W = W 1 r1 + . . . + W k rk , where W j is a sum of loans in j-th group. Risk managers take different τ 1 and τ 2 , such as 6, 12 or 24 months. The problem is how to choose the period of transition to the problematic state and the reconstruction period. In some approaches value of τ 1 depends on j and a period of an overdue.

We propose another approach taking into account time structure of possible losses. Let's define a risk of loans in j-th group as maximum of a sum of problematic loans in future for loans of this group. Thus, the risk for j-th loans group equals

r j = max t∈0,...,T 1 (1 + ρ) t p [t] jm , (16)
where ρ is a month discount factor.

For a new loan we have

r 0 = max t∈0,...,T 1 (1 + ρ) t p [t] 0m . (17)
We can take into account that transition probabilities are known incompletely and instead of (16) use its quantile:

q j (α) : P{ max t∈0,...,T 1 (1 + ρ) t p [t] jm ≤ q j (α)} ≥ α, (18)
where P(A) is a probability of a random event A. The quantile may be calculated using the confidence approach or the simulation method (see Sect. 4).

Example 2. Let's consider a scheme with repayment (Scheme 3) and calculate the value of reserves for new loans using the proposed approach. Is proposed that we may estimate statistical moments (6) of the transition probabilities p ij based on a previous data and estimates (4), [START_REF] Hanson | Confidence intervals for probabilities of default[END_REF] or their modifications. We take a possible value of transitions probabilities matrix

P = {p ij }, cal- culate h(t) = 1 (1 + ρ) t p [t]
0m , t = 0, . . . , T, and find their maximum r 0 (t) which depended on matrix P . In Fig. 3 one can see changes of h(t) (the discounted sum of problem loans) for a fixed matrix P . Values of h(t) decreasing for t ≥ 20.

We generate many times transition probabilities according normal distribution with the given statistical moments, calculated the discounted sum of problem loans h(t) = h(t, P ), and its maximum r 0 = r 0 (P ), then we obtain estimates of a mean value and quantile (18).

In Fig. 4 one can see a histogram for risk estimates r 0 = r 0 (P ). For considered data mean value of risk r0 = 0.039, q 0.95 = 0.065.

Conclusion

We studied the discrete Markov chain model for loan portfolio. It is proposed that transition probabilities are unknown and estimated during the process. We proposed methods to estimate system state probabilities in future. Obtained results apply to forecast credit portfolio shares and to define necessary reserves.

Fig. 1 .

 1 Fig. 1. Graph of the Markov chain for the simplified scheme.

Fig. 2 .

 2 Fig. 2. Histogram of the share of problematic loans and normal density.

Fig. 3 .

 3 Fig. 3. Dependence of the discounted sum of problem loans h(t) on time.

Fig. 4 .

 4 Fig. 4. Histogram of risk r0(P ) and normal density

⋆ The investigations are partially supported by the Russian Foundation for Basic Research, project 10-01-00672a.


Peter Benner
email: benner@mpi-magdeburg.mpg.de

Sabine Hein
email: sabine.hein@mathematik.tu-chemnitz.de

MPC/LQG for Infinite-Dimensional Systems Using Time-Invariant Linearizations

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

Introduction

The control of nonlinear processes is a fundamental problem in engineering. A usual strategy for computer-aided control consists in pre-computing, in an offline phase, an optimal trajectory and control input, and in the implementation (online phase) to endow the system with a feedback controller in order to compensate for external disturbances and deviations from the optimized trajectory. A successful strategy for designing a nonlinear control scheme for complex dynamical systems, whose global optimization is impossible in real time, is model predictive control (MPC), see, e.g., [START_REF] Camacho | Model Predictive Control[END_REF][START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF]. In this approach, the behavior of the dynamical process is predicted on a small (local) time horizon and then optimized for a certain time interval using an auxiliary problem for which the computational solution of the local optimization problem is feasible in real-time. The control strategy is then applied for a small time step, the process is advanced in time, and for the next time step, prediction and optimization are repeated for the new state of the system based on new available measurements. Under certain conditions, this process converges to the optimal solution of the global control problem if the time steps and prediction/optimization horizons tend to zero, see, e.g., [START_REF] Camacho | Model Predictive Control[END_REF][START_REF] Grüne | Nonlinear Model Predictive Control: Theory and Algorithms[END_REF] and references therein. If used as a feedback control scheme, the "optimization" goal is to minimize the deviation from the desired trajectory so that stabilization, i.e., convergence to zero, becomes the goal.

If the state is not fully available in the prediction step, one is faced with the problem of incomplete observations. This requires to include a state estimator in the prediction/optimization step. If the local optimization problem is solved via an auxiliary linear-quadratic optimal control (LQ) problem (based on a suitable linearization of the nonlinear system) without control and state constraints, the optimal state estimator in a least-squares sense is given by the Kalman(-Bucy) filter [START_REF] Kalman | New results in linear filtering and prediction theory[END_REF], and the solution of the local LQ problem is obtained by the linearquadratic Gaussian (LQG) design, consisting of a combination of the Kalman filter and the linear-quadratic regulator (LQR). This requires the solution of two Riccati equations. Depending on whether the linearization is time-invariant or time-varying, these are the algebraic or differential Riccati equations (ARE or DRE, respectively) -see, e.g., [START_REF] Datta | Numerical Methods for Linear Control Systems[END_REF] or any other textbook on linear control design. In [START_REF] Ito | Receding horizon control with incomplete observations[END_REF], this MPC/LQG approach was suggested for finite-dimensional problems. Here, we will extend this idea to infinite-dimensional systems.

In the following we consider the control problem min

T f 0 y(t), Q(t)y(t) Y + u(t), R(t)u(t) U dt + G(x(T f )), (1)
subject to the semi-linear stochastic system with additive unmodeled disturbance

ẋ(t) = f (x(t)) + Bu(t) + F v(t), t > 0, x(0) = x 0 + η, (2)
u(t) ∈ U, x(t) ∈ X , where v(t) is an unknown Gaussian disturbance process with covariance V and η denotes the noise in the initial condition. Since in many applications the state is not completely available we introduce the output function (simulating measurements)

y(t) = Cx(t) + w(t), y ∈ Y,
where w(t) is a measurement noise process which will also be assumed to be Gaussian with covariance W . If (2) is an ordinary differential equation then we have a finite-dimensional problem with X = R n , Y = R p and U = R m . In the case of a partial differential equation (PDE), the problem is infinite-dimensional and X , Y, U are appropriate Hilbert spaces. Here, B, F, Q, R are linear operators on these Hilbert spaces, f is a nonlinear map, and . , . are inner products on the respective Hilbert spaces. The outline of the paper is as follows: in the next section, we will briefly sketch the MPC/LQG control design. In Section 3, we will then provide the necessary theoretical background to solve the local LQG problems for infinite-dimensional systems and briefly discuss a framework for the numerical approximation of the solution of the AREs to be solved in a practical implementation of the infinite-dimensional controller. Note that we have demonstrated the efficiency of the suggested MPC/LQG approach for the stabilization of the noisy Burgers equation in [START_REF] Benner | MPC for the Burgers equation based on an LGQ design[END_REF] and for a 3D reaction-diffusion system in [START_REF] Benner | Model predictive control for nonlinear parabolic differential equations based on a linear quadratic Gaussian design[END_REF]. Concluding remarks are provided in Section 4.

The MPC/LQG controller

Given a reference trajectory and control (x(t), ū(t)) obtained, e.g., from an offline optimization procedure, in the following we design a model predictive feedback control strategy. This MPC/LQG approach is based on a linearization of (2) on small intervals to obtain a linear time-invariant (LTI) or time-varying (LTV) problem. Due to space restrictions, we will concentrate here on the LTI case. For the general strategy in the LTV case, we refer to [START_REF] Benner | Model predictive control based on an LQG design for time-varying linearizations[END_REF][START_REF] Hein | MPC/LQG-Based Optimal Control of Nonlinear Parabolic PDEs[END_REF] We solve this linear problem on a small interval by using an LQG design. Note that we write M * to denote the adjoint operator corresponding to the linear operator M and the derivative of f from ( 2) is to be understood as the Fréchet derivative. With these preliminaries, the strategy is the following:

(1) Prediction and optimization step on

[t i , t i + T p ], t i + T p ≤ T f : linearize (2) around a given set point x to obtain A = f ′ (x(t i ))
and the linear state equation

ż(t) = Az(t) + B ũ(t) + F v(t), z(t i ) = x(t i ) -x(t i ), y(t) = Cx(t) + w(t), with z(t) = x(t) -x(t) and ũ(t) = u(t) -ū(t). Then solve the ARE 0 = XA + A * X -XBR -1 B * X + C * QC (3)
in order to obtain

X ⋆ and K = -R -1 B * X ⋆ .
(

) Implementation step on [t i , t i + T c ], T c ≤ T p : solve the filter ARE (FARE) 0 = AΣ + ΣA * -ΣC * W -1 CΣ + F V F * , 2
where V , W are the covariance matrices of the noise processes. Feed the real system on [t i , t i + T c ] with

u(t) = ū(t) -K(x(t) -x(t)),
and obtain the "measurement" y(t) by solving the nonlinear system on [t i , t i + T c ]. Estimate the state by x(t) by solving

ż(t) = Aẑ(t) + B ũ(t) + L y(t) -C x(t) , ẑ(t) = x(t) -x * (t), using the estimator gain L = Σ ⋆ C * W -1 .
(3) Receding Horizon Step: update t i+1 = t i + T c and go to the first step.

Note that the solutions of the AREs (3) and ( 4) are linear selfadjoint operators on D(A), the domain of A, and D(A * ), respectively. Some remarks are in order:

Remark 1. For the finite-dimensional case, if G in (1) is selected as a control Lyapunov function, Ito and Kunisch established the asymptotic stability and estimated the performance for the receding horizon synthesis in [START_REF] Ito | On asymptotic properties of receding horizon optimal control[END_REF]. Analogous results for the LTV case in the slightly more general MPC setting are obtained in [START_REF] Benner | Model predictive control based on an LQG design for time-varying linearizations[END_REF][START_REF] Hein | MPC/LQG-Based Optimal Control of Nonlinear Parabolic PDEs[END_REF].

Remark 2. Solving the AREs (3) and ( 4) yields an LQG controller for an infinite time-horizon. Therefore, we can also consider this scheme as an MPC scheme with infinite prediction and optimization horizon.

In the following, we will discuss an appropriate setting in which this procedure is well-posed and can be approximated using an appropriate numerical scheme. Note that using efficient numerical algorithms, large-scale AREs resulting from discretized PDE control problems can be solved in a reasonable time-scale, see [START_REF] Benner | Numerical solution of large Lyapunov equations, Riccati equations, and linear-quadratic control problems[END_REF]. Whether or not this is real-time feasible depends on the control horizon T c of the process. Further advances in computer hardware and improvements of the numerical algorithms will certainly allow real-time solution of AREs for moderately fast processes with medium-fine granularity of the discretization in the near future.

Infinite-dimensional LQG theory

Consider the following nonlinear optimal control problem:

min J (u) := x T f , Gx T f X + ∞ 0 x(t), C * QCx(t) X + u(t), Ru(t) U dt, (5)
subject to ẋ(t) = f (x(t)) + Bu(t) + F v(t), t > 0,

y(t) = Cx(t) + w(t), t > 0, x(0) = x 0 + η.
Following the infinite-dimensional LQG theory derived in [START_REF] Curtain | Infinite Dimensional Linear System Theory[END_REF] and denoting the set of linear maps from M to N by L(M, N ), we will assume the following:

Assumption 1 -X , Y, U are Hilbert spaces, f : D(f ) ⊆ X → X is a nonlinear map; -B ∈ L(U , X ), F ∈ L(U , X ), C ∈ L(X , Y), G ∈ L(X ); -Q ∈ L(Y), R, R -1 ∈ L(U )
, all self-adjoint and nonnegative and ν, Rν ≥ α||ν|| 2 for all ν ∈ U and some α > 0; -x 0 ∈ X , η is a zero mean Gaussian random variable on X with covariance Σ 0 , -v(t) and w(t) are Wiener processes (Gaussian with zero mean) on the Hilbert spaces U and Y with incremental covariance operators V ∈ L(U ) and W, W -1 ∈ L(Y), respectively.

Assuming that f (x) is Fréchet-differentiable and linearizing on small intervals [t i , t i + T p ] around a stationary operating point x, we obtain the stochastic LTI problem in differential form on [t i , t i + T p ]:

dz(t) = Az(t)dt + B ũ(t)dt + F dv(t), t i < t < t i + T p , dỹ(t) = Cz(t)dt + dw(t), t i < t < t i + T p , (6)
z(t i ) = z ti , with z(t) := h(t) = x(t) -x(t), ũ(t) = u(t) -ū(t) and f (x + h)(t) ≈ f (x(t)) + Ah(t),
where A := f ′ (x(t i )) is the Fréchet derivative of f , evaluated at x(t i ).

Note that the linear system ( 6) is only a local approximation to the original nonlinear system, but nevertheless we will use it to solve the LQG problem on an infinite horizon. The so obtained control is then only applied locally on the control horizon [t i , t i + T c ], then the prediction horizon is shifted by T c , and a new linearization based on x(t i + T c ), leading to a new LQG problem, is used.

To avoid problems of existence and uniqueness of the stochastic evolution equation ( 6), we use its integral form on [t i , t i + T p ]:

z(t) = S t-ti z(t i ) + t ti S t-s B ũ(s) ds + t ti S t-s F dv(s), t i ≤ s ≤ t ≤ t i + T p , (7)
ỹ(t) = t ti Cz(s) ds + w(t), t i < t ≤ t i + T p , z(t i ) = z ti ,
where S t is a strongly continuous semigroup on X generated by A on [t i , t i + T p ] (see, e.g., [START_REF] Engel | One-Parameter Semigroups for Linear Evolution Equations[END_REF] for the notion of semigroups and their properties).

A direct consequence of results from [START_REF] Curtain | Infinite Dimensional Linear System Theory[END_REF] is then the following theorem which yields the solution to the MPC/LQG/LTI problem on [t i , t i + T p ] for T p = ∞: Theorem 1. Under the Assumptions 1., the optimal control and corresponding estimated state for the minimization problem (5) subject to [START_REF] Benner | Numerical solution of large Lyapunov equations, Riccati equations, and linear-quadratic control problems[END_REF]

on [t i , t i + T p ] are given by u * (t) = u r (t) -R -1 B * Π ∞ x * (t) -x(t) , x * (t) = S ti x(t i ) + t ti S t-s Σ ∞ C * W -1 dy(s) + t ti S t-s (f (x(s)) -Ax(s)) ds,
where S t is the strongly continuous semigroup generated by

A -BR -1 B * Π ∞ -Σ ∞ C * W -1 C,
and Π ∞ and Σ ∞ are the unique nonnegative, self-adjoint solutions of the ARE

0 = A * Π + ΠA -ΠBR -1 B * Π + C * QC,
and the FARE

0 = AΣ + ΣA * -ΣC * W -1 CΣ + F V F * .
Note again that the global solution of the LQG problem which is computed for T p = ∞ in the above theorem is only used locally on [t i , t i + T c ]. In this way, we can implement the infinite-dimensional MPC/LQG controller for any infinite dimensional system satisfying the assumptions made in this section.

Checking the Assumptions 1 for a practical problem is tedious. In [START_REF] Hein | MPC/LQG-Based Optimal Control of Nonlinear Parabolic PDEs[END_REF], as an example the control problem for the Burgers equation

x t (t, ξ) = ν x ξξ (t, ξ) -x(t, ξ) x ξ (t, ξ) + b(ξ)u(t), on (0, T f ] × (0, 1) x(t, 0) = x(t, 1) = 0, t ∈ (0, T f ], x(0, ξ) = x 0 (ξ), ξ ∈ (0, 1),
is considered and it is shown that linearization about a set point leads to an infinite-dimensional LTI system satisfying the Assumptions 1.

For a practical implementation, it is of course necessary to discretize the infinite-dimensional system and to work with a finite-dimensional approximation. If the infinite-dimensional problem is defined via a PDE, this can be achieved using a spatial semi-discretization based on finite differences or finite elements. As both A and its adjoint A * appear in the formulation of the problem and their discretizations A h and A T h are used to define the finite-dimensional AREs that need to be solved to obtain the approximate feedback and estimator gain matrices K h and L h , standard convergence results for finite element discretizations are not sufficient. The necessary conditions for dual convergence are stated in [START_REF] Banks | The linear regulator problem for parabolic systems[END_REF] for linear parabolic equations with distributed control and are generalized in [START_REF] Benner | Linear-quadratic regulator design for optimal cooling of steel profiles[END_REF] to boundary control problems. Both papers only consider the LQR problem, the extension to the LQG case is rather straightforward and is executed in [START_REF] Hein | MPC/LQG-Based Optimal Control of Nonlinear Parabolic PDEs[END_REF]Section 8.2.6]. The numerical solution of the resulting large-scale AREs is the computational bottleneck of the suggested control approach. The effective solution of large-scale AREs and associated LQR problems is discussed, e.g., in [START_REF] Banks | A numerical algorithm for optimal feedback gains in high dimensional linear quadratic regulator problems[END_REF][START_REF] Benner | Solving large-scale control problems[END_REF][START_REF] Benner | Numerical solution of large Lyapunov equations, Riccati equations, and linear-quadratic control problems[END_REF][START_REF] Morris | Solution of algebraic Riccati equations arising in control of partial differential equations[END_REF]. The main idea is to apply a Newton-type method to the quadratic nonlinear systems of equations defined by the AREs and to solve the Newton steps by effective iterative methods. It should be noted that for LQG design as discussed here, the actual solution operators/matrices are not necessary as one is only interested in the feedback and estimator gain matrices K and L. Note that the Newton iteration for AREs can be re-written in such a way that one directly iterates on approximations to these operators rather than on approximations to the ARE solutions. This saves a significant amount of workspace and computational effort and is thus recommended in the context of the suggested MPC/LQG scheme. For an efficient variant of the Newton-Kleinman iteration suitable for large-scale AREs, see [START_REF] Benner | A Galerkin-Newton-ADI Method for Solving Large-Scale Algebraic Riccati Equations[END_REF]. Numerical examples demonstrating the effectiveness of the proposed MPC/LQG feedback control design for PDE-constrained optimization problems are shown in [START_REF] Benner | MPC for the Burgers equation based on an LGQ design[END_REF] for the Burgers equation and in [START_REF] Benner | Model predictive control based on an LQG design for time-varying linearizations[END_REF] for a bilinear 3D reaction-diffusion system.

Conclusions

We have presented a framework for model predictive control of infinite-dimensional nonlinear systems subject to stochastic perturbations based on an LQG design to implement the optimization step. This includes the state estimation using a Kalman filter. Linearization about the set point leads to an LTI system. We have focused here on the necessary theoretical ingredients to render this step well-posed. Sufficient conditions for convergence of a numerical approximation scheme to implement the LQG design in a computational procedure can be derived, but are not detailed here due to space restrictions. Though stabilization properties of the nonlinear infinite-dimensional MPC/LQG controller have not been shown yet, numerical experiments in [START_REF] Benner | MPC for the Burgers equation based on an LGQ design[END_REF][START_REF] Benner | Model predictive control based on an LQG design for time-varying linearizations[END_REF] illustrate the good performance of this control scheme. Further improvements can be obtained if one allows for time-varying linearizations in the optimization step, i.e., linearization around the reference trajectory. The treatment of this case is similar to the LTI case and will be described, together with further numerical experiments, in a forthcoming detailed publication. A convergence and stabilization proof of the infinite-dimensional design based on ideas presented in [START_REF] Ito | On asymptotic properties of receding horizon optimal control[END_REF][START_REF] Ito | Receding horizon optimal control for infinite dimensional systems[END_REF][START_REF] Ito | Receding horizon control with incomplete observations[END_REF] is in progress. Further investigations are necessary in order to make the approach real-time feasible. This may require algorithmic improvements in the Riccati solvers or the inclusion of a model reduction strategy in the prediction and optimization step.


Marina Blizorukova
email: msb@imm.uran.ru

On an Algorithm for Dynamic Reconstruction in Systems with Delay in Control

Keywords:

We discuss a problem of the dynamic reconstruction of unknown input controls in nonlinear vector equations. A regularizing algorithm is proposed for reconstructing these controls simultaneously with the processes. The algorithm is stable with respect to informational noises and computational errors.

1 Introduction. Problem statement. Consider a controlled system described by the following equation ẋ(t) = f 1 (t, u t (s), x t (s)) + f 2 (t, x t (s))u(t) [START_REF] Osipov | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF] with the initial state

u t0 (s) = u 0 (s) ∈ C([-τ u m , 0]; R n1 ), x t0 (s) = x 0 (s) ∈ C([-τ x n , 0]; R n2 ). ( 2
)
Here t is time from a fixed interval T = [t 0 , ϑ] (t 0 < ϑ < +∞); x(t) = (x 1 (t), . . . , x n2 (t)) is the phase state of the system; u(t) = (u 1 (t), . . . , u n1 (t)) is a control; the symbols x t (s) and u t (s) mean the functions x t (s) = x(t + s) for s ∈ [-τ x n , 0] and u t (s) = u(t + s) for s ∈ [-τ u m , 0], respectively. We assume that initial state (2) is Lipschitz. For simplicity, we assume also that the initial state x 0 (s), u 0 (s) is fixed and known. The control u = u(t) = (u 1 (t), . . . , u n1 (t)) is called an admissible control if its components u i (t), i ∈ [1 : n 1 ], are Lebesgue measurable functions on the interval T and values u(t) belong to a given compact set P from Euclidean space R n1 for almost all t ∈ T . The set of all admissible controls is denoted by P (•). Therefore, P (•) = {u(•) ∈ L 2 (T ; R n1 ) : u(t) ∈ P for a. a. t ∈ T }. By the trajectory (or the solution) x(•) of equation [START_REF] Osipov | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF] with initial state [START_REF] Maksimov | Dynamical Inverse Problems of Distributed Systems[END_REF] corresponding to some admissible control u(•), we call absolutely continuous on T function x = x(t) satisfying (1) for a.a. t ∈ T .

Condition 1. The elements of matrix function

f 2ij (t, x t (s)) = f 2ij (t, x(t), x(t -τ x 1 ), . . . , x(t -τ x n )), i ∈ [1 : n 2 ], j ∈ [1 : n 1 ],
and vector-valued function

f 1i (t, u t (s), x t (s)) = = f 1i (t, u(t -τ u 1 ), . . . , u(t -τ u m ), x(t), x(t -τ x 1 ), . . . , x(t -τ x n )), i ∈ [1 : n 2 ] satisfy the Lipschitz conditions |f 2ij (t 1 , x (1) 0 , x (1) 1 , . . . , x (1) n ) -f 2ij (t 2 , x (1)
0 , x

(2) 1 , . . . , x (2) n )| ≤ (3) ≤ C 1 (|t 2 -t 1 | + n j=0 |x (1) j -x (2) j |), |f 1i (t 1 , u (1)
1 , . . . , u (1) m , x

0 , x

, . . . , x (1) n ) -f 1i (t 2 , u (1) 1
1 , . . . , u (2) m , x

(2) 0 , x (2) 1
, . . . , x (2) n )| ≤ d 1 (|t 2 -t 1 | + m i=1 |u (1)
i -u

(2)

i | + n j=0 |x (1) j -x (2) j |). (4)
In this case, under this condition for any pair, i.e., for initial state (2) and the control u(•) ∈ P (•), there exists a unique solution of equation [START_REF] Osipov | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF].

Let u(•) be an admissible control realizing during the given time interval T ; x(•) be the real motion generated by this control. We assume that the phase states x(τ i ) of the system are inaccurately measured at frequent enough time moments

τ i ∈ T in the process. Measurement results ξ h (τ i ) ∈ R n2 satisfy the inequalities |ξ h (τ i ) -x(τ i )| ≤ h. (5)
Here, the quantity h ∈ (0, 1) specifies the measurement error.

In the present paper, we construct an algorithm that reconstructs the control u(•) on the basis of the current information ξ h (•) in real time. Since the exact reconstruction is impossible due to the error of measurements ξ h (•) we require that the algorithm should generate some approximation. Namely, it is required to construct an algorithm allowing us, on the basis of the inaccurate measurements ξ h (•), and in real time, to form the admissible control v h (•) such that the meansquare deviation of v h (•) from u(•):, i.e.,

|v h (•) -u(•)| 2 L2(T ) = ϑ t0 |v h (t) -u(t)| 2 dt, (6)
is arbitrarily small for the sufficiently small measurement error h. Since the measurements are inaccurate it is in general impossible to identify u(t) precisely, therefore the problem is to approximate the input by some function v h (t).

Here and below, the symbol | • | stands for both the Euclidean norm and the corresponding matrix norm and for the modulo of a number. In what follows, we set τ u m = τ x n = τ for simplicity, and by ξ h (•) we denote the function

ξ h (t), t ∈ [t 0 -τ, ϑ] such that ξ h (t) = x 0 (t -t 0 ) for t ∈ [t 0 -τ, t 0 ), ξ h (t) = ξ h (τ i ) for t ∈ [τ i , τ i+1 ), i ∈ [0 : d -1], where τ i = τ h,i , d = d h , ξ h (τ i ) satisfies (5).
The suggested solution outline is the following ( [START_REF] Osipov | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF][START_REF] Maksimov | Dynamical Inverse Problems of Distributed Systems[END_REF][START_REF] Maksimov | Dynamical reconstruction of inputs for constraction semigroup systems: the boudary input case[END_REF][START_REF] Maksimov | Dynamical state and control reconstruction for a phase field model. Dynamics of continuous, discrete and impulsive systems[END_REF][START_REF] Osipov | Methods of Dinamical Reconstruction of Inputs of Controlled Systems[END_REF][START_REF] Maksimov | Lyapunov function method in input reconstruction problems of systems with aftereffect[END_REF]). An auxiliary control system (model M ) described by equation of the form

ẇ(t) = F (t, ξ h t (s), v h t (s)), w t0 (s) = w 0 (s), t ∈ T (7)
is associated with the real dynamical system (1). Here the vector w ∈ R n2 characterizes state of the model, the form of function F is corrected below, vector v h is control action. After that, the problem of reconstruction of input u(•) is replaced by the problem of positional control of the model. This process is realized on the time interval T in such a way that control v h (•) "approximates" appropriately u(•). First, one takes a uniform net

∆ = {τ i } m i=0 , τ i+1 = τ i + δ, δ > 0, i ∈ [0 : m], τ 0 = 0, τ m = T with the step δ. Then, on the interval t ∈ [τ i , τ i+1 ) the model is acted upon the controls v h i = V h (τ i , w τi (s), ξ h τi (s)) (8)
calculated at the moment τ i by use of some rule, which hereinafter we shall identify with mapping V h . Thus, the controls in the model are realized by the method of feedback control. Its value on the interval [τ i , τ i+1 ] depends on the measurement results ξ h (•) corresponding to the phase state x(•) of the system (1) and state w of the model (7). The described process forms the piece-wise function

v h (t) = v h i , t ∈ [τ i , τ i+1 )
in real time synchro with the motion of real system [START_REF] Osipov | Inverse Problems for Ordinary Differential Equations: Dynamical Solutions[END_REF]. Thus, to solve the problem above, we should specify a model and a control law for this model.

Algorithm for solving the problem

As a model, we take the following system of linear ordinary differential equation

ẇ(t) = f 1 (τ i , v h τi (s), ξ h τi (s)) + f 2 (τ i , ξ h τi (s))v h i + 2(ξ h (τ i ) -w(τ i )), (9)
w ∈ R n2 , t ∈ [τ i , τ i+1 ), τ i = τ h,i , v h t0 (s) = u 0 (s), with the initial state w(t 0 ) = ξ h (t 0 ). The solution of this equation w(•) = w(•; t 0 , w t0 (s), v h (•)
) is understood in the sense of Caratheodory. So, the righthand side of equation of the model (7) has the form

F (t, ξ h t (s), v h t (s)) = f 1 (τ i , v h τi (s), ξ h τi (s)) + f 2 (τ i , ξ h τi (s))v h i + +2(ξ h (τ i ) -w(τ i )), t ∈ [τ i , τ i+1 ).
Introduce the following notation: ∆ (j) = [t j , t j+1 ], t j = t 0 +τ x 1 j; the symbol l stands for the integer part of the number τ /τ x 1 ; j * = max{j : t j < ϑ},

g j (h) = h (1/3) j , j ∈ [1 : j * ].
Fix a partition of the interval T with a step δ = δ(h) depending on the measurement error h, i.e.,

∆ h = {τ h,i } d h i=0 , τ i = τ h,i , τ h,0 = t 0 , τ h,d h = ϑ, (10)
(for simplicity, we assume that τ i -τ i-1 = δ = δ(h)). Without loss of generality, we can suppose that the partition ∆ h is chosen in such a way that t j ∈ ∆ h . Define the law of forming the control v h i in the model (for

τ i ∈ [t j , t j+1 ) ∩ T ) by the relations V h (τ i , w τi (s), ξ h τi (s)) = V j (τ i , w τi (s), ξ h τi (s)) = arg min{2(l i , f 2 (τ i , ξ h τi (s))v) + α j |v| 2 : v ∈ P }. ( 11
)
Here α j is a parameter, j ∈ [0 :

j * ], l i = w(τ i ) -ξ h (τ i ). Condition 2.
Let n 2 ≥ n 1 , and let there exists a number c * > 0 such that the matrix f 2 (t, x t (s)) has a minor of order n 1 with the property: the n 1 × n 1dimensional matrix f2 (t) = f2 (t, x t (s)) corresponding to this minor satisfies the inequality | f2 (t)u| ≥ c * |u| for each t ∈ T and all u ∈ R n1 . We choose the parameter α j which plays the role of the regularizer, as follows:

α 0 = Ch 2/3 , α j = Cg 2/3 j (h), j ≥ 1, C = const > 0. ( 12
)
Let us describe the algorithm for solving the problem above. Before the initial moment the value h and the partition ∆ = ∆ h with diameter δ = δ(h) are fixed. The work of the algorithm starting at time t = 0 is decomposed into m h -1 steps. At the i-th step carried out during the time interval δ i = [τ i , τ i+1 ), τ i = τ h,i , the following actions take place. First, at time moment τ i vector v h i is calculated by formula (11). Then the control v h (t) = v h i is fed onto the input of the model (9). After that, we transform the state w τi (s) of the model into w τi+1 (s). The procedure stops at time ϑ.

The following theorem is true.

Theorem 1. Let δ = δ(h) ≤ h. Then the inequalities ν (j) ≡ |v h (•) -u(•)| 2 L2(∆ (j-1) ;R n 1 ) ≤ c j g j (h), j ∈ [1 : j * ], are valid. Here, v h (t) = u(t) for t ∈ [t 0 -τ, t 0 ], v h (t) = u 0 (-τ ) for t ∈ [t 0 -τ - τ u 1 , t 0 -τ ).
The proof of the theorem is based on auxiliary statements, which are used in forthcoming considerations. Introduce two systems

ṗ(t) = f 1 (t) + f 2 (t)u 1 (t), t ∈ T, q(t) = F 1 (t) + F 2 (t)u 2 (t), where p(t), q(t) ∈ R n , f 1 (•), F 1 (•) ∈ L 2 (T ; R n ), f 2 (•) ∈ L 2 (T ; R n×r ), F 2 (•) ∈ L 2 (T ; R n×r ), u 1 (•), u 2 (•) ∈ L 2 (T ; R r ), |u l (•)| L∞(T ;R r ) ≤ K, l = 1, 2.
Introduce the notation:

∆ (j) * = [t * j , t * j+1 ] ∩ T , t * j = t 0 + τ * j, j ∈ [0 : j 0 ], ∆ (-1) = [t 0 -τ * , t 0 ], τ * = const ∈ (0, ϑ -t 0 ), j 0 = max{j : t * j ≤ ϑ}.
Let r ≤ n and let there exists a number c > 0 such that the matrix f 2 (t) has a minor of order r such that the r × r-matrix f2 (t) corresponding to this minor satisfies the following inequality: | f2 (t)u| ≥ c|u| for each t ∈ T and all u ∈ R r .

It is easy to verify the following lemmas.

Lemma 1. Let the function t → ( f2 (t)) -1 u 1 (t) be a function of bounded variation on T and let the conditions

|f 1 (•) -F 1 (•)| 2 L2(∆ (j) * ;R n ) ≤ a (j) 1 , |f 2 (•) -F 2 (•)| 2 L2(∆ (j) * ;R n×r ) ≤ a (j) 2 , |p(t) -q(t)| 2 + αj t t * j {|u 2 (ν)| 2 -|u 1 (ν)| 2 } dν ≤ a (j) 3 t ∈ [t * j , t * j+1 ], |p(t * j ) -q(t * j )| 2 ≤ a (j)
4 , αj = const ∈ (0, +∞) be true. Then the inequality

|u 1 (•) -u 2 (•)| 2 L2(∆ (j) * ;R r ) ≤ K j { 4 l=1 (a (j) l ) 1/2 + α1/2 j } + a (j) 3 / αj is valid.
Lemma 2. The bunches of solutions of systems (1) and (9) are bounded in the space W 1,∞ (T ;

R n2 ) = {x(•) ∈ L 2 (T ; R n2 ); ẋ(•) ∈ L 2 (T ; R n2 )}.
We use the relation

ε j (t) = |x(t) -w(t)| 2 + α j t tj {|v h (ν)| 2 -|u(ν)| 2 } dν, j ∈ [0 : j * ], t ∈ T.
Lemma 3. The following inequalities

ε j (t) ≤ b j , t ∈ ∆ (j) ∩ T, j ∈ [0 : j * ],
are valid, where

b j = |x(t j ) -w(t j )| 2 + c (1) j (h + δ) + c (2) j j k=j-l ν (k) , c (1)
j , c

(2) j are some constants, which can be explicitly written.

Proof. Fix τ i ∈ ∆ (j) . Then for t ∈ ∆ (j) ∩ δ i = [τ i , τ i+1 ], we obtain

ε j (t) ≤ ε j (τ i ) + 4 j=1 Λ ji (t), (13)
where

Λ 1i (t) = 2(s i , t τi {f 1 (ν, u ν (s), x ν (s))-f 1 (τ i , v h ν (s), ξ h τi (s))} dν), s i = x(τ i )-w(τ i ), Λ 2i (t) = 2(s i , t τi {f 2 (ν, x ν (s))u(ν) - -f 2 (τ i , ξ h τi (s))v h i } dν) + α j t τi {|v h (ν)| 2 -|u(ν)| 2 } dτ, Λ 3i (t) = -2(t -τ i )(s i , ξ h (τ i ) -w(τ i )), Λ 4i (t) = (t -τ i ) t τi | ẇ(τ ) -ẋ(τ )| 2 dτ.
By virtue of lemma 2, we have

Λ 4i (t) ≤ K (j) * (t -τ i ) 2 , t ∈ δ i . ( 14
) Note that v h (τ i + s) = v h (t + s) for s ≥ t 0 -τ i , t ∈ [τ i , τ i+1
] and in addition

|ξ h (τ i + s) -x(t + s)| ≤ K * (h + t -τ i ) for τ i + s ≥ t 0 -τ. (15)
Taking into account lemma 2, as well as the Lipschitz property of the functions u 0 (s) and x 0 (s), inequalities (4) and the relation

|ξ h (τ i + s) -x(t + s)| ≤ K * (h + t -τ i ) for τ i + s ≥ t 0 -τ, (16)
we obtain for t ∈ δ i the estimate

t τi |f 1 (ν, u ν (s), x ν (s)) -f 1 (τ i , v h ν (s), ξ h τi (s))| dν ≤ ≤ K (j) 1 (t -τ i )(h + t -τ i ) + K (j) 2 (t -τ i ) 1/2 m k=1 t-τ u k τi-τ u k |u(ν) -v h (ν)| 2 dν 1/2 .
Here, τ x 0 = 0. In this case, the inequality

Λ 1i (t) ≤ 2(t -τ i )|x(τ i ) -w(τ i )| 2 + K (j) 3 {(t -τ i )(h + t -τ i ) 2 + + m k=1 t-τ u k τi-τ u k |u(ν) -v h (ν)| 2 dν} (17)
holds for t ∈ δ i . In view of (5), we have

Λ 3i (t) ≤ -2(t -τ i )|x(τ i ) -w(τ i )| 2 + K (j) 4 h(t -τ i ), t ∈ δ i . (18)
Moreover, from ( 5), (3), and (16), we derive

|f 2 (ν, x ν (s))u(ν) -f 2 (τ i , ξ h τi (s))u(ν)| ≤ K 0 (h + ν -τ i ) for ν ∈ [τ i , τ i+1 ]. In this case, Λ 2i (t) ≤ K (j) 5 (t -τ i )(h + t -τ i )+ + t τi {2(l i , f 2 (τ i , ξ h τi (s)){v h i -u(ν)} + α j {|v h i | 2 -|u(ν)| 2 }} dν.
The rule for forming the control v h i (11) and the last inequality imply

Λ 2i (t) ≤ K (j) 5 (t -τ i )(h + t -τ i ). (19)
Finally, taking into account ( 13)-( 19), we conclude that for t ∈ ∆ (j) ∩ δ i

ε j (t) ≤ ε j (τ i ) + K (j) 6 δ(h + δ) + K (j) 3 m k=1 t-τ u k τi-τ u k |u(ν) -v h (ν)| 2 dν, i.e., for t ∈ ∆ (j) = [t j , t j+1 ], ε j (t) ≤ ε j (t j ) + K (j) 7 (h + δ) + K (j) 8 tj+1-τ u 1 tj -τ |u(ν) -v h (ν)| 2 dν.
Note that τ = lτ u 1 + γ, γ ≥ 0. Therefore, t j+1 -τ u 1 = t j , t j-l-1 ≤ t j -τ ≤ t j-l . In this case, for t ∈ ∆ (j) we have

ε j (t) ≤ ε j (t j ) + K (j) 7 (h + δ) + K (j) 9 j k=j-l ν (k) .
Here, constants K (j) k , k ∈ [0 : 9] are written explicitly. Thus, one can assume that c

(1)

j = K (j) 7 and c (2) j = K (j)
9 . The lemma is proved. Lemma 4. Let δ ≤ h and values α j be given by (12). Then the inequalities

ν (j) ≤ c j g j (h), ( 20
) b j ≤ c (0) j g j (h) (21)
are valid.

Proof. For simplicity, set t j * +1 = ϑ. By virtue of lemma 3, we have for t ∈ ∆ (j)

|x(t) -w(t)| ≤ ε j (t) + α j t tj {|v h (ν)| 2 + |u(ν)| 2 } dν 1/2 ≤ b j + α j ρ A 1/2 , (22)
where ρ A = 2τ * d 2 (P ) and d(P ) = sup{|u| : u ∈ P }. Taking into account the inclusion t j ∈ ∆ h , we conclude that for any j ∈ [0 : j * ], one can specify the number i = i j (h) such that t j = τ ij (h) . Introduce the notation

̺ j ≡ |f 1 (•) - F 1 (•)| 2 L2(∆ (j) ;R n 2 )
. In this case, by virtue of lemma 2, as well as of ( 4) and ( 16), we obtain

̺ j ≤ d (1) j i=ij+1(h)-1 i=ij (h) τi+1 τi {δ 2 + h 2 + γ h (ν) + γ h i (ν) + |ξ h (τ i ) -w(τ i )| 2 } dν,
where

γ h (ν) = m k=1 |u(ν -τ u k ) -v h (ν -τ u k )| 2 , γ h i (ν) = n k=0 |x(ν -τ x k ) -ξ h (τ i -τ x k )| 2 .
Note that

tj+1 tj γ h (ν) dν ≤ d (2) j tj t j-l-1 |u(ν) -v h (ν)| 2 dν = d (2) j j k=j-l ν (k) , (23)
tj+1 tj γ h i (ν) dν ≤ d (3)
j (h 2 + δ 2 ). (24)
In addition,

ν (k) = 0 k ∈ [-l : 0]. (25)
Therefore, combining inequalities ( 22)-(24), we obtain the estimates

̺ j ≤ d (5)
j {h 2 + δ 2 + j k=j-l ν (k) + b j + α j }, j ∈ [0 : j * ]. (26)
One can easily see that the following estimates also hold:

|f 2 (•) -F 2 (•)| 2 L2(∆ (j) ;R n 2 ×n 1 ) ≤ d (5) j (h 2 + δ 2 ), j ∈ [0 : j * ]. (27)
Here d

(1)

j -d (5)
j are some constants, which can be explicitly written. By lemma 3, (22), and (25), for δ ≤ h, we have the inequalities

ε 0 (t) ≤ b 0 ≤ c * 0 h, t ∈ ∆ (0) , (28)
with initial conditions x 0 (s) = y 0 (s) = 1, u(s) = 0 for s ∈ [-1, 0] and control u(t) = t was considered. The solution x(t) = {x 1 (t), x 2 (t)} of system (33) was calculated analytically. During the experiment, we assumed that ξ h (τ i ) = x 1 (τ i ) + h. As a model, we took the system (9), which has the form ẇ(0) (t) = 2ξ h 1 (τ i -1) + v h i + 2(ξ h 1 (τ i ) -w (0) (τ i )) for t ∈ [τ i , τ i+1 ) ẇ(1) (t) = ξ h 2 (τ i -1) + ξ h 1 (τ i ) + v h (τ i -1) + 2(ξ h 2 (τ i ) -w (1)

(τ i )), (34)
with the initial condition w (0) (s) = w (1) (s) = 1, for s ∈ [-1, 0]. Here v h (τ i ) = v h i for t ∈ [τ i , τ i+1 ), i ≥ 0, v h (s) = 0 for s ∈ [-1, 0). The controls v h i in model (34) were calculated by the following formula (see (11)) v h i = arg min{2l i v + α j |v| 2 : |v| ≤ K}, where l i = w (0) (τ i ) -ξ h 1 (τ i ). In figures 1 and 2 the results of calculations are presented for the case when δ = 10 -4 , α 0 = Ch 2/3 , α 1 = Ch 2/9 , C = 0.2, K = 10.

Fig. 1

 1 corresponds to the case when h = 0.001, fig.2-h = 0.02. In these figures the solid (dashed) lines represent the model control v h (•) (the real control u(•)). The equations were solved by the Euler method with step δ.

  Fig. 1. Fig. 2.

Foundation for Basic Research (12-01-00175-a), by the Ural-Siberian Integration Project (12-C-1-1017), and by the Program for support of leading scientific schools of Russia (6512.2012.1).

Taking into account (25)-(28), for h ∈ (0, 1), we obtain

By virtue of condition 1, one can use lemma 1. Set p = x, q = w, u

4 = c * h 2/3 , α0 = α 0 = ch 2/3 , we have

It means that inequality (20) holds for j = 1. Further, by using ( 29) and (30), we deduce that

Inequality (21) for j = 1 is also verified. It follows from (22) that

Consequently, in view of relations (31), as well as of the rule for definition b j , we have the inequality

Setting a

(5)

= b j-1 + ρ A α j-1 , j ∈ [1 : j * ] for j ≥ 1 in lemma 1 and taking into account inequalities (32), we obtain

Here, we used lemma 3 and inequalities ( 27), (28), and (31)) for choosing values a

i . Now, to proof inequalities (20) and ( 21), one can use the proof by induction. The lemma is proved.

Example

The algorithm was tested by a model example. The following system ẋ1 (t) = 2x 1 (t -1) + u(t) ẋ2 (t) = x 2 (t -1) + x 1 (t) + u(t -1), t ∈ T = [0, 2], (33)


Salah- Eddine Rebiai
email: rebiai@hotmai.com

Exponential stability of the system of transmission of the wave equation with a delay term in the boundary feedback

Keywords: Wave equation, transmission problem, time delays, boundary stabilization, exponential stability

We consider a system of transmission of the wave equation with Neumann feedback control that contains a delay term and that acts on the exterior boundary. First, we prove under some assumptions that the closed-loop system generates a C0-semigroup of contractions on an appropriate Hilbert space. Then, under further assumptions, we show that the closed-loop system is exponentially stable. To establish this result, we introduce a suitable energy function and use multiplier method together with an estimate taken from [3] (Lemma 7.2) and compactnessuniqueness arguments.

Introduction

It is by now well-known that certain infinite-dimensional second-order systems are not robust with respect to arbitarily small delays in the damping. This lack of stability robustness was first shown to hold for the one-dimensional wave equation ( [START_REF] Datko | An example on the effect of time delays in boundary feedback stabilization of wave equations[END_REF]). Later, further examples illustrating this phenomenon were considered in [START_REF] Datko | Not all feedback stabilized hyperbolic systems are robust with respect to small time delays in their feedbacks[END_REF]: the two-dimensional wave equation with damping introduced through Neumann-type boundary conditions on one edge of a square boundary and the Euler-Bernoulli beam equation in one dimension with damping introduced through a specific set of boundary conditions on the right end point.

Recently, Xu et al [START_REF] Xu | Stabilization of wave systems with input delay in the boundary control[END_REF] established sufficient conditions that guarantee the exponential stability of the one-dimensional wave equation with a delay term in the boundary feedback. Nicaise and Pignotti [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] extended this result to the multi-dimensional wave equation with a delay term in the boundary or internal feedbacks. The same type of result was obtained by Nicaise and Rebiai [START_REF] Nicaise | Stabilization of the Schrödinger equation with a delay term in boundary feedback or internal feedback[END_REF] for the Schrödinger equation.

Motivated by the references [START_REF] Xu | Stabilization of wave systems with input delay in the boundary control[END_REF], [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] and [START_REF] Nicaise | Stabilization of the Schrödinger equation with a delay term in boundary feedback or internal feedback[END_REF]; we investigate in this paper the problem of exponential stability for the system of transmission of the wave equation with a delay term in the boundary feedback.

Let Ω be an open bounded domain of R n with a boundary Γ of class C 2 which consists of two non-empty parts Γ 1 and

Γ 2 such that Γ 1 ∩ Γ 2 = ∅. Let Γ 0 with Γ 0 ∩ Γ 1 = Γ 0 ∩ Γ 2 = ∅
be a regular hypersurface of class C 2 which separates Ω into two domains Ω 1 and Ω 2 such that Γ 1 ⊂ ∂Ω 1 and Γ 2 ⊂ ∂Ω 2 . Furthermore, we assume that there exists a real vector field h ∈ (C 2 (Ω)) n such that: (H.1) The Jacobian matrix J of h satisfies

Ω J(x)ζ(x).ζ(x)dΩ ≥ α Ω |ζ(x)| 2 dΩ,
for some constant α > 0 and for all

ζ ∈ L 2 (Ω; R n ); (H.2) h(x).ν(x) ≤ 0 on Γ 1 ; (H.3) h(x).ν(x) ≥ 0 on Γ 0 .
where ν is the unit normal on Γ or Γ 0 pointing towrds the exterior of Ω or Ω 1 .

Let a 1 , a 2 > 0 be given. Consider the system of transmission of the wave equation with a delay term in the boundary conditions:

y ′′ (x, t) -a(x)∆y(x, t) = 0 in Ω × (0, +∞), (1)
y(x, 0) = y 0 (x), y ′ (x, 0) = y 1 (x, 0) in Ω, (2)
y 1 (x, t) = 0 on Γ 1 × (0, +∞), (3)
∂y 2 (x, t) ∂ν = -µ 1 y ′ 2 (x, t) -µ 2 y ′ 2 (x, t -τ ) on Γ 2 × (0, +∞), (4)
y 1 (x, t) = y 2 (x, t), on Γ 0 × (0, +∞), (5)
a 1 ∂y 1 (x, t) ∂ν = a 2 ∂y 2 (x, t) ∂ν on Γ 0 × (0, +∞), (6)
y ′ 2 (x, t -τ ) = f 0 (x, t -τ ) on Γ 2 × (0, τ ). (7)
where:

-a(x) = a 1 , x ∈ Ω 1 a 2 , x ∈ Ω 2 -y(x, t) = y 1 (x, t), (x, t) ∈ Ω 1 × (0, +∞) y 1 (x, t), (x, t) ∈ Ω 2 × (0, +∞) -∂.
∂ν is the normal derivative. -µ 1 and µ 2 are positive real numbers.

τ is the time delay y 0 , y 1 , f 0 are the initial data which belong to suitable spaces.

In the absence of delay, that is µ 2 = 0, Liu and Williams [START_REF] Liu | The exponential stability of the problem of transmission of the wave equation[END_REF] have shown that the solution of (1)-( 6) decays exponentially to zero in the energy space

H 1 Γ1 (Ω) × L 2 (Ω) provided that a 1 > a 2 (8)
and {Ω, Γ 0 , Γ 1 , Γ 2 } satisfies (H.1), (H.2), (H.3), and (H.4) h(x).ν(x) ≥ γ > 0.

The purpose of this paper is to investigate the stability of problem ( 1) -(7) in the case where both µ 1 and µ 2 are different from zero. To this end, assume as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] that

µ 1 > µ 2 . ( 9
)
and define the energy of a solution of ( 1) -( 7) by

E(t) = 1 2 Ω |y ′ (x, t)| 2 + a(x) |∇(y(x, t)| 2 dx+ ξ 2 Γ2 1 0 |y(x, t -τ ρ)| 2 dρ dσ(x), (10)
where

a 2 τ µ 2 < ξ < a 2 τ (2µ 1 -µ 2 ), (11)
We show that if {Ω, Γ 0 , Γ 1 , Γ 2 } satisfies (H.1), (H.2) and (H.3), then there is an exponential decay rate for E(t).The proof of this result combines multipliers technique and compactness-uniqueness arguments.

The main result of this paper can be stated as follows.

Theorem 1. Assume (H1), (H.2), (H.3), ( 8) and [START_REF] Xu | Stabilization of wave systems with input delay in the boundary control[END_REF]. Then there exist constants M ≥ 1 and ω > 0 such that

E(t) ≤ M e -ωt E(0).
Theorem 1 is proved in Section 3. In Section 2, we investigate the wellposedness of system ( 1) -(7) using semigroup theory.

2 Well-poseness of problem ( 1) -(7)

Inspired from [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] and [START_REF] Nicaise | Stabilization of the Schrödinger equation with a delay term in boundary feedback or internal feedback[END_REF], we introduce the auxilliary variable z(x, ρ, t) = y(x, tτ ρ). With this new unknown, problem (1) -( 7) is equivalent to

y ′′ (x, t) -a(x)∆y(x, t) = 0 in Ω × (0, +∞), (12)
y(x, 0) = y 0 (x), y ′ (x, 0) = y 1 (x) in Ω, (13)
y(x, t) = 0 on Γ 1 × (0, +∞), (14) ∂z(x, ρ, t) ∂t + 1 τ ∂z(x, ρ, t) ∂ρ = 0 on Γ 2 × (0, +∞) (15) ∂y 2 (x, t) ∂ν = -µ 1 y ′ 2 (x, t) -µ 2 z(x, 1, t) on Γ 2 × (0, +∞), (16)
y 1 (x, t) = y 2 (x, t) on Γ 0 × (0, +∞), (17)
a 1 ∂y 1 (x, t) ∂ν = a 2 ∂y 2 (x, t) ∂ν on Γ 0 × (0, +∞), ( 18
) z(x, 0, t) = y ′ (x, t) on Γ 2 × (0, +∞) (19) z(x, ρ, 0) = f 0 (x, -τ ρ) on Γ 2 × (0, 1) (20)
Now, we endow the Hilbert space

H =H 1 Γ1 (Ω) × L 2 (Ω) × L 2 (Γ 2 ; L 2 (0, 1))
with the inner product

  u v z   ;   u v z   = Ω (a(x)∇u(x)∇u(x)+v(x)v(x)) dx+ξ Γ2 1 0 z(x, ρ)z(x, ρ)dρ dσ(x)
and define a linear operator in H by

D(A) = {(u, v, z) T ∈ H 2 (Ω 1 , Ω 2 , Γ 1 ) × H 1 Γ1 (Ω) × L 2 (Γ 2 ; H 1 (0, 1); ∂u ∂ν = -µ 1 v -µ 2 z(., 1), v = z(., 0) on Γ 2 } (21) A   u v z   =   v a(x)∆u -τ -1 ∂z ∂ρ   (22)
The spaces used for the definition of H and D(A) are

H 1 Γ1 (Ω) = {u ∈ H 1 (Ω) : u = 0 on Γ 1 } H 2 (Ω 1 , Ω 2 , Γ 1 ) = {u i ∈ H 2 (Ω i ) : u = 0 on Γ 1 , u 1 = u 2 and a 1 ∂u 1 ∂ν = a 2 ∂u 2 ∂ν on Γ 0 }
Then we can rewrite ( 12) -(20) as an abstract Cauchy problem in

H d dt Y (t) = AY (t) Y (0) = Y 0 (23)
Y (.) ∈ C([0, +∞); H) if Y 0 ∈ H, Y (.) ∈ C([0, +∞); D(A)) ∩ C 1 ([0, +∞); H) if Y 0 ∈ D(A). Proof. Let Y =   u v z   ∈ D(A). Then AY, Y = Ω a(x)∇u(x).∇v(x)dx + Ω (a(x)∆u(x))v(x)dx - ξ τ Γ2 1 0 z ρ (x, ρ)z(x, ρ)dρdΓ (24)
Applying Green's first theorem, we obtain

Ω (a(x)∆u(x))v(x)dx = a 1 Γ1 v(x) ∂u(x) ∂ν dΓ -a 1 Ω1 ∇u(x).∇v(x)dx + a 2 Γ2 v(x) ∂u(x) ∂ν dΓ -a 2 Ω2 ∇u(x).∇v(x)dx = a 2 Γ2 v(x){-µ 1 v(x) -µ 2 z(x, 1)}dΓ - Ω a(x)∇u(x).∇v(x)dx (25)
Integrating by parts in ρ, we get

Γ2 1 0 z ρ (x, ρ)z(x, ρ)dρdΓ = 1 2 Γ2 {z 2 (x, 1) -z 2 (x, 0)}dΓ (26)
Inserting ( 25) and ( 26) into (24) results in

AY, Y = -a 2 µ 1 Γ2 v 2 (x)dΓ -a 2 µ 2 Γ2 v(x)z(x, 1)dΓ - ξ 2τ Γ2 z 2 (x, 1)dΓ + ξ 2τ Γ2 v 2 (x)dΓ
from which follows using the Cauchy-Schwarz inequality

AY, Y ≤ -(a 2 µ 1 - a 2 µ 2 2 + ξ 2τ ) Γ2 v 2 (x)dΓ -( ξ 2τ - a 2 µ 2 2 ) Γ2 z 2 (x, 1)dΓ (27) (27) implies that AY, Y ≤ 0
Thus A is dissipative. Now we show that for a fixed λ > 0 and (g, h, k) T ∈ H, there exists Y = (u, v, z) T ∈ D(A) such that

(λI -A)Y = (g, h, k) T or equivalently λu -v = g (28) λv -a(x)∆u = h (29) λz + 1 τ z ρ = k (30)
Suppose that we have found u with the appropriate regularity, then we can determine z. Indeed, from ( 21) and (30) we have

z ρ (x, ρ) = -λτ z(x, ρ) + τ k(x, ρ) z(x, 0) = v(x)
The unique solution of the above initial value problem is z(x, ρ) = e -λτ ρ v(x) + τ e -λτ ρ ρ 0 e λτ s k(x, s)ds and in particular z(x, 1) = λe -λτ u(x) + z 0 (x), x ∈ Γ 2 where z 0 (x) = -e -λτ g(x) + τ e -λτ 1 0 e λτ s k(x, s)ds By ( 28) and ( 29), the function u satisfies

λ 2 u -a(x)∆u = h + λg (31)
Problem (31) can be reformulated as

Ω (λ 2 u -a(x)∆u)wdx = Ω (h + λg)wdx, w ∈ H 1 Γ1 (Ω) (32)
Using Green's first theorem and recalling (21), we express the right-hand side of (32) as follows

Ω (λ 2 u -a(x)∆u)wdx = Ω (λ 2 uw + a(x)∇u.∇w)dx + a 2 Γ2 {µ 1 (λu -g)w +µ 2 (λe -λτ u(x) + z 0 (x))w}dΓ
Therefore (32), can be rewritten as

Ω (λ 2 uw + a(x)∇u.∇w)dx + a 2 Γ2 (µ 1 + µ 2 e -λτ )λuwdΓ = Ω (h + λg)wdΓ +a 2 µ 1 Γ2 gwdΓ -a 2 µ 2 Γ2 z 0 wdΓ, ∀w ∈ H 1 Γ1 (Ω). ( 33
)
Since the left-hand side of (33) is coercive on H 1 Γ1 (Ω), the Lax-Milgram Theorem guarantees the existence and uniqueness of a solution y ∈ H 1 Γ1 (Ω) of (31). If we consider w ∈ D(Ω) in (28), then y is a solution in D ′ (Ω) of

λ 2 u -a(x)∆u = h + λg (34)
and thus ∆u ∈ L 2 (Ω). Combining (33) together with (34), we obtain after using Green's first theorem

a 2 Γ2 (µ 1 + µ 2 e -λτ )λuwdΓ + a 2 Ω ∂u ∂ν wdΓ = a 2 µ 1 Γ2 gwdΓ -a 2 µ 2 Γ2 z 0 wdΓ which implies that ∂u(x) ∂ν = -µ 1 v(x) -µ 2 z(x, 1)
So, we have found (u, v, z) T ∈ D(A) which satisfies (28) -(30). Thus, by the Lumer-Phillips Theorem (see for instance [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF], Theorem 1.4.3), generates a strongly continuous semigroup of contractions on H.

Proof of Theorem 1

We prove Theorem 1 for smooth initial data. The general case follows by a standard density argument. We proceed in several steps.

Step 1.

Since

E(t) = 1 2 (y, y ′ , z) 2 H
Then, we deduce from the proof of Proposition 1 that E(t) is non-increasing and

d dt E(t) ≤ -C Γ2 {y 2 (x, t) + y ′2 (x, t)}dΓ (35)
where

C = min{a 2 µ 1 - a 2 µ 2 2 + ξ 2τ , ξ 2τ - a 2 µ 2 2 } Step 2. Set E(t) = E(t) + E d (t)
where

E(t) = 1 2 Ω {a(x) |∇y(x, t)| 2 + |y ′ (x, t)| 2 }dx
and

E d (t) = ξ 2τ Γ2 1 0 |y ′ (x, t -τ ρ)| 2 dρdΓ
E d (t) can be rewritten via a change of variable as

E d (t) = ξ 2τ 2 t+τ t Γ2 y ′2 (x, s -τ )dΓ ds (36)
From (36), we obtain

E d (t) ≤ C 1 T 0 Γ2 y ′2 (x, s -τ )dΓ ds (37)
for 0 ≤ t ≤ T and T large enough.

Step 3.

By applying energy methods (multiplier 2h.∇y +(divh-α)y) (see the appendix) to problem ( 1) -(7), we obtain for all T > 0.

T 0

E(t)dt ≤ C 2 (E(0) + E(T )) + C 3 T 0 Γ2 {( ∂y(x, t) ∂ν ) 2 + y ′2 (x, t)}dΓ dt + C 4 T 0 Γ2 |∇ σ y(x, t)| 2 dΓ dt + C 5 T 0 Ω |y(x, t)| 2 dΩdt (38)
where ∇ σ y is the tangential gradient of y.

Step 4.

We eliminate the tangential gradient from (38) by using the following estimate due to Lasiecka and Triggiani (Lemma 7.2 in [START_REF] Lasiecka | Uniform stabilization of the wave equation with Dirichlet or Neumann feedback control without geometrical conditions[END_REF])

T -ǫ ǫ Γ2 |∇ σ y(x, t)| 2 dΓ dt ≤ C 6 { T 0 Γ2 {( ∂y(x, t) ∂ν ) 2 + y ′2 (x, t)}dΓ dt + y 2 L 2 (0,T ;H 1/2+δ (Ω)) }
where ǫ and δ are arbitrary positive constants. We obtain

T 0 E(t)dt ≤ C 2 (E(0) + E(T )) + C 7 T 0 Γ2 {( ∂y(x, t) ∂ν ) 2 + y ′2 (x, t)}dΓ dt + C 8 y 2 L 2 (0,T ;H 1/2+δ (Ω)) (39)
Step 5.

We differentiate E(t) with respect to t and apply Green's first theorem. We obtain

d dt E(t) = a 2 Γ2 y ′ (x, t) ∂y(x, t) ∂ν dΓ dt (40)
From (40), we get via the Cauchy-Schwarz inequality

E(0) ≤ E(T ) + a 2 2 T 0 Γ2 {y ′2 (x, t) + ( ∂y(x, t) ∂ν ) 2 }dΓ dt (41)
Insertion of (41) into (39) yields

T 0 E(t)dt ≤ 2C 2 E(T ) + C 9 T 0 Γ2 {( ∂y(x, t) ∂ν ) 2 + y ′2 (x, t)}dΓ dt + C 8 y 2 L 2 (0,T ;H 1/2+δ (Ω)) (42)
Step 6.

Since E(t) is non-increasing and E(t) = E(t) + E d (t), then (42) together with (37) implies that

T E(T ) ≤ 2C 2 E(T ) + C 9 T 0 Γ2 {( ∂y(x, t) ∂ν ) 2 + y ′2 (x, t)}dΓ dt + C 8 y 2 L 2 (0,T ;H 1/2+δ (Ω)) + T C 1 T 0 Γ2 y ′2 (x, t -τ )dΓ dt (43)
Thus invoking again the identity E(t) = E(t) + E d (t) and recalling the boundary condition (4), we obtain from ( 43)

E(T ) ≤ C 10 T 0 Γ2 {y ′2 (x, t)+y ′2 (x, t-τ )}dΓ dt+C 11 y 2 L 2 (0,T ;H 1/2+δ (Ω)) (44)
for T large enough.

Step 7.

We drop the lower order term on the right-hand side of (44) by a compactnessuniqueness argument to obtain

E(T ) ≤ C 12 T 0 Γ2
{y ′2 (x, t) + y ′2 (x, tτ )}dΓ dt (45)

Step 8.

The estimate (45) together with (35) yields

E(T ) ≤ C 13 1 + C 13 E(0) (46)
The desired conclusion follows now from (46) since the system (1) -( 7) is invariant by translation.

Appendix: Sketch of Proof of (38)

We multiply both sides of (1) by 2h.∇y+(divh-α)y and integrate over (0, T )×Ω. We obtain

Proposition 1 .

 1 where Y (t) = (y, y ′ , z) T and Y 0 = (y 0 , y 1 , f 0 (., -.τ )) T The operator A defined by (21) and (22) generates a strongly continuous semigroup on H. Thus, for every Y 0 ∈ H, problem (23) has a unique solution Y whose regularity depends on the the initial datum Y 0 as follows:

  h.∇y + (divhα)y ′ y} T 0 dΩ -

  after using the boundary conditions (3) and[START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF]. Identity (47) is used together with (H.1), (H.2), (H.3) and (8) to obtain estimate (38).

					∂y 1 ∂ν	2	h.νdΓ dt -(a 1 -a 2 )	0	T	Γ0	|∇y 1 |	2 h.νdΓ dt -
	(a 1 -a 2 ) 2 a 2	0	T		Γ0	∂y 1 ∂ν	2	h.νdΓ dt +	0	T	Γ2	|y ′ 2 |	2 h.νdΓ dt +
	2a 2	0	T	Γ2	∂y 2 ∂ν	2	h.∇y 2 dΓ dt -a 2	0	T	Γ2	|∇y 2 |	2 h.νdΓ dt +
	a 2	0	T	Γ2	∂y 2 ∂ν	2	(divh -α)dΓ dt	(47)


Alexander Tarasyev
email: tarasiev@iiasa.ac.atanastasy.ousova@gmail.com

Anastasia Usova

Nonlinear Stabilizers in Optimal Control Problems with Infinite Time Horizon

Keywords: optimal control, nonlinear control system, nonlinear stabilizer, economic systems

In optimal control problems with infinite time horizon, arising in models of economic growth, there are essential difficulties in analytical and even in numerical construction of solutions of Hamiltonian systems. The problem is in stiff properties of differential equations of the maximum principle and in non-stable character of equilibrium points connected with corresponding transversality conditions. However, if a steady state exists and meets several conditions of regularity then it is possible to construct a nonlinear stabilizer for the Hamiltonian system. This stabilizer inherits properties of the maximum principle, generates a nonlinear system with excluded adjoint variables and leads its trajectories to the steady state. Basing on the qualitative theory of differential equations, it is possible to prove that trajectories generated by the nonlinear stabilizer are close to solutions of the original Hamiltonian system, at least locally, in a neighborhood of the steady state. This analysis allows to create stable algorithms for construction of optimal solutions.

Introduction

This paper deals with optimal control problems with infinite time horizon basing on economical growth models which is relied on classical constructions of growth theory (see [START_REF] Shell | Applications of Pontryagin's Maximum Principle to Economics[END_REF], [START_REF] Solow | Growth Theory: An Exposition[END_REF]). Also it includes ideas of a SEDIM model [START_REF] Sanderson | The SEDIM Model: Version 0.1[END_REF] describing the role of different economic factors such as the demographic ones in a country's economic development. Another technique in the background (see [START_REF] Ayres | The Economic Growth Engine: How Energy and Work Drive Material Prosperity[END_REF], [START_REF] Krasovskii | Conjugation of Hamiltonian Systems in Optimal Control Problems[END_REF]) considers capital and useful work as the key drivers of economic growth and uses optimal control theory to design past and future growth trajectories.

The research of optimal control problems uses as basis the Pontryagin's maximum principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF] for the problem with infinite time horizon (see [START_REF] Aseev | The Pontryagin Maximum Principle and Optimal Economic Growth Problems[END_REF], [START_REF] Balder | An existance result for optimal economic growth problems[END_REF], [START_REF] Krasovskii | Conjugation of Hamiltonian Systems in Optimal Control Problems[END_REF]). We investigate properties of the maximized Hamiltonian function and provide analysis of existence of steady states in domains of specific control regimes and focus attention on the domain corresponding to the transient control regimes of investment. We consider linearized Hamiltonian system in this domain. Special attention is given to the Jacobi matrix which has two negative and two positive eigenvalues that is the steady state has the saddle character. According to the results of the qualitative theory of differential equations [START_REF] Ph | Ordinary Differential Equations[END_REF] the trajectory of the nonlinear Hamiltonian dynamics converges to the steady state tangentially to the plane generated by eigenvectors corresponding to negative eigenvalues of the Jacobi matrix. This analysis provides the important information about the growth rates of optimal synthetic trajectories.

A novelty of the proposed solution is based on the idea of creating of nonlinear stabilizers built on the feedback principle (see [START_REF] Krasovskii | Control under Lack of Information[END_REF], [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF]) which lead the system from any current position to a steady state. The constructed nonlinear stabilizer generates the dynamic system closed in phase variables and having the property of local stability. Also we construct solutions of the Hamiltonian system and the stabilized Hamiltonian system in a steady state neighborhood and compare behavior of these trajectories. Simulated optimal trajectories of nonlinear Hamiltonian systems are obtained numerically by the implicit Runge-Kutta method.

Two-sectors economical growth model and optimal control problem

The model.

The model is based on analysis of the Gross Domestic Product (GDP) dynamics which is denoted by symbol Y . It is supposed that changes of GDP depend on three production factors: capital stock K, labor L (or it can be named as human capital) and useful work U . The production function F describes the relation between these factors and GDP (Y ), that is

Y = F [K, L, U ].
It is assumed that the production function F has the property of homogeneity of degree one, i.e.

F [αK, αL, αU ] = αF [K, L, U ] ∀α > 0.
This model includes also a parameter P (t) denoting the number of workers in a country at time t. According to the Sanderson model [START_REF] Sanderson | The SEDIM Model: Version 0.1[END_REF] we assume that labor L is proportional to the number of workers P with coefficient E. This coefficient has the sense of labor efficiency of one worker. Hence, we have the following equality: L(t) = E(t)P (t). Due to the homogeneous property of the production function we introduce relative variables: k = K/P, l = E = L/P, u = U/P, y = Y /P.

It is supposed that the total number of workers P has exponential growth trend Ṗ (t) = ρP (t), ρ > 0.

(

) 1
Here ρ is a positive predefined constant denoting the relative growth rate. It should be mentioned that the considered dynamics for the labor force is quite adequate for the US statistical data in the period from 1900 to 2005. Parameter ρ is small enough and equal to (approx.) 10

Investments in capital constitute a part of GDP (Y ). Hence, it can be written as follows: S(t) = s(t)Y (t), where function s(t) may take any value in the range from zero to the positive constant a s which is less than one, i.e. 0 s(t) a s < 1. Changes in labor are described by the equation:

L(t) = bR(t), (3)
F [k, l, u] = F [k, l, u] = f (k, l).
Based on equations (1), ( 2) and (3), one can evaluate dynamics of relative variables k and l.

Let functions C(t) and c(t) describe the total consumption level in a country and the consumption level per one worker, respectively. It is assumed that the closed economical system is considered in which GDP (Y ) is spent on consumption (C) and investments in capital stock (S) and human capital (R): Y (t) = C(t) + S(t) + R(t), or in relative variables: y(t) = c(t) + (s(t) + r(t))y(t). Hence one can easily calculate consumption per one worker

c(t) = (1 -s(t) -r(t))y(t) ≈ (1 -s(t))(1 -r(t))y(t). (4)
Optimal control problem.

Let us consider investments s and r as control variables. It is supposed that the utility function of the growth process is described by an integral consumption index discounted on the infinite horizon. We use the consumption index of the logarithmic type, rather common for the theory of endogenous growth (see [START_REF] Tarasyev | Optimal Dynamics of Innovation in Models of Economic Growth[END_REF]). Let us note that the utility of such type is closely related to the notion of entropy in thermodynamics, mechanics and

dynamic systems J = +∞ ∫ t0 e -λt ln c(t) dt.
Here, parameter λ is the discounting factor.

It should be mentioned that the following equality d ln c(t) = dc(t) c(t) determines relative growth of the consumption c(t) (4) per one worker. In fact, the introduced utility function presents the summary growth of the relative consumption adjusted to the value of money depreciation.

Problem 1. The optimal control problem presumes maximization of the utility function

J = +∞ ∫ t0 e -λt (ln (1 -s(t)) + ln (1 -r(t)) + ln f (k(t), l(t))) dt over trajectory (k(•), l(•), s(•), r(•)) of the system { k(t) = s(t)f (k(t), l(t)) -(δ + ρ)k(t) l(t) = br(t)f (k(t), l(t)) -ρl(t) with control parameters (s(•), r(•)) subject to constraints 0 ≤ s(t) ≤ a s < 1, 0 ≤ r(t) ≤ a r < 1, 0 ≤ a s + a r < 1, (5)
and phase variables (

k(•), l(•)) satisfying initial conditions k(t 0 ) = k 0 , l(t 0 ) = l 0 .
The production function y = f (k, l) meets the following conditions P F 1 . For all positive values of phase variables k and l function f (k, l) is positive with its partial derivatives, i.e.

f (k, l) > 0, f k > 0, f l > 0. P F 2 .
For all positive values of phase variables k and l function f (k, l) is a strictly concave function in phase variables, i.e. f kk < 0, f kk f ll -f 2 kl > 0. Here we use the following notations for the first and second order derivatives of the production function f = f (k, l)

f k = ∂f (k, l) ∂k , f l = ∂f (k, l) ∂l , f kl = ∂ 2 f (k, l) ∂k∂l , f kk = ∂ 2 f (k, l) ∂k 2 , f ll = ∂ 2 f (k, l) ∂l 2 .
Let us note that the problem 1 can be solved within the optimal control theory for problems with infinite horizon (see [START_REF] Aseev | The Pontryagin Maximum Principle and Optimal Economic Growth Problems[END_REF], [START_REF] Krasovskii | Conjugation of Hamiltonian Systems in Optimal Control Problems[END_REF]).

Model analysis

Model analysis is based on the Pontryagin maximum principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF] for problems with infinite time horizon [START_REF] Aseev | The Pontryagin Maximum Principle and Optimal Economic Growth Problems[END_REF].

Hamiltonian function.

We investigate properties of the Hamiltonian function H = H(t; k, l; s, r; ψ 1 , ψ 2 ) which is defined by the equality:

H(t; k, l; s, r; ψ 1 , ψ 2 ) = e -λt (ln (1 -s) + ln (1 -r) + ln f (k, l)) + + ψ 1 (sf (k, l) -(δ + ρ)k) + ψ 2 (brf (k, l) -ρl). (6)
Let us formulate the main property of the Hamiltonian function [START_REF] Krasovskii | Control under Lack of Information[END_REF].

Proposition 1. The Hamiltonian function H(t; k, l; s, r; ψ 1 , ψ 2 ) is concave in control variables s and r.

It is convenient to introduce new variables for excluding the exponential time term: ψ 1 = ψ 1 e λt , ψ 2 = ψ 2 e λt and H = He λt . Substituting new variables to the Hamiltonian function [START_REF] Krasovskii | Control under Lack of Information[END_REF] we get the expression:

H(k, l; s, r; ψ 1 , ψ 2 ) = ln (1 -s) + ln (1 -r) + ln f (k, l) + + ψ 1 (sf (k, l) -(δ + ρ)k) + ψ 2 (brf (k, l) -ρl). (7)
Since control variables s and r satisfy to restrictions (5), the optimal control has the following structure:

s 0 =            0, (k, l, ψ 1 ) ∈ ∆ 1 s = {(k, l, ψ 1 ) : ψ 1 f (k, l) 1} ; 1 - 1 ψ 1 f (k, l) , (k, l, ψ 1 ) ∈ ∆ 2 s = { (k, l, ψ 1 ) : 1 ψ 1 f (k, l) 1 1 -a s } ; a s , (k, l, ψ 1 ) ∈ ∆ 3 s = { (k, l, ψ 1 ) : ψ 1 f (k, l) 1 1 -a s } ; r 0 =            0, (k, l, ψ 2 ) ∈ ∆ 1 r = {(k, l, ψ 2 ) : bψ 2 f (k, l) 1} ; 1 - 1 bψ 2 f (k, l) , (k, l, ψ 2 ) ∈ ∆ 2 r = { (k, l, ψ 2 ) : 1 bψ 2 f (k, l) 1 1 -a r } ; a r , (k, l, ψ 2 ) ∈ ∆ 3 r = { (k, l, ψ 2 ) : bψ 2 f (k, l) 1 1 -a r } . (8)
Substituting values of optimal control to the Hamiltonian function H(•) in [START_REF] Krasovskii | Game-Theoretical Control Problems[END_REF] we obtain the maximized Hamiltonian: H(k, l; ψ 1 , ψ 2 ) = H(k, l; s 0 , r 0 ; ψ 1 , ψ 2 ). There exist nine domains D ij = ∆ i s ∩ ∆ j r (i, j = 1, 2, 3) of definition of the maximized Hamiltonian function. These domains are determined by the structure of optimal controls. Let us discuss important properties of the maximized Hamiltonian function.

Proposition 2. The maximized Hamiltonian function H(k, l; ψ 1 , ψ 2 ) is a smooth function in variables k, l and ψ 1 , ψ 2 in domains D ij (i, j = 1, 3) and on boundaries between these domains. Proposition 3. The maximized Hamiltonian is a strictly concave function in phase variables k, l for all positive values of conjugate variables ψ 1 and ψ 2 , if the following matrix is negatively defined:

∂f (k, l) =   -f f k f l f k f kk f kl f l f lk f ll   , ∀ (k, l, ψ 1 , ψ 2 ) ∈ D 22 , ψ 1 > 0, ψ 2 > 0.
Necessary and sufficient conditions of optimality.

Let us mention that for the control problem 1 all conditions of the existence theorem (see [START_REF] Aseev | The Pontryagin Maximum Principle and Optimal Economic Growth Problems[END_REF], [START_REF] Balder | An existance result for optimal economic growth problems[END_REF]) are fulfilled. Moreover, one can formulate necessary [START_REF] Aseev | The Pontryagin Maximum Principle and Optimal Economic Growth Problems[END_REF] and sufficient [START_REF] Krasovskii | Conjugation of Hamiltonian Systems in Optimal Control Problems[END_REF] conditions of optimality for problems with infinite horizon in the form of the Pontryagin maximum principle. It should be noted that properties 2 and 3 ensure sufficiency of necessary optimality conditions [START_REF] Krasovskii | Conjugation of Hamiltonian Systems in Optimal Control Problems[END_REF].

Qualitative analysis.

Firstly, we construct the Hamiltonian system and investigate the existence of steady states. Due to the structure of the optimal control (s 0 (t), r 0 (t)) in [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF] the Hamiltonian system has different form in each domain D ij (i, j = 1, 3). The special attention is given to domain D 22 with the transient control regime, where both controls are not constant.

In the domain D 22 = ∆ 2 s ∩ ∆ 2 r the Hamiltonian system has the following form:

                     k = f (k, l) -(δ + ρ)k - k z 1 = H 1 , l = bf (k, l) -ρl - l z 2 = H 2 , ż1 = ( λ -f k (k, l) + f (k, l) k ) z 1 -b k l f k (k, l)z 2 + k f (k, l) f k (k, l) -1 = H 3 , ż2 = - l k f l (k, l)z 1 + ( λ -bf l (k, l) + b f (k, l) l ) z 2 + l f (k, l) f l (k, l) -1 = H 4 , (9)
where new adjoint variables z 1 and z 2 are defined as follows: z 1 = kψ 1 and z 2 = lψ 2 and symbols H i denotes functions

H i = H i (k, l, z 1 , z 2 ), i = 1, 4.
Let us suppose that the Hamiltonian system has a steady state P * with coordinates P * = (k * , l * , z * 1 , z * 2 ). In this case conjugate coordinates z * 1 and z * 2 of the steady state can be found from the first two equations of the Hamiltonian system, namely

z * 1 = k * f (k * , l * ) -(δ + ρ)k * , z * 2 = l * bf (k * , l * ) -ρl * . (10)
Further, we construct the linearized Hamiltonian system in a neighborhood of the steady state. Let symbol A = {α ij } 4 i,j=1 denotes the matrix of the linearized Hamiltonian system, where

α i1 = ∂H i (P * ) ∂k , α i2 = ∂H i (P * ) ∂l , α i3 = ∂H i (P * ) ∂z 1 , α i4 = ∂H i (P * ) ∂z 2 , i = 1, 4.

Nonlinear stabilizer

A nonlinear stabilizer is constructed under the following assumptions A 1 . It is assumed that matrix A has two real negative λ 1 and λ 2 and two real positive λ 3 and λ 4 eigenvalues.

This assumption means that the steady state P * = (k * , l * , z * 1 , z * 2 ) has the saddle character.

Let the symbols h i = {h ij } 4 j=1 , i = 1, 4 denote eigenvectors corresponding to eigenvalues λ i , i = 1, 4, respectively.

A 2 . It is supposed that first two coordinates of eigenvectors h 1 and h 2 corresponding to negative eigenvalues λ 1 and λ 2 meet the restriction h 11 h 22 ̸ = h 12 h 21 .

Construction of nonlinear stabilizer.

Idea of construction of the nonlinear stabilizer is based on results of the qualitative theory of differential equations (see [START_REF] Ph | Ordinary Differential Equations[END_REF]). Namely, the trajectory of the nonlinear Hamiltonian dynamics converges to the steady state tangentially to the plane generated by eigenvectors corresponding to negative eigenvalues of the Jacobi matrix. Let us describe the algorithm of construction of the nonlinear stabilizer. 1. To build the plane π generated by two eigenvectors h 1 and h 2 corresponding to two negative eigenvalues λ 1 and λ 2 , so that the steady state P * belongs to this plane π. 2. To extract conjugate variables z 1 and z 2 from equations of the plane. 3. To substitute the obtained relations of extraction instead of conjugate variables into control functions s 0 (t) and r 0 (t) corresponding to domain D 22 .

As a result, the algorithm provides construction of control s(t) and r(t) which is called nonlinear stabilizer. Let us consider each step in details.

Plane construction. Any vector v located in the plane π can be expressed through eigenvectors h 1 and h 2 in the following way: v = ν 1 h 1 +ν 2 h 2 . Hence, the plane π generated by two eigenvectors h 1 and h 2 and containing the equilibrium point P * can be written as follows:

k -k * = ν 1 h 11 + ν 2 h 21 , l -l * = ν 1 h 12 + ν 2 h 22 , z 1 -z * 1 = ν 1 h 13 + ν 2 h 23 , z 2 -z * 2 = ν 1 h 14 + ν 2 h 24 . (11)
Due to assumption A 2 coefficients ν 1 and ν 2 can be found from the first two equations [START_REF] Solow | Growth Theory: An Exposition[END_REF].

Extraction of conjugate variables. Conjugate variables z 1 and z 2 can be extracted from the equations [START_REF] Solow | Growth Theory: An Exposition[END_REF] of the plane π. As a result, we obtain

z 1 = z 1 (k, l) = z * 1 + γ 11 (k -k * ) + γ 12 (l -l * ), ( 12
)
z 2 = z 2 (k, l) = z * 2 + γ 21 (k -k * ) + γ 22 (l -l * ),
where γ 11 = - It should be mentioned that the following equalities take place

h
z 1 (k * , l * ) = z * 1 , z 2 (k * , l * ) = z * 2 . ( 13
)
Nonlinear stabilizer. The only thing left is to substitute expressions [START_REF] Tarasyev | Optimal Dynamics of Innovation in Models of Economic Growth[END_REF] into relations [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF] for optimal controls in the domain D 22 . Finally, we get the following structure of the nonlinear stabilizer:

s(k, l) = 1 - k z 1 (k, l)f (k, l) , r(k, l) = 1 - l bz 2 (k, l)f (k, l) . (14)
Substituting expressions for conjugate variables [START_REF] Tarasyev | Optimal Dynamics of Innovation in Models of Economic Growth[END_REF] to the first two equations of the Hamiltonian system (9) we get the stabilized Hamiltonian system:

k = f (k, l) -(δ + ρ)k - k z 1 (k, l) , l = bf (k, l) -ρl - l z 2 (k, l) . ( 15
)
Properties of the nonlinear stabilizer. Let us indicate main properties of the constructed nonlinear stabilizer.

Proposition 4. The nonlinear stabilizer (14) generates the nonlinear system (15) having the steady state with coordinates (k * , l * ) which are the same as the first two coordinates at the steady state of the original Hamiltonian system [START_REF] Sanderson | The SEDIM Model: Version 0.1[END_REF].

Proof of this propositions is based on the property (13) of the representation of adjoint variables z 1 = z 1 (k, l), z 2 = z ( k, l) in the plane π and relations [START_REF] Shell | Applications of Pontryagin's Maximum Principle to Economics[END_REF] for conjugate coordinates z * 1 , z * 2 of the steady state. Let us consider the linearized Hamiltonian system with Jacobi matrix A projected on subspace π. We substitute representation z 1 = z 1 (k, l) and z 2 = z 2 (k, l) (12) of conjugate variables into the first two equations of the linearized dynamics and collect similar terms

k = a 11 (k -k * ) + a 12 (l -l * ), l = a 21 (k -k * ) + a 22 (l -l * ), ( 16
)
where a ij = α ij + α 13 γ 1j + α 14 γ 2j , i, j = 1, 2.

Proposition 5. The matrix of the linearized stabilized system is the same as the matrix A = {a ij } 2 i,j=1 of the linearized Hamiltonian system projected on plane π.

In order to prove this proposition it is necessary to linearized stabilized Hamiltonian system (15) at the steady state (k * , l * ) neighborhood.

Let the symbol A denote the matrix of the linearized stabilized system (16). The next important question deals with eigenvalues of the stabilized Hamiltonian system (15). Proposition 6. The linearized stabilized Hamiltonian system (15) has two real negative eigenvalues coinciding with eigenvalues λ 1 and λ 2 and the following eigenvectors

h 1 = (h 11 , h 12 ), h 2 = (h 21 , h 22 ). ( 17
)
Proof. Basing on property 5 one can assert that the linearized stabilized Hamiltonian system coincides with the linearized Hamiltonian system (16) projected on plane π. For the Jacobi matrix A evaluated at the steady state P * the following equalities are fulfilled

A h i = λ i h i , i = 1,
) = = ( α 11 h i1 + α 12 h i2 + α 13 h i3 + α 14 h i4 α 21 h i1 + α 22 h i2 + α 23 h i3 + α 24 h i4 ) = ( λ i h i1 λ i h i2 ) = λ i h i , i = 1, 2.

⊓ ⊔

The following theorem collects all obtained results. The proof of the theorem follows directly from properties of the nonlinear stabilizer.

Remark 1. The constructed nonlinear stabilizer generates the nonlinear system which is closed with respect to phase variables. The solution of the obtained stabilized system approximates optimal trajectories of the original Hamiltonian system in a neighborhood of the steady state, since a trajectory of the nonlinear Hamiltonian system that tends to the equilibrium point is tangent to the plane formed by two eigenvectors corresponding to negative eigenvalues. One can use this fact to estimate the growth rates of optimal trajectories The growth rates are determined by values of negative eigenvalues. Trajectories of the system (15) generated by the nonlinear stabilizer (14) and the original Hamilton system (9) are calculated numerically by the Runge-Kutta method. Figure 1.(a) demonstrates phase trajectories k(l) as a solutions of the stabilized (15) and Hamiltonian (9) dynamics. One can see that these trajectories almost coincide with each other especially at the vicinity of the steady state. Optimal trajectories of the capital stock k(t) and labor efficiency l(t) and its stabilized solutions are depicted at figures 1.(b) and 1.(c) respectively. In the steady state neighborhood optimal trajectories are very close to its stabilized solutions.

  Numerical simulations. The calculations are carried out on the basis of the data on the US economy in the period of 1900 to 2005. The data values are normalized with respect to the data values of 1900. The production function of the Cobb-Douglas type is used: f (k, l) = µ k α l β . The calibration procedure for the model parameters provides the following values: µ = 2.19942, α = 0.31, β = 0.09, λ = 0.03, δ = 0.2, ρ = 0.013, b = 0.31, a s = 0.3, a r = 0.2, k 0 = 1, l 0 = 1. The Hamiltonian system has the equilibrium point P * with coordinates k * = 5.75, l * = 5.2, z * 1 = 1.8188, and z * 2 = 2.9684. The control parameters at the equilibrium point take the values s * = 27.95 and r * = 3.79. All four eigenvalues of the matrix A calculated at the equilibrium point P * = (k * , l * , z * 1 , z * 2 ) are real numbers; two of them are positive, and the other two are negative: λ 1 = -0.268, λ 2 = -0.094, λ 3 = 0.124, λ 4 = 0.298.

Fig. 1 :

 1 Fig. 1: Stabilized and optimal graphs (a) of phase trajectories, k(l); (b) of the capital stock, k(t); (c) of the labor efficiency, l(t).

  4. Moreover eigenvectors h 1 , h 2 are located at the plane π. Thus, for coordinates of these vectors are valid relationsh i3 = γ 11 h i1 + γ 12 h i2 , h i4 = γ 21 h i1 + γ 22 h i2 , i = 1, 2.Using these facts let us check the following equalitiesI h i = λ i h i , i = 1, 2. i1 + α 12 h i2 + α 13 (γ 11 h i1 + γ 12 h i2 ) +α 14 (γ 21 h i1 + γ 22 h i2 ) α 21 h i1 + α 22 h i2 + α 23 (γ 11 h i1 + γ 12 h i2 ) + α 24 (γ 21 h i1 + γ 22 h i2 )

	=	(	α 11 h
	A h i =	(	a 11 a 12 a 21 a 22	) (	h i1 h i2	)	=

  Theorem 1. Under assumptions A 1 and A 2 for the Hamiltonian system[START_REF] Sanderson | The SEDIM Model: Version 0.1[END_REF] constructed in domain D 22 and linearized in a neighborhood of the steady state P * the nonlinear stabilizer (14) exists and generates the nonlinear dynamical system (15) which is closed with respect to the phase variables k, l and has the following properties 1. the steady state of the closed system (15) has coordinates (k * , l * ) coinciding with the phase coordinates k and l of the steady state P

* of the original Hamiltonian system (9); 2. the system (15) is stabilized at the steady state P * ; 3. the eigenvectors h 1 and h 2 of the linearized closed system (16) generated by the nonlinear stabilizer are evaluated by formulas (17).

⋆ The research is supported by the Russian Fund of Basic Research (Grants 11-01-0042-a, 11-01-12088-ofi-m-2011, 11-01-12112-ofi-m-2011), by the Program for the Sponsorship of Leading Scientific Schools (Grant NSCH-64508.2010.1), by the Program of the Presidium of RAS "Dynamic Systems and Control Theory", by the Program of the Presidium of RAS No. 38Π (Project 12-Π-1-1038), by the Project of the Ural Branch of RAS "Socio-Economic Development of Regions: Forecasting and Optimal Control" (Grant 12-Π-7-1001), and the International Institute for Applied System Analysis (IIASA).