Моделирование систем

Моделирование многоканальных систем массового обслуживания

Разбить на страницы
Показывать лекцию целиком

Теоретическая часть

Многоканальные системы массового обслуживания — это системы с параллельно включенными приборами обслуживания. Для них принято использовать символику Кендалла, которая состоит из основных четырех позиций вида $$A/B/m/K/М$$, где $$А$$ — закон поступления требований в систему, $$В$$ — закон обслуживания требований, m — число параллельно функционирующих приборов (каналов, узлов) обслуживания, $$K$$ — допустимое число требований в системе, т. е. число требований в очереди плюс число требований, принятых на обслуживание. В системе $$A/B/m/K/М$$ последний символ $$М$$ — число источников нагрузки [2]. При этом параметры системы постоянны, а процессы в системах рассматриваются с позиций теории размножения и гибели.

Предметом данной лабораторной работы будут системы с простейшим пуассоновским потоком, экспоненциальным обслуживанием и, возможно, экспоненциальном уходом из очереди "нетерпеливых" требований. При сделанных условиях исследуемые системы будут иметь следующие обозначения: $$M/M/m/K$$, $$M/M/m/K$$ с нетерпеливыми требованиями, $$M/M/m/K/M$$, где первые буквы $$М$$ означают пуассоновский поток требований и экспоненциальное обслуживание, а в последнем примере крайняя правая буква $$М$$ — это число требований, формируемых конечным числом источников нагрузки. В частности, может быть система с отказами, т. е. $$M/M/m$$ с отказами.

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

Для составления уравнений Колмогорова можно использовать мнемоническое правило [4]. Сначала определим понятие потока вероятности: это такой поток, который переводит систему из одного состояния в другое соседнее и определяется как произведение вероятности $$P_{i}(t)$$ $$i$$ -го состояния, из которого происходит переход, на интенсивность потока событий (интенсивность поступления требований или интенсивность обслуживания). Теперь приведем мнемоническое правило составления уравнений Колмогорова: производная вероятности любого состояния равна сумме потоков вероятности, переводящих систему в это состояние, минус сумма всех потоков вероятности, выводящих систему из этого состояния [4].

Прежде чем применять правило Колмогорова, целесообразно изобразить размеченный граф состояний заданной системы. На рис. 2.1 показан пример размеченного графа состояний для системы M/M/3/5 ( $$m = 3$$, $$K = 5$$, $$\lambda$$ — интенсивность входного потока требований, $$\mu$$ — интенсивность обслуживания одним прибором).

(рис 2.1) Размеченный граф состояний системы М/М/3/5

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

$$\frac{dP_{0}}{dt}=-\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=\lambda P_{0}-(\lambda +\mu) P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=\lambda P_{1}-(\lambda +2\mu) P_{2}+3\mu P_{3};\\ \frac{dP_{3}}{dt}=\lambda P_{2}-(\lambda +3\mu) P_{3}+3\mu P_{4};\\ \frac{dP_{4}}{dt}=\lambda P_{3}-(\lambda +3\mu) P_{4}+3\mu P_{5};\\ \frac{dP_{5}}{dt}=\lambda P_{4}-3\mu P_{5};\\$$

Для решения систем дифференциальных уравнений типа (2.1) обычно задают естественные начальные условия:

$$P_{0}(0)=1,\qquad P_{i}(0)=0,\qquad i=\overline{1,5}.$$

Если в системе (2.1) производные приравнять нулю, то можно будет получить соотношения для расчета стационарных вероятностей состояний системы. При этом следует использовать нормировочное условие для $$n = 5$$:

$$\sum\limits_{k=0}^{5} P_{k}(t)=1.$$

Практическая часть

2.1. Пример моделирования системы типа М/М/M/K

Система $$M/M/m/K$$ — это система с пуассоновским входящим потоком требований, с экспоненциальным законом обслуживания в m приборах, с допустимым числом требований в системе, не превышающим $$K$$, которое не менее, чем заданное количество приборов обслуживания. Параметры системы постоянны, т. е. $$\lambda = const$$, $$\mu = const$$. Систему $$M/M/m/K$$ называют многоканальной системой массового обслуживания с ограниченной длиной очереди.

Пример 1. Проинтегрируйте систему массового обслуживания $$М/М/4/7$$ при естественных граничных условиях и параметрах $$\lambda = 1.23$$, $$\mu = 0.678$$. Рассчитайте операционные характеристики системы [2].

Для решения примера сначала составим дифференциальные уравнения в соответствии с мнемоническим правилом Колмогорова:

$$\frac{dP_{0}}{dt}=-\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=\lambda P_{0}-(\lambda +\mu) P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=\lambda P_{1}-(\lambda +2\mu) P_{2}+3\mu P_{3};\\ \frac{dP_{3}}{dt}=\lambda P_{2}-(\lambda +3\mu) P_{3}+4\mu P_{4};\\ \frac{dP_{4}}{dt}=\lambda P_{3}-(\lambda +4\mu) P_{4}+4\mu P_{5};\\ \frac{dP_{5}}{dt}=\lambda P_{4}-(\lambda +4\mu) P_{5}+4\mu P_{6};\\ \frac{dP_{6}}{dt}=\lambda P_{5}-(\lambda +4\mu) P_{6}+4\mu P_{7};\\ \frac{dP_{7}}{dt}=\lambda P_{6}-4\mu P_{7};\\$$

Естественные начальные условия:

$$P_{0}(0)=1,\qquad P_{i}(0)=0,\qquad i=\overline{1,7}.$$

Программный код решения примера в MATLAB:

function MMmK;
clc,close 
% Параметры системы
L = 3.52;
M = 0.678;
m = 4;
K = 7; 
%%%%%%%%%%%%%%%%%%%%%%%%
global A
A = [-L,M,0,0,0,0,0,0;
L,-(L+M),2*M,0,0,0,0,0;
0,L,-(L+2*M),3*M,0,0,0,0;
0,0,L,-(L+3*M),4*M,0,0,0;
0,0,0,L,-(L+4*M),4*M,0,0;
0,0,0,0,L,-(L+4*M),4*M,0;
0,0,0,0,0,L,-(L+4*M),4*M;
0,0,0,0,0,0,L,-4*M];
%% Численное интегрирование дифф. уравнений
P0 = [1;zeros(length(A)-1,1)];
T = [0,20];
[t,P] = ode23(@cmo, T, P0);
 
%% Построение диаграммы вероятностей состояний
%% line(t,P,'linew',2) %% с различными цветами
line(t,P(:,1),'linew',2, 'color','r') %% Po
line(t,P(:,2), 'linew',2,'lines','--') %% P1
line(t,P(:,3), 'linew',2,'lines','-.') %% P2
line(t,P(:,4), 'linew',2,'lines',':') %% P3
line(t,P(:,5), 'marker','o', 'color', 'm') %% P4
line(t,P(:,6), 'marker','h', 'color','k') %% P5
line(t,P(:,7), 'marker','p','color','r') %% P6
line(t,P(:,8), 'marker','>') %% P7
grid on
N = length(A)-1;
arr = [0:N]';
str = num2str(arr);
legend(strcat('\bf\itP\rm\bf_', str, '(\itt\rm\bf)'));
title(sprintf('%s Вероятности состояний системы M/M/%d/%d', '\bf\fontsize{12}',m, K));
xlabel('\bf\it\fontsize{12} -  -  -  -  -  -  -  -   t   -  -  -  -  -  -  -  -')
ylabel('\bf\fontsize{12}\itP\rm\bf(\itt\rm\bf)');
set(gca, 'fontweight','bold', 'fontsize',10)

fprintf('\n Стационарные вероятности:\n');
for J = 1 : length(A)
    fprintf('\tP%d = %f\n', J-1, P(end,J));
end
 fprintf('\n\t ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ:\n');
Pnot = P(end,end);

fprintf(' Вероятность отказа Pnot = %f\n', P(end,end));
Q = 1 - Pnot;
fprintf(' Относительная пропускная способность Q = %f\n', Q);
Ab = L*Q;
fprintf(' Абсолютная пропускная способность A = %f\n', Ab);
Pq = sum(P(end, m+1:end));
fprintf(' Вероятность наличия очереди Pq = %f\n', Pq);
Ps = sum(P(end, m:end));
fprintf(' Вероятность загрузки всех каналов обслуживания Ps = %f\n', Ps);
 
Ns = [0:length(A)-1]*P(end,:)';
fprintf(' Среднее количество требований в системе Ns = %f\n', Ns);
fprintf(' Среднее время пребывания требования в системе Ts = %f\n', Ns/L);
Nq = [0:(K-m)]*P(end,m:K)';
fprintf(' Средняя длина очереди Nq = %f\n', Nq);
fprintf(' Среднее время пребывания требования в очереди Tq = %f\n', Nq/L);
 
 
function f = cmo(t,P);
%% Функция описания правых частей 
%% дифференциальных уравнений
global A
f = A*P;

Результат выполнения программы

Стационарные вероятности:
	P0 = 0.004352
	P1 = 0.022604
	P2 = 0.058629
	P3 = 0.101614
	P4 = 0.131597
	P5 = 0.171220
	P6 = 0.221807
	P7 = 0.288176

	 ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ
 Вероятность отказа: Pnot = 0.288176
 Относительная пропускная способность: Q = 0.711824
 Абсолютная пропускная способность: A = 2.505620
 Вероятность наличия очереди: Pq = 0.812800
 Вероятность загрузки всех каналов обслуживания: Ps = 0.914415
 Среднее количество требований в системе: Ns = 5.175269
 Среднее время пребывания требования в системе: Ts = 1.470247
 Средняя длина очереди: Nq = 1.139458
 Среднее время пребывания требования в очереди: Tq = 0.323710

Среднее время пребывания требования в системе было рассчитано по формуле Литтла. Для расчета среднего времени пребывания требования в очереди использовался аналог формулы Литтла [2].

Диаграмма с вероятностями состояний системы $$М/М/4/7$$ показана на рис. 2.2.

(рис 2.2) Вероятности состояний системы М/М/4/7

Задание 1

  • Рассчитайте нормировочное условие для трех значений времени.
  • Определите вероятность того, что в системе будет не более $$k$$ требований. Значение $$k$$ выбирайте случайно по равномерному закону из интервала целых чисел $$[1; 6]$$ (см. $$help randperm$$ ).
  • Примените диалоговое окно $$inputdlg$$ (см. $$help inputdlg$$ ) для ввода параметров системы.
  • Напишите программу для произвольно задаваемых параметров системы $$\lambda$$, $$\mu$$ и значений $$m$$, $$K$$. Предусмотрите также ввод длительности интервала интегрирования.
  • 2.2. Пример моделирования системы типа М/М/M/K с ограниченным временем ожидания

    Пример 2. Получите операционные характеристики системы $$M/M/4/7$$ с ограниченным временем ожидания при следующих параметрах системы: $$\lambda = 3.52 \mbox{ }с^{–1}$$ — интенсивность входного потока требований, $$\mu = 0.678 \mbox{ }с^{–1}$$ — интенсивность обслуживания каждым прибором, $$\nu = 1.3\mbox{ }с^{–1}$$ — интенсивность ухода из очереди "нетерпеливых" требований. Моделирование проведите в MATLAB и GPSS/PC.

    Размеченный граф состояний системы показан на рис. 2.3.

    (рис 2.3) Граф состояний системы с ограниченным временем ожидания

    На рис. 2.3 вокруг состояния 4 проведена воображаемая окружность для пояснения применения мнемонического правила составления дифференциальных уравнений Колмогорова. В соответствии с мнемоническим правилом получаем следующую систему дифференциальных уравнений:

    $$\frac{dP_{0}}{dt}=-\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=\lambda P_{0}-(\lambda +\mu) P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=\lambda P_{1}-(\lambda +2\mu) P_{2}+3\mu P_{3};\\ \frac{dP_{3}}{dt}=\lambda P_{2}-(\lambda +3\mu) P_{3}+4\mu P_{4};\\ \frac{dP_{4}}{dt}=\lambda P_{3}-(\lambda +4\mu) P_{4}+(4\mu +\nu) P_{5};\\ \frac{dP_{5}}{dt}=\lambda P_{4}-(\lambda +4\mu+\nu) P_{5}+(4\mu +2\nu)P_{6};\\ \frac{dP_{6}}{dt}=\lambda P_{5}-(\lambda +4\mu+2\nu) P_{6}+(4\mu+3\nu) P_{7};\\ \frac{dP_{7}}{dt}=\lambda P_{6}-(4\mu+3\nu) P_{7};\\$$

    где $$P_{i},\mbox{ } i=\overline{0,7}$$ — вероятности состояний системы.

    Для интегрирования системы дифференциальных уравнений зададим естественные граничные условия:

    $$P_{0}(0)=1,\qquad P_{i}(0)=0,\qquad i=\overline{1,7}.$$

    Программный код решения примера в системе MATLAB:

    function MMmKv
    %% Система M/M/4/7 с ограниченным временем ожидания в очереди
    % Параметры системы
    L = 3.52;
    M = 0.678;
    v = 1.3;
    m = 4; 
    K = 7;
    
    % Матрица коэффициентов системы дифференциальных уравнений
    global A
    A = [-L, M, zeros(1,6);
        L, -(L+M), 2*M, zeros(1,5);
        0, L, -(L+2*M), 3*M, zeros(1,4);
        0, 0, L, -(L+3*M), 4*M, zeros(1,3);
        zeros(1,3), L, -(L+4*M), (4*M+v), 0, 0;
        zeros(1,4), L, -(L+4*M+v), (4*M+2*v), 0;
        zeros(1,5), L, -(L+4*M+2*v), (4*M+3*v);
        zeros(1,6), L ,-(4*M+3*v)];
    % Решение дифференциальных уравнений
    T = [0,20]; % Время интегрирования
    P0 = [1;zeros(length(A)-1,1)]; % Начальные условия
    [t, P] = ode23(@cmo, T, P0);
    
    %% Построение диаграмм вероятностей состояний
    % line(t,P,'linew',2) %% сполошные линии с различными цветами
    line(t,P(:,1),'linew',2, 'color','r') %% Po
    line(t,P(:,2), 'linew',2,'lines','--') %% P1
    line(t,P(:,3), 'linew',2,'lines','-.') %% P2
    line(t,P(:,4), 'linew',2,'lines',':') %% P3
    line(t,P(:,5), 'marker','o', 'color', 'm') %% P4
    line(t,P(:,6), 'marker','h', 'color','k') %% P5
    line(t,P(:,7), 'marker','p','color','r') %% P6
    line(t,P(:,8), 'marker','>') %% P7
    grid on
    N = length(A)-1;
    arr = [0:N]';
    str = num2str(arr);
    legend(strcat('\bf\itP\rm\bf_', str, '(\itt\rm\bf)'));
    title(sprintf('%s Вероятности состояний системы M/M/%d/%d', '\bf\fontsize{12}',m, K));
    xlabel('\bf\it\fontsize{12} -  -  -  -  -  -  -  -   t   -  -  -  -  -  -  -  -')
    ylabel('\bf\fontsize{12}\itP\rm\bf(\itt\rm\bf)');
    
    set(gca, 'fontweight','bold', 'fontsize',10)
    
    Pcm = P(end,:); % Стационарные вероятности
    fprintf('\n Стационарные вероятности:\n');
    for J = 1 : length(A)
        fprintf('\tP%d = %f\n', J-1, Pcm(J));
    end
    fprintf('\n\t ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ:\n');
    Pnot = P(end,end);
    fprintf(' Вероятность отказа Pnot = %f\n', P(end,end));
    Q = 1 - Pnot;
    fprintf(' Относительная пропускная способность Q = %f\n', Q);
    Ab = L*Q;
    fprintf(' Абсолютная пропускная способность A = %f\n', Ab);
    Pq = sum(P(end, m+1:end));
    fprintf(' Вероятность наличия очереди Pq = %f\n', Pq);
    Ps = sum(P(end, m:end));
    fprintf(' Вероятность загрузки всех каналов обслуживания Ps = %f\n', Ps);
    Ns = [0:length(A)-1]*P(end,:)';
    fprintf(' Среднее количество требований в системе Ns = %f\n', Ns);
    fprintf(' Среднее время пребывания требования в системе Ts = %f\n', Ns/L);
    Nq = [0:(K-m)]*P(end,m:K)';
    fprintf(' Средняя длина очереди Nq = %f\n', Nq);
    fprintf(' Среднее время пребывания требования в очереди Tq = %f\n', Nq/L);
    
    function f = cmo(t,P);
    %% Функция описания правых частей 
    %% дифференциальных уравнений
    global A
    f = A*P;

    Результат выполнения программы

    Стационарные вероятности:
    	P0 = 0.007887
    	P1 = 0.040947
    	P2 = 0.106289
    	P3 = 0.183963
    	P4 = 0.238701
    	P5 = 0.209568
    	P6 = 0.138707
    	P7 = 0.073939
    
    	 ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ
     Вероятность отказа: Pnot = 0.073939
     Относительная пропускная способность: Q = 0.926061
     Абсолютная пропускная способность: A = 3.259736
     Вероятность наличия очереди: Pq = 0.660915
     Вероятность загрузки всех каналов обслуживания: Ps = 0.844877
     Среднее количество требований в системе: Ns = 4.157868
     Среднее время пребывания требования в системе: Ts = 1.181213
     Средняя длина очереди: Nq = 1.073958
     Среднее время пребывания требования в очереди: Tq = 0.305102

    Диаграмма вероятностей состояний системы $$М/М/4/7$$ с ограниченным временем ожидания показана на рис. 2.4.

    (рис 2.4) Вероятности состояний системы М/М/4/7 с ограниченным временем ожидания

    Задание 3

  • Постройте зависимость среднего времени пребывания в очереди от интенсивности входного потока, изменяя интенсивность от $$0.Х$$ до $$2.Х$$, где $$Х$$ — номер компьютера (1, 2, 3, ...), за которым выполняется лабораторная работа.
  • Напишите программу для анализа системы с параметрами, вводимыми с клавиатуры пользователем, т. е. $$\lambda$$, $$\mu$$, $$\nu$$, $$m$$, $$K$$. Предусмотрите также ввод интервала интегрирования дифференциальных уравнений.
  • Для решения примера в системе GPSS/PC подготовим данные о функциях распределения входного потока, экспоненциальном обслуживании и экспоненциальном уходе из очереди "нетерпеливых" требований.

    Программа в MATLAB формирования данных:

    clear all,clc
    L = 3.52; %% Интенсивность входного потока
    M = 0.678; %% Интенсивность обслуживания
    v = 1.3; %% Интенсивность ухода из очереди
    %% Функции экспоненциального распределения
    x = 0 : 0.2 : 20;
    F = 1 - exp(-L*x); 
    F2 = 1 - exp(-M*x);
    F3 = 1 - exp(-v*x);
    %%% Запись в текстовый файл MC2.txt
    fid = fopen('MC2.txt', 'w');
    %% Для блока generate
        fprintf(fid, 'puas  function RN1,C50\r\n');
    for J = 1 : length(x)
        if J <= 50
            if mod(J, 5) 
    fprintf(fid, '%g,%g/', F(J), x(J));    
            
            else
    fprintf(fid, '\r\n');  
    fprintf(fid, '%g,%g/', F(J), x(J));    
             end
        end
    end
    fprintf(fid, '\r\n;-----------------------------------\r\n');
    %% Для блоков advance
    fprintf(fid, '\r\nexpM  function RN2,C50\r\n');
    for J = 1 : length(x)
        if J <= 50
            if mod(J, 5) 
    fprintf(fid, '%g,%g/', F2(J), x(J));    
            
            else
    fprintf(fid, '\r\n');  
    fprintf(fid, '%g,%g/', F2(J), x(J));  
             end
        end
    end
    
    fprintf(fid, '\r\n;-----------------------------------\r\n');
    %% Для функции распределения ухода из очереди
    fprintf(fid, '\r\nexpM  function RN3,C50\r\n');
    for J = 1 : length(x)
        if J <= 50
            if mod(J, 5) 
    fprintf(fid, '%g,%g/', F3(J), x(J));    
            
            else
    fprintf(fid, '\r\n');  
    fprintf(fid, '%g,%g/', F3(J), x(J));  
             end
        end
    end
    fclose(fid);

    GPSS-программа приводится ниже. В ней задано условие обработки 1000 требований.

    Программный код решения примера в системе GPSS/PC:

    simulate
    puas  function RN1,C50
    0,0/0.505397,0.2/0.755368,0.4/0.879004,0.6/
    0.940155,0.8/0.970401,1/0.98536,1.2/0.992759,1.4/0.996419,1.6/
    0.998229,1.8/0.999124,2/0.999567,2.2/0.999786,2.4/0.999894,2.6/
    0.999948,2.8/0.999974,3/0.999987,3.2/0.999994,3.4/0.999997,3.6/
    0.999998,3.8/0.999999,4/1,4.2/1,4.4/1,4.6/
    1,4.8/1,5/1,5.2/1,5.4/1,5.6/
    1,5.8/1,6/1,6.2/1,6.4/1,6.6/
    1,6.8/1,7/1,7.2/1,7.4/1,7.6/
    1,7.8/1,8/1,8.2/1,8.4/1,8.6/
    1,8.8/1,9/1,9.2/1,9.4/1,9.6/
    1,9.8/
    ;-----------------------------------
    expM  function RN2,C50
    0,0/0.126808,0.2/0.237536,0.4/0.334223,0.6/
    0.418649,0.8/0.492369,1/0.556741,1.2/0.612949,1.4/0.662031,1.6/
    0.704888,1.8/0.742311,2/0.774988,2.2/0.803521,2.4/0.828436,2.6/
    0.850192,2.8/0.869189,3/0.885777,3.2/0.900261,3.4/0.912909,3.6/
    0.923953,3.8/0.933596,4/0.942017,4.2/0.949369,4.4/0.95579,4.6/
    0.961396,4.8/0.966291,5/0.970566,5.2/0.974298,5.4/0.977558,5.6/
    0.980403,5.8/0.982888,6/0.985058,6.2/0.986953,6.4/0.988607,6.6/
    0.990052,6.8/0.991314,7/0.992415,7.2/0.993377,7.4/0.994217,7.6/
    0.99495,7.8/0.995591,8/0.99615,8.2/0.996638,8.4/0.997064,8.6/
    0.997437,8.8/0.997762,9/0.998045,9.2/0.998293,9.4/0.99851,9.6/
    0.998699,9.8/
    ;-----------------------------------
    expV  function RN3,C50
    0,0/0.228948,0.2/0.405479,0.4/0.541594,0.6/
    0.646545,0.8/0.727468,1/0.789864,1.2/0.837974,1.4/0.87507,1.6/
    0.903672,1.8/0.925726,2/0.942731,2.2/0.955843,2.4/0.965953,2.6/
    0.973748,2.8/0.979758,3/0.984392,3.2/0.987966,3.4/0.990721,3.6/
    0.992845,3.8/0.994483,4/0.995746,4.2/0.99672,4.4/0.997471,4.6/
    0.99805,4.8/0.998497,5/0.998841,5.2/0.999106,5.4/0.999311,5.6/
    0.999469,5.8/0.99959,6/0.999684,6.2/0.999756,6.4/0.999812,6.6/
    0.999855,6.8/0.999888,7/0.999914,7.2/0.999934,7.4/0.999949,7.6/
    0.999961,7.8/0.99997,8/0.999977,8.2/0.999982,8.4/0.999986,8.6/
    0.999989,8.8/0.999992,9/0.999994,9.2/0.999995,9.4/0.999996,9.6/
    0.999997,9.8/
    
    ***************************************
    T_s table m1,0,2,6
    **************************************
    10  GENERATE fn$puas
    20  test LE  q1,3,exit1
    30  test L  qx1,fn$expV,exit2
    40  QUEUE 1
    50  TRANSFER  ALL,CHAN1,CHAN4,3
    
    **************************************
    100 CHAN1  SEIZE zet1
    110  assign 1,zet1
    120  transfer  ,come
    **************************************
    200 chan2   SEIZE zet2
    210  assign 1,zet2
    220  transfer  ,come
    **************************************
    300 chan3   SEIZE zet3
    310  assign 1,zet3
    320  transfer  ,come
    **************************************
    400 chan4   SEIZE zet4
    410  assign 1,zet4
    
    **************************************
    510 come  DEPART 1
    520  advance fn$expM
    530  release p1
    540  tabulate T_s
    550 work   terminate 1
    560 exit1  terminate
    570 exit2  terminate
    
    **************************************
    START 1000
    ;end

    В программе использованы следующие стандартные числовые атрибуты: $$q1$$ — длина первой очереди, $$qx1$$ — среднее время пребывания в первой очереди без нулевых входов, $$m1$$ — время пребывания транзакта в модели. Освобождение того или иного устройства блоком $$release$$ выполнено на основе записи в первый параметр транзактов имени устройства, например, $$assign 1,zet4$$.

    Программа реализует следующий алгоритм. Поступающие по пуассоновскому закону требования (транзакты в системах GPSS) попадают в очередь на обслуживание. Если в очереди меньше или равно 3 транзактам, то они ожидают начала обслуживания. В противном случае они покидают систему по метке $$exit1$$. Проверка длины очереди осуществляется блоком $$TEST$$ в строке 20. Время ожидания соотносится со случайным временем ожидания "нетерпеливых" требований —— со значением функции распределения $$expV$$. Если оно меньше, то требования поступают в устройства обслуживания, в одно из свободных устройств, начиная с первого (реализуется блоком $$TRANSFER$$ в режиме $$ALL$$ ). Поступающим в устройства обслуживания транзактам в первый параметр записывается имя устройства ( $$zet1, … , zet4$$ ), после чего они отправляются по метке $$come$$ в блок $$DEPART$$, который освобождает одно место в очереди под номером 1. Далее происходит обслуживание транзактов по экспоненциальному закону блоком $$ADVANCE$$ в строке 520. Вывод транзакта из работающего устройства осуществляется блоком $$RELEASE$$ по первому параметру, в котором записано имя устройства. Обслуженные транзакты выводятся из системы блоком $$TERMINATE$$ с меткой $$work$$.

    Для получения распределения времени пребывания в модели использован оператор таблицы $$table$$ с меткой $$T_s$$. Сбор статистики для введенной таблицы осуществляется блоком $$TABULATE$$.

    Результат выполнения программы система GPSS/PC оформляет в виде файла стандартного отчета, который имеет следующий вид:

    GPSS/PC Report file REPORT.GPS.  (V 2, # 37349)  01-08-2010 11:16:14   page 1
    
         START_TIME    END_TIME  BLOCKS    FACILITIES  STORAGES   FREE_MEMORY
                0          416     23           4          0        137072
    
     LINE       LOC          BLOCK_TYPE       ENTRY_COUNT   CURRENT_COUNT   RETRY
      10         1            GENERATE            14227              0         0
      20         2            TEST                14227              0         0
      30         3            TEST                13127              0         0
      40         4            QUEUE                1002              0         0
      50         5            TRANSFER             1288              1         0
      100       CHAN1         SEIZE                 317              0         0
      110        7            ASSIGN                317              0         0
      120        8            TRANSFER              317              0         0
      200       CHAN2         SEIZE                 265              0         0
      210        10           ASSIGN                265              0         0
      220        11           TRANSFER              265              0         0
      300       CHAN3         SEIZE                 232              0         0
      310        13           ASSIGN                232              0         0
      320        14           TRANSFER              232              0         0
      400       CHAN4         SEIZE                 187              0         0
      410        16           ASSIGN                187              0         0
      510       COME          DEPART               1001              0         0
      520        18           ADVANCE              1001              1         0
      530        19           RELEASE              1000              0         0
      540        20           TABULATE             1000              0         0
      550       WORK          TERMINATE            1000              0         0
      560       EXIT1         TERMINATE            1100              0         0
      570       EXIT2         TERMINATE           12125              0         0
    
    
    FACILITY   ENTRIES  UTIL.   AVE._TIME AVAILABLE  OWNER PEND INTER RETRY DELAY
     ZET1          317  0.776        1.02      1         0   0     0     1     0
     ZET2          265  0.675        1.06      1     14161   0     0     1     0
     ZET3          232  0.557        1.00      1         0   0     0     1     0
     ZET4          187  0.507        1.13      1         0   0     0     1     0
    
    
    QUEUE        MAX  CONT. ENTRIES ENTRIES(0) AVE.CONT. AVE.TIME  AVE.(-0) RETRY
      1            4     1    1002       765     0.72        0.30      1.27    0
    
    
    TABLE            MEAN      STD.DEV. RETRY  RANGE           FREQUENCY  CUM.%
     T_S              1.34        1.58    0 
                                              -         0         374      37.40
                                            0 -         2         448      82.20
                                            2 -         4         122      94.40
                                            4 -         6          42      98.60
                                            6 -         8          14     100.00
    
    
    XACT_GROUP         GROUP_SIZE       RETRY
     POSITION                  0           0

    По файлу стандартного отчета определим следующие операционные характеристики: среднее время пребывания в очереди и среднее время пребывания в модели (в системе). Эти характеристики обозначим $$Tq$$ и $$Ts$$ соответственно.

    Значение $$Tq$$ расположено в поле $$AVE.TIME$$ статистики очереди. Оно равно $$Tq = 0.30$$.

    Значение Ts есть среднее значение таблицы $$T_S$$, т. е. в поле $$MEAN$$. Оно равно $$Ts = 1.34$$.

    Поскольку очередь имеет нулевые входы – поле $$ENTRIES(0)$$, для расчета средней длины очереди $$Nq$$ используем аналог формулы Литтла [2], т. е.

    $$Nq = \lambda \cdot Tq = 3.52\cdot 0.30 = 1.0560.$$

    Для определения среднего количества требований $$Ns$$ в системе используем формулу Литтла:

    $$Ns = \lambda \cdot Ts = 3.52\cdot 1.34 = 4.7168.$$

    Задание 4

  • По данным файла стандартного отчета определите вероятность отказа, относительную пропускную способность, абсолютную пропускную способность.
  • Если требуется, то внесите изменения в программу для определения ве-роятности наличия очереди и вероятности загрузки всех каналов обслуживания.
  • Внесите изменения в программу для определения вероятностей состояния системы.
  • Полученные результаты сведите в таблицу, в которую также внесите операционные характеристики системы $$М/М/4/7$$ с ограниченным временем ожидания, полученные с помощью анализа в MATLAB.
  • 2.3. Пример моделирования системы типа М/М/M/K/M

    Система $$M/M/m/K/M$$ — это система с пуассоновским входящим пото-ком требований, с экспоненциальным законом обслуживания в m приборах, с допустимым числом требований в системе, не превышающим K, и с ограниченным числом источников нагрузки, которые создают поток из $$M$$ требований. Общее число K требований в системе заключено в интервале $$m \le K \le M$$, где $$M$$ — число требований, формируемых конечным числом источников нагрузки.

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

    $$\lambda_{k}= \left \begin{cases} \lambda(M-k), 0\le k \le K-1, \lambda=const;\\ 0\mbox{ }\mbox{ в остальных случаях}\\ \end{cases} \right.\\ \\ \\ \mu_{k}= \left \begin{cases} k\mu, 0\le k \le m, \mu=const;\\ m\mu, m\le k \le K.\\ \end{cases} \right.$$

    Диаграмма интенсивностей переходов для системы $$M/M/m/K/M$$ будет представлять собой конечный размеченный граф состояний, который показан на рис. 2.5.

    (рис 2.5) Граф состояний системы M/M/m/K/M

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

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

    Запишем следующее выражение для определения стационарных вероятностей состояний от $$0$$ до $$m – 1$$:

    $$p_{k}=p_{0}\rho^{k}C_{M}^{k},\mbox{ }0\le k\le m-1,$$

    где:

    $$C_{M}^{k}=\frac{M!}{k!(M-k)!}$$ — биномиальные коэффициенты;

    $$\rho = \lambda / \mu$$,

    $$p_{0}=\frac{1}{\sum\limits_{k=0}^{m-1} \rho^{k}C_{M}^{k}+\sum\limits_{k=m}^{K} \rho^{k}C_{M}^{k}\frac{k!m^{m-k}}{m!}} $$ — вероятность нулевого состояния.

    Расчет вероятностей состояний от $$m$$ до $$K$$:

    $$p_{m+r}=p_{0}\rho^{m+r}C_{M}^{m+r}\frac{(m+r)!}{m!m^{r}},\mbox{ }r=0,1,2,...,K-m.$$

    Если в последнем выражении сделать замену $$m+r=k$$, $$r=k-m$$, то получим

    $$p_{k}=p_{0}\rho^{k}C_{M}^{k}\frac{k!m^{m-k}}{m!},\mbox{ }m\le k\le K.$$

    Вероятность нулевого состояния была определена из нормировочного условия:

    $$\sum\limits_{k=0}^{K} p_{k}=1.$$

    Рассмотрим случай, когда система $$M/M/m/K/M$$ работает в режиме чистых потерь, т. е. когда параметры системы удовлетворяют условию $$M\ge K=m$$. Расчет вероятностей состояний будет определяться при $$K = m$$:

    $$p_{k}=\frac{\rho^{k}C_{M}^{k}}{\sum\limits_{i=0}^{m}\rho^{i}C_{M}^{i}},\mbox{ }k=0,1,2,...,m.$$

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

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

    $$\frac{dP_{0}}{dt}=-M\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=M\lambda P_{0}-[(M-1)\lambda +\mu] P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=(M-1)\lambda P_{1}-[(M-2)\lambda +2\mu] P_{2}+3\mu P_{3};\\ .............................\\ \frac{dP_{m}}{dt}=[M-(m-1)]\lambda P_{m-1}-[(M-m)\lambda +m\mu] P_{m}+m\mu P_{m+1};\\ .............................\\ \frac{dP_{K-1}}{dt}=[M-(K-2)]\lambda P_{K-2}-\{[M-(KI-1)]\lambda +m\mu\} P_{K-1}+m\mu P_{K};\\ \frac{dP_{K}}{dt}=[M-(K-1)]\lambda P_{K-1}-m\mu P_{K}.$$

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

    $$P_{0}(0)=1,\qquad P_{j}(0)=0,\qquad j=1,2,...,K.$$

    Для стационарного режима рассмотрим ряд операционных характеристик в достаточно общем виде.

    Среднее число требований в системе:

    $$N_{ср}=\sum\limits_{k=0}^{K}kp_{k}.$$

    Среднее время пребывания одного требования в системе определим по формуле Литтла:

    $$T=\frac{N_{ср}}{\lambda}.$$

    Средняя длина очереди:

    $$N_{q}=\sum\limits_{k=m}^{K}kp_{k+m}.$$

    Среднее время пребывания требования в очереди определим по формуле Литтла:

    $$T_{q}=\frac{N_{q}}{\lambda}.$$

    Вероятность отказа в обслуживании соответствует вероятности того, что в системе находится $$K$$ требований (максимально допустимое число):

    $$p_{отк}=p_{K}=p_{0}\rho^{K}C_{M}^{K}\frac{K!m^{m-K}}{m!}.$$

    Относительная пропускная способность:

    $$Q=1-p_{отк}.$$

    Абсолютная пропускная способность:

    $$A=\lambda Q.$$

    Таким образом, операционные характеристики рассчитываются по известным стационарным вероятностям и параметрам системы.

    Практическая часть

    Пример 3. Для системы массового обслуживания $$М/М/4/7/9$$ постройте переходные вероятности состояний и определите для нее операционные характеристики при следующих параметрах: $$\lambda = 3.52 с^{–1}$$, $$\mu = 0.678 с^{–1}$$. Среда программирования — MATLAB.

    Запишем с буквенными коэффициентами дифференциальные уравнения для системы $$М/М/4/7/9$$:

    $$\frac{dP_{0}}{dt}=-M\lambda P_{0}+\mu P_{1},\\ \frac{dP_{1}}{dt}=M\lambda P_{0}-[(M-1)\lambda +\mu] P_{1}+2\mu P_{2},\\ \frac{dP_{2}}{dt}=(M-1)\lambda P_{1}-[(M-2)\lambda +2\mu] P_{2}+3\mu P_{3},\\ \frac{dP_{3}}{dt}=(M-2)\lambda P_{2}-[(M-3)\lambda +3\mu] P_{3}+4\mu P_{4},\\ \frac{dP_{4}}{dt}=(M-3)\lambda P_{3}-[(M-4)\lambda +4\mu] P_{4}+4\mu P_{5},\\ \frac{dP_{5}}{dt}=(M-4)\lambda P_{4}-[(M-5)\lambda +4\mu] P_{5}+4\mu P_{6},\\ \frac{dP_{6}}{dt}=(M-5)\lambda P_{5}-[(M-6)\lambda +4\mu] P_{6}+4\mu P_{7},\\ \frac{dP_{7}}{dt}=[M-(K-1)]\lambda P_{6}-4\mu P_{7},\\$$

    где:

    $$М = 9; K = 7; \lambda = 3.52 с^{–1}; \mu = 0.678 с^{–1}.$$

    Полученную систему представим в матричном виде:

    $$\frac{dP}{dt}=AP,$$

    где матрица $$А$$ размера $$8\times 8$$ составляется из коэффициентов при вероятностях в правой части уравнений.

    Для интегрирования дифференциальных уравнений примем естественные граничные условия.

    Программный код решения примера в MATLAB:

    function MMmKM;
    clc,close all
    % Параметры системы M/M/4/7/9
    L = 2.2;
    M = 3.678;
    m = 4;
    K = 7;
    N = 9; 
    % Матрица коэффициентов А
    global A
    A = [-N*L,M,zeros(1,6);
    N*L,-[(N-1)*L+M],2*M,zeros(1,5);
    zeros(1,1),(N-1)*L,-[(N-2)*L+2*M],3*M,zeros(1,4);
    zeros(1,2),(N-2)*L,-[(N-3)*L+3*M],4*M,zeros(1,3);
    zeros(1,3),(N-3)*L,-[(N-4)*L+4*M],4*M,zeros(1,2);
    zeros(1,4),(N-4)*L,-[(N-5)*L+4*M],4*M,zeros(1,1);
    zeros(1,5),(N-5)*L,-[(N-6)*L+4*M],4*M;
    zeros(1,6),(N-6)*L,-4*M];
     
    T = [0,2]; %% интервал интегрирования
    P0 = [1,zeros(1,length(A)-1)]; %% начальные условия
    [t,P] = ode23(@cmo,T,P0);
     
    %% Построение диаграмм вероятностей состояний
    % line(t,P,'linew',2) %% сплошные линии с различными цветами
    line(t,P(:,1),'linew',2, 'color','r') %% Po
    line(t,P(:,2), 'linew',2,'lines','--') %% P1
    line(t,P(:,3), 'linew',2,'lines','-.') %% P2
    line(t,P(:,4), 'linew',2,'lines',':') %% P3
    line(t,P(:,5), 'marker','o', 'color', 'm') %% P4
    line(t,P(:,6), 'marker','h', 'color','k') %% P5
    line(t,P(:,7), 'marker','p','color','r') %% P6
    line(t,P(:,8), 'marker','>') %% P7
    grid on
    Na = length(A) - 1;
    arr = [0:Na]';
    str = num2str(arr);
    legend(strcat('\bf\itP\rm\bf_', str, '(\itt\rm\bf)'));
    
    title(sprintf('%sСистема M/M/%d/%d/%d; %s%g; %s%g',...
    '\bf\fontsize{11}',m,K,N,'\lambda = ',L,'\mu = ',M));
    xlabel('\bf\it\fontsize{12} -  -  -  -  -  -  -  -   t   -  -  -  -  -  -  -  -');
    ylabel('\bf\fontsize{12}\itP\rm\bf(\itt\rm\bf)');
    set(gca, 'fontweight','bold', 'fontsize',10)
     
    Pcm = P(end,:); % Стационарные вероятности
    fprintf('\n Стационарные вероятности системы M/M/%d/%d/%d:\n', m,K,N);
    for J = 1 : length(A)
        fprintf('\tP%d = %f\n', J-1, Pcm(J));
    end
    fprintf('\n ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ СИСТЕМЫ M/M/%d/%d/%d:\n', m,K,N);
    Pnot = P(end,end);
    fprintf(' Вероятность отказа Pnot = %f\n', P(end,end));
    Q = 1 - Pnot;
    fprintf(' Относительная пропускная способность Q = %f\n', Q);
    Ab = L*Q;
    fprintf(' Абсолютная пропускная способность A = %f\n', Ab);
    Pq = sum(P(end, m+1:end));
    fprintf(' Вероятность наличия очереди Pq = %f\n', Pq);
    Ps = sum(P(end, m:end));
    fprintf(' Вероятность загрузки всех каналов обслуживания Ps = %f\n', Ps);
    Ns = [0:length(A)-1]*P(end,:)';
    fprintf(' Среднее количество требований в системе Ns = %f\n', Ns);
    fprintf(' Среднее время пребывания требования в системе Ts = %f\n', Ns/L);
    Nq = [0:(K-m)]*P(end,m:K)';
    fprintf(' Средняя длина очереди Nq = %f\n', Nq);
    fprintf(' Среднее время пребывания требования в очереди Tq = %f\n', Nq/L);
     
    function f = cmo(t,P)
    % М-функция описания правых частей дифференциальных уравнений:
    global A
    f = A*P;

    Результат выполнения программы

    Стационарные вероятности системы M/M/4/7/9:
    	P0 = 0.013187
    	P1 = 0.070956
    	P2 = 0.169841
    	P3 = 0.236851
    	P4 = 0.212652
    	P5 = 0.158847
    	P6 = 0.095048
    	P7 = 0.042619
    
     ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ СИСТЕМЫ M/M/4/7/9
     Вероятность отказа: Pnot = 0.042619
     Относительная пропускная способность: Q = 0.957381
     Абсолютная пропускная способность: A = 2.106238
     Вероятность наличия очереди: Pq = 0.509166
     Вероятность загрузки всех каналов обслуживания: Ps = 0.746017
     Среднее количество требований в системе: Ns = 3.634652
     Среднее время пребывания требования в системе: Ts = 1.652115
     Средняя длина очереди: Nq = 0.815489
     Среднее время пребывания требования в очереди: Tq = 0.370677

    Зависимости вероятностей состояний от времени показаны на рис. 2.6.

    (рис 2.6) Изменение вероятностей состояний во времени системы М/М/4/7/9

    Задание 5

  • Рассчитайте стационарные вероятности по аналитическим формулам (см. теоретическую часть). Сравните результаты.
  • Постройте зависимость среднего времени пребывания в очереди от интенсивности входного потока, изменяя интенсивность от $$0.Х$$ до $$5.Х$$, где $$Х$$ — номер компьютера (1, 2, 3, ...), за которым выполняется лабораторная работа.
  • Напишите программу для анализа системы с параметрами, вводимыми с клавиатуры пользователем, т. е. $$\lambda, \mu, m, K, М$$. Предусмотрите также ввод интервала интегрирования дифференциальных уравнений.
  • Контрольные вопросы

  • Что называется простейшим пуассоновским потоком?
  • По какому закону распределены интервалы времени между требованиями в простейшем пуассоновском потоке?
  • Что такое символика Кендалла, применяемая для систем массового обслуживания?
  • Что определяет собой точка на графике функции распределения случайной величины?
  • Какие случайные величины в системе массового обслуживания являются дискретными, а какими непрерывными?
  • Что означает несовместность событий?
  • Каким условиям должны отвечать начальные условия при решении дифференциальных уравнений относительно вероятностей состояния системы массового обслуживания?
  • Какова длина очереди в системе $$M/M/m/K/M$$?
  • Через какой тип данных могут выражаться интенсивность входного потока требований и интенсивность обслуживания требований в системах массового обслуживания?
  • Страницы:

    Теоретическая часть

    Многоканальные системы массового обслуживания — это системы с параллельно включенными приборами обслуживания. Для них принято использовать символику Кендалла, которая состоит из основных четырех позиций вида $$A/B/m/K/М$$, где $$А$$ — закон поступления требований в систему, $$В$$ — закон обслуживания требований, m — число параллельно функционирующих приборов (каналов, узлов) обслуживания, $$K$$ — допустимое число требований в системе, т. е. число требований в очереди плюс число требований, принятых на обслуживание. В системе $$A/B/m/K/М$$ последний символ $$М$$ — число источников нагрузки [2]. При этом параметры системы постоянны, а процессы в системах рассматриваются с позиций теории размножения и гибели.

    Предметом данной лабораторной работы будут системы с простейшим пуассоновским потоком, экспоненциальным обслуживанием и, возможно, экспоненциальном уходом из очереди "нетерпеливых" требований. При сделанных условиях исследуемые системы будут иметь следующие обозначения: $$M/M/m/K$$, $$M/M/m/K$$ с нетерпеливыми требованиями, $$M/M/m/K/M$$, где первые буквы $$М$$ означают пуассоновский поток требований и экспоненциальное обслуживание, а в последнем примере крайняя правая буква $$М$$ — это число требований, формируемых конечным числом источников нагрузки. В частности, может быть система с отказами, т. е. $$M/M/m$$ с отказами.

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

    Для составления уравнений Колмогорова можно использовать мнемоническое правило [4]. Сначала определим понятие потока вероятности: это такой поток, который переводит систему из одного состояния в другое соседнее и определяется как произведение вероятности $$P_{i}(t)$$ $$i$$ -го состояния, из которого происходит переход, на интенсивность потока событий (интенсивность поступления требований или интенсивность обслуживания). Теперь приведем мнемоническое правило составления уравнений Колмогорова: производная вероятности любого состояния равна сумме потоков вероятности, переводящих систему в это состояние, минус сумма всех потоков вероятности, выводящих систему из этого состояния [4].

    Прежде чем применять правило Колмогорова, целесообразно изобразить размеченный граф состояний заданной системы. На рис. 2.1 показан пример размеченного графа состояний для системы M/M/3/5 ( $$m = 3$$, $$K = 5$$, $$\lambda$$ — интенсивность входного потока требований, $$\mu$$ — интенсивность обслуживания одним прибором).

    (рис 2.1) Размеченный граф состояний системы М/М/3/5

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

    $$\frac{dP_{0}}{dt}=-\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=\lambda P_{0}-(\lambda +\mu) P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=\lambda P_{1}-(\lambda +2\mu) P_{2}+3\mu P_{3};\\ \frac{dP_{3}}{dt}=\lambda P_{2}-(\lambda +3\mu) P_{3}+3\mu P_{4};\\ \frac{dP_{4}}{dt}=\lambda P_{3}-(\lambda +3\mu) P_{4}+3\mu P_{5};\\ \frac{dP_{5}}{dt}=\lambda P_{4}-3\mu P_{5};\\$$

    Для решения систем дифференциальных уравнений типа (2.1) обычно задают естественные начальные условия:

    $$P_{0}(0)=1,\qquad P_{i}(0)=0,\qquad i=\overline{1,5}.$$

    Если в системе (2.1) производные приравнять нулю, то можно будет получить соотношения для расчета стационарных вероятностей состояний системы. При этом следует использовать нормировочное условие для $$n = 5$$:

    $$\sum\limits_{k=0}^{5} P_{k}(t)=1.$$

    Практическая часть

    2.1. Пример моделирования системы типа М/М/M/K

    Система $$M/M/m/K$$ — это система с пуассоновским входящим потоком требований, с экспоненциальным законом обслуживания в m приборах, с допустимым числом требований в системе, не превышающим $$K$$, которое не менее, чем заданное количество приборов обслуживания. Параметры системы постоянны, т. е. $$\lambda = const$$, $$\mu = const$$. Систему $$M/M/m/K$$ называют многоканальной системой массового обслуживания с ограниченной длиной очереди.

    Пример 1. Проинтегрируйте систему массового обслуживания $$М/М/4/7$$ при естественных граничных условиях и параметрах $$\lambda = 1.23$$, $$\mu = 0.678$$. Рассчитайте операционные характеристики системы [2].

    Для решения примера сначала составим дифференциальные уравнения в соответствии с мнемоническим правилом Колмогорова:

    $$\frac{dP_{0}}{dt}=-\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=\lambda P_{0}-(\lambda +\mu) P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=\lambda P_{1}-(\lambda +2\mu) P_{2}+3\mu P_{3};\\ \frac{dP_{3}}{dt}=\lambda P_{2}-(\lambda +3\mu) P_{3}+4\mu P_{4};\\ \frac{dP_{4}}{dt}=\lambda P_{3}-(\lambda +4\mu) P_{4}+4\mu P_{5};\\ \frac{dP_{5}}{dt}=\lambda P_{4}-(\lambda +4\mu) P_{5}+4\mu P_{6};\\ \frac{dP_{6}}{dt}=\lambda P_{5}-(\lambda +4\mu) P_{6}+4\mu P_{7};\\ \frac{dP_{7}}{dt}=\lambda P_{6}-4\mu P_{7};\\$$

    Естественные начальные условия:

    $$P_{0}(0)=1,\qquad P_{i}(0)=0,\qquad i=\overline{1,7}.$$

    Программный код решения примера в MATLAB:

    function MMmK;
    clc,close 
    % Параметры системы
    L = 3.52;
    M = 0.678;
    m = 4;
    K = 7; 
    %%%%%%%%%%%%%%%%%%%%%%%%
    global A
    A = [-L,M,0,0,0,0,0,0;
    L,-(L+M),2*M,0,0,0,0,0;
    0,L,-(L+2*M),3*M,0,0,0,0;
    0,0,L,-(L+3*M),4*M,0,0,0;
    0,0,0,L,-(L+4*M),4*M,0,0;
    0,0,0,0,L,-(L+4*M),4*M,0;
    0,0,0,0,0,L,-(L+4*M),4*M;
    0,0,0,0,0,0,L,-4*M];
    %% Численное интегрирование дифф. уравнений
    P0 = [1;zeros(length(A)-1,1)];
    T = [0,20];
    [t,P] = ode23(@cmo, T, P0);
     
    %% Построение диаграммы вероятностей состояний
    %% line(t,P,'linew',2) %% с различными цветами
    line(t,P(:,1),'linew',2, 'color','r') %% Po
    line(t,P(:,2), 'linew',2,'lines','--') %% P1
    line(t,P(:,3), 'linew',2,'lines','-.') %% P2
    line(t,P(:,4), 'linew',2,'lines',':') %% P3
    line(t,P(:,5), 'marker','o', 'color', 'm') %% P4
    line(t,P(:,6), 'marker','h', 'color','k') %% P5
    line(t,P(:,7), 'marker','p','color','r') %% P6
    line(t,P(:,8), 'marker','>') %% P7
    grid on
    N = length(A)-1;
    arr = [0:N]';
    str = num2str(arr);
    legend(strcat('\bf\itP\rm\bf_', str, '(\itt\rm\bf)'));
    title(sprintf('%s Вероятности состояний системы M/M/%d/%d', '\bf\fontsize{12}',m, K));
    xlabel('\bf\it\fontsize{12} -  -  -  -  -  -  -  -   t   -  -  -  -  -  -  -  -')
    ylabel('\bf\fontsize{12}\itP\rm\bf(\itt\rm\bf)');
    set(gca, 'fontweight','bold', 'fontsize',10)
    
    fprintf('\n Стационарные вероятности:\n');
    for J = 1 : length(A)
        fprintf('\tP%d = %f\n', J-1, P(end,J));
    end
     fprintf('\n\t ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ:\n');
    Pnot = P(end,end);
    
    fprintf(' Вероятность отказа Pnot = %f\n', P(end,end));
    Q = 1 - Pnot;
    fprintf(' Относительная пропускная способность Q = %f\n', Q);
    Ab = L*Q;
    fprintf(' Абсолютная пропускная способность A = %f\n', Ab);
    Pq = sum(P(end, m+1:end));
    fprintf(' Вероятность наличия очереди Pq = %f\n', Pq);
    Ps = sum(P(end, m:end));
    fprintf(' Вероятность загрузки всех каналов обслуживания Ps = %f\n', Ps);
     
    Ns = [0:length(A)-1]*P(end,:)';
    fprintf(' Среднее количество требований в системе Ns = %f\n', Ns);
    fprintf(' Среднее время пребывания требования в системе Ts = %f\n', Ns/L);
    Nq = [0:(K-m)]*P(end,m:K)';
    fprintf(' Средняя длина очереди Nq = %f\n', Nq);
    fprintf(' Среднее время пребывания требования в очереди Tq = %f\n', Nq/L);
     
     
    function f = cmo(t,P);
    %% Функция описания правых частей 
    %% дифференциальных уравнений
    global A
    f = A*P;

    Результат выполнения программы

    Стационарные вероятности:
    	P0 = 0.004352
    	P1 = 0.022604
    	P2 = 0.058629
    	P3 = 0.101614
    	P4 = 0.131597
    	P5 = 0.171220
    	P6 = 0.221807
    	P7 = 0.288176
    
    	 ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ
     Вероятность отказа: Pnot = 0.288176
     Относительная пропускная способность: Q = 0.711824
     Абсолютная пропускная способность: A = 2.505620
     Вероятность наличия очереди: Pq = 0.812800
     Вероятность загрузки всех каналов обслуживания: Ps = 0.914415
     Среднее количество требований в системе: Ns = 5.175269
     Среднее время пребывания требования в системе: Ts = 1.470247
     Средняя длина очереди: Nq = 1.139458
     Среднее время пребывания требования в очереди: Tq = 0.323710

    Среднее время пребывания требования в системе было рассчитано по формуле Литтла. Для расчета среднего времени пребывания требования в очереди использовался аналог формулы Литтла [2].

    Диаграмма с вероятностями состояний системы $$М/М/4/7$$ показана на рис. 2.2.

    (рис 2.2) Вероятности состояний системы М/М/4/7

    Задание 1

  • Рассчитайте нормировочное условие для трех значений времени.
  • Определите вероятность того, что в системе будет не более $$k$$ требований. Значение $$k$$ выбирайте случайно по равномерному закону из интервала целых чисел $$[1; 6]$$ (см. $$help randperm$$ ).
  • Примените диалоговое окно $$inputdlg$$ (см. $$help inputdlg$$ ) для ввода параметров системы.
  • Напишите программу для произвольно задаваемых параметров системы $$\lambda$$, $$\mu$$ и значений $$m$$, $$K$$. Предусмотрите также ввод длительности интервала интегрирования.
  • 2.2. Пример моделирования системы типа М/М/M/K с ограниченным временем ожидания

    Пример 2. Получите операционные характеристики системы $$M/M/4/7$$ с ограниченным временем ожидания при следующих параметрах системы: $$\lambda = 3.52 \mbox{ }с^{–1}$$ — интенсивность входного потока требований, $$\mu = 0.678 \mbox{ }с^{–1}$$ — интенсивность обслуживания каждым прибором, $$\nu = 1.3\mbox{ }с^{–1}$$ — интенсивность ухода из очереди "нетерпеливых" требований. Моделирование проведите в MATLAB и GPSS/PC.

    Размеченный граф состояний системы показан на рис. 2.3.

    (рис 2.3) Граф состояний системы с ограниченным временем ожидания

    На рис. 2.3 вокруг состояния 4 проведена воображаемая окружность для пояснения применения мнемонического правила составления дифференциальных уравнений Колмогорова. В соответствии с мнемоническим правилом получаем следующую систему дифференциальных уравнений:

    $$\frac{dP_{0}}{dt}=-\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=\lambda P_{0}-(\lambda +\mu) P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=\lambda P_{1}-(\lambda +2\mu) P_{2}+3\mu P_{3};\\ \frac{dP_{3}}{dt}=\lambda P_{2}-(\lambda +3\mu) P_{3}+4\mu P_{4};\\ \frac{dP_{4}}{dt}=\lambda P_{3}-(\lambda +4\mu) P_{4}+(4\mu +\nu) P_{5};\\ \frac{dP_{5}}{dt}=\lambda P_{4}-(\lambda +4\mu+\nu) P_{5}+(4\mu +2\nu)P_{6};\\ \frac{dP_{6}}{dt}=\lambda P_{5}-(\lambda +4\mu+2\nu) P_{6}+(4\mu+3\nu) P_{7};\\ \frac{dP_{7}}{dt}=\lambda P_{6}-(4\mu+3\nu) P_{7};\\$$

    где $$P_{i},\mbox{ } i=\overline{0,7}$$ — вероятности состояний системы.

    Для интегрирования системы дифференциальных уравнений зададим естественные граничные условия:

    $$P_{0}(0)=1,\qquad P_{i}(0)=0,\qquad i=\overline{1,7}.$$

    Программный код решения примера в системе MATLAB:

    function MMmKv
    %% Система M/M/4/7 с ограниченным временем ожидания в очереди
    % Параметры системы
    L = 3.52;
    M = 0.678;
    v = 1.3;
    m = 4; 
    K = 7;
    
    % Матрица коэффициентов системы дифференциальных уравнений
    global A
    A = [-L, M, zeros(1,6);
        L, -(L+M), 2*M, zeros(1,5);
        0, L, -(L+2*M), 3*M, zeros(1,4);
        0, 0, L, -(L+3*M), 4*M, zeros(1,3);
        zeros(1,3), L, -(L+4*M), (4*M+v), 0, 0;
        zeros(1,4), L, -(L+4*M+v), (4*M+2*v), 0;
        zeros(1,5), L, -(L+4*M+2*v), (4*M+3*v);
        zeros(1,6), L ,-(4*M+3*v)];
    % Решение дифференциальных уравнений
    T = [0,20]; % Время интегрирования
    P0 = [1;zeros(length(A)-1,1)]; % Начальные условия
    [t, P] = ode23(@cmo, T, P0);
    
    %% Построение диаграмм вероятностей состояний
    % line(t,P,'linew',2) %% сполошные линии с различными цветами
    line(t,P(:,1),'linew',2, 'color','r') %% Po
    line(t,P(:,2), 'linew',2,'lines','--') %% P1
    line(t,P(:,3), 'linew',2,'lines','-.') %% P2
    line(t,P(:,4), 'linew',2,'lines',':') %% P3
    line(t,P(:,5), 'marker','o', 'color', 'm') %% P4
    line(t,P(:,6), 'marker','h', 'color','k') %% P5
    line(t,P(:,7), 'marker','p','color','r') %% P6
    line(t,P(:,8), 'marker','>') %% P7
    grid on
    N = length(A)-1;
    arr = [0:N]';
    str = num2str(arr);
    legend(strcat('\bf\itP\rm\bf_', str, '(\itt\rm\bf)'));
    title(sprintf('%s Вероятности состояний системы M/M/%d/%d', '\bf\fontsize{12}',m, K));
    xlabel('\bf\it\fontsize{12} -  -  -  -  -  -  -  -   t   -  -  -  -  -  -  -  -')
    ylabel('\bf\fontsize{12}\itP\rm\bf(\itt\rm\bf)');
    
    set(gca, 'fontweight','bold', 'fontsize',10)
    
    Pcm = P(end,:); % Стационарные вероятности
    fprintf('\n Стационарные вероятности:\n');
    for J = 1 : length(A)
        fprintf('\tP%d = %f\n', J-1, Pcm(J));
    end
    fprintf('\n\t ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ:\n');
    Pnot = P(end,end);
    fprintf(' Вероятность отказа Pnot = %f\n', P(end,end));
    Q = 1 - Pnot;
    fprintf(' Относительная пропускная способность Q = %f\n', Q);
    Ab = L*Q;
    fprintf(' Абсолютная пропускная способность A = %f\n', Ab);
    Pq = sum(P(end, m+1:end));
    fprintf(' Вероятность наличия очереди Pq = %f\n', Pq);
    Ps = sum(P(end, m:end));
    fprintf(' Вероятность загрузки всех каналов обслуживания Ps = %f\n', Ps);
    Ns = [0:length(A)-1]*P(end,:)';
    fprintf(' Среднее количество требований в системе Ns = %f\n', Ns);
    fprintf(' Среднее время пребывания требования в системе Ts = %f\n', Ns/L);
    Nq = [0:(K-m)]*P(end,m:K)';
    fprintf(' Средняя длина очереди Nq = %f\n', Nq);
    fprintf(' Среднее время пребывания требования в очереди Tq = %f\n', Nq/L);
    
    function f = cmo(t,P);
    %% Функция описания правых частей 
    %% дифференциальных уравнений
    global A
    f = A*P;

    Результат выполнения программы

    Стационарные вероятности:
    	P0 = 0.007887
    	P1 = 0.040947
    	P2 = 0.106289
    	P3 = 0.183963
    	P4 = 0.238701
    	P5 = 0.209568
    	P6 = 0.138707
    	P7 = 0.073939
    
    	 ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ
     Вероятность отказа: Pnot = 0.073939
     Относительная пропускная способность: Q = 0.926061
     Абсолютная пропускная способность: A = 3.259736
     Вероятность наличия очереди: Pq = 0.660915
     Вероятность загрузки всех каналов обслуживания: Ps = 0.844877
     Среднее количество требований в системе: Ns = 4.157868
     Среднее время пребывания требования в системе: Ts = 1.181213
     Средняя длина очереди: Nq = 1.073958
     Среднее время пребывания требования в очереди: Tq = 0.305102

    Диаграмма вероятностей состояний системы $$М/М/4/7$$ с ограниченным временем ожидания показана на рис. 2.4.

    (рис 2.4) Вероятности состояний системы М/М/4/7 с ограниченным временем ожидания

    Задание 3

  • Постройте зависимость среднего времени пребывания в очереди от интенсивности входного потока, изменяя интенсивность от $$0.Х$$ до $$2.Х$$, где $$Х$$ — номер компьютера (1, 2, 3, ...), за которым выполняется лабораторная работа.
  • Напишите программу для анализа системы с параметрами, вводимыми с клавиатуры пользователем, т. е. $$\lambda$$, $$\mu$$, $$\nu$$, $$m$$, $$K$$. Предусмотрите также ввод интервала интегрирования дифференциальных уравнений.
  • Для решения примера в системе GPSS/PC подготовим данные о функциях распределения входного потока, экспоненциальном обслуживании и экспоненциальном уходе из очереди "нетерпеливых" требований.

    Программа в MATLAB формирования данных:

    clear all,clc
    L = 3.52; %% Интенсивность входного потока
    M = 0.678; %% Интенсивность обслуживания
    v = 1.3; %% Интенсивность ухода из очереди
    %% Функции экспоненциального распределения
    x = 0 : 0.2 : 20;
    F = 1 - exp(-L*x); 
    F2 = 1 - exp(-M*x);
    F3 = 1 - exp(-v*x);
    %%% Запись в текстовый файл MC2.txt
    fid = fopen('MC2.txt', 'w');
    %% Для блока generate
        fprintf(fid, 'puas  function RN1,C50\r\n');
    for J = 1 : length(x)
        if J <= 50
            if mod(J, 5) 
    fprintf(fid, '%g,%g/', F(J), x(J));    
            
            else
    fprintf(fid, '\r\n');  
    fprintf(fid, '%g,%g/', F(J), x(J));    
             end
        end
    end
    fprintf(fid, '\r\n;-----------------------------------\r\n');
    %% Для блоков advance
    fprintf(fid, '\r\nexpM  function RN2,C50\r\n');
    for J = 1 : length(x)
        if J <= 50
            if mod(J, 5) 
    fprintf(fid, '%g,%g/', F2(J), x(J));    
            
            else
    fprintf(fid, '\r\n');  
    fprintf(fid, '%g,%g/', F2(J), x(J));  
             end
        end
    end
    
    fprintf(fid, '\r\n;-----------------------------------\r\n');
    %% Для функции распределения ухода из очереди
    fprintf(fid, '\r\nexpM  function RN3,C50\r\n');
    for J = 1 : length(x)
        if J <= 50
            if mod(J, 5) 
    fprintf(fid, '%g,%g/', F3(J), x(J));    
            
            else
    fprintf(fid, '\r\n');  
    fprintf(fid, '%g,%g/', F3(J), x(J));  
             end
        end
    end
    fclose(fid);

    GPSS-программа приводится ниже. В ней задано условие обработки 1000 требований.

    Программный код решения примера в системе GPSS/PC:

    simulate
    puas  function RN1,C50
    0,0/0.505397,0.2/0.755368,0.4/0.879004,0.6/
    0.940155,0.8/0.970401,1/0.98536,1.2/0.992759,1.4/0.996419,1.6/
    0.998229,1.8/0.999124,2/0.999567,2.2/0.999786,2.4/0.999894,2.6/
    0.999948,2.8/0.999974,3/0.999987,3.2/0.999994,3.4/0.999997,3.6/
    0.999998,3.8/0.999999,4/1,4.2/1,4.4/1,4.6/
    1,4.8/1,5/1,5.2/1,5.4/1,5.6/
    1,5.8/1,6/1,6.2/1,6.4/1,6.6/
    1,6.8/1,7/1,7.2/1,7.4/1,7.6/
    1,7.8/1,8/1,8.2/1,8.4/1,8.6/
    1,8.8/1,9/1,9.2/1,9.4/1,9.6/
    1,9.8/
    ;-----------------------------------
    expM  function RN2,C50
    0,0/0.126808,0.2/0.237536,0.4/0.334223,0.6/
    0.418649,0.8/0.492369,1/0.556741,1.2/0.612949,1.4/0.662031,1.6/
    0.704888,1.8/0.742311,2/0.774988,2.2/0.803521,2.4/0.828436,2.6/
    0.850192,2.8/0.869189,3/0.885777,3.2/0.900261,3.4/0.912909,3.6/
    0.923953,3.8/0.933596,4/0.942017,4.2/0.949369,4.4/0.95579,4.6/
    0.961396,4.8/0.966291,5/0.970566,5.2/0.974298,5.4/0.977558,5.6/
    0.980403,5.8/0.982888,6/0.985058,6.2/0.986953,6.4/0.988607,6.6/
    0.990052,6.8/0.991314,7/0.992415,7.2/0.993377,7.4/0.994217,7.6/
    0.99495,7.8/0.995591,8/0.99615,8.2/0.996638,8.4/0.997064,8.6/
    0.997437,8.8/0.997762,9/0.998045,9.2/0.998293,9.4/0.99851,9.6/
    0.998699,9.8/
    ;-----------------------------------
    expV  function RN3,C50
    0,0/0.228948,0.2/0.405479,0.4/0.541594,0.6/
    0.646545,0.8/0.727468,1/0.789864,1.2/0.837974,1.4/0.87507,1.6/
    0.903672,1.8/0.925726,2/0.942731,2.2/0.955843,2.4/0.965953,2.6/
    0.973748,2.8/0.979758,3/0.984392,3.2/0.987966,3.4/0.990721,3.6/
    0.992845,3.8/0.994483,4/0.995746,4.2/0.99672,4.4/0.997471,4.6/
    0.99805,4.8/0.998497,5/0.998841,5.2/0.999106,5.4/0.999311,5.6/
    0.999469,5.8/0.99959,6/0.999684,6.2/0.999756,6.4/0.999812,6.6/
    0.999855,6.8/0.999888,7/0.999914,7.2/0.999934,7.4/0.999949,7.6/
    0.999961,7.8/0.99997,8/0.999977,8.2/0.999982,8.4/0.999986,8.6/
    0.999989,8.8/0.999992,9/0.999994,9.2/0.999995,9.4/0.999996,9.6/
    0.999997,9.8/
    
    ***************************************
    T_s table m1,0,2,6
    **************************************
    10  GENERATE fn$puas
    20  test LE  q1,3,exit1
    30  test L  qx1,fn$expV,exit2
    40  QUEUE 1
    50  TRANSFER  ALL,CHAN1,CHAN4,3
    
    **************************************
    100 CHAN1  SEIZE zet1
    110  assign 1,zet1
    120  transfer  ,come
    **************************************
    200 chan2   SEIZE zet2
    210  assign 1,zet2
    220  transfer  ,come
    **************************************
    300 chan3   SEIZE zet3
    310  assign 1,zet3
    320  transfer  ,come
    **************************************
    400 chan4   SEIZE zet4
    410  assign 1,zet4
    
    **************************************
    510 come  DEPART 1
    520  advance fn$expM
    530  release p1
    540  tabulate T_s
    550 work   terminate 1
    560 exit1  terminate
    570 exit2  terminate
    
    **************************************
    START 1000
    ;end

    В программе использованы следующие стандартные числовые атрибуты: $$q1$$ — длина первой очереди, $$qx1$$ — среднее время пребывания в первой очереди без нулевых входов, $$m1$$ — время пребывания транзакта в модели. Освобождение того или иного устройства блоком $$release$$ выполнено на основе записи в первый параметр транзактов имени устройства, например, $$assign 1,zet4$$.

    Программа реализует следующий алгоритм. Поступающие по пуассоновскому закону требования (транзакты в системах GPSS) попадают в очередь на обслуживание. Если в очереди меньше или равно 3 транзактам, то они ожидают начала обслуживания. В противном случае они покидают систему по метке $$exit1$$. Проверка длины очереди осуществляется блоком $$TEST$$ в строке 20. Время ожидания соотносится со случайным временем ожидания "нетерпеливых" требований —— со значением функции распределения $$expV$$. Если оно меньше, то требования поступают в устройства обслуживания, в одно из свободных устройств, начиная с первого (реализуется блоком $$TRANSFER$$ в режиме $$ALL$$ ). Поступающим в устройства обслуживания транзактам в первый параметр записывается имя устройства ( $$zet1, … , zet4$$ ), после чего они отправляются по метке $$come$$ в блок $$DEPART$$, который освобождает одно место в очереди под номером 1. Далее происходит обслуживание транзактов по экспоненциальному закону блоком $$ADVANCE$$ в строке 520. Вывод транзакта из работающего устройства осуществляется блоком $$RELEASE$$ по первому параметру, в котором записано имя устройства. Обслуженные транзакты выводятся из системы блоком $$TERMINATE$$ с меткой $$work$$.

    Для получения распределения времени пребывания в модели использован оператор таблицы $$table$$ с меткой $$T_s$$. Сбор статистики для введенной таблицы осуществляется блоком $$TABULATE$$.

    Результат выполнения программы система GPSS/PC оформляет в виде файла стандартного отчета, который имеет следующий вид:

    GPSS/PC Report file REPORT.GPS.  (V 2, # 37349)  01-08-2010 11:16:14   page 1
    
         START_TIME    END_TIME  BLOCKS    FACILITIES  STORAGES   FREE_MEMORY
                0          416     23           4          0        137072
    
     LINE       LOC          BLOCK_TYPE       ENTRY_COUNT   CURRENT_COUNT   RETRY
      10         1            GENERATE            14227              0         0
      20         2            TEST                14227              0         0
      30         3            TEST                13127              0         0
      40         4            QUEUE                1002              0         0
      50         5            TRANSFER             1288              1         0
      100       CHAN1         SEIZE                 317              0         0
      110        7            ASSIGN                317              0         0
      120        8            TRANSFER              317              0         0
      200       CHAN2         SEIZE                 265              0         0
      210        10           ASSIGN                265              0         0
      220        11           TRANSFER              265              0         0
      300       CHAN3         SEIZE                 232              0         0
      310        13           ASSIGN                232              0         0
      320        14           TRANSFER              232              0         0
      400       CHAN4         SEIZE                 187              0         0
      410        16           ASSIGN                187              0         0
      510       COME          DEPART               1001              0         0
      520        18           ADVANCE              1001              1         0
      530        19           RELEASE              1000              0         0
      540        20           TABULATE             1000              0         0
      550       WORK          TERMINATE            1000              0         0
      560       EXIT1         TERMINATE            1100              0         0
      570       EXIT2         TERMINATE           12125              0         0
    
    
    FACILITY   ENTRIES  UTIL.   AVE._TIME AVAILABLE  OWNER PEND INTER RETRY DELAY
     ZET1          317  0.776        1.02      1         0   0     0     1     0
     ZET2          265  0.675        1.06      1     14161   0     0     1     0
     ZET3          232  0.557        1.00      1         0   0     0     1     0
     ZET4          187  0.507        1.13      1         0   0     0     1     0
    
    
    QUEUE        MAX  CONT. ENTRIES ENTRIES(0) AVE.CONT. AVE.TIME  AVE.(-0) RETRY
      1            4     1    1002       765     0.72        0.30      1.27    0
    
    
    TABLE            MEAN      STD.DEV. RETRY  RANGE           FREQUENCY  CUM.%
     T_S              1.34        1.58    0 
                                              -         0         374      37.40
                                            0 -         2         448      82.20
                                            2 -         4         122      94.40
                                            4 -         6          42      98.60
                                            6 -         8          14     100.00
    
    
    XACT_GROUP         GROUP_SIZE       RETRY
     POSITION                  0           0

    По файлу стандартного отчета определим следующие операционные характеристики: среднее время пребывания в очереди и среднее время пребывания в модели (в системе). Эти характеристики обозначим $$Tq$$ и $$Ts$$ соответственно.

    Значение $$Tq$$ расположено в поле $$AVE.TIME$$ статистики очереди. Оно равно $$Tq = 0.30$$.

    Значение Ts есть среднее значение таблицы $$T_S$$, т. е. в поле $$MEAN$$. Оно равно $$Ts = 1.34$$.

    Поскольку очередь имеет нулевые входы – поле $$ENTRIES(0)$$, для расчета средней длины очереди $$Nq$$ используем аналог формулы Литтла [2], т. е.

    $$Nq = \lambda \cdot Tq = 3.52\cdot 0.30 = 1.0560.$$

    Для определения среднего количества требований $$Ns$$ в системе используем формулу Литтла:

    $$Ns = \lambda \cdot Ts = 3.52\cdot 1.34 = 4.7168.$$

    Задание 4

  • По данным файла стандартного отчета определите вероятность отказа, относительную пропускную способность, абсолютную пропускную способность.
  • Если требуется, то внесите изменения в программу для определения ве-роятности наличия очереди и вероятности загрузки всех каналов обслуживания.
  • Внесите изменения в программу для определения вероятностей состояния системы.
  • Полученные результаты сведите в таблицу, в которую также внесите операционные характеристики системы $$М/М/4/7$$ с ограниченным временем ожидания, полученные с помощью анализа в MATLAB.
  • 2.3. Пример моделирования системы типа М/М/M/K/M

    Система $$M/M/m/K/M$$ — это система с пуассоновским входящим пото-ком требований, с экспоненциальным законом обслуживания в m приборах, с допустимым числом требований в системе, не превышающим K, и с ограниченным числом источников нагрузки, которые создают поток из $$M$$ требований. Общее число K требований в системе заключено в интервале $$m \le K \le M$$, где $$M$$ — число требований, формируемых конечным числом источников нагрузки.

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

    $$\lambda_{k}= \left \begin{cases} \lambda(M-k), 0\le k \le K-1, \lambda=const;\\ 0\mbox{ }\mbox{ в остальных случаях}\\ \end{cases} \right.\\ \\ \\ \mu_{k}= \left \begin{cases} k\mu, 0\le k \le m, \mu=const;\\ m\mu, m\le k \le K.\\ \end{cases} \right.$$

    Диаграмма интенсивностей переходов для системы $$M/M/m/K/M$$ будет представлять собой конечный размеченный граф состояний, который показан на рис. 2.5.

    (рис 2.5) Граф состояний системы M/M/m/K/M

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

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

    Запишем следующее выражение для определения стационарных вероятностей состояний от $$0$$ до $$m – 1$$:

    $$p_{k}=p_{0}\rho^{k}C_{M}^{k},\mbox{ }0\le k\le m-1,$$

    где:

    $$C_{M}^{k}=\frac{M!}{k!(M-k)!}$$ — биномиальные коэффициенты;

    $$\rho = \lambda / \mu$$,

    $$p_{0}=\frac{1}{\sum\limits_{k=0}^{m-1} \rho^{k}C_{M}^{k}+\sum\limits_{k=m}^{K} \rho^{k}C_{M}^{k}\frac{k!m^{m-k}}{m!}} $$ — вероятность нулевого состояния.

    Расчет вероятностей состояний от $$m$$ до $$K$$:

    $$p_{m+r}=p_{0}\rho^{m+r}C_{M}^{m+r}\frac{(m+r)!}{m!m^{r}},\mbox{ }r=0,1,2,...,K-m.$$

    Если в последнем выражении сделать замену $$m+r=k$$, $$r=k-m$$, то получим

    $$p_{k}=p_{0}\rho^{k}C_{M}^{k}\frac{k!m^{m-k}}{m!},\mbox{ }m\le k\le K.$$

    Вероятность нулевого состояния была определена из нормировочного условия:

    $$\sum\limits_{k=0}^{K} p_{k}=1.$$

    Рассмотрим случай, когда система $$M/M/m/K/M$$ работает в режиме чистых потерь, т. е. когда параметры системы удовлетворяют условию $$M\ge K=m$$. Расчет вероятностей состояний будет определяться при $$K = m$$:

    $$p_{k}=\frac{\rho^{k}C_{M}^{k}}{\sum\limits_{i=0}^{m}\rho^{i}C_{M}^{i}},\mbox{ }k=0,1,2,...,m.$$

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

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

    $$\frac{dP_{0}}{dt}=-M\lambda P_{0}+\mu P_{1};\\ \frac{dP_{1}}{dt}=M\lambda P_{0}-[(M-1)\lambda +\mu] P_{1}+2\mu P_{2};\\ \frac{dP_{2}}{dt}=(M-1)\lambda P_{1}-[(M-2)\lambda +2\mu] P_{2}+3\mu P_{3};\\ .............................\\ \frac{dP_{m}}{dt}=[M-(m-1)]\lambda P_{m-1}-[(M-m)\lambda +m\mu] P_{m}+m\mu P_{m+1};\\ .............................\\ \frac{dP_{K-1}}{dt}=[M-(K-2)]\lambda P_{K-2}-\{[M-(KI-1)]\lambda +m\mu\} P_{K-1}+m\mu P_{K};\\ \frac{dP_{K}}{dt}=[M-(K-1)]\lambda P_{K-1}-m\mu P_{K}.$$

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

    $$P_{0}(0)=1,\qquad P_{j}(0)=0,\qquad j=1,2,...,K.$$

    Для стационарного режима рассмотрим ряд операционных характеристик в достаточно общем виде.

    Среднее число требований в системе:

    $$N_{ср}=\sum\limits_{k=0}^{K}kp_{k}.$$

    Среднее время пребывания одного требования в системе определим по формуле Литтла:

    $$T=\frac{N_{ср}}{\lambda}.$$

    Средняя длина очереди:

    $$N_{q}=\sum\limits_{k=m}^{K}kp_{k+m}.$$

    Среднее время пребывания требования в очереди определим по формуле Литтла:

    $$T_{q}=\frac{N_{q}}{\lambda}.$$

    Вероятность отказа в обслуживании соответствует вероятности того, что в системе находится $$K$$ требований (максимально допустимое число):

    $$p_{отк}=p_{K}=p_{0}\rho^{K}C_{M}^{K}\frac{K!m^{m-K}}{m!}.$$

    Относительная пропускная способность:

    $$Q=1-p_{отк}.$$

    Абсолютная пропускная способность:

    $$A=\lambda Q.$$

    Таким образом, операционные характеристики рассчитываются по известным стационарным вероятностям и параметрам системы.

    Практическая часть

    Пример 3. Для системы массового обслуживания $$М/М/4/7/9$$ постройте переходные вероятности состояний и определите для нее операционные характеристики при следующих параметрах: $$\lambda = 3.52 с^{–1}$$, $$\mu = 0.678 с^{–1}$$. Среда программирования — MATLAB.

    Запишем с буквенными коэффициентами дифференциальные уравнения для системы $$М/М/4/7/9$$:

    $$\frac{dP_{0}}{dt}=-M\lambda P_{0}+\mu P_{1},\\ \frac{dP_{1}}{dt}=M\lambda P_{0}-[(M-1)\lambda +\mu] P_{1}+2\mu P_{2},\\ \frac{dP_{2}}{dt}=(M-1)\lambda P_{1}-[(M-2)\lambda +2\mu] P_{2}+3\mu P_{3},\\ \frac{dP_{3}}{dt}=(M-2)\lambda P_{2}-[(M-3)\lambda +3\mu] P_{3}+4\mu P_{4},\\ \frac{dP_{4}}{dt}=(M-3)\lambda P_{3}-[(M-4)\lambda +4\mu] P_{4}+4\mu P_{5},\\ \frac{dP_{5}}{dt}=(M-4)\lambda P_{4}-[(M-5)\lambda +4\mu] P_{5}+4\mu P_{6},\\ \frac{dP_{6}}{dt}=(M-5)\lambda P_{5}-[(M-6)\lambda +4\mu] P_{6}+4\mu P_{7},\\ \frac{dP_{7}}{dt}=[M-(K-1)]\lambda P_{6}-4\mu P_{7},\\$$

    где:

    $$М = 9; K = 7; \lambda = 3.52 с^{–1}; \mu = 0.678 с^{–1}.$$

    Полученную систему представим в матричном виде:

    $$\frac{dP}{dt}=AP,$$

    где матрица $$А$$ размера $$8\times 8$$ составляется из коэффициентов при вероятностях в правой части уравнений.

    Для интегрирования дифференциальных уравнений примем естественные граничные условия.

    Программный код решения примера в MATLAB:

    function MMmKM;
    clc,close all
    % Параметры системы M/M/4/7/9
    L = 2.2;
    M = 3.678;
    m = 4;
    K = 7;
    N = 9; 
    % Матрица коэффициентов А
    global A
    A = [-N*L,M,zeros(1,6);
    N*L,-[(N-1)*L+M],2*M,zeros(1,5);
    zeros(1,1),(N-1)*L,-[(N-2)*L+2*M],3*M,zeros(1,4);
    zeros(1,2),(N-2)*L,-[(N-3)*L+3*M],4*M,zeros(1,3);
    zeros(1,3),(N-3)*L,-[(N-4)*L+4*M],4*M,zeros(1,2);
    zeros(1,4),(N-4)*L,-[(N-5)*L+4*M],4*M,zeros(1,1);
    zeros(1,5),(N-5)*L,-[(N-6)*L+4*M],4*M;
    zeros(1,6),(N-6)*L,-4*M];
     
    T = [0,2]; %% интервал интегрирования
    P0 = [1,zeros(1,length(A)-1)]; %% начальные условия
    [t,P] = ode23(@cmo,T,P0);
     
    %% Построение диаграмм вероятностей состояний
    % line(t,P,'linew',2) %% сплошные линии с различными цветами
    line(t,P(:,1),'linew',2, 'color','r') %% Po
    line(t,P(:,2), 'linew',2,'lines','--') %% P1
    line(t,P(:,3), 'linew',2,'lines','-.') %% P2
    line(t,P(:,4), 'linew',2,'lines',':') %% P3
    line(t,P(:,5), 'marker','o', 'color', 'm') %% P4
    line(t,P(:,6), 'marker','h', 'color','k') %% P5
    line(t,P(:,7), 'marker','p','color','r') %% P6
    line(t,P(:,8), 'marker','>') %% P7
    grid on
    Na = length(A) - 1;
    arr = [0:Na]';
    str = num2str(arr);
    legend(strcat('\bf\itP\rm\bf_', str, '(\itt\rm\bf)'));
    
    title(sprintf('%sСистема M/M/%d/%d/%d; %s%g; %s%g',...
    '\bf\fontsize{11}',m,K,N,'\lambda = ',L,'\mu = ',M));
    xlabel('\bf\it\fontsize{12} -  -  -  -  -  -  -  -   t   -  -  -  -  -  -  -  -');
    ylabel('\bf\fontsize{12}\itP\rm\bf(\itt\rm\bf)');
    set(gca, 'fontweight','bold', 'fontsize',10)
     
    Pcm = P(end,:); % Стационарные вероятности
    fprintf('\n Стационарные вероятности системы M/M/%d/%d/%d:\n', m,K,N);
    for J = 1 : length(A)
        fprintf('\tP%d = %f\n', J-1, Pcm(J));
    end
    fprintf('\n ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ СИСТЕМЫ M/M/%d/%d/%d:\n', m,K,N);
    Pnot = P(end,end);
    fprintf(' Вероятность отказа Pnot = %f\n', P(end,end));
    Q = 1 - Pnot;
    fprintf(' Относительная пропускная способность Q = %f\n', Q);
    Ab = L*Q;
    fprintf(' Абсолютная пропускная способность A = %f\n', Ab);
    Pq = sum(P(end, m+1:end));
    fprintf(' Вероятность наличия очереди Pq = %f\n', Pq);
    Ps = sum(P(end, m:end));
    fprintf(' Вероятность загрузки всех каналов обслуживания Ps = %f\n', Ps);
    Ns = [0:length(A)-1]*P(end,:)';
    fprintf(' Среднее количество требований в системе Ns = %f\n', Ns);
    fprintf(' Среднее время пребывания требования в системе Ts = %f\n', Ns/L);
    Nq = [0:(K-m)]*P(end,m:K)';
    fprintf(' Средняя длина очереди Nq = %f\n', Nq);
    fprintf(' Среднее время пребывания требования в очереди Tq = %f\n', Nq/L);
     
    function f = cmo(t,P)
    % М-функция описания правых частей дифференциальных уравнений:
    global A
    f = A*P;

    Результат выполнения программы

    Стационарные вероятности системы M/M/4/7/9:
    	P0 = 0.013187
    	P1 = 0.070956
    	P2 = 0.169841
    	P3 = 0.236851
    	P4 = 0.212652
    	P5 = 0.158847
    	P6 = 0.095048
    	P7 = 0.042619
    
     ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ СИСТЕМЫ M/M/4/7/9
     Вероятность отказа: Pnot = 0.042619
     Относительная пропускная способность: Q = 0.957381
     Абсолютная пропускная способность: A = 2.106238
     Вероятность наличия очереди: Pq = 0.509166
     Вероятность загрузки всех каналов обслуживания: Ps = 0.746017
     Среднее количество требований в системе: Ns = 3.634652
     Среднее время пребывания требования в системе: Ts = 1.652115
     Средняя длина очереди: Nq = 0.815489
     Среднее время пребывания требования в очереди: Tq = 0.370677

    Зависимости вероятностей состояний от времени показаны на рис. 2.6.

    (рис 2.6) Изменение вероятностей состояний во времени системы М/М/4/7/9

    Задание 5

  • Рассчитайте стационарные вероятности по аналитическим формулам (см. теоретическую часть). Сравните результаты.
  • Постройте зависимость среднего времени пребывания в очереди от интенсивности входного потока, изменяя интенсивность от $$0.Х$$ до $$5.Х$$, где $$Х$$ — номер компьютера (1, 2, 3, ...), за которым выполняется лабораторная работа.
  • Напишите программу для анализа системы с параметрами, вводимыми с клавиатуры пользователем, т. е. $$\lambda, \mu, m, K, М$$. Предусмотрите также ввод интервала интегрирования дифференциальных уравнений.
  • Контрольные вопросы

  • Что называется простейшим пуассоновским потоком?
  • По какому закону распределены интервалы времени между требованиями в простейшем пуассоновском потоке?
  • Что такое символика Кендалла, применяемая для систем массового обслуживания?
  • Что определяет собой точка на графике функции распределения случайной величины?
  • Какие случайные величины в системе массового обслуживания являются дискретными, а какими непрерывными?
  • Что означает несовместность событий?
  • Каким условиям должны отвечать начальные условия при решении дифференциальных уравнений относительно вероятностей состояния системы массового обслуживания?
  • Какова длина очереди в системе $$M/M/m/K/M$$?
  • Через какой тип данных могут выражаться интенсивность входного потока требований и интенсивность обслуживания требований в системах массового обслуживания?
  • Вернуться к учебному плану