Обчислення характеристик немарківських систем масового обслуговування чисельними методами

Нехай потік задається послідовністю випадкових величин, які визначають моменти надходження заявок. Легко помітити, що потік можна визначити, якщо задавати послідовність інтервалів між моментами надходження сусідніх заявок.

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

Потік з обмеженою післядією, для якого при , називається рекурентним потоком із запізненням або процесом відновлення із запізненням. Потік з обмеженою післядією, для якого при , називається рекурентним або процесом відновлення. В подальшому, розглядаючи процеси відновлення, будемо мати на увазі, що при .

Для простого процесу відновлення, який характеризується функцією розподілу F(x), позначимо через число відновлень (заявок, вимог), які відбулися до моменту t, а через H(t) - функцію відновлення, яка визначається рівністю

. (1.2.4.4.1)

Позначимо через момент k-го відновлення (). Тоді, за визначенням процесу відновлення, , - послідовність незалежних в сукупності випадкових величин із законом розподілу . Виразимо ймовірність через ймовірність F(x).

Помітимо, що має місце рівність , яка випливає із визначення випадкових величин та . Надалі, подію представимо у вигляді суми двох несумісних подій: . Тоді

. (1.2.4.4.2)

Підставляючи рівність (1.2.4.4.2) в співвідношення (1.2.4.4.1), отримуємо

, (1.2.4.4.3)

де - n-кратна згортка закону F(t), при чому

. (1.2.4.4.4)

Розподіл випадкової величини можна представити у вигляді n-кратної згортки закону F(t), так як - незалежні в сукупності випадкові величини і .

Тепер, використовуючи рівність (1.2.4.4.3), введемо інтегральне рівняння, яке пов'язує функцію відновлення H(t) та функцію розподілу F(x) випадкових величин .

Із визначення згортки двох функцій розподілу випливає рекурентне співвідношення

. (1.2.4.4.5)

Використовуючи рівності (1.2.4.4.3) та (1.2.4.4.4), запишемо

.

Але із рівності (1.2.4.4.3) маємо

.

Випливає шукане інтегральне рівняння відновлення

. (1.2.4.4.6)

Якщо останній інтеграл взяти по частинам і враховувати рівності , то інтегральне рівняння відновлення можна записати в іншому вигляді:

. (1.2.4.4.7)

Аналогічні рівняння можна отримати для функції відновлення процесу відновлення із запізненням. При введені співвідношення (1.2.4.4.2) не була використана та обставина, що випадкова величина і ті що зосталися , розподілені за одним і тим самим законом. Тому для процесу відновлення із запізненням залишається справедливою рівність (1.2.4.4.2). Тоді

(1.2.4.4.8)

Так як та - незалежні в сукупності , розподілені за законом при та при , то розподіл їх суми можна представити у вигляді згортки закону і (n-1)-го закону , звідси при або

, (1.2.4.4.9)

де задається рекурентним співвідношенням (1.2.4.4.5) при , або

, (1.2.4.4.10)

де визначається рекурентним співвідношенням (1.2.4.4.5) при .

Використовуючи рівності (1.2.4.4.8) та (1.2.4.4.9), запишемо

,(1.2.4.4.11)

, (1.2.4.4.12)

де H(t) - функція відновлення простого процесу відновлення, який характеризується функцією .

Використовуючи рівності (1.2.4.4.8) та (1.2.4.4.10), отримуємо

(1.2.4.4.13)

або

. (1.2.4.4.14)

Рівняння відновлення можуть бути розв'язані в термінах перетворення Лапласа-Стілтьєса.

Позначимо через перетворення Лапласа-Стілтьєса функції H(t). Тоді, використовуючи властивість цього перетворення, яке відноситься до інтегралу згортки, запишемо

або

. (1.2.4.4.15)

Для процесу відновлення із запізненням

. (1.2.4.4.16)

Для будь-якого моменту визначимо функцію

, (1.2.4.4.17)

якщо ця похідна існує. Функція називається щільністю відновлення. Враховуючи рівність (1.2.4.4.3) та визначення (1.2.4.4.17), можна величину трактувати як імовірність появи моменту відновлення (неважливо якого за підрахунком) в інтервалі . Неважко також отримати рівняння відновлення для щільності відновлення і розв'язку цих рівнянь в термінах перетворень Лапласа-Стілтьєса.

Сформулюємо елементарну теорему відновлення.

Теорема 1.2.4.4.1: Нехай H(t) - функція відновлення простого процесу відновлення, який визначається розподілом F(t); функція відновлення процесу відновлення із запізненням, який визначається розподілом та . Тоді та існують і

. (1.2.4.4.18)

Також наведу узлову теорему відновлення, яка часто використовується в різних задачах масового обслуговування та теорії надійності.

Теорема 1.2.4.4.2: Нехай функція Q(t)є невід'ємною функцією, не зростаюча та інтегрована в границях . Тоді при

.

Розділ 2 Наближене обчислення характеристик СМО M|G|1 з необмеженою чергою

2.1 Обчислення ймовірностей станів системи. Метод Філона

На рисунку зображена одна з реалізацій процесу - процес немарківський, але в моменти закінчення обслуговування заявок на приладі «пам'ять» про минуле втрачається, це так звані марківські моменти, значення процесу в ці моменти утворюють ланцюг Маркова, , який і називається вкладеним, де - момент закінчення обслуговування n-ї заявки.

рис.1

Для системи M|G|1 з необмеженою чергою Хінчіним і Поллячеком була знайдена твірна функція для ймовірностей станів процесу

(2.1.1).

Нехай , де - число заявок в системі в момент часу t, тобто - граничні ймовірності станів процесу обслуговування.

Формула Поллячека-Хінчіна має вигляд:

(2.1.2).

В цій формулі , (2.1.3).

Ці фігурують в матриці ймовірностей переходу за один крок вкладеного ланцюга Маркова в процес обслуговування .

На рисунку зображена одна з реалізацій процесу - процес немарківський, але в моменти закінчення обслуговування заявок на приладі «пам'ять» про минуле втрачається, це так звані марківські моменти, значення процесу в ці моменти утворюють ланцюг Маркова, який і називається вкладеним.

Матриця ймовірностей переходу за один крок для вкладеного ланцюга Маркова наведена в 1.2.1.10.

, що мають вигляд (2.1.3) - це відомі величини, вони виражаються через відому інтенсивність вхідного потоку л і відому функцію розподілу часу обслуговування однієї заявки одним приладом.

Оскільки з однієї сторони - це сума ряду в комплексній області, а ліва частина в формулі Поллячека-Хінчіна - аналітична функція, то на основі теореми Лорана можна записати

(2.1.4).

На основі цієї формули я будую обчислювальну процедуру для знаходження .

Спочатку криволінійний інтеграл з допомогою підстановки перетворюю на звичайний інтеграл, який можна обчислювати за допомогою відомих квадратурних формул (формули прямокутника, трапецій, Сімпсона). Нажаль, ці формули практично не застосовні для великих номерів к, оскільки підінтегральна функція у визначеному інтегралі містить у вигляді множників швидкоосцилюючі функції типу sinkx, coskx.

А саме:

Підставляючи (2.1.3):

,

де скориставшись розподілом Ерланга k-го порядку, отримала, що:

.

Розглянемо окремо чисельник:

Позначимо

і

.

Тоді отримаємо:



Також ввівши умовні позначення

.

Таким же чином введемо умовні позначення і для знаменника:

.

Тоді наш інтеграл матиме вигляд:

.

Необхідно позбавитися від ірраціональності в знаменнику дробу, для цього ми чисельник і знаменник помножимо на вираз, спряжений до знаменника, і розкривши дужки отримаємо:

Другий інтеграл буде рівним нулю, тому залишається, що

.

Далі перемножимо многочлени та винесемо та,тоді інтеграл матиме вигляд:

,

де

Для такого випадку добре пристосований метод Філона.

Метод Філона для наближеного обчислення інтегралів від тригонометричних функцій.

Наближене обчислення інтегралу вигляду , коли к не є малою величиною, пов'язано із значними труднощами, так як звичайні квадратурні формули, наприклад Сімпсона, потребують в наслідок швидкого коливання функції ділення проміжку інтегрування на дуже маленькі частини.

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

Розділимо проміжок інтегрування на 2n рівних частин з інтервалом h (або на n частин з інтервалом 2h), так, що

. (2.1.5)

Для зручності позначимо:

(2.1.6)

(2.1.7)

(2.1.8)

де s - ціле число.

Нехай в інтервалі , іншими словами, в можна з достатньою точністю функцію f(p) замінити параболічною дугою

(2.1.9)

Тоді знайдемо:

(2.1.10)

Продиференцюємо (2.1.9) за р, підставляючи замість В і С знайдені значення (2.1.10) поклавши потім , отримуємо:

(2.1.11)

Розглянемо інтеграл

.

Інтегрування його по частинам дає:

Так як згідно з (2.1.9) , то ми можемо написати:

(2.1.12)

Підставляючи сюди замість С його значення (2.1.10), користуючись позначенням (2.1.6) та маючи на увазі, що

знаходимо:

(2.1.13)

Користуючись надалі (2.1.11) та перетворюючи щойно вказаним способом , отримуємо:

(2.1.14)

Додаючи (2.1.13) та (2.1.14) на підставі (2.1.12), отримуємо

(2.1.15)

Оскільки , то коефіцієнт при в цьому виразі дорівнює

Аналогічно, так як , то коефіцієнт при в (2.1.15) дорівнює:

.

Таким чином, якщо мі покладемо:

(2.1.16)

І помітивши, що , то рівняння (2.1.15) можна написати в наступному вигляді:

Просумуємо тепер для , ми отримаємо формулу:

(2.1.17)

де - це сума усіх парних ординат кривої , яка знаходиться між а та в, окрім половини першої й останньої ординат; - це сума всіх непарних ординат, а величини б,в,г знаходяться з формули (2.1.16). Формула (2.1.17) заміняє формулу Сімпсона для розглянутого інтегралу, вона працює, навіть коли величина к велика, за умови, що інтеграл підібраний таким чином, що f(p) апроксимується з достатньою точністю параболічною дугою.

Для інтегралу, який має також виведена формула:

(2.1.18)

де - сума всіх парних ординат кривої , яка знаходиться між а та в, окрім половини першої й останньої ординат; - це сума всіх непарних ординат, а величини б,в,г мають попередній зміст.

Приклад:

Нехай , , тоді

Для інтегралу (2.1.17):

> I[1]=h*(A*(f(b)*sin(k*b)-f(a)*sin(k*a))+B*C[2*s]+Y*C[2*s-1]);

Для інтегралу (2.1.18):

> I[2]=h*((-A)*(f(b)*cos(k*b)-f(a)*cos(k*a))+B*S[2*s]+Y*S[2*s-1]);

Програмна реалізація цього прикладу наведена в додатку В.

2.2 Метод обернення перетворення Лапласа за допомогою ряду Фур'є за синусами. Приклад

Цей метод заснований на двох припущеннях. По-перше припускається, що зображення існує для . Це завжди можна зробити, якщо розглядати замість F(p) зображення F(p+a) для достатньо великих а. Останнє рівносильне множенню оригінала на . По-друге, для нашої задачі відоме зображення і за ним мі відновлюємо оригінал.

Інтеграл Лапласа перетворюють за допомогою підстановки , де д - довільне додатне число.

Тоді де і

. (2.2.1)

Функцію розкладають в ряд Фур'є за синусами непарних кратних дуг: , коефіцієнти якого знаходяться способом:

.

Для того, щоб виразити значення через значення зображення F(p) , зробимо наступним чином. В інтегралі в правій частині (2.2.1) покласти і урахувати, що

,

то після ряду перетворень одержують лінійну систему рівнянь з трикутною матрицею для визначення коефіцієнтів :

, (2.2.2)

де N - кількість членів ряду Фур'є , якими будемо апроксимувати функцію . Вибір значення д обумовлюється величиною проміжку, для якого необхідно обчислити значення оригіналу (д треба вибирати малим для великих t і, навпаки, великим для малих t).

Приклад:

Нехай , . Знаючи перетворення відновимо оригінал.

Отримуємо:

> c[0]:=12*2/(Pi*3^3);

> c[1]:=48*2/(9^3*Pi)-c[0];

> c[2]:=192*2/(Pi*15^3)-2*c[0]-3*c[1];

> c[3]:=768*2/(Pi*21^3)-5*c[0]-9*c[1]-5*c[2];

> c[4]:=3072*2/(Pi*27^3)-14*c[0]-28*c[1]-20*c[2]-7*c[3];

> c[5]:=12288*2/(Pi*33^3)-42*c[0]-90*c[1]-75*c[2]-35*c[3]-9*c[4];

Sum(c[k]*sin((2*k+1)*27.99845),k=0..5)=sum(c[k]*sin((2*k+1)*27.99845),k=0..5);

Програмна реалізація цього прикладу наведена в додатку В.

2.3 Алгоритм обчислення характеристик M|G|1 з необмеженою чергою методом статистичного моделювання.

Розглянемо систему M|G|1 з необмеженою чергою.

На рисунку 2 представлений алгоритм моделювання роботи СМО методом Монте-Карло:



рис. 2

Де на 1-й осі моделюємо моменти надходження заявок в систему масового обслуговування.

На 2-й осі - обслуговування заявок, 3-а вісь - результуюча, на ній зображено, яка кількість заявок знаходиться в СМО в кожен момент часу.

Нехай в систему надходить простіший потік заявок з інтенсивністю б. Щільності вважаємо відомими.

Поставимо у відповідність наш прилад обслуговування в комірку ЕОМ, в яку будемо записувати момент звільнення приладу . В початковий момент, очевидно, .

Момент надходження першої заявки розіграємо за формулою:

.

Час обслуговування ф знайдемо з рівняння:

.

Прилад буде зайнятий до моменту . Тому у відповідній комірці потрібно замінити на і обчислити нове .

Нехай тепер долю (k-1)-ї заявки, яка надійшла в момент , ми вже з'ясували. Тоді розіграємо момент надходження k-ї заявки:

та перевіряємо, чи вільний прилад. Для цього достатньо перевірити нерівність

а) Якщо нерівність виконується, то прилад, для якого , починає обслуговування k-ї. Тривалість обслуговування ф розігрується за формулою:

;

Потім обчислюється і нове значення . Після цього можна перейти до розгляду наступної (k+1)-ї заявки.

б) Якщо нерівність не виконується, то прилад в момент зайнятий і заявка стає в чергу. Кількість місць в черзі не обмежена, тому можна перейти до розгляду (k+1)-ї заявки.

Розрахунок продовжується до тих пір, поки момент надходження заявок не стане більше, ніж t_скін.

І повторимо такий дослід N разів (з різними числами).

Таким чином мі розіграли моменти надходження заявок в систему, час обслуговування заявок приладом.

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

Розділ 3 Описання програмного забезпечення

Що таке Maple? Одним зі світових лідерів в комп'ютеризації математичних обчислень (у тому числі символьних) являється корпорація Waterloo Maple Inc. (Канада), що випускає програмний продукт Maple. Остання версія охоплює майже усю математику, починаючи з елементарної математики і закінчуючи спеціальними математичними розділами. Maple - математичний windows-додаток, який дозволяє розв'язувати задачі з цього широкого діапазону за мінімальний час.

Чи складно працювати в Maple? Виявляється, що не дуже. Програми рішень основних математичних задач і геометричних побудов, які складені авторами Maple, надаються користувачеві тільки іменами, в круглих дужках після яких вводяться необхідні дані. Враховуючи залежність від даних, їх називають вбудованими функціями. Задача користувача - вибудовувати з них і операторів потрібні послідовності і задавати дані. Втім, часто виявляється, що досить скористатися однією із вбудованих функцій, тим більше, що їх в Maple більше 3000.

Чи можна довіряти отриманим в Maple результатам? Так, можна. Наведені далі приклади зайвий раз підтверджують це. Більше того, Maple - перша комп'ютерна програма, що пройшла тестування з результатом 100 %. Комп'ютерна програма Maple незамінна як для перевірки остаточних і проміжних результатів, що отримують аналітично - без комп'ютера, так і для пошуку методів розв'язків. Сам пакет широко поширений в університетах провідних наукових держав, дослідницьких центрах і компаніях. При цьому пакет Maple розвивається, вбираючи в себе нові уміння, нові розділи математики і забезпечуючи краще середовище для роботи.

Maple складається з ядра - процедур, що написані на мові С і вкрай оптимізованих, бібліотеки, написаної на Maple -мові, та інтерфейсу. Ядро виконує більшість базисних операцій. Бібліотека містить безліч команд - процедур, що виконуються в режимі інтерпретації. Програмуючи власні процедури, користувач може поповнювати ними стандартний набір і, таким чином, розширювати можливості Maple.

Алгебраїчні перетворення:

Вбудовані функції елементарних перетворень:

simplify - спростити,

expand - розкрити дужки,

factor - розкласти на множники,

normal - привести до спільного знаменника,

combine - перетворити степені (чи тригонометричні вирази),

collect - привести подібні члени.

Після ключового слова в дужках вводиться аналітичний вираз або
його ім'я - ідентифікатор, а також параметри, частина яких або усі можуть бути відсутні - бути необов'язковими. Наприклад, застосовуючи collect, щоб не було повідомлення про помилку, необхідно вказати змінну, по степенях якої приводяться подібні члени - обов'язковий параметр. У simplify може
бути додана вбудована функція assume - прийняти, яка задає умови, при
яких відбувається спрощення, - необов'язковий параметр. Вбудована функція combine також має необов'язкові параметри.

> simplify((a^3-b^3)/(a-b));

> expand((a-b)*(a^2+a*b+b^2));

> factor(a^3-b^3);

> normal(y/x+1/(x^2));

> collect(x^2+3*x^2+4*x+4*x+y,x);

> simplify(2*a/sqrt(a^2),assume(a<0));

> combine(x^(1/2)*x^(3/2));

Оператори циклу:

Часто буває необхідним циклічне повторення виконання виразів задане число разів, або до тих пір, поки виконується задана умова.

В Maple є декілька операторів циклу, і для їх запису використовуються службові слова for, from, by, to, while, do, end do. Тілом всіх операторів циклу є послідовність команд, які зосереджені між do та end do. Почнемо з перелічного циклу, який є практично в усіх алгоритмічних мовах:

for VAR from VAL1 by VAL2 to VAL3 do EXPR end do;

Тіло циклу EXPR виконується при кожному значенні параметру VAR , який змінюється від VAL1 з кроком VAL2 до тих пір, поки не стане більше VAL3. Відмічу, що якщо крок зміни VAL2 дорівнює одиниці, то оператор циклу допускає скорочену форму:

for VAR from VAL1 to VAL3 do EXPR end do;

> for h from 1 by 2 to 20 do h end do;

Математичний аналіз:

В Maple для деяких математичних перетворень існує по дві команди: одна прямого, а друга - відкладеного значення, при чому назви цих команд складаються з однакових літер. Команди прямого значення, як правило, починаються з маленької літери і виконується негайно. Відкладені команди часто починаються з великої літери, за звичай в тому випадку, коли існує команда-синонім прямої дії.

Для операцій сумування виконуються команди

Sum(expr, var=var1...var2); та sum(expr, var=var1...var2);

Де expr - вираз, який залежить від змінної сумування var, а var1...var2 - границі сумування. Границі сумування можуть бути скінченими, нескінченими або взагалі відсутні. Така команда також може використовуватись для сумування рядів.

> Sum(3*(k-1),k=1..n-1)=sum(3*(k-1),k=1..n-1);

> Sum(3*(k-1),k=1..n)+Sum((8+3*k)/2,k=0..n)=sum(3*(k-1),k=1..n)+sum((8+3*k)/2,k=0..n);

Для операцій інтегрування виконуються команди

Для такої мети в пакеті передбачено декілька команд, які знаходяться в різних бібліотеках. В стандартній бібліотеці знаходяться процедури int(expr,par) та Іnt(expr,par), які в залежності від параметрів par можуть використовуватись для знаходження невизначених інтегралів, аналітичного чи числового обчислення визначених, власних і не власних інтегралів. Для обчислення визначених інтегралів виконується команда

іnt(expr,var=a..b);

Де expr - інтегрований вираз, var - змінна інтегрування, а (a;b) - проміжок інтегрування, його кінці можуть приймати значення нескінченності. Для обчислення подвійних, потрійних і т. д. інтегралів необхідно застосовувати цю команду декілька разів.

> Int(6-x-5/x,x=1..5)=int(6-x-5/x,x=1..5);

Висновки

Для деяких класичних систем масового обслуговування теоретичні результати для обчислення характеристик системи через вхідні параметри (інтенсивність вхідного потоку, інтенсивність обслуговування, число каналів) задаються аналітично, в основному це системи, процеси обслуговування в яких описуються марківськими процесами. Для немарківських СМО результати для обчислення характеристик, як правило, знаходяться в термінах твірних функцій, в термінах перетворень Лапласа та інших перетворень.

Для практичного використання результатів необхідно вміти знаходити за відомими перетвореннями потрібні характеристики.

Я в дипломній роботі на прикладі СМО M|G|1 з необмеженою чергою, показала деякі підходи до наближеного обчислення обернених перетворень.

І. Для системи M|G|1 з необмеженою чергою Хінчіним і Поллячеком була знайдена твірна функція для ймовірностей станів процесу

Нехай (1), де - число заявок в системі в момент часу t, тобто - граничні ймовірності станів процесу обслуговування.

Формула Поллячека-Хінчіна має вигляд:

.

В цій формулі , (4) .

Ці фігурують в матриці ймовірностей переходу за один крок вкладеного ланцюга Маркова в процес обслуговування .

На рисунку (1) зображена одна з реалізацій процесу - процес не Марківський, але в моменти закінчення обслуговування заявок на приладі «пам'ять» про минуле втрачається, це так звані марківські моменти, значення процесу в ці моменти утворюють ланцюг Маркова, який і називається вкладеним.

- момент закінчення обслуговування n-ї заявки,

- вкладений ланцюг Маркова.

Матриця ймовірностей переходу за один крок для вкладеного ланцюга Маркова має вигляд (1.2.1.10)

- це відомі величини, вони виражаються через відому інтенсивність вхідного потоку л і відому функцію розподілу часу обслуговування однієї заявки одним приладом.

Оскільки з однієї сторони - це сума ряду в комплексній області, а ліва частина в формулі Поллячека-Хінчіна - аналітична функція, то на основі теореми Лорана можна записати

.

На основі цієї формули я побудувала обчислювальну процедуру для знаходження .

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

Нажаль, ці формули практично не застосовні для великих номерів к, оскільки підінтегральна функція у визначеному інтегралі містить у вигляді множників швидкоосцилюючі функції типу sinkx, coskx. Для такого випадку добре пристосований метод Філона.

ІІ. Для деяких характеристик СМО результати отримані в термінах перетворень Лапласа або Лапласа-Стілтьєса.

рис.3

Наприклад, для СМО M|G|1 знайдено перетворення Лапласа-Стілтьєса функції розподілу періоду зайнятості системи. Під періодом зайнятості розуміють випадкову величину Х, яка має наступний зміст, що це час безперервної зайнятості обслуговуючого приладу. Ця характеристика відіграє важливу роль при дослідженні питання, чи існує в СМО стаціонарний процес, чи ні? Якщо пропускна здатність СМО буде не достатньо високою, то величина черги з плином часу буде необмежено збільшуватись, що в реальних СМО недопустимо.

Нехай - функція розподілу періоду зайнятості,

- перетворення Лапласа-Стілтьєса функції розподілу періоду зайнятості.

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

Інтеграл Лапласа перетворюють за допомогою підстановки , де д - довільне додатне число. Тоді

і

.

Функцію розкладають в ряд Фур'є: для наближеного обчислення беруть скінчену кількість доданків.

Після ряду перетворень отримуємо, що , де с_к - знаходяться із рекурентних співвідношень

.

ІІІ. Ці ж самі характеристики можна оцінити, застосовуючи чисельні алгоритми методу Монте-Карло.

На рисунку (2) представлений алгоритм моделювання роботи СМО методом Монте-Карло:

Де на 1-й осі моделюємо моменти надходження заявок в СМО.

На 2-й осі - обслуговування заявок, 3-а вісь - результуюча, на ній зображено, яка кількість заявок знаходиться в СМО в кожен момент часу. Для того щоб оцінити, наприклад, ймовірність треба розіграти одну досить довгу траєкторію процесу і по ній знайти суму довжин відрізків де число заявок в системі дорівнює нулю, а потім це число розділити на Т. Аналогічно оцінюємо і інші імовірності.

Таким чином я дослідила наближене обчислення характеристик для системи M|G|1 з необмеженою чергою, склала програмну реалізацію для методу перетворення Лапласа та методу Філона, завдяки чому полегшується обчислювальний апарат.

Список літератури

1. Г.І. Івченко, В.А. Каштанов, И.Н. Коваленко «Теорія масового обслуговування»

2. В.І.Крилов, Н.С. Скобля «Методи наближеного перетворення Фур'є і обернення перетворення Лапласа»

3. М.Л. Краснов, А.І. Кисилев, Г.И. Макаренко «Функції комплексної змінної. Операційне числення. Теорія стійкості»

4. К. Дж. Трнатер «Інтегральні перетворення в математичній фізиці»

5. О.А. Сдвіжков «Математика на комп'ютері: MAPLE»

6. В.С. Королюк, Турбін А.Ф. «Ніпівмарківські процеси і їх додатки»

7. А. Кофман, Р. Крюон «Масове обслуговування»

8. Т.Сааті «Елементи теорії масового обслуговування та її додатки»

9. О.Я. Хінчин «Праці за математичної теорії масового обслуговування»

10. С.М. Броді, І.А. Погосян «Вкладені стохастичні процеси в теорії масового обслуговування»

11. В.С. Королюк, О.Ф. Турбін «Напівмарківські процеси та їх додатки»

12. Дж. Кендалл «Стохастичні процеси, які зустрічаються в теорії черг, та їх аналіз методом вкладених ланцюгів маркова»

13. В.Л.Смит «Теорія відновлення і суміжні з нею питання»

14. О. С. Вентцель «Дослідження операцій»

Додаток А

Обчислення характеристик немарківських СМО чисельними методами

І. Наближене обчислення імовірнісних розподілів за заданою твірною функцією.

СМО M|G|1 з необмеженою чергою:

М - (markoven) інтервали між надходженнями заявок розподілені за показниковим законом,

G - (general) час обслуговування заявок розподілений за довільним законом,

1 - система має один пристрій для обслуговування

- число заявок в системі в момент часу ,

(1)

(2)

(3)

- формула Поллячека-Хінчіна.

, (4)

 , (5)

- функція розподілу часу обслуговування однієї заявки.

рис.1

- момент закінчення обслуговування n-ої заявки,

- вкладений ланцюг Маркова.

Матриця ймовірностей переходу за один крок для вкладеного ланцюга Маркова має вигляд:

(6)

 (7)

,

.(8)

Метод Філона для наближеного обчислення інтегралів

(9)

Розділимо проміжок інтегрування на 2n рівних частин з інтервалом h:

(10)



ІІ. Наближене обчислення характеристик СМО за допомогою оберненого перетворення Лапласа.

рис.3

(14)

- функція розподілу періоду зайнятості,

 (15)

- перетворення Лапласа-Стілтьєса функції розподілу періоду зайнятості.

Інтеграл Лапласа

(16)

,

де д - довільне додатне число.

,(17)

. (18)

Розвиток в ряд Фур'є:

 (19)

- знаходяться із рекурентних співвідношень:

. (20)

Наближене обчислення обернення перетворення Лапласа:

, (21)

Додаток В:

Програмна реалізація для 2.2:

Інтеграл Лапласа:

> F(p):=Int(exp^(-pt)*f,t=0..infinity);

Перетворюють за допомогою підстановки е^wt=cosv, де w - довільне додатне число, тоді:

> f(t):=f(-1/w*ln(cos(v)))=y(v);

> w*F(p)=Int(((cos(v))^(p/w-1))*(sin(v)*y(v)),v=0..Pi/2);

Функцію розкладають в ряд Фур'є за синусами непарних кратних дуг:

> y(v)=sum(c[k]*sin((2*k+1)*v),k=0..infinity);

> sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..n):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> restart;

Проведемо розрахунки для w=3, t=0,1:

> k:=0: n:=0: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=0: n:=1: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=1: n:=1: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=0: n:=2: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=1: n:=2: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=2: n:=2: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);



> k:=0: n:=3: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=1: n:=3: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=2: n:=3: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=3: n:=3: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=0: n:=4: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=1: n:=4: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);



> k:=2: n:=4: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=3: n:=4: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=4: n:=4: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=0: n:=5: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=1: n:=5: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=2: n:=5: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);



> k:=3: n:=5: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=4: n:=5: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> k:=5: n:=5: w:=3:

>sum(((2*k+1)/(2*n+1))*(factorial(2*n+1)/(factorial(n-k)*factorial(n+k+1)))*c[k],k=0..5):=(4^(n+1)/Pi)*w*F((2*n+1)*w);

> restart;

> c[0]:=12*2/(Pi*3^3);

> c[1]:=48*2/(9^3*Pi)-c[0];

> c[2]:=192*2/(Pi*15^3)-2*c[0]-3*c[1];

> c[3]:=768*2/(Pi*21^3)-5*c[0]-9*c[1]-5*c[2];

> c[4]:=3072*2/(Pi*27^3)-14*c[0]-28*c[1]-20*c[2]-7*c[3];

> c[5]:=12288*2/(Pi*33^3)-42*c[0]-90*c[1]-75*c[2]-35*c[3]-9*c[4];

>Sum(c[k]*sin((2*k+1)*27.99845),k=0..5)=sum(c[k]*sin((2*k+1)*27.99845),k=0..5);

Програмна реалізація для 2.1:

Введіть необхідні данні:

> restart;

> a:=0;

> b:=2*3.14159;

> n:=1000;

> h:=(b-a)/(2*n);

> k:=10;

> Q:=k*h;

> f(p):=1;

> f(b):=1;

> f(a):=1;

> for j from 1 to 2*n do f(j*h):=1 end do:

Підрахунок:

> A:=(Q^2+Q*sin(Q)*cos(Q)-2*sin(Q)^2)/Q^3;

> B:=2*(Q*(1+cos(Q)^2)-2*sin(Q)*cos(Q))/Q^3;

> Y:=4*(sin(Q)-Q*cos(Q))/Q^3;

Для інтегралу з coskx:

>

> for j from 1 to n do C[j]:=f(2*j*h)*cos(k*(2*j)*h) end do:

>

> C[2*s]:=sum(C[i],i=1..n)+0.5*(f(a)*cos(k*a)+f(b)*cos(k*b));

>

>

> for j from 0 to n-1 do C[j]:=f((2*j+1)*h)*cos(k*(2*j+1)*h) end do:

>

> C[2*s-1]:=sum(C[i],i=0..n-1);

> I[1]=h*(A*(f(b)*sin(k*b)-f(a)*sin(k*a))+B*C[2*s]+Y*C[2*s-1]);

Для інтегралу з sinkx:

> for j from 1 to n do S[j]:=f(2*j*h)*sin(k*(2*j)*h) end do:

> S[2*s]:=sum(C[i],i=1..n)+0.5*(f(a)*sin(k*a)+f(b)*cos(k*b));

> for j from 0 to n-1 do S[j]:=f((2*j+1)*h)*sin(k*(2*j+1)*h) end do:

> S[2*s-1]:=sum(C[i],i=0..n-1);

> I[2]=h*((-A)*(f(b)*cos(k*b)-f(a)*cos(k*a))+B*S[2*s]+Y*S[2*s-1]);