ON ROBUST ON-LINE PARAMETER IDENTIFICATION TECHNIQUES IN ENVIRONMENTAL MODELING

14th IFAC Symposium on System Identification, Newcastle, Australia, 2006

ON ROBUST ON-LINE PARAMETER IDENTIFICATION TECHNIQUES IN ENVIRONMENTAL MODELING Vyacheslav Maksimov 1

Institute of Mathematics and Mechanics, Ekaterinburg, Russia

Abstract: In our paper we would like to discuss about one method of auxiliary controlled models and the application of this method to solving some problems of identification for differential equations. Problems of determining some parameters through the information on equation’s solutions are often called reconstruction (identification) problems. Therewith it is assumed that the input information (results of measurements of current phase states of a dynamical system) forthcomes in the process. As to unknown parameters, they should be reconstructed in the process too. One of the methods of solving similar problems was suggested in (Kryazhimskii and Osipov, 1983; Kryazhimskii and Osipov, 1984; Osipov and Kryazhimskii, 1995; Maksimov, 1995; Maksimov, 2002a; Maksimov, 2002b). This method based on the ideas of the theory of ill-posed problems actually reduces an identification problem to a control problem for an auxiliary dynamical system—the model (Krasovskii and Subbotin, 1988). Regularization of the problem under consideration is locally realized during the process of choosing a positional control in the system-model. The method mentioned above was applied to a number of problems described by some classes of ordinary differential equations as well as by equations with distributed parameters. Different system’s characteristics varying in time were under reconstruction, for example, unknown discontinuous inputs, initial and boundary data, distributed disturbances, coefficients of an elliptic operator and so on. In the present paper, we illustrate this method with a model example. Keywords: Robust identification, Environmental modeling

ment of complex industrial projects. In recent years, a considerable number of mathematical models describing the process of interaction between the climate and human activities has been designed.

1. THE MODEL The modern period in history is characterized by a rapidly growing complexity in the interrelations between the human activities and the processes in the environment. Large-scale positive and negative impacts of various human-driven processes and their growing dynamism have led to the necessity of the active and systematic treatment of actual and virtual responses in the environment. Studying different aspects of ecological modeling has become a critical stage in manage-

One of such models is a dynamical model connecting main economic and climatic indices suggested in (Nordhaus, 1994). This model is oriented to developing an economic strategy directed to deceleration of global warming. The main goal of the analysis of the model is to provide the means for tackling the following question: whether the reduction of emissions of greenhouse gases is justified from the economical viewpoint or not. The model takes into account global processes: it is assumed that the structure of economy is the same for all countries; the climate change is characterized by the average value of the temperature

1

Work was supported in part by RFBR (project 06-01-00359), by the International Institute for Applied Systems Analysis, by the Program for Support of Leading Scientific Schools of Russia, by the Ural-Siberian Integration Project, and by the Program on Basic Research of Presidium of the Russian Acad. Sci. (project “Catastrophe risk estimation: Theoretical research”).

744

T˙ (t) = c1 T (t) + c2 T1 (t) + c3 F (t) T˙1 (t) = c4 (T (t) − T1 (t)) M˙ 1 (t) = βE(t) − δM M1 (t)

on Earth’s surface and so on. This model contains three types of parameters. 1) Constant parameters (their list is presented in tables 2.3 and 2.4 on page 21 of work (Nordhaus, 1994)). 2) Functions that are considered (for simplicity of the analysis) as exogenous with respect to the model and are a priori given. 3) Inner functions that are connected to one another and to exogenous parameters by means of some algebraic and differential equations. The list of these functions is presented in table 2.3. (see (Nordhaus, 1994)), and the model equations are presented in table 2.2. Let us give the list of functions: μ(t) is the rate of emissions reduction with respect to uncontrollable emissions; E(t) is emissions of greenhouse gases GHGs (CO2 (carbonic acid gas) and chlorine-fluorine carbons only); M1 (t) = (M (t) − 590) is the excess of mass of GHGs in the atmosphere in comparison with the pre-industrial period; T (t) is the average atmospheric temperature (on Earth’s surface); T1 (t) is the average deep-ocean temperature; I(t) is the gross investment; K(t) is the capital stock; F (t) is the atmospheric radiative forcing from GHGs; O(t) is the forcing of exogenous GHGs (i.e., of gases, which are considered as uncontrollable; there are all GHGs, except CO2 (carbonic acid gas) and chlorine-fluorine carbons); A(t) is the level of technology; σ(t) is the ratio of GHGs emissions to global output; L(t) is the population at time t, also equal to labor inputs, Ω(t) is the output scaling factor due to emissions control and to damages from climate change, Q(t) is the gross world product.

˙ K(t) = −δK K(t) + I(t),

(1)

t ∈ [0, ϑ],

where t is time, ϑ is a terminal time moment,   M1 (t) + O(t), F (t) = 4, 1 · log2 1 + 590 E(t) = (1 − μ(t))σ(t)Q(t), Q(t) =

(1 − b1 μ(t)b2 ) . (1 + θ1 T (t)θ2 )A(t)K(t)γ L(t)1−γ

An initial state of Σ, {T (0), T1 (0), M1 (0), K(0))}, is assumed to be known and a priori given. Functions μ(t) and I(t) are considered as control parameters determining a strategy of global control of climate and economy. Our attention to the rather simplified model described in (Nordhaus, 1994) is not casual. It has been actively investigated in the recent years. In particular, it served as the basis for creation of the Dynamic Integrated model of Climate and the Economy (the DICE model). The latter is still under improvement and is used by some specialists for analyzing implications of joining the Kyoto Protocol by different countries. 2. DYNAMICAL INVERSE PROBLEM At that, the direct problem is solved, namely, possible strategies (rules of forming μ(t) and I(t)) are specified, and system’s dynamics is computed. The comparative analysis of simulation results for different structures is performed. In addition, the analysis of sensitivity of results with respect to some model parameters is fulfilled. Our aim differs from that of (Nordhaus, 1994). We dwell on an inverse problem. It consists in the following. Let us assume that some function I(t) is known. The latter implies that function K(t) is also known. Neglecting small values (b1 = 0,0686, ϑ1 = 0,00144), we transform the system (1) to the form

Schematically, the connections between the inner functions can be pictured in the following way:

T˙ (t) = c1 T (t) +c2 T1 (t) +  M1 (t) c5 · log2 1 + + c3 O(t), 590 T˙1 (t) = c4 (T (t) − T1 (t)), (2) ˙ M1 (t) = E1 (t)(1 − μ(t)) − δM M1 (t),

T ∗∗ −→ Ω ↑ ↑  F μ Q ↑ ↓  ↑  K ∗ ←− I M1∗ ←− E

˙ K(t) = −δK K(t) + I(t),

t ∈ [0, ϑ],

Here the functions marked by asterisk are solutions of linear differential equations of the first order; the function T (t) is a solution of a linear differential equation of the second order.

where E1 (t) = βσ(t)A(t)K(t)γ L(t)1−γ . Hereinafter, we consider the system Σ of the form (2).

If we pass from the discrete model suggested by the authors to the “continuous” one, then the equations of the model Σ take the form:

The problem under consideration may be formulated in the following way. At frequent enough time moments τi ∈ Δ = {τi }m i=0 , τi+1 = τi + δ, τ0 = 0,

745

τm = ϑ, values of the average atmospheric temperature T (τi ) and the average deep-ocean temperature T1 (τi ) are inaccurately measured. Results of measureh h , ξ2i } ∈ R2 ) satisfy the ments (vectors ξih = {ξ1i inequalities h 2 h 2 | + |T1 (τi ) − ξ2i | ≤ h2 , |T (τi ) − ξ1i

Dynamical inverse problem. It is required to indicate differential equations of the model M w˙ h (t) = f1 (ξih , wh (τi ), uhi ), t ∈ δh,i = [τi,h , τi+1,h ),

(3)

wh (0) = w0h ,

where h ∈ (0, 1) is a level of informational noise. It is required to design an algorithm allowing us to reconstruct (synchro with the process) the unknown excess of mass of GHGs M1 (t) and rate of emissions reduction μ(t). This is the meaningful formulation of dynamical reconstruction problems being investigated in the present paper.

- U 6

- M

wh (t) ∈ R4 ,

and the rule of choosing controls uhi at the moments τi being a mapping of the form U : {τi , ξih , wh (τi )} → uhi = {uhi1 , uhi2 } ∈ R2 (5) such that the convergences ϑ

The scheme of algorithms for solving such problems is given in Figure 1. uh (t)

(4)

τi = τi,h ,

0

wh (t) -

|uh1 (t) − M1 (t)|2 dt → 0,

ϑ

(6) |uh2 (t) − μ(t)|2 dt → 0

0 h

ξ (t)

Σ 

μ(t)

take place whereas h tends to 0. Here (see Figure 1) uh (t) = {uh1 (t), uh2 (t)}, uh1 (t) = uhi1 , uh2 (t) = uhi2 for t ∈ δh,i .

Fig. 1. According to this scheme (see also (Kryazhimskii and Osipov, 1983; Kryazhimskii and Osipov, 1984; Osipov and Kryazhimskii, 1995; Maksimov, 1995; Maksimov, 2002a; Maksimov, 2002b)), an auxiliary dynamical system M (the model) is introduced. This model functioning on the time interval [0, ϑ] has an input uh (t) and an output w h (t). The process of synchronous feedback control of the systems Σ and M is organized on the interval [0, ϑ]. This process is decomposed into (m − 1) identical steps. At the i-th step carried out during the time interval δi = [τi , τi+1 ) the following actions are fulfilled. First, at the time moment τi according to the chosen rule U the functions

3. THE ALGORITHM FOR SOLVING THE INVERSE PROBLEM Let us turn to the description of the algorithm for solving the inverse problem. From the above, it is necessary to indicate the model (4) and the strategy U (5) providing the convergences (6). Let a restriction on the rate of emissions reduction be known, i.e., we know a number f > 0 such that |μ(t)| ≤ f

are calculated by measurements ξih and wh (τi ). Then (till the moment τi+1 ) the control

max xμ (t) ≤ K,

0≤t≤ϑ

sup x˙ μ (t) ≤ K,

τi ≤ t < τi+1 ,

0≤t≤ϑ

is fed onto the input of the model M . The values h and wh (τi+1 ) are the results of the work of the ξi+1 algorithm at the i-th step.

τ0,h = 0,

M1 (t) ∈ [a1 , a2 ].

(7)

Here xμ is the Euclidean norm of the vector xμ . Let the following condition be fulfilled.

Thus, the inverse problem may be formulated as follows. In the sequel, a family of partitions h Δh = {τi,h }m h=0 ,

t ∈ [0, ϑ].

From now on, it is assumed that we know numbers K, a1 , a2 ∈ (0, +∞), a1 < a2 , such that each solution xμ (t), xμ (t) = {T (t), T1 (t), M1 (t), K(t)}, of the equation (2) satisfies the following conditions

uh (t) = uhi = U (τi , ξih , wh (τi )),

u = uh (t),

for all

Condition 1. Functions σ(t), A(t) and L(t) are con˙ tinuous. Function Q(t) is differentiable and |Q(t)| ≤ c < +∞.

τi+1,h = τi,h + δ(h), τmh ,h = ϑ,

Introduce some function α(h): (0, 1) → R+ = {r ∈ R : r ≥ 0} with the properties:

of the interval [0, ϑ] with diameters δ = δ(h) is assumed to be fixed.

α(h) → 0, (h1/6 + ω(h))α−1 (h) → 0 as h → 0.

746

Estimate the variation of the value

Here ω(h) = ωE (δ) + δ; ωE (δ) is the modulo of continuity of function E1 (t). It is easily seen that the modulo of continuity ωM (δ) of function M1 (t) equals Cδ, where C is some constant that can be written in an explicit form. The function α = α(h) plays the role of a parameter of regularization (Tikhonov and Arsenin, 1977). Let the model (4) have the form

ε(t) = |w2h (t) − M1 (t)|2 + t α(h) {|uh2 (τ )|2 − |μ(τ )|2 } dτ. 0

h h + c2 ξ2i + w˙ h (t) = c1 ξ1i

It is easily seen that for t ∈ δi = [τi , τi+1 ) the following inequality is true:

uh 4,1c3 log2 (1 + 1i ) + c3 Q(τi ) 590 h h (8) w˙ 1h (t) = c4 (ξ1i − ξ2i )

t ε(t) ≤ ε(τi ) + δ(h)

w˙ 2h (t) = E1 (τi )(1 − uh2i ) − δM uh1i t ∈ δi = [τi , τi+1 ),

τi

t

w˙ 3h (t) = −δK K(τi ) + I(τi ),

t

μi (τ ) dτ + α(h)

τi = τi,h , τi

and the rule U of forming the control uhi = {uhi1 , uhi2 } is as follows h

uh1i = 590(2πi − 1),

uh2i =

⎧ (1) ⎨ βi α−1 (h), if |βi(1) | ≤ α(h)f ⎩

(1)

Here

Consider the value μi (t). We have for t ∈ δi

(9)

μi (t) =

(10)

t ∈ δi ,

(13)

(2)

w1h (0) = T1 (0),

(E1 (τi ) − E1 (t)) ,  (2)

λ2i (t) = 2βi μ(t) − uh2i E1 (τi ), λ1i (t) = 2βi

w3h (0) = K(0).

(2)

(E1 (t) − E1 (τi )) μ(t),  (2)

λ4i (t) = 2δM βi M1 (t) − uh1i ,

λ3i (t) = 2βi

(2)

βi

= w2h (τi ) − M1 (τi ).

Estimate each term in the right-hand part of the equality (13). From (7) and (11) it follows that

b2 = log2 (1 + a2 /590),

λ1i (t) ≤ C1 ω(δ),

= E1 (τi )(w2h (τi ) − uh1i ),

h ). ρi = −8,2c3 (wh (τi ) − ξi1

m h −1

Let the family of partitions Δh have the diameters

i=0

δ(h) ≤ h.

ϑ C3

Then the following theorem is true.

t ∈ δi .

τi+1

λ4i (t) dt ≤ τi

|M1 (t) − uh1 (t)| dt ≤ C4 h1/6 .

(14)

In addition, the estimate m h −1 i=0

Proof. It is easily seen that from results of [2] the following inequality follows

m h −1 i=0

|uh1 (t) − M1 (t)|2 dt ≤ Ch1/3 .

λ3i (t) ≤ C2 ω(δ),

0

Theorem 1. Let E1 (t) > 0, t ∈ [0, ϑ]. Then the convergences (6) take place under choosing the model equation in the form of (4), (8) and the strategy U in the form of (5), (9), (10).

ϑ

λji (t),

where

b1 = log2 (1 + a1 /590), (1)

4 j=1

⎧ ⎪ ρ /(2h2/3 ), if ρi (2h2/3 )−1 ∈ [b1 , b2 ] ⎪ ⎨ i πih = b1 , if ρi (2h2/3 )−1 < b1 ⎪ ⎪ ⎩ b2 , if ρi (2h2/3 )−1 > b2 ,

βi

{|uhi2 |2 − |μ(τ )|2 } dτ,

μi (t) = 2(w2h (τi ) − M1 (τi ))(w˙ 2h (t) − M˙ 1 (t)).

The initial state of the model is

w2h (0) = M1 (0),

(12)

τi

f sign βi , otherwise.

wh (0) = T (0),

|w˙ 2h (τ ) − M˙ 1 (τ )|2 dτ +

(11)

τi+1

τi+1

λ6i (t) dt ≤ τi

(h)

Ci |u1i − M1 (τi )| dt ≤ C5 (h1/6 + δ),

τi (h)

|Ci | ≤ C5 ∀i ∈ [0 : mh − 1], h ∈ (0, 1),

0

747

is valid. Then we have

The initial conditions for the system are as follows:

λ2i (t) ≤ λ5i (t) + λ6i (t), λ5i (t) =

(1) 2βi (μ(t)

−

uh2i ),

T (0) = 1,

T1 (0) = 0.5,

M1 (0) = 50,

K(0) = 10.

λ6i (t) = 2E1 (τi )(uh1i − M1 (τi ))(μ(t) − uh2i ). The upper panel shows the samples of M1 (thick line) and uh1 (thin line); the samples of μ (thick line) and its estimate uh1 (thin line) are in the plot below.

Note that the control uh2i of the form (10) is calculated by the following formula: (1)

uh2i = arg min{−2βi u + αu2 : |u| ≤ f }.

The equation was solved using the Euler method with step δ. The results of numerical experiments show that the mean-square convergence (6) takes place under “reduction” of parameters h and δ or of one of them. This is seen, for example, from Fig. 2, 3 that correspond to the fixed h = 10−3 and different δ. The presence of “jumps” of the uh2 (t) curve in Fig. 2 matches the fact that uh2 (t) converges to μ(t) in L2 [0, 2] metric (see (6)), not in the continuous metric.

Therefore, in virtue of (3) we obtain τi+1

{|uhi2 |2 − |μ(τ )|2 } dτ =

λ5i (τi+1 ) + α(h) τi τi+1



 (1) −2βi uhi2 + α(h)|uhi2 |2 −

(15)

τi



(1) −2βi μ(τ ) + α(h)|μ(τ )|2 dτ ≤ 0. Taking into account (12)–(15) and the inequality δ(h) ≤ h, we have for all i ∈ [1 : mh ] the following estimate ε(τi ) ≤ C(h1/6 + ω(δ(h))) ≤ C(h1/6 + ω(h)). Further argument corresponds to the standard scheme (see, for example, (Kryazhimskii and Osipov, 1983; Osipov and Kryazhimskii, 1995; Maksimov, 2002b)). The theorem is proved.

4. RESULTS OF COMPUTER MODELING The algorithm described above was tested on computers. System (2) evolves on the interval [0, 2] and it is subject to the external function μ(t) which is equal to 1 + 0.5t. Our point of view is that we do not know μ(t) and we want to recover it from a finite number of corrupted samples of the evolution of {T, T1 }.

Fig. 2. δ = 0.001

We discuss now the result of simulation, see Fig. 2, 3. In these figures, the results of computer modeling of the dynamic inverse problem are presented for the following case (see (1)): c1 = c2 = c3 = 1,

σ = 1 + 0.5t,

c4 = 0.5,

O(t) = 5 sin(t),

δK = 0.65,

L(t) = 1,

δM = 0.0833,

μ(t) = 1 + 0.5t,

β = 0.1,

I(t) = 1 = 0.15t2 ,

γ = 0.1,

A(t) = 2t1/2 .

The parameters of the algorithm are the following: f = 2,

a1 = 30,

Fig. 3. δ = 0.0001

a2 = 60.

748

5. REFERENCES Krasovskii, N.N. and A.I. Subbotin (1988). GameTheoretical Control Problems. Springer Verlag. Kryazhimskii, A.V. and Yu.S. Osipov (1983). On modeling of control in a dynamical system. Izv. Akad. Nauk USSR. Tech. Cybern. (2), 51–60. In Russian. Kryazhimskii, A.V. and Yu.S. Osipov (1984). On positional calculation of normal controls in dynamical systems. Probl. Control and Inform. Theory 13(6), 425–436. Maksimov, V.I. (1995). On the reconstruction of a control through results of observations. In: Proceedings of the Third European Control Conference. Rome, Italy. pp. 3766–3771. Maksimov, V.I. (2002a). Dynamical estimation of an input in nonlinear differential systems. In: Proc. of the 15th IFAC World Congress on Automatic Control. Vol. D. Barcelona, Spain, July 21–26, 2002. pp. 199–203. Maksimov, V.I. (2002b). Dynamical inverse problems for distributed systems. VSP. Nordhaus, W.D. (1994). Managing the global commons. The economics of climate change. MIT Press. Osipov, Yu.S. and A.V. Kryazhimskii (1995). Inverse problems of ordinary differential equations: dynamical solutions. Gordon and Breach. Tikhonov, A.N. and V. Arsenin (1977). Solutions of ill-posed problems. Wiley. New York.

749