Моделирование при проектировании информационно-управляющих систем

W(j)=A()ej(),(2.18)

где A()=|W(j)| ()=argW(j).

Функция А=А() называется амплитудной характеристикой системы и представляет собой отношение амплитуды гармонического сигнала y(t) к амплитуде гармонического сигнала U(t).

Функция =() называется частотной характеристикой системы, показывает, на сколько выходной сигнал y(t) при данной частоте сдвинут по фазе относительно входного сигнала U(t).

W(j) можно представить в виде

W(j)=Р()+jQ()

где Р()=ReW(j) Q()=ImW(j).

Р() и Q()-соответственно вещественная и мнимая частотные характеристики системы с передаточной функцией (р).

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

Примем аппроксимацию передаточной функции (р) в виде

где, l=0,1,2 - порядок астатизма объекта.

Определение k и ai осуществляется методом наименьших квадратов по инверсной с, которая при l=0 имеет вид

-1(р)=bnрn+…+b1р+b0, где ,(2.19)

n-порядок аппроксимации.

Очевидно, если определим b0 и bi, то параметры (р) будут определены по формулам:

k=1/b0, Ai=bi/b0=bik

Подставим в (2.19) р=j и учитывая, что n5, запишем

W-1(j)=U()+jV()=(b0-b22+b44)+j(b1-b32+b54),

где U()=b0-b22+b44,(2.20)

U()/=b1-b32+b54.

Если частотная характеристика задана прямоугольными координатами, то

U()+jU()=,

где Re()=Р(), Im()=Q()

или тогда.(2.21)

Если частотная характеристика задана полярными координатами, то

U()+jU()=

или

Формулы (2.20) и (2.21) выведены для объектов с самовыравниванием, когда l=0. Аналогичные формулы можно вывести для l=1 и l=2.

Уравнения регрессии (2.20) сходны между собой, поэтому их можно записать в общем виде



=с1+с2i2+с3i4.

Обозначение - это величина, которая определяется через экспериментальные значения i и неизвестные коэффициенты Сj. Индекс i показывает, что соответствующие величины относятся к i-й точке частотной характеристики.

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

,(2.22)

можно вычислить значения коэффициентов Сj, Zi - экспериментальная величина, определяемая по координатам i-й точки частотной характеристики.

Если положить Zi=Ui, то согласно (2.20)

b0=С1, b2=-С2, b4=С3.(2.23)

Если положить Zi=Ui/i, то согласно (2.20)

b1=С1, b3=-С2, b5=С3.(2.24)

Взяв частные производные в (2.22) с учетом (2.20), получим систему нормальных уравнений для определения коэффициентов Сj:

(2.25)

Т.к. мы рассматриваем систему с порядком аппроксимации n=5, то для вычисления b0,b1,…,b5 необходимо сделать следующее. При Zi=Ui находим С1,С2,С3 из системы (2.25), а затем из формулы (2.23) b0,b2 и b4. При Zi=Vi/i вычислим С1,С2,С3 из системы (2.25), а затем из (2.24) определим.

Зная b0,b1,…,b5, определим параметры k и ai передаточной функции (р).

Эксперимент проходит по схеме, представленной на рис.2.4.

Индикаторы фиксируют амплитуды Ai1 и Ai2 - входного и выходного токов.

Регистрирующий прибор записывает U(t) и Y(t) в переходном и устанавливаемом режимах:

,

где k1 и k2 - масштабные коэффициенты.

тогда

Р(i)=A(i)сos(i),(2.26)

Q(j)=A(i)sin(i).

Подставляем (2.26) в формулы (2.21) и получаем U() и V() для определения С1,С2,С3 через Zi в уравнении (2.25)

Р(i)Q(i)U(i)V(i)/iZiС1,С2,С3b1,b2,b3,b4,b5.

Рис.2.4

3. Стохастические модели объектов

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

Если входные параметры объекта, смена состояний объекта или его выходные параметры описываются случайными распределениями, то эти объекты относятся к классу стохастических. При моделировании поведения данных объектов применяется аппарат теории вероятностей, а для идентификации параметров моделей применяется аппарат математической статистики. Рассмотрим виды моделей, которые могут быть применены для описания стохастических объектов.

3.1.1 Определение случайных функций

Течение случайного процесса описывают некоторой функцией (), где -- аргумент функции со значениями из множества . Функцию (), наблюдаемую в некотором опыте, соблюдая определенный комплекс условий, называют выборочной функцией или реализацией случайного процесса.

Если множество произвольно, то вместо термина «случайный процесс» удобнее пользоваться термином «случайная функция». Название «случайный процесс» применимо в тех случаях, когда параметр интерпретируется как время. Если аргумент случайной функции является пространственной переменной, то функцию называют случайным полем.

Определение. Моделью случайного процесса называют случайную функцию (), заданную на множестве , принимающую действительные значения и описываемую семейством распределений [5]

F12... n(х1,х2,...,хn), i, i=1,2,...,n, n=1,2,...,

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

F12...nn+1,...,n+1(х1,х2,...,хn,...,+,...,+)= F12.. n(х1,х2,...,хn)

F12...n(х1,х2,...,хn) =Fi1i2...in(хi1,хi2,...,хin),

где i1,i2,...,in -- любая перестановка индексов 1, 2,..., n.

Набор функций F12...n(х1,х2,...,хn) называется конечномерными распределениями случайной функции.

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

Последовательность функций 1(),2(),…,n() возможно заменить векторной функцией (), компонентами которой служат случайные функции i(), (i=1,2,…,n).

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

Если f12...n(х1,х2,...,хn) -- плотность функций распределения F12,...n(х1,х2,...,хn), то

.

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



1(), 2(), … n(), i0 >, i=1,n, n=1,2,...,

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

где M -- символ математического ожидания, u1,u2,...,uk -- вещественные числа.

Если существует плотность конечномерного распределения, то модель в виде характеристической функции является преобразованием Фурье плотности распределения.

3.1.2 Корреляционные функции

Исчерпывающую характеристику модели стохастического объекта в виде случайной функции в широком смысле дает семейство конечномерных распределений. Однако решение многих теоретико-вероятностных задач зависит только от небольшого числа параметров, характеризующих входящие в задачу распределения. Наиболее важными числовыми характеристиками распределений являются их моменты. В теории случайных функций роль моментов распределений играют моментные функции.

Определение. Модель случайной функции (i), i в виде моментной функции задается отношением

Mj1,j2,…,jn(1,2,...,n)=M{[(1)]j1,...,[(n)]jn},

если математическое ожидание в правой части равенства имеет смысл при всех i, i=1,n. Величина q=j1+j2+...+jn называется порядком моментной функции.

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

при u1=u1=…=un=0.

Кроме моментных функций в качестве моделей часто рассматривают центральные моменты функции

mj1,j2,…,jn(1,2,...,n)=M{[(1)-m(1)]j1,...,[(n)-m(n)]jn},

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

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

m()=m1(1)=M(),

R1(1,2)=m1(1,2)=M{[(1)-m(2)][(2)-m(2)]}.

Функции m() называются средним значением, а R1(1,2) - корреляционной функцией.

При 1=2= корреляционная функция дает дисперсию () величины (), R1(1,2)=2().

Величину

называют коэффициентом корреляции случайных величин (1) и (2).

3.1.3 Классификация моделей случайных процессов

Случайные процессы делятся на следующие широкие классы: гауссовы процессы; процессы с независимыми приращениями; стационарные в широком смысле; марковские процессы.

Модели на базе гауссовых случайных функций. Важную роль во многих прикладных вопросах играют случайные функции, конечномерные распределения которых являются гауссовыми (нормальными). Определение многомерного гауссового распределения следующее.

Определение. Случайный вектор =(1,2,...,n) имеет гауссово (нормальное) распределение, если характеристическая функция распределения представима в виде

(u)=M{eхр[j(u,)]}=eхр[j(m,u)-0,5R(u,u)],

где m=(m1m2...mn), u=(u1u2...un) - векторы, R - неотрицательно-определенная вещественная симметричная матрица, R=||rij||, i,j=1,n. Здесь (, ) обозначает скалярное произведение векторов и , так, что

Модель процессов с независимыми приращениями. Пусть T - конечный отрезок T=[0,a] или T=[0,].

Определение. Случайный процесс {(t), tT} со значениями в евклидовом пространстве Rn называется процессом с независимыми приращениями, если для любых n, таких, что 0<t1<t2<...<tn, случайные векторы (0), (t1)-(0),...,(tn)-(tn-1) - взаимно независимы.

Вектор (0) называется начальным состоянием (значением) процесса, а его распределение -- начальным распределением процесса. Чтобы задать процесс с независимыми приращениями в широком смысле, достаточно задать начальное распределение Р0(B)=Р{(0)B} и набор распределений Р(t,н,B) - распределений вектора Р{(t+н)-(t)}B.

Процесс с независимыми приращениями называется однородным, если распределения вектора (t+н)-(t) не зависят от t, Р(t,н,B)=Р(н,B).

Модель процессов, стационарных в широком смысле. Стационарные процессы - это такие процессы, теоретико-вероятностные характеристики которых не изменяются со временем.

Пусть T=[0,a] или T=[0,).

Определение. Модель случайного процесса (в широком смысле) {(t), tT} со значениями в Rn называется стационарной, если для любого n и любых t1,t2,...,tт, таких, что tk+tT, (k=1,n), совместное распределение случайных векторов, описывающих случайный процесс (t1+t),...,(tn+t), не зависит от t.

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

Определение. Случайный процесс (t), t>0 со значениями в пространстве Rn называют процессом, стационарным в широком смысле, если M[(t)]2< и M[(t)]=m=сonst, M[(t)-m][(s)-m]=R(t-s), (t>s), где R(t) -- непрерывная матричная функция.

Функцию R(t) называют корреляционной (матричной) функцией процесса (t).

В качестве примера стационарных в широком смысле процессов можно рассмотреть колебания со случайными параметрами.

3.1.4 Модели марковских процессов

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

Представим себе систему, которая может находиться в разных состояниях. Возможные состояния образуют множество Х, называемое фазовым пространством. Пусть система эволюционирует во времени. Состояние системы в момент времени t обозначим через хt. Если хtB, где BХ, то будем говорить, что система в момент времени t находится во множестве B.

Предположим, что эволюция системы носит стохастический характер, т.е. состояние системы в момент времени t не определяется однозначно через состояние системы в моменты времени s, предшествующие t, где s<t, а является случайным и описывается теоретико-вероятностными законами.

Пусть Р(s,х,t,B) - вероятность события хtB (s<t), при условии, что хs=х. Функцию Р(s,х,t,B) называют вероятностью перехода рассматриваемой системы. Под системой без последействия понимают систему, для которой вероятность попадания в момент времени t во множество B, при полностью известном движении системы до момента времени s (s<t), по-прежнему равна Р(s,х,t,B) и, таким образом, зависит только от состояния системы в последний момент времени.

Обозначим через Р(s,х,u,y,t,B) условную вероятность события хtB при гипотезах хs=х, хu=y (s<u<t). В соответствии с общими свойствами условных вероятностей имеет место равенство

.(3.1)

Для системы без последствия естественно предположить, что

Р(s,х,u,y,t,B)=Р(u,y,t,B).

Тогда равенство (3.1) примет вид

.(3.2)

Соотношение (3.2) называется уравнением КолмогороваЧепмена. Это уравнение определяет модель марковского процесса.

Пусть {Х,B}-некоторое измеримое пространство. Функцию Р(х,B), хХ, BB, удовлетворяющую условиям:

а) Р(х,B) при фиксированном х является мерой на B и Р(х,Х)=1;

б) при фиксированном B Р(х,B) является B - измеримой функцией от х будем называть стохастическим ядром.

Пусть I - некоторый конечный или бесконечный полуинтервал (отрезок).

Семейство стохастических ядер {Рst(х,B)=Р(s,х,t,B), s<t, (s,t)II}, удовлетворяющих уравнению Колмогорова-Чепмена (3.2), будем называть марковским семейством стохастических ядер.

Определение. Моделью марковского процесса в широком смысле называется совокупность следующих объектов:

- измеримое пространство {х, B};

- полуинтервал I (отрезок) действительной оси;

- марковское семейство стохастических ядер {Рst(х,B), s<t, (s,t)II}.

Семейство ядер Рst(х,B)=Р(s,х,t,B) называют вероятностью перехода марковского процесса, пространство {х,B} - фазовым пространством системы, точка множества I интерпретируется как моменты времени, а величина Рst(х,B)=Р(s,х,t,B) - как условная вероятность того, что система в момент времени t окажется во множестве B, если в момент времени s она находилась в точке х фазового пространства (s<t).

Классификация марковского процесса. Дискретные случайные процессы, обладающие марковскими свойствами, называются цепями Маркова. В фазовом пространстве простейшими марковскими процессами являются процессы со счетным или конечным числом состояний. В фазовых пространствах выделяются следующие классы марковского процесса.

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

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

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

Марковские процессы в конечномерном пространстве, аппроксимируемые на малых промежутках времени произвольным процессом с независимыми приращениям.

3.2 Методы имитации случайных факторов

3.2.1 Датчики случайных чисел

Базовой последовательностью случайных чисел, используемой в ЭВМ для формирования случайных элементов различной природы, с различными функциями распределения, является совокупность случайных чисел с равномерным законом распределения.

Строго говоря, с применением ЭВМ получить последовательность случайных величин (СВ) с равномерным распределением не представляется возможным. Поэтому, если считать, что число разрядов ЭВМ равно k, а случайное число Х будет сформировано согласно формуле [6]

где i=0, рi=0,5; i=1, рi=1-0,5, то множество, состоящее из чисел ={i/(2k-1)}, принимает значение i/(2k-1) (i=0,1,2,...,2k-1) с вероятностью Р=1/2k.

Такое распределение чисел Х называется квазиравномерным в интервале [0,1], причем математическое ожидание и дисперсия определяются следующими соотношениями:

.(3.3)

Из формулы (3.3) видно, что математическое ожидание M[х] точно совпадает с генеральным средним для равномерного распределения в интервале [0,1], а дисперсия при k асимптотически стремится к дисперсии для равномерного распределения, равной 1/12.

Практически при k>15 обеспечивается требуемая точность в имитационных исследованиях.

При выводе выражения (3.3) предполагалось, что Х формируется на основе случайных чисел , принимающих значения (0,1) с вероятностью Р=1/2, для чего в машине должен существовать случайный генератор, дающий строго случайные последовательности чисел с соответствующим распределением. В ЭВМ такого генератора нет и случайные числа вырабатываются программным путем, в силу чего они, строго говоря, не являются случайными, т.к. формируются на основе вполне детерминированных преобразований, поэтому их называют псевдослучайными. Однако, если при моделировании число обращений к программному датчику случайных чисел оказывается меньше периода, измеряемого числом различных случайных чисел, то такая периодичность программного датчика не оказывает существенного влияния на результаты моделирования.

Методы получения псевдослучайных квазиравномерных чисел программным путем можно разбить на две основные группы: а) аналитические; б) методы перемешивания.

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

Хr=(Хr-1, Хr-2,..., Х0).

Простейшим примером может служить метод вычетов, в котором используется следующее рекуррентное соотношение:

Хi+1=bХi(mod M),

где выражение bХi(mod M) означает остаток от деления произведения bХi на число M; Хi+1 - очередное случайное число; Хi - предыдущее случайное число; b - некоторая константа; M - число, определяющее значение получаемых случайных чисел.

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

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

Проверочные тесты для проверки качества серии квазиравномерных псевдослучайных чисел следующие.

Рассмотрим тест частот. Отрезок [0,1] разбивается на m (обычно 10-20) равных интервалов. Полученные эмпирические частоты ni/N (i=), сравнивают с теоретическими вероятностями 1/m. Согласие проверяется по критерию 2, т.к. стохастическая формула [7,8]

подчиняется распределению 2 с (m-1) степенями свободы, где N -- объем выборки.

Рассмотрим тест пар частот. Рассматриваются последовательные пары случайных чисел. Квадрат [0,1]х[0,1] делится на m2 частей. Каждая пара случайно попадает в одно из m2 делений квадратной таблицы.

Пусть дана серия чисел х1,х2,...,хn. Если пары образовать в виде (х1,х2), (х3,х4)..., то пары взаимно независимы, эмпирические частоты (а их число равно m2) сравниваются с теоретическими вероятностями равномерного распределения 1/m2. Функция

распределена по закону 2 с (m2-m) степенями свободы, где nij -- число попаданий в (i,j)-ю клетку таблицы mхm, N/2 -- объем выборки пар случайных чисел [7,8].

Более сложная ситуация возникает, если пары образовать в виде (х1,х2)(х2,х3)... . Этот метод образования пар более выгодный, т.к. полнее использует выборку чисел, но из-за зависимостей пар случайная величина 2определится по формуле

]

и будет иметь распределение 2 с (m2- m) степенями свободы.

3.2.2 Имитация случайных событий

Пусть в результате эксперимента должно наступить с вероятностью рi одно из несовместимых событий A1, A2,..., An, которые образуют полную группу событий.

Разбиваем отрезок [0,1] на n частей длиной Р1,Р2,...,Рn, при этом точки деления отрезка имеют следующие координаты:

Пусть теперь х -- очередное число от генератора случайных чисел.

Если lk-1х<lk, то считаем, что произошло событие Ak. Действительно

Р(Ak)=Р(lk-1<х<lk)=lk-lk-1=Рk.

3.2.3 Имитация непрерывных случайных величин

[6]. Существует несколько методов, основанных на преобразованиях равномерно распределенных случайных чисел в числа с заданным законом распределения.

Метод обратной функции. Если х - случайная величина, равномерно распределенная на отрезке [0,1], то случайная величина y, являющаяся решением уравнения

, (3.4)

плотность распределения f(y). Этот метод позволяет вывести правило генерирования случайного числа, имеющего произвольное непрерывное распределение f(y):

- вырабатывается случайное число х генератором случайной равномерной последовательности;

- случайное число yi, имеющее распределение f(y), находится из решения уравнения (3.4).

Графическая иллюстрация метода обратных функций приведена на рис.3.1.

Метод ступенчатой аппроксимации. Зависимость плотности распределения f(y) представляется графически в интервале изменения y от a до b. Если случайная величина задана на [0,], то произведем усечение распределения с заданной точностью. Разобьём отрезок [a,b] на n частей таких, что

где ai (i=0,n) -- координата точки разбиения.



Рис3.1

Тогда вероятность того, что случайная величина y попадет в один из интервалов определится по формуле

Т.е. попадание на любой отрезок равновероятно.

На каждом из интервалов функция f(y) аппроксимируется так, чтобы значение f(y) в каждом интервале было постоянным. Тогда координата случайной точки может быть представлена как yi=ai+сi, где сi- - расстояние от левого конца интервала. В силу ступенчатой аппроксимации сi является равномерно распределенной случайной величиной на интервале [ai,ai+1] (i=0, n-1).

Правило имитации в этом случае сводится к следующему:

- получаем два числа х1 и х2 от генератора равномерно распределенных чисел;

- с помощью х1 находим индекс i=[nх1] интервала, где [nх1];

- целая часть числа nх1, причем [nх1]<nх1;

- с помощью числа х2 находим сi= х2(ai+1-ai);

- находим случайное число, имеющее интересующий нас закон распределения f(y), по формуле 

.

Использование предельных теорем. Для получения нормального закона распределения используется свойство сходимости независимых величин к нормальному распределению. Для получения нормального распределения чисел с параметрами my=0, y=1 удобен искусственный прием, основанный на центральной предельной теореме теории вероятностей.

Для этого в качестве исходных чисел возьмем n равномерно распределенных на отрезке [1,-1] чисел, получаемых из интервала [0,1] по правилу zi=2хi-1.

Сформируем величину z согласно следующей формуле:

.

По центральной предельной теореме при достаточно большом значении n величина z может считаться нормально распределенной с параметрами

нормирование величины z и получим

(3.5)

Величина y будет иметь нормальное распределение с my=0, y=1.

Установлено, что при n>8 формула (3.5) дает хорошие результаты.



3.2.4 Имитация марковского процесса

Рассмотрим марковский процесс с дискретным временем перехода из одного состояния в другое.

Пусть система находится в некотором множестве состояний zi и (ziZ), i=1,2,…,n. Динамика процесса будет заданной, если задан граф перехода либо матрица переходов ||Рij||, где Рij - вероятность перехода из состояния zi в состояние zj в некоторый дискретный момент времени kt (k=1,2,3,...).

Рij не зависит от времени, т.е. процесс является однородным. Начальные условия пусть заданы в виде распределения: Р1(0), Р2(0),...,Рn(0), где Рi(0) - вероятность нахождения процесса в i-ом состоянии при t=0.

Моделирование такого процесса основано на принципе имитации системы случайных событий.

Определяем числовые границы

Выбирается начальное случайное число х0, и определяется номер r начального состояния z0, для которого будет справедливо условие

выбирается случайное число х1, которое также сравнивается с границами

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

Рассмотрим моделирование марковского процесса, если переход из состояния в состояние может произойти в любой момент времени. Марковский процесс в этом случае определяется:

- распределением вероятностей состояния Рi(0) в момент t0;

- набором интенсивностей ij переходов процесса из состояния zi в состояние zj.

Один из возможных способов моделирования сводится к следующему.

Пусть z(t)=zk, и известен момент времени tk, в который системы должна покинуть состояние zk.

Для определения следующего состояния zk+1 и момента перехода в следующее состояние tk+1 генерируется ряд независимых случайных величин zkzi по показательному закону распределения с параметрами k-ой строки матрицы ||ki ||. Затем определяется

.

Следующий момент изменения состояния будет tk+1=tk+zi, а состояние, в которое переходит система в момент tk+1, будет zi, где zi - состояние, при котором zi было минимальным.

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



4. Алгоритмизация процессов функционирования систем

4.1 Моделирующие алгоритмы

Для моделирования любого объекта, заданного при помощи математической модели, а также в виде последовательности процедур, имитирующих отдельные элементарные процессы, необходимо построить соответствующий моделирующий алгоритм. Структура программы вычислений, составленная применительно к типу ЭВМ, зависит от вида алгоритма и от характеристик ЭВМ. Моделирующий алгоритм необходимо записать в таком виде, который бы отражал в первую очередь особенности его построения без излишних второстепенных деталей.

Создание моделирующего алгоритма - этап исследования, когда уже решены все вопросы выбора математического аппарата для исследования.

Необходимо сделать запись алгоритма независимо от характеристик ЭВМ. Способы представления моделирующего алгоритма следующие: запись алгоритмов при помощи операторных схем; запись в языках программирования; использование методов прикладных программ.

Применительно к имитационному моделированию это называется: операторные схемы моделирующих алгоритмов (ОСМА); языки программирования; универсальные имитационные модели.

ОСМА содержит последовательность операторов, каждый из которых изображает достаточно крупную группу элементарных операций. Эта запись не содержит развернутых схем счета, но достаточно полно отражает логическую структуру моделирующего алгоритма. ОСМА не учитывает особенности системы команд. Это происходит при построении программы.

Требования к операторам: оператор должен иметь ясный смысл, связанный с природой моделируемого процесса; любой оператор может быть выражен последовательностью элементарных операций.

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

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

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

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

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

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

Для изображения операторной схемы моделирующих алгоритмов удобно пользоваться арифметическими и логическими операторами.

Арифметические операторы производят действия, связанные с вычислениями. Обозначаются A14 - арифметический оператор №14.

Свойство арифметического оператора состоит в том, что после выполнения изображенных им операций передается действие другому оператору. - передача управления от А14 к А16 (графически отображается стрелкой).

Логические операторы предназначены для проверки справедливости заданных условий и выработки признаков, обозначающих результат проверки.

Свойство логического оператора состоит в том, что после его реализации управление передается одному из двух операторов алгоритма, в зависимости от значения признака, вырабатываемого логическим оператором. Обозначается в виде Рi, а графически в виде круга или ромба, внутри которого символически записывается условие.

Изображение передачи управления - Р352212. Если условие выполняется, то управление передается оператору №22, если нет -- то оператору №12.

Для операторов всех классов обозначение передачи управления оператора, следующему непосредственно за ним, опускается.

Передача управления данному оператору от других операторов обозначается 16,14A18. Оператору A18 управление передается от операторов №16 и №14..

Обозначение оператора, обозначающего окончание вычислений, - Я.

Пример. Рассмотрим решение уравнения х2+рх+q= 0,

Введем операторы:

A1 -- вычисление р/2;

A2 -- вычисление р2/4-q;

A3-- вычисление ;

Р4 -- проверка условия D0;

A5 -- определение действительных корней х12=-(р/2)R;

A6 -- определение мнимых корней х12=-(р/2)jR;

Я -- окончание вычислений и выдача (х1,х2).

Операторная схема алгоритма

A1 A2 A3 Р46 A57 A6, 5Я7.

Операторную схему алгоритма можно заменить рисунком алгоритма, вид которого показан на рис.4.1.

Операторные схемы алгоритмов позволяют перейти от схематического изображения алгоритма к его записи в виде формулы.

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

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

Важнейшие типы операторов следующие. Вычислительные операторы (операторы счета) описывают сколь угодно сложную и громоздкую группу операторов, если она удовлетворяет требованиям, предъявляемым к операторам алгоритма (подготовленность исходных данных, передача управления только одному оператору в операторных схемах моделирующего алгоритма). Обозначаются Ai.

Операторы формирования реализаций случайных процессов решают задачу преобразования случайных чисел стандартного вида в реализации случайных процессов с заданными свойствами. Обозначаются i.

Операторы формирования неслучайных величин формируют различные константы и неслучайные функции времени. Обозначаются Fi.

Счетчики подсчитывают количества различных объектов, обладающих заданными свойствами. Обозначаются Ki.

Рис.4.1

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

Процесс функционирования сложной системы можно рассматривать как последовательную смену ее состояний, описываемых характеристиками z1(t),z2(t),...,zn(t) в n-мерном фазовом пространстве. Задачей моделирования является построение функций zi(t), а также вычисление некоторых величин, зависящих от этих функций. Математическая модель связывает характеристики состояний системы zi(t) с ее параметрами и временем при начальных условиях для t0 : zi(t0).

Пусть существует сложная система с детерминированными характеристиками. Преобразуем соотношения математической модели к виду, удобному для вычисления значений zi(t+t) по известным значениям zi(t), при t<t. Выделим ячейку для фиксации текущего времеи t и назовем ее часами (таймером). При t0 определим , при t0+t определим zi(t0+t) и т.д. Если шаг t0, то получим приближенные значения zi(t).

Рассмотрим сложную систему со стохастическими параметрами.

Состояние zi(t) и соотношения математической модели определяют распределение вероятностей величин zi(t+t), состояния также могут быть случайными и задаваться соответствующими распределениями вероятностей.

Структура моделирующего алгоритма для таких систем та же. Но вместо состояния z(t+t) необходимо вычислить распределение вероятностей для возможных состояний.

В соответствии с заданным распределением вероятностей выбирается одно из состояний . Затем при (t0+t) вычисляется условное распределение вероятностей состояний при условии . По жребию определяется состояние zi(t0+t) и т.д.

Принцип построения моделирующего алгоритма, позволяющий определить последовательные состояния сложной системы через некоторые интервалы времени, иногда называют "способ t-моделирования" (неэкономичен с точки зрения расхода машинного времени).

При рассмотрении некоторых сложных систем можно обнаружить неравномерность состояний системы в заданном интервале времени t.

Выделяются два типа состояний: обычные состояния, в которых система находится почти все время; особые состояния, характерные для системы в некоторые изолированные моменты времени, совпадающие с моментами поступления в систему входных сигналов от внешней среды, выхода характеристики zi(t) на границу области существования и т.д. Координаты zi(t) в эти моменты времени могут изменяться скачком.

Очевидно, что моделирующие алгоритмы, построенные по принципу t-моделирования, оказываются не эффективными.

Для данных систем моделирующие алгоритмы строятся по способу "особых состояний". Они отличаются от принципа t только тем, что включают в себя процедуру определения момента времени, соответствующего следующему особому состоянию по известным характеристикам данного или предыдущего состояния.

При моделировании обработки заявок в системах массового обслуживания строится моделирующий алгоритм по способу последовательной проводки заявок. Идея этого способа состоит в последовательном воспроизведении истории отдельных заявок в порядке поступления их в систему: алгоритм обращается к сведениям о других заявках лишь в том случае, если это необходимо для решения задачи о дальнейшем порядке обслуживания данной заявки. Оператор имеет сложную логическую структуру, экономичен по машинному времени.

4.3 Фиксация и обработка результатов моделирования

При реализации моделирующего алгоритма на ЭВМ вырабатывается информация о состоянии систем, которая представляет собой исходный материал для определения приближенных искомых величин.

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

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

Оценка Р(A) вероятности Р(A) события A



Р*(A)=m/N,

где m -- число случаев наступления событий A, N -- число реализаций.

Если событие принимает значения в некоторой области величин, то область значений n случайной величины разбивается на отрезки так, что n={n1,n2…nm} Оценка вероятностей возможных k-ых значений случайной величины определяется

Р*k(A)=mk/N,

где mk- число значений случайной величины в интервале nk.

Среднее значение случайной величины

,

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

Оценкой S2* дисперсии случайной величины определится

Эта формула неудобна. Известна упрощенная формула, согласно которой

т.е. для определения S2* достаточно накапливать значения 

и .

Для оценки корреляционного момента K случайной величины и с возможными значениями хk и yk применяется формула

эта формула преобразуется к виду

требующему запоминания небольшого количества величин.

Иногда искомыми величинами являются математическое ожидание и корреляционные функции случайного процесса, например х(t).

Интересующий интервал (0,T) разбивается на части с шагом t. Накапливают значения хk(t) реализаций случайного процесса х(t) для фиксированных моментов времени ti. Затем вычисляют оценки для математического ожидания по формуле

Оценки для корреляционной функции вычисляются по формуле

где u и s «пробегают» все значения ti. Или 

4.4 Точность. Количество реализаций

Если х*(t) - результат измерения некоторой величины х(t), то текущая погрешность дискретизации определится: (t)=х(t)-х*(t).

Выбор критерия оценки (t) зависит от назначения величины х(t). Известны следующие критерии.

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

Критерий применим, если известны априорные сведения о сигнале в форме условие Липшица , где l -некоторая константа.

Среднеквадратичный критерий приближения определяется по формуле

.

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

Интегральный критерий как мера отклонения х(t) от х*(t) имеет вид

.

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

Выбор количества реализаций зависит от того, какие требования предъявляются к результатам моделирования. Пусть для оценки параметра a, оцениваемого по результатам моделирования хi, выбирается величина х*, являющаяся функцией от хi, х* будет отличаться от a в силу случайных факторов, т.е.

|a-х*|< ,(4.1)

где - точность оценки. Вероятность того, что неравенство (4.1) выполняется, называется достоверностью точности оценки х*, т.е.

Р(|a-х*|< )=. (4.2)

Воспользуемся сформулированным принципом (4.2) для определения точности результатов методом статистического моделирования.

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

Количество наступления события A в реализации процесса является случайной величиной, принимающей значение х1=1 с вероятностью р и значение х2=0 с вероятностью 1-р.

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

M[]=х1р+х2(1-р)=р.(4.3)

Это совпадает с вероятностью наступления события A.

Дисперсия определится



D[]=[х1-M[]]2р+[х2-M[]]2(1-р)=р(1-р).(4.4)

Оценкой вероятности р является частость m/N наступления события A при N реализациях. Частость m/N можно представить в виде

(4.5)

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

Из формул (4.3), (4.4) и (4.5) можно определить математическое ожидание и дисперсию частости m/N

(4.6)

В силу центральной предельной теоремы вероятностей частость m/N при N имеет распределение, близкое к нормальному. Поэтому для каждого значения достоверности можно выбрать из таблиц нормального распределения такую величину t, что точность будет равна

(4.7)

Подставим в (4.7) значение D из (4.6), тогда

(4.8)

Из (4.8) можно определить количество реализаций N, необходимых для получения оценки m/N с точностью и достоверностью :



(4.9)

В практике моделирования вероятность р обычно неизвестна. Для определения N поступают следующим образом. Выбирают N0=50-100. По результатам N0 реализаций определяют m/N, а затем окончательно выбирают N, принимая р = m/N0.

Другим случаем является оценка по результатам моделирования среднего значения некоторой случайной величины. Пусть случайная величина имеет среднее значение A и дисперсию 2.

В реализации с номером i она принимает значение хi. В качестве оценки для среднего значения (мат.ожидания) A используется среднее арифметическое

В силу центральной предельной теоремы при Nх будет иметь приблизительно нормальное распределение с математическим ожиданием А и дисперсией 2/N. Поэтому точность определится

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

(4.10)

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



5. Моделирование систем с использованием типовых математических схем

5.1 Модели систем массового обслуживания

5.1.1 Общие сведения

Моделирование объекта с применением математического языка систем массового обслуживания (СМО) предусматривает в процессе формализации выделение понятий: заявка (требование), поток заявок, прибор обслуживания, очередь на обслуживание, дисциплины выбора на обслуживание, закон обслуживания, поток обслуженных заявок, поток потерянных заявок.

Известна классификация, которая произведена, исходя из характеристик СМО [9,10]. СМО классифицируют следующим образом.

По потокам заявок СМО делятся на СМО с однородным потоком и приоритетные СМО.

По дисциплинам обслуживания СМО делятся на СМО с дисциплиной FIFO (первый пришел - первый обслуживается), СМО с дисциплиной LIFO (последний пришел - первый обслуживается), СМО со случайным выбором на обслуживание.

Исходя из того, каким временем на ожидание располагает заявка, СМО делятся на СМО с отказами, если эта величина времени ожидания равна нулю, смешанные СМО, если время ожидания является конечной величиной (СМО с ограниченной очередью), СМО с ожиданием, если время ожидания является бесконечной величиной (с бесконечной очередью).

По количеству и структурному расположению приборов обслуживания СМО делятся на одноканальные, если имеется один прибор обслуживания, n-канальные, если имеются n параллельно расположенных приборов, m-фазные СМО, если имеются последовательно расположенные m приборов, СМО смешанной структуры. Классификация может быть осуществлена, исходя из математических законов, описывающих математические модели потока входных заявок и времени обслуживания.

Моделирование систем с применением схем СМО предусматривает определение выходных параметров и параметров состояния, которые могут быть представлены как показатели эффективности СМО.

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

5.1.2 Модель входного потока заявок

Входной поток заявок характеризуется начальным моментом времени t0, моментами времени ti поступления i-ых заявок, случайными величиными i - интервалами времени между заявками, i=ti-ti-1. Модель потока в общем виде представляет собой конечномерную функцию распределения вероятностей:

F(х1,х2,...,хn)=Р{1<х1, 2<х2,..., n<хn}.

Если i - величины детерминированные, то имеем дело с равномерным потоком заявок. Можно задать для каждого i плотности распределения fi(х). В том случае, когда плотность совместного распределения будет определяться как f(х1,х2,...,хn)=f1(х)f2(х)...fn(х), получим поток Пальма с ограниченным последействием.

Известны три характеристики для классификации входных потоков:

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

- стационарный поток, если вероятность поступления k-заявок за интервал времени (t0,t) не зависит от выбора момента t0;

- поток без последействия, если вероятность появления k-заявок внутри некоторого интервала не зависит от появления заявок до момента начала этого интервала.

Простейший поток (поток Пуассона) удовлетворяет всем трем условиям. Для этого потока вероятность поступления k-событий за время t определится

Функция распределения времени поступления между двумя заявками определяется экспоненциальным распределением - A(t)=1--t.

Наиболее часто применяется при моделировании экспоненциальное распределение и распределение Эрланга. Функция распределения плотности вероятности интервалов между заявками для эрланговского потока r-го порядка определится

Если r=0, то получаем экспоненциальное распределение. Эрланговские распределения описывают модели потоков с последействием.

5.1.3 Модель времени обслуживания

Моделями времени обслуживания могут служить функция и плотность распределения вероятности длительности обслуживания.

При исследовании прибора обслуживания необходимо определить эмпирическую плотность распределения длительности обслуживания, а затем ее аппроксимировать известными теоретическими распределениями. Наиболее часто применяемые - нормальное, постоянное, экспоненциальное распределения и распределение Эрланга.

5.1.4 Модель Эрланга

Определенный интерес при моделировании СМО представляет подход, при котором исследуются изменения в системе за сколь угодно малый отрезок времени. Составляются уравнения в частных приращениях, от которых затем осуществляется переход к дифференциальным уравнениям. Рассмотрим вывод дифференциальных уравнений, известных как модель Эрланга [9,10].

Будем рассматривать одноканальную СМО с бесконечной очередью, с ожиданием, пуассоновсим потоком заявок и экспоненциальным временем обслуживания. Поток ординарный, простейший, функция распределения интервалов между заявками является экспоненциальной. Модель смены состояний можно представить в виде графа, приведенного на рис.5.1.

Рис.5.1

Составим уравнения Эрланга в частных приращениях, которые будут отображать те изменения, которые произошли в системе за сколь угодно малое время t.

Из графа состояний (см. рис.5.1) мы видим, что следует выделить начальное состояние, когда число заявок в СМО n=0 и состояния с числом заявок в СМО n1.

Вероятность Р0(t+t) того, что СМО к моменту t+t останется в нулевом состоянии, определится из анализа полной группы событий:

- в момент времени t система была в нулевом состоянии и за время t заявки не поступали;

- в момент времени t система была в единичном состоянии (в СМО была одна заявка) и за время t обслуживание заявки окончилось.

Вероятность Р0(t+t) определится

Р0(t+t)=Р0(t)(1-t)+Р1(t)t,(5.1)

1-t -- вероятность непоступления заявки в СМО за время t, t - вероятность окончания обслуживания заявки за время t.

Вероятность Рn(t+t) того, что к моменту времени t+t система будет в n-м состоянии, определится из рассмотрения следующей полной группы событий:

- в момент времени t в системе было n-1 заявок и за время t поступила заявка;

- в момент времени t система была в n-м состоянии и за время t заявки в СМО не поступили и обслуживание не окончено;

- в момент времени t в системе была n+1 заявка и за время t обслуживание заявки было окончено.

Вероятность Рn(t+t) определится

Рn(t+t)=Рn-1(t)t+Рn(t)[1-(+)t1+Рn+1(t)t, (5.2)

t - вероятность поступления заявки за время t; 1-(+)t - вероятность непоступления заявки в СМО и неокончания обслуживания заявки за время t.

Уравнения (5.1) и (5.2) представляют собой модель рассматриваемой СМО в виде уравнений Эрланга в частных приращениях.

От уравнений в частных приращениях перейдем к дифференциальным уравнениям.

Для этого Рn(t) из правой части перенесем в левую, разделим каждую часть на t и определим предел при t 0. Получим уравнения:

(5.3)

Уравнения (5.3) представляют собой модель исследуемой СМО в виде дифференциальных уравнений Эрланга для нестационарного случая.

Так как поток заявок, поступающих в систему, отвечает условиям стационарности, то значение производных можем приравнять к нулю. Получим модель СМО в виде уравнений Эрланга для стационарного режима

Р1=Р0, n=0, (1+)Рn=Рn+1+Рn-1, n1,(5.4)

где /= - коэффициент использования системы.

Решение системы уравнений (5.4) будет иметь следующий вид:

Рn=nР0, Р0 =(1-), Рn = n (1-),

где Рn - вероятность того, что в СМО будет n-заявок.

Затем могут быть определены такие характеристики СМО, как математическое ожидание числа заявок в СМО, математическое ожидание числа заявок в очереди и другие.



5.1.5 Исследование модели пуассоновского процесса с помощью производящих функций

[10]. Будем считать, что на вход СМО поступает пуассоновский поток заявок с интенсивностью и вероятностью Рn(t) того, что за время t в СМО поступит n-заявок. Делаем предположение, что при сколь угодно малом отрезке t вероятность поступления заявки определится через t . Вероятность непоступления заявки определится как 1-t. Поток является ординарным. Можно записать уравнение в частных приращениях. Вероятность того, что к моменту времени t+t в системе не будет заявок, определится через вероятность того, что в системе в момент времени t не было заявок, и за отрезок времени t заявки в систему не поступили:

Р0(t+t)=Р0(t)(1-t).(5.5)

Вероятность того, что к моменту времени 1+t в СМО будет n заявок, определится как вероятность того, что в момент t в СМО было n заявок, и за время t заявка не поступила, или к моменту времени t в СМО была n-1 заявок, и за время t поступила еще одна заявка:

Рn(t+t)=Рn(t)(1-t)+ Рn-1(t) t.(5.6)

После проведения преобразований уравнений (5.5) и (5.6), аналогичных преобразованиям уравнений (5.1), (5.2), получим дифференциальные уравнения

(5.7)



Рассмотрим решение уравнений (5.7) с применением производящих функций. Производящая функция Р(z,t) для функции Рn(t) определится

Вероятность Рn(t) получим из производящей функции после того, как продифференцируем ее n раз и положим z=0.

При решении уравнения в частных приращениях начало отсчета времени выбирается произвольно даже после того, как в систему поступило i заявок. Будем считать, что при t=0 в СМО есть i-заявок. В этом случае Рn(0)=0, если ni и Рn(0)=1, если n=i. Таким образом,

Если умножим дифференциально-разностное уравнение (5.5) на zn, а дифференциально-разностное уравнение (5.6) на z0 и просуммируем по всем значениям n, так что

то получим, что сумма в левой части равна



а сумма первых членов правой части равна - Р(z,t). Просуммировав вторые члены правой части по n, получим

.

Если в правой части выделить множитель z, то всю сумму можно записать в виде zР(z,t). Таким образом, система приводится к линейному дифференциальному уравнению для производящей функции, которое имеет вид

Решение этого уравнения при постоянном значении z (поскольку оно не зависит от t) имеет вид

Р(z,t)=Сe(z-1)t.

Допустим, что к моменту t=0 не поступило ни одного требования, тогда Р(z,0)=1, так как i=0. Таким образом, С=1 и

Р(z,t)=e(z-1)t.

Как говорилось выше, Рn(t) определится

...