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

function twophase;
clc,close all
% ВХОДНЫЕ ДАННЫЕ
L = 1.7;
M = 0.67;
global A
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;
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
	 Относительная пропускная способность системы: 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