Вычисление показателя Ляпунова обобщенных линейных систем с показательным распределением элементов переходной матрицы Текст научной статьи по специальности «Математика»

ВЫЧИСЛЕНИЕ ПОКАЗАТЕЛЯ ЛЯПУНОВА ОБОБЩЕННЫХ ЛИНЕЙНЫХ СИСТЕМ С ПОКАЗАТЕЛЬНЫМ РАСПРЕДЕЛЕНИЕМ ЭЛЕМЕНТОВ ПЕРЕХОДНОЙ МАТРИЦЫ*

Н. К. Кривулин

С.-Петербургский государственный университет,

канд. физ.-мат. наук, доцент, [email protected]

1. Введение. При исследовании различных систем в экономике, технике и других областях находят применение обобщенные линейные динамические модели вида

2 (к) = А(к)2 (к — 1),

где А(к) — случайная матрица, 2 (к) — вектор состояний, а умножение понимается в смысле идемпотентной алгебры [1-3] с операциями вычисления максимума и сложения.

Во многих случаях представляет интерес определение средней скорости роста вектора состояний (показателя Ляпунова) системы

А = Иш \\г(к)||1/к, к——

где символ нормы и степень понимаются в смысле идемпотентной алгебры.

Однако вычисление значения предела часто оказывается довольно трудоемкой задачей даже для весьма простых моделей систем невысокого порядка. Известные результаты в этой области включают решения, полученные в работах [4-8] для систем второго порядка. В работе [4] изучена система с матрицей, элементы которой независимы и имеют экспоненциальное распределение с единичным средним. Решение сводится к построению последовательности одномерных распределений для разности компонент вектора системы с последующим решением интегрального уравнения для нахождения предельной плотности, которая используется для получения результата в виде А = 407/228.

Подход, предложенный в работе [5], опирается на построение последовательности двумерных распределений вероятностей вектора состояний и ее анализ при помощи преобразования Лапласа соответствующих функций распределения. Для системы с матрицей, у которой диагональные элементы распределены экспоненциально с параметром м, а недиагональные элементы — экспоненциально с параметром V, получен результат (см. также [6]) А = Р(м, v)/Q(м, V), где

Р (м, V) = 160м10 + 1776м9V + 8220^2 + 21378^3 + 35595^4 + 41566^5 +

+ 35595^6 + 21378^7 + 8220^8 + 1776^9 + 160v10,

Q(м, V) = 16^(м + v)(8м8 + 80м7 V + 321м^2 + 690м^3 + 880м*^4+

+ 690^3 V5 + 321^6 + 80^7 + 8v8).

Кроме того, в [5] показано, что для системы второго порядка с симметричной матрицей, у которой диагональные элементы имеют экспоненциальное распределение с параметром м, а недиагональные элементы равны нулю, выполняется А = 5/(4м). Этот

* Работа выполнена при финансовой поддержке РФФИ (грант №09-01-00808).

© Н. К. Кривулин, 2009

результат был затем обобщен в работе [7] на случай, когда диагональные элементы имеют распределения с различными параметрами м и т:

_ м4 + /л3т + м2г2 + /лт3 + т4 /Хт(/Х + т)(/.42 + Г2)

В работе [8] рассматривается система с матрицей второго порядка, элементы которой имеют экспоненциальные распределения с произвольными параметрами. Задача нахождения средней скорости роста сводится к ряду алгебраических вычислений, включая решение алгебраической системы уравнений и вычисление значения некоторого линейного функционала от полученного решения.

Целью настоящей работы является определение средней скорости роста вектора состояний для ряда систем второго порядка с матрицами, некоторые элементы которых имеют экспоненциальные распределения и независимы, а другие элементы являются нулевыми константами. Используется подход, развитый в [7], который опирается на построение и анализ некоторой последовательности одномерных функций распределения. Средняя скорость роста затем находится как математическое ожидание случайной величины, которая определяется предельным распределением этой последовательности.

2. Линейная стохастическая динамическая система. Рассмотрим динамическую систему, эволюция которой при всех к = 1, 2, . . . описывается векторным уравнением в полукольце с операциями вычисления максимума и сложения

2 (к) = А(к)2 (к — 1),

где

а«=(а *)■ 2«=(Уч)’ 2<»>=(°).

Векторное уравнение можно представить в виде системы скалярных уравнений

ж(к) = акх(к — 1) ® вку(к — 1),

У(к) = 7кх(к — 1) ® 5ку(к — 1)

или с использованием знаков обычных арифметических операций в виде

ж(к) = шах(х(к — 1) + ак, у(к — 1) + вк), у(к) = шах(х(к — 1) + 7к, у(к — 1) + 4).

Пусть каждая последовательность {ак}, {вк}, |7к} и } состоит из независимых случайных величин, которые имеют общее экспоненциальное распределение. Случайные величины ак, вг, 7т и 5п являются независимыми при любых к,/,ш, п.

Средняя (асимптотическая) скорость роста вектора состояний, которую часто называют показателем Ляпунова системы [5], определяется в идемпотентной алгебре как предел

А = Иш ||г(к)||1/к к——

или, с использованием обычных обозначений, как предел

А= Ит —та х(х(к),у(к)). к—к

Применяя эргодическую теорему из [9], нетрудно показать, что когда элементы матрицы А(к) неотрицательны и имеют конечные средние, этот предел с вероятностью 1 существует и является конечным числом. При этом также выполняется равенство

Ит —Ета х(х(к),у(к)) = X.

к——ж к

3. Вычисление средней скорости роста вектора состояний. Для решения задачи вычисления А применим замену переменных, предложенную в [4],

X (к) = ж(к) — ж(к — 1), У (к) = у (к) — ж(к),

после которой уравнения системы приобретают вид

X (к) = шах(ак,вк + У (к — 1)),

У (к) = шах(7к, 4 + У (к — 1)) — шах(ак, вк + У (к — 1)).

Обозначим функции распределения вероятностей величин ак, вк, 7к и ^ через

^а(£) = шах(0,1 — е-м4), ^3 (4) = шах(0,1 — е-^),

(4) = шах(0,1 — е-ст4), ^5 (4) = шах(0,1 — е-т4),

а их функции плотности через /а, /, /7 и /5.

Введем функции распределения

Фк(4) = Р{Х(к) < 4}, Фк(*) = Р{У(к) < 4}.

Функцию Фк можно представить в виде

р ОО

Фк(4) = -^а(£) / Фк-1(4 — и)/в (и)^и.

0

Предположим, что последовательность функций Фк при к ^ то сходится к некоторой функции распределения Ф. Тогда, применяя теорему Лебега об ограниченной сходимости, приходим к выводу, что последовательность Фк сходится к функции распределения некоторой случайной величины X:

Ф(4) = ^а(*п Ф(4 — и)/в (м)^и. (1)

0

Из сходимости распределений следует, что последовательность ЕХ(к) сходится к ЕХ.

Наконец, учитывая, что ж(к) = X(к) + • • • + X(1), а также принимая во внимание неравенство Е шах(0, У (к)) < Е(7к + #к) < то, получим

1 1 1 „к

А= Ит — Е тах(ж(/г), г/(/г)) = Ит — Е тах(0, У (к)) + Ит — ЕХ (г) = ЕХ.

к к к к к к

г=1

4. Система с матрицей с нулями на диагонали. Рассмотрим динамическую систему с матрицей

0 вк ^

А(к) = ' •» 0

Скалярные уравнения системы в обычных обозначениях принимают вид

х(к) = шах(х(к — 1), у(к — 1) + вй), у(к) = шах(х(к — 1) + 7Й, у(к — 1)).

После перехода к новым переменным X (к) и У (к) получим уравнения

X (к) = шах(0, У (к — 1) + вй),

У (к) = шах(7й, У (к — 1)) — шах(0, У (к — 1) + вй).

Функция Ф&(£) = Р{Х(к) < имеет вид

рЖ

Фй (4)=/ 1 (4 — и)/а(и)йи.

ио

Применяя формулу полной вероятности, функцию Ф^ = Р{У(к) < 4} запишем так:

рЖ р Ж

Фй(4) = / / Р{У(к) <4|вй = «,7й = «}/в(и)/7Н^и^.

7о 7о

Рассмотрим условную вероятность

Р{У (к) < 4|вй = и, 7й = V} = Р{шах(«, У (к — 1)) — шах(0, У (к — 1) + и) < 4}.

После преобразования выражения под знаком вероятности имеем

_1 (4), если V < 4, и < —4,

1, если V < 4, и > —4,

0, если V > 4, и < —4,

1 — Фй_1^ — и — 4), если V > 4, и > —4.

Подстановка в уравнение для Ф^ и замена переменных приводит к рекуррентному уравнению

т . . | (1 — аь_ 1)е^4, если 4 < 0,

Фй(4) = < й / _

|^1 — ай _ 1е , если 4 > 0,

где числа определяются выражением

/*Ж /* ОО

а^ = ^а / / Фй^ — и)е__ ст-иЙи^у.

оо

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

Из двух последних равенств имеем

1 2^ + а

а-к = ~пак-1 +

2 2(^ + а)'

Найденному рекуррентному соотношению соответствует некоторое отображение, заданное на множестве неотрицательных чисел, которое является сжатием. Следовательно, последовательность а^ сходится при к ^ то к неподвижной точке отображения а.

Переходя к пределу, получаем уравнение для определения а, решение которого дает

2^ + а

а = —---------

3(г/ + а)

Последовательность функций Ф& тогда тоже сходится к некоторой функции, которая имеет вид

{ если г < О,

ф(*) = \ 3(,у+<т9) , ’ _ , “

I 1 — ч?~г~те > если ь > о.

^ 3(^+ст) ’ ^

Ясно, что Ф является функцией распределения некоторой случайной величины. Найдем функцию распределения Ф исходя из (1). Предположим сначала, что а = V. После вычисления интеграла имеем

фм = 1 _

6^ — а) 3(^ — а2)

Окончательно получаем

4^2 + 7vа + 4а2

/*Ж

Л = 4ЙФ(4) =

о

6vа(v + а)

Нетрудно проверить, что этот результат остается верным и в случае, когда а = V.

5. Система с матрицей с нулевой строкой. Предположим, что матрица системы имеет вид

Л(к)=( а в*

Рассмотрим скалярные уравнения

х(к) = шах(х(к — 1) + а, у(к — 1) + вй), у(к) = шах(х(к — 1), у(к — 1)),

которые после замены переменных приобретают вид

X (к) = шах(а, У (к — 1) + вй),

У (к) = шах(0, У (к — 1)) — шах(ай, У (к — 1) + вй).

Учитывая, что У (к) < 0 при всех к, имеем У (к) = —X (к). Следовательно, анализ системы уравнений можно свести к анализу уравнения

X(к) = шах(ай, вй — X(к — 1)).

Функция распределения Ф^ (4) = P{X(к) < 4} имеет вид

Фй(4) = (1 — е м*)(1 — ай _ 1е иЬ)

при условии, что

Ж

ай = V / Фй(и)е _^“Йи.

о

Из двух последних соотношений имеем

2(М + 2^^" М + V'

Переходя к пределу при к то, получим уравнение для а, решением которого

является

2 /х(/4 + 2г/) а (/X + г/)(3/х + 4г/)'

Последовательность функций Фд сходится к функции

*(*) = С1 - -—) (і - )

которая, очевидно, является функцией распределения некоторой случайной величины. Вычисляя математическое ожидание, получаем

_ 2/х4 + 7ц3V + 10/42г/2 + П/хг/3 + 4г/

/хг/(/х + и)2 (3/х + 4г/) '

6. Система с матрицей с нулевым элементом на диагонали. Рассмотрим динамическую систему с матрицей

а№>=( а »

Запишем динамическое уравнение в виде системы скалярных уравнений

х(к) = шах(х(к - 1) + а, у(к - 1) + ДО, у(к) = шах(х(к - 1) + 7Й, у(к - 1)).

После замены переменных получим уравнения

X (к) = шах(ай, У (к - 1) + ДО,

У (к) = шах(7й, У (к - 1)) - шах(ай, У (к - 1) + Дк).

По формуле полной вероятности имеем

/•^О /*Ж /*Ж

Фй(і)=/ / / Р{У(к) <і|«к = и,вк = «,7к = т}/а(и)/а^)/7(т)йи^т.

^0 70 70

Найдем условную вероятность

Р{У(к) < і |ай = и, вк = V, 7к = ад} =

О, если и < т - і, V < -і,

1 - Фд _і(-V + т - і), если и < т - і, V > -і,

Фк _і(и + і), если и > т - і, V < -і,

_ 1, если и > т - і, V > -і.

С учетом обозначений

р СО рО

а& = ма / / Фк (и + v)e_(м+ст)u_мv Л^,

и0 -/0

л ОО /* оо

6к = / / Фд(-и + v)e_vu_стvЛ^,

^0 ^0

л О /* о

Ск = / / Фд(-и + v)e_v“_(м+ст)vЙи^у

00

/0 ./0

функция Фд принимает вид

Ф (Л - Iай-1в^4 + (1 - 6й_1)е^ - (ай-1 - сй-1)е(^+^)4, если 4 < 0, к |^1 — (6к-1 — ек-1)е-ст4, если 4 > 0.

Имеем рекуррентные соотношения

_ иа I, , иа , 17

ай — — 7-1--------, г, \bk-l + т---;---гт--, ~ ,сц +

(м + а)(м + 2а) (м + а)(м + 2а) м + а’

1 ^(м + 2^ + 2а) 2^ + а

-а-к-1 — + 777-;----77-----. Л чсй-1 +

(м + ^)(м + 2^)(^ + а) 2 ^ 2(^ + а)(м + 2^^х 2(^ + а)’

V2 а а(/х + 2г/ + 2а)

й (/х + г/)(/х + 2г/)(/х + V + а) й 1 2(/л + 2а) (/х + г/ + а) й 1 +

2^а а(м + 2^ + а)

+ 7—, „ ч „ чСй-1 +

(м + 2^)(м + 2а) 2(м + а)(м + V + а)

Введем обозначение Vк = (ак, 6к,ск)т. Определим матрицу

/ о

С =

_______(Л£_______ __(Л£__

(^+<т)(м+2(Т) (М+<7 )(/^+2<7 )

__________1/^(7__________ _ _1 ^/(М+2^'+2(7 )

(^+і/)(//+2і/)(і/+(т) 2 2(і/+(7 )(іі-\-2и)

. _____________іу2а______________________<7(іл-\-2і/-\-2<7 ) ________2 у а______

\ (ц-\-іу)(ц-\-2іу)(ц-\-іу-\-<7) 2{ц-\-2<т){ц-\-і/-\-<т) {ц-\-2і/){ц-\-2<т)

и вектор

/

И+а

4

у-\-<т

2 (>+<т) <т(іи+2і'+<т)

\ 2(/і+<т)(м+і'+<т)

Теперь рекуррентные соотношения можно представить в виде векторного уравнения

V к = к_і + 4, (2)

которое задает некоторое отображение трехмерного пространства в себя.

Пусть указанное отображение является сжимающим. Тогда последовательность векторов Vк сходится к неподвижной точке отображения V = (а, 6, с)т. Учитывая взаимно однозначное соответствие между последовательностями функций Фк и векторов V к, заключаем, что последовательность Фк также сходится к некоторой функции Ф. Неподвижная точка отображения определяется уравнением

(I - = 4. (3)

Проверка сходимости последовательности Vк и решение уравнения (3) в общем случае требуют слишком громоздких вычислений. Рассмотрим частные случаи, в которых решение задачи может быть получено в достаточно компактной форме.

6.1. Случай а = V = м- При условии а = V = м

О -1/6 1/6 \ / 1/2

С = ( 1/12 -1/2 5/12 I , 4 = I 3/4

1/18 -5/18 2/9 у \ 1/3

Вычисление евклидовой нормы матрицы дает ||С||2 = 397/648 < 1, откуда следует, что последовательность векторов V к сходится к неподвижной точке V. После решения уравнение (3) имеем координаты V в виде

а - 123 6- 165 с- 69

0 ~ 278’ ““ 278’ С~ 278'

Последовательность функций Фк тогда сходится к функции

{ 118 , 27 2 ,

---ем---------е м , если і < О,

139 48 139 - ’

1------е если і > О,

139 ’ ’

которая является функцией распределения некоторой случайной величины.

Найдем функцию распределения

89 48

ф(г) = (1 - 11-----------------

139 139

Вычисление математического ожидания дает результат

439

Л —

278м

6.2. Случай а = м, V = м- Матрица С и вектор д принимают вид

С=

______^ /г_____ 1 ^(3^+2^) | /г+2^

(ц+и)2(ц+2и) 2 2(/х+1')(м+21') 1 4 2(/х+1')

_________^___________ 3(л+2и 2у

\ (/л+1^)(/л+21х)(2/л+1х) 6(2/х+1') 3(ц+2и) ) \ 2(2^1+!/)

Нетрудно проверить, что усу2 < 1. Решение уравнения (3) имеет вид (^ + v)(17м3 + 50^,^ + 44^2 + 12v3)

6 =

36м4 + 147м^ + 215м2v2 + 130мv3 + 28v4 ’

12 + 67/л31/ + 134/х2г/2 + 95/хг/3 + 22г/4

36/х4 + 147/х3г/ + 215/х2г/2 + 130/хг/3 + 28И ’

(ц + г/) (Юг/3 + ЗІ/хг/2 + 22/і2і/ + б/х3)

36/х4 + 147/х3г/ + 215/х2г/2 + 130/хг/3 + 28И '

с

Теперь можно записать предельную функцию Ф, а затем найти функцию

Ф(*) = (1 - е-*“) I 1----—(с - Ъ)е-^~

\ ^ - V

V

2

— — —

/х + г/ г/(2/х + г/) ,

а, + —------------------о — ------------—-----------------с е

ч2 (^ + v)(yU, + 2v) 2(^ — V) (^ + 2v)(м — V)

Наконец, вычисляя математическое ожидание, получаем

Л 48/л5 + 238/х4г/ + 495/л3г^2 + 581/х2г/3 + 326/хг/4 + 68г/5 ““ 2/хг/(36 /х4 + 147/х3г/ + 215/х2г/2 + 130 /хг/3 + 28И)

6.3. Случай а = V. Запишем матрицу и вектор:

С:

/ П__________^^___________ \ / V

' (ц+и)(ц+2и) (ц+и)(ц+2и) ' '

2 1 |Ц+41/

2(^+1')(^+2!') 2 4(^+2^) I ^ Я.

_______________ и()1+Аи) 2 и2 '

\ (м+^)(м+2^)2 2(ц+2и)2 (>+2и)2

\ 2(м+^)(м+2^)

Заметим, что матрица О — вырожденная, причем для всех к выполняется

_ М + ^ 2г/ г/(/х - г/)

к 2(/х + 2г/) к /х + 2г/ к 2(/х + г/)(/х + 2г/)

Введем вектор Vк = (ак,Ьк)Т. Подставляя выражение (4) для сд в уравнение (2), приходим к уравнению V к = Н V к_ 1 + г, где

____ци __________//2^__ \ / 1/(^2 + 7^1/+8^2) \

2(/х+2|/) (м+^)(м+2^)2 | г = I 2(ц+и)([1+2и)2 \

ц,3 + 6ц,2 1у+5ц1у2 — Ау3 ц(ц+Зу) I ’ I 3,ц2 + 14|Ш/+20^2 I '

8(ц+1;)(ц+21;)2 2([1+2и)2 / \ 4(ц+2и)2 /

Это уравнение задает отображение двумерного пространства в себя. Как нетрудно проверить, уНу2 < 1. Учитывая, что отображение является сжатием, последовательность v'k сходится к неподвижной точке V' = (а, 6)т, которая определяется уравнением (I — Н)v' = г. Решение уравнения имеет вид

^ _ г/(6/х4 + 57/х3г/ + 173/х2г/2 + 192/лг^3 + 64г/4)

12/х5 + 97/х4г/ + 286/х3г/2 + 397/х2г/3 + 256/хИ + 64г/5 ’

_ 6/х6 + 55/хбг/ + 197/х4г/2 + 367/х3г/3 + 391/х2г/4 + 240г/5/х + 64г/6 “ (/х + г/) (12/х5 + 97/х4г/ + 286/х3г/2 + 397/х2г/3 + 256/хИ + 64г/5) '

Переходя к пределу в (4), найдем предел последовательности ск в виде

v(3/u5 + 27^ + 99^2 + 183^3 + 176^4 + 64v5)

с = ---------------------------------------------------------------------

(/X + V)(12/х5 + 97/х4г/ + 286/х3г/2 + 397/х2г/3 + 256/хИ + 64г/5)'

Запишем функцию распределения:

ф(() =(1 - е"‘> (' - Й - (У+Х+^)а+^ - с)<г“() •

45

1

Вычисляя математическое ожидание, получаем Л = Р(м, v)/Q(м, V), где

Р(м, V) = 15/ + 152м7V + 624м^2 + 1382м^3 + 1838м^4 + 1592м^5 + + 973м^6 + 384мv7 + 64v8, ф(м, V) = м^(м + v)2(12м5 + 97м^ + 286м^2 + 397м^3 + 256м^4 + 64v5).

7. Система с матрицей с нулевым элементом ниже диагонали. Рассмотрим динамическую систему с матрицей

Ясно, что для решения задачи вычисления Л можно применить подход, использованный выше для анализа системы с матрицей с нулевым элементом на диагонали. Учитывая, что при этом общий порядок решения сохраняется без существенных изменений, описание решения здесь опускается.

7.1. Случай т = v = м. При таком условии выполняется

515 Л —------.

352м

7.2. Случай v = м- Имеет место результат Л = P(м,т)/Q(m,t), где

P (м, т) = 288м8 + 1048м7т + 1936м6т2 + 2688м5 т3 + 3012м4т 4+

+ 2226м3т5 + 941м2т6 + 204мт7 + 17т8,

Q(m, т) = 2мт (144м7 + 524м6т + 968м5т2 + 1200м4т3 + 910м3т4+

+ 387м2т5 + 84мт6 + 7т7).

7.3. Случай т = м. Результат записывается в виде Л = P(м, v)/Q(m, v), где

P (м, v) = 256м10 + 2112м9 v + 8044m8v2 + 19355м7v3 + 32167m6v4+

+ 36887m5v5 + 28709m4v6 + 14854m3v7 + 4912m2v8 + 944mv9 + 80v10,

Q(M, v) = 2mv(m + v)(192m8 + 1344m7v + 4047m6v2 + 6770m5v3+

+ 6799m4v4 + 4216m3v5 + 1600m2v6 + 344mv7 + 32v8).

Литература

1. Baccelli F., Cohen G., Olsder G. J., Quadrat J.-P. Synchronization and linearity: An algebra for discrete event systems. Chichester: Wiley, 1992.

2. Маслов В. П., Колокольцов В. Н. Идемпотентный анализ и его применение в оптимальном управлении. М.: Физматлит, 1994.

3. Литвинов Г. Л., Маслов В. П., Шпиз Г. Б. Идемпотентный функциональный анализ. Алгебраический подход // Матем. заметки. 2001. Т. 69. Вып. 5. С. 758-797.

4. Olsder G. J., Resing J. A. C., De Vries R. E., Keane M. S., Hooghiemstra G. Discrete event systems with stochastic processing times // IEEE Trans. Automat. Contr. 1990. Vol. 35, N 3. P. 299-302.

5. Jean-Marie A. Analytical computation of Lyapunov exponents in stochastic event graphs // Performance Evaluation of Parallel and Distributed Systems. Solution Methods: Proc. 3rd QMIPS Workshop. Amsterdam: CWI, 1994. P. 309-341. (CWI Tracts, Vol. 106.)

6. Кривулин Н. К. Вычисление скорости роста вектора состояний для одной модели обобщенной линейной стохастической системы // Вестн. С.-Петерб. ун-та. Сер. 1. 2007. Вып. 3. С. 91-99.

7. Кривулин Н. К. Скорость роста вектора состояний обобщенной линейной стохастической системы с симметричной матрицей // Зап. науч. семин. ПОМИ. 2007. Т. 341. С. 134-141.

8. Кривулин Н. К. О вычислении скорости роста вектора состояний обобщенной линейной стохастической системы второго порядка // Вестн. С.-Петерб. ун-та. Сер. 1. 2008. Вып. 1. С. 38-48.

9. Kingman J. F. C. Subadditive ergodic theory // Ann. Probab. 1973. Vol. 1. P. 883-909. Статья поступила в редакцию 30 октября 2008 г.