Оптимизация и параллельная реализация статистического моделирования диффузионных процессов

Для описания различных явлений в физике и технике широко используются математические модели, в основе которых лежит понятие стохастического диффузионного процесса. Подобные модели применяются, например, при расчете концентрации загрязнений в водной и воздушной средахпри исследовании оценок устойчивости и надежности радиотехнических, и телекоммуникационных систем (см., например, [26, 23, 19, 18]). В общем случае, стохастическая диффузионная задача может рассматриваться как вычисление функционала, представляющего собой математическое ожидание некоторой функции от случайной траектории диффузионного процесса. Как известно, математические модели диффузионных процессов определяются системами стохастических дифференциальных уравнений (СДУ) [26, 23] и поэтому диффузионная задача состоит в вычислении функционала от траекторий решения СДУ.

Для численной оценки подобных функционалов обычно применяется метод статистических испытаний или метод Монте-Карло (см., например, [2, 15]). Для этого решение СДУ аппроксимируется с помощью подходящей дискретной стохастической схемы, обеспечивающей сходимость численного решения к точному с некоторым порядком в среднеквадратическом смысле (см., например, [26, 19]). Соответственно, значение функционала Ф от точных диффузионных траектории аппроксимируется значением функционала (р от приближенных траекторий. В общем виде метод Монте-Карло заключается в оценке функционала ср выборочным средним от независимых реализаций (i, i = 0,1,., iV случайной оценки (, являющейся некоторой функцией от приближенных траекторий:

1 N м г= 1.

При этом предполагается, что у случайной величины? существуют и конечны ее математическое ожидание Е£, и дисперсия.

Практически важной проблемой при практическом использовании метода Монте-Карло является построение случайных оценок с малой трудоемкостью (см., например, [2, 15]). Вообще говоря, трудоемкость случайной оценки следует определять в зависимости от существа конкретной задачи. Разберем это положение подробнее. Пусть t (C) — среднее время ЭВМ, затрачиваемое на моделирование одного выборочного значения случайной оценки Как известно, случайная N величина N" 1 ^ Q при достаточно большом N имеет нормальное расi= 1 пределение и при заданном уровне доверия р выполняется следующее неравенство:

DC.

N' где константа 7 зависит от величины р. На основе этого неравенства, абсолютная статистическая погрешность статистической оценки определяется величиной.

DC, а относительная статистическая поN.

DC.

В случае малых значении функгрешность — величинои,. (ЕСг • N ционала ip = ЕС <С 1 точность расчетов целесообразно определять величиной относительной погрешности, и поэтому в этом случае под трудоемкостью *5© случайной оценки следует понимать среднее количество вычислительной работы, необходимой для достижения заданной относительной погрешности (см., например, [2]), а именно:

Диссертационная работа посвящена проблемам вычисления именно малых функционалов и поэтому в ней по существу будет использоваться это определение трудоемкости. Подчеркнем, что целью исследований является построение случайных оценок диффузионных функционалов с малой трудоемкостью, а не только с малой дисперсией.

В диссертационной работе разработаны методы улучшения статистического моделирования для широкого класса диффузионных задач. При этом детально исследовались методы оптимизации случайных оценок следующих практически важных функционалов: во-первых, вероятности не достижения границы заданой области диффузионным процесом в определенном интервале времени и, во-вторых, полной концентрации траекторий диффузионного процесса в заданной точке в определенном интервале времени [15]. В первом случае особенностью задачи является большая величина временного интервала, а во втором случае — большое расстояние от источника траекторий до заданной точки, что в обоих случаях приводит к весьма малой величине искомых функционалов. Прямое моделирование для статистических оценок для таких функционалов трудоемко, поэтому встает вопрос об улучшении вычислений.

Для решения вопроса уменьшения трудоемкости статистических оценок диффузионных функционалов существует целый рад подходов (см., например, [2, 15, 24, 21, 31]). В применении к решаемым проблемам наиболее подходящими представляются методы оптимизации случайных оценок многократных интегралов и случайных оценок решений интегральных уравнений второго рода [2, 15]. Эти соображения основаны на том, что функционал от приближенных траекторий диффузионных процессов можно либо выразить в виде многократного интеграла, либо представить в виде решения интегрального уравнения второго рода.

В диссертационной работе для уменьшения трудоемкости случайной оценки вероятности не достижения границы заданной области диффузионным процессом использовались, с одной стороны, метод расщепления и, с другой стороны, весовое моделирование на основе принципа выборки по важности [2, 15]. При этом проделан сравнительный анализ целесообразности применения этих методов оптимизации.

Для задачи расчета полной концентрации траекторий диффузионного процесса в заданной точке в определенном интервале времени использовалось представление концентрации в виде решения интегрального уравнения второго рода и далее для уменьшения трудоемкости оценки решения этого уравнения применялась комбинация метода расщепления и моделирования по ценности [2, 15].

Во многих случаях эффективной методикой ускорения расчетов на ЭВМ служит их распараллеливание [4, 15]. При этом методы Монте-Карло находятся, безусловно, в более «привилегированном» положении по отношению, к примеру, к конечно-разностным алгоритмам вычислительной математики. Это прямо следует из их структурыметоды Монте-Карло, по сути, заключаются в моделировании независимых реализаций случайных процессов. Поэтому алгоритмически задача эффективного распределения независимых реализаций по независимым процессорам может быть с успехом решена. Как известно, основное машинное время многопроцессорного вычислительного комплекса уходит на обмен данными между процессорами и на время простоя процессоров [4]. При распараллеливании методов Монте-Карло подобные временные затраты сведены к минимуму: здесь время, затраченное на распределение заданий по процессорам и на финальное осреднение независимых результатов практически не играет роли, особенно если выполняется моделирование большого числа независимых реализаций случайных оценок [15]. При этом, структура методов Монте-Карло позволяет добиться обратно пропорциональной зависимости величины трудоемкости случайных оценок от числа процессоров (при условии, что используемые процессоры имеют одинаковую производительность).

В связи с задачей распараллеливания статистического моделирования возникает, без преувеличений, основополагающая проблемаиспользование специально приспособленных для этого генераторов псевдослучайных чисел. Традиционно, в основе моделирования случайных процессов лежит использование равномерно распределенных на единичном интервале случайных величин, для численной реализации которых можно использовать конгруэнтный генератор псевдослучайных чисел [2, 15]. Естественно, в параллельной системе необходимо использовать обоснованные и проверенные статистическими тестами методы распределения псевдослучайных чисел по независимым процессорам.

Здесь, для полноты освещения проблемы, следует упомянуть об одном недостаточно обоснованном, с точки зрения автора диссертационной работы, подходе. В работах [27, 22] для распределения псевдослучайных чисел по процессорам на каждом из них предлагается использовать свой множитель конгруэнтного генератора, т. е. по сути, разные генераторы на каждом из процессоров. Однако, например, одна только проблема обоснования некоррелированности получаемых последовательностей псевдослучайных чисел и подбора соответствующих параметров генераторов представляется весьма сложной.

Для моделирования независимых реализаций оценки на каждом независимом процессоре можно использовать свою, непересекающуюся с другими, подпоследовательность исходной последовательности псевдослучайных чисел [2, 15]. При этом, для практического применения распараллеливания необходимо быть уверенным в том, что для конкретного генератора указанные подпоследовательности псевдослучайных величин в совокупности удовлетворяют необходимым статистическим тестам.

В диссертационной работе предлагается модифицированный генератор псевдослучайных чисел, равномерно распределенных в единичном интервале, с «астрономически» длинным периодом. В его основе лежит предварительно протестированный в работе [22] генератор. Предлагаемый в диссертационной работе модифицированный генератор является пригодным для широкого круга практических задач и, в особенности, для вычисления диффузионных функционалов, где на каждом процессоре велика потребность в большом количестве псевдослучайных чисел. Для построенного модифицированного генератора были успешно реализованы тесты многомерной равномерности и его предлагается использовать в составе описанной далее системы MONC. Отметим, с другой стороны, что предлагаемый генератор является 128-разрядным, т. е. в перспективе он может применятся на будущих мощных компьютерах с 128-ю двоичными разрядами машинного слова.

В диссертационной работе для распараллеливания методов статистического моделирования предлагается универсальная программная система — система MONC (сокращение от Monte Carlo), позволяющая распределять вычисления независимых реализаций статистических оценок (случайные испытания) по независимым вычислительным процессорам и, фактически, специально уменьшать величину t© среднего времени моделирования одного выборочного значения случайной оценки и, следовательно, уменьшать ее трудоемкость. С помощью системы MONC возможно моделировать на ЭВМ весьма большое число реализаций, необходимое из соображений точности оценок. Подчеркнем особо, что данная система может быть использована для широкого круга задач статистического моделирования, не ограничиваясь, конечно, только лишь диффузионными задачами.

Далее следует краткое содержание диссертационной работы по главам.

В первой главе предлагаются методы улучшения статистического моделирования для широкого класса диффузионных задач. При этом детально исследуются методы оптимизации случайных оценок следующих практически важных функционалов: во-первых, вероятности не достижения границы заданой области диффузионным процесом в определенном интервале времени и, во-вторых, полной концентрации траекторий диффузионного процесса в заданной точке в определенном интервале времени. Для уменьшения трудоемкости случайной оценки вероятности не достижения границы заданной области диффузионным процессом применяются, с одной стороны, метод расщепления и, с другой стороны, весовое моделирование на основе принципа выборки по важности. Для задачи расчета полной концентрации траекторий диффузионного процесса в заданной точке в определенном интервале времени используется представление концентрации в виде решения интегрального уравнения второго рода и далее для уменьшения трудоемкости оценки решения этого уравнения применяется комбинация метода расщепления и моделирования по ценности.

В параграфе 1.1 вводится система стохастических дифференциальных уравнений (СДУ), рассматривается класс функционалов от траекторий диффузионных процессов, рассматривается метод Эйлера для моделирования приближенных траекторий, делаются замечания о его общности при рассмотрении проблем оптимизации, вводится определение функционала от приближенных траекторий, представленного в рекуррентном виде, рассматриваются понятия дискретной и непрерывной постановок задач вычисления диффузионных функционалов.

В параграфе 1.2 рассматривается проблема вычисления методом Монте-Карло вероятности не достижения границы заданной области диффузионным процессом. Рассматривается приближение к точной диффузионной вероятности, которое будет использоваться в дальнейших исследованиях.

В параграфе 1.2.1 вводится аналоговая оценка вероятности не достижения границы заданной области диффузионным процессом.

Приводятся формулы для ее трудоемкости в рамках непрерывной и дискретной постановок задачи.

В параграфе 1.2.2 рассматривается применение метода многократного расщепления траекторий в точках временного интервала для улучшения аналоговой оценки. Приводятся известные соотношения о дисперсии расщепленной оценки, о количестве условно-независимых траекторий. На основе приближения к точной диффузионной вероятности делается вывод о выборе точек расщепления на временном интервале. Приводятся соотношения о выигрыше в трудоемкости от применения метода расщепления. Выводятся формулы трудоемкости расщепленной оценки в рамках непрерывной и дискретной постановок диффузионной задачи. Предлагается методика предварительной оценки параметров метода расщепления.

В параграфе 1.2.3 рассматривается весовая оценка вероятности не достижения границы заданной области диффузионным процессом. На основании принципа выборки по важности рассматриваются «идеальная» весовая оценка и «идеальная» плотность перехода. Доказывается теорема о величине дисперсии весовой оценки при использовании приближения к «идеальной» плотности перехода. На основании доказанной теоремы выводятся формулы трудоемкости весовой оценки в рамках непрерывной и дискретной постановок диффузионной задачи. Предлагается использовать приближение к точной диффузионной вероятности вместо неизвестной функции «ценности» в «идеальной» переходной плотности весовой оценки. Для улучшения весовой оценки в рамках дискретной постановки задачи предлагается использовать приближение к точной диффузионной вероятности, соответствующей увеличенной области (относительно заданной области). В качестве приближения к точной диффузионной вероятности предлагается использовать первую собственную функцию решения уравнения Колмогорова (в двумерном случае — это функция Бесселя индекса 0).

В параграфе 1.2.4 приведены численные результаты решения методических двумерных задач. Для весовой оценки на основе метода исключения описывается алгоритм моделирования траекторий, соответствующих новой переходной плотности. Расщепленная и весовая оценки исследуются в рамках непрерывной и дискретной постановок. Приведены результаты для весовой оценки, соответствующей увеличенной области. Делаются выводы о том, что при небольших по значению векторов сноса в рамках непрерывной постановки задачи весовая оценка особенно эффективнав рамках дискретной постановки задачи для СДУ с большими по значению векторами сноса предпочтительно применять расщепленную оценку, а для СДУ с небольшими по значению векторами сноса весовая оценка, соответствующая увеличенной области эффективнее расщепленной.

В параграфе 1.2.5 приведены исследования порядка детерминированной погрешности аналоговой оценки функционала методом зависимых испытаний.

В параграфе 1.3 рассматривается проблема расчета полной концентрации траекторий диффузионного процесса в заданной точке в определенном интервале времени. Рассматривается оценка по времени для концентрации.

В параграфе 1.3.1 рассматривается аналоговая оценка полной концентрации траекторий диффузионного процесса в выбранной точке в заданном интервале времени. Аналоговая оценка рассматривается в рекуррентном виде, на основе которого выводится интегральное уравнение второго рода для концентрации, которое рассматривается в расширенном фазовом пространстве, получаемом включением переменной t. Выводятся соотношения для дисперсии аналоговой оценки.

В параграфе 1.3.2 рассматривается весовая оценка полной концентрации траекторий диффузионного процесса в заданной точке в определенном интервале времени. Для улучшения оценки решения интегрального уравнения предлагается использовать метод моделирования по ценности. Сформулирована теорема о величине дисперсии весовой оценки при использовании приближения к «идеальной» переходной плотности. На основании этой теоремы выводятся формулы трудоемкости весовой оценки. Предлагается использовать приближение к точной диффузионной концентрации вместо неизвестной функции «ценности» в переходной плотности «идеальной» весовой оценки. В качестве приближения предлагается использовать точное решение для стационарной задачи для СДУ с постоянными коэффициентами сноса и диффузии.

В параграфе 1.3.3 на основе метода исключения описывается алгоритм моделирования весовой оценки. Делаются выводы об эффективности метода исключения. На основе этих выводов предлагается переходить на весовое моделирование только внутри ограниченного шара Ощ. Для повышении эффективности весового моделирования предлагается использовать метод расщепления. Переход к расщеплению предлагается делать при попадании траектории в шар содержащий шар Пру. Предлагается способ выбора радиуса шара 0,$ на основе приближенного представления решаемой проблемы, а именно, в качестве приближения рассмотривается задача о вычислении вероятности попадания траекторий диффузионного процесса в заданную область. Для определения количества условно-независимых траекторий в акте расщепления предлагается использовать специальные предварительные расчеты.

В параграфе 1.3.4 приведены численные результаты решения методических задач с известным решением. Для организации численных экспериментов применялась система MONC. Исследования проводились для двух различных шагов интегрирования для нескольких значений вектора сноса. Делается вывод о том, что при не слишком больших значениях вектора сноса применение комбинированной расщепленно-весовой" оценки существенно улучшает расчеты сравнительно с прямым моделированием.

В параграфе 1.3.5 приведены исследования порядка детерминированной погрешности аналоговой оценки концентрации методом зависимых испытаний.

Во второй главе рассмотрены различные аспекты параллельной реализации алгоритмов метода Монте-Карло на многопроцессорных вычислительных системах. Предлагается модификация длиннопери-одного «конгруэнтного» программного генератора псевдослучайных чисел, для которой успешно реализованы тесты многомерной равномерности. В приложении к главе приведены соответствующие вычислительные программы на языке Fortran 90.

В параграфе 2.1 рассмотрены различные аспекты параллельной реализации алгоритмов метода Монте-Карло на многопроцессорных вычислительных системах.

В параграфе 2.2 рассмотрены вопросы параллельной реализации генераторов псевдослучайных чисел. На основе конгруэнтного генератора (метода вычетов) вводятся понятия распределения псевдослучайных чисел по испытаниям (lf-генератор) и по процессорам (bf-генератор). Приводятся сведения о комбинированном использовании «физических» генераторов и генераторов псевдослучайных чисел.

В параграфе 2.3 описывается исследование параллельной реализации конкретного конгруэнтного 128-битного генератора. Для выбранной длины «прыжка» bf-генератора описывается методика статистического исследования модифицированного генератора на многомерную равномерность. Для применяемой методики статистической проверки по критерию х2 Для одномерного случая предлагается использовать оптимальное в смысле максимума мощности критерия X2 число классов. Показано, что модифицированный генератор удовлетворяет статистическим тестам на многомерную равномерность вплоть до восьмой размерности.

В параграфе 2.4 приводится сравнение численных результатов решения конкретных диффузионных задач при использовании предложенного генератора и известного конгруэнтного генератора с параметрами г = 40, М = 517. Показано, что для рассмотренных задач используемые генераторы статистически эквивалентны. Кроме того, проделанные в параграфе 1.3.4 расчеты показали, что модифицированный 128-битный конгруэнтный генератор позволяет получить удовлетворительное согласие статистической оценки для полной концентрации с известным точным ее значением. Это, по сути, является еще одной успешной проверкой предложенного модифицированного генератора.

В третьей главе предлагается универсальная программная система — система MONC (сокращение от Monte Carlo), позволяющая распределять вычисления независимых реализаций статистических оценок (случайные испытания) по независимым вычислительным процессорам. Описываются методологические принципы системы MONC, приводятся инструкции по ее использованию и методика организации численного эксперимента в сети персональных компьютеров.

В параграфе 3.1 описываются методологические принципы системы MONC. Предлагается специальная методика организации численного эксперимента в сети персональных компьютеров. Приводятся соотношения, позволяющие оценить производительность системы MONC.

В параграфе 3.2 описываются функциональные возможности системы MONC. Приводится схема взаимодействия компьютеров в составе системы, описываются функции выделенных компьютеров -" Серверов" и «Клиента» .

В параграфе 3.3 приводятся требования к Проекту пользователя: к программе, которую должен написать пользователь и которая предназначена для выполнения с использованием системы MONCк специальным файлам, необходимым для взаимодействия программы пользователя и системы.

В параграфе 3.4 описывается пример решения конкретной диффузионной задачи с использованием системы MONC. На этом примере показаны универсальные принципы организации численного эксперимента с использованием системы. Показано, что требования к Проекту пользователя не являются излишне ограничительными и поэтому система MONC является достаточно универсальной. Показано, что применение системы MONC при решении задач методом Монте-Карло предоставляет существенный выигрыш в быстродействии сравнительно с «непараллельными» вычислениями.

В приложении 1 приводятся необходимые для использования модифицированного 128-битного генератора вычислительные программы, константы и инструкции для использования.

В приложении 2 приводятся сведения по инсталляции и использованию системы MONC.

В параграфе 2.1 приложения 2 описывается служебное программное обеспечение системы MONC, приводятся описания необходимых файлов настройки системы и инструкции по инсталляции.

В параграфе 2.2 приложения 2 описывается интерфейс системы MONC, приводятся инструкции по ее использованию.

Основные результаты диссертационной работы опубликованы в работах [5]-[11]. Дистрибутивы программ и материалы, касающиеся предложенного генератора псевдослучайных чисел и системы MONC, опубликованы для широкого использования также на Web-странице автора http://osmf.sscc.ru/~mam.

На защиту выносятся следующие основные результаты.

1.) Разработана методика улучшения вычислений методом Монте.

Карло вероятности не достижения границы^{заданной} области диффузионным процессом за определенное время на основе метода расщепления.

2.) Разработана методика улучшения вычислений методом Монте-Карло вероятности не достижения границы области диффузионным процессом за определенное время на основе принципа выборки по важности.

3.) Разработана методика улучшения вычислений методом Монте-Карло полной концентрации траекторий диффузионного процесса в заданной точке в выбранном интервале времени на основе представления концентрации в виде решения интегрального уравнения второго рода и соответствующего применения комбинации метода расщепления и моделирования по ценности.

4.) Предложен и проверен статистическими тестами многомерной равномерности модифицированный 128-битный конгруэнтный генератор псевдослучайных чисел с «астрономически» длинным периодом, пригодный для распределения псевдослучайных чисел по испытаниям и по независимым процессорам.

5.) Разработана и внедрена система параллельных вычислений MONC, предназначенная для организации распределения независимых испытаний в методе Монте-Карло по независимым процессорам с финальным осреднением результатов расчетов.

Результаты, изложенные в диссертационной работе, докладывались и обсуждались на семинаре Отдела статистического моделирования в физике ИВМиМГ СО РАН, на конференциях молодых ученых ИВМиМГ СО РАН (1999;2001 гг.), на 4-м международном семинаре по математическому моделированию в г. С.-Петербурге (2001 г.), на международной конференции РаСТ-2001 в г. Новосибирске (2001 г.).

Пользуясь случаем, автор выражает искреннюю благодарность своему научному руководителю чл.-корр. РАН Г. А. Михайлову за руководство всей работой, за постоянное к ней внимание и за плодотворные обсуждения задач и результатов, а также д.ф.-м.н. С. М. Пригарину за полезные советы по написанию вычислительной программы генератора псевдослучайных чисел.

Заклю чение.

Сформулируем основные результаты диссертационной работы:

1.) Разработана методика улучшения вычислений методом Монте-Карло вероятности не достижения границы заданной области диффузионным процессом за определенное время на основе метода расщепления.

2.) Разработана методика улучшения вычислений методом Монте-Карло вероятности не достижения границы области диффузионным процессом за определенное время на основе принципа выборки по важности.

3.) Разработана методика улучшения вычислений методом Монте-Карло полной концентрации траекторий диффузионного процесса в заданной точке в выбранном интервале времени на основе представления концентрации в виде решения интегрального уравнения второго рода и соответствующего применения комбинации метода расщепления и моделирования по ценности.

4.) Предложен и проверен статистическими тестами многомерной равномерности модифицированный 128-битный конгруэнтный генератор псевдослучайных чисел с «астрономически» длинным периодом, пригодный для распределения псевдослучайных чисел по испытаниям и по независимым процессорам.

5.) Разработана и внедрена система параллельных вычислений MONC, предназначенная для организации распределения независимых испытаний в методе Монте-Карло по независимым процессорам с финальным осреднением результатов расчетов.