Модифицированный алгоритм статистического моделирования систем со случайным периодом квантования Текст научной статьи по специальности «Математика»

УДК 519.676

Т.А. Аверина

МОДИФИЦИРОВАННЫЙ АЛГОРИТМ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ

СИСТЕМ СО СЛУЧАЙНЫМ ПЕРИОДОМ КВАНТОВАНИЯ

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

Численные методы, стохастические дифференциальные уравнения, системы со случайной структурой.

T.A. Averina

MODIFIED ALGORITHM OF STATISTICAL MODELING OF SYSTEMS WITH A RANDOM QUANTIZATION PERIOD

The algorithm for statistical modeling of systems with a random quantization period, which can be described as a system with random structure with distributed transitions has been constructed. The offered algorithm is based on numerical methods of the solution to the stochastic differential equations and uses the modified maximum cross-section method when the intensity of transition depends on a vector of state.

Numerical methods, stochastic differential equations, systems with random structure.

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

большое место занимают локационные системы, предназначенные для измерения параметров, определяющих пространственное положение объекта, таких, как угловые координаты, дальность, углы ориентации. Если в течении времени эти параметры изменяются и система измеряет их с заданной степенью точности, она называется следящей. Радиоизотопные локационные следящие системы используются в космической технике в задачах стыковки сближающихся аппаратов, в робототехнике, при управлении атомными энергетическими установками и в ряде других случаев. Системы со случайным периодом квантования можно описать как системы со случайной структурой с распределенными переходами [1]. Характерными особенностями таких систем являются: структурная неопределенность (смена структуры в случайные моменты времени в процессе функционирования) и стохастичность процессов в них.

Система со случайной структурой характеризуется вектором состояния Y(t) и номером структуры Ц^=1,.„, S0; S0 - число детерминированных структур. Номер L(t) является

случайным дискретным скалярным процессом, принимающим целочисленные значения. Достаточно общая математическая модель динамической непрерывной нелинейной стохастической системы со случайной структурой записывается как задача Коши для стохастических дифференциальных уравнений (СДУ). Векторное уравнение для фиксированной 1-й структуры имеет вид СДУ в смысле Стратоновича:

dY (0 = a а ^, t)dt + а{1 ^, t)dW ^), Y (^) = Y0, l = 1,..., S0. (1)

Для каждой 1-й структуры вектор состояния системы Y(t) является непрерывным случайным процессом размерности п (1); W(t) - m(1) -мерный стандартный винеровский процесс; a(1) (У^) - п (1)-мерная вектор-функция; а('1 ^, t) - матричная функция размера п (1) хm(1). Начальное состояние системы задается случайным вектором Y0. Случайный дискретный

процесс L(t) может быть произвольным немарковским, марковским или условно марковским, зависящим от вектора Y(t).

Пусть процесс L(t) является условно марковским процессом и зависимость от Y(t) проявляется статистически: моменты перехода из одного состояния в другое случайным образом зависят от изменения фазовых координат. Переходы от одной структуры к другой могут происходить при любых значениях Y(t), но с различной вероятностью. Системы, обладающие подобными свойствами, называются системами с распределенными переходами или системами с условной марковской структурой при распределенных переходах [1]. Условные вероятности перехода из 1-й структуры к г-й для малых временных интервалов Лt выражаются через условные интенсивности переходов У1г ^, t) и имеют вид:

plг (г, t + Лt I ^ t, Y) = У1г ^, t)Лt + o(Лt), l Ф г; (2)

N0

pH (^ t + Лt 11, t, Y) = 1 - уи ^, t)Лt + о^), уи ^, t) = ^ у1г ^, t),

г=1ф1

где о^) является малой величиной порядка (Л)2, У1г > 0. Компоненты функций поглоще-* * ния vlг и восстановления ulr имеют вид:

v;г(^0 = У1г(Y,t)p;(l)(Y,t), l,г = 1,.,So; l Ф г,

Щг ^, t) = | Уг ^',0Р;(г)(Y',t)q1г (^11 Y',t)dY', (3)

u;l (Y, t)=v;l (Y, t) = 0,

где qlr ^, 11 Y', t) - условная функция плотности вероятности восстановления г-й реализаций из 1-й; p1*(l) (У, t) - условная функция плотности фазовых координат системы при L(t)=l (звездочка означает, что речь идет о непоглощенных реализациях, т.е. отсутствует нормировка).

Вид функции qlr в (3) определяется физическим содержанием задачи и характеризует

начальные условия при восстановлении процесса в r-м состоянии при переходе из l-го состояния. Могут иметь место различные условия восстановления, определяемые функциями qlr. В частности, если

qr (Y, 11 Y', t) = S(Y - Y'), (4)

то восстановление точное ("жесткое", "без потерь"), т.е. конечные условия процесса в l-м состоянии совпадают с начальными в r-м состоянии. Если условная плотность восстановления

qlr (Y, 11 Y', t) = S(Y -Yr)(t)), (5)

то имеют место несвязанные условия восстановления. При восстановлении процесс всегда начинается с заданной функции времени y(r) (t). Если условная плотность вероятности не зависит от предыдущего состояния

qlr (Y, 11 Y', t) = ¥{r)(Y), (6)

то имеет место общий случай процесса с несвязанными условиями восстановления.

В общем случае условная функция qlr (Y, 11 Y', t) может иметь произвольный закон распределения фазового вектора Y(t) в состоянии r при заданном Y'(t) в состоянии l.

Так как каждая структура описывается системой стохастических дифференциальных уравнений, поэтому численный алгоритм должен в себя включать: численное решение СДУ, а также моделирование моментов смены структуры и номера новой структуры. Численные методы решения СДУ в смысле Стратоновича можно найти в работе [3]. В рассматриваемом случае распределение моментов смены структуры определяется интенсивностями переходов. Так как интенсивности переходов зависят от вектора состояния, то моделирование моментов смены, согласно алгоритму из [2], осуществляется по методу максимального сечения.

В работе [4] построен новый экономичный способ моделирования последовательности независимых дискретных случайных величин с помощью лишь одного случайного числа, равномерно распределенного в интервале (0,1). Использование этого способа моделирования, а также лемм, доказанных в работе [5], позволило построить более экономичную модификацию метода максимального сечения.

Построим модифицированный алгоритм решения систем со случайной структурой при распределенных переходах, используя модифицированный метод максимального сечения, построенный в работе [5]. Применение этого метода требует выполнения следующих условий: vu (Y, t) < vJ = const, i = 1, к, S0, i * l.

Модифицированный алгоритм численного моделирования перехода из l-го состояния для систем с условной марковской структурой при распределенных переходах:

Пусть в момент tk система находилась в l-м состоянии и вектор состояния равен Yk.

Моделируем вспомогательную случайную величину а1 - равномерную на интервале (0,1) и заводим счетчик, полагая z=1.

Моделируем возможный момент выхода из l-го состояния: tk+1 = tk + т, где т - случайная величина с плотностью распределения p(x) = vm exp(-vlmx), vm = '^uMV’u (по формуле т = - ln а/ vlm , а - равномерная на интервале (0,1) случайная величина).

Вычисляем возможный r -й номер новой структуры, распределенный с вероятностью

pr (x)=vm /vm, r * l.

На интервале [tk, tk+1 ] численным методом для СДУ решаем уравнения (1) для l-й структуры, находим Yk+1 - вектор состояния системы в момент времени tk+1 (шаг численного метода должен быть согласован с интенсивностью перехода, например h < 0.1/ vj1).

tk := tk+1 , Yk := Yk +1 ; Z := z(1 vlr (Yk , tk )/vlr ) .

Проверяем условие смены структуры: если 1 - а1 > z, то идем на 8); иначе идем на 3).

Меняем номер структуры на r-й; вычисляем Yk согласно заданной условной плотности восстановления qlr.

Следует отметить, что данный алгоритм построен для произвольных систем с условной марковской структурой при распределенных переходах. В более простом, частном случае, когда интенсивности перехода постоянны vlr (Y, t) = vJ = const (такие системы называются системами с независимой марковской структурой при распределенных переходах), в алгоритме будет отсутствовать пункт 7, так как условие, проверяемое в этом пункте, будет выполнено автоматически.

Обозначим Yt - кусочно-линейный процесс, полученный по значениям Yk, а также введем обозначения для функционалов от решения

f (h) := f ( Yt), J(h) := Ef (h), J := Ef (Y).

Для оценки некоторого функционала J от решения (1), численным алгоритмом моделируется N траекторий процесса Y(t) и величина J(h) оценивается средним арифметическим полученных выборочных значений f(h):

Jn (h) = f (i)(h)/N, EJn (h) = J(h).

Погрешность оценки JN (h) определяется величиной

E I J - JN (h) l<l J - J(h) I + E I J(h) - JN (h) l< Chp +JVf (h) / VN,

где E - математическое ожидание; V - дисперсия, а p - порядок слабой сходимости используемого численного метода решения стохастических дифференциальных уравнений (в пятом пункте алгоритма). Таким образом, для уравнивания статистической и детерминированной погрешностей целесообразно полагать N = O(h-2p).

Испытание построенного алгоритма проведем на примерах из [1], для которых удалось записать аналитические формулы для математического ожидания решения. В работе [1] в качестве примеров систем со случайным периодом квантования приведены скачкообразные СДУ с пуассоновской случайной мерой, зависящей от вектора состояния (в этом случае в СДУ (1) добавится пуассоновская составляющая, заменяющая условную функцию плотности вероятности восстановления). Чтобы продемонстрировать влияние пуассоновской составляющей на поведение решения СДУ, рассмотрим простейшее уравнение первого порядка

dy(t) = jev(dffx dT), (7)

г

где характеристическая мера П, задающая пуассоновскую случайную меру v, определяется через неотрицательную функцию п следующим образом:

П(Я, t, y(t-)) = \п(6,t, y(t- ))d6, B еГ. (8)

B

Обозначим

№(t, У) = f П(в, t, y)d6, h(6, t, y) =п(в, t, y)/ ju(t, y) (9)

функции, характеризующие пуассоновскую случайную меру v. Пусть эти функции не зависят от y. Тогда, задавая по-разному функцию h(e), можно получить скачкообразные процессы с различными размерами скачков. Так если

h(e) = S(e-X), (10)

то все скачки одинаковы. При Л > 0 процесс y(t) возрастающий, а при Л < 0 - убывающий.

При

N N

к(в) = £рШ(в-я,), £р! = !• (Ц)

г=1 г=1

где Л - некоторые числа, можно получить модель случайного процесса с дискретным числом состояний, но более сложной структуры.

Рассмотрим СДУ (7), когда функции (9), характеризующие пуассоновскую меру V, зависят от у.

Пример 1. Рассмотрим СДУ (7) при

Ш(в-Л), I/ у < 0]

к(в,у) = 1*,* ^ п\, Л> 0. (12)

[Ш(в + Л), I/ у > 0

Получаем процесс с двумя возможными состояниями: у0 и у1 = у0 + Л, если - Л < у0 < 0; у0 и у1 = у0 - Л, если 0 < у0 < Л.

Пример 2. Рассмотрим СДУ (10), когда к(в, у) задана как в примере 1, но функция ц также зависит от у :

| а, I/ у < 0 ]

Му) = 1 ' Л, Л> 0, а,в> 0. (13)

[в, г/ у > 0

В этом случае, если: -Л < у0 < 0, получаем процесс у(г) с двумя возможными состояниями у0 и у1 = у0 + Л. При а < в в среднем процесс у(г) будет чаще принимать значение у0, а при а > в - значение у1. При в = 0 процесс у (г), приняв значение у1, далее не меняется. Точные значения первого т(г) и второго d (г) моментов имеют вид

т(г) = у0 Р1 (г) + у Р2(г), d (г) = у0 р^г) + у2 р2(г), (14)

где

Р1(г) = -0-+ (р,(0)------а-)е~'а'в, р2(г) = -в + (р1(0)----------а_)е-а+в . (15)

а+в а + в а+в а+в

При приближенном вычислении математического ожидания Е% случайной величины

Ь с конечной дисперсией = о2 по формуле

1 N

Ь = N £ % (16)

при заданном уровне доверия (1 - £) имеет место соотношение

о 1 ^(£)

Р(! Ь - Е% \<7(£)^) | е~у,Чу = 1 - £ , (17)

^ 42ж _/(£)

где у(е) - константа, определяемая выбором величины £. При £ = 0.003 имеем у(е) = 3 , а при £ = 0.3 получаем у(е) = 1.

Рассмотренные примеры являются примерами систем со случайным периодом квантования с интенсивностями переходов, зависящими от вектора состояния системы. В уравнении (7), соответствующем уравнению модели (1) отсутствуют коэффициенты сноса и диффузии. Поэтому при численной реализации построенного алгоритма, пункт 5 отсутствует, и погрешность численного решения состоит только из статистической составляющей.

Оба примера были просчитаны предложенным алгоритмом, использующим модифицированный метод максимального сечения. Оценивалось математическое ожидание решения в узлах временной сетки с шагом к = 0.1 на интервале [0,Т]. Моделировалось N = 106 траекторий и полагалось р1 (0) = 1. При численных расчетах использовался «генератор» псевдо-

случайных чисел RAND [6] с модулем 240 и множителем 517. Длина периода данного датчика составляет 2 38. Данный датчик рекомендуется для численных расчетов, в которых требуется последовательность псевдослучайных чисел не длиннее, чем 237. Расчеты проводились на РС Intel Celeron, 2.02 ГГц, 768 МБ.

Ниже приводятся точные стационарные значения оцениваемого функционала в момент времени T = 20. В таблице с результатами численных экспериментов указаны полученные абсолютные погрешности 8 оценки математического ожидания в последней точке, и приведены доверительные интервалы для оценки математического ожидания при у(е) = 1. Из таблицы видно, что полученные оценки, попадают в доверительные интервалы.

Точные стационарные значения оцениваемого функционала

m

8+V D#(20) 4n

Пример 1 2, = ?0 =1 0 = 0 0.00032 + 0.001

Пример 1 4, = ^0 = 3 0 = 1 0.00064 + 0.00224

Пример 2 2, = ^0 =1 а = 5, в II 0 -1/3 0.00011 + 0.001

Пример 2 2, = ^0 = 3 а = 10, в = 10 1/3 0.00021 + 0.00191

Временные графики m(t), m(t) не приводятся ввиду их сильного совпадения.

Основные достоинства построенного алгоритма состоят в следующем:

1) использование модифицированного метода максимального сечения позволило уменьшить время счета за счет уменьшения числа обращений к «генератору» псевдослучайных чисел;

2) кроме уменьшения трудоемкости вычислений, уменьшение количества используемых значений псевдослучайных чисел снизило конструктивную размерность алгоритма, связанную с многомерной равномерностью используемых псевдослучайных чисел;

3) статистическое соответствие оценок метода из [2] и модифицированного метода является дополнительным критерием удовлетворительности используемого «генератора» псевдослучайных чисел.

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

Работа выполнена при финансовой поддержке РФФИ (проекты N° 09-01-00798, N° 11-01-00282).

ЛИТЕРАТУРА

1. Артемьев В.М. Управление дискретными системами со случайным периодом квантования/ В.М. Артемьев, А.В. Ивановский. М.: Энергоатомиздат, 1986. 97 с.

2. Averina T.A. Algorithm of statistical simulation of dynamic systems with conditional Markovian change of structure // Proc. Intern. Conf. Comput. Math. Novosibirsk: ICM&MG Publisher, 2002. V.1. P. 196-200.

3. Аверина Т.А. Новое семейство численных методов решения стохастических дифференциальных уравнений / Т.А. Аверина, С.С. Артемьев // ДАН СССР. 1986. Т. 288, № 4. С.777-780.

4. Аверина Т.А. Новые алгоритмы статистического моделирования неоднородных пуассоновских ансамблей // Журнал вычислительной математики и математической физики. 2010. Т.50, №1. С. 16-23.

5. Аверина Т. А. Алгоритмы точного и приближенного статистического моделирования пуассоновских ансамблей / Т. А. Аверина, Г.А. Михайлов // Журнал вычислительной математики и математической физики. 2010. Т.50, № 6. С. 1005-1016.

6. Ермаков С.М. Статистическое моделирование / С.М. Ермаков, Г.А. Михайлов. М.: Наука, 1982. 320с.

Аверина Татьяна Александровна -

кандидат физико-математических наук, старший научный сотрудник лаборатории «Численного анализа стохастических дифференциальных уравнений (СДУ)» Института вычислительной математики и математической геофизики Сибирского отделения РАН; доцент кафедры «Вычислительной математики» Новосибирского государственного университета

Статья поступила в редакцию 10.07.11, принята к опубликованию 15.07.11