Моделирование систем. Моделирование систем 2е издание, исправленное
 Скачать 5.72 Mb.
|
|
Моделирование систем 2-е издание, исправленное Афонин В.В. Федосин С.А. Национальный Открытый Университет "ИНТУИТ" 2016 УДК 004.94(076.5) ББК 16 А94 Моделирование систем / Афонин В.В., Федосин С.А. - M.: Национальный Открытый Университет "ИНТУИТ", 2016 (Основы информационных технологий) ISBN 978-5-9963-0352-6 Курс состоит из теоретической и практической частей. В каждой работе приводятся необходимые теоретические сведения для выполнения и понимания практических примеров. В практической части рассматриваются характерные примеры и их полная программная реализация в таких средах, как MATLAB, GPSS/PC. Большая часть программного кода реализована в MATLAB. Данный курс охватывает следующие разделы: системы массового обслуживания (многофазные и многоканальные системы обслуживания), статистическое моделирование на ЭВМ (включая метод Монте-Карло, интервальные методы оценки параметров вероятностных законов и метод максимального правдоподобия), дисперсионный анализ, планирование машинных экспериментов с моделями систем (пассивный и активный эксперименты), регрессионный анализ в случае неполного ранга наблюдений и идентификации линейных непрерывных систем управления. (c) ООО "ИНТУИТ.РУ", 2010-2016 (c) Афонин В.В., Федосин С.А., 2010-2016 Моделирование многофазных систем массового обслуживания Цель работы: практически освоить методы моделирования двухфазных и трехфазных систем массового обслуживания с нулевой вместимостью блоков ожидания в программных средах MATLAB и GPSS/PC с целью получения операционных характеристик. Теоретическая часть Анализ многофазных систем массового обслуживания основан на теоретическом материале, взятом из [20] Система массового обслуживания может представляться в виде много- фазной модели, когда каждое требование в ней последовательно обслуживается во всех фазах (приборах обслуживания). При этом если очереди перед каждой фазой не допускаются, то система будет называться системой с нулевой вместимостью блоков ожидания [20] Работа двухфазной системы обслуживания состоит в следующем. Каждая фаза может быть занята на обслуживание либо свободна. Поскольку перед фазой очередь не допускается, то принимается, что первая фаза обслуживания заблокирована, если обслуживание требования в данной (первой) фазе завершено, а вторая фаза не готова к приему требования по той причине, что в ней не закончено обслуживание. Принимается также, что если первая фаза занята, то очередное входящее требование получает отказ. В системе могут быть следующие состояния: "фаза свободна", "фаза занята", "фаза заблокирована", которые обозначают как , , соответственно. 1. Двухфазная система обслуживания Схема двухфазной системы обслуживания показана на рис. 1.1 Рис. 1.1. Модель двухфазной системы обслуживания Если состояние первой фазы обозначить символом , а состояние второй фазы — символом , то множество состояний двухфазной системы обслуживания будет следующим: (1.1) Примем, что входной поток требований имеет пуассоновское распределение, а обслуживание в фазах осуществляется по экспоненциальному вероятностному закону. Рассматривая вероятности переходов из одного состояния в другое во времени, можно получить следующие дифференциальные уравнения относительно вероятностей состояний двухфазной системы: (1.2) Система (1.2) — это система однородных обыкновенных дифференциальных уравнений с постоянными коэффициентами. Ее можно представить в матричном виде: где: — вектор размера с элементами , , , , ; — матрица коэффициентов размера , которая имеет следующий вид: Для интегрирования системы (1.2) будем задавать естественные начальные условия, т. е. такие, когда в начальный момент времени, равный нулю, вероятность отсутствия требований в системе равна единице, а остальные вероятности в начальный момент времени равны нулю: (1.3) 2. Трехфазная система обслуживания В трехфазной системе каждая из фаз может быть свободной (символ ) либо занятой (символ ), а фазы 1 и 2 могут быть к тому же заблокированы (символ ). Если состояние первой фазы обозначить символом , состояние второй фазы — символом , а состояние третьей фазы — символом , то возможные состояния трехфазной системы будут следующими: (1.4) В соответствии с возможными состояниями (1.4) трехфазной системы можно получить следующую систему дифференциальных уравнений 13- го порядка относительно вероятностей состояний : (1.5) Система (1.5) — это система обыкновенных однородных дифференциальных уравнений с постоянными коэффициентами. Ее также можно представить в матричном виде: где: — вектор состояний системы 13-го порядка; — матрица постоянных коэффициентов размера Для решения системы (1.5) будем использовать естественные начальные условия, т. е. (1.6) Отметим также, что система линейных дифференциальных уравнений будет иметь нетривиальное решение, если ее характеристический определитель будет равен нулю. В матричном виде получаем так называемое характеристическое уравнение (1.7) где: — единичная матрица того же размера, что и матрица коэффициентов ; — скалярная в общем случае комплексная переменная, относительно которой решается характеристическое уравнение. Если действительная часть корней характеристического уравнения (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 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 Переходные процессы вероятностей состояний показаны на рис. 1.3 Рис. 1.3. Переходные процессы вероятностей состояний в двухфазной системе Задание 1 1. Вычислите по формуле Литтла среднее количество требований в системе. 2. Проверьте нормировочное условие для вероятностей состояний системы в трех точках временного интервала интегрирования. 3. В соответствии с номером компьютера, на котором выполняется лабораторная работа (1, 2, 3 и т. д.), увеличьте интенсивность входного потока на этот номер. Выполните два предыдущих пункта задания. 4. Результаты расчета операционных характеристик запишите в текстовый файл с именем compX.txt, где — номер компьютера, на котором выполняется лабораторная работа. 5. Напишите программу с различными интенсивностями обслуживания в фазах обслуживания. Рассчитайте также операционные характеристики двухфазной системы обслуживания. 2. Пример моделирования в системе GPSS/PC Для моделирования в системе GPSS/PC необходимо подготовить данные по функциям распределения случайных величин для имитации пуассоновского входного потока и экспоненциального обслуживания. В обоих случаях используется экспоненциальное распределение, поскольку интервалы времени между требованиями в пуассоновском потоке распределены по экспоненциальному закону. В рассматриваемом случае параметры пуассоновского потока и экспоненциального обслуживания заданы и равны соответственно , Для подготовки данных по функциям распределения случайных величин использован следующий программный код: clear all,clc L = 1.7; %% Интенсивность входного потока M = 0.67; %% Интенсивность обслуживания %% Функции экспоненциального распределения x = 0 : 0.2 : 20; F = 1 - exp(-L*x); F2 = 1 - exp(-M*x); %%% Запись в текстовый файл input.txt fid = fopen('input.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 fclose(fid); В программе функция применена для формирования заданного количества массива строк с данными для выбора случайных чисел, распределенных по экспоненциальному закону. Сформированные массивы будут представлены в GPSS-программе. GPSS-программа приводится ниже. В ней задано условие обработки 500 требований. Программный код решения примера в системе GPSS/PC: simulate puas function RN5,C50 0,0/0.28823,0.2/0.493383,0.4/0.639405,0.6/ 0.743339,0.8/0.817316,1/0.869971,1.2/0.907449,1.4/0.934125,1.6/ 0.953112,1.8/0.966627,2/0.976246,2.2/0.983093,2.4/0.987966,2.6/ 0.991434,2.8/0.993903,3/0.995661,3.2/0.996911,3.4/0.997802,3.6/ 0.998435,3.8/0.998886,4/0.999207,4.2/0.999436,4.4/0.999598,4.6/ 0.999714,4.8/0.999797,5/0.999855,5.2/0.999897,5.4/0.999927,5.6/ 0.999948,5.8/0.999963,6/0.999974,6.2/0.999981,6.4/0.999987,6.6/ 0.99999,6.8/0.999993,7/0.999995,7.2/0.999997,7.4/0.999998,7.6/ 0.999998,7.8/0.999999,8/0.999999,8.2/0.999999,8.4/1,8.6/ 1,8.8/1,9/1,9.2/1,9.4/1,9.6/ 1,9.8/ ;------------------------------ expM function RN6,C50 0,0/0.12541,0.2/0.235092,0.4/0.331019,0.6/ 0.414916,0.8/0.488291,1/0.552465,1.2/0.60859,1.4/0.657677,1.6/ 0.700608,1.8/0.738154,2/0.770992,2.2/0.799712,2.4/0.82483,2.6/ 0.846798,2.8/0.866011,3/0.882815,3.2/0.897511,3.4/0.910364,3.6/ 0.921605,3.8/0.931437,4/0.940035,4.2/0.947556,4.4/0.954133,4.6/ 0.959885,4.8/0.964916,5/0.969316,5.2/0.973164,5.4/0.976529,5.6/ 0.979473,5.8/0.982047,6/0.984299,6.2/0.986268,6.4/0.98799,6.6/ 0.989496,6.8/0.990813,7/0.991965,7.2/0.992973,7.4/0.993854,7.6/ 0.994625,7.8/0.995299,8/0.995889,8.2/0.996404,8.4/0.996855,8.6/ 0.99725,8.8/0.997595,9/0.997896,9.2/0.99816,9.4/0.998391,9.6/ 0.998593,9.8/ tab1 table m1,0,5,12 tab2 table mp2,0,5,10 ******** Basic Program ******* 5 generate fn$puas 10 savevalue FULL+,1 ;--- Probability P00 ----------- 15 test E f1,0,met1 20 test E f2,0,met1 25 savevalue P_00+,1 ;--- Probability P01 ----------- 30 met1 test E f1,0,met2 35 savevalue P_01+,1 ;--- Probability P10 ----------- 40 met2 test E f2,0,met3; 45 savevalue P_10+,1 ;--- Probability P11+Pb1--------- 50 met3 test NE f1,0,met4 55 test NE f2,0,met4 60 test NE q2,0,met4 65 savevalue P_11b1+,1 70 met4 assign 2,1 75 gate NU 1,exitNOT ;---------- 1st phase --------------- 80 seize 1 85 advance fn$expM 90 release 1 95 mark 2 ;---------- Delay ------------------ 100 queue 2 105 gate NU 2 110 depart 2 120 tabulate tab2 ;---------- 2nd phase --------------- 125 seize 2 130 advance fn$expM 135 release 2 ;------------------------------------ 140 tabulate tab1 145 terminate 1 150 exitNOT savevalue exit0+,1 155 terminate start 500 ;end Комментарии к программе Введение таблицы с аргументом (системный стандартный числовой атрибут) — временем пребывания в модели транзакта, обрабатываемого программой в данный момент, позволяет табулировать величину и найти среднее значение времени пребывания транзакта в модели. Введение таблицы с аргументом (стандартный числовой атрибут транзакта) — значением времени, равного разности относительного модельного времени и содержимого 2-го параметра текущего транзакта. Табулируется время пребывания транзакта в заблокированном состоянии. Отметка времени осуществляется блоком В программе использованы стандартные числовые атрибуты , , которые определяют состояния соответствующих устройств. Если устройство свободно, то величина стандартного числового атрибута равна 0, и 1 — во всех остальных случаях. Результат выполнения программы система GPSS/PC оформляет в виде файла стандартного отчета, который имеет следующий вид: GPSS/PC Report file REPORT.GPS. (V 2, # 37349) 01-08-2010 10:28:19 page 1 START_TIME END_TIME BLOCKS FACILITIES STORAGES FREE_MEMORY 0 675 30 2 0 135296 LINE LOC BLOCK_TYPE ENTRY_COUNT CURRENT_COUNT RETRY 5 1 GENERATE 3092 0 0 10 2 SAVEVALUE 3092 0 0 15 3 TEST 3092 0 0 20 4 TEST 501 0 0 25 5 SAVEVALUE 159 0 0 30 MET1 TEST 3092 0 0 35 7 SAVEVALUE 501 0 0 40 MET2 TEST 3092 0 0 45 9 SAVEVALUE 982 0 0 50 MET3 TEST 3092 0 0 55 11 TEST 2591 0 0 60 12 TEST 1768 0 0 65 13 SAVEVALUE 1453 0 0 70 MET4 ASSIGN 3092 0 0 75 15 GATE 3092 0 0 80 16 SEIZE 501 0 0 85 17 ADVANCE 501 1 0 90 18 RELEASE 500 0 0 95 19 MARK 500 0 0 100 20 QUEUE 500 0 0 105 21 GATE 500 0 0 110 22 DEPART 500 0 0 120 23 TABULATE 500 0 0 125 24 SEIZE 500 0 0 130 25 ADVANCE 500 0 0 135 26 RELEASE 500 0 0 140 27 TABULATE 500 0 0 145 28 TERMINATE 500 0 0 150 EXITNOT SAVEVALUE 2591 0 0 155 30 TERMINATE 2591 0 0 FACILITY ENTRIES UTIL. AVE._TIME AVAILABLE OWNER PEND INTER RETRY DELAY 1 501 0.800 1.08 1 3090 0 0 0 0 2 500 0.740 1.00 1 0 0 0 0 0 QUEUE MAX CONT. ENTRIES ENTRIES(0) AVE.CONT. AVE.TIME AVE.(-0) RETRY 2 19 0 500 207 2.95 3.98 6.79 0 TABLE MEAN STD.DEV. RETRY RANGE FREQUENCY CUM.% TAB1 6.06 5.22 0 - 0 50 10.00 0 - 5 233 56.60 5 - 10 120 80.60 10 - 15 67 94.00 15 - 20 24 98.80 20 - 25 6 100.00 TAB2 3.98 5.03 0 - 0 207 41.40 0 - 5 137 68.80 5 - 10 97 88.20 10 - 15 40 96.20 15 - 20 18 99.80 20 - 25 1 100.00 XACT_GROUP GROUP_SIZE RETRY POSITION 0 0 GPSS/PC Report file REPORT.GPS. (V 2, # 37349) 01-08-2010 10:28:19 page 2 SAVEVALUE VALUE RETRY FULL +3092 0 P_00 +159 0 P_01 +501 0 P_10 +982 0 P_11B1 +1453 0 EXIT0 +2591 0 Расчет операционных характеристик двухфазной системы по файлу стандартного отчета будем на окончании помечать буквой . Расчет вероятности отказа в обслуживании определим как отношение числа требований, получивших отказ (число 2591 в сохраняемой ячейке ), к общему числу (число 3092 в сохраняемой ячейке ): Относительная пропускная способность и вероятность того, что вновь поступившее требование будет принято на обслуживание 1-й фазой, равны между собой и определяются как противоположная вероятность отказа в обслуживании, т. е. Расчет эффективной частоты поступления требований в систему: Абсолютная пропускная способность двухфазной системы: Полное время пребывания требования в системе определяется как среднее значение таблицы с именем : Среднее время блокировки 1-й фазы определяется как среднее значение таблицы с именем : Среднее время обработки одного требования в двухфазной системе можно определить двумя путями: как сумму значений полей для устройств 1 и 2 ( ) и как разность между и : Среднее число требований в двухфазной системе определим по формуле Литтла: Финальные (стационарные) вероятности определим как частное от деления числа требований, отвечающим заданным условиям (успешное событие), к общему числу требований, поступивших в систему. Для этого в программе были созданы ячейки сохраняемых величин: для вероятности — ячейка , для вероятности — ячейка , для вероятности — ячейка , для суммы вероятностей (поскольку они равны между собой) — ячейка . В соответствии с результатами, приведенными в файле стандартного отчета, получим ; ; ; Проведем контрольную проверку. Для этого сложим вероятности (которые являются несовместными): Поскольку полученный результат практически равен единице, можно считать, что расчет вероятностей произведен успешно. Задание 2 1. Выполните прогон программы с датчиками случайных чисел в соответствии с номером компьютера, за которым выполняется лабораторная работа, т. е. , и т. д. Сравните полученные операционные характеристики с соответствующими характеристиками, определенными при моделировании двухфазной системы в MATLAB. 2. Видоизмените программу так, чтобы можно было построить функцию распределения времени блокировки 1-й фазы. 3. Напишите программу с различными интенсивностями обслуживания в фазах обслуживания. Рассчитайте также операционные характеристики двухфазной системы обслуживания. Сравните с операционными характеристиками, полученными при моделировании в системе MATLAB. 4. Смоделируйте в MATLAB и в GPSS/PC (GPSS World) трехфазную систему обслуживания и рассчитайте операционные характеристики по параметрам примера для двухфазной системы, т. е. , Контрольные вопросы |