Линейные детерминированные и стохастические системы

1 Линейные детерминированные и стохастические системы Лекция 15. Линейное преобразование случайных процессов. Качество систем при случайных воздействиях. Стохастические дифференциальные уравнения. Белый шум. Преобразование стационарного случайного процесса линейной системой. Точность LTI моделей систем при воздействии случайных возмущений. Системы дифференциальные уравнений со случайной правой частью. Стохастический интеграл и формула Ито. Линейные стохастические дифференциальные уравнения. Приложение. Стохастические интегралы Ито и Стратоновича. Численное моделирование линейных стохастических систем. Функции среды Matlab для моделирования случайных величин. Белый шум. Белым шумом nt называют стационарный случайный процесс с постоянной спектральной плотностью S n ( )  S0 . Тогда корреляционная функция определяется из следующих соотношений: Rn (t )      it it  S n ( )e d  S0  e d  2S0 (t ) . Таким образом, корреляционная функция белого шума с точностью до постоянного множителя представляет собой  - функцию Дирака. Очевидно, что дисперсия белого шума Dn   . Физический процесс типа белого шума реализовать невозможно, так как мощность этого процесса должна быть бесконечно большой. Однако такой процесс является удобной математической моделью, с помощью которой можно аппроксимировать входной сигнал системы в определенном диапазоне частот. Понятие приближения случайного процесса к белому шуму процесса достаточно условно. Один и тот же случайный процесс может считаться белым шумом по отношению к одной динамической системе и «цветным» по отношению к другой, имеющей более широкую полосу пропускания. Если полоса пропускания системы не достаточно узкая по отношению к полосе равномерности спектра возмущающего процесса  (t ) , то аппроксимация этого процесса белым шумом будет не адекватна постановке задачи стохастической динамики. В этих случаях, необходим учет свойств «цветности» процесса  (t ) . Преобразование стационарного случайного процесса линейной системой. Спектральная плотность установившегося решения обыкновенного дифференциального уравнения с постоянными коэффициентами. Не вдаваясь в понятие производной стационарной случайной функции  t , примем ее существование. Примем также, что ее производная t также будет стационарной, а ее 2 корреляционная функция будет определяться равенством R   d 2 R (t ) t  дифференцирование соотношения R (t )   S ( )e it dt 2 . Тогда d по параметру t позволяет получить   равенство: R (t )   2 S ( )eit d . С другой стороны, вследствие стационарности процесса t ,   имеет место соотношение R (t )   S ( )e it d . Отсюда получим следующее равенство:  S ( )   2 S ( ) . Для того, чтобы существовала спектральная плотность S ( ) необходимо, чтобы выполнялось  следующее неравенство:  2 S ( )eit d   , то есть свойства произведения  2 S ( ) позволяют  судить о дифференцируемости процесса  t . Как дальше будет показано, если интеграл от этого произведения ограничен, то процесс дифференцируем. Если же интеграл будет не ограничен, то процесс  t не имеет производной / 1 /. Для производных высшего порядка, в случае их существования, будут справедливо следующее соотношение: S ( n ) ( )   2n S ( ) . Рассмотрим теперь следующее линейное выражение: X t  bm d mt dt m  bm 1 d m 1t dt m 1  ...  b0t  Pm ( d )t . Отсюда dt получим следующее выражение для корреляционной функции процесса X t :  R X (t )   | bm (i ) m  bm 1 (i ) m 1  ...  b0 |2 S ( )eit d . Тогда несложно получить соотношение,  определяющее спектральную плотность процесса X t : S X ( ) | bm (i ) m  bm 1 (i ) m 1  ...  b0 |2 S ( ) | Pm (i ) |2 S ( ) . Рассмотрим теперь дифференциальное уравнение вида: Pn ( d d ) X t  Pm ( ) t . Если dt dt предположить, что  t стационарная случайная функция и решение уравнения также обладает свойством стационарности, то можно сразу записать спектральную плотность S X ( ) через спектральную плотность S ( ) : | Pn (i ) | X t | Pm (i ) | S (t ) или S X ( )  2 Уравнение Pn ( 2 | Pm (i ) |2 | Pn (i ) |2 S ( ) . d d ) X t  Pm ( ) t будет иметь стационарное решение в том случае, когда момент dt dt времени, для которого отыскивается решение этого уравнения, достаточно велик. То есть, это такой момент времени, когда все переходные процессы, связанные с начальными условиями затухнут (решение устойчиво). 3 Частотная характеристика стационарного линейного преобразования. Пусть стационарная случайная функция  t имеет спектральное разложение вида: t   e it dZ (d ) , где Z (d ) некоторая случайная мера с некоррелированными приращениями,  имеющая спектральную функцию F ( ) . Определение.  Если случайная функция  t , t  R , допускает представление вида: t  e it ( ) Z (d ) , где  комплексная функция ( ) удовлетворяет условию   | ( ) | 2 dF ( )   , то говорят, что  случайная функция  t получена из случайной функции  t с помощью стационарного линейного преобразования. Функция ( ) называется частотной характеристикой этого преобразования. Несложно  убедиться, что тогда будут справедливы следующие соотношения: D   | ( ) | 2 dF ( ) ,  F ( )    |  (v ) | 2 dF (v); , v  R . Если же, спектральная функция случайной функции  t будет  иметь спектральную плотность S ( ) , то случайная функция  t будет иметь спектральную плотность S ( ) и, эти спектральные плотности связаны равенством: S ( ) | ( ) |2 S ( ) .  Процесс  t будет стационарным, если выполняется неравенство  | ( ) | 2 dF ( )   . Тогда  корреляционная функция процесса  t определяется с помощью следующего соотношения:  e R (t )  it | ( ) |2 dF ( ) .  Рассмотрим теперь стационарное линейное преобразование во временной области. Пусть задана  функция w(t ) такая, что  | w(t ) | dt   . Рассмотрим следующее с.к. интегральное преобразование   стационарной центрированной случайной функции t   w(t  s) s ds . Очевидно, что будут  справедливы следующие соотношения: M {t }    w(t  s)M { s }ds  const ,  M {| t | }  2      w(t  s1 ) R (s2  s1 )w(t  s2 )ds1ds2  D (  | w(s) | ds)    2 . 4  Если функция  t имеет спектральное разложение  t  e it dZ (d ) , то тогда будет также  выполняться следующее соотношение: t     ( )   w(t  s)(  e it Z (d ))ds   e it ( ) Z (d ) , где    it  e w(s)ds . Причем будет справедливо неравенство  Следовательно, преобразование t     | ( ) | 2 dF ( )  M {| t |2 }   .    w(t  s) s ds является стационарным линейным  преобразованием с частотной характеристикой ( )   e  it w( s)ds .  Пусть теперь стационарный процесс  t является m раз с.к. дифференцируем. Пусть также существует n раз с.к. дифференцируемый процесс  t , связанный с  t линейным n дифференциальным уравнением вида:  akt(k )  k 0 m  bk t(k ) . Найдем условия, при которых процесс k 0  t будет стационарным случайным процессом, а также определим частотную характеристику преобразования  t в  t . Обозначим ( )  B(i ) . Если существует спектральная плотность S ( ) , то тогда, как было A(i ) показано выше, будет справедливо соотношение S ( ) | ( ) |2 S ( ) . По построению случайный процесс  t существует и является стационарным, если его дисперсия будет конечной, то есть M {| t |2 }   B(i )  | A(i ) | 2 dF ( )   . Данное соотношение будет выполняться, если полином  A(s) не будет иметь положительных и чисто мнимых корней. Теперь рассмотрим процесс прохождения через физически реализуемый линейный фильтр с частотной характеристикой ( ) стационарного процесса nt с постоянной спектральной плотностью, равной S n ( )  1 . Спектральная плотность S ( ) процесса  t на выходе фильтра, тогда определяется следующим образом: S ( ) | ( ) |2 , где ( )   e  it w( s)ds , w(t ) -  импульсная функция фильтра, удовлетворяющая условию абсолютной интегрируемости   | w(t ) | dt   .  Можно поставить обратную задачу: когда спектральную плотность некоторого стационарного процесса можно представить в виде результата прохождения белого шума nt со спектральной плотностью S n ( )  1 через физически реализуемый фильтр, то есть в дробно-рациональной форме. Ответ дается следующей теоремой. 5 Теорема. Для того, чтобы неотрицательная функция S ( ) допускала представление S ( ) | ( ) |2 ,  необходимо и достаточно, чтобы выполнялось условие   ln( S ( )) 1  2 d   . Точность LTI моделей систем при воздействии случайных возмущений. При воздействии на систему случайных процессов, как управляющих, так и возмущающих рассматривают вероятностные показатели качества. К ним обычно относят следующие оценки: - средний квадрат ошибки; - вероятность превышения сигналом ошибки заданного значения; - показатель накопления ущерба; - показатель правдоподобия полученной оценки ошибки; - число выходов ошибки из заданного интервала и другие. Кроме того, часто используют понятие критерия эффективности, как вероятность получения заданного результата. В данном разделе будем рассматривать только оценки ошибки функционирования замкнутой линейной стационарной системы при воздействие на нее случайного сигнала. Обычно применяют следующие показатели: - среднеквадратичная ошибка, которая определяется с помощью следующего соотношения: De   e2  M {(e(t )  me (t )) 2 } ; - вероятность превышения ошибкой заданного значения P{| e(t )  me (t ) |  } . Общая ошибка состоит из ошибки, обусловленной неточной «отработкой» системой полезного сигнала, и ошибки, обусловленной отработкой сигнала случайной помехи. В частотной области это соотношение можно графически изобразить следующим образом. 6 Характеристики LTI модели системы в переходном режиме. Пусть на вход линейной стационарной системы с импульсной переходной функцией w(t ) , где w(t )  L1{Wy, v (s)} , yt - выход системы, подается случайный процесс vt  u(t )  t , где u (t ) неслучайная функция, а  t - случайная составляющая. Физическая реализуемость системы  заключается в том, что  | w(t ) | dt   и w(t )  0; t  0 ,. В данном случае реализуемость системы  заключается в том, что система устойчива. Рассмотрим, для простоты, случай, когда u (t )  0 . Процесс yt связан с процессом  t с помощью следующего линейного оператора t yt   w(t   ) d . Таким образом, осуществляется линейное преобразование случайного процесса  t . Пока переходные процессы в системе не завершились, такое преобразование будет нестационарным. Зная импульсную характеристику w(t ) , математическое ожидание M { t } и корреляционную функцию R (t1 , t 2 ) входного сигнала, можно найти математическое ожидание M { yt } и корреляционную функцию R y (t1 , t 2 ) выходного сигнала с помощью следующих соотношений: t t1 t 2 00 M { yt }  m y (t )   w(t   )m ( )d , R y (t1 , t 2 )    w(t1   1 ) w(t 2   2 ) R ( 1 , 2 )d 1d 2 . Если, входной случайный процесс является стационарным в широком смысле, то для корреляционной функции выходного сигнала будет справедливо следующее соотношение: R y (t1 , t 2 )  t1 t 2   w(t1  1 )w(t2   2 ) R (1   2 )d1d 2 . 0 0 Характеристики сигнала ошибки, такие как: M {et } и Re (t1 , t 2 ) , будут определяться по аналогичным формулам. Однако импульсная функция будет определяться передаточной функцией We, (s) . Дисперсия ошибки будет определяться функцией De (t )  Re (t , t ) . Пример. Рассмотрим LTI модель системы, заданную передаточной функцией W0 ( s)  2( s  1) . На вход системы s  2s  3 подается смесь полезного сигнала с помехой: vt  u (t )  t , где u(t )  0.2 cos(0.5t ) - полезный сигнал, а помеха t представляет собой стационарный случайный процесс с моментами: 2 m (t )  0.05 , R (t )  0.01e |t | cos(t ) . Вычислить первый и второй моменты выхода системы на интервале [0,15] с. Математическое ожидание выхода Корреляционная функция равна t yt равно M { yt }   w(t   )u ( )d   w(t   )m ( )d . R y (t1 , t 2 )  t1 t 2   w(t1  1 )w(t2   2 ) R (1   2 )d1d 2 , а дисперсия 0 0 R y (t , t ) . t 7 Файл E74_TransRn. Характеристики LTI модели системы в установившийся режим системы. Для анализа процессов в устойчивой LTI модели системе при установившемся режиме будем предполагать, что начальный момент подачи входного воздействия происходит при t0   . Тогда t выход системы будет определяться соотношением yt   w(t   ) d ,  а математическое 8 ожидание выходного сигнала соотношением: M { yt }  t    w(t   )M { }d   w(s)m (t  s)ds . Такое линейное преобразование стационарного случайного процесса  t будет порождать на выходе также стационарный процесс yt , который будет иметь постоянные первый и второй моменты.   Следовательно, будет справедливо равенство: m y (t )  m w( s )ds . Соответственно, корреляционная функция выхода будет удовлетворять соотношению:  R y (t )    w(1 )w( 2 ) R (t  1   2 )d1d 2 . 00 Так как линейное преобразование является стационарным, то спектральная плотность S ( ) случайного процесса  t на входе системы будет связана со спектральной плотностью S y ( ) сигнала на выходе следующим равенством: S y ( ) | W y, (i ) |2 S ( ) , где W y , (i ) частотная характеристика линейного преобразования. Ковариационная функция и дисперсия выходного сигнала, при этом, удовлетворяют следующим соотношениям:  R y (t )   S y ( )e it d     | Wy, (i ) | 2 S ( )e  it d , D y     S y ( )d   | Wy, (i ) |  2 S ( )d .  Таким образом, выходной сигнал линейной системы представляет собой стационарный случайный процесс тогда и только тогда, когда рассматривается устойчивая стационарная линейная система в установившемся режиме при стационарном случайном входном сигнале. Обозначим среднее значение систематической ошибки в установившемся режиме от случайного задающего воздействия vt  u (t )  t , где u(t )  0; | u(t ) |  , как meu (t ) . Очевидно, что систематическая ошибка будет удовлетворять соотношению: meu (t )  eu (t )  me, (t ) , где eu (t ) составляющая систематической ошибки от неслучайной функции u (t ) , а me, (t ) - математическое ожидание установившейся ошибки от случайной функции  t . Функция me, (t ) будет определяться   формулой me, (t )  w(t   )m ( )d . Если на систему будет также воздействовать случайное стационарное возмущение f t , то тогда значение систематической ошибки будет равно meu (t )  eu (t )  me, (t )  me, f (t ) , где функция me, f (t )  будет определяться формулой me, f (t )   w(t   )m f ( )d . При этом для средних значений составляющих установившейся ошибки от случайных воздействий можно использовать те же выражения, с использованием коэффициентов ошибок, что и для детерминированного варианта функционирования системы: me, (t )  C ,0 m (t )  C ,1m (t )  C ,2 m (t )  ... , где C ,0  Wev (0), C , k  1 d kWev ( s) |s  0 ; k! ds k me, f (t )  C f ,0 m f (t )  C f ,1m f (t )  C ,2 m f (t )  ... , где C f ,0  Wef (0), C f , k k 1 d Wef ( s)  |s  0 . k! ds k 9 Если система обладает астатизмом r - го порядка, то коэффициенты ошибок C , k , при k  1,2,.., r  1 будут удовлетворять равенствам C ,0  0,..., C , r 1  0 , а коэффициент C , r можно определить с помощью соотношения: C , r  Wev ( s) sr |s  0 . То есть, упрощенными формулами можно пользоваться до первого, отличного от нуля, коэффициента включительно. Это же утверждение будет справедливо и для случайного возмущения f t , если система обладает по отношению к нему соответствующей степенью астатизма. Рассмотрим случай, когда внешние воздействия стационарны и не коррелированны. Тогда формулы для дисперсии и среднеквадратической установившейся ошибки имеют следующий вид:   Se ( )d  De  Re (0)   | Wev (i ) | 2 S ( )d ,  Se ( )d   | Wev (i ) | 2 S ( )d .   De  Re (0)       Если управляющее воздействие или возмущающее воздействие не коррелированны, то дисперсия ошибки будет равна сумме дисперсий: De  De  Def . Соответственно, формула для вычисления среднеквадратичной ошибки примет вид:  e  De  Def . Пример. Рассмотрим LTI модель разомкнутой системы, заданную следующей структурной схемой: 2s  1  s 1 , W2 ( s)  . Задающее воздействие имеет вид: g (t )  0.5t . s(0.5s  1) s 2  2s  3 Возмущение является стационарным случайным процессом с математическим ожидание m f  0.2 и Здесь W1 ( s)  2 корреляционной функцией R f (t )   2e  |t | cos(t ) , где  2  0.09,  1,   2 . Определить систематическую и среднеквадратическую ошибки в установившимся режиме. Очевидно, что имеют место следующие соотношения meg  Cg ,0 g (t )  Cg ,1g (t ), g(t )  ...  0 и mef  C f 0 m f . Так как система является астатической первого порядка (типа 1) , то отсюда получим, что C g ,0  0; C f 0  0 , а C g ,1  Weg ( s) s |s  0 . Таким образом, me  C g ,1  0.5 . Спектральная плотность 1 возмущения будет равна S f ( )  2  e  it   2 2   2   2 R f (t )dt   ( 2   2   2 ) 2  4 2 2 также вычислена численно с использованием функции fft (.) среды Matlab. Дисперсия ошибки тогда будет  определяться соотношением De  Def   | Wef (i ) |  равна  e  De . / 1 / и, может быть 2 S f ( )d . Среднеквадратическая ошибка будет 10 Файл E75_StdRn. Понятие формирующего фильтра стационарного случайного процесса. Подадим на вход устойчивой линейной стационарной системы входной стационарный случайный процесс, спектральная плотность которого зададим, для простоты, в следующем виде: S ( )  1 . Тогда, для установившегося режима системы, можно записать следующую формулу S y ( ) | W (i ) |2 . Таким образом, спектральная плотность выходного сигнала представляет собой дробнорациональную функцию, которая является действительной, четной и неотрицательной функцией параметра (частоты)  . Такую функцию можно представить в виде S y ( )  S y ( ) S y ( ) , где функция S y ( ) содержит нули и полюса функции S y ( ) , расположенные в верхней полуплоскости; функция S y ( ) содержит нули и полюса функции S y ( ) , расположенные в нижней полуплоскости. 11 Если ввести комплексный параметр s  i , то функция S y (s ) будет содержать нули и полюса рациональной функции S y ( ) в правой части s – плоскости, а функция S y (s) в левой части s – плоскости. Это следует из того, что для действительных значений  справедливо следующее равенство: S y (i )  S y (i ) . Такая процедура называется факторизацией спектральной плоскости. Указанным свойством обладает рациональная передаточная функция устойчивой линейной стационарной минимально-фазовой системы. Тогда, учитывая равенство S y ( ) | W (i ) |2 , можно записать, что W (i )  S y (i ) . Таким образом, разложив спектральную плотность формируемого сигнала на комплексно-сопряженные множители, можно определить передаточную функцию формирующего фильтра. Пусть заданы следующие равенства: W ( s)  bm s m  bm 1s m 1  ...  b0 an s n  an 1s n 1  ...  a0 , S y ( ) | W (i ) |  2 Bm 2m  Bm 1 2( m 1)  ...  B0 An 2n  An 1 2( n 1)  ...  A0 Тогда несложно получить следующие зависимости: Пример. Задана спектральная плотность стационарного случайного процесса t в виде: S ( )  2   2  4  2 2 2   4 . Провести факторизацию спектральной плотности и построить передаточную функцию формирующего фильтра. Факторизация спектральной плотности имеет вид: S ( )  S ( ) S ( )  передаточная функция формирующего фильтра имеет вид: W (s)  s (s   ) 2 . (  i ) (  i ) 2  (  i ) (  i ) 2 . Отсюда 12 Системы дифференциальные уравнений со случайной правой частью. Введенные понятия с.к. производной и с.к. интеграла позволяют дать корректное описание линейного дифференциального уравнения со случайными возмущениями в правой части и случайными начальными условиями: t  a(t )t  b(t )t , t  t0 ,t 0  v , где t ,  t - с.к. непрерывные случайные функции, t - с.к. производная функции  t , a(t ), b(t ) - непрерывные неслучайные функции, а v - некоторая случайная величина с известным распределением. Случайная функция  t является решением дифференциального уравнения при t  t0 , если будет t t   t0 t0 удовлетворять уравнению: t  v  a( ) d  b( ) d , где все интегралы понимаются в с.к. смысле. Для явного построения решения, введем вспомогательную неслучайную функцию  (t , t0 ) ,  (t , t0 )  a(t ) (t , t0 ), t  t0 . Тогда случайная  (t0 , t0 )  1  удовлетворяющую следующим соотношениям:  t  функция t   (t , t0 )v   (t ,  )b( ) d удовлетворяет исходному дифференциальному уравнению t0 со случайной правой частью. Рассмотрим теперь векторный случай, когда система состоит из n линейных дифференциальных уравнений, записанных в матричной форме: t  A(t )t  B(t )t , t  t0 ,t 0  v , где t - с.к. производная n - мерной случайной функции t  R n ,  t  R m - с.к. непрерывная при t  t0 случайная функция, v  R n - случайный вектор, а матричные функции A(t )  R n n , B(t )  R n m имеют неслучайные кусочно - непрерывные компоненты . t  Общее решение уравнения t  A(t )t  B(t )t ,t 0  v имеет вид: t   (t , t0 )v   (t ,  ) B( ) d , t0 где  (t , t0 ) - матричная переходная функция, удовлетворяющая системе однородных    A(t ) дифференциальных уравнений  , t  t0 . Интеграл в правой части решения понимается  (t 0 , t0 )  E в с.к. смысле и применяется к случайной векторной подынтегральной функции по каждому ее компоненту. Анализ стационарных линейных систем в пространстве состояний при воздействии белого шума. Рассмотрим линейную стационарную систему, описываемую следующим уравнением: x  Ax  Bu; y  Cx; x  R n , u  R q , y  R p , x(t0 )  xt 0 . Пусть входное воздействие u(t )  (u1, u2 ,..., uq )T представляет собой стационарную случайную векторную функцию. В этом случае запись уравнений состояния будем понимать в среднем квадратичном смысле. В соответствии с условием задания случайных функций, такое входное воздействие будет полностью характеризоваться многомерным совместным распределением вида Fut1 ,ut 2 ,..., utk ( z1, z2 ,...., zkq ) . Для задания стационарных в широком смысле случайных функций достаточно рассмотреть первые два момента такой векторной функции и ее компонентов: вектор 13 средних значений mu  M {u (t )}  (mu1 , mu 2 ,..., mu q )T и матрицу авто- и взаимно- корреляционных функций Ru (t )  ( Ru i u j (t )  M {ui0u 0j }; i, j  1,2,..., q , где M {ui0u 0j }  M {(u  mu )(u  mu )T } . Общее решение уравнений состояния можно записать в следующем виде : t xt   (t , t 0 )v   (t , ) B( )u ( )d , где  (t , t0 ) - матричная переходная функция. Интеграл в правой t0 части решения понимается в с.к. смысле и применяется к случайной векторной подынтегральной функции. Для постоянной матрицы A переходная функция имеет вид  (t , t0 )  e A(t  t 0 ) . Если перейти к центрированным значениям вектора состояния x 0  x(t )  mx (t ) и вектора входного воздействия u 0  n(t )  mu (t ) , то уравнение состояния системы можно записать в следующем виде: x 0  Ax 0  Bu 0 . Тогда, предполагая, что вектор состояния x(t ) также является стационарным процессом, можно записать следующее соотношение: M {x 0 ( x 0 )T }  D x (t )  M {Ax 0 ( x 0 )T  Bu 0 ( x 0 )T }  M {x 0 ( x 0 )T AT  x 0 (u 0 )T BT }   ADx (t )  Dx (t ) AT  BM {u 0 ( x 0 )T }  M {x 0 (u 0 )T }BT . Учитывая, что входной сигнал и начальные условия обычно не коррелированны, можно записать следующее соотношение: M {u ( x ) }  M {u 0 T ( xt00 )T } T (t , t0 )  t  M {u t ( )(u ( )) }B G (t   )d   Ru ( ) BT G T (t   )d . T T T t0 t0 Здесь G(t   )   (t , t0 ) 1 ( , t0 ) / 4 /. Рассмотрим случай, когда входной сигнал представляет собой белый шум с независимыми компонентами, то есть u 0  N 0 (t ) и RN (t )  V (t ) (t ) , где V (t ) диагональная матрица интенсивностей векторного единичного белого шума N 0 (t ) . Тогда, с учетом данного допущения, t  получим: M {u 0 ( x 0 )T }  V ( ) ( ) BT GT (t   )d  t0 1 1 1 [V ( ) BT ]  t  V (t ) BT . Множитель 2 2 2 учитывает тот факт, что момент приложения  - функции совпадает с верхним пределом интегрирования. Аналогично можно показать, что M {x 0 (u 0 )T }  1 BV (t ) . 2 Тогда уравнение для соотношения M {x 0 ( x 0 )T } может быть записано в виде: D x (t )  ADx (t )  Dx (t ) AT  BV (t ) BT . Данное матричное уравнение называется дисперсионным уравнением и является уравнением Риккати. Так как матрица Dx (t ) является симметричной матрицей, то оно эквивалентно системе n(n  1) скалярных уравнений. 2 Таким образом, учет того, что белый шум является случайным процессом с независимыми приращениями, то есть не дифференцируемым случайным процессом, привел к тому, что в уравнении для корреляционной функции Dx (t ) появилось дополнительное слагаемое BV (t ) BT . 14 Стохастический интеграл и формула Ито. Как уже обсуждалось выше (см.лекция 14), довольно часто в теории случайных процессов b встречаются стохастические интегралы  f (t , t )dt , где  t - случайная функция. Такие интегралы a вводятся для недифференцируемых процессов  t ( процессов с некоррелированными или ортогональными приращениями), так как в противном случае, когда процесс дифференцируем, b достаточно рассмотреть интеграл  f (t , t )t dt . a Пусть wt - стандартный винеровский процесс с некоррелированными (ортогональными) приращениями с нулевым математическим ожиданием, определен на отрезке числовой прямой t  (a, b) . Можно показать, что для такого процесса существует неубывающая функция F (t ) , такая, что приращение t   s , s  t; t , s  T имеет дисперсию, равную F (t )  F (s) . Процесс с некоррелированными приращениями будет непрерывен в среднем квадратичном (соответственно справа или слева) тогда и только тогда, когда непрерывна (справа или слева ) функция F (t )  M {| wt  wt 0 |2 } . Определим ступенчатые функции, имеющие вид: f (t )  f 0 при t  t1 , f (t )  f1 при t1  t  t 2 , …, f (t )  f n при t  t n . Положим J ( f n )  n  f k  (t k 1 k 1 M {J ( f n )}  0 и M {| J ( f n ) |2 }    t k ) . Легко заметить, что будут справедливы соотношения: n  | f k |2 [ F (tk 1 )  F (tk )] . Тогда интеграл k 0 b  f (t )dwt определяется с a помощью следующего предельного перехода J ( f )  l.i.m J ( f n ) . Такое равенство будет n  выполняться почти наверное, а интеграл J ( f ) называется стохастическим интегралом Ито. При b  этом будут справедливы следующие равенства: M {| J ( f ) |2 }  | f (t ) |2 dF (t ) и a ____ b ___ M {J ( f ) J ( g )}   f (t ) g (t ) dF (t ) . a Без ограничения общности можно считать, что M { f (t )}  0 на интервале t  (a, b) , так как в ~ ~ противном случае J ( f )  J (m f )  J ( f ) , где m f (t )  M { f (t )} , а f (t )  f (t )  m f (t ) - центрированная случайная функция. В приведенной схеме значение интеграла будет зависеть от точки отсчета значения функции f k на каждом k - ом интервале. В интеграле Ито точка отсчета на каждом интервале выбирается как 15 t t ) левый конец отрезка, то есть t k . Если выбирать в качестве точки середину отрезка ( k 1 k , то 2 такой интеграл будет называться стохастическим интегралом Стратоновича, который, в общем случае, будет иметь отличающееся значение от интеграла Ито. Следует отметить, что свойства интеграла Ито отличаются от привычных свойств интеграла Римана-Стилтьеса. Например, формула интегрирования по частям в общем случае не имеет места. Стохастический дифференциал. Формула Ито. При проведенном построении интеграла Ито не строится аппарат дифференцирования. Тем не менее, оказывается возможным получить вариант правила дифференцирования сложной функции, называемой формулой Ито / 2 /. На основе полученного построения стохастического интеграла, введем следующее понятие стохастического дифференциала, которое определим с помощью следующего соотношения : t2 t2 t1 t1  t 2   t1   a(t )dt   b(t )dwt , где a(t ), b(t ) - неслучайные функции, измеримые по L2 (a, b) . Тогда стохастическим дифференциалом процесса  t называют следующее выражение: dt  a(t )dt  b(t )dwt . Теорема. Формула Ито. Пусть процесс  t имеет стохастический дифференциал dt  a(t )dt  b(t )dwt , а неслучайная функция u (t , x) определена на интервале t  (a, b) , x  R , является непрерывной и имеет  (t , x) . Тогда процесс t  u(t , t ) также имеет непрерывные производные u't (t , x), ux (t , x), uxx стохастический дифференциал и справедлива следующая формула (формула Ито): 1  (t ,  t )b 2 (t )]dt  u (t ,  t )b(t )dwt . dt  [ut (t ,  t )  u x (t ,  t )a(t )  u xx 2 Следует отметить отличие правил дифференцирования Ито от правила обычного дифференцирования. Пусть  (t ) является дифференцируемой неслучайной функцией, а  (t )  u(t , (t )) . Тогда можно записать следующее соотношение: d (t )  ut (t ,  (t ))dt  u (t ,  (t ))d (t ) . Очевидно, что оба правила дифференцирования будут  (t , t )  0 . совпадать в случае, когда uxx Пример. Обычно процесс с независимыми ортогональными приращениями wt (винеровский процесс) связывают с броуновским движением частиц жидкости. Уравнения такого движения (уравнение диффузии) появляется в случае большого трения, когда ускорением частицы  x можно пренебречь, то есть x  t . Здесь  t представляет собой случайные силы, возникающие при тепловом возбуждении молекул. Если предположить, как обычно, что величины t коррелированны по Гауссу M {t }  0; M {tt }   (t  t ) , то получим M {x(t )}  0 и M {x 2 (t )}  t . Это означает, что квадрат расстояния до начала координат линейно зависит от времени, если на частицу действуют случайные силы, в виде белого гауссова шума (в отличии от закона x 2 (t )  t 2 при постоянной силе). Можно также показать, что данный результат будет выполняться и тогда, когда ускорением не пренебрегают. 16  t  Q(t )t , Если процесс  t является белым шумом, то процесс wt , связанный с  t уравнением w где Q(t ) детерминированная векторная функция (матрица), является процессом с ортогональными приращениями. При этом, если wt является стандартным винеровским процессом, то такой белый шум называется нормальным (гауссовым) белым шумом. Стохастические дифференциальные уравнения. Рассмотрим стохастическую систему, которая описывается с помощью уравнения dxt  a( xt , t )dt   ( x, t )dwt на интервале t  [0, T ] с начальным условием xt 0  v . Здесь wt  R m стандартный винеровский процесс. Данное уравнение может быть представлено в скалярном виде: dxk , t (t )  ak ( xt , t )dt  обозначения | a( x, t ) | 1 m  ki ( x(t ), t )dwi,t , xk ,t 0  vk , k  1,2,..., n . Введем 2 i 1 n n m k 1 k 1i 1  ak2 ( x, t ) и ||  ( x, t ) || 2   ki2 ( x, t ) . Рассмотрим условия существования и единственности решения приведенной стохастической системы. Теорема / 3 /. Пусть случайная величина v не зависит от {wt , t  [0, T ]}, M {| v |2 }   , а коэффициенты уравнения dxt  a( xt , t )dt   ( x, t )dwt непрерывны по переменным t  [0, T ], x  R n . Пусть также: - найдется такая величина K   , что при всех t  [0, T ], x  R n будет выполняться соотношение | a( x, t ) |2  ||  ( x, t ) || 2  K (1 | x |2 ) ; - найдется такая величина C   , что при всех t  [0, T ], x, y  R n будет выполняться соотношение | (a( x, t )  a( y, t ) | 2  ||  ( x, t )   ( y, t ) || 2  C || x  y || 2 . Тогда на интервале t  [0, T ] существует и единственно ( почти наверное) непрерывное решение xt уравнения dxt  a( xt , t )dt   ( x, t )dwt , xt 0  v , причем M {| xt |2 }  L(1  M {| v |2 }) , где константа L зависит лишь от T , K . Рассмотрим теперь следующее стохастическое дифференциальное уравнение: dt  A(t )t dt  u(t )dt  B(t )dwt ;t 0  v , где A(t )  R nn , B(t )  R n q , u(t )  R n непрерывные неслучайные функции. Можно показать, что процесс, определяемый равенством: t t t0 t0  t   (t , t0 )v   (t , )u ( )d   (t , ) B( )dw будет являться решением стохастического дифференциального уравнения, если функция  (t , t0 ) будет удовлетворять следующей системе:  (t , t0 )  A(t ) (t , t0 ), t  t0 .   (t0 , t0 )  E  Второй интеграл в равенстве, определяющем решение, понимается в смысле Ито. 17 Уравнения моментов для линейных стохастических систем. Рассмотрим линейную стохастическую систему следующего вида dxt  A(t )dt  u(t )dt  B(t )dwt , где xt 0  v , wt - случайный процесс с независимыми приращениям, mw  0; Rw (t , t ' )  Q (t  t ' ) ; а детерминированные матричные функции A(t ), u(t ), B(t ) являются непрерывными . Случайная величина v является независимой от процесса wt . Решением этого уравнения, как уже t t   t0 t0 отмечалось, является процесс  t   (t , t 0 )v   (t , )u ( )d   (t , ) B( )dw , если матричная  (t , t0 )  A(t ) (t , t0 ), t  t0 .  (t0 , t0 )  E  Пусть процесс xt в момент t  [0, T ] имеет среднее значение mx (t )  M {x(t )} и дисперсионную функция  (t , t0 ) , удовлетворяет системе уравнений:  матрицу Dx (t )  M {( x(t )  mx (t ))( x(t )  mx (t ))T } , а начальное значение xt 0  v не зависит от wt . Тогда функции mx (t ), Dx (t ) будут являться решениями следующих систем обыкновенных  x (t )  A(t )mx (t )  u(t ) , mx (0)  mv ; дифференциальных уравнений: m D x (t )  A(t ) Dx (t )  Dx (t ) AT (t )  B(t )QBT (t ) , Dx (0)  Dv , где Dv  M {(v  M {v})(v  M {v})T } дисперсионная матрица вектора x(0)  v . Данные уравнения называются уравнениями метода моментов. Пример. Пусть скалярная случайная функция xt удовлетворяет уравнению dxt  adt  bdwt , xt  0  v, a  0 . t Решением этого уравнения будет процесс xt  e at mv  be at  e  a w( )d . Уравнения моментов будут иметь вид:  x (t )  amx (t ) , mx (0)  mv и D x (t )  2aDx (t )  b 2 , Dx (0)   v2 . Отсюда получим следующие m соотношения mx (t )  mv e at , Dx (t )  e 2at v2  mx (t )  0, Dx (t )  b2 (1  e 2at ) . Следовательно, при t   получим 2|a| b2 2 для любых mV , Dv   v . 2|a| Определение. Система дифференциальных уравнений dxt  Adt  Du(t )dt  Bdwt , xt 0  v с постоянными матрицами A, D, B называется асимптотически устойчивой, если все корни {k } уравнения det(E  A)  0 лежат в левой полуплоскости, то есть Re(k )  0, k  1,2,..., n Утверждение. Если система dxt  Adt  Dudt  Bdwt , xt 0  v асимптотически устойчива и u  u(t )  const , то существуют mx (t )  mx и Dx (t )  Dx при t   . Предельные значения mx , Dx удовлетворяют стационарным уравнениям метода моментов Am x  Du  0 и ADx (t )  Dx (t ) AT  BBT  0 18 Рассмотрим теперь линейное стохастическое уравнение dxt  Ax (t )dt  Bdwt , xt 0  v , где mw  0; Rw (t , t ' )   (t  t ' ) . Тогда величина Dx (t ) удовлетворяет уравнению моментов D x (t )  ADx (t )  Dx (t ) AT  BBT , Dx (0)  Dv . Можно показать, что Dx (t )  0 при любом Dv  0 тогда и только тогда, когда rank{B, AB,..., An1 B}  n , то есть линейная стохастическая система является вполне управляемой. Различные формы представления стохастических систем. Во многих, практически важных приложениях, стохастические системы описываются уравнением вида xt  a( xt , t )   ( xt , t )t , xt 0  v . Здесь xt - производная в среднеквадратическом (с.к.) смысле случайного процесса xt  R n ;  t  R m - с.к. – непрерывная при t  0 случайная функция; a( x, t ), ( x, t ) - непрерывные неслучайные функции; v - случайная величина. Данная запись является, по сути, записью следующего соотношения: t t xt  v   a( x( ), )d    ( x( ), ) d , где все интегралы понимаются в с.к. – смысле. Если полоса пропускания системы значительно меньше полосы частот равномерности спектральной плотности стационарного возмущающего процесса  t , то этот процесс может считаться белым шумом по отношению к данной системе, то есть имеющим постоянную спектральную плотность в диапазоне частот существенном для динамики системы. Таким образом, понятие белошумности случайного процесса достаточно условно. Один и тот же случайный процесс может считаться белошумным по отношению к одной динамической системе и «цветным» по отношению к другой, имеющей более широкую полосу пропускания. Для формализации этого утверждения, была предложена параметризованная модель белошумного случайного процесса. Если ввести малый параметр   0 и построить процесс  t  1   , где t  2 t  R m - стационарный гауссов случайный процесс с корреляционной функцией R ( ) , такой, что   R (t , t   )d  E , и спектральной плотностью S ( ) , то процесс  t будет иметь корреляционную  функцию R ( ) и спектральную плотность S  ( ) , такие, что при   0 будут выполняться соотношения R ( )  1  2 R (  )   ( ) E , S ()  S ( 2)  S (0) . 2  19 Таким образом, стохастическая система, для параметризованных случайных возмущающих процессов, принимает вид xt  a( x, t )   ( x, t ) процесс  t  1     t , x(t0 )  v . При достаточно малом значении  2 может считаться белошумным по отношению к данной динамической системе, t  1 2 так как полоса равномерности спектра возмущающего процесса 1   t неограниченно расширяется 2 при   0 , а полоса пропускания системы остается постоянной. Если полоса пропускания системы не достаточно узкая по отношению к полосе равномерности спектра возмущающего процесса  (t ) , то белошумная аппроксимация этого процесса будет не адекватна постановке задачи стохастической динамики. В этих случаях, необходим учет свойств «цветности» процесса  t . Для этого спектральную плотность процесса подвергают процедуре факторизации, то есть представляют в следующем виде: S ( )  (i )(i )  B(i ) B(i ) , где A(i ) A(i ) B(s), A(s) - многочлены конечных степеней порядка nB , n A ; nB  n A . Причем многочлен A(s) является гурвицевым. Отсюда можно построить уравнения соответствующего формирующего фильтра: zt  Dzt  Gnt ; t  C T zt , где nt - процесс типа белого шума; ( s)  C T ( sE  D) 1 G  B( s ) . A( s) Таким образом, достаточно универсальной математической моделью, описывающей стохастическую динамику разнообразных систем, будет следующая: xt  a( xt , t )   ( xt , t )t , x(t0 )  v . Однако, необходимо отметить, что случайные процессы типа белого шума, имеют бесконечную дисперсию. То есть, говорить о с.к. – смысле полученного соотношения надо с большой осторожностью, так как все математические преобразования придется выполнять с обобщенными функциями. Поэтому такую систему обычно приводят к форме стохастического уравнения Ито при проведении аналитических вычислений, то есть к форме dxt  a( xt , t )   ( xt , t )dwt , x(t0 )  v . Решение такого уравнения представляется в виде соотношения t t xt (t )  v   a( x ,  )d    ( x ,  )dw , где первый интеграл в правой части понимается в с.к. – смысле, а второй – является интегралом Ито. Приложение. Стохастические интегралы Ито и Стратоновича Рассмотрим предварительно скалярный стандартный винеровский процесс wt с нулевым математическим ожиданием и коэффициентом диффузии (дисперсией)  2  1 . Пусть интеграл t I (t , s)   w dw ,0  s  t   от винеровского процесса по этому же винеровскому процессу. Если бы s функция wt была бы непрерывно дифференцируема на интервале [ s, t ] , то тогда бы имело место равенство I (t , s)  0.5( wt2  ws2 ) . Определим численную аппроксимацию интеграла I (t , s) N  [wt t  0 следующим образом: I (t , s)  lim k 1 k 1  (1   ) wt k ](wt k 1  wt k ) , где t  ts , N 20 tk  s  (k  1)t; k  1,2,..., N  1 и  [0,1] , параметр определяющий точку взятия значения wt на интервале. Тогда, проведя преобразования, можно записать: I (t , s)  lim [ t   N 1 N 2 1 N 1 1 ( wt k 1  wt2k )  (  )  ( wt k 1  wt k ) 2 ]  ( wt2  ws2 )  (  ) lim  ( wt k 1  wt k ) 2 .  2 k 1 2 k 1 2 2 t  0 k 1 Величины ( wt k 1  wt k ) , k  1,2,..., N  1 представляют независимые случайные величины, распределенные по нормальному закону с нулевым математическим ожиданием и дисперсией, равной t (см.свойства винеровского процесса в лекции 14). Тогда M {(wt k 1  wt k ) 2 }  t и D{(wt k 1  wt k ) 2 }  M {[(wt k 1  wt k ) 2  t ]2 }  (t ) 2 M {[( величина ( M {( wt k 1  wt k wt k 1  wt k t t wt k 1  wt k t ) 2  1]2 } . Так как случайная ) 2 распределена по закону  2 с одной степенью свободы, то получим: ) 2 }  1 и D{( wt k 1  wt k t ) 2 }  2 . Отсюда найдем, что D{(wt k 1  wt k ) 2 }  2(t ) 2 . Введем случайные величины  k  ( wt k 1  wt k ) 2  t; k  1,2,..., N . Данные величины являются, очевидно, как функции независимых переменных и имеют следующие первый и второй моменты: M { k }  0; D{ k }  2(t ) 2  N lim  (wt k 1  wt k ) 2  lim t  0 k 1 2(t  s) 2 N2 . Тогда найдем, что N  (t  k )  t  s  lim t  0 k 1 N N k  t  s , так как M { k }  0 и t  0 k 1 k 1 N D{ k }  N 2(t ) 2  0 , а любая случайная величина с нулевым математическим ожиданием и t  0 k 1 нулевой дисперсией есть (почти наверное) нуль. Таким образом, можно записать следующее 1 2 1 ( wt  ws2 )  (  )(t  s) . То есть получили однопараметрической 2 2 семейство интегралов с параметром  [0,1] . Обычно на практике используют два типа: соотношение: I (t , s)  I  (t , s)  1 1 1 1 I 0 (t , s)  ( wt2  ws2 )  (t  s);  0 и I 0.5 (t , s)  ( wt2  ws2 );  . Первый интеграл I 0 (t , s) 2 2 2 2 называют стохастическим интегралом Ито, а второй – стохастическим интегралом Стратоновича. Несмотря на то, что интеграл I 0.5 (t , s) совпадает с интегралом от неслучайной функции w  w(t ) , интеграл Ито более удобен при анализе стохастических моделей. Определение. Пусть wt , t T  [0, ) - n - мерный винеровский процесс, выходящий из точки 0, и f ( wt , t ) матричная функция размера n  n , определенная для всех t  T . Если для любых s, t  T , таких, что t  s , и для любого  [0,1] существует предел N  f [(wt N  lim k 1 k 1  (1   ) wt k ), t k 1  (1   )t k ](wt k 1  wt k ) , где t  ts , N tk  s  (k  1)t; k  1,2,..., N  1 , то его называют стохастическим интегралом функции f (wt , t ) по 21 винеровскому процессу wt и обозначают I  (t , s) . Если   0 , то предел называют стохастическим интегралом Ито, а при значении   0.5 - стохастическим интегралом Стратоновича. Чтобы не путать интегралы Ито и Стратоновича вводят следующие обозначения: t I 0 (t , s)   f ( wt , t )dwt для интеграла Ито, а для интеграла Стратоновича используют s t дополнительный символ ‘*’ - I 0.5 (t , s)   f (wt , t )dwt . s Теорема. Интеграл Ито I 0 (t , s) имеет нулевое математическое ожидание и дисперсию t D{I 0 (t , s)}   M {[ f ( wt , t )]2 }dt . s Доказательство. По определению интеграла можно записать  lim N  M { f (wt N   k 1 k N  f (wt N  M {I 0 (t , s)}  M { lim k 1 k , t k )(wt k 1  wt k )}  , t k )}M {(wt k 1  wt k )}  0 , так как математическое ожидание винеровского процесса равно нулю и случайные величины wt k и ( wt k 1  wt k ) независимы. t Равенство D{I 0 (t , s)}   M {[ f ( wt , t )]2 }dt доказывается аналогично, с учетом соотношения s M {(wt  ws ) 2 }  t  s . Теорема. Интегралы Ито и Стратоновича связаны следующим равенством: t 1 f ( x, t ) I 0.5 (t , s)  I 0 (t , s)   | x  wt dt . 2 s x Численное моделирование линейных стохастических систем. Некорректность применения численных методов для обыкновенных дифференциальных уравнений для стохастических дифференциальных уравнений. Большинство численных методов для обыкновенных дифференциальных уравнений строятся с помощью разложения в ряд Тейлора. Рассмотрим скалярное стохастическое дифференциальное уравнение вида: xt  a( xt , t )  b( xt , t )nt , где xt  R - решение уравнения; a, b : R  [0, T ]  R ; nt - случайный процесс в виде гауссова белого шума, nt dt  dwt . Запишем это же уравнение в форме стохастического дифференциального уравнения Ито: dxt  a( xt , t )  b( xt , t )dwt , где wt - стандартный винеровский процесс. Применим обычную идеологию численных методов к данным стохастическим уравнениям / 6 /. Согласно формуле Тейлора при значении s  t имеем: 22 dxt d 2 xt ( s  t ) 2 xs  xt  (s  t )  2  O((s  t ) 3 ) , где запись  t  O( s  t ) ) понимается как: dt 2 dt lim s t M { t2 } ( s  t )  C  const   . Используя первое уравнение, найдем: d (s  t ) 2 [a( xt , t )  b( xt , t )nt ]  O((s  t ) 3 ) . Правая часть dt 2 данного соотношения неопределенна. Это связано с тем, что белый шум nt не является xs  xt  a( xt , t )(s  t )  b( xt , t )(s  t )nt  дифференцируемым в каком – либо смысле случайным процессом. Рассмотрим теперь применение этой же идеологии в отношении уравнения Ито. Используем для случайных функций a( xt , t ), b( xt , t ) такие же правила дифференцирования, как и для неслучайных функций, то есть: da( xt , t )  a( xt , t ) a( xt , t ) b( xt , t ) b( xt , t ) dxt  dt , db( xt , t )  dxt  dt , как это x t x t делается в обычных алгоритмах численного решения. Запишем решение уравнения в интегральной s s   t t форме: xs  xt  a( x , )d  b( x , )dw , где: a( x ,  ) a( x ,  ) a( x , ) a( xs , s)  a( xt , t )   [a( x ,  )  ]d   b( x ,  ) dw , x  x t t s s b( x , ) b( x , ) b( x , ) b( xs , s)  b( xt , t )   [a( x , )  ]d   b( x , ) dw . x  x t t s s Подставляя данные выражения в разложение Тейлора, получим: s xs  xt  a( xt , s)(s  t )  b( xt , t )  dw  b( xt , t ) t s b( xt t ) dw 1 dw  x t t ss  [a( xt , t ) b( xt , t ) b( xt , t ) dxt  ] d  1 dw  O(( s  t ) 2 ) . x t t t Однако, поскольку для случайных процессов a( xt , t ), b( xt , t ) , правила дифференцирования должны выполняться по формуле Ито, то с вероятностью 1 получим разложение: s s s b( xt t ) a( xt , t ) xs  xt  a( xt , s)( s  t )  b( xt , t )  dw  b( xt , t ) dw 1 dw b( xt , t ) dw 1 d   x t t x  t t t s s b( xt , t ) b( xt , t ) 1 2  2b( xt , t )  [a( xt , t ) dxt   b ( xt , t ) ]  d  1 dw  O((s  t ) 2 ) . Таким образом, 2 x t 2 x t t приходим к выводу, что при применении обычных правил дифференцирования будут отсутствовать часть слагаемых, которые имеют порядок O( s 3  t) 2 ) . То есть формальное применение правил дифференцирования для неслучайных функций к случайным процессам приводит либо к неопределенности, либо к неверному результату. 23 Моделирование линейных стохастических систем. Рассмотрим линейную систему, описываемую стохастическим дифференциальным уравнением вида: dxt  A(t ) xt dt  B(t ) f (t )dt  G(t )dwt ; xt t 0  v , где xt  R n - случайный процесс, являющийся решением уравнения с начальным случайным условием xt  t 0  v ; wt  R m - винеровский случайный процесс с независимыми компонентами wt  (w1,t ,..., wm,t )T , wk ,t  R , p M {wt }  0; M {wt wtT }   2w  diag{ w2 k }m k 1 ; f (t )  R - детерминированное возмущение; A(t ), B(t ), G(t ) - соответствующие матричные функции, удовлетворяющие условиям существования и единственности решения стохастического уравнения. Приведенную систему можно рассматривать, как систему уравнений Ито или Стратоновича, поскольку можно показать, что в силу своей структуры оно имеет одно и то же решение, в каком бы смысле она не рассматривалась. Поэтому, часто применяемая в технической литературе форма описания линейной стохастической системы в виде: xt  A(t ) xt  B(t ) f (t )  G(t ) wnt ; xt t 0  v с белым шумом nt с единичной дисперсионной матрицей в правой части не приводит к неоднозначности математического определения решений системы. Рассмотрим стационарный случай системы, то есть когда ее модель описывается стохастическим уравнением: dxt  Axt dt  Bf (t )dt  Gdwt ; xt  t 0  v . В этом случае решение будет определяться с вероятностью 1 в виде: xt  e A(t  t 0 ) t t t0 t0 v   e A(t  s ) Bf ( s)ds   e A(t  s ) Gdws . Для численного решения выберем регулярную сетку t k 1  t k  h; t0  0; k  0,1,2,... , где h - шаг сетки. Введем следующие h h A( h  s ) ˆ  e Ah , ( Bˆ f )  e A( h  s ) Bf (kh  s)ds , nˆ Gdw( kh  s ) . Тогда непрерывную обозначения A k k 1   e  ˆ x  ( Bˆ f )  nˆ , где систему можно аппроксимировать дискретной моделью вида: xk 1  A k k k 1 случайные величины xk , ( Bˆ f ) k , nˆ k 1 сохраняют свое значение на интервале h . Векторный гауссов дискретный случайный процесс n̂k имеет дисперсионную матрицу: h D{nˆ k 1}  Dnˆ (h)   e A( h  s )G 2wG T e A T (h  s) ds . Несложно показать, что дисперсионная матрица Dnˆ (t ), t [0, h] является решением уравнения: D nˆ (t )  ADnˆ (t )  Dnˆ (t ) AT  G 2wGT ; Dnˆ (0)  0 . Таким образом, для моделирования случайной составляющей решения надо на каждом шаге моделировать гауссов случайный процесс с нулевым средним и дисперсионной матрицей Dnˆ (h) . Если шаг является переменным, то надо вычислять, при каждом его изменении, новую дисперсионную матрицу, что является вычислительно трудоемким. Вычисление систематической составляющей решения линейной стохастической модели. Оценка систематической составляющей сводится к вычислению интеграла h Z h   e A( h  s ) Bf (kh  s)ds . Представим векторную функцию f (kh  s), s [0, h] с точностью  24 следующей полиномиальной функцией: f (kh  s)  L  q j (k )s j  r (k , s) , где | r (k , s) | s L 1 для j 0 величин k  0,1,2,... ; а векторные коэффициенты q j (k ), j  1,2,..., L могут быть различными способами. Тогда можно записать Z h  L  Z h, j Bq j (k )  Rh, L 1 , где j 0 h h || Z h, L 1 ||||  e A( h  s ) Br (k , s)ds ||  || Z h, L 1 ||  || B || , Z h, j   s j e A( h  s ) ds . Для вычисления интегралов Z h, j можно использовать следующие очевидные рекуррентные соотношения: Z h,0  A1 (e Ah  E ), Z h, j   A1 (h j E  jZ h, j 1 ); j  1,2,..., L . Отсюда найдем: j h l Al ); j  1 или, отбрасывая остаточный член Rh, L 1 , получим l ! l 1 Z h, j   A ( j 1) j!( AZ h,0   следующее приближенное представление: Z h, j   A ( j 1) L  j 0 j hl Al ). l ! l 0 j!(e Ah   При кусочно – постоянной аппроксимации ( L  0) выполняется соотношение q0 (k )  f (kh)  f k . Тогда Z h  A1 (e Ah  E ) Bf k  Bˆ f k . То есть ( Bˆ f ) k  Bˆ f k , где Bˆ  A1 (e Ah  E ) B . При линейной интерполяции ( L  1 ) получим: q0 (k )  f k , q1 (k )  Z h  A1 (e Ah  E ) Bf k  1 ( f k 1  f k ) . Отсюда найдем h A2 Ah (e  E  hA) B( f k 1  f k ) h Вычисление случайной составляющей решения линейной стохастической модели. h Рассмотрим дисперсионную матрицу Dnˆ (h)  M {nˆ k 1nˆ kT1}   e A( h  s ) G 2wG T e A T (h  s) ds . Так как эта матрица является квадратной, симметричной и невырожденной, то ее SVD разложение (см.приложение, лекции 4 и 5) можно записать в виде: Dnˆ (h)  M {nˆk 1nˆkT1}  T2DT 1 , где матрица  D  diag (1, 2 ,.., n ) , k , k  1,2,..., n - собственные значения матрицы Dnˆ (h) , T - матрица ортонормированных собственных векторов матрицы Dnˆ (h) . Очевидно, что, так как матрица T является ортогональной. Представим вектор nˆ k 1 в виде nˆ k 1  T D nk 1 , где nk 1 - случайный вектор – столбец, такой, что его компоненты n j , k 1 , j  1,2,.., n являются случайными величинами с нулевым средним с единичной дисперсией. Поэтому, чтобы смоделировать случайный вектор nˆ k 1 надо найти матрицу Dnˆ (h) , ее собственные значения и собственные векторы, а также смоделировать n независимых гауссовых случайных величин с нулевым средним и единичной дисперсией, для построения вектора nˆ k 1 . Для нахождения матрицы Dnˆ (h) можно использовать дифференциальное уравнение D nˆ (t )  ADnˆ (t )  Dnˆ (t ) AT  G 2wGT ; Dnˆ (0)  0 . 25 Приближенный метод моделирования линейных стохастических систем. Вычисление дисперсионной матрицы Dnˆ (h) является алгоритмически достаточно трудоемким. Поэтому достаточно часто используют приближенный метод численного решения стохастических уравнений. Для этого белый шум nt , (dwt   wnt dt ) заменяется на промежутке [kh, (k  1)h]; k  0,1,2,... кусочно постоянным случайным процессом nt с достаточно малым шагом дискретизации  . Запишем следующее приближенное представление стохастической системы в следующей форме: ~ xt  A~ xt  Bu (t )  Gnt ; yt  C~ xt ; xt t 0  v , где nt - кусочно постоянный на h ; N  1 , случайный процесс с независимыми значениями N на этих промежутках. То есть имеют место соотношения: nt  zk , r  const при значении промежутках [kh, (k  1)h]; k  0,1,2,...;  t  [r , (r  1) ) ; M {zk , r }  0 ; M {zk , r zkT, r }  1   w . Корреляционная функция случайного процесса |t | 1  (1  ) 2w , | t |  . Таким образом, при   0 корреляционная функция nt имеет вид: Rn (t )     , | t |    процесса nt будет стремиться к  - функции Дирака, а спектральная плотность к постоянной величине. Если выполняется условие max   2 , где max - максимальная частота, наблюдаемая 8 в системе, то будет справедливо соотношение 0.95 2w  S n ( )   2w . Таким образом, на интервале длиной h должно быть не менее 8 подинтервалов  постоянства процесса nt . Построим дискретное представление линейной стохастической модели на регулярной сетке h h xk 1  e Ah ~ xk   e A( h  s ) Bu (kh  s)ds   e A( h  s )Gdn( kh  s ) . Оценим [kh, (k  1)h]; k  0,1,2,...; в виде: ~ стохастический интеграл, используя кусочно постоянное приближение процесса nt на более мелкой сетке с шагом  : h e A( h  s ) N Gdn( kh  s )   j  j 1( j 1) e A( h  s ) dsGn( kN  j 1)  1  N A1 (e A  E )  e A( h  j )G  z k , j . Представляя j 1 дисперсионную матрицу  w  T wT 1 в виде SVD представления, получим: h A( h  s ) Gdn( kh  s )  e 1  N ~ ~ z k , j - вектор-столбец A1 (e A  E )  e A( h  j )G~ z k , j , где G   GT w , ~ j 1 независимых стандартных гауссовых случайных величин. Представляя систематическую составляющую с помощью кусочно – постоянной аппроксимации, h найдем, что e A( h  s ) Bu (kh  s)ds  A1 (e Ah  E ) Bu k . Тогда можно записать следующую рекуррентную формулу для приближенного моделирования линейной стохастической системы: N 1 1 A ~ ~ xk 1  e Ah ~ xk  A1 (e Ah  E ) Bu k  A (e  E )  e A( h  j )G~ z k , j ; yk  C~ xk , где N  8 .  j 1 26 Алгоритм моделирования можно записать в следующей форме. 1. Определение временного интервала моделирования [0, t max] и задание шагов размеров сеток h  t max h , k  1,2,..., M ;   , j  1,2,..., N ; N  8 M N 2. Вычисление с заданной точностью матриц: e A , e Ah , A1 (e Ah  E ) B, 1  ~ A1 (e A  E ); G 3. k  1; ~ x1  xt 0  v ; u1  u(t  0) . 4. Численное моделирование N случайных независимых гауссовых векторов ~ z k , j ; j  1,2,..., N размерности m . ~ 5. Вычисление с заданной точностью матриц e A( h  j )G, j  1,2,..., N . 6. Вычисление ~ xk 1  e Ah ~ xk  A1 (e Ah  E ) Bu k  1  N ~ A1 (e A  E )  e A( h  j )G~ z k , j ; yk  C~ xk j 1 7. Если k  M , то k  k  1 и переход к шагу 4. 8. Конец работы алгоритма Расчет моментных характеристик решений линейных стохастических уравнений. Пусть задана линейная стохастическая система dxt  Axt dt  Bf (t )dt  Gdwt ; yt  Cxt ; xt t 0  v , где xt  R n - случайный процесс, являющийся решением уравнения с начальным случайным условием xt  t 0  v ; yt  R q - выход системы; wt  R m - винеровский случайный процесс с независимыми p компонентами wt  (w1,t ,..., wm,t )T , wk ,t  R , M {wt }  0; M {wt wtT }   2w  diag{ w2 k }m k 1 ; f (t )  R - детерминированное возмущение; A, B, G, C - постоянные вещественные матрицы соответствующего размера. Введем обозначения mx (t )  M {xt }; xt0  xt  mx (t ) , Rx (t , s)  M {xt0 xs0 }; Dx (t )  Rx (t , t ) . Используя уравнения моментов можно записать следующие соотношения:  x (t )  Am x (t )  Bf (t ); mx (0)  M {v} ; m D x (t )  ADx (t )  Dx (t ) AT  G 2wGT ; Dx (0)  M {(v  M {v})(v  M {v})T }  Dv ;  Rx (t , s)  AR x (t , s); Rx ( s, s)  Dx ( s), t  s и Rx (t , s)  RT (s, t ), t  s . t Если f (t )  f  const и матрица A является гурвицевой, то приведенные соотношения имеют стационарное решение, определяемое уравнениями: Am x ()  Bf (t )  0 ; ADx ()  Dx () AT  G 2wGT  0 ; A(t  s )  Dx (), t  s e Rx (t , s)  Rx (t  s)   AT ( s  t )  D (  ) e ,t  s  x Таким образом, процесс xt является асимптотически стационарным. Моментные характеристики выхода yt системы определяются соотношениями: M { yt }  m y (t )  Cmx (t ) ; M { yt0 ( yt0 )T }  D y (t )  CD y (t )C T ; R y (t , s)  CR x (t , s)C T . В установившемся режиме, если матрица A является гурвицевой, соотношения принимают вид: m y (t )  Cmx () , 27 D y (t )  CD y ()C T , R y (t  s)  CR x (t  s)C T . Спектральная плотность стационарного выходного процесса yt определяется формулой 1 S y ( )  2   R y (t )e  it dt  C (iE  A) 1 G 2wGT (iE  AT ) 1 C T , где C (sE  A) 1G w - матричная  передаточная функция системы от входа nt , (dwt   wnt dt ) в виде белого шума с единичной дисперсионной матрицей. Функции среды Matlab для моделирования случайных величин. Для формирования массива элементов, распределенных по равномерному закону, используется функция rand (.) . Синтаксис функции rand (n) формирует массив размера n  n , элементами которого являются случайные величины, распределенные по равномерному закону в интервале (0, 1). Синтаксис rand (m, n) формирует массив размера m  n , элементами которого являются случайные величины, распределенные по равномерному закону в интервале (0, 1). Команда X  rand (size( A)) формирует массив соразмерный с матрицей A , элементами которого являются случайные величины, распределенные по равномерному закону в интервале (0, 1). Синтаксис функции rand без аргументов формирует одно случайное число, подчиняющееся равномерному закону распределения в интервале (0, 1), которое изменяется при каждом последующем вызове. Алгоритм генерации равномерно распределенных случайных чисел основан на линейном конгруентном методе, при котором вычисление следующего случайного числа реализовано согласно соотношению: seed  (77  seed )(mod( 231  1) . Синтаксис rand (' seed ' ) возвращает текущее значение базы (начального значения) генератора случайных чисел, а команда rand (' seed ' , x0) присваивает базе (начальному значению) генератора случайных чисел значение x0 . Этот результат может оказаться иным и зависит от версии системы и предыстории сеанса работы. Для формирования массива элементов, распределенных по нормальному закону, используется функция randn(.) . Синтаксис функции randn(n) формирует массив размера n  n , элементами которого являются случайные величины, распределенные по нормальному закону с математическим ожиданием 0 и среднеквадратическим отклонением 1. Синтаксис randn(m, n) формирует массив размера m  n , элементами которого являются случайные величины, распределенные по нормальному закону с математическим ожиданием 0 и среднеквадратическим отклонением 1. Команда X  randn(size( A)) формирует массив соразмерный с матрицей A , элементами которого являются случайные величины, распределенные по нормальному закону с математическим ожиданием 0 и среднеквадратическим отклонением 1. Синтаксис randn без аргументов формирует одно случайное число, распределенное по нормальному закону с математическим ожиданием 0 и среднеквадратическим отклонением 1, которое изменяется при каждом последующем вызове. Этот результат может оказаться иным и зависит от версии системы и предыстории сеанса работы. 28 Список функций пакета Statistics Toolbox Matlab для генерации псевдослучайных чисел. Генерация псевдослучайных чисел по заданному закону распределения  betarnd - Бета распределение  binornd - Биномиальное распределение  chi2rnd - Функция распределения хи-квадрат  exprnd - Экспоненциальное распределение  frnd - Распределение Фишера  gamrnd - Гамма распределение  geornd - Геометрическое распределение  hygernd - Гипергеометрическое распределение  iwishrnd - Обратная матрица случайных чисел распределения Уишарта  lognrnd - Логнормальное распределение  mvnrnd - Многомерное нормальное распределение  mvtrnd - Многомерное распределение Стьюдента  nbinrnd - Отрицательное биномиальное распределение  ncfrnd - Смещенное распределение Фишера  nctrnd - Смещенное распределение Стьюдента  ncx2rnd - Cмещенное хи-квадрат распределение  normrnd - Нормальное распределение  poissrnd - Распределение Пуассона  random - Параметризованная функция генерации псевдослучайных чисел  raylrnd - Распределение Релея  trnd - Распределение Стьюдента  unidrnd - Дискретное равномерное распределение  unifrnd - Непрерывное равномерное распределение  weibrnd - Распределение Вейбулла  wishrnd - Матрица случайных чисел распределения Уишарта Функция normrnd(.) генерации псевдослучайных чисел по нормальному закону Синтаксис normrnd(MU , SIGMA) предназначен для генерации псевдослучайного числа по нормальному закону для каждой пары параметров MU (математического ожидания) и SIGMA (среднего квадратического отклонения). Размерность векторов или матриц параметров MU и SIGMA должна быть одинаковой. Скалярный параметр увеличивается до размера остальных входных аргументов. Команда R  normrnd(MU , SIGMA, m) позволяет получить вектор псевдослучайных чисел из m элементов, распределенных по нормальному закону для параметров MU и SIGMA , где m - вектор размерностью 1x2 определяющий размерность матрицы R . Синтаксис normrnd(MU , SIGMA, m, n) позволяет получить матрицу псевдослучайных чисел с размерностью m  n элементов, распределенных по нормальному закону для параметров MU , SIGMA . Функция mvnrnd(.) генерации псевдослучайных чисел по многомерному нормальному распределению. Команда R  mvnrnd(MU , SIGMA) генерирует матрицу n  d псевдослучайных чисел распределенных по многомерному нормальному закону с параметрами математического ожидания MU и матрицей корреляции SIGMA . Размерность матрицы MU равна n  d . Функция mvnrnd(.) генерирует каждый ряд R на основе соответствующего ряда значений MU . Матрица SIGMA 29 должна быть квадратной и положительно определенной. SIGMA может быть задана матрицей с размерностью d  d или 3-х мерным массивом с размерностью d  d  N . Если матрица SIGMA задана как 3-х мерный массив, то каждый ряд R генерируется с использованием страницы массива SIGMA , то есть R(i, :) генерируется с использованием значений MU (i,:) и SIGMA(:, :, i) . Если массив MU задан как вектор с размерностью 1  d , то функция mvnrnd(.) формирует матрицу размерности соответствующей размерности массива SIGMA . Команда R  mvnrnd(MU , SIGMA, cases ) генерирует матрицу псевдослучайных чисел с размерностью cases  d , распределенных по многомерному нормальному распределению для вектора средних MU с размерностью 1  d и матрицы SIGMA с размерностью d  d . Список использованных источников к лекции 15. 1. Свешников А.А. Прикладные методы теории случайных функций. Москва, Наука, 1968 – 463 с. 2. Оксендаль Б. Стохастические дифференциальные уравнения. Москва, Мир, 2003 – 408 с. 3. Миллер Б.М., Панков А.Р. Теория случайных процессов в примерах и задачах. Москва, Физматлит, 2007 – 320 с. 4. Методы классической и современной теории автоматического управления. Том.2: Статистическая динамика и идентификация систем автоматического управления / Под ред. Пупкова К.А., Егупова Н.Д. Москва, МГТУ им. Н.Э. Баумана, 2004 – 640 с. 5. Волков И.К., Зуев С.М., Цветкова Г.М. Случайные процессы. Москва, МГТУ им. Н.Э. Баумана, 1999 – 446 с. 6. Кузнецов Д.Ф. Стохастические дифференциальные уравнения: теория и практика численного решения. СПб. Издательство Политехнического университета, 2010 – 816 с.