CC BY
58
УДК 681.3.06
В. Н. ЗАДОРОЖНЫЙ
Омский государственный технический университет
О КАЧЕСТВЕ ПРОГРАММНЫХ ГЕНЕРАТОРОВ СЛУЧАЙНЫХ ЧИСЕЛ
Рассматривается проблема качества датчиков случайных чисел, встроенных в языки программирования. Разрабатываются показатели, упрощающие анализ сходимости статистических оценок математического ожидания. Учитывается случай зависимых испытаний.
Ключевые слова: тестирование датчиков случайных чисел, статистическое моделирование.
1.Введение
Для реализации любых случайных факторов, имеющих заданные вероятностные свойства, в имитационном моделировании (ИМ) используются датчики базовой случайной величины (БСВ). По определению, БС В должна иметь равномерное распределение на интервале и = (0, 1), и, кроме того, каждое очередное значение на выходе датчика БСВ должно быть случайной величиной (сл.в.), статистически независимой от предшествующих значений. Идеальный датчик БСВ (ИД) - это датчик, точно удовлетворяющий требованию равномерности распределения на и и требованию независимости. Использование соответствующих теорем теории вероятностей позволяет путем преобразования идеальных БСВ получать идеальные реализации других сл.в. (или любых случайных объектов), обладающих требуемыми вероятностными характеристиками.
Существует множество статистических тестов, применяемых для проверки равномерности распределения и независимости чисел на выходе датчика БСВ. Хорошо известны и методы улучшения свойств недостаточно точных датчиков БСВ, такие, например, как метод Макларена-Марсальи [1], в котором числа на выходе одного датчика случайным образом перемешиваются с помощью другого. Опубликованная более тридцати лет назад библиография [2] содержит 491 ссылку на литературу по датчикам случайных чисел. В книге [3], вышедшей двумя годами позднее, Дж. Клейнен приходит, однако, к выводу, что «все еще нет высококачественного генератора псевдослучайных чисел». Одновременно Лыоис [4] и Марсалья [5] указывают: «Многие широко распространенные генераторы фактически непригодны». Актуальность проблемы качества датчиков случайных чисел не снижается и позже [6], и до настоящего времени [7,8].
Плохой датчик отбраковать достаточно легко, т.к. его недостатки проявляются в самых простых тестах. Но чем лучше датчик, тем сложнее выявляются его недостатки, если таковые (достойные внимания) имеются. Ермаков и Михайлов [9] отмечают, что в отношении датчиков «чрезвычайно важна проверка с помощью решения типовых задач, допускающих независимую оценку результатов аналитическими или численными методами. Можно сказать, что представление о надежности псевдослучайных чисел создается в процессе их использования с тщательной проверкой результатов всегда, когда это возможно».
Очевидно и то, что качество проверок, проводимых в отношении датчиков, требует не меньшего внимания, чем качество самих датчиков. Этим обусловливается критическое отношение автора настоящей статьи к большинству публикуемых проверок, послужившее отправной точкой к исчислению показателей, которые наглядно характеризуют сходимость статистических оценок при использовании ИД. Рассчитанные показатели сходимости приводятся в статье в качестве справочных данных, полезных для корректной интерпретации результатов «штатных» статистических экспериментов. Полезность предлагаемых показателей иллюстрируется примеромтого, как они позволяют глубже разобраться в качестве датчиков системы GPSS World и, одновременно, обнаружить проблемы в части корректности программной реализации этой системы.
2.Сходимость оценок
в независимых испытаниях
Рассмотрим пример из [8], в котором приводится наблюдаемая при использовании реального датчика «зависимость от числа испытаний N погрешности 8 расчета числа п методом статистических испытаний через долю точек, попавших во вписанный в квадрат круг» (табл. 1). Немонотонное убывание погрешности по числу испытаний приводится здесь как демонстрация того, что «реальные датчики равномерных псевдослучайных чисел отнюдь не идеальны». Но, заметим, при использовании ИД погрешность тоже убывает немонотонно!
При использовании ИД исходом каждого независимого испытания является попадание случайной точки в круг, с вероятностью точно р = 7t/4, или с вероятностью 1 — р, не попадание. В результате N испытаний вычисляется оценка р вероятности р, определяющая оценку 4р числа п. И, поскольку поведение оценки 4р однозначно определяется поведением оценки р, то сведем исследуемый вопрос о монотонности к анализу последней.
Обозначим через z, е {0,1} индикатор исхода/-го испытания: сл.в. z, принимает значение 1, если в i-м испытании произошло событие А, и значение 0 в противном случае. Оценка p = ^z,/N вероятностир имеет математическое ожидание M(p) = M(J]z(//v)=p и
дисперсию сг2 = /N) = р(1 - p)/N, а распределе-
ние вероятностей сл.в. р при больших N с высокой точностью описывается нормальным законом. Ошиб-
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК N*2 (ВО). 200?
ка х = р — р также распределена нормально с дисперсией ст2=р (1 — р)/1Ч, а ее математическое ожидание М(х) = 0.
Это - общеизвестные свойства оценки р, получаемой при использовании ИД. Исходя из них, можно определить вероятность О того, что при числе испытаний К, получится оценка более точная, чем при числе испытаний Ы,. Эта вероятность
О = Р(|х2|<|х,|), где х2, х, — независимые нормальные сл.в. с нулевыми средними и со стандартными отклонениями ст2 = ^/р(1 — р) / АГ2, ст, = л/р(1-р)/Л/| соответственно. Следовательно,
О = Р(|х2|<|х,|) =
= Р( — х2 < — х,, х2 < 0, х, < 0) +
+ Р(-х2 < -х,, х2 < 0, х, £ 0) +
+ Р(х2 < - х,, х2 £ 0, х, < 0) +
+ Р(х2 <х,, х2 £ 0, х, 2:0) =
= 4Р(х2 <х(, х2 £ 0, х, £ 0),
где все четыре вероятности равны в силу симметричности нормального распределения независимых х2 и хг Далее, переходя отх2 и х, к центрированным нормированным (стандартным) нормальным сл.в. х2 и х,, получаем:
0= 4Р(ст2х2<ст,х,, х2 £0,х, >0) =
= 4Р х2< — лг,, х2£0,х,£0
V а2
= 4Р(х2< ^N^/N1X1, х2£0,х,£0) = = 4)р0(к-'^0т,
(1)
где к = <т,/сг2 = ^Л/2/Л/, — коэффициент сжатия доверительного интервала, Г0 и Р0 — плотность вероятностей и, соответственно, дополнительная функция распределения стандартной нормальной сл.в. (Р0 =1-Р0).
Значения О, вычисленные по формуле (1) методами численного интегрирования, приводятся для некоторых М2/1М, в табл. 2.
Сопоставим данные табл. 1 и значения С? в табл. 2:
— в табл. 1 имеется 9 случаев удвоения числа опытов, и в 6 из них погрешность снизилась; согласно табл. 2, при удвоении числа испытаний погрешность оценки уменьшается с вероятностью О» 0,6;
— из 5 случаев увеличения числа опытов в 2,5 раза в табл. 1 находим 2 или 3 случая снижения погрешности (для N./14, =2,5 вероятность О ~ 0,64);
— из 12 случаев увеличения числа опытов вдесятеро имеем 9 случаев снижения погрешности (для ИД соответствующая вероятность О» 0,8);
— погрешность оценки в табл. 1 снизилась во всех 9 случаях увеличения числа опытов в 100 раз (для ИД при М2/М, = 100вероятность О»0,94).
Таким образом, становится очевидным, что немонотонное убывание погрешности, представленное в табл. 1, вполне возможно и при использовании ИД. При этом как сам факт немонотонного убывания погрешности в табл. 1, так и статистические параметры этого немонотонного убывания на самом деле свидетельствуют о хорошем качестве характеризуемого этой таблицей реального датчика. По данным табл. 1, этот реальный датчик неотличим от ИД.
На самом деле монотонная сходимость в ходе эксперимента гарантируется, как правило, не в отношении самой оценки, а в отношении каких-либо ее вероятностных характеристик, например, в отноше-
Таблица 1
Результаты испытаний датчика БСВ из [8]
N 5 N 5 N 5
1 тыс. 9.8е-2 50 тыс. 2.8е-3 2 млн -1.4е-4
2 тыс. 8.6е-2 100 тыс. 8.3е-3 5 млн -І.Ос-5
5 тыс. 5.6е-3 200 тыс. 4.3с-3 10 млн 2.9с-4
10 тыс. 18е-2 500 тыс. 4.3е-3 20 млн 4.2е-5
20 тыс. 2.2е-3 1 млн 2.0е-3 50 млн 1 2е-4
Таблица 2
Вероятность О уточнения оценки при независимых испытаниях
Ы2/Ы, к 0 Ы2/Ы, к 0 N2/^1, к О Ы2/Ы, к 0
1 1,000 0,5000 5 2,236 0,7323 40 6,325 0,9002 250 15,81 0,9598
2 1,414 0,6082 10 3,162 0,8050 50 7,071 0,9106 400 20,00 0,9682
2,5 1,581 0,6410 20 4,472 0,8600 100 10,00 0,9365 500 22,36 0,9715
4 2,000 0,7048 25 5,000 0,8743 200 14,14 0,9551 1000 31,62 0,9799
Таблица 3
Вероятность Ос уточнения оценки в продляемой серии независимых испытаний
Ы2/Ы, к <}с М2/Ы, к <2с N2/14, к <5с ы2/ы, к <3с
Ы2/Ы, 1 1 к і 1 <5*0,5 5 2,236 0,7500 40 6,325 0,9013 250 15,81 0,9599
2 1,414 0,6476 10 3,162 0,8128 50 7,071 0,9114 400 20,00 0,9682
2,5 1,581 0,6749 20 4,472 0,8631 100 10,00 0,9369 500 22,36 0,9715
4 2,000 0,7272 25 5,000 0,8766 200 14,14 0,9552 1000 31,62 0,9799
нии длины доверительного интервала, имеющего заданную надежность. Так, для рассмотренного эксперимента по расчету числа п «правило трех сигм» с вероятностью 0,9973 гарантирует, что при N = 50 млн оценка 4р числа я отличается от точного значения р не более чем на 3-(4^р[\-р)/» 0,0007. Сама же статистическая оценка сходится к точному значению стохастично, и всегда имеются шансы при небольшом N получить «чересчур точное» ее значение, которому проиграет «естественная точность» оценок, получаемых при дальнейшем увеличении N.
Согласимся, однако, что в реальных экспериментах, где ни вычисляемая величина, ни погрешность ее оценки не известны, и где поэтому классические методики типа «правила трех сигм» рекомендуют ориентироваться лишь на вероятностные характеристики оценки, мы все же нередко этими методиками не пользуемся. Использовать их не всегда удобно, зачастую они требуют приложения больших усилий и затрат времени, да и не всегда выполняются вероятностные условия, требуемые для их пригодности. И тогда, учитывая только поведение самой оценки, мы принимаем интуитивное решение, достаточно ли много испытаний выполнено для достижения текущих целей. Приведенные в табл. 2 показатели сходимости могут помочь лучше ориентироваться в подобных ситуациях. Нетрудно видеть, что вывод этих показателей элементарно распространяется на оценки математического ожидания любых сл.в. г.
Теперь заметим следующее: формула (1) и табл. 2 получены в предположении, что погрешность х2 формируется в Ы2 новых испытаниях, независимых от тех Ы, испытаний, в которых формировалась погрешность х,. В то же время при статистическом моделировании часто продляют одну и ту же серию испытаний и наблюдают за изменением оценки, вычисляемой по всей длине серии. В этом случае Ы2 исходов продленной серии испытаний включают в себя 14, первоначально полученных исходов, поэтому оценки х2 и х2 оказываются зависимыми, и характер сходимости оценок получается несколько иным.
3. Сходимость оценок в продляемой
серии независимых испытаний
Пусть по выборкег1,...,2м вычисляется оценка р, математического ожидания р сл.в. г, после чего выборка продляется и вычисляется оценка р2 =(г, +... + г„1 + г^... + г^)/М2. Найдем вероятность Ос = Р(|х2|<|х,|) того, что погрешность х2 = р2-р оценки р2 по абсолютной величине меньше погрешности х, =р,-р оценки р,. Обозначим дисперсию сл.в. г{ через а2.
Чтобы перейти от события |х2|<|х,|, определенного на зависимых нормально распределенных сл.в. х, и х2, к эквивалентному по вероятности событию, определенному на независимых стандартных сл.в., представим р2 в виде р2 =(г, + ... + гм )/М2 + (г„ ,г..+2ц )/Л/2 = р, ■(Л/,/Л/2) + р'-(Л/2 -Л^,)/Л/2,' где р' = {г^... + гыу /(Л/2-Л/,) — независимая оценка для раскрытая»). Погрешность можно тогда выразить через погрешность х' = р'-р «скрытой» оценки следующим образом:
N. . (А/,-Л/,)
’Р = ТГ"Р| + —-
х2=р2-
-р-р =
N... . (ЛГ2 - /V.).., .
Ы2 ' ЛГ2
(2)
Рис. 1. Зависимость О и Ос от N/N1, для двух стратегий увеличения числа опытов
независимыми нормально распределенными сл.в. с нулевыми средними и со среднеквадратичными отклонениями а/уЩ и стЛ/Л/2 -Л/, соответственно.
С учетом (2) искомая вероятность Ос может быть представлена в виде:
Ос = Р(|х2|<|х,|) =
= Р
/V,
<м
(3)
Раскрывая прямые скобки абсолютной величины с перебором четырех сочетаний знаков + и - (по аналогии с выводом формулы (1)), приводя затем в полученных неравенствах подобные члены, и переходя от х, и х' к стандартным нормальным сл.в. х, и х', сведем (3) к интегралу:
ос = г\ р0(7*7м-*)-р0(-
где, как легко видеть, погрешности X, и х' являются
где параметр к по-прежнему представляет собой ко-эффициент сжатия доверительного интервала:
ктщщ.
Рассчитанные по формуле (4) значения Ос приведены в табл. 3.
Сопоставляя данные таблиц 3 и 2, легко заметить, что вероятность уточнения оценки в стратегии продления испытаний немного выше, чем в стратегии их возобновления. Это объясняется появлением положительной корреляции между ошибками х, и х2, которая при случайном получении «чересчур малого» значения |х,| влечет снижение в среднем и значения |х2|. Разумеется, маргинальные распределения ошибки х2, а, следовательно, и безусловные средние значения |х2| для двух стратегий одинаковы. Но стратегия продления экономичнее, т.к. требует выполнить не 1Ч2, как стратегия возобновления, а только Ь12 — Ы, новых испытаний. Графики зависимости О и Ос от Ы2 /Ы, приведены на рис. 1.
Максимальная разность ДО = Ос —О » 0,044 достигается при Ы2 /Ы,» 1,43, а максимальное отношение 0с/0 я 1,08 — приЫ2/Ы,» 1,33.
Используя (2), нетрудно вычислить коэффициент г = у/ыуы.; =к~' корреляции между ошибками х, и Х2 в стратегии продления испытаний, и записать уравнение регрессии:
М(х21 х, = 0 = М (х2) + г^.(< - М(х,)) = ,
<7, Л/2
характеризуемой остаточной дисперсией
ОМСКИЙ НАУЧНЫЙ ВЕСТНИК (>0). 200» ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК N>2 (ВО).
0(х2|х, =0 = £>(х2Н1-Г) = —
ст2 (Л/2-/У,)
N. N.
которая, как легко видеть, при небольшом относительном удлинении выборки оказывается во много раз меньшей, чем безусловная дисперсия 0(х2) = = а71Ч2.
Разумеется, предположение о нормальном распределении ошибок х,, х' и х2 ориентирует применение полученных формул на случаи достаточно больших Ы, и (Ы2— Ы,), достигающих хотя бы нескольких десятков или сотен.
Доверительные интервалы
в зависимых испытаниях
Случай зависимых испытаний имеет место в следующей проверке (7):
«Для системы массового обслуживания М/М/1 теоретическая средняя длина очереди я = р2/ (1 — р) = = 0,92/(1 —0,9) = 8,1. Среднее время ожидания... должно быть равно ц/Х и в данном случае (1=1) численно совпадать со средней длиной очереди... [далее приводятся результаты моделирования, часть которых (табл. 4), слишком медленно сходится к точным значениям. — Прим. авт.]. Результаты, связанные с ожиданием, оставляют желать лучшего. Причина может быть только одна: недостаточно качественный генератор псевдослучайных чисел. Увеличение числа проведенных испытаний в общем приближает результаты к ожидаемым, но даже при 500 тыс. испытаний погрешность составляет около 17 %. Судя по динамике версий, вРБЭ/Ш проходит процесс интенсивной отладки. Можно надеяться, что встроенные в систему датчики равномерно распределенных чисел будут совершенствованы. В этом направлении можно приложить и собственные усилия».
Судя по всему, сомнений в состоятельности рассматриваемых оценок в работе [7] не возникает, и речь идет только о скорости убывания их дисперсии и смещения, слишком низкой по экспертной оценке автора (как будет видно далее, - правильной). Ноне будем спешить соглашаться с тем, что «причина может быть только одна», и «просчитаем» еще несколько не упомянутых причин замедленной сходимости оценок. Первая причина замедления состоит в том, что для расчета оценки используется выборка зависимых реализаций сл.в., между которыми имеется положительная корреляция. Действительно, если некоторая заявка попадает в длинную очередь и ее вре-
мя ожидания оказывается большим, то и следующая за ней заявка, скорее всего, попадет в длинную очередь, и, следовательно, ее время ожидания также окажется большим. Например, корреляция между реализациями wl,...1wN времени ожидания V/ изменяет дисперсию оценки йг математического ожидания М(ш) следующим образом:
■ _ 1 V
™ /V 1 ’
IV 1.1
г Я ЛЫ N
ч'->
. _1_ /V2
_1_
Ы2
М
/-І >-/♦! N-1 N
-ЛГ-АГ(нг)
ХМ(иг2)+2£ X
Ы /-1 /-/♦!
(5)
где а1 — дисперсия выборочного значения — не зависит от номера заявки і (мы рассматриваем стационарную последовательность реализаций .......,у/ы),
| N-1 N
7 =—X 2 гч — усредненный по выборке коэффи-” /-І >-/♦! циент корреляции между \У( и
Гу — СОГГ(УГі, ) =
^Щ]Щ)
М(иг(иг )-М2[мг)
Множитель (1 + 27) имеет смысл коэффициента увеличения дисперсии оценки.
Значения г() (рассчитанные путем ИМ, использующего качественные датчики, при длине выборки N = 500 000), в рассматриваемой системе М/М/1 медленно убывают с ростом расстояния б = |1—^ между элементами и выборки (рис. 2 слева). Корреляция г() существенно зависит от коэффициента загрузки системы р: при р= 0,75 она убывает на порядок быстрее (рис. 2 справа). Неплохая для целей приближенного расчета аппроксимация экспериментальных данных дается экспоненциальным представлением:
согг(н> 1,Ч> /+*)
М|М|1 загрузка 0.9
согг(н> і,И> і +,)
М|М|1 загрузка 0.75
202
Рис. 2. Корреляция времени ожидания заявок в системе М/М/1
Часть результатов моделирования системы М/М/1 в проверке [7]
Показатель Теория Число испытаний
50 000 200 000 500 000
Средняя длина очереди 8.100 5.132 5.757 6.690
Среднее время ожидания 8.100 5.148 5.785 6.703
Таблица 5
Результаты моделирования системы М/М/1 на GPSS World версии 4.3.2
Показатель Число испытаний
1 млн 10 млн 100 млн 1 млрд 2 млрд
Оценка средней длины очереди 8.332 8.061 8.087 8.043 7.799
Фактическая погрешность оценки 0.232 0.039 0.013 0.057 0.301
Допустимый уровень погрешности 0.404 0.128 0.040 0.013 0.009
Таблица 6
Результаты выполнения теста «расчет числа п» на GPSS World версии 4.3.2
Показатель Число испытаний
1 тыс. 10 тыс. 100 тыс. 1 млн 10 млн 100 млн 1 млрд 10 млрд
Абсолютная погрешность 0,0064 0,002 0,0054 4 10-* 9,9 10-* 8,4 Ю-5 3.5I0-* 2 10-6
Допустимая погрешность 0,16 0,049 0,016 5 10 3 1,6-Ю-3 510-1 1,6- Ю-* 5-Ю-’
Г,-'
fo.QSe-0'0055^, р = о,9 [о,9е"0'004Ч'"У|, /3 = 0,75
(6)
Соответствующие формулам (6) графики функций гч = т(|1—]|) = ф) изображены на рис. 2 сплошными линиями.
Используя в (5) экспоненциальную аппроксимацию, можно продолжить вычисление дисперсии О(Н') следующим образом:
1 N-l N . N-l N-I < ЛМ N-I
f—ZS-i - ІНІ.',,. - *11» -
1.1 y.M /V /e, /V /e, Ja|
= ;^Zr(sHN_s) * JJ ]a0e-a,(N-s)ds =
N
. N-1 , N-I
= — J Na0e"“ds- — I a0se~a,ds =
Л/
N
I
N-\
= a0 J e'a,ds-~j j se"a,ds ~
»a0 fe'“ds-lim— f se"asds = o0 fe'“ds = ^e" 0 J N-мо ft J 0 J л
l-r(l)
-n-nni
r(l)
l-r(l)'
(7)
Погрешность грубого приближения (7) можно охарактеризовать следующими сопоставлениями. Расчет величины 7 с помощью ИМ системы М/М/1 прир= 0,75 дает оценку 7 я 19, а по формуле (7) 7 я 16 (здесь г(1)» 0,943). При р= 0,9 найденные в ИМ значения г(э) охватили лишь часть значимого диапазона значений в (рис. 2) и остальная его часть в проверочном расчете была заменена аппроксимацией (6). При этом получено 7 я 70, а по грубому приближению (7) имеем 7 я 1 Ю (здесьг(1)я 0,991).
Подставляя (7) в (5), получаем простое приближе-
ние:
D(vir) = £L(l+27)«
сгЧі+грГ
N{\-r( 1),
(8)
где коэффициенты а0 и а соответствуют общему виду приближения для коэффициента корреляции г|( = ф) в экспоненциальной форме (6). Переход от коэффициента 7, усредненного по выборке конечной длины, к его пределу, соответствующему бесконечной выборке, выполнен в (7) в связи с естественным предположением, что ее длина N с большим запасом превышает значения 5, при которых корреляция ф) может заметно влиять на время ожидания заявок. Соответственно, часть интеграла, добавляемая при замене /Уна ос, пренебрежимо мала. Предпоследнее приближенное равенство в (7) соответствует переходу к грубой экспоненциальной аппроксимации функции ф), определяемой параметрами а0 = 1 и а = 1 — г( 1); последнее приближенное равенство в (7) учитывает, что коэффициент корреляции г(1) близок к единице.
При р= 0,9 найденные путем ИМ параметры а2 (дисперсия времени ожидания) и г( 1) равны примерно 82 и 0,991 соответственно.
Используем (8) для построения доверительного интервала оценки \у при значениях N. перечисленных в табл. 4. Подставляя эти N в (8) при известных г( 1) я 0,991 и ст2я82, находим, что при N = 50 тыс. Э(й') я 0,62, при N = 200 тыс. Э(й') я 0,32, и при N = 500 тыс. Б( й) я 0,22. Согласно правилу трех сигм доверительные отклонения оценки V/ ОТ ММ при этих значениях N составляют, соответственно, 3 0,6 = 1,8, 3 0,3 = 0,9 и 3 0,2 = 0,6. Если же О(й') рассчитывать по формуле (5) при 7 я 70, то получим для N = 50 тыс., 200 тыс. и 500 тыс. доверительные отклонения 1,4,0,7 и 0,5 соответственно. В табл. 4 отклонения оценок IV от 8,1 примерно вдвое превышают эти интервалы, каким бы из двух рассмотренных способов они ни рассчитывались. Таким образом, объяснить представленное в табл. 4 замедление сходимости оценок только влиянием корреляции нельзя, что подтверждает правильность цитированной выше экспертной оценки («результаты, связанные с ожиданием, оставляют желать лучшего»).
Правда, существует еще один фактор, замедляющий сходимость в зависимых испытаниях — это переходный процесс (ПП), который возникает в связи с фиксированным (например, пустым) начальным состоянием моделируемой системы, влияющим на ее ближайшее будущее. Однако, как показывают специальные имитационные эксперименты с системой М/М/1, при р= 0,9 влияние ПП на скорость сходимости здесь существенно ниже, чем влияние поло-
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ ОМСКИЙ НАУЧНЫЙ ВЕСТНИК N>2 (80). 2009
жительной корреляции между наблюдениями W| или L между q, (рис. 3). ш Слева на рис. 3 (верхняя кривая) показан график изменения средней длины q очереди, который получен усреднением ансамбля из 10 тысяч независимых реализаций процесса изменения очереди. Нижняя кривая - это график среднего значения оценки О средней длины очереди (оценка О вычисляется усреднением процесса по ходу модельного времени). График получен усреднением 10 тыс. независимых реализаций оценки О. На рис. 3 справа этот график показан в другом масштабе времени. Как видим, ПП протекает достаточно быстро. Сравнивая в этих экспериментах смещение оценки, вызванное ПП, и отклонение, вызванное дисперсией (5) (при 7 =70), находим, что уже при N = 500+600 смещение уверенно входит в пределы доверительного интервала. В целом оно убывает как Р1, т.е. значительно быстрее, чем среднеквадратичная составляющая погрешности, убывающая лишь как t-0,5. Поэтому при N = 5+10 тыс. смещение, вызванное ПП, становится не значимым.
Заметим, что при ИМ этой же системы на GPSS World студенческой версии 4.3.2 (copyright 2000) при тех же N, которые перечислены в табл. 4, получаются вполне приемлемые результаты, например, 8,702,8,362 и 8,727 соответственно. Хотя, как теперь легко подсчитать, результат 8,727 лежит на границе доверительного интервала. В связи с этим проведем для данной версии GPSS более полную проверку сходимости оценки w и рассмотрим ее результаты (табл. 5).
Как видим, при N = 1 млрд погрешность результатов моделирования оказалась существенно выше допустимой, а при N = 2 млрд — уже совершенно неприемлемой. И нам, пожалуй, пришлось бы поставить под сомнение качество работы датчика в дальних зонах N... если бы не одно «но»: слишком уж плохо он должен там себя повести, чтобы настолько резко испортить оценку.
Реабилитация датчика GPSS World
Проверочный эксперимент, отражаемый таблицей 5, имеет комплексный характер. В нем проверяется не только датчик случайных чисел, но и другие, скрытые от пользователя, инструменты системы моделирования, в том числе — алгоритм расчета среднего времени ожидания заявки. Есть основания полагать, что последние два столбца таблицы 5 сигнализируют не о проблеме датчиков, а о проблемах, связанных с корректностью использования машинной арифметики механизмами GPSS, скрытыми от пользователя.
Первое основание — это результат специального тестирования, свободного от использования скрытых
алгоритмов расчета показателей. Результаты этого тестирования приведены в табл. 6. Ее анализ показывает, что во всем диапазоне проверенных значений N абсолютная ошибка лежит в пределах доверительного интервала, т.е. датчик GPSS ведет себя как ИД. Это же показывают и перемены точности, статистика которых соответствует критериям, исчисленным в табл. 1.
Второе основание: в ходе работы над статьей в рассматриваемой версии GPSS действительно обнаружены дефекты машинной арифметики, проявляющиеся при большом числе испытаний N. После того, как число сгенерированных транзактов переваливает за 2 млрд номера транзактов - целые положительные числа — становятся отрицательными. Когда блок GENERATE при счетчике входов 2 147 483 646 пытается создать новый транзакт, возникает прерывание и выдается номер транзакта «XN: — 2 147 481 501» (отрицательный). Дополнительные проверки показывают, что отрицательные номера появляются гораздо раньше. Вот текст модели, позволяющий проверить вышесказанное:
GENERATE (Exponential (1,0,1))
QUEUE 1
SEIZE 1
DEPART
ADVANCE
RELEASE
TERMINATE
GENERATE
TERMINATE
1
(Exponential(l,0,0.9))
1
10000000000
1
Двоичный вид 1111111111111111111111111111111
«следующего» значения 2 147 483 647 счетчика, запрещенного к реализации и вызвавшего прерывание, обеспечивает контроль переполнения разрядной сетки, но, увы, запоздалый. Номера транзактов уже задолго до этого момента превышают свои пределы.
Выводы
1. При использовании ИД статистика колебаний точности подчиняется строгим вероятностным соотношениям. Вероятность повышения точности оценки при увеличении числа независимых испытаний с ^ до N2 рассчитывается численными методами и несколько различается в двух разных стратегиях наращивания числа испытаний. При этом стратегия продления серии независимых испытаний оказывается эффективнее, чем стратегия их возобновления.
2. В случае зависимых испытаний доверительные интервалы для оценки V/ математического ожидания М^) можно рассчитывать по методике, продемонстрированной на примере системы М/М/1 и позволившей найти простое приближенное выражение для
204
Рис. 3. Переходные процессы в системе М/М/1 при р = 0,9
коэффициента 1 + 2Г увеличения дисперсии оценки. Для вычисления этого коэффициента достаточно вычислить коэффициент корреляции г( 1) между двумя соседними элементами W, и wl+, выборки.
3. В системе М/М/1 при коэффициенте загрузки
0,9 имеем г(1)« 0,991, что приводит к увеличению дисперсии оценки w приблизительно в 150+200 раз и, соответственно, к такому же замедлению сходимости оценки.
4. Применение найденного коэффициента увеличения дисперсии подтвердило его адекватность как первого приближения, полезного для контроля точности зависимых испытаний и для анализа результатов проверки генераторов. Его применение позволило косвенным образом провести диагностику корректности программной реализации системы моделирования GPSS World версии 4.3.2 и обнаружить проблемы, связанные с использованием машинной арифметики.
Библиографический список
1. Харин Ю.С., Степанова М.Д. Практикум на ЭВМ по математической статистике : для мат. спец. ун-тов. — Минск, 1987. — 304 с.
2. Nance R.E. and Overstreet С. (1972). A bibliography on random number generation : Computing Rev., 13, P. 495-508.
3. Клейнен Дж. Статистические методы в имитационном моделировании : пер с англ. / под ред. Ю.П.Адле-
ра и В.Н. Варыгина. — М. : Статистика, 1978. — Вып.1 — 221 с.; вып. 2 — 335 с.
4. Lewis P.A.W. Large-Scale Computer-Aided Statistical Mathematics. Naval Postgraduate Shool, Monterey, Calif. / / Proc. Computer Science and Statistics: 6th Annual Symp. Interface, Western Periodical Co., Hollywood, Calif. 1972.
5. Marsaglia G. The structure of linear congruental sequences. — In: Applications of Number Theory to Numerical Analysis, S.K.Zaremba (ed), Academic, New York, 1972.
6. Соболь И.М. Метод Монте-Карло. — М., 1978. - 64 с.
7. Оценка системы моделирования GPSS World. Режим доступа: http://gpss.ru/paper/ryzhikov/index_w.html.
8. Рыжиков Ю.И., Соколов Б.В., Юсупов P.M. Проблемы теории и практики имитационного моделирования // Имитационное моделирование. Теория и практика : материалы 3-й Всероссийской конференции.— СПб: ФГУП ЦНИИТС, 2007. - Том 1, Том 2. Режим доступа: http://gpss.rU/immod07/doklad/6.html.
9. Ермаков С.М., Михайлов Г.А. Статистическое моделирование. — 2-е изд., доп. — М. : Наука, 1982. — 296 с.
ЗАДОРОЖНЫЙ Владимир Николаевич, кандидат технических наук, доцент, кафедра «Автоматизированные системы обработки информации и управления».
644050, г. Омск, пр. Мира, 11
Дата поступления статьи в редакцию: 26.05.2009 г.
© Задорожный В.Н.
yw681 3 06 В. Н. ЗАДОРОЖНЫЙ
Е. С. ЕРШОВ
Омский государственный технический университет
ГРАДИЕНТНЫЙ МЕТОД И ПРОГРАММА ОПТИМИЗАЦИИ СЕТЕЙ С ОЧЕРЕДЯМИ
Предлагается новый эффективный аналитико-имитационный метод оптимизации немарковсиих сетей массового обслуживания. Осуществляется программная реализация метода в среде АпуЦодх-б. Экспериментально оцениваются скорость сходимости и точность метода. Даются практические рекомендации по его применению.
Ключевые слова: сеть массового обслуживания, аналитико-имитационное моделирование.
1. Введение
Сеть массового обслуживания (СеМО) или, иначе, сеть с очередями (queueingNetwork), в которой циркулируют статистически однородные заявки, называется однородной сетью (1 ]. Узлами сети являются системы массового обслуживания (СМО) [2]. Заявки поступают в СМО, обслуживаются и передаются в соответствии с переходными вероятностями на входы других СМО или на выход сети (рис. 1—3). В терминах СеМО решаются многие задачи проектирования информационно-вычислительных систем (ИВС) [1-4]. При этом заявки интерпретируются как передаваемые
сообщения или как пользовательские запросы, обрабатываемые ресурсами ИВС. Теоретические методы анализа марковских сетей разработаны достаточно полно [ 1 — 3]. Успешно развиваются и методы синтезамар-ковских сетей, необходимые для оптимизации ИВС. Так, в [ 1) задача оптимизации производительности однородной замкнутой марковской сети сводится к системе л нелинейных алгебраических уравнений, эффективно решаемых численными методами.
В то же время для адекватной постановки задач оптимизации ряда ИВС и других систем приходится учитывать, что распределения интервалов поступле-
ОМСКИЙ НАУЧНЫЙ ВЕСТНИК N*2 (SO). 2009 ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
