CC BY
183
№5(23)2009
А. А. Емельянов
Стохастические сетевые модели массового обслуживания
Коммуникационные системы и сетевая организация сферы услуг стали причинами нового этапа развития двух направлений исследования операций — имитационного моделирования и теории массового обслуживания. В статье рассматриваются стохастические модели сетей массового обслуживания.
В реальной действительности существуют многочисленные объекты, в которых требования последовательно проходят несколько систем массового обслуживания (как однока-нальных, так и многоканальных), например:
• сети связи, компьютерные сети, телекоммуникационные системы;
• разветвленные или распределенные системы сферы услуг;
• образовательные системы, реализующие недетерминированные учебные планы (например, дополнительного образования или профессиональной подготовки в заочной или дистанционной форме, но с очной сессией).
Эти системы образуют более или менее сложную сеть, характеризуемую, с одной стороны, структурой, т. е. связями, существующими между системами, составляющими сеть, а с другой стороны, свойствами самих систем. Понятно, что аналитическое исследование подобной сети в общем виде чрезвычайно затруднительно. Однако некоторые типы сетей, как мы увидим в настоящей статье, оказываются удивительно простыми.
Рассмотрим сначала частный случай сети, для которого примем исходные предположения только в отношении составляющих систем; мы будем предполагать, что в этих системах могут образовываться неограниченные очереди, т. е. в каждой системе массового обслуживания, составляющей сеть, длительность ожидания и длина очереди не ограничиваются.
Линейные стохастические сети. Сети, которые мы исследуем в настоящей статье1, мо-
гут быть построены путем соединения конечного числа М систем массового обслуживания и внешнего источника заявок следующим образом. Требования, выходящие из системы / (/ = 1,2,...,М) с постоянной вероятностью 6Ц (у = 1,2,..., М), поступают в систему } (если }=/,2,...,М) или покидают сеть (у =0). С другой стороны, требования от внешнего источника поступают в сеть; вероятность того, что какое-нибудь из этих требований поступит в систему у (у = 1,2,..., М), равна d0у. Очевидно, должно иметь место равенство
£ dj¡, / = 0, 1,2,..., М,
у=0
(1)
причем d00 =0 (предположение d00 ^0 не представляет практического интереса).
Мы будем называть такую сеть линейной стохастической сетью; этот термин представляется оправданным, так как вероятность поступления требования в систему у в течение интервала (■,■ + dt) является линейной комбинацией с постоянными коэффициентами вероятностей выхода требований из разных систем сети d у.
В дальнейшем для упрощения изложения мы будем рассматривать внешний источник как нулевую систему сети; особенность этой системы состоит в том, что она содержит бесконечное число требований.
Структура сети может быть представлена графом, пример которого показан на рис. 1. Не следует смешивать такой граф, который мы будем называть графом передач сети,
1 Первые публикации по сетям этого вида с несколько иными обозначениями принадлежат А. Кофману и Р. Крюо-ну [1].
103
IS
is
«о
i
о
S3
о
«о §
«3 € а
5
о € QJ
№5(23)2009
с графом, представляющим матрицу перехода марковского процесса, который используется при изучении отдельных систем массового обслуживания.
Сеть
Внешний источник
Рис. 1. Структура сети в виде графа
Рассмотрим произвольное требование, которое последовательно проходит системы с номерами 1,2,...,I,...(/=0,1,2,...,М); из определения сети непосредственно вытекает, что последовательность таких переходов образует некоторую реализацию цепи Маркова, описывающей действие сети, со стохастической матрицей
0 d0i d02 • d0 M
di0 du di2 • d1M
D = d20 d2i d22 • d0 M
dM 0 dM1 dM2 ' dMM,
где выполняются условия (1) и справедливо следующее неравенство:
0 < dij< 1, /,} = 0,1,2,..., М.
Матрица Э называется матрицей передач. Понятия графа и матрицы сети поясняются двумя примерами на рис. 2 и рис. 3. Очереди
на рисунках имеют условные обозначения queue; узлы обслуживания — обозначения serve (система обозначений Actor Pilgrim). Сеть, показанная на рис. 2, в качестве узлов имеет одну систему массового обслуживания с неограниченным числом каналов (поэтому очередь не показана), и одну систему массового обслуживания с одним каналом.
Рис. 2. Сеть с двумя системами массового обслуживания Эта сеть имеет матрицу передач
D =
0 1 0' 0 0 1
0 1-0 0
где 0 — вероятность возврата заявки из сети в источник.
Соответствующий граф сети показан на рис. 3.
Рис. 3. Граф двухузловой сети
Более сложная сеть показана на рис. 4. Она имеет в качестве узлов три системы массового
104
№5(23)2009
обслуживания с очередями, у которых нет ограничений на длину. Причем системы 1 и 2 имеют по одному каналу обслуживания, а система 3 — два канала.
queue н ^serve^ i®j
Система 2
Рис. 4. Сеть с тремя системами массового обслуживания Матрица передач для этой сети имеет вид
D =
0 do, 0 do31
0 0 d,2 0
d20 d21 d22 d23
0 d3, d32 0
Интенсивности потоков заявок. Структура сети. Рассмотрим суммарную интенсивность потока Xу в системе у, т. е. среднее число заявок, поступающих в эту систему за единицу времени в установившемся режиме. Элементы Xу образуют вектор R. Среднее число требований, покидающих систему, очевидно, равно среднему числу поступающих требований; с другой стороны, предполагается, что вероятность того, что требование, покидающее систему /, направится в систему у, не зависит от предыдущего пути этого требования и состояния сети. Следовательно, если существует установившийся режим, то
Xj =ЕХ-'dj, j=0,1,2,...,M,
(2)
где X0 — суммарная интенсивность источника.
Далее используем матричные выражения и методы линейной алгебры. Введем в рассмотрение динамическую матрицу
В=О -1,
где I — матрица, в которой элементы главной диагонали равны 1, а остальные элементы равны нулю.
0
где do2 + doз =1; ^2 =1;
d2o + d2l + d22 + d2з = 1; dзl + dз1 = 1.
Соответствующий граф сети показан на рис. 5.
После этого получим уравнение В' R =0,
(3)
где В' — матрица, транспонированная по отношению к В;
R — вектор, имеющий М +1 элементов, удовлетворяющих условию (2).
В результате получена система М +1 линейных однородных уравнений (3) относительно X, (включая X0), (/'=0,1,2,...,М). Чтобы эта система имела нетривиальное решение R ^0, необходимо, чтобы ее определитель равнялся нулю:
Рис. 5. Граф трехузловой сети
i—1 d,0 d20 dM 0
d0, d,i -1 d2i dM1
B '| = d02 d,2 d22 -1 ' dM2 =0.
dM 0 dM1 dM2 ' dMM — 1
105
№5(23)2009
Это условие всегда удовлетворяется; если прибавить к какому-либо столбцу определителя сумму всех других столбцов, получится столбец, содержащий только нули, а это и показывает, что данный определитель равен нулю.
Итак, система (3) допускает бесконечное число решений. Однако, на самом деле, X0 известно заранее, поэтому X, ( >0) определяются единственным образом. Выражаясь более строго, это имеет место, если ранг системы в точности равен М.
Пример. Рассмотрим граф сети, имеющей три узла массового обслуживания. Значения интенсивностей потоков показаны на рис.6.
Рис. 6. Граф сети и значения интенсивностей потоков
Имеем
§ §
«о
i
о
о
«о §
«3 € а
5
о €
I
ш §
» а
i
.5
О =
0 1 0 0
1 1
0 0
2 2
0 1 0 3
4 4
1 0 1 1
2 4 4;
откуда
1 0
1 -
0 -
2 1
1 -
0 - - -2 4 4
Таким образом, система уравнений (3) может быть представлена в следующем виде:
-X 0 X 0
-X, +
2 X
IX2
4
-3 X2 42
+2 х.
1
Н—X2 4
- 3 Xз
4
0 0 0 0.
Решение этой системы: 9 8
^ = 7 0, ^ = 7 ^^ 0, Xз =2X 0.
(4)
Вернемся к нашему исследованию.
Чтобы определить, при каких условиях ранг матрицы В равен М, необходимо рассмотреть более тщательно структуру марковской цепи, определяемой матрицей О.
Цепь Маркова называется неразложимой [1], если каждое ее состояние может быть достигнуто из любого другого состояния посредством определенного числа переходов (здесь мы будем говорить, что требование находится в состоянии к, если оно находится в системе к). Если для процесса, который мы исследуем, указанное свойство не выполняется, это означает, что множество систем может быть разделено на р + 1 классов Б0, 5!, 52,..., Б,,..., Бр-1, Бр, таких, что выполняются следующие два условия.
1. При , = 0, 1, 2, ..., р - 1 , заявка, находящаяся в одном из состояний класса Б,, может попасть в другое состояние только из того же класса. Состояния класса Б, назовем возвратными.
2. При , = р заявка, которая находится в состоянии Бр, может покинуть Бр, но она не может туда возвратиться, так как вероятность того, что она не попадет в Бр после п переходов при п^га стремится к 1. Состояния класса Бр назовем невозвратными. Класс Бр может быть, в частности, пустым, т. е. может не существовать невозвратных состояний.
Предположим сначала, что система 0 (источник) является невозвратной, т. е. составляет часть класса Бр; тогда требования, поступающие из нее, неограниченно накапливаются в сети, не имея возможности покинуть последнюю; при этом установившийся режим
106
2
В
№5(23)2009
становится невозможным. Мы отбросим этот случай и будем предполагать, что система 0 является частью класса Б / (/^р), например 50, но тогда в установившемся режиме никакое требование не будет проходить через систему из класса Бр; даже если вначале этот класс содержал требования, они постепенно покинут его и не смогут возвратиться в прежнее состояние. Мы ограничимся исследованием установившегося режима; предположим, что невозвратные состояния в сети отсутствуют, т. е. все возможные состояния являются возвратными.
При этих условиях число требований, находящихся в классе Б/ , постоянно, и никакой переход в том или ином направлении из класса Б / в классе Б0 невозможен. Это соответствует системе массового обслуживания с ограниченным числом клиентов. Но если классы Б/ полностью изолированы друг от друга, они образуют независимые сети, которые могут быть исследованы отдельно.
Следовательно, мы можем без нарушения общности задачи предположить, что цепь Маркова является неразложимой.
Если X0 ^0, сеть будет называться разомкнутой.
Если X0 =0, сеть будет называться замкнутой; в этом случае она содержит постоянное число требований т, зависящее от начальных условий.
Перейдем к расчету интенсивностей потоков требований для разомкнутой и замкнутой сетей.
Потоки в разомкнутой сети (теорема без доказательства). Если цепь Маркова является неразложимой (как мы предположили), то ранг системы уравнений (3) равен М; следовательно, для каждого значения X0 существует единственное решение
X / =а / X 0, (5)
где а / — постоянные коэффициенты.
Пример был показан на рис. 6. Формула (4) подтверждает правильность выражения (5).
Потоки в замкнутой сети. Для замкнутой сети обязательно все d¡0 =0, '/ = 1,2,..., М, иначе в установившемся режиме сеть не содер-
жала бы ни одного требования, так как состояние 0 является невозвратным.
Матрица О в данном случае приводится к виду
О,
d11 d21
d12 d22
<1М М
<М1 с^31 ••• dMM,
а формула (2) преобразуется к виду
X
М
=Е X ,
(6)
где у = 1,2,.,М;
107
Определитель матрицы В С, транспонированной по отношению к Вс = Ос — I, всегда равен нулю, а величины X/ неопределенны,так как ни одна из них заранее не известна. В этом нет ничего удивительного, потому что число требований в сети зависит от начальных условий, которые не учитываются в случае разомкнутой сети. На том же основании, что и в случае открытой сети, ранг системы уравнений равен (М — 1). И в этом случае справедливы линейные соотношения вида
X/ =а(1 ^, (7)
для любых/,I = 0,1,2,.,М,/'^I,
где а (1) — постоянные коэффициенты, если выразить все X/ через заранее выбранную неизвестную интенсивность X1.
Однако решения эти соотношения не дают. Для получения решения необходимо учитывать следующее условие баланса: в любой момент времени суммарное число заявок во всех узлах сети равно общему числу заявок в сети, которое не меняется во времени. Это обстоятельство существенно усложняет практические расчеты, а во многих случаях они становятся невозможными. В таких случаях приходится прибегать либо к приближенным оценкам, либо к имитационному моделированию, которое дает довольно точные результаты (при наличии моделирующего пакета).
Максимальная интенсивность в разомкнутой сети. До сих пор мы не принимали в расчет ограничений, которые возникают, если учитывать свойства систем массового об-
№5(23)2009
о
о
«о §
«3 € а
5
о €
служивания, составляющих сеть. Для того чтобы существовал установившийся режим, каждая система сети должна удовлетворить условию ненасыщения X, <м, (X, <м, — суммарная интенсивность обслуживания системы ,, иными словами, среднее число требований, обслуживаемых в единицу времени, когда прибор (или приборы) непрерывно занят (заняты))2, т. е. ф, <1.
Следовательно, с учетом формулы (5), значение интенсивности внешнего источника должно удовлетворять неравенствам3
а,Xо <М,,
где , = 1,2,...,М, т. е. X0 < тт
I=1,2,., М
а,
Это неравенство является необходимым условием существования установившегося режима в разомкнутой сети.
Максимальная интенсивность в замкнутой сети. Как и в случае открытой сети, должно быть X, <м, (/ = 1,2,.,М).
Но здесь совокупность X,зависит от одного параметра (например, от X|). Предположим, что в уравнении I=1 (7), а величина X, возрастает, начиная с некоторого малого значения е. Для упрощения формул обозначим а(1) как а,. Тогда, начиная с некоторого значения X,, одно из вышеуказанных неравенств больше не будет удовлетворяться; это произойдет в системе ¡,
Мм
для которой отношение — минимально, при
а
а! =1.
Само собой разумеется, этот индекс , нисколько не зависит от выбора X! в качестве параметра. Мы можем, следовательно, предположить, что м 1 при к=2,3,.,М, т.е. что
а к
именно система с индексом 1 является узким местом сети; тогда в установившемся режиме будет выполняться неравенство X1 <м 1, каково бы ни было число требований в сети.
Заметим, что при заданном значении X1 отношение
I Xl ф1 =— М1
будет наибольшим значением коэффициента использования
ф = X,- = X1а¡ , М, М,
Величину ф1 назовем коэффициентом насыщения сети.
Отметим два важных момента.
Во-первых, очевидно, X1 будет тем больше, чем больше будет число требований в сети; если мы устремим это число к бесконечности, требования будут накапливаться в системе 1, для которой коэффициент использования ближе всего к 1.
Во-вторых, число требований в других системах не будет бесконечно увеличиваться, так как коэффициент использования этих систем имеет предел, меньший единицы. Тогда системы 2,3,.,М образуют в пределе разомкнутую сеть, по отношению к которой система 1 образует внешний источник.
Действительно, разомкнутая сеть, как мы ее определили, — это сеть, нагруженная стационарным источником; но система 1, в которой число содержащихся требований стремится к бесконечности, прибор (или приборы) которой непрерывно занят и интенсивность выходного потока требований которой равна М1, не зависит от остальной сети, и, в частности, от числа требований в системах 2,3,.,М. Заметим, что условие существования установившегося режима для этой разомкнутой сети удовлетворяется4, так как
М,
М1 <—.
а
Если число требований в сети конечно и равно М, величина коэффициента насыще-
§ »
а
2 Если система , содержит бесконечное число приборов, то м ; бесконечно, и условие ненасыщения удовлетворяется при любом X/.
3 Мы предполагаем здесь, что для всех систем интенсивность обслуживания не зависит от интенсивности потока.
4 Коэффициенты а, в этом случае аналогичны коэффициентам, определяемым выражением (5), так как система 1 здесь играет роль системы 0 в теории открытых сетей.
108
№5(23)2009
ния ф,, очевидно, зависит от Ы, но равным образом она зависит также и от свойств систем, составляющих сеть.
Время пребывания в разомкнутой сети.
Предположим, что нам известно для каждой системы массового обслуживания , среднее время пребывания и, (включающее продолжительность обслуживания) для установившегося режима сети. Пусть и — среднее время пребывания в сети, т. е. среднее значение промежутка времени, который проходит с момента поступления заявки в сеть (в очередь одной из ее систем), до момента, когда заявка покидает сеть (ту же или другую систему).
Рассмотрим снова цепь Маркова, описываемую матрицей О, и обратимся сначала к случаю открытой сети. Пусть ш, — вероятность состояния , в установившемся режиме, т. е. вероятность того, что рассматриваемая заявка находится в системе , после произвольного, но очень большого числа переходов. Вероятности ш, удовлетворяют соотношению
В' W=0,
(8)
где №=
ш, ш 2
ш м
м
Е ш,
Сравнение выражений (8) и (3) показывает, что ш, пропорциональны X,, т. е.
=а,.
ш, X,-ш 0 X 0
Но ш, может быть истолкована как относительная частота прохождений через состояние (т. е. через систему) ,; иначе говоря, если за достаточно большой промежуток времени Т заявка проходит Тш 0 раз через систему 0, она проходит в среднем Тш, раз через систему ¡. Среднее число прохождений через систему ,, соответствующее одному прохождению через систему 0, равно, следовательно,
Тш ш Тш 0 ш 0
=а,-.
Таким образом, требование, поступающее из источника, проходит в среднем а, раз через систему , перед тем как покинуть сеть; среднее время пребывания в системе равно, следовательно, а и , а среднее время пребывания в сети вычисляется как:
= Еа и,.
(9)
/(' ) =
= Еа ■ 1)и
(10)
где ), очевидно, должны быть вычислены для той же основной системы }, т. е.
,(1)
-X.
X,-'
и 1) =Еа|1 и, = — ЕХ,и ,
X,
(11)
V
-X¡d¡j,
(12)
Время пребывания в замкнутой сети.
Рассуждая аналогичным образом, можно вычислить среднее значение промежутка времени между моментом выхода требования из произвольной системы * и моментом первого возвращения того же требования в систему , для замкнутой сети. Пусть это среднее значение равно и('). Тогда получим
Если прибавить к и(') среднюю продолжительность и, пребывания в системе ,, то получится средняя продолжительность полного цикла с момента поступления требования в систему * до его следующего поступления в ту же систему
Как можно было ожидать, продолжительность цикла, измеренная на входе (или на выходе) различных систем сети, обратно пропорциональна интенсивности потока требований, поступающих в эти системы.
Измерение среднего потока в сети. Возвратимся к уравнениям (2) или (6) и рассмотрим сеть в установившемся режиме. На практике часто бывает трудно оценить интенсивности X, и коэффициенты передачи d¡j статистическими методами; зато можно измерить величины
109
м
0
№5(23)2009
представляющие средний поток требований, проходящих в единицу времени из системы , в систему ] (т. е. интенсивность поступления заявок в систему ] из системы ,).
Вводя эту величину в уравнения (2) и (6), получим
М
X, = Е . (13)
1 = 0
Из (12) и (13) выводим значение
а- = 1=1
и1 \ М ■
Xí М-,
1 = 0
Таким образом, зная поток Ц, нетрудно получить X, и а.
Покажем теперь, что в каждой системе средний поток Ц «консервативен», т. е. справедливо соотношение
(14)
§ §
«о
i
о
о
«о §
«3 € а
5
о €
I
ш §
» а
i
.5
ММ
Е-,1 = . ,'=0 к=0
На основании (13) имеем
М
Е Ц =X 1;
с другой стороны,
М М М
Е -1к = Е^^ 1 ^¡к ^ 1 Е а1к ^ 1 , к=0 к = 0 к = 0
что и подтверждает справедливость равенства (14). Если величины Ц получены в результате измерений, при помощи этого соотношения можно проверить, что система находится в установившемся режиме, и что в графе, изображающем сеть, не пропущена ни одна дуга. Физически свойство «консервативности» очевидно.
Линейные показательные сети. Если в разомкнутой сети внешний источник выдает в сеть пуассоновский поток заявок, а в качестве систем массового обслуживания выступают одно- или многоканальные системы массового обслуживания, с временем обслуживания, распределенным по экспоненциальному закону, и не имеющие ограничений на длину очереди, то формулы расчета среднего времени пребывания в сети становятся пригодными
для реальных расчетов: (9) — для разомкнутых сетей, (10) и (11) — для замкнутых.
С другой стороны, каждое значение и,
1
состоит из времени обслуживания ts¡ =—
М,
и длительности ожидания в очереди w¡ =<м-, (5,), являющейся функцией от числа каналов обслуживания Б,:
и
— + wi (Б,), М,
(15)
где Б, — интенсивность обслуживания в канале.
Ниже приведены соответствующие расчетные формулы для определения w¡ [2, 3]: 1) одноканальная система Б, =1
wi
М,(1-Р,)
-,, = 1,2,____ М,
(16)
где Р — загрузка системы:
, М,
2) многоканальная система Б, >1
Б,
W¡ =-—-~ Р0 (Б,),
М, • Б; • Б;!
Б,
(17)
где р0 (Б,) — вероятность того, что в системе , нет ни одной заявки, вычисляемая по формуле
Р0 (Б )=
Б Р, -+ 1+Р,+
Б,! 1- Р,^ 1
Б;
р2 р[ р Б-1 1 1
! 3! (Б, -1)!
В результате видно, что формулы (9), (10), (11), (15), (16) и (17) позволяют выполнять практические расчеты в линейных показательных сетях (при обязательном соблюдении условия баланса в случае замкнутой сети).
Замечание. Формула (17) содержит факториалы Б,! как в числителе, так и в знаменателе. Но современная вычислительная техника (64-разрядный компьютер) позволяет точно вычис-
110
0
№5(23)2009
лить только 20! (для 20 каналов в узле обслуживания). При увеличении числа каналов компьютер будет вычислять факториал только приближенно с использованием плавающей арифметики. Например, 100! (100 каналов в узле) будет содержать значительную погрешность, которая увеличится в разы, что приведет к ошибочным (иногда — к абсурдным результатам). Поэтому, когда S, >>20 для определения значения w¡ единственным инструментом является имитационная модель [4]. Пример фрагмента простейшей модели, позволяющей получать w¡ при произвольных распределениях интервала поступления показан ниже:
top(2): queue ("Очередь", none, 3) ; place; top(3): serve ("Обслуживание", 10 0 , none, tria, 5, 10, 15, 4); place.
Автоматически измеренное время пребывания в вершине top(2) этой модели — задержка в очереди, которая обслуживается 100-ка-нальным обслуживающим узлом top(3); причем интервал времени обслуживания распределен по треугольному закону (tria) с параметрами 5, 10 и 15.
Однако если прибегнуть к имитационному моделированию, то можно получить весьма точные результаты вообще для любых сетей без каких-либо математических расчетов: для этого нужна только мощная инструментальная система имитационного моделирования, например, Actor Pilgrim [4]. Сложность применения методов имитационного моделирования заключается в необходимости приобретения моделирующего пакета и в получении специальных знаний и подготовки для создания таких моделей.
В качестве заключения: единственность установившегося режима.
1. Если определены законы, которым подчинен ход процесса массового обслуживания, то установившийся режим будет зависеть только от начальных условий. Эта зависимость может быть исследована на основе обратных уравнений [1] рассматриваемого процесса. В статье не приведены соответствующие выкладки, так как они потребовали бы введения громоздких обозначений и преоб-
разований. Следует отметить, что структура сети влияет только на расчет интенсивности потоков. Если последние вычислены, параметры установившегося режима определяются только в зависимости от свойств составляющих систем.
2. Классические методы теории массового обслуживания базируются на общих результатах теории вероятностей: марковские и полумарковские процессы, вложенные цепи Маркова (метод Кендалла). Они имеют два существенных практических достоинства:
• возможность исследования очередей с приоритетами, причем источники заявок могут быть как конечными, так и бесконечными (широко применяется для анализа технических автоматизированных систем);
• возможность исследования линейных стохастических сетей (применяется при проектировании сетей связи и телекоммуникационных систем широкого назначения, в которых из-за их сложности начинают работать «законы больших чисел», в результате чего гипотеза о линейности сети начинает выполняться).
3. Принципиальный недостаток классической теории массового обслуживания: нет возможности анализа систем в случаях, когда потоки заявок не являются простейшими. Преодоление этого недостатка оказалось возможным с применением методов, известных под названием «диффузная аппроксимация». При их разработке сделан решительный отход от классического марковского аппарата. Существуют точные математические приемы описания такого непрерывного процесса в виде дифференциальных уравнений, имеющих аналитические решения. Рассмотрение этих методов предполагается в последующих статьях.
СПИСОК ЛИТЕРАТУРЫ
1. КофманА., Крюон Р. Массовое обслуживание. Теория и приложения. М.: Мир, 1965.
2. КофманА. Методы и модели исследования операций. М.: Мир, 1966.
3. Емельянов А. А. Модели процессов массового обслуживания / Прикладная информатика. 2008. №5(17).
4. Емельянов А. А., Власова Е. А., Дума Р. В, Емельянова Н.З. Компьютерная имитация экономических процессов: учебник. М.: Маркет ДС, 2009.
111
