Финаев В.И. Аналитические и имитационные модели

Подождите немного. Документ загружается.

101

Подпрограммы

WWOD и WIWOD предназначены для

реализации интерфейса имитационной программы. В

подпрограмме

WWOD осуществляется ввод элементов

массива

Р0[I], массива Р[I,J], а также параметров функций

распределения

(t). Определяется T=0 и TZ. В

подпрограмме

OРRZ0 (см. блок 3) определяется начальное

состояние, в подпрограмме

OРRZ определяются состояния

в процессе моделирования смены состояний.

Подпрограмма

STAT1 предназначена для набора

статистических данных о частотах пребывания

марковского процесса в состояниях.

В подпрограмме

OРRTAU определяется время

пребывания марковского процесса в состоянии,

определенном в подпрограмме

OРRZ. Подпрограмма

STAT2 предназначена для набора статистических данных о

времени пребывания марковского процесса в состояниях.

Рассмотрим реализацию подпрограмм

OРRZ0, OРRZ,

OРRTAU.

Подпрограммы

OРRZ0, OРRZ имеют такой же вид, как

и аналогичные подпрограммы для дискретной цепи

Маркова (см. рис. 4.24 и рис. 4.25). Выходным параметром

подпрограммы

OРRZ является индекс J состояния z

(T) в

текущем такте моделирования

Схема алгоритма подпрограммы

OРRTAU будет иметь

вид алгоритмов генерации случайной величины,

рассмотренных в разд. 4.5, метод обратных функций, метод

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

теорем. О схемной реализации подпрограмм

STAT1 и

STAT2 будет сказано ниже.

Существует способ моделирования вложенной цепи

Маркова, который называется схемой развернутой

рекуррентной имитацией. Для этого способа

математическая модель определена:

102

- распределением вероятностей состояния

(0) в момент

;

- матрице

Λ, имеющей вид

11 12 1n

21 22 2n

n1 n2 nn

...

... ... ... ...

...

λλ λ

, (4.14)

где

- интенсивностей переходов марковского процесса

из состояния

в состояние z

Способ моделирования сводится к следующему.

Пусть в момент

t состояние марковского процесса

z(t)=z

. Известен момент времени t

, в который марковский

процесс должен покинуть состояние

. Для определения

следующего состояния

i+1

и момента t

i+1

перехода в

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

случайных величин

по показательному закону

распределения с параметрами

i-й строки матрицы ||λ

||.

Затем определяется минимальное значение из полученной

совокупности независимых случайных величин

, т.е.

ij ki

k=1,2,...,n

τ =minτ

. (4.15)

Следующий момент изменения состояния будет

i+1

+τ

, а состояние, в которое переходит система в

момент

i+1

, будет z

, где z

- состояние, при котором

было минимальным из всех значений совокупности

независимых случайных величин

{

}

. Физическая

интерпретация такого подхода к моделированию

марковского процесса с непрерывным временем позволяет

разрабатывать имитационные модели реальных систем.

Схема алгоритма моделирования вложенной цепи

Маркова по данному способу приведена на рис. 4.29.

103

WWOD

Начало

T=T+1

OPRZ0

MINTAU

WIWOD

Конец

T<TZ

OPRTAU

STAT1

Рис. 4.29

Отличие алгоритма по схеме развернутой рекуррентной

имитации (см. рис. 4.29) от алгоритма, приведенного на

рис. 4.28, состоит в блоках 4 и 5, в которых реализована

генерация случайных величин

и формулы (4.15).

На рис. 4.30 приведена схема алгоритма подпрограммы

OРRTAU для данного способа. На рис. 4.31 приведена

схема алгоритма подпрограммы

MINTAU.

104

TAU[J]=-(1/λ

)lnP

СЕN(P)

Начало

Конец

J=J+1

J<n

J=0

Рис. 4.30

В подпрограмме OРRTAU (см. рис. 4.30) входным

параметром является индекс

I предыдущего состояния

марковского процесса.

Для каждого значения

J в блоке 4 по методу обратных

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

TAU[J], отождествляющей собой описанную выше

величину

. Таким образом, в данной подпрограмме будет

сформирован массив

TAU[J].

В подпрограмме

MINTAU (см. рис. 31) входным

параметром являются элементы массива

TAU[J], из

которых необходимо выбрать минимальный элемент.

В блоках 1 – 5 осуществляется выбор минимальной

величины (

MIN) из массива TAU[J]. Происходит это

следующим образом. При переменной массива

J=1 первый

105

элемент массива

TAU[1] определяется, как минимальный

элемент

MIN (см. блок 1).

MIN=TAU[J]

Начало

Конец

J=J+1

TAU[J]<MIN

J=1, MIN=TAU[1]

J<n

J=J+1

TAU[J]=MIN

J=0

Рис. 4.31

Затем второй элемент массива TAU[2] сравнивается со

значением

MIN (см. блок 3). Если второй элемент массива

TAU[2] меньше значения MIN, то идентификатору MIN

присваивается значение

TAU[2] (см. блок 4). Если условие

TAU[2]<MIN не выполняется, то рассматривается условие

TAU[3]<MIN для элемента массива TAU[3] и т.д.

В блоках 6 – 8 определяется индекс

J элемента массива

TAU[J], равного числу MIN. Это индекс J определяет

индекс состояния

, в который переходит марковский

процесс.

106

5. ОБРАБОТКА РЕЗУЛЬТАТОВ

МОДЕЛИРОВАНИЯ НА ЭВМ

5.1. Выбор числа опытов

При разработке имитационных моделей для

исследования случайных объектов существует задача

выбора числа опытов (объема выборки). Это непростая

задача, т.к. во-первых, необходимо обосновать

достоверность результатов моделирования и связать

достоверность с точностью, а во-вторых, существуют

события, вероятность появления которых очень мала (

Р→0)

или, наоборот очень велика (

Р→1). Для обоснования

объема выборки при имитационном моделировании

применяют аналитические подходы [8,12,13].

В инженерной практике известны критерии для оценки

погрешности. Непрерывную (аналоговую) величину

х(t) в

ряде инженерных задач рассматривают как дискретную,

как показано на рис. 5.1.

(

)

x(t)

(t)

x(t)

Рис. 5.1

Если х

(t) - результат измерения непрерывной величины

х(t), то для любого момента t текущая погрешность

107

дискретизации этой непрерывной величины определится:

δ(t)=х(t)-х

(t). Выбор критерия оценки δ(t) зависит от

назначения величины

х(t). Известны следующие критерии.

Критерий наибольшего отклонения имеет вид

)t(x)t(x)t(max

−=δ≤δ

Δ∈

. Критерий применим, если

известны априорные сведения о сигнале в форме условия

Липшица

)tt(l)t(x)t(x

−≤−

, где l -некоторая константа.

Среднеквадратичный критерий приближения

определяется по формуле

()

dt)t(x)t(x

dt)t(

σ=−

=δ

=σ

∫∫

∗

Среднеквадратичный критерий применим для функций,

интегрируемых в квадрате. Использование

среднеквадратичного критерия связано с усложнениями,

например аппаратуры измерения, по сравнению с

критерием наибольшего отклонения.

Интегральный критерий как мера отклонения

х(t) от х

(t)

имеет вид

dt)t(

μ≤=δ

=μ

∫

Если моделируются случайные процессы, то выше-

названные критерии не применимы.

Выбор количества реализаций зависит от того, какие

требования предъявляются к результатам моделирования.

Пусть для оценки случайной величины

A, оцениваемой

по результатам моделирования

х, выбирается величина х

являющаяся функцией от

х. Значения х

будут отличаться

от

A в силу случайных факторов, т.е. можно связать

точность оценки, теоретическое значение случайной

величины

A и её статистическую оценку х

в виде формулы

|A-х*|<ε, (5.1)

108

где

ε - точность оценки.

В силу того, что каждый результат

моделирования

случайной величины

A также является случайной

величиной, то вероятность того, что неравенство (5.1)

выполняется, будет достоверностью точности оценки

случайной величины

A, т.е. справедлива формула

Р(|A-х*|<ε)=α. (5.2)

Воспользуемся сформулированным критерием (5.2) для

определения точности результатов методом

статистического моделирования. Пусть цель

моделирования - вычисление вероятности

р появления

события

A. Количество ε наступления события A в

реализации процесса является случайной величиной,

принимающей значение

=1 с вероятностью р и значение

=0 с вероятностью 1-р. Пример случайного появления

события

A показан на рис. 5.2.

А А А А

Рис. 5.2

Математическое ожидание случайной величины ε

определится по формуле

M[ε]=х

р+х

(1-р)=р. (5.3)

Если

=1, как появление события A, а х

=0, как не

появление события, то значение

M[ε] совпадает с

вероятностью

р наступления события A.

Дисперсия определится по формуле

D[ε]=[х

-M[ε]]

р+[х

-M[ε]]

(1-р)=р(1-р). (5.4)

При выполнении имитационного моделирования

оценкой вероятности

р является частость m/N наступления

события

A при N реализациях, m - число испытаний, в

109

которых событие

A наступило. Частость m/N определяется

формулой

∑

p* = m N = 1 N ε , ε [0,1]

n=1

(5.5)

где ε

- количество наступлений событий A в реализации с

номером

Из формул (5.3), (5.4) и (5.5) можно определить

математическое ожидание и дисперсию частости

m/N

m m p(1- p)

M[ ]=p; D[ ]=

NNN

. (5.6)

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

частость

m/N при NÆ∞ имеет распределение, близкое к

нормальному. Поэтому для каждого значения

достоверности

α (вероятности) можно выбрать из таблиц

нормального распределения такую величину (значение

случайной величины)

, что точность ε будет равна

.]Nm[D

tε =

(5.7)

Подставив в формулу (5.7) значение D из формулы (5.6),

получим

)p1(p

tε

−

(5.8)

Из формулы (5.8) можно определить количество

реализаций

N, необходимых для получения оценки m/N с

точностью

ε и достоверностью α:

)p1(p

−

(5.9)

В формуле (5.9) неизвестны величины N и р, т.к.

вероятность

р определяется, исходя из ее оценки m/N, а N -

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

моделирования для определения

N поступают следующим

образом. Выбирают

=50-100. По результатам N

110

реализаций определяют

m/N

, т.е. осуществляют

примерную оценку вероятности

р=m/N

. Затем по формуле

(5.9) окончательно выбирают

N, принимая р=m/N.

Другим случаем является оценка по результатам

моделирования среднего значения некоторой случайной

величины. Пусть непрерывная случайная величина

A имеет

среднее значение

и дисперсию σ

. В реализации с

номером

i случайная величина A принимает значение х

. В

качестве оценки для среднего значения (математического

ожидания)

A используется среднее арифметическое

∑

x= x

i=1

В силу центральной предельной теоремы при NÆ∞

будет иметь приблизительно нормальное распределение с

математическим ожиданием

и дисперсией σ

/N, поэтому

точность определится по формуле

ε =t σ N

Число реализаций определится по формуле

22 2

N=t σε

. (5.10)

Так как в формуле (5.10) неизвестными являются число

реализаций

N и среднеквадратичное отклонение σ

, то

также выбирают

=50-100. По результатам N

реализаций

определяют оценку дисперсии, а затем по формуле (5.10)

окончательно выбирают

N. Количество реализаций N в

формуле (5.9) зависит от

р, а в формуле (5.10) от σ

Целесообразно так строить моделирующий алгоритм,

чтобы методом моделирования оценивались параметры

величин, имеющих возможно меньшую дисперсию, или

вероятности случайных событий, не близкие к

0,5.

Вероятности не должны быть также близки к

0 или 1, т.к. в