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

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

2-е издание, исправленное

Афонин В.В.
Федосин С.А.

Национальный Открытый Университет “ИНТУИТ”
2016

2

УДК 004.94(076.5)
ББК 16
А94
Моделирование систем / Афонин В.В., Федосин С.А. - M.: Национальный Открытый Университет
“ИНТУИТ”, 2016 (Основы информационных технологий)
ISBN 978-5-9963-0352-6

Курс состоит из теоретической и практической частей. В каждой работе приводятся необходимые
теоретические сведения для выполнения и понимания практических примеров. В практической части
рассматриваются характерные примеры и их полная программная реализация в таких средах, как
MATLAB, GPSS/PC. Большая часть программного кода реализована в MATLAB.
Данный курс охватывает следующие разделы: системы массового обслуживания (многофазные и
многоканальные системы обслуживания), статистическое моделирование на ЭВМ (включая метод
Монте-Карло, интервальные методы оценки параметров вероятностных законов и метод
максимального правдоподобия), дисперсионный анализ, планирование машинных экспериментов с
моделями систем (пассивный и активный эксперименты), регрессионный анализ в случае неполного
ранга наблюдений и идентификации линейных непрерывных систем управления.

(c) ООО “ИНТУИТ.РУ”, 2010-2016
(c) Афонин В.В., Федосин С.А., 2010-2016

3

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

Цель работы: практически освоить методы моделирования двухфазных и трехфазных
систем массового обслуживания с нулевой вместимостью блоков ожидания в
программных средах MATLAB и GPSS/PC с целью получения операционных
характеристик.

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

Анализ многофазных систем массового обслуживания основан на теоретическом
материале, взятом из [20].

Система массового обслуживания может представляться в виде много-фазной модели,
когда каждое требование в ней последовательно обслуживается во всех фазах
(приборах обслуживания). При этом если очереди перед каждой фазой не допускаются,
то система будет называться системой с нулевой вместимостью блоков ожидания [20].

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

1. Двухфазная система обслуживания

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

Рис. 1.1.  Модель двухфазной системы обслуживания

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

(1.1)

Примем, что входной поток требований имеет пуассоновское распределение, а

4

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

(1.2)

Система (1.2) — это система однородных обыкновенных дифференциальных уравнений
с постоянными коэффициентами. Ее можно представить в матричном виде:

где:

 — вектор размера 
 с элементами 
, 
, 
, 
, 
 ;

 — матрица коэффициентов размера 
, которая имеет следующий вид:

Для интегрирования системы (1.2) будем задавать естественные начальные условия, т.
е. такие, когда в начальный момент времени, равный нулю, вероятность отсутствия
требований в системе равна единице, а остальные вероятности в начальный момент
времени равны нулю:

(1.3)

2. Трехфазная система обслуживания

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

(1.4)

5

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

(1.5)

Система (1.5) — это система обыкновенных однородных дифференциальных
уравнений с постоянными коэффициентами. Ее также можно представить в матричном
виде:

где:

 — вектор состояний системы 13-го порядка;

 — матрица постоянных коэффициентов размера 
.

Для решения системы (1.5) будем использовать естественные начальные условия, т. е.

(1.6)

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

(1.7)

где:

 — единичная матрица того же размера, что и матрица коэффициентов  ;

 — скалярная в общем случае комплексная переменная, относительно которой

решается характеристическое уравнение.

6

Если действительная часть корней характеристического уравнения (1.7) будет
отрицательной, то решение системы дифференциальных уравнений (1.5) с начальными
условиями (1.6) будет устойчивым, т. е. будет стремиться к установившимся
значениям.

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

Схема трехфазной системы обслуживания показана на рис. 1.2.

Рис. 1.2.  Модель трехфазной системы обслуживания

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

Рассмотрим моделирование двухфазной системы массового обслуживания в
программных средах MATLAB и GPSS/PC. Сначала выполним моделирование
двухфазной системы в MATLAB. Примем, что параметр входного простейшего потока
равен 
, а интенсивность обслуживания каждой фазы равна 
, т. е. 
, 

. В основе моделирования стоит задача интегрирования системы

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

1. Пример моделирования в системе matlab

Программа моделирования оформлена в виде 
 -функции без аргументов и без

возвращаемого значения. Это позволит в одном 
 -файле использовать

дополнительные 
 -функции.

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

function twophase;
clc,close all
% ВХОДНЫЕ ДАННЫЕ
L = 1.7;
M = 0.67;
global A

7

A = [-L, M, 0, 0, 0;
    0, -(L+M), M, 0, M;
    L, 0, -M, M, 0;
    0, L, 0, -2*M, 0;
    0, 0, 0, M, -M];

% ИНТЕРВАЛ ИНТЕГРИРОВАНИЯ
T = [0, 20];
% НАЧАЛЬНЫЕ УСЛОВИЯ
P0 = [1;zeros(length(A)-1,1)];
% ВЫЗОВ РЕШАТЕЛЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
[t,P] = ode45(@faza,T,P0);

% ПОСТРОЕНИЕ ДИАГРАММ ВЕРОЯТНОСТЕЙ СОСТОЯНИЙ
figure(1);
line(t,P(:,1),'linew',2, 'color','r'),
line(t,P(:,2), 'marker','o','color', [0,102,102]/255),
line(t,P(:,3), 'lines','--','linew',2, 'color', 'k'),
line(t,P(:,4), 'lines','-.','linew',2, 'color', 'b'),
line(t,P(:,5), 'linestyle',':','linew',2, 'color','m'),

 
str='\bfВероятности состояний двухфазной системы обслуживания:';
 
title(sprintf('%s %s = %g; %s = %g',str,'\lambda',L,'\mu', M));
legend('P_0_0(t)','P_0_1(t)','P_1_0(t)','P_1_1(t)','P_b_1(t)');
xlabel('\bf -  -  -  -  -  -  t -  -  -  -  -  - ');
ylabel('\bf P(t)');
%%% ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ СИСТЕМЫ
% Ncp - СРЕДНЕЕ ЧИСЛО ТРЕБОВАНИЙ В СИСТЕМЕ
k =[0, 1, 2];
Pk = [P(end,1);P(end,2)+P(end,3);P(end,4)+P(end,5)];
Ncp = k*Pk;

fprintf('\n\t ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ ДВУХФАЗНОЙ СИСТЕ-МЫ:\n');
fprintf('\t Среднее число требований в системе: Ncp = %f\n',Ncp);

% Pf1 - ВЕРОЯТНОСТЬ НАЧАЛА ОБСЛУЖИВАНИЯ 
Pf1 = P(end,1)+P(end,2);
fprintf('\t Вероятность начала обслуживания: Pf1 = %f\n',Pf1);

% Pomk - ВЕРОЯТНОСТЬ ОТКАЗА В ОБСЛУЖИВАНИИ  
Pomk = 1-Pf1;
fprintf('\t Вероятность отказа: Pomk = %f\n',Pomk);

%  Q - ОТНОСИТЕЛЬНАЯ ПРОПУСКНАЯ СПОСОБНОСТЬ 
Q = 1-Pomk;

8

fprintf('\t Относительная пропускная способность системы: Q = %f\n',Q);

% Lef - ЭФФЕКТИВНАЯ ЧАСТОТА ПОСТУПЛЕНИЯ ТРЕБОВАНИЙ В СИСТЕМУ
Lef = L*Pf1;
fprintf('\t Эффективная частота поступления требований в систему: Lef = %f\n',Lef);

% Ab - АБСОЛЮТНАЯ ПРОПУСКНАЯ СПОСОБНОСТЬ 
Ab = Lef*Q;
fprintf('\t Абсолютная пропускная способность системы: Ab = %f\n',Ab)

% Ts - ПОЛНОЕ ВРЕМЯ ПРЕБЫВАНИЯ ТРЕБОВАНИЙ В СИСТЕМЕ  
Ts = Ncp/Lef;
fprintf('\t Полное время пребывания требования в системе: Ts = %f\n',Ts);
% Tcp1 - СРЕДНЕЕ ВРЕМЯ ОБСЛУЖИВАНИЯ ОДНОГО ТРЕБОВАНИЯ
Tcp1 = 2*(1/M);
fprintf('\t Среднее время обслуживания требования в системе: Tcp = %f\n',Tcp1);
% Td - СРЕДНЕЕ ВРЕМЯ ЗАДЕРЖКИ 
Td = Ts - Tcp1;
fprintf('\t Среднее время задержки в обслуживании: Td = %f\n',Td);

% Финальные (стационарные) вероятности
fprintf('\n\t Финальные (стационарные) вероятности:\n');
k = 0;
for J = 1 : length(P0)
    if J == 1 | J == 2
fprintf('\t P0%d = %f\n ',J-1,  P(end, J));
    
    elseif J > 2 & J < length(P0)
fprintf('\t P1%d = %f\n ', k,  P(end, J));
k = k + 1;
   
    else
fprintf('\t Pb1 = %f\n', P(end,end));
   end
end
 

% М-функция описания правых частей дифференциальных уравнений
function f = faza(t,P)
global A
f = A*P;

Результаты моделирования по определению операционных характеристик

ОПЕРАЦИОННЫЕ ХАРАКТЕРИСТИКИ ДВУХФАЗНОЙ СИСТЕМЫ
 
 Среднее число требований в системе: Ncp = 1.345670

 
 Вероятность начала обслуживания: Pf1 = 0.224859

 
 Вероятность отказа: Pomk = 0.775141

9

 
 Относительная пропускная способность системы: Q = 0.224859

 
 Эффективная частота поступления требований в систему: Lef = 0.382260

 
 Абсолютная пропускная способность системы: Ab = 0.085955

 
 Полное время пребывания требования в системе: Ts = 3.520298

 
 Среднее время обслуживания требования в системе: Tcp = 2.985075

 
 Среднее время задержки в обслуживании: Td = 0.535224

 
 

          Финальные (стационарные) вероятности:
 
 P00 = 0.063570

  
 P01 = 0.161289

  
 P10 = 0.365900

  
 P11 = 0.204626

  
 Pb1 = 0.204615

Переходные процессы вероятностей состояний показаны на рис. 1.3.

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

Задание 1

1. Вычислите по формуле Литтла среднее количество требований в системе.
2. Проверьте нормировочное условие для вероятностей состояний системы в трех

точках временного интервала интегрирования.

3. В соответствии с номером компьютера, на котором выполняется лабораторная

работа (1, 2, 3 и т. д.), увеличьте интенсивность входного потока на этот номер.
Выполните два предыдущих пункта задания.

10