METHOD FOR FINDING A SOLUTION TO A LINEAR NONSTATIONARY INTERVAL ODE SYSTEM Текст научной статьи по специальности «Математика»

Вестник СПбГУ. Прикладная математика. Информатика... 2021. Т. 17. Вып. 2 UDC 517.926.7+517.977.58 MSC 34В05, 49М05

Method for finding a solution

to a linear nonstationary interval ODE system*

A. V. Fominyh

St. Petersburg State University, 7-9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation

Institute of Problems in Mechanical Engineering, Russian Academy of Sciences, 61, Bolshoy pr. V. O., St. Petersburg, 199178, Russian Federation

For citation: Fominyh A. V. Method for finding a solution to a linear nonstationary interval ODE system. Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2021, vol. 17, iss. 2, pp. 148-165. https://doi.org/10.21638/11701/spbul0.2021.205

The article analyses a linear nonstationary interval system of ordinary differential equations so that the elements of the matrix of the system are the intervals with the known lower and upper bounds. The system is defined on the known finite time interval. It is required to construct a trajectory, which brings this system from the given initial position to the given final state. The original problem is reduced to finding a solution of the differential inclusion of a special form with the fixed right endpoint. With the help of support functions, this problem is reduced to minimizing a functional in the space of piecewise continuous functions. Under a natural additional assumption, this functional is Gateaux differentiable. For the functional, Gateaux gradient is found, necessary and sufficient conditions for the minimum are obtained. On the basis of these conditions, the method of the steepest descent is applied to the original problem. Some numerical examples illustrate the constructed algorithm realization.

Keywords: linear nonstationary interval system of ordinary differential equations, differential inclusion, support function, the steepest descent method.

1. Introduction. Finding a solution of interval systems of ordinary differential equations is an important problem. This is because in many real processes the object is governed by the system in the presence of uncertainty. It can be of a different nature: structural parameters of the system, errors, material properties, perturbations. However, the change interval of these uncertainties is usually known heuristically. Such interval systems are considered in this paper.

One basic approach, which can be effectively applied to the problem considered, is using standard interval methods for reachable set estimation. For linear models an enclosure of the exact solution set is obtained (see, for example, [1^3]). Let us also note also some works devoted to estimating attainability sets in nonlinear case [4-7].

One problem, which is closely connected with the considered in this paper, is numerical solutions of fuzzy differential equations. Since in the work [8] the authors introduced the notion of fuzzy numbers and defined the basic arithmetic operations with them and Chang and Zadeh [9] introduced the fuzzy derivative, this mathematical field began to develop rapidly. Elementary fuzzy calculus was constructed by the authors in [10]. Later fuzzy calculus was applied to the differential equations. Kandel and Byatt [11] applied

* The main results of this paper (Sections 3 and 5) were obtained in IPME RAS and supported by Russian Science Foundation (project N 20-71-10032). © St. Petersburg State University, 2021

the concept of fuzzy differential equation to the analysis of fuzzy dynamical problems. The important Cauchy problem for fuzzy differential equations was considered by many-authors, for example, Kaleva [12]. One of the main concepts of fuzzy calculus is Hukuhara derivative [13].

A scheme based on the classical Euler method was used in the article [14] for numerical solution of the fuzzy systems. The papers [15-17] also used Euler method and Runge^ Kutta scheme for such problems. Solving numerically the fuzzy differential equation by the Taylor method is considered in the following paper [18]. In [19] modified Taylor series method is considered. In [20] the approximate solution of fuzzy differential equations system with the homotopy analysis method, introduced in [21], is found. The work [22] seeks for the "minimal" solution of the fuzzy system, based on the pseudo inverse. All the mentioned papers deal with the known right-hand side, but with the interval initial condition. The paper [23] proposes a new interval analysis method for the dynamic response of nonlinear systems with uncertain-but-bounded parameters using Chebyshev polynomial series. This paper deals with the interval right-hand side and shows on examples, that as a rule Chebyshev polynomial series are more effective, than Taylor polynomial series.

Although there are many papers considering the systems with some uncertainties, much less works are devoted to the solution of boundary value problems. The paper [24] considers the problem of the existence of solutions to the two-point boundary value problem. The tolerable solution set for interval linear systems is also estimated in [25]. The work [26] considers questions connected with the construction of two-sided estimates for solution sets of systems of nonlinear differential equations with interval parameters. In work [27] numerical methods for constructing the boundaries of the solution set of boundary-value problems for ordinary differential equations (ODE) with interval coefficients are investigated. The paper [28] develops a new solution method, which constructs a solution as a fuzzy set of real vector-functions (the crisp coefficients of the system are considered, while forcing functions and initial conditions are fuzzy). The papers [29] and [30] present a method for modelling and simulation of uncertain dynamical systems, using differential inclusions and fuzzy arithmetic transformation respectively. Real paper is devoted to the boundary value problem with the linear nonstationary nonhomogeneous right-hand side. The initial problem is reduced to the boundary value solution of a differential inclusion, which is sought for in the same fashion as in the work [31]. See some advantages of the proposed scheme in Remarks 1 and 7 below.

Consider the linear nonstationary interval system of ODE

x = A(t)x + g(t), t e [0,T],

(1)

where

A

( [ân an] [a21 Ô21]

ai2] Â22]

[aln ai„] ^

• [a2n a2n]

\ [anl a"i] [Sn2 a"2] • • • [ann ann] J

( \

(2)

with the boundary conditions

and

x(0) = xo x(T ) = xt ■

(3)

(4)

In formula (1) T is the given finite moment of time, x is a n-dimensional continuous vector-function of the phase coordinates with piecewise continuous and bounded derivative in the interval [0, T], In formula (2) «¿ -(i) and «¿¿(i), «¿¿(t) s' «¿¿(i) Vt G [0, T], i,j = l,n, are given continuous functions and g (t) and <?j(i), g (t) < ^¿(t) Vt G [0, T], i = l,n, are given continuous functions. Further, for brevity, we will write aLipaij, g.,gi instead of Sjj(i), a>ij(t), g^t), ~§i{t) respectively, i, j = l,n, t G [0, T]. In formulas (3), (4) xq,xt G 1" are the given vectors.

It is required to find the vector-function x* G Cn[0, T], which is a solution of system (1) and satisfies conditions (3), (4). We assume, that there exists such a solution.

Here Cn [0, T] is the space of n-dimensional vector-functions, continuous in [0,T], with derivative from the space Pn[0, T]; Pn[0, T] is the space of n-dimensional vector-functions, piecewise continuous and bounded in [0, T], In the sequel, we also need the space ^[0, T] of square-summable in [0, T] n-dimensional vector-functions.

If t0 G [0, T) is the point of discontinuity of the vector-function X, then we suppose that X(t0) is a right-handed derivative of the vector-function X at the point t0, X(T) is

xT

Remark 1. The stated problem may be considered as the following control problem. It is required to find the "controls" «¿¿(t) G P[0,T] (such as «¿ -(t) < «¿¿(i) < «¿¿(i) Vt G [0,T], i,j = T~n) and Vi(t) G P[0,T\ {git) < Vi(t) < Vt G [0,T], i = T^n), which bring system (1) (with the matrix A(t) consisting of the elements «¿¿(t), i,j = l,n, and the vector g(t) consisting of the elements Vi(t), i = 1, n) from initial point (3) to final state (4) in the time T. However, ^n such formulation it is required to find n2 + n control functions. On the other hand, from the construction of the method described in the article, one will see that if we write the system of this problem in form (1), (2) and apply the

n

of the considered system.

For the arbitrary set F C Rn define the support function (of the vector ^ G Rn) as

follows: c(F, = sup (f, where (a, b) is the inner product of the vectors a,b G Rn. f eF

Let us rewrite the system (1) in the form of the differential inclusion

X G F(x,t), t G [0,T], (5)

here

F(x, t) = A(t)X + g(t).

Apparently, the mapping F(x, t) is a convex compact set from Rn for every moment of time t G [0, T] and for every phase point x G Rn. One can also see that F(x; t) is continuous.

Now we can reformulate the original problem as follows. It is required to find the vector-function x* G Cn[0,T], which is a solution of differential inclusion (5) and satisfies conditions (3), (4).

Remark 2. Instead of trajectories from the space Cn [0,T] with derivatives from the space Pn[0,T] in the paper one can consider absolutely continuous in the interval [0, T] trajectories with measurable and almost everywhere bounded in [0, T] derivatives respectively (which must satisfy the differential inclusion almost everywhere on [0, T]). The choice of the space of solutions in the article is explained by the possibility of their practical construction. Note that under the assumptions made there exists [32] even a continuously differentiable ("classical") solution of the Cauchy problem (3), (5) at least in the vicinity xo

2. Reduction to the unconstrained minimization problem. Further, for brevity, we sometimes write F instead of F(x, t). Since for all t G [0, T] and for all x G 1" the multivalued mapping F(x, t) is a convex, closed and bounded set in 1", inclusion (5) may be rewritten as follows [33]:

x(t), 4) < c(F(x(t),t), 4) W G S, Vt G [0, T],

where S is the unit sphere in 1" with the center in the origin.

FF n-dimensional parallelepiped with the the center c(x) = (ci(x),..., cn(x))/:

j=i

and with the half-sides

Si~ 9,

' 2 "" 1*3-1 +

j=i

So the support function of the set F can be expressed [33] by the formula

n n

c(F (x,t)^) = YJ Ci(x)^i +Y; ri(*Wi\ ^ e (6)

i=1 i=1

We have for every tp e S (if Xi(t) ^ 0, i = 1, n)

= ¿ ^pi ^ + slgnXi (7)

i j=1 j=1

So we see that the support function of the set F is continuously differentiate in the phase coordinates x, if Xj ^ 0 V« = 1, n.

Denote z(t) = X(t), z e Pn[0,T], then from (3) we get

t

<t)= X0 + j z(t)dr.

Put

1(ф, z, t) = (z(t), ф) - c(F(x(t), t), ф), (8)

h(z,t) = ггтхтгАх{0,£(ф, z,t)} (9)

^es

and construct the functional

T

Ф) = \J h?{z{t),t)dt. (10)

0

Remark 3. The magnitude h(z(t),t) at the every fixed t G [0,T] is the Euclidean distance from the point z(t) to the set F(x(t),t), and functional (10) is the half-squared deviation (in the L"[0,T]-norm) of the trajectory z(t) from the set F(x,t).

Consider the set

Q = {z e Pn[0,T] | p(z) = 0}. It is not difficult to see that for functional (10) the following relations are valid:

v(z) = 0(z e Q), if (x(t),^) < c(F(x(t),tw e S, Vt e [0,T], (11) ¥>(z) > 0 (z e Q), otherwise, so inclusion (1) is satisfied, if ^(z) = 0.

Introduce the functional

X(z) = i + J z(t)dt -

One can see that boundary conditions (3), (4) are satisfied, if x(z) = 0. Construct the functional

I (z) = <p(z) + x(z). (12)

z*

t

x* (t) = xo + J z* (t)dr o

is the solution of the initial problem. So, the problem of finding an approximate solution of the original problem has been reduced to the minimization of functional (12) in the space Pn[0,T].

I

Suppose the coordinate Xi(t), i = 1, n, t G [0, T], vanishes only at the finite number of time

moments Uj G [0, T], i = 1 ,n. Denote such a set of trajectories X[0,T]. If x G X[0,T],

, , . . dc(F (x,t),^) . . . . , . ,

then the derivative --- exists and is continuous m (x, w.t) with the exception

ox _

of the finite number of points Uj, i = 1, n.

Suppose the optimal trajectory x* belongs to X [0, T].

Remark 4. The assumption made is presumably not burdensome, clS it IS violated, for example, in the case when there exists such a "coordinate" xj, i G {1, Idots, n}, which remains in the zero position (and has zero "speed") in some interval [ti, t2] C [0, T], t2 > t\. Such a case does not seem natural for control systems considered on the finite interval [0,T].

Let us obtain some crucial sufficient conditions, when this case for sure doesn't occur for n = 1 and n = 2 for homogeneous equations. Consider the case n = 1. We have the equation

x{t) = [aii(i) an(t)]x(t), x(0) = xo, t G [0, T],

Here the explicit solution is x(t) = x(0)exp^ u(t)d^j =0 Vt G [0,T], if x(0) = 0

(where u(t) G [an(i) ou(i)] Vt G [0,T]).

Consider the case n = 2. Suppose x1 (t) =0 Vt G [t1 ,t2] C [0,T^ t2 > t1. Then Vt G [t1,t2] we have

0 = [al2{t) a12(t)]x2(t), xi(0) = x0i, x2(t) = [a22(t) a22(t)]x2(t), x2(0) = x02.

If there exists such a point t* G [ti,¿2], that 0 ^ [a12(t*) ai2(t*)\, then from the first equation we conclude, that x2 (t*) = 0. The considered (second) linear homogeneous equation has the zero "initial" condition, so its solution vanishes on the whole interval: x2(t) =0 Vt G [t1,t2], But in this case we have xi(t1) = 0 x2(t1) = 0, so if we consider the original system on the whole interval [0, T] with this "initial" condition, we see, that it has the only solution x1(t) = 0 x2(t) =0 Vt G [0, T], what is impossible, if x21 + x01 = 0. If we argue in a completely analogous fashion in the following case x2(t) = 0 t G [t3,t4] C [0,T], t4 > t3 (suppose there exists such a point t* G [t3,t4], that 0 ^ [a21(t*) a2i(t*)]), then we will come to the analogous conclusion (that this is also impossible, if x01 + x01 = 0).

Let us study the differential properties of functional (12) for the support function of a more general form (later we will return to the specific support function given by formula (6), see Remark 5 below).

Assume that the support function c{F(x, t),^) of the arbitrary set F(x,t) (which is convex compact set from 1" for all t G [0,T] and for all x G M") is differentiate in

dc(F (x,t),i>)

the phase variable x and that the vector-function---is continuous in (x, ib.t).

dx

Then for all x,y G C"[0, T], and for all ^ G S,t G [0, T] the following relation holds:

c{F(x(t) + ay(t), t),i>) - c(F(x(t),t), ^) =

jdc(F(x(t),t),VQ dx

o(a,x(t),y(t),^,t)

a

Let v G Pn [0, T}. Fix the point z G Pri [0, T]. Put

^ 0, when a ^ 0.

t

za(t) = z(t) + av(t), y(t) = J v(t)dT. (14)

0

Using well-known properties of support functions [32] and equations (13), (14), let us calculate

£(^,za,t) = £(ip,z,t) + aH1(^,z,v,t) + o(a,^,t), (15)

o(a, ^, t)

^ 0, when a ^ 0,

where

t

Н\(ф, z, v, t) = (ф,v(t))-(J v(t )d

0

dc(F(x(t),t),t)\ dx /

Using relations (8), (9), find

h(za, t) = h(z, t) + aH(z, v,t) + o(a, t),

o(a, t) v ' ->• 0, when a 4. 0,

a

here

H(z,v,t) = m_ax Hi(ip, z,v,t), if тглх£(ф, z,t) > 0,

ipeR(z,t)

a

H(z,v,t) = 0, if max £(^,z,t) < 0,

^es

H(z,v,t) = max max{0, Hi(ip, z, v, t)}, if max. z,t) = 0,

ipeR(z,t)

with

R(z, t) = \ ip{z, t) € S | max{0, i(ip, z, t)} = max max{0, i(ip, z, t)} ^.

L ^(t)es J

Due to the structure of functional (8) it is easy to see, that in the case £(^, z,t) > 0 the maximum of the function max{0, £(^, z,t)} = £(^, z, t) is reached at the only element ^*(z,t) G S. Actually, in this case the point z does not belong to the set F. From the properties of the support function, it is known [33], that at every fixed t G [0, T] the value h(z(t),t) is the Euclidean distance from the point z(t) to the set F(x(t),t), that is

h(z(t),t)= \\z(t) - f (x(t),t)\\Rn ,

where f (x(t),t) is the projection of the point z(t) on the set F(x(t),t), which is unique as F(x(t),t) is the convex compact. Hence it is clear, that the maximum of the expression

h(z(t),t) = e(^,z(t),t) = (z(t),^) - c(F(x(t),t), ^) = \\z(t) - f (x(t),t)\\

is reached at the vector

r(z(t),t) = (z(t) - f (x(t),t))/\\z(t) - f (x(t),t)\\,

which is unique because of the uniqueness of the vector f(x(t),t). So in this case the set R(z,t) consists of the only element ^*(z,t).

Now it is not difficult to obtain the expansion

T

<p(za) = <(z) + a J h(z(t),t)H (z (t),v(t),t)+ o(a,t)dt. (16)

o

It is obvious, that if z G Q, then the functional < is Gateaux differentiable, and its Gateaux gradient vanishes.

Consider the case z G Q- Denote

. Фes

so in this case

(z,t) = ^max£(^,z,t)J > 0

T

<fi(z) = \ J w(z,t)dt,

0

as

w(z,t) = h2(z,t) = [maxmax{0,i(^,z,t)Y^ 2

o2(j,*t

(max £^,z,t)) = £2^*(z,t),z,t),

v ^gs у

because z ф i. As it has been already noted, in this case this maximum is reached at the only element ^*(z,t). Then the function ф* (z, t) is continuous in z at the fixed t [34],

w

therefore, in view of the continuity of the function ^c(F(x,t), ip) ^ conciUf[e) that the

dw(z t)

function —— is continuous in 2 at the fixed t. Further, by Lagrange's mean value dz

theorem there exists such a number 6 = 6(t) G [0,that (note that the function ^*(z,t) is in fact continuous in (z,t) [34], so all the integrals below exist)

/ dw(z + 6av,t)

wiz + av, t) — wiz, t) = a(---, v

\ dz

/ dw(z,t) \ / dw(z + 6av,t) dw(z,t) \ dz ' / \ dz dz

V Vt e [0,T].

Then we have

+ av) — <p(z) 1

dw(z(t),t) )dt dz

dw(z(t) + 6av(t),t) dw(z(t),t) , v(t),-------)at

dz

dz

<

1

4тах|Ж|| /

dw(z(t) + Oav(t), t) dw(z(t),t)

dz

dz

dt.

(17)

The multiplier may] l|v(t)|| ^s bounded, as v G P"[0,T]. The norm under the integral tends to zero, when a I 0 at the every fixed t G [0, T] due to the continuity of the function

dw(z,t) ,

dz

zt

dC(F(x,t), ^)

function ---'—!—. Then by Lebesgue's dominated convergence theorem expression

dx

(17) tends to zero, when a I 0, what proves Gateaux differentiability of the functional ip in the case z G Note that this reasoning also means that in formula (16) we in fact have

that — ->• 0, when a 4- 0, uniformly in t e [0, T],

a

Using the rule of the maximum function differentiation [34], calculate

dw(z,t) dl(ip*,z,t) / dc(F(x, t), ip*)\ ,

= 2madz =2 h(z,t)^----). (18)

Using classical variation, it is easy to calculate Gateaux gradient of the functional

i

V\(z) = x0 + j z(t)dt — xT.

(19)

function

d)

dc(F(x, t), ip) / N _

- holds true for support function (6), if Xj 0 Vi = l,n. By the

dx

assumption made this condition is violated only at the finite number of points tij. So we can divide the interval [0,T] in finite number of the intervals [0,t11],... and considered all the integrals separately on these intervals. Then on each of the corresponding open

2

a

dc(F (x,t),p)

intervals (0,tn),... the existence and the continuity of the vector-function---

holds true for support function (6) as on these intervals each coordinate Xj, i = 1, n, retains its sign. So for all such intervals we can precisely repeat the previous proof, made for the whole interval [0,T]. Hence incorporating these obvious minor changes in the previous reasoning (see formulas (15), (17) and (18), and making use of (19) we finally get the following theorem.

Theorem 1. Let x G X[0,T]. The functional I is Gateaux different-table at the point z and, its Gateaux gradient at the point z is given by the formula

t

T

+ xo + J z(t)dt - XT. (20)

o

dC(F(X,t), p*)

In (20) the expression--7—-- is given by the formula (7) (with ip* instead of tp in

dx

this formula).

The following theorem formulates the known minimum conditions of a differentiate functional.

Theorem 2. Let x* G X[0,T], For the point z* to minimize the functional I, it is necessary to have

n h( * * T{h(*(\ ,dc(F(x*(T),T),r(z*(r),T)) , ,

0„ = h{z ,t)ip [z ,t) - / h{z (t),t)-—-dr +

t

T

+ x0 + J z*(t)dt - xT, (21)

o

here 0n is the zero element of the space Pn[0,T]. It is apparent, that if I(z*) = 0, then the condition (21) is also sufficient.

4. The steepest descent method. Describe the well-known steepest descent

I

z1 G Pn[0,T]. Let the point zk G Pn[0,T] be already constructed. In order to apply

this method correctly, we still assume that xk G X [0,T]. If minimum condition (21) holds,

zk I

Otherwise let us put

zk+1 = zk - YkVI(zk),

where the vector-function xk (t) = x0 + / zk(r)dr, ^^d ^^e magnitude jk is the solution

o

of the following one-dimensional minimization problem:

minI(zk - yVI(zk)) = I(zk - YkVI(zk)). (22)

Due to (22) I(zk+1) < I(zk). If the sequence {zk} is finite, then its last point is the

I

Let the functional VI be uniformly continuous and bounded on the ball in Ln[0,T] with the center in the origin and the radius r' > sup HzH^o-y]. We suppose that the

Lebesgue set Z1 = {z G Pn[0,T] | I(z) < I(z1)} is bounded in the Ln[0,T] norm. If the

{zk}

|| VI (zk)\\L2 [0 T ] ^ 0, when k

5. Numerical examples. Let us consider some examples of the implementation of the proposed algorithm. In both examples the process was interrupted at the k*-th iteration, provided that \\VI(zk*T] ^ e* = 10~5. Example 1. Consider the interval system

x1 = [1, 3]x1 + [-4, -2]x2,

x2 = [0, 2]x1 + [-3, 0]x2, where t G [0,1]. There are the boundary conditions

x(0) = (1, 2)', x(1) = (-5.45, -0.75)'.

It is required to find a solution of this differential inclusion, which satisfies the given boundary conditions. In this example

c(F,P) = (2x1 - 3x2)p1 + (x1 - 1.5x2)p2 + (|x1| + |x2|)|p1| + (|x11 + 1.5lx2\)lp2\. Put z1 = (1,1)', then x1 = (1+1, 2 +1)'.

Four iterations were made according to the proposed algorithm. As a result the point x4(t) = 1 - 5.3358t - 1.4325t2 + 0.3183t3,

x4(t) = 2 - 1.8884t - 1.2197t2 + 0.4775t3 - 0.1194t4

was constructed (it is written out with precision up to the 4-th sign), see Fig. 1, a, b.

x*

Besides, the interval system is also satisfied, as (it can be easily checked) £(p, z*, t) < 0 Vt G

[0,1] (see (2) and (11)), see Fig. 1, c, d. Note, that due to the structure of the functional (21)

x*

are satisfied exactly. (This is because on some iteration! (including the last iteration) the solution of the problem (22) was obtained exactly.) So 11VI(z4)||L2 [0 T] = 0. Example 2. Consider the interval system

x1 = -x1 + [t, 2t]x2,

x2 = x1 - 2x2 + [-t2, 0], where t G [0,1]. There are the boundary conditions

x(0) = (2, 3)', x(1) = (1.5,1)'.

Fig. 1. Solution of example 1

It is required to find a solution of this differential inclusion, which satisfies the given boundary conditions.

In this example

c(F, ^) = (-x1 + 1.5tx2+ (x1 - 2x2 - 0.5t2+ 0.5tx1|^11 + 0.5t21

Put z1 = (1,1)', then x1 = (2 + t, 3 + t)f.

Seven iterations were made according to the proposed algorithm. As a result the point x7(t) = 2 - 1.0002t3 + 2A732t2 - 1.9937t - 0.0253t4 + 0.0425t5 + 0.0014t6 + 0.0022t7,

x2(t) = 3 - 1.1473t3 + 3.003t2 - 4.0125t + 0.2522t4 - 0.1034t5 + 0.0125t6 - 0.0044t7

was constructed (it is written out with precision up to the 4th sign), see Fig. 2, a, b. As it has been mentioned, \ (z7)\\L2 T] ^ 10_5.

One can make sure, that the boundary conditions for the trajectory x* are satisfied. Besides, the interval system is also satisfied, as (it can be easily checked) £(^,z*,t) < 0 Vt G [0,1] (see (2) and (11)) (with the set accuracy 10"5), see Fig. 2, c, d. Note the boundary conditions are satisfied exactly.

Remark 6. The choice of the accuracy £* is based on the balance between carrying out a not very big number of iterations and making compliance with sufficient for applications accuracy. The calculations were performed symbolically in the Maple 12.0 package. Note that various other known methods may be applied to minimize functional (12). The steepest descent method has been used for simplicity. So further studies are of interest for improving the computing efficiency of the realization of the concepts proposed in the article (this improvements can be carried out with the help of using discretization and also some other methods, such as conjugate gradients, to increase the effectiveness of the minimization of the given functional).

Fig. 2. Solution of example 2

Remark 7. As it has been noted in Introduction, the methods of estimating the attainability set may be effectively used while solving examples 1, 2. For the problems considered, the bounds provided by such methods are exact [1-3]. Note that solution of the considered problems exists if and only if the vector xT lies within these bounds at t = T. However, the proposed method of this paper gives an approach for solving more general optimization problem (see below Remark 8), while standard interval methods can't be applied to such problems straightforward.

Remark 8. Consider the following more general problem. It is required to minimize the functional

J(z) = J f (x,z,t)dt

subject to constraints (1), (3), (4), where we have x(t) = x0 + z(r)dr and the function

J 0

f (x, z, t) is continuous in (x, z, t) and continuously differentiate in x and in X. Note that under these assumptions the functional J(z) is Gateaux differentiate and its Gateaux

z

= dfi^llA + 1 df{x{r)z{T),r)

dz J dx

t

Using the results obtained in this article it is easy to give two rather straightforward approaches in attempt to solve this problem.

One approach for solving this problem is to construct the functional

Ф^ ) = J (z) + A^

0

T

J h2(z, t)dt.

In the case of the free right end (i. e., when constraint (4) is absent) the recent work [37] contains some natural sufficient conditions for penalty functional $(z) to be completely-exact. Speaking informally, it means that there exists the finite number A* such that the problem of minimizing functional J(z) under constraints (1), (3) is equivalent in some sense to the unconstrained minimization of the functional $(z) for all A > A*. We don't give the details here; see the paper [37] for rigorous formulations. This result is more of theoretical interest, as the functional $(z) is non-differentiable, besides its structure seems to be too complicated for applying methods of even nonsmooth optimization.

Another approach for solving this problem (including the case of the fixed right end) is considering the unconstrained minimization of the functional

^(z) = J (z) + nI (z)

with the increasing value of the parameter ^ to satisfy the constraints (1), (3), (4) with the desired accuracy (see formula (12)). Although the functional ^(z) is Gateaux differentiate, it is known that there are computational difficulties occurring while minimizing this functional caused by its "ravine" structure, when the value of ^ is rather big. We also don't consider this functional in detail here. See [38] for the known theorems on the connection

J(z)

the unconstrained minimization of the functional ^(z) with the increasing value of the parameter as well as some techniques aimed at overcoming the noted computational difficulties. We give the following simple example of solving a real physical problem with the help of this approach. As previously, we use the steepest descent method for finding the stationary points of the functional ^(z) in this example.

Example 3. Consider the interval system

x1 = x2,

x2 = [-1, -0.5]x1 + [-2,-1]x2 - 9.8, where t G [0,1]. There is the initial condition

x(0) = (1,1)'

and the right end is free. It is required to find both solutions of this differential inclusion, which satisfy the given initial condition and minimize as well as maximize the functional

i

jm =

л.

Such a statement of the problem admits the following physical interpretation. Consider a weight suspended on an ideal spring in a gravity field. In addition to the gravity force and the force from the the spring, the friction force ¡3v from the environment acts on the weight, which we assume to be proportional to the weight speed v, where 3 is the

coefficient of resistance of the environment. Write down the Newton's second law for the

y

my = -ky - 3v - mg,

here m is the weight mass; k is the spring stiffness; y is the coordinate measured from the state, in which the spring is not stretched in this oscillatory system. Now suppose that the

spring stiffness к as well as the environment resistance coefficient в are the functions of time, which can vary (but lie in the known intervals at every moment of time) depending on changes in the environment and the spring properties (for example while heating, etc.).

It is now interesting to find the minimal and the maximal possible under these conditions

2

^HV

values of the average kinetic energy —— of the system during the considered time interval.

If we now put m = 1, g = 9.8, k(t) e [0.5,e(t) e [1, 2], t e [0,1] and denote x1 = y, x2 = v = y, we will get the differential inclusion and the cost functional above. In this example we have

c(F,ф) = Х2Ф1 + (-0.75x1 - 1.5x2 - 9.8)ф2 + (0.25\x1\ + 0.5\х2\)\ф2\.

Put ц = 85, z1 = (1,0)', then x1 = (1 +1,1)'. The results in this example are written out with precision up to the 6th sign.

Example 3.1. To minimize the functional ^(z), 15 iterations were made according to the proposed algorithm. As a result the point

x15(t) = 1 - 0.001341t5 - 0.219599t4 + 1.801792t3 - 5.645627t2+

+ 0.006182t6 + t, 0 < t < 0.092, 0.998539 + 0.531280t5 - 1.741569t4 + 4.040669t3 - 6.193228t2-- 0.049752t6 + 1.048461t - 0.065932t7 + 0.031599t8, 0.092 < t < 0.624, 0.061205t6 + 3.775047t3 - 6.51987t2 - 0.975091t4-

- 0.004879t8 + 0.931305 + 1.344730t, 0.624 <t < 1, x25 (t) = 1 + 0.037092t5 - 0.0066704t4 - 0.878399t3+

+ 5.405375t2 - 11.291254t, 0 < t < 0.092, 0.252796tr - 0.461527t6 - 0.298511t5 + 2.656402t4 - 6.966278t3+ + 12.122007t2 - 12.386457t + 1.048461, 0.092 < t < 0.624,

- 0.039036tr - 13.03974t + 0.367231t5 - 3.900365t3+

+ 11.325142t2 + 1.344730, 0.624 <t < 1,

was constructed (see Fig. 3, I, a). The process was terminated at this iteration as 11 V^(z15) 11[0 T] < 2.5 • 10-2. So the obtained point z* = z15 is at least the stationary (with the set precision) point of the functional ^(z).

One can make sure that the initial condition for the trajectory x* is satisfied and we have I(z*) = p(z*) = 1.1 • 10-4 (so the interval system is satisfied (according to criterion (11)) with the set accuracy 5 • 10-4 and, as can be easily verified, the point-wise violation of the interval system does not exceed the value 1.5 • 10-2) and J(z*) = 7.711839. Fig. 3, I, b shows that the trajectory x*(t) lies practically on the boundary of the permissible region.

Example 3.2. To maximize the functional ^(z), 32 iterations were made according to the proposed algorithm. As a result the point

x32(t) = 1 - 6.40899t2 + 4.100813t3 - 1.555853t4+

+ 0.618022t5 +1, 0 < t < 0.0843,

Fig. 3. Solution of examples 3.1 (J) and 3.2 (II)

1.001158 - 0.099432t4 - 5.893454t2 + 1.830437t3 + 0.295029t5+

+ 0.95809t - 0.203586t8 + 0.554507t7 - 0.618585t6, 0.0843 < t < 0.562,

1.050245 - 0.25779t4 - 5.571737t2 - 0.009198t7 + 1.847321t3+

+ 0.020355t6 + 0.001292t8 + 0.720453t, 0.562 <t < 1,

x22(t) = 1 + 3.090111t4 - 6.223412t3 + 12.30244t2-

- 12.81798t, 0 < t < 0.0843,

0.95809 + 3.881548t6 - 1.628687t7 - 3.711512t5 + 1.475147t4-

- 0.397728t3 + 5.491312t2 - 11.786907t, 0.0843 < t < 0.562,

0.720453 - 0.064387t6 + 0.010339t7 + 0.122129t5 - 11.143474t -

- 1.03116t3 + 5.541962t2, 0.562 <t < 1,

was constructed (see Fig. 3, II, a). The process was terminated at this iteration as ||V^(z32)||l2 j0Tj < 2.5 • 10"2. So the obtained point z* = z32 is at least the stationary (with the set precision) point of the functional ^(z).

One can make sure that the initial condition for the trajectory x* is satisfied and we have I(z*) = ^(z*) = 3.65 • 10"4 (so the interval system is satisfied (according to criterion (11)) with the set accuracy 5 • 10"4 and, as can be easily verified, the point-wise violation of the interval system does not exceed the value 2.5 • 10-2) and J(z*) = 14.085116. Fig. 3, II, b shows that the trajectory X* (t) lies practically on the boundary of the permissible region.

So we obtain the lower and the upper bounds of the system average kinetic energy change with a possible variation of the parameters k and p in the specified intervals: we have 7.711839 < J(z) < 14.085116.

6. Conclusion. Thus, in this paper the problem of finding the solution of a linear nonstationary interval system of differential equations with the given boundary conditions

has been considered. With the help of the support functions, the original problem is reduced to minimizing some functional. For it, Gateaux gradient is found, necessary and sufficient conditions for a minimum are obtained. On the basis of these conditions, a numerical method (the method of the steepest descent) for solving the initial problem is constructed. The numerical examples illustrate implementation of this algorithm. The possibility of applying the developed in the article method to finding the optimal (in the sense of the integral functional) solution of the considered interval system is noted, and this approach is demonstrated on a simple physical example.

References

1. Durieu C., Walter E., Polyak B. Multi-input multi-output ellipsoidal state bounding. Journal of Optimization Theory and Applications, 2001, vol. Ill, no. 2, pp. 273-303.

2. Harrison G. W. Dynamic models with uncertain parameters. Proceedings of the First International Conference on Mathematical Modeling. Ed. by X. J. R. Avula. Rolla, University of Missouri Publ., 1977, vol. 1, pp. 295-304.

3. Maksarov D. G., Norton J. P. Computationally efficient algorithms for state estimation with ellipsoidal approximation. International Journal of Adaptive Control and Signal Processing, 2002, vol. 16, no. 6, pp. 411-434.

4. Jaulin L. Nonlinear bounded-error state estimation of continuous-time systems. Automatica, 2002, vol. 38, no. 6, pp. 1079-1082.

5. Scott J. K., Barton P. I. Bounds on the reachable sets of nonlinear control systems. Automatica, 2013, vol. 49, pp. 93-100.

6. Singer A. B., Barton P. I. Bounding the solutions of parameter dependent nonlinear ordinary differential equations. SIAM J. Sci. Comput., 2006, vol. 27, no. 6, pp. 2167-2182.

7. Kieffer M., Walter E., Simeonov I. Guaranteed nonlinear parameter estimation for continuous-time dynamical models. Proceedings of l/tth IFAC Symposium on System Identification, 2006, pp. 843-848.

8. Dubois D., Prade H. Towards fuzzy differential calculus. Pt 3. Differentiation. Fuzzy Sets and Systems, 1982, no. 8, pp. 225-233.

9. Chang S. L., Zadeh L. A. On fuzzy mapping and control. IEEE Trans. Systems Man Cybernet, 1972, vol. 2, pp. 30-34.

10. Goetschel R., Voxman W. Elementary fuzzy calculus. Fuzzy Sets and Systems, 1987, no. 24, pp. 31-43.

11. Kandel A., Byatt W. J. Fuzzy differential equations. Proceedings of International Conference on Cybernetics and Society, 1978, pp. 1213-1216.

12. Kaleva O. The Cauchy problem for fuzzy differential equations. Fuzzy Sets and Systems, 1990, no. 35, pp. 389-396.

13. Stefanini L., Bede B. Generalized Hukuhara differentiability of interval-valued functions and interval differential equations. Nonlinear Analysis, 2009, vol. 71, pp. 1311-1328.

14. Ma M., Friedman M., Kandel A. Numerical solutions of fuzzy differential equations. Fuzzy Sets and System-s, 1990, no. 105, pp. 133-138.

15. Pederson S., Sambandham M. Numerical solution to hybrid fuzzy systems. Mathematical and Computer Modelling, 2007, vol. 45, pp. 113-144.

16. Jayakumar T., Maheskumar D., Kanagarajan K. Numerical solution of fuzzy differential equations by Runge^Kutta method of order five. Applied Mathematical Sciences, 2012, vol. 6, pp. 2989-3002.

17. Bandyopadhyay A., Kar K. Fuzzy continuous dynamical system: a multivariate optimization technique. Annual meeting of the North American Fuzzy Information Processing Society. Berkeley, CA, USA, 2007, pp. 1-7.

18. Abbasbandy S., Viranloo T. A. Numerical solutions of fuzzy differential equations by Taylor method. Computational Methods in Applied Mathematics, 2002, vol. 2 (2), pp. 113-124.

19. Morozov A. Yu., Reviznikov D. L. Modifikacija metodov reshenija zadachi Koshi dlja sistem obyknovennyh differencial'nyh uravnenij s interval'nymi parametrami [Modification of methods for solving the Cauchy problem for systems of ordinary differential equations with interval parameters]. Proceedings of MAI Institute, 2016, vol. 89, pp. 1-17. (In Russian)

20. Hashemi M. S., Malekinagad J., Marasi H. R. Series solution of the system of fuzzy differential equations. Advances in Fuzzy Systems, 2012, pp. 1-16.

21. Liao S. Homotopy analysis method: a new analytical technique for nonlinear problems. Communications in Nonlinear Science and Numerical Simulation, 1997, no. 2, pp. 95-100.

22. Mosleh M., Otadi M. Minimal solution of fuzzy linear system of differential equations. Neural Computing & Applications, 2012, vol. 21, pp. 329-336.

23. Wu J., Zhang Yu., Chen L., Luo Zh. A Chebyshev interval method for nonlinear dynamic systems under uncertainty. Applied Mathematical Modelling, 2013, vol. 37, pp. 4578-4591.

24. Chen M., Wu C., Xue X., Liu G. On fuzzy boundary value problems under uncertainty. Information Sciences, 2008, no. 178, pp. 1877-1892.

25. Sokolova S. P., Kuzmina E. A. Dinamicheskie svojstva interval'nyh sistem [Dynamic properties of interval systems]. SPIIRAS Proceedings, 2008, vol. 7, pp. 215-221. (In Russian)

26. Rogoza A. A. Ob odnom podhode к postroeniju dvustoronnih ocenok mnozhestv reshenij nelinejnyh difFerencial'nyh uravnenij s interval'nymi parametrami na osnove proekcionnyh metodov [On an approach to constructing two-sided estimates of solution sets of non-linear differential equations with interval parameters on the basis of projective methods]. Differential Equations and Control Processes, 2015, vol. 2, pp. 19-58. (In Russian)

27. Dobronec B. S., Strelnikova E. A. Reshenija kraevyh zadach dlja ODU s interval'nymi neo-predelennostjami [Solutions of boundary value problems for ordinary differential equations with interval uncertainties]. Computing Technologies, 2005, vol. 10, pp. 16-21. (In Russian)

28. Gasilov N. A., Fatullayev A. G., Amrahovb S. E. Solution method for a non-homogeneous fuzzy linear system of differential equations. Applied Soft Computing, 2018, vol. 70, pp. 225-237.

29. Hullermeier E. An approach to modelling and simulation of uncertain dynamical systems. international Journal of Uncertainty, Fuzziness and Knowledge-Based Systems, 1997, vol. 5, no. 2, pp. 117-137.

30. Hanss M. The transformation method for the simulation and analysis of systems with uncertain parameters. Fuzzy Sets and Systems, 2002, no. 130, pp. 277-289.

31. Fominyh A. V. A numerical method for finding the optimal solution of a differential inclusion. Vestnik St. Petersburg University. Mathematics, 2018, no. 51, pp. 397-406.

32. Blagodatskih V. I., Filippov A. F. Differencial'nye vkljuchenija i optimal'noe upravlenie [Differential inclusions and optimal control]. Proceedings of the Steklov Institute of Mathematics, 1985, vol. 169, pp. 194-252. (In Russian)

33. Blagodatskih V. I. Vvedenie v optimal'noe upravlenie [Introduction to optimal control]. Moscow, Vysshaya shkola Publ., 2001, 239 p. (In Russian)

34. Bonnans J. F., Shapiro A. Perturbation analysis of optimization problems. New York, Springer Science+Business Media Press, 2000, 601 p.

[]

Publ., 1977, 752 p. (In Russian)

36. Penot J. P. On the convergence of descent algorithms. Computational Optimization and Applications, 2002, vol. 23, no. 3, pp. 279-284.

37. Dolgopolik M. V., Fominyh A. V. Exact penalty functions for optimal control problems. I: Main theorem and free-endpoint problems. Optimal Control Applications & Methods, 2019, vol. 40, no. 6, pp. 1018-1044.

[]

(In Russian)

Received: November 29, 2020. Accepted: April 05, 2021.

Author's information:

Alexander V. Fominyh — PhD in Physics and Mathematics, Associate Professor; alexfomster@mail.ru

Метод нахождения решения линейной нестационарной интервальной

*

А. В. Фоминых

Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9

* Основные результаты статьи (разделы 3 и 5) были получены в ИПМаш РАН при финансовой поддержке Российского научного фонда (проект № 20-71-10032).

Институт проблем машиноведения Российской академии наук, Российская Федерация, 199178, Санкт-Петербург, Большой проспект В. О., 61

Для цитирования: Fominyh А. V. Method for finding a solution to a linear nonstationary interval ODE system // Вестник Санкт-Петербургского университета. Прикладная математика. Информатика. Процессы управления. 2021. Т. 17. Вып. 2. С. 148-165. https://doi.org/10.21638/11701/spbul0.2021.205

Исследуется линейная нестационарная интервальная система обыкновенных дифференциальных уравнений, в которой элементами матрицы ее правой части являются интервалы с известными нижними и верхними границами. Система определена на известном конечном интервале времени. Требуется построить траекторию, которая приведет эту систему из заданного начального положения в заданное конечное состояние. Исходная задача — нахождение решения дифференциального включения специального вида с закрепленным правым концом. С помощью опорных функций такая задача сводится к минимизации функционала в пространстве кусочно-непрерывных функций. При естественном дополнительном предположении данный функционал дифференцируем по Гато. Для функционала найден градиент Гато, получены необходимые и достаточные условия минимума. На основе этих условий к исходной задаче применяется метод наискорейшего спуска. Некоторые численные примеры иллюстрируют реализацию построенного алгоритма.

Ключевые слова: линейная нестационарная интервальная система обыкновенных дифференциальных уравнений, дифференциальное включение, опорная функция, метод наискорейшего спуска.

Контактная информация:

Фоминых Александр Владимирович — канд. физ.-мат. наук, доц.; alexfomster@mail.ru