close

Вход

Забыли?

вход по аккаунту

?

Алгоритмы численного решения стохастических дифференциальных систем с переключаемой диффузией.

код для вставкиСкачать
Управление большими системами. Выпуск 36
УДК 519.6+004.94
ББК 22.193
АЛГОРИТМЫ ЧИСЛЕННОГО РЕШЕНИЯ
СТОХАСТИЧЕСКИХ ДИФФЕРЕНЦИАЛЬНЫХ
СИСТЕМ С ПЕРЕКЛЮЧАЕМОЙ ДИФФУЗИЕЙ
Черных Н. В.1, Пакшин П. В.2
(Арзамасский политехнический институт (филиал)
Нижегородского государственного технического университета им Р.Е. Алексеева)
Рассматриваются математические модели сложных систем в
виде стохастических дифференциальных уравнений с марковскими переключениями диффузионной составляющей. Предлагается обобщение известных численных схем Тейлора для аппроксимации решения таких уравнений. Представлены результаты
численного моделирования в среде Scilab.
Ключевые слова: стохастические системы, диффузионные процессы, марковские переключения, схемы Тейлора, сходимость,
устойчивость.
1. Введение
Стохастические дифференциальные уравнения служат фундаментом для многих разделов прикладных наук, т.к. успешно
моделируют системы, функционирующие при наличии случайных возмущений. Они являются одним из основных объектов
современной теории управления. Большое значение, которое
приобрели стохастические дифференциальные уравнения, объяс-
1
Надежда Валентиновна Черных, аспирантка (nadezdacher@mail.ru).
Павел Владимирович Пакшин, доктор физико-математических наук,
профессор (pakshinpv@gmail.com).
2
106
Математическая теория управления
няется также их тесной связью с уравнениями математической
физики.
Один из известных и перспективных подходов к численному
интегрированию стохастических дифференциальных уравнений
Ито основан на стохастических аналогах формулы Тейлора.
Теоретически это разложение позволяет строить методы как
угодно высокого порядка точности (при соответствующих предположениях о коэффициентах системы уравнений).
Важнейшей отличительной особенностью стохастических
аналогов формулы Тейлора для решения стохастических дифференциальных уравнений Ито является присутствие в них, так
называемых, повторных стохастических интегралов в форме Ито
или Стратоновича, которые являются функционалами сложной
структуры относительно компонент векторного винеровского
процесса [3]. Аналогом обычных степеней винеровского процесса
служат в теории Ито многочлены Эрмита, а обычную экспоненту
от стохастического интеграла дублирует экспоненциальный
супермартингал. Проблема совместного численного моделирования
(исходя из среднеквадратического критерия сходимости) совокупностей повторных стохастических интегралов Ито является сложной
как с теоретической, так и с вычислительной точки зрения.
В последнее время внимание исследователей привлекают математические модели в виде стохастических дифференциальных
уравнений (СДУ) со скачкообразными изменениями диффузионной составляющей. Они получили название моделей с переключаемой диффузией и описывают сложные системы, которые
могут испытывать резкие изменения структуры и параметров,
вызванные возможными отказами, перерывами в поступлении
информации и воздействиями внешней среды. Такие модели
получают все более широкое распространение в современной
теории управления и информации.
Исследование скачкообразных систем началось с работ
I. YA. Kats, N.N. Krasovskii, E.A. Lidskii (1960–1961), а затем
T. Kazangey и D.D. Sworder (1971) представили скачкообразную
систему, где макроэкономическая модель народного хозяйства
использовалась, чтобы изучить эффект федеральной политики по
107
Управление большими системами. Выпуск 36
стабилизации
жилищного
сектора.
Известна
работа
Г.Н. Мильштейна (1972) «Mean square stability of linear systems
driven by Markov chain». A.S. Willsky и B.C. Levy рассматривали
скачкообразные системы для моделирования систем электроэнергии и управления солнечным тепловым центральным приемником. M. Mariton в своей работе «Jump Linear Systems in Automatic
Control» поясняет, что скачкообразные системы появились как
удобная математическая структура для моделирования разнообразных проблем в различных областях, таких как задачи слежения, контроля допустимых ошибок и производственных процессов. [11]
Как правило, рассматриваются переключения диффузии по
закону марковской цепи с конечным числом состояний – марковские переключения.
При компьютерном моделировании таких систем возникает
необходимость численного решения СДУ с марковскими переключениями.
Известны работы по изучению управляемости, устойчивости,
стабилизации таких систем. Вопросы численного решения СДУ с
марковскими переключениями изучены мало.
В работах [9, 10] рассмотрена известная численная схема
Эйлера (для стохастических систем) в применении к решению
стохастического дифференциального уравнения с марковскими
переключениями, приведены доказательства некоторых результатов частного характера.
В данной работе рассматривается возможность применения
численных схем, основанных на стохастических аналогах формулы Тейлора для аппроксимации решений таких уравнений.
Теоретические выводы подтверждаются проведенными компьютерными экспериментами в среде SCILAB.
2. Предварительные сведения
Пусть (Ω, F, P) – вероятностное пространство, Ft, t0 ≤ t ≤t0 + T
– неубывающее семейство σ – подалгебр F, (ωr(t), Ft), r = 1,…,d –
независимые винеровские процессы.
108
Математическая теория управления
Рассмотрим стохастическое дифференциальное уравнение
Ито
d
(1) dX  a(t , X )dt    r (t , X )dr (t ),
r 1
где X, a, σr - векторы размерности n. [4]
Предполагается, что функции a(t, x) и σr(t, x) определены и
непрерывны при t  [t 0 , t 0  T ] , x   n и удовлетворяют условию
Липшица при всех t [t0 , t0  T ] , x   n , y   n :
(2)
a(t , x)  a(t , y )   (t , x )   (t , y )  K x  y . [4]
Будем использовать далее следующие обозначения:
|x| означает евклидову норму вектора x, xy – скалярное произведение векторов x и y; К – положительная константа.
Xt,x(t) или просто X(t) - решение уравнения (1), удовлетворяющее
начальным данным Xt,x(0) = x. Разобьем промежуток [t0, t0 + T]
точками деления tk на N равных частей, так что tk+1 – tk = h, k = 0,
1, …, N - 1, t0 + T = tN, h = T/N. Приближение к X(tk) будем обозначать X (t k ) , где X 0  X (t0 ) . Далее пусть X – Ftk -измеримая
случайная величина и E|X|2 < ∞; X t k , x (t ) означает решение уравнения (1) для tk ≤ t ≤ t0 + T, удовлетворяющее начальным данным
при t = tk.
Пусть (X(t), F(t)), t0 ≤ t ≤ t0 + T – некоторое решение уравнения (1), с начальным условием, удовлетворяющим неравенству
E|X(t0)|2 < ∞.
Определим одношаговую аппроксимацию
X t , x (t  h) ,
t0 ≤ t ≤ t + h ≤ t0 + T, которая формируется в зависимости от x, t, h,
и 1 ( )  1 (t ),...,d ( )  d (t ) : t    t  h :
(3)
X t , x (t  h)  x  f (t , x, h; i ( )  i (t ), i  1,...d , t    t  h) ,
где функция f определяется конкретным методом нахождения
решения.
На основании одношаговой аппроксимации рекуррентно построим приближение
X k , Ft k , k = 0, …, N, tk+1 - tk = hk+1,


tN = t0 + T:
109
Управление большими системами. Выпуск 36
(4)
X 0  X 0  X t0 , X k 1  X tk , X k (t k 1 ) =
 X k  f tk , X k , hk 1; i ( )  i (t k ), i  1,...d , t k    t k 1  .
Для простоты считаем, что tk+1 - tk = h = T / N. [4]
Теорема 1. Пусть одношаговая аппроксимация X t , x (t  h)
имеет порядок точности ρ1 для математического ожидания
отклонения и порядок точности ρ2 для среднеквадратичного
отклонения, т. е., при любых t0 ≤ t ≤ t0 + T - h, x   n выполняются неравенства
(5)
E X t , x t  h   X t , x (t  h)   K (1  x )1/ 2 h 1
2
2
(6)  E X t , x t  h   X t , x (t  h) 


и пусть
(7)  2  1 / 2, 1   2  12 .
1/ 2
2
 K (1  x )1 / 2 h 2
Тогда при любых N и k = 0, 1, …, N выполняется неравенство
1


1
2 2
2 2  2  12
(8)  E X t0 , X 0 (t k )  X t0 , X 0 (tk )   K 1  E X 0
h
,


т.е. порядок точности метода, построенного с использованием
одношаговой аппроксимации X t , x (t  h) , равен    2  12 .
Постоянные K не зависят от X0 и N. [4]
3. Постановка задачи
Пусть (Ω, F, P) – вероятностное пространство, Ft,
t0 ≤ t ≤ t0 + T – неубывающее семейство σ – подалгебр F, ωr(·),
r = 1, …, d – независимые винеровские процессы. Пусть
M = {1, …, m} – конечное множество.
Рассмотрим стохастическое дифференциальное уравнение с
марковскими переключениями в форме
d
(9) dX (t )  a (  (t ), X (t ))dt    r (  (t ), X (t )) dr (t ) ,
r 1
где  (t ) - однородный марковский процесс со счетным множеством состояний M,  (0)  u0 , X (0)  x0 .
110
Математическая теория управления
(10) P  t  h   l  (t )  u, x ( s),  ( s), s  t   qul t h  oh , u  l ,
где x(t )   n , a, :  n  M   n и  , :  n  M   nn .
Переходная функция такого процесса определяется набором
функций P (t , u, l )  pul (t ) , образующих стохастическую матрицу
P(t) переходных вероятностей ( pul (t )  0,  pul (t )  1 ), где числа
l
pul имеют смысл вероятностей перехода из u в l за интервал h при
условии, что процесс β(t) на этом интервале покинул состояние u.
Эволюция процесса  (t ) целиком описывается с помощью
чисел qul, ql, образующих матрицу интенсивности переходов
P(t )  P(0)
Q (t )  (qul (t ))   mm , Q (t )  lim
, где для каждого t
t 0
t
m
qul (t )  0 при u  l , quu   qu ,
 qul (t )  0 для каждого u  M ; [1]
l 1
Предполагается, что функции a(β(t), x(t)) и σr(β(t), x(t)) определены и непрерывны при t  [t0 , t0  T ] , x   n и удовлетворяют
условию Липшица при всех t [t0 , t0  T ] , x   n ,
uM :
(11) a (u, x)  a(u , y )   (u , x )   (u , y )  K x  y ,
y n ,
а также условию:
(12) a (u, x)   (u , x )  K (1  x ) ,
которое накладывает ограничения на скорость изменения компонент функций по x.
При выполнении всех описанных выше условий уравнение
(9) имеет единственное, непрерывное решение Xu, x на интервале
t ≥ 0 для каждого начального значения. [8]
Поставим задачу численного решения уравнения (9), т. е. нахождения аппроксимации его решения на рассматриваемом
временном интервале.
При такой постановке проблемы, когда матрица переходов Q
не зависит от x, а зависит только от t, ее можно свести к известной задаче решения стохастического дифференциального уравнения вида (1) с помощью схем, описанных в [2, 4, 7] на совокуп111
Управление большими системами. Выпуск 36
ности случайных подынтервалов [0, t1 ), [t1 , t1  t 2 ),... ,где tk – случайные моменты переключения марковской цепи, которые можно
определить заранее. Так как следующий момент преключения tk и
состояние β(tk) есть случайные величины, распределение которых
зависит лишь от (β(tk), tk), можно предварительно сгенерировать
совокупность (β(tk), tk) правых частей уравнения (1) и решать
уравнение (1) известными методами. При этом точность аппроксимации, которая зависит от выбора шага на отдельных подынтервалах, будет достижима и на всем рассматриваемом временном интервале.
Но можно предложить для решения уравнения (9) другой
подход.
Для стохастического дифференциального уравнения вида (1)
Е. Платен предложил очень простой вывод (использующий лишь
формулу Ито) разложения решения Xt, x(t + h) по степеням h и
интегралам, зависящим от приращений  r ( )   r (t ) , где
t    t  h , r = 1, …, d. Это разложение в детерминированном
случае представляет собой формулу Тейлора для Xt, x(t + h) по
степеням h в окрестности точки (t, x). [2, 4, 7]
Будем использовать разложение Платена, подробно рассмотренное в [4], для решения уравнения (9), заменив функцию f(t, x)
функцией f(β(t), x(t)) с переключаемой компонентой.
Пусть Xu, x(s) = X(s) - решение уравнения (9) для всех u  M ;
f(β, x), где β = β(t), X = X(t) – достаточно гладкая функция (скалярная или векторная). Согласно формуле Ито имеем для
t0  t    t0  T
(13) f (  ( ), X ( ))  f (  , x) 
d 

r 1 t
t
    r f  (1 ), X 1 d r 1    Lf  (1 ), X 1 d1 ,
где операторы
 r , r = 1, …, d и L определены как:
 

(14)  r    r ,  ,
x 

112
Математическая теория управления
2
   1 q n n
(15) L   a,       ri  rj i j ,
x x
 x  2 r 1 i 1 j 1
Применим формулу (13) к функциям  r f и Lf, а затем полученные выражения для  r f  ( ), X ( )  и Lf  ( ), X  
подставим в (13). Имеем
d
s
s
r 1
t
t
(16) f  ( s ), x ( s)   f    r f  d r    Lf  d 
d s d 
       s  r f  (1 ), X (1 ) d s (1 )d r ( ) 
r 1 t  s 1 t
d s 
     L r f  (1 ), X 1 d1 d r   
r 1 t  t
d s 
      r Lf  (1 ), X (1 ) d r (1 ) d 
r 1 t  t
s 

    L2 f (  (1 ), X (1 ))d1 d ,
tt
где, например,  r f вычисляется в точке (β, x).
Поступая так дальше можно получить разложения для
f(β(t + h), X(t + h)), где роль степеней выполняют случайные
величины вида (которые не зависят от Ft)
t h
(17) I i1 ,...i j ( h) 

1
 j 2
 d i j ( )  di j 1 (1 )  ...  d i1 ( j 1 ),
t
t
t
t
где i1, …, ij принимают значения из множества чисел 0, 1, …, d,
t    1  2  ...   j 1  t  h и под d0 (r ) понимается dr .
Интеграл (17) соответствует j-степени параметра h в детерминированном разложении.
Очевидно, что EI i1 ,...,i j  0 , если хотя бы одно из ik ≠ 0,
k = 1, …, j, и EI i1 ,..., i j  O(h j ) , если все ik = 0, k = 1, …, j.
Рассмотрим примеры вычисления нескольких интегралов
Ито. Для этого будем использовать формулы (18) из теории
интегралов Ито ([4, 5, 6]):
113
Управление большими системами. Выпуск 36
t
(18)   ( ) d ( ) 
0
 2 (t )  t
;
2
3
t
2
  ( )d ( ) 
0
0
th
I1 (h ) 
t
 d ( )   (t  h)   (t )   (h);
t

th
t
t
t
th
t
t
th
d ( )  d1 (1 ) 

t h
 d  d1   d 
t

th

th
d 0 ( )  d 0 (1 ) 

I1,1 (h) 
0
 (t )
   ( )d ;
3
0
 d  h ;
I 0, 0 ( h ) 
t
t
th
(19) I 0 (h) 
t
 d ( )  t (t )    ( )d ;

t
 ( )d ( ) 
t
th


h2
;
2
 2 ( h)  h
;
2
t h
 d0 ( )  d (1 )   d  d (1 )    ( )d ;
I1, 0 (h) 
t
th
t
t
t h

t

t
t h
 d ( )  d0 (1 )   d ( )  d   d ( ) 
I 0,1 (h) 
t
t h
t
t
t
t
 h    ( )d  h (h)  I1, 0 ;
t
t h
I1,1,1 (h) 

t

1
t
t

t h


d ( )  d1 (1 )  d 2 (2 ) 
d ( )  1 (1 )d1 (1 ) 
t
t
t h

t
 2 ( )  
d ( ) 
2
3
 ( h) 1
 h (h) ;
6
2
В формулах (19) все рассмотренные повторные интегралы

th
Ито выражаются через случайные величины  (h) ,
  ( )d .
t
Вернемся к разложению (16) и выпишем формулу, которая
получается путем ряда непосредственных подстановок:
114
Математическая теория управления
d
(20) f  (t  h), X (t  h)   f    r f
r 1
d
r 1 i 1
t
d
   s  i  r f
r 1 i 1 s 1
1

t
t
d
t
th

 d  dr (1 )   Lf  dr ( ) d1 
   r Lf
r 1
t

t h
t
 d  d1   ,
t
d
t
 dr ( ) di (1 )  d s (2 ) 

t h
 L2 f
 d 
t
t h
d
d
t
 dr ( ) di (1 ) 
   i  r f
d

t h
d r ( )  Lf

t h
d
t h
r 1
t
t
где f  f (  (t ), x (t )) ,
t
d
d
d t  h   1 2
(21)     
   (  (   j  s  i  r f 
t t t
  (3 ), X (3 ) d j (3 ) d j (3 ) d s (2 ) di (1 ) d r ( ) 
r 1 i 1 s 1 j 1 t
d
d t  h    1
 
     L i  r f ( (2 ), X (2 ))d2 di (1 ) dr ( ) 
1
t t
d d t  h    1
        i L r f  (2 ), X (2 ) di (2 ) d1 d r ( ) 

r 1 i 1 t  t  t
d d t  h    1
        i  r Lf  (2 ), X (2 ) di (2 ) d r (1 ) d 

r 1 i 1 t  t  t
r 1 i 1 t
d
d
d t  h    1 2
 
     (  L s Li  r f ( (3 ), X (3 ))d3 )d s (2 )) 
t t t
 di (1 ) dr ( ) 
r 1 i 1 s 1 t
d t  h    1

     L  r f (  (2 ), X (2 ))d2 d1 dr ( ) 
r 1 t
2
t t
d t  h    1

     L r Lf  (2 ), X (2 )d2 dr (1 )d 
r 1 t
t t
115
Управление большими системами. Выпуск 36
d t  h    1

2
      r L f  (2 ), X (2 ) dr (2 )d1 d 
r 1 t

t t
t  h    1
3
     L f  (2 ), X (2 ) d2 d1 d.
t t
Рассмотрим леммы, доказанные в [4], применительно к разложению (20).
Лемма 1. Справедливо соотношение
 j 2ik 
1

2 2
2 
(22) E I i1 ,...,i j
 O h k 1  .




Другими словами, при подсчете порядка малости интеграла (17)
нужно руководствоваться правилом: d вносит в порядок
малости единицу, а d r ( ) , r = 1, …, d – одну вторую.
Лемма 2. Пусть
t



(23)  i j ... i1 f (  (t ), x (t ))  K 1  x
2

1
2
.
Тогда величина
(24) I i1 ,...,i j ( f , h) 

t h

1
 j 2
t
t
 di ( ) di (1 )  ...   i ... i f  ( j 1 ), X ( j 1 )di ( j 1 )
j
t
j 1
t
j
1
1
удовлетворяет неравенству
j
2

2

 ( 2ik )
(25) E I i1 ,...,i j ( f , h)  K 1  E X (t ) h k 1
,
т. е., в частности, ее порядок малости тот же, что и у величины I i 1 ,..., i j ( h) . Далее, если хотя бы один из индексов ik,
k = 1, …, j отличен от нуля, то
j
(26) EI i1 ,...,i j ( f , h)  0 ,
 ik2  0 .
k 1
Из этой леммы вытекает, что каждое слагаемое, составляющее ρ, имеет не более чем второй порядок малости. Более того,
математическое ожидание всех членов второго порядка малости
116
Математическая теория управления
и порядка малости 2,5 из ρ равно нулю согласно (26). Поэтому
|Eρ| = O(h3) Разумеется, это верно, если для всех подынтегральных функций из ρ выполняется, например, (23).
Леммы 1, 2 справедливы для разложения (20) в силу того,
что в формуле (20) интегралы Ито остались в том же виде, в
котором они используются в разложении Платена (см. [4, 7]).
На основании этих лемм непосредственно подсчитывается
порядок каждого слагаемого в разложениях вида (20).
Легко увидеть [4], что можно получить как разложения, содержащие все члены включительно до какого-то полуцелого
порядка малости, так и до какого-то целого порядка.
Возьмем теперь в (20), (21) в качестве f(β, x) вектор x. В этом
случае  r f   r , Lf = a. Поэтому
d
th
(27) X u , x (t  h)  x   r  d r ( )  ah 
r 1
d
t
t h
d
   i r
r 1 i 1
 i ( )  i (t ) dr ( ) 
t
t h
d
  L r
r 1
d
t h
r 1
t
 (  t )dr ( )    r a  r ( )  r (t )d 
t
h2





(

)


(
t
)
d

(

)
d

(

)

La
.
s
i
1
r
  s 1
2
r 1 i 1 s 1
t t
В формуле (27) все коэффициенты σr, a, Λi σr, Lσr, Λr a,
ΛS Λi σr, L a вычисляются в точке (β, x), а остаток ρ равен
d d d d t  h   1 2
(28)        (  (   j  s  i r (  (3 ), X (3 )) 
r 1 i 1 s 1 j 1 t  t t t
 d j (3 ))d 3 ( 2 )) d i (1 ) d r ( ) 
d
d
d
   s  i r
d
t h 
d t  h   1
 
   (  L i r (  (2 ), X (2 ))d2 ) di (1 ) dr ( ) 
t t
d d t  h   1
     (   i L r (  (2 ), X (2 ))di (2 )) d1 d r ( ) 
r 1 i 1 t  t t
r 1 i 1 t
117
Управление большими системами. Выпуск 36
d
d t  h   1
 
   (   i  r a(  (2 ), X (2 ))di (2 )dr (1 ) d 
t
r 1 i 1 t
d
d
t
d t  h   1 2
   (  (  L s  i r (  (3 ), X (3 ))d3 ) ds (2 )) 
 
t t t
 d i (1 ) d r ( ) 
r 1 i 1 s 1 t
d t  h   1
   (  L  r (  (2 ), X (2 ))d2 )d1 dr ( ) 
r 1
2

t t
d t  h   1
     (  L r a(  (2 ), X (2 ))d2 ) d r (1 )  
r 1 t  t t
d t  h   1
     (   r La(  (2 ), X (2 ))d r (2 )) d1 d 
r 1 t  t t
t  h   1
    (  L2 a (  (2 ), X (2 ))d2 ) d1 d.
t t t
В связи с формулами (27), (28) рассмотрим следующую одношаговую аппроксимацию:
t
d
(29)
 r  r (t  h)   r (t )   ah 
X u ,x (t  h)  x  
r 1
d
t h
d
   i r
r 1 i 1
d
th
 i ( )  i (t )dr ( )   L r  (  t )dr ( ) 
r 1
t
t
th
d
  r a
r 1
 r ( )  r (t )d 
t


h2




(

)


(
t
)

d

(

)
d

(

)

La
,
s
i
1 
r
  s 1
2
r 1 i 1 s 1
t t

где a = a(β(t), x(t)), σr = σr(β(t), x(t)).
Пользуясь леммами 1 и 2, можно установить (при выполнении условия (23) для соответствующих функций), что для одношаговой аппроксимации (29)
d
d
d
th 
    s  i r
E  O (h 3 ),
118
2
 
E   O h4 ,
Математическая теория управления
т. е. ρ1 = 3, ρ2 = 2, и этот метод имеет порядок сходимости равный
3 / 2 (на основании Теоремы 1).
Общий результат построения одношаговых аппроксимаций
обоснован и сформулирован в [4].
Все многократные интегралы Ито вида (17) имеют множителем ∆ω(h) (см. формулы (19)), т. е. величины, способные принимать сколь угодно большие значения. Тем не менее, схемы, основанные на разложении Платена, обеспечивают хорошую
аппроксимацию решений стохастических дифференциальных
уравнений, что подтверждается компьютерными экспериментами, например в [7].
Обобщим теорему Платена (приведенную в [4]) на рассматриваемый случай решения стохастического дифференциального
уравнения с марковскими переключениями (9). Все выводы
теоремы Платена остаются в силе для разложения (27), так как
выполняются леммы 1, 2.
Теорема 2 (обобщение теоремы Платена).
Пусть
X u , x (t  h) содержит все слагаемые вида  i j ... i1 f  I i1 ,...,i j ,
где f  x , до целого порядка m включительно. Пусть все функции
j
2  ik
 i j ... i1 f (  (t ), x (t )) , где f  x , 
 m  1 , удовлетворяют
2
k 1
неравенству (23). Тогда среднеквадратичный порядок точности
метода, основанный на такой аппроксимации, равен m.
Пусть
X u , x (t  h )
содержит все слагаемые вида
 i j ... i1 f  I i1 ,..., i j , где f  x , до полуцелого порядка m 
чительно
Пусть
и
все
слагаемое
функции
m
L a
t h

m 1
t
t
t
1
2
вклюm 1
h
m
 d  d1...  dm  L a (m  1)! .
 i j ... i1 f (  (t ), x (t )) ,
где
f x,
j
2  ik
 m  2 , удовлетворяют неравенству (23). Тогда сред2
k 1
неквадратичный порядок точности метода, основанный на
такой аппроксимации, равен m  12 .

119
Управление большими системами. Выпуск 36
Для нахождения решения уравнения (9) предложим следующую численную схему, построенную на одношаговой аппроксимации (29), т. е. среднеквадратичного порядка точности 3 / 2.
(30) X t k 1  X tk  a(  t k , X t k )h   (  tk , X tk )  k 
1
2
  (  tk , X tk ) (  t k , X t k )  k   h  a(  t k , X t k ) (  t k , X t k )Z k 
2
1
1

  a(  t k , X t k )a(  t k , X t k )   2 (  tk , X tk ) a(  tk , X tk ) h 2 
2
2

1


  a(  t k , X tk ) (  tk , X tk )   2 (  t k , X t k ) (  tk , X t k )  k h  Z k  
2


1
2
  (  tk , X tk )  (  tk , X tk ) (  t k , X t k )   (  t k , X t k ) 
2
1

2
   k   h  k ,
3






t k 1 s2
где    (t ) , x  x (t ) , Z k 
  d s ds2 ,
1
 k - есть N(0, h) гаус-
tk t k
совски распределенных приращений винеровского процесса ω на
подынтервалах t k  t  t k 1 , X t0  x0 .
Рассмотрим теперь модель переключения, зависящего от фазового состояния.
Вернемся к уравнению (9), где β(t) - однородный марковский
процесс со счетным множеством состояний M, β(0) = u0, X(0) = x0,
(31) P  t  h   l  (t )  u, x ( s),  ( s), s  t   qul x h  oh , u  l ,
где x(t )   n , a, :  n  M   n и  , :  n  M   nn .
(32) Пусть Q  :  n   mm - ограниченная и непрерывная функция. Q x   (qul ( x ))   mm для каждого x, qul ( x )  0 при u  l ,
m
 qul ( x)  0 для каждого u  M .
l 1
Такая состояние - зависимая модель находит значительно
большее применение в управлении и оптимизации. Если матрица
Q зависит от состояния x(t), то случай становится значительно
120
Математическая теория управления
более сложным. Одна из главных трудностей - та, что из-за непрерывной состояние - обусловленной связи, β(t) и x(t) становятся
зависимы. β(t) – является цепью Маркова только для фиксированного x, то есть фактически не является марковской. [9, 10]
В отличие от обычных диффузионных процессов, моделируемых стохастическими дифференциальными уравнениями,
распределение такой переключаемой диффузии имеет смешанный
характер.
Суть задачи в следующем - рассмотреть пару процессов β(t)
и x(t) совместно; т. е. процесс с двумя компонентами рассматривается как марковский. [9, 10]
Пусть переключаемая диффузия (β(t), x(t)) задана уравнением (9). Пусть первое начальное значение (β(0), x(0)) = (u0, x0),
также задано другое начальное значение (β(0), y(0)) = (u0, y0) для
y ≠ x. Для состояние - зависимого случая  u0 , x0 (t )   u0 , y0 (t )
бесконечно часто, даже при том, что первоначально
 u0 , x0 (0)   u0 , y0 (0)   , где верхний индекс показывает начальное значение зависимости.
Будем считать, что β(t) – стохастически непрерывный процесс, для которого  (t  h)   (t ) при h → 0. В соответствии с
q
требованием сепарабельности будем считать, что β(t) не может
менять свое состояние за «нулевое время» более одного раза. В
этом случае траектории процесса будут кусочно – постоянными,
т. е. временной интервал разбивается на подынтервалы
[0, t1), [t1, t1 + t2),…, на которых β(t) постоянно, tk – случайные
моменты переключения марковской цепи. [1]
Далее будем рассматривать последовательность  (t k ) как
дискретно - временной стохастический процесс, аппроксимирующий β(t) в соответствующем значении.
Используя функцию Q(x) можно построить матрицу переходных вероятностей P  e Q ( x ) h . Используя разложение функции
в ряд Тейлора I  Q ( x )  O( 2 ) и отбросив O (2 ) ввиду ограниченности и непрерывности Q(x), будем моделировать матрицу
121
Управление большими системами. Выпуск 36
переходных вероятностей как P  I  Q (x ) , где I – единичная
матрица.
Будем рассматривать пару процессов β(t) и x(t) совместно как
марковскую цепь следующим образом.
Значения X tk 1 генерируются рекурсивно согласно (30), используя предыдущие значение X t k , и одновременно генерируются значения  t k 1 , также используя значение X t k ( x  X t k в матрице переходных вероятностей P  I  Q(x ) ).
При всех вышеописанных условиях Леммы 1, 2 и обобщенная теорема Платена выполняются для случая марковского процесса, зависящего от фазового состояния, т. е. схема (30) будет
обеспечивать нахождение приближенного решения уравнения (9)
с матрицей переходов Q(x).
Подобный подход используется в [9, 10], где представлены
леммы и теоремы частного характера (рассматриваемая в этих
работах схема Эйлера является частным случаем схемы (30)).
Однако, при численном решении в модели с переключением,
которое зависит от фазового состояния, моменты переключения и
состояния после него, тоже определяются приближенно. При
этом возникает ряд вопросов по поводу определения того, как
приближенное решение отличается от точного. Возможно, нужно
ввести какую-то меру отклонения приближенного распределения
моментов скачков от точного распределения. Один из возможных
подходов - сделать замену меры, то есть превратить все скачки в
пуассоновские. Тогда проблема переходит в другую плоскость,
где надо оценивать, как отклоняется эта мера от точной, и насколько будут отличаться вычисленные средние значения функционалов от траекторий и точные значения. Эти вопросы пока
остаются открытыми.
4. Пример
Рассмотрим процесс Ито X = {Xt, t ≥ 0}, удовлетворяющий
линейному стохастическому уравнению
122
Математическая теория управления
(33) dX t  f (  t , X t ) X t dt  g (  t , X t ) X t dt
на временном интервале [0, T] при T = 1, X0 = 1.
Пусть βt - однородная марковская цепь, M = {1, 2} - число
состояний марковской цепи, Q – матрица интенсивности переходов, P – матрица переходных вероятностей.
f(βt, Xt) и g(βt, Xt) - непрерывные, ограниченные, дважды
дифференцируемые функции, удовлетворяющие условиям (11),
(12).


1

(34) X t  X 0 exp  f (  t , xt )  g 2 (  t , xt ) t  g (  t , xt )t 
2



есть решение уравнения (33) для t [0, T ] и данного винеровского
процесса   {t , t  0} .
Перепишем схему (30) в виде
(35) Yn 1  Yn  a(  n , Yn ) n   (  n , Yn ) I (1) 
  (  n , Yn ) (  n , Yn ) I (1,1)  a(  n , Yn ) (  n , Yn ) I (1,0 ) 
1


  a(  n , Yn )a(  n , Yn )   2 (  n , Yn ) a(  n , Yn )  I (0 ,0 ) 
2


1 2
 (a(  n , Yn ) (  n , Yn )   (  n , Yn ) (  n , Yn )) I (0 ,1) 
2
2
  (  n , Yn )  (  n , Yn ) (  n , Yn )   (  n , Yn )  I 1,1,1 ,


где a(βt, xt) = f(βt, xt) xt, σ(βt, xt) = g(βt, xt) xt для рассматриваемого
уравнения (33).
Для численного моделирования интегралов Ито введем в
рассмотрение независимые N(0, ∆) – распределенные случайные
величины
 n 1

1
3 
1
(36) 1   n2  n и  2  12 n2   ( ( )   ( )) d   n  n  .


2
 n

С помощью этих величин получим:
 n 1
1

1
3  1
(37)  n  2n1 ,
 ( ( )   ( ))d   n2  2 1  12  2  .


n
Тогда формулы (19) можно переписать в виде:
123
Управление большими системами. Выпуск 36
1
1
 n2   n ;
I 0, 0   2n ;
2
2


1
1
I 1,0  
 n  n   1 
 2 ; I 0,1   n  n  I 1, 0  ;
2
3 

1
1 1

I 1,1,1    n2   n   n ; где  n  2n1 ;
23

Зададим начальные значения Y0 = X0, β0 = u0 и будем рекурсивно генерировать 100 значений Yn с равным значением шага ∆
согласно (35), где Δn - есть длина временного интервала дискретизации t0 = τ0 < τ1 < … < τn < … < τN = T на временном интервале
[t0, T];
Для сравнения будем использовать (34), чтобы определить
соответствующие значения точного решения, используя ту же
примерную траекторию винеровского процесса ωt на подынтервалах τn < t < τn + 1, а именно
n


1

(38) X n1  X 0 exp  f (  n , X n )  g 2 (  n , X n )  n  g ( n , X n ) i 
2

i 1


Рассмотрим результаты численного решения уравнения (33),
выбирая различные варианты задания матрицы переходов P,
шага дискретизации ∆, значений функций f(βt, xt), g(βt, xt).
Случай 1.
I 1   n ;
I 1,1 


Пусть f(βt, xt) принимает два значения – {α1, α2}, соответствующие первому и второму состоянию марковской цепи. Положим α1 = 1,5, значение α2 будет изменяться в ходе эксперимента.
g(βt, xt) принимает два значения - {λ1, λ2}. Положим λ1 = 0,1,
значение λ2 будет изменяться в ходе эксперимента.
Будем полагать, что начальному моменту времени соответствуют значения α1 и λ1.
 0,2 0,8 
 ;
1. P = 
 0,4 0,6 
124
α2 = 10α1;
λ2 = 50λ1;
Математическая теория управления
Рис. 1. ∆ = 0,01 ( t  [0;0,1] )
Рис. 2. ∆ = 0,0001 ( t  [0;0,01] )
125
Управление большими системами. Выпуск 36
Рис. 3. Марковская цепь
В данном случае значения коэффициентов меняются не
очень значительно, поэтому устойчивая аппроксимация решения
уравнения наблюдается с размером шага дискретизации ∆ = 0,01.
 0,3 0,7 
 ; α2 = 100α1;
2. P = 
 0,6 0,4 
λ2 = 200λ1;
Рис. 4. ∆ = 0,001 ( t  [0;0,01] )
126
Математическая теория управления
Рис. 5. ∆ = 0,00001 ( t  [0;0,0001] )
В данном примере оба коэффициента увеличиваются значительно (в 100 и 200 раз), поэтому выбор меньшего шага дискретизации заметно повышает точность вычислений.
Случай 2.
f(βt, xt) принимает два значения - {α1, α2}, соответствующие
первому и второму состоянию марковской цепи. g(βt, xt) принимает два значения - {λ1, λ2}. Будем полагать, что начальному
моменту времени соответствуют значения α1 и λ1.
 0,2 0,8 
 ; α1 = cos x + sin x;
3. P = 
α2 = 1 + cos x;
 0,4 0,6 
λ1 = 0,5x;
λ2 = 0,2x;
127
Управление большими системами. Выпуск 36
Рис. 6. ∆ = 0,01 ( t  [0;0,1] )
Рис. 7. Марковская цепь
 0,6 0,4 
 ;
P = 
 0,4 0,6 
α1 = cos x + sin x;
α2 = 1 + cos x;
λ1 = 0,5x;
λ2 = 2x;
4.
128
Математическая теория управления
Рис. 8. ∆ = 0,01 ( t  [0;0,1] )
Рис. 9. ∆ = 0,001 ( t  [0;0,01] )
129
Управление большими системами. Выпуск 36
Рис. 10. Марковская цепь
Варианты 3, 4 численного моделирования отличаются матрицей переходов P и значением λ2.
Случай 3.
  50 cos 3 x 50 cos 3 x 
 ; P = I + Q ∆;
5. Q = 
 10 sin x 
 10 sin x
α1 = 2 + sin x;
λ1 = 0,5;
130
α2 = 1 + sin x cos x;
λ2 = 0,1;
Математическая теория управления
Рис. 11. ∆ = 0,01 ( t  [0;0,1] )
Рис. 12. Марковская цепь
131
Управление большими системами. Выпуск 36
Рис. 13. ∆ = 0,001 ( t  [0;0,01] )
Рис. 14. Марковская цепь
132
Математическая теория управления
  5 cos 2 x
5 cos 2 x 
 ; P = I + Q ∆;
Q = 
2
2 
 10 cos x  10 cos x 
α1 = 20 + cos2 x;
α2 = sin x + 5cos 2x;
λ1 = 0,2;
λ2 = 3;
6.
Рис. 15. ∆ = 0,001 ( t  [0;0,01] )
Рис. 16. Марковская цепь
133
Управление большими системами. Выпуск 36
Рис. 17. ∆ = 0,0001 ( t  [0;0,001] )
Рис. 18. Марковская цепь
134
Математическая теория управления
Эксперимент показывает, что в случае состояние - зависимой
модели переключения предложенная численная схема также
позволяет найти приближенное решение уравнения (33).
α1 = 20 + cos2 x; α2 = sin x + 5 cos 2x;
λ1 = 0,2;
λ2 = 0,3;
Этот вариант данных отличается от предыдущего значением λ2.
Рис. 19. ∆ = 0,01 ( t  [0;0,1] )
Рис. 20. Марковская цепь
135
Управление большими системами. Выпуск 36
Рис. 21. ∆ = 0,001 ( t  [0;0,01] )
Рис. 22. Марковская цепь
136
Математическая теория управления
Рис. 23. ∆ = 0,0001 ( t  [0;0,001] )
Рис. 24. Марковская цепь
137
Управление большими системами. Выпуск 36
  5 cos 2 x
5 cos 2 x 
;
7. P = 
2
2 
 10 cos x  10 cos x 
α1 = 3 + cos x;
α2 = 20 – sin x;
λ1 = 0,5;
λ2 = 0,1;
Рис. 25. ∆ = 0,01 ( t  [0;0,1] )
Рис. 26. Марковская цепь
138
Математическая теория управления
Рис. 27. ∆ = 0,001 ( t  [0;0,01] )
Рис. 28. Марковская цепь
139
Управление большими системами. Выпуск 36
Рис. 29. ∆ = 0,0001 ( t  [0;0,001] )
Рис. 30. Марковская цепь
140
Математическая теория управления
5. Заключение
Теорему Платена можно обобщить на рассматриваемый случай диффузионной системы с марковскими переключениями, что
подтверждает сходимость схем, основанных на стохастических
аналогах формулы Тейлора к решению стохастических дифференциальных уравнений с марковскими переключениями при
наложении определенных условий.
Возникающие вопросы устойчивости схем можно решить,
учитывая, что во все слагаемые в той или иной степени входит
множителем h. При достаточно малом h устойчивость схем будет
обеспечена. Эксперименты показали, что для многих случаев
подходящим является шаг дискретизации h = 0,01, но для случаев
сильного колебания значений функций a(βt, Xt), σ(βt, Xt) в уравнении (9) устойчивость будет обеспечивать выбор более малого
шага h = 0,001, h = 0,0001.
В случае состояние - зависимой модели переключений остается открытым вопрос – что считать точным решением и как
определять меру отклонения от получаемых приближенных
значений.
Также необходимо отметить, что для практических применений важнее оказывается не поиск решения, как реализации единственной траектории, хотя бы и с заданной точностью, а скорее
нахождение значений функционалов от случайных траекторий,
таких как математическое ожидание интегральных функционалов
и/или вероятностей пребывания в заданном множестве. Нахождение этих величин с заданной точностью методом Монте-Карло
требует предварительной оценки числа итераций, которое также
должно быть связано с точностью аппроксимации.
141
Управление большими системами. Выпуск 36
Литература
БОРОВКОВ А.А. Теория вероятностей: Учеб. пособие для
вузов.- 2-е изд, перераб. и доп. - М.: Наука. Гл. ред. физ.-мат.
лит., 1986. – 432 с.
2. КУЗНЕЦОВ Д.Ф. Стохастические дифференциальные уравнения: теория и практика численного решения. - Спб.: Издательство Политехнического университета, 2007. – 800 с.
3. КУЗНЕЦОВ Д.Ф. Стохастические дифференциальные уравнения: теория и практика численного решения // Дифференциальные уравнения и процессы управления. -2008. - № 1.
4. МИЛЬШТЕЙН Г.Н. Численное интегрирование стохастических дифференциальных уравнений. - Свердловск: Изд-во
Уральского университета, 1988. – 224 с.
5. МИЛЛЕР Б.М., ПАНКОВ А.Р. Теория случайных процессов
в примерах и задачах. – М.: ФИЗМАТЛИТ, 2002. – 320 с.
6. ОКСЕНДАЛЬ Б. Стохастические дифференциальные уравнения. Введение в теорию и приложения: Пер. с англ. – М.:
Мир, ООО «Издательство АСТ», 2003. – 408 с.
7. KLOEDEN P. E., PLATEN E., SCHURZ H. Numerical Solution
of SDE Through Computer Experiments. – Berlin: Springer –
Verlag, 1994.
8. MAO X. Stability of stochastic differential equations with Markovian switching // Stoch. Proce. Appl. – 1999. – Vol. 79. – 45-67.
9. YIN G., MAO X., YUAN C., AND CAO D. Approximation
methods for hybrid diffusion systems with state-dependent
switching processes: numerical algorithms and existence and
uniqueness of solutions// SIAM Journal on Mathematical Analysis – 2010. – Vol. 41, №6. - p. 2335-2352.
10. YIN G., ZHU C. Stochastic modeling and applied probability.
Hybrid switching diffusions. Properties and applications. //
Springer Science + Business Media, LLC 2010.
11. YUAN C., LYGEROS J. Stochastic markovian switching hybrid
processes // Project IST-2001-38314, COLUMBUS. Design of
Embedded Controllers for Safety Critical Systems, June 17, 2004.
1.
142
Математическая теория управления
ALGORITHMS FOR NUMERICAL SOLUTION OF
STOCHASTIC DIFFERENTAL SYSTEMS WITH
SWITCHING DIFFUSION
Nadezda Chernykh, post-graduate student, Arzamas Polytechnic
Institute of R.E. Alekseev, Nizhny Novgorod State Technical University, 19, Kalinina Street, Arzamas, 607227, Russia,
(nadezdacher@mail.ru).
Pavel Pakshin, Dr.Sci., professor, Arzamas Polytechnic Institute of
R.E. Alekseev, Nizhny Novgorod State Technical University, 19,
Kalinina Street, Arzamas, 607227, Russia, (pakshinpv@gmail.ru)
Abstract: Mathematical models are considered of hybrid systems in
the form of stochastic differential equations with Markovian switching of the diffusion component. An extension of Taylor schemes for
numerical approximation of their solutions is proposed. Modeling
results in SCILAB are presented to demonstrate efficiency of the
obtained algorithms.
Keywords: stochastic systems, diffusion process, Markovian
switching, Taylor schemes, convergence, stability.
Статья представлена к публикации
членом редакционной коллегии А. В. Добровидовым
143
Документ
Категория
Без категории
Просмотров
6
Размер файла
1 828 Кб
Теги
диффузией, решение, дифференциальной, алгоритм, система, стохастических, переключаемых, численного
1/--страниц
Пожаловаться на содержимое документа