CC BY
8
ОБ ИМИТАЦИОННОМ МОДЕЛИРОВАНИИ СЛУЧАЙНЫХ ВЕЛИЧИН С ПОМОЩЬЮ ИНТЕНСИВНОСТИ1
Зверкина Г. А.2
(ФГБУН Институт проблем управления им. В.А. Трапезникова РАН, Москва) Кошелев А. А.3 (Московский государственный университет им. М.В. Ломоносова, Москва)
Как известно, случайные величины обычно моделируются с помощью функции, обратной к функции распределения. Однако часто при моделировании поведения сложных технических моделей исследователь сталкивается с изменением параметров моделируемых случайных величин. Часто в этих случаях меняется интенсивность завершения определенного случайного периода. Следовательно, важно иметь возможность моделировать случайную величину в условиях (случайного) изменения интенсивности этой случайной величины. В данной статье предлагается новый метод моделирования случайной величины по её интенсивности. Текст снабжён иллюстративным материалом результатов численных экспериментов, подтверждающих эффективность предложенного метода. Естественно, этот текст - только начало работы. В дальнейшем будут уточнены технические детали эффективности применения предложенного метода.
Ключевые слова: интенсивность случайной величины, численное моделирование, вопросы оптимизации моделирования.
1. Введение
Как хорошо известно, значительное количество известных моделей в теории массового обслуживания, в теории надёжности и теории сетей массового обслуживания описывается линейчатыми процессами (см. [3]).
1 Работа выполнена при частичной финансовой поддержке РФФИ, грант №20-01-00575А
Авторы признательны М.П.Фархадову за ценное обсуждение содержания статьи.
2Галина Александровна Зверкина, к.ф.-м.н., доцент (zverkina@gmail.com).
3Александр Анатольевич Кошелев, (koshelev030698@yandex.ru).
Как правило, большинство таких процессов являются эрго-дическими, т.е. с течением времени распределение состояния такого процесса (слабо) сходится к некоторому стационарному инвариантному (относительно соответствующих переходных вероятностей) распределению.
Обычно для стационарного распределения можно вычислить некоторые важные характеристики, однако в значительном количестве случаев теоретическое определение скорости сходимости распределения изучаемой системы к стационарному распределению оценить или вычислить невозможно.
Поэтому для получения оценок скорости сходимости распределения исследуемой системы к стационарному распределению может быть использовано имитационное моделирование поведения системы с помощью программных средств.
Напомним данное в [3] определение линейчатого процесса (с некоторой модификацией для того, чтобы кусочно-линейный процесс был марковский):
Определение 1. Кусочно-линейным марковским процессом
называется случайный процесс Xt = (и(1),у(1)), определяемый следующим образом. Пространство состояний процесса X - это множество пар , где г - элемент конечного или счётного
пространства, а щ - вектор = (ь\,... ьщ), где |г| ^ 0 является «размером» базового состояния г; ^ О. Поведение процесса Xt описывается так.
Пусть Xt = (г,у), у = (уг,.. .ущ). Тогда с вероятностью \iiXt) х рц(Х^ёИ за время произойдёт случайный переход Xt в базовое состояние ]. После перехода в состояние ] новое значение компоненты случайно и определяется измеримой по у функцией распределения (ф.р.)
В((0)(х\у) = Р (ф + ё*) < х\у(г) = у, и СО = I, и (* + = з} .
Вероятность того, что за малое время Ъ произойдёт более одного случайного перехода, есть о(Ъ). При этом если V= = V(£), то + = у(1) + аг где а.г = (ац,... ,ацц) -
вектор с неотрицательными компонентами, если на интервале (¿, £ + не было случайных переходов (см. [3]).
Таким образом, поведение кусочно-линейного процесса определяется случайными величинами (сл.в.), для которых обычно известна интенсивность окончания того или иного периода. Кроме того, часто бывает известна функция распределения периодов кусочно-линейного процесса - в этом случае может применяться классический метод моделирования случайных величин с помощью обратной функции к функции распределения.
То есть случайная величина £ с функцией распределения Р(х) = Р {£ ^ х} может быть промоделирована как
(1) £ = р-1(и),
где
р-1(у) = ^{% е к : Р(х) > у}, у е М,
а и - равномерно распределённая случайная величина на интервале [0; 1).
Такое моделирование одной системы надёжности проводилось в работе [4]; при этом было обнаружено, что начиная с некоторого момента за счет накопления ошибки точность вычисления падает и становится примерно равной около 10-4, что демонстрирует статья [4]. Случайные величины в работе [4] моделировались с точностью до 15 значащих цифр по схеме, описанной выше по формуле (1).
Естественно, в дальнейшем следует учитывать тот факт, что имитационное моделирование непрерывных случайных величин фактически является моделированием дискретных случайных величин с шагом, определяемым используемым программным продуктом.
В значительном количестве изучаемых моделей ТМО и смежных областей распределения работы/восстановления/интервала между заявками и пр. известны. Но не всегда даже для известных распределений можно найти явное выражение для обратной функции Р-1( ) - например, для часто встречающегося в приложениях гамма-распределения. Также в ряде сложных сто-
хастических систем «поведение» элементов системы может существенно зависеть от общего состояния системы (т.е. состояний других элементов системы) и не может быть определено до начала моделирования поведения сложной системы. В этом случае будем рассматривать случайную величину (длину периода восстановления) как случайную величину, заданную с помощью (случайно) меняющейся интенсивности окончания того или иного периода (отказа, ремонта, появления требования в системе, обслуживания требования и др.) меняется в зависимости от полного состояния изучаемой системы (см., например,
Поэтому целью настоящей работы является представление алгоритма моделирования случайной величины с помощью её интенсивности.
2. Случайные величины и интенсивности
Пусть при исследовании поведения некоторого случайного процесса нам необходимо учесть поведение некоторого объекта, на поведение которого может влиять состояние всей исследуемой системы.
Пусть в некоторый момент времени £ известно, что еще не закончился период работы (ремонта, ожидания появления нового требования и пр.). Если распределение этого периода £ имеет непрерывную функцию распределения Р(в), то вероятность того, что на бесконечно малом промежутке времени (Ь, Ь + ДЬ) произойдет окончание интересующего нас периода, имеет вид
[2, 8, 9, 10, 11, 12, 13, 14]).
Р {£ е (М + Д1) | £>П = Р {{ < г + Д^ - Р {£<*}_ ^(г + Дг) - ^(¿)
Р {С > Ц
1 - р(г)
При Д£ ^ 0 имеем
(2) Р {{ е + Дг) | £ > г} = Х(г)Дг + о(Дг),
где
№
1п(1 - ^(I)),
Р(¿) — функция распределения, р(Ь) - плотность распределения. Функция Х(Ь) называется интенсивностью случайной величины (см., например, [7]). Мы предполагаем, что интенсивность может зависеть от времени и от некоторых дополнительных параметров, связанных с поведением других элементов сложной системы. Заметим, что
Р(5) = 1-ехр J , р(в) = А(«)ехр J \(и)йи
Отсюда видно, что п.в. (почти всюду) положительная случайная величина с абсолютно непрерывным распределением, являющаяся периодом некоторого процесса восстановления, определяется однозначно либо функцией распределения, либо плотностью, либо интенсивностью - подробнее см. [1, 6].
3. Идея предлагаемого метода моделирования случайных величин
В данной работе предлагается для моделирования поведения сложных систем (массового обслуживания, надёжности и пр.), поведение которых определяется кусочно-линейными процессами, использовать не «классический» метод моделирования случайных величин по схеме (1), а опираясь на формулу (2).
А именно, рассматривается последовательность независимых одинаково распределённых величин {^¿}, равномерно распределённых на [0; 1).
Время с момента Ь = 0 до произвольного достаточно большого времени Т разбивается на интервалы точками 0 < т\ < < Т2 < ... < тп = Т .В рассматриваемый момент времени £ = Тг < Т предполагается, что время пребывания некоторого параметра рассматриваемой сложной системы (обслуживания, надёжности и пр.) в данном состоянии равно величине в ^ 0; интенсивность окончания этого периода равна А(з). Тогда изменение этого состояния (т.е. окончание пребывания в данном состоянии)
на интервале времени (¿; ¿+Д£) (где = Тг+1 — ъ) близко к значению А(,§)Д1 Поэтому с вероятностью Р, = Р {Л(s)Дí > Щ в момент Тг+\ закончится рассматриваемый период времени пребывания исследуемого параметра в данном состоянии.
Таким образом, общая схема моделирования случайной величины по её интенсивности выглядит так:
1. Выбираем точки 0 <т1 < т2 < ... <тп = Т.
2. Если выполнено неравенство А(0)(т1 — 0) > и0, то принимаем, что £ = 0; если это не так, то переходим к следующему этапу.
3. Если выполнено неравенство Л(г1)(г2 — т-\) > и1, то принимаем, что £ = т1; если это не так, то переходим к следующему этапу.
4. Повторяем описанную в п.3 процедуру для ^ последовательно до тех пор, пока не случится событие
Чп)(п+1 — Тг) > Ыг.
Полагаем £ = п. Напомним, что Р {Л < а} = а при а € [0; 1].
Очевидно, что для любой сл.в. £ с конечным математическим ожиданием эта процедура конечна с вероятностью 1. Заметим, что при
(3)
У \(з) ёз =
о
описанная процедура конечна с вероятностью 1, но применение её в случае бесконечности Е £ нецелесообразно.
Данный метод удобнее классического (формула (1)) тем, что можно моделировать положительные случайные величины не зависимо от того, насколько сложно (и возможно) вычисление обратной функции к функции распределения случайной величины.
С
Например, для функции интенсивности А(з) = невоз-
3 83 + 1
можно выписать функцию распределения, поскольку интеграл [ С л
аз не вычисляется в явном виде, но сама подын-
У 3 я3 + 1 38
тегральная функция соответствует формуле (3), т.е. задаёт собственную функцию распределения сл.в.
Основной проблемой данного метода является то, что от выбор интервалов An = тг+i — тг в зависимости от значения величины А(т¿) может существенно влиять на точность результатов предложенного метода моделирования случайных величин. Авторы предполагают в дальнейшем более детально исследовать этот вопрос.
Также важно и качество используемого в программе генератора случайных чисел.
4. Результаты моделирования
Ниже представлены результаты моделирования случайных величин по интенсивности как в случае, когда она имеет экспоненциальный момент, так и в случае, когда случайная величина имеет лишь конечное число степенных моментов.
Для этого рассматривались случайные величины с постоянной интенсивностью Ai(s) = Ао = const, что соответствует экспоненциальному распределению F1(s) = 1 — exp(—A0s). Также
рассматривалась переменная интенсивность A2(s) =-, кото-
1 + ^ 1
рая соответствует функции распределения F2(s) = 1 — —-,
(1 + s)
которую далее мы будем называть полиномиальным распределением.
Предложенный выше метод моделирования случайных величин был реализован на языке программирования C в программе mod.c (см. Приложение). С помощью стандартного генератора псевдослучайных целых неотрицательных чисел rand() был сконструирован генератор случайных чисел с равномерным распределением:
#define frand() ((double) rand() / (RAND_MAX + 1.0))
Для того чтобы функция frand() при каждом запуске не давала одинаковые последовательности случайных чисел, она
была дополнена функцией srand с аргументом, равным значению функции time(NULL), которая возвращает текущее время в секундах с 1 января 1970 года [5].
Далее была реализована функция modulus, моделирующая случайную величину по методу, описанному выше, с использованием интенсивности. Чтобы убедиться в корректности работы программы, а также подобрать оптимальный шаг А, были рассмотрены и смоделированы две случайные величины с функциями распределения F1(s) и F2(s) с различными параметрами. В качестве критерия достоверности результатов моделирования сравнивались выборочные и теоретически вычисленные моменты случайных величин, а также гистограммы полученных выборок промоделированных случайных величин и плотностей соответствующих распределений.
Далее обозначаем — случайная величина, распределённая по показательному закону Fi(s), а £2 — по полиномиальному закону F2(s).
Приведем основные вероятностные характеристики этих случайных величин: функцию распределения F(s), плотность распределения p(s), интенсивность A(s), математическое ожидание E£, дисперсию В£. Для показательного распределения:
Fi(s) = 1 - е
F\(s) = 1 - e-X°s,
s, s ^ 0 и, s < 0,
Ai(s) = Aq,
где Aq = const > 0, а для полиномиального
т - 1 в^ - 2 1
= С—1' = С2 - 3С + 2 - (С - 1)2'
где С = const > 2.
5. Приложения
5.1. Приложение 1
По выборке, полученной моделированием для каждого из распределений, были подсчитаны математические ожидания и дисперсии для различных шагов Д. Выборки состояли из 100000 сгенерированных чисел. Далее приведем результаты вычислений. В качестве параметров распределений были взяты А0 = 2, С = 8.
А0 = 2, С = 8, Д = 0,1
Эксперимент N E£i E6 D6
1 0,401235 0,200112 0,028346 0,004858
2 0,399677 0,200059 0,028995 0,004950
3 0,400694 0,201088 0,028509 0,004851
А0 = 2, С = 8, Д = 0,01
Эксперимент N E£i E&
1 0,487648 0,243396 0,132102 0,024996
2 0,488431 0,240377 0,130432 0,024451
3 0,487839 0,244927 0,132280 0,025362
А0 = 2, С = 8, Д = 0,001
Эксперимент N E£i E6 D6
1 0,496119 0,245180 0,141784 0,027176
2 0,498286 0,249556 0,141512 0,027039
3 0,495537 0,245814 0,141822 0,027003
А0 = 2, С = 8, Д = 0,0001
Эксперимент N E£i E&
1 0,466462 0,219541 0,139288 0,025180
2 0,468898 0,218354 0,138768 0,025246
3 0,468516 0,218312 0,139447 0,025323
Управление большими системами. Выпуск 94 _А0 = 2, С = 8, А = 0,00001
Эксперимент N ка юа Е6
1 0,327335 0,107316 0,112233 0,014785
2 0,327366 0,107663 0,112095 0,014777
3 0,328139 0,107909 0,112430 0,014857
Из эксперимента видно, что для выбранных значений Ао = 2 и С = 8 оптимальный шаг для моделирования случайной величины А = 0,001.
В результате моделирования было замечено, что величина шага А может существенно влиять на результат численного эксперимента.
5.2. Приложение 2.
Ниже приводятся гистограммы промоделированных случайных величин с указанным в заготовке распределением (голубые столбики) в сравнении со значениями плотности распределения (100000 реализаций, шаг моделирования А = 0,001).
Рис. 1. Гистограмма промоделированных случайных величин
с распределением Р2 (голубые столбики) в сравнении со значениями плотности распределения (100000 реализаций)
Рис. 2. Гистограмма промоделированных случайных величин
с распределением Р\ (голубые столбики) в сравнении со значениями плотности распределения (100000 реализаций)
5.3. Приложение 3
Текст программы mod.c
i #include <stdio.h>
2#include <stdlib.h> 3#include <math.h> 4#include <time.h>
5
6#define lambdaO 2.0 7#define C 8.0 8#define deltal 0.001 9#define delta2 0.001
io#define frand() ((double) rand() / (RAND_MAX + 1.0))
ii
i2void sort(double *p0, int n); i3double lambda(double t, int num);
i4void modulus(double *mass, double *mE, double *disp,
double delta, int num); i5void nuller(double *mass, int n);
i6
i7void sort(double *p0, int n)
i8 {
19 double first;
20 int j, i;
21 for (i = 0; i < n-1; i++)
for(j = i+1; j < n; j++)
if(p0[i] > p0[j])
24 {
25 first = p0[i]; p0[i] = p0[j];
27 p0[j] = first;
28 }
29 }
30
31 double lambda(double t, int num)
32 {
33 double l;
if(num == 1)
35 {
36 l = lambda0; //1 --- Exponential distribution
37 }
if(num == 2)
39 {
40 l=C/(1+t);//2 --- Polynomial distribution
41 }
42 return l;
43 }
44
45 void modulus(double *mass, double *mE, double *disp,
double delta, int num)
46 {
47 int i;
double randnum, t = 0, tNew = 0, sumt = 0, tsq = 0;
49 for (i = 0; i < 100000; i++)
50 {
while(1)
52 {
53 randnum = frand();
if(delta * lambda(t, num) > randnum)
55 {
56 tNew = t;
57 mass[i] = t;
58 break;
59 }
60 t += delta;
61 }
62 sumt += tNew;
63 tsq += t*t;
64 t = 0;
65 }
66 *тЕ = эш^ / 100000;
67 = (tsq / 100000 - *тЕ * *тЕ) * 100000 / 99999;
68 ге^т;
69 }
70
71 void пи11ег^оиЬ1е *тазз, int п)
72 {
73 int ^
74 £ог^ = 0; i < п; i++)
75 {
76 тазэ^] = 0;
77 }
78 }
79
80 int main(void)
81 {
82 doub1e гv1[100000], гv2[100000], техр[30], тро1[16];
83 doub1e х, тЕ = 0, disp = 0, dt1 = 0.1, 1е^1 = 0, эит = 0;
84 int i, j;
85 Е11Е *£1, *£2;
86 £1 = fopen("out1.txt", "ш");
87 £2 = fopen("out2.txt", "ш");
88 sгand(time(NULL));
89 modu1us(гv1, &тЕ, &disp, de1ta1, 1);
90 pгintf("Exponentia1 distгibution\n de1ta = %1£\П М[ exp] = M[exp_theoг] = %1£\П D[exp] = %1f\t Б[ exp_theoг] = %1£\п", de1ta1, тЕ, 1.0 / 1ambda0, disp,
1.0 / pow(1ambda0, 2));
91 modu1us(гv2, &тЕ, &disp, de1ta2, 2);
92 pгintf("Po1ynomia1 distribution\n de1ta = %1£\П М[ po1] = %1f\t M[po1_theoг] = %1£\П D[po1] = %1f\t D[ po1_theor] = %1£\п", de1ta2, тЕ, 1.0 / (С - 1), disp,
2.0 / (С*С -3*С+2) -1.0/ pow(C - 1, 2));
93
94 soгt(гv1, 100000);
95 soгt(гv2, 100000);
96 nu11er(mexp, 30);
97 nu11eг(mpo1, 16);
98 j = 0;
99 1eft1 = dt1;
foг(i = 0; i < 100000 && j < 30; i++)
101 {
102 1£(^1[1] > 1е^1)
103 {
104 1е^1 +=
105 j++;
106 }
107 шехр^] += 1;
108 }
109 j =0;
110 1е^1 = dt1;
£ог(1 = 0; 1 < 100000 && j < 16; 1++)
112 {
113 1£(^2[1] > 1eft1)
114 {
115 1е^1 += dt1;
116 j++;
117 }
118 mpo1[j] += 1;
119 }
120
121 1е^1 = 0.05;
122 ^г^ = 0; j < 30; j++)
123 {
124 fprintf(f1, "%1f\t %1f\t %1f\t %1f\n", 1eft1, (
mexp[j] / 100000) / dt1, 1ambda0 * exp(-1ambda0 * ( 1eft1)), 1 - exp(-1ambda0 * (1е^1)));
125 1eft1 += dt1;
126 }
127
128 1е^1 = 0.05;
129 ^г^ = 0; j < 16; j++)
130 {
131 fprintf(f2, "%1f\t %1f\t %1f\t %1f\n", 1eft1, (
mpo1[j] / 100000) / dt1, С / pow(1 + 1eft1, С + 1), 1 - 1 / pow(1 + 1eft1, С));
132 1е^1 += dt1;
133 }
134
135 ^^эе^^;
136
137 ге^т 0;
138 }
Литература
1. БОРОВКОВ А.А. Вероятностные процессы в теории массового обслуживания. - М.: Физматлит, 1972. - 368 с.
2. ВЕРЕТЕННИКОВ А.Ю. О скорости сходимости к стационарному распределению в системах обслуживания с одним прибором // Автоматика и телемеханика. - 2013. - №10. -С. 23-35. Engl. tranls.: On the rate of convergence to the stationary distribution in the single-server queuing system // Automation and Remote Control. - 2013. - Vol. 74, Iss. 10. -P. 1620-1629.
3. ГНЕДЕНКО Б.В., КОВАЛЕНКО И.Н. Введение в теорию массового обслуживания. - М.: Наука, 1966.
4. ЗВЕРКИНА Г.А. Об экспоненциальной сходимости коэффициента готовности // Управление большими системами. -2021.-Вып. 90.-С. 5-35.
5. КЕРНИГАН Б.У., РИЧЧИ Д.М. Язык программирования С. 2-е изд. - М.: ООО «И.Д. Вильямс», 2015. - 288 с.
6. КОКС Д., СМИТ В. Теория восстановления. - М.: Советское радио, 1967.
7. ТАРТАКОВСКИЙ А.Г. О последовательном оценивании интенсивности процессов восстановления // Проблемы передачи информации. - 1985. - №21:1. - С. 40-47; Problems Information Transmission. - 1985. - Vol. 21:1. - P. 30-36.
8. VERETENNIKOV A. On convergence rate for Erlang-Sevastyanov type models with infinitely many servers // Theory of Stochastic Processes. - 2017. - No. 1. - P. 88-102.
9. VERETENNIKOV A. On mean-field GI/GI/1 queueing model: existence and uniqueness // Proc. of the Int. Conf. "Analytical and Computational Methods in Probability Theory and Its Applications" (ACMPT-2017), 23-27 October 2017 / Ed.: A.V. Lebedev. - RUDN, Moscow, Russia. - P. 182-186.
10. VERETENNIKOV A. On mean-field (GI/GI/1) queueing model: existence and uniqueness // Queueing Systems. - 2020. -Vol. 94(3). - P. 243-255.
11. VERETENNIKOV A. On Polynomial Recurrence for Reliability System with a Warm Reserve // Markov Processes and Related Fields. - 20i9. - Vol. 25. - P. 745-76i.
12. VERETENNIKOV A. On Recurrence and Availability Factor for Single-Server System With General Arrivals // Reliability: Theory and Applications (RT& A). - 20i6. - Vol. ii, No. 3(42). - P. 49-58.
13. VERETENNIKOV A.YU. On the rate of convergence for infinite server Erlang-Sevastyanov's problem // Queueing Systems. - 20i4. - Vol. 76, Iss. 2. - P. i8i-203.
14. VERETENNIKOV A.YU., ZVERKINA G.A. Simple Proof of Dynkin's Formula for Single-Server Systems and Polynomial Convergence Rates // Markov Processes and Related Fields. -20i4. - Vol. 20. - P. 479-504.
ON SIMULATION OF RANDOM VARIABLES USING INTENSITY
Galina Zverkina, V.A. Trapeznikov Institute of Control Sciences of RAS, Moscow, Cand.Sc., assistant professor (zverkina@gmail.com). Alexandr Koshelev, Moscow State University, Moscow, student (koshelev030698@yandex.ru).
Abstract: As is well known, random variables are usually modelled using the inverse function of the distribution function. However, often when modelling the behaviour of complex technical models, a researcher is faced with a change in the parameters of the simulated random variables. Often in these cases, the intensity of the completion of a certain random period changes (or "failure rate", or "hazard rate"). Therefore, it is important to be able to simulate a random variable under conditions of a (random) change in the intensity of this random variable. This paper proposes a new method for modelling a random variable by its intensity. The text is provided with illustrative material of the results of numerical experiments confirming the efficiency of the proposed method. Naturally, this text is just the beginning of work. In the future, the technical details of the effectiveness of the application of the proposed method will be clarified.
Keywords: random variable intensity, numerical simulation, discussion of simulation optimization.
УДК 519.2, 519.6
ББК 22.171, 22.176, 22.19
Б01: 10.25728/^.2021.94.2
Статья представлена к публикации членом редакционной коллегии Я.И. Квинто.
Поступила в редакцию 15.11.2021. Дата опубликования 30.11.2021.
