Портал GPSS.RU

Ю. И. Рыжиков

ОПЫТ ПОНИЖЕНИЯ ДИСПЕРСИИ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ


 

 

1.         Введение


Имитационное моделирование является важнейшим инструментом поисковых исследований в области анализа сложных систем и верификации аналитических и чис-ленных расчетных методик. Известный недостаток его – статистическая погрешность результатов – в условиях резко возросшего быстродействия ЭВМ перестал быть крити-чески важным. Однако остаются задачи, в которых проблема понижения дисперсии ре-зультатов моделирования остается актуальной:
1) анализ сильно загруженных систем массового обслуживания (СМО);
2) расчет структурно сложных систем с сильно зависящими от состояния вы-ходными показателями (например, те же СМО с отказами каналов);
3) расчет вероятностей редких событий (число испытаний, требуемое для полу-чения оценки с заданной относительной точностью, обратно пропорционально искомой вероятности и квадрату допустимой ее погрешности).
В таких случаях приходится прибегать к методам понижения дисперсии (МПД) результатов единичного испытания [1, 6], в которых можно найти ссылки на дополни-тельные источники. К сожалению, во всех этих источниках не обсуждаются или не рас-крываются критически важные детали технологии МПД-моделирования, без которых воспользоваться МПД не удается.
В докладе рассмотрены следующие вопросы реализации упомянутых технологий:
• логарифмическая шкала выдач; • раздельные датчики псевдослучайных чисел (ДСЧ);
• выбор оценок с наименьшей дисперсией;
• дополняющие ДСЧ;
• использование контрольных переменных;
• условная имитация.
Подчеркнем, что общей идейной основой МПД является использование дополнительных знаний о моделируемой системе. Предложения по МПД, не привлекающие таких знаний, могут «отметаться с порога» – по аналогии с проектами «вечных двигателей».


 

 

2.         Логарифмические выдачи

Исследования в области МПД, как правило, связаны с прогрессивным наращи-ванием числа испытаний – предпочтительно в ходе одного прогона с выдачей результа-тов при числе испытаний, соответствующем точкам логарифмической шкалы.
Типичное построение имитационной модели – цикл типа while, в теле которого стоит управляющее выбором вариантов (прибытие заявки, завершение обслуживания, поломка прибора и т. п.) разветвление. Непосредственно перед концом упомянутого цикла можно вставить фрагмент следующего вида:
if (n.gt.n1) then
k=mod(n,m)
if (k.eq.0) then
k=n/m
if (k.eq.1 .or. k.eq.2 .or. k.eq.5 .or. k.eq.10)then
print *, 'n = ',n
/*вывод нужных результатов>*/
n1=n
if (k.eq.10) m=10*m
end if
end if
end if

Здесь n – текущее число испытаний, а m – размер «текущей порции» (для начала, например, 1000). Работа с переменной n1 (ее стартовое значение 1) нужна для исклю-чения дублирования вывода, поскольку отработка некоторых событий (например, при-бытие заявки с постановкой ее в очередь) не связана с увеличением n.


 

 

3.         Раздельные датчики


Принципиальным условием применения МПД является независимость случай-ных величин каждого типа (для модели СМО – интервалов между заявками, длительно-стей обслуживания, типов заявок, узлов-преемников в сетях обслуживания и т. п.). Не-зависимость достигается использованием для их генерации раздельных ДСЧ. Для удоб-ства реализации датчики должны быть однотипны. Как правило, это генераторы муль-типликативного вида с общим множителем, но различными начальными значениями. Требуется гарантировать непересекаемость генерируемых серий, но конструктивные рекомендации на этот счет в литературе отсутствуют.
Упомянутую проблему можно решить, если разделить период генератора (для хороших ДСЧ это где m – разрядность датчика) на требуемое число серий и взять в качестве начальных значений числа с соответствующими номерами. Для получения этих чисел можно применить алгоритм ускоренного получения случайного числа, ос-нованный на двоичном разложении его номера k:
subroutine fastrand(k,p)
integer k,p
integer x0,z
integer j,j1,j2
data x0 /57539/
j=k
p=x0
z=1220703125
do while(j.gt.0)
j1=j/2
j2=j-2*j1
if (j2.eq.1) then
p=p*z
if (p.le.0) p=(p+2147483647)+1
end if
z=z*z
if (z.le.0) z=(z +2147483647)+1
j=j1
end do
end
Здесь как очередной «полуфабрикат» p, так и степень множителя z копятся по модулю датчика . Возведение в степень выполняется на каждом шаге, а домножение на p – в среднем на половине шагов, что дает для трудоемкости получения k-го числа оценку умножений. Заметим, что двоичный логарифм миллиарда не превосхо-дит 30. Алгоритм применим к любым начальным значениям и множителям и (после очевидной модификации) – к любому модулю датчика.
Идея этого алгоритма некогда была предложена А. П. Ершовым для реализации возведения в степень в Альфа-трансляторе; применительно к имитации она использо-валась Ю. Г. Полляком [4], но была изложена неотчетливо и не связывалась с разнесе-нием серий чисел. Правильность приведенной программы подтверждается совпадением результатов прямой и ускоренной генерации для заданного номера.

Сформированные целые числа для получения делятся на модуль датчика или умножаются на обратную ему величину. На большой выигрыш в быстро-действии от этой замены рассчитывать не стоит. Согласно [3], трудоемкость умноже-ния вещественных чисел по отношению к присваиванию целых оценивается в 5 еди-ниц, деление – в 9, а вычисление стандартных функций – в 150. Следовательно, для по-казательно распределенных чисел выигрыш составит менее 2%.
Из полученного набора начальных установок по степени согласия результатов моделирования с эталонными можно выбрать наилучшую комбинацию датчиков.


 

 

4.         Выбор оценок с наименьшей дисперсией


В данном разделе уместно привести основные формулы для системы M/M/1:
• среднее время ожидания
• дисперсия времени ожидания
• средняя длина очереди
• дисперсия длины очереди

С их помощью легко установить относительный выигрыш от оценки средней длины очереди через среднее время ожидания по формуле Литтла. Ре-зультаты расчета сведены в табл. 1. Таким образом, реальный вы-игрыш от подобного пересчета возможен только при умеренной загрузке.
Еще один часто обсуждаемый способ этого класса – определение среднего вре-мени пребывания заявки в системе через статистическое среднее времени ожидания и известное априорно среднее время обслуживания. Он может дать заметный эффект лишь в тех случаях, когда дисперсия чистого времени обслуживания соизмерима с дис-персией времени ожидания. В нашем случае дисперсия обслуживания равна 1. Сравни-вая ее с числами из четвертого столбца, убеждаемся в незначительности выигрыша.


 

 

5.         Дополняющие ДСЧ


Большие надежды связываются с усреднением результатов, полученных с чис-лами , и с их дополнениями до единицы. При их совместном использовании в од-ном прогоне можно в полтора раза уменьшить трудоемкость обращений к ДСЧ – в за-висимости от состояния «кувыркающейся единицы»:
с=1
...........
if (c.eq.1) then
u=random()
else
u=1.0-u
end if
dt =-log(u)/lam ! для EXP закона
c=-c
Попутно отметим полезность этого приема для программирования любых аль-тернирующих процессов.
Для пары прогонов можно использовать дополнение начальной установки до модуля датчика – числа (проверено!) будут тождественными. Коэффициент кор-реляции между парными числами оказывается очень близким к -1. Известно [7] обобщение метода дополняющих переменных, в котором проводится m параллельных опытов со случайными аргументами и

В процессе моделирования используются не сами , а , полученные с их помощью через обратные функции распределения требуемых случайных величин. В ряде источников утверждается, что при этом модуль отрицательной корреляции (и со-ответственно, дисперсия среднего) будут уменьшаться из-за нелинейности упомянуто-го преобразования. Моделирование дало –0,6445 для показательного закона, –0,9475 – для закона Релея и практически –1 – для треугольного закона, что и предсказывалось в [6, c. 150]. Заметим, что преобразование нелинейно для всех перечисленных законов, в том числе, треугольного. Следовательно, сохранению –1 мешает не нелинейность, а асимметрия плотности распределения.
Наконец, известен еще один способ введения отрицательной корреляции откликов: «перекоммутация» датчиков, используемых для генерации интервалов между заявками и длительностей обслуживания соответственно. В табл. 2 для модели СМО M/M/1 c точным значением среднего времени ожида-ния w = 9,0000 сопоставляются результаты прямого моделирования и «глобальных» ус-реднений «конечных» результатов прогона. В этой таблице:
• вариант 1 – результаты стандартного моделирования;
• варианты 2, 3 и 4 – полусуммы с результатами , дополняющей начальной установкой и перекоммутацией датчиков соответственно;
• варианты 5 и 6 – результаты «мультиметодов» [7] второго и четвертого порядка.

Обращает на себя внимание очень медленная сходимость результатов «сольного» моделирования, что объясняется большой дисперсией времени ожидания (99,0). Чтобы довести среднеквадратическое отклонение статистического среднего до 0,01, требуется примерно (100/0,01)2 = 108 испытаний. Эта оценка хорошо согласуется с результатами моделирования. Усреднение результатов (случаи с дополняющими и перекрестными датчиками, в особенности, варианты 2 и 4) дает заметно лучшие результаты. Эффект усреднения вполне окупает удвоение общего числа испытаний, в чем легко убедиться по следующей строке варианта 1 из той же таблицы. «Мультиметоды» дают весьма неустойчивую картину, причем увеличение m не гарантирует повышения точности.
Отметим, наконец, неожиданный и нигде не обсуждавшийся результат: локальное усреднение, т. е. чередование основных и «дополняющих» переменных датчика в ходе одного прогона, хотя и уменьшает дисперсию результатов, но дает очень большое смещение оценок.

Анализ показал, что при упомянутом чередовании среднее логарифмически пре-образованных чисел остается единичным, а их дисперсия . Следовательно, второй начальный момент сгенерированного распределения длительности обслуживания с единичным средним составит 1,5 вместо необходимого при экспоненциальном законе 2,0. Согласно формуле Полячека–Хинчина

, среднее время ожидания уменьшится на четверть (при коэффициенте загрузки 0,9 – с 9 до 6,75), что и подтвердил эксперимент по 200 тыс. реализаций с «кувыркающимся» датчиком длительности обслуживания – было получено 6,7067.


 

 

6.         Контрольные переменные


Метод контрольных переменных состоит в уточнении статистического среднего случайной величины через одновременно наблюдаемую другую с точно известным математическим ожиданием. Запишем соответствующие формулы из [1, 6] применительно к средним длительности ожидания и длины очереди:

Здесь «крышками» помечены статистические оценки средних, а q – точное матожидание длины очереди. Оптимальный множитель c* дается формулой и уменьшает дисперсию в [ ]-1 раз (в этой формуле вычитаемое – квад-рат коэффициента корреляции). Корреляционный момент получался накапливанием произведений (W-w)(Q-q) в момент выборки заявки из очереди. Случайное время ожи-дания W определялось как разность текущего значения таймера и вносимого в паспорт заявки момента постановки заявки в очередь, а длина очереди Q – как число заявок в очереди после выбираемой на обслуживание. Делать это можно, так как распределение числа заявок в очереди перед прибытием меченой заявки совпадает с распределением числа прибывших за время ожидания (закон сохранения стационарной очереди, [5]). Дисперсии W и Q определялись стандартным способом.
Результаты применения описанного подхода (табл. 3), который можно считать «самым умным», оказались заметно точнее остальных – в особенности, при умеренном числе испытаний. Уже при N=10 тыс. испытаний погрешность среднего времени ожидания составила около 2%. Выигрыш f в уменьшении дисперсии асимптотически оказывается четырехкратным.
Поправочный коэффициент c* с увеличением N очень быстро стабилизируется и может быть определен по умеренному числу испытаний. Кроме того, с ростом N уменьшается поправочная разность (напомним, что в нашем примере q=8,1). С учетом обоих этих обстоятельств метод контрольных переменных позволяет ограничиться небольшим N.

К сожалению, выбор контрольных переменных весьма ограничен, а к попыткам сконструировать их искусственно приходится относиться критически. В частности, ра-бота с суммарной трудоемкостью очереди к удовлетворительным результатам не при-вела. Вероятной причиной является то, что здесь фактически нет новой контрольной переменной – обсуждаемая сумма является частью статистики по основной перемен-ной. Видимо, так же обстоит дело и с предлагаемой в [1] средней из 100 предыдущих задержек.


 

 

7.         Условная имитация


Условная имитация, как и контрольные переменные, комбинирует имитацию и дополнительные знания о моделируемом процессе, но является гораздо более гибким инструментом понижения дисперсии.

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

Применяя к его слагаемым формулу для дисперсии произведения независимых случайных величин, имеем для дисперсии искомой оценки
Точное знание вероятностей слоев cводит к нулю взвешенную сумму по-слойных вторых моментов длительности ожидания и значительно уменьшает Эта формула подсказывает организацию математического эксперимента: определить (предпочтительно аналитически) вероятности и – отдельно для каждого слоя – получить имитацией первые и вторые моменты времени ожидания. Успех расслоения обусловливается двумя требованиями:
• стратифицирующая переменная должна быть коррелирована с откликом (результатом);
• должны быть известны вероятности попадания в слои.
Наиболее естественной областью применения расслоенных выборок к имитации СМО является учет переменных внешних условий (изменение интенсивностей потоков, длительностей обслуживания, числа действующих каналов и т. п.). В частности, представляется целесообразным моделировать n-канальную СМО с ненадежными каналами как обычную систему с числом каналов от 1 до n и затем взвесить результаты с учетом распределения числа исправных каналов, полученного расчетом или моделированием процесса отказов и воccтановлений. Ограничением этого подхода является требование существования стационарного режима в каждом слое.

7.2. Общий случай
Пусть искомая характеристика СМО является математическим ожиданием функции двух независимых случайных величин X и Y:

причем аналитический вид g(x,y) неизвестен. Зафиксируем значение одного из случайных факторов и получим с помощью имитационной модели условный отклик
Тогда конечный результат можно получить численным интегрированием согласно
Поскольку влияние фактора Y учитывается аналитически, дисперсия будет меньше, чем в случае полного моделирования по всему пространству . Проиллюстрируем эту идею на аналитических моделях систем с очередями. Рас-смотрим расчет среднего времени ожидания заявки в системе M/G/1 для равномерного на интервале [0,3; 0,9] распределения времени обслуживания через соответствующую систему с регулярным обслуживанием (при переходе к имитационным моделям в по-следнем случае дисперсия результатов будет меньшей). По формуле Полячека-Хинчина для упомянутого равномерного распределения =0,6, =0,39 и при =1 среднее вре-мя ожидания w=0,4875. Моделирование системы M/D/1 при длительности обслужива-ния 0,3; 0,6 и 0,9 и 5000 испытаний дало соответственно 0,00627; 0,415 и 2,919. Под-ставляя эти результаты в формулу численного интегрирования по Симпсону, получаем

Непосредственное моделирование равномерно распределенной длительности обслуживания дало 0.439, то есть заметно худший результат. Разумеется, утроение чис-ла испытаний даст выигрыш в точности и на исходной модели (было получено 0,471). Однако работоспособность предложенного подхода несомненна. Реализация этой идеи имеет существенное ограничение: все частные варианты обсчитываемых имитационных моделей должны иметь докритический коэффициент загрузки . Например, при равномерно распределенной на отрезке [0,3;1,5] длительности обслуживания и =1 коэффициент загрузки составит 0.9, однако прогоны модели M/D/1 для будут бессмысленны из-за отсутствия стационарного режима. Комбинация аналитики и имитации эффективна и в случае периодических процессов обслуживания, методы аналитического расчета которых отсутствуют. Чтобы встроить в датчик интервалов между заявками циклическое изменение интенсивности входящего потока, возможны (в зависимости от системы программирования) следующие варианты:
• объявить таймер модели в общем блоке;
• оформить датчик как внутреннюю процедуру или ее аналог;
• сделать таймер параметром датчика.
Ясно, что прямая работа с таким датчиком увеличит дисперсию времени ожидания. С другой стороны, можно прогнать модель для интенсивностей потока, соответствующих узлам квадратурной формулы (для составной формулы Симпсона – 5 узлов от минимума до максимума с шагом фазы ). Тогда

(стандартный множитель перед суммой равен h/3, сумма делится на длину интервала интегрирования по фазе). Опыт подтвердил хорошую точность даже при малом числе испытаний.
В рассмотренной сумме наибольшую роль играют два последних слагаемых – из-за гиперболического роста ожидания по коэффициенту загрузки и коэффициента «4» при четвертом слагаемом. Они же обладают наибольшей дисперсией. Поэтому ограниченный общий лимит испытаний следует распределять между точками пропорционально произведениям квадратов весовых коэффициентов на вторые моменты длительности ожидания, определенные по данным пробных прогонов.

7.3. Определение вероятностей редких событий Условная имитация является единственным средством получения достаточно точных оценок вероятностей редких событий (например, отказа высоконадежных систем). Для решения такой задачи искомую вероятность записывают как сумму произведений условных вероятностей отказа или других неблагоприятных событий (предполагается, что эти вероятности существенно больше искомой). Затем вероятности звеньев цепочек определяются на серии «условных» моделей, а при возможности – аналитическими методами. Наконец, они подставляются в исходную формулу. К сожалению, принцип оценки потенциального выигрыша по числу испытаний в немногих затрагивающих данный вопрос статьях [2] четко не формулируется.
Ключевой идеей оценки является обратная пропорциональность требуемого числа испытаний искомой вероятности и квадрату ее допустимой относительной погрешности. Пусть искомая вероятность
причем составляющие имеют порядок 0,01. Тогда сама  имеет порядок , и для ее прямого определения с относительной погрешностью 0,1 требуется испытаний, где k – функция доверительной вероятности.
При перемножении их относительные погрешности будут складываться, и для расчета каждой из них потребуется испытаний. Выигрыш по числу испытаний будет в раз (для сохранения доверительной вероятности придется пойти на некоторое его уменьшение).


 

 

8.         Выводы


1. Стандартная технология моделирования не дает трех верных знаков среднего времени ожидания сильно загруженной системы M/M/1 даже при числе испытаний, измеряемом миллионами. В таких случаях необходимо применение МПД.
2. Тестирование всех обсуждаемых в литературе МПД показало их работоспособность и в то же время – существенно различную эффективность.
3. Исследования в области технологий имитационного моделирования требуют получения и представления результатов в логарифмической шкале. Предлагаемая организация цикла моделирования эту проблему решает.
4. Приведенная программа ускоренного прогона мультипликативного ДСЧ позволяет легко получать стартовые значения однородных датчиков, разнесенных сколь угодно далеко по длине периода.
5. Ожидания от выбора оценок с наименьшей дисперсией оказались преувеличенными. В частности, выигрыш от расчета средней длины очереди через среднее время ожидания при загрузке свыше 0,7 оценивается приблизительно в 2-p.
6. Метод дополняющих переменных должен применяться только к прогонам в целом. Предпочтительна его версия с перекрестным использованием датчиков.
7. При реальном использовании новой информации эффективен метод контрольных переменных. Он требует умеренного числа испытаний.
8. Формула для дисперсии условного матожидания позволяет оценить потенциальный выигрыш от комбинированного использования аналитических вероятностей и условных имитаций. Предложенная комбинация имитации и численного интегрирования дает возможность быстрого и точного анализа периодических режимов, но требует существования стационарных режимов во всем диапазоне обсчитываемых условий.
9. Приведенная схема оценки эффективности расчета вероятностей редких событий методом условной имитации четко выделяет компоненты эффекта и указывает пути их вычисления.


 

 

9.         Литература


1. Кельтон В., Лоу А. Имитационное моделирование/Пер. с англ. –СПб.: Питер; Киев: Изд. группа БХВ, 2004. – 847 с.
2. Коваленко И. Н., Кривуца В. Г., Кузнецов Н. Ю. Опыт практического применения методов статистического моделирования в теории надежности//Кибернетика.–1987. № 5, – С. 111–117.
3. Мик Б., Хит П., Рашби Н. Практическое руководство по программированию/пер. с англ. –М.: Радио и связь, 1986. – 168 с.
4. Полляк Ю. Г. Вероятностное моделирование на электронных вычислительных машинах. М.: Сов. радио, 1971. – 400 с.
5. Рыжиков Ю. И. Теория очередей и управление запасами. –СПб.: Питер, 2001. – 376 с.
6. Рыжиков Ю. И. Имитационное моделирование. Теория и технологии. –СПб.: КОРОНА принт. – 2004. – 384 с.
7. Nelson B. L. Antithetic-Variable Splitting for Steady-State Simulations//Europian J. of Operational Research. – 1988. V. 36. – P. 360–370.

 


Распечатано с портала GPSS.RU (c) Ю. И. Рыжиков , 2005 г.