ПРОГРАММНЫЙ КОМПЛЕКС ВЕРОЯТНОСТНОГО АНАЛИЗА СИСТЕМ С КУСОЧНО-ЛИНЕЙНОЙ СТРУКТУРОЙ Текст научной статьи по специальности «Математика»

Электронный журнал «Труды МАИ». Выпуск № 72

www.mai.ru/science/trudy/

УДК 28.29.51

Программный комплекс вероятностного анализа систем с кусочно-линейной структурой Кан Ю.С., Травин А.А*.

Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993,

Россия

*e-mail:Dron-mail2001 @mail. т

Аннотация

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

Ключевые слова: вероятность, квантиль, линейные неравенства, выпуклые многогранники

1. Введение

Задачи вычисления вероятностного и квантильного критериев возникают при оценке терминальной точности систем управления летательными

аппаратами [1,2] и составляют предмет изучения специального раздела прикладной математики, находящегося на стыке теорий стохастического программирования и стохастического оптимального управления и известного под названием "вероятностный анализ". Вероятностный критерий по-существу моделирует вероятность достижения цели (целей) управления, а квантильный -гарантированный с заданной вероятностью результат принимаемого решения.

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

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

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

допускают естественное распараллеливание вычислений, что позволяет им успешно конкурировать с алгоритмами из [3].

Указанный подход применен в [4] для решения задач вероятностного анализа в случае квадратичной функции случайных параметров. Ниже он применяется для кусочно-линейной функции гауссовского вектора, размерность которого не выше трех. Предлагаемые ниже численные процедуры основаны на численном решении вспомогательной задачи вероятностного анализа - определении вероятностной меры выпуклого многогранника. Отметим, что такая задача возникает, как вспомогательная, на одном из этапов алгоритма решения задачи стохастического линейного программирования с квантильным критерием [5]. Указанные процедуры реализованы в виде программного комплекса, описание которого также приводится в настоящей статье.

2. Постановка задачи

Пусть -п-мерный случайный вектор, распределенный по нормальному закону N(Оп, 1п), где Оп -п-мерный вектор из нулей, 1п - единичная п х п матрица. Рассмотрим кусочно-линейную функцию потерь вида

Ф(й = тах [ат{ + Ь }, (1)

где а1 - детерминированный и-мерный вектор, Ь1 - детерминированная константа. Будем предполагать, что параметры функции (1) таковы, что она

достигает своего минимума в точке z0, и множество [г: Ф(¿) <р} является ограниченным выпуклым многогранником в Я" для любого р >Ф(г0) . Предположим также, что

же8п[г: Ф(г) = Ф(¿0)} = 0, (2)

где mesn - мера Лебега борелевских множеств в Я".

Если обозначить п = Ф(£), то вероятностный критерий, определенный согласно [3] выражением

Р(р) = Р(Ф(£ < р), (3)

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

Квантильный критерий для ае (0,1) определим выражением [3]:

ра = тт[р: Р(р) > а}. (4)

Величина (а является квантилью уровня а распределения случайной величины пи подлежит оценке. Отметим, что в силу сделанных допущений

Ра е (Р(¿0), +»).

3. Вспомогательные результаты

Из сделанных допущений легко следует [3], что Р(р) непрерывна и

строго возрастает по ре [Ф(г0), . Поэтому ра можно найти как единственный корень уравнения

Г (ф) = а (5)

относительно ф. Пусть Гк" (ф), Гк+ (ф) - некоторые последовательности функций распределения, поточечно сходящихся к Г(ф) в диапазоне фе [Ф(г0), , причем

Г"(ф) < Г(ф) < Г/ (ф), Уфе [Ф(го), +-). (6)

Теорема 1 [4]. Пусть для некоторых хп, уп е Я1 справедливы неравенства

Г- (хп) > а, Г+ (уп) <а. (7)

Тогда функция квантили фа, определенная выражением (4), удовлетворяет неравенству

уп <фа< хп. (8)

Более того, если функция Г(х) непрерывна на некотором отрезке [а, Ь] и фае (а,Ь) - единственный корень уравнения Г(ф)-а = 0, то для каждого п можно подобратьхп, уп, удовлетворяющие (8) и такие, что

хп - уп ^ 0 при п . (9)

Рассмотрим следующий алгоритм [4], порождающий последовательности хп и уп .

1. Выбираем точность £ и полагаем к = 1, п = 1, где к- счетчик, отвечающий за точность Ек - ¥к оценки функции распределения Р, а п- шаг алгоритма. На первом шаге полагаем а1 = а , Ь1 = Ь.

2а + Ь

2. Делим отрезок [ап, Ьп ] на три части и находим точки —"3—п = сп и

2Ь" + ап_ = ^п. Вычисляем значение функций ¥к, Рк в точках сп, (1п и переходим

3

к шагу 3.

3. Проверяем, где находятся вычисленные на шаге 2 значения и уменьшаем отрезок путем присвоения старых границ новым. Здесь возможны 6 вариантов расположения значений (см. рис. 1):

3.1. Если Гк (сп) <а< ¥- (^п) , то значению ап присваиваем значение сп , а значению Ьп значение с1п и увеличиваем п на единицу.

3.2. Если Рк(сп) >а , то значение ап не меняем, а значению Ьп присваиваем значение сп и увеличиваем п на единицу.

3.3. Если Гк (сп) < а и Гк (^п) < а < Гк (^п), то значение Ьп не меняем, а значению ап присваиваем значение сп и увеличиваем п на единицу.

3.4. Если ) >а и Рк(сп) <а< ¥к (сп) , то значение ап не меняем, а значению Ьп присваиваем значение dn и увеличиваем п на единицу.

3.5. Если (йп) <а, то значение Ьп не меняем, а значению ап присваиваем значение dп и увеличиваем п на единицу.

п п п ип II

3.5 3.6

Рис.1. Иллюстрация к алгоритму, порождающему последовательности хп и уп

3.6. Если Г- (сп) <а< (сп) и ) <а< Г+ (dп) , то увеличиваем к и

увеличиваем п на единицу.

Зацикливание алгоритма ввиду бесконечного повторения шага 3.6

произойти не может в силу предположения о том, что фа - единственный

7

корень уравнения Р (ф) -а = 0. Поэтому за конечное число шагов переходим в один из случаев 3.1-3.5.В результате получаем новый отрезок меньшей длины.

При этом в силу предложенной процедуры на каждом шаге алгоритма

2

неравенство (8) остается выполненным и хп+1 - уп+1 < хп - уп), откуда и следует (9).

4. Проверяем условие окончания алгоритма: I ап - Ьп 1< £ . Если условие выполнено, то переходим к шагу 5, иначе повторяем шаги 2-4.

а„ + Ь

5.Окончательно, для оценки фа, вычисляем ф

п_п_

а

2

4. Оценка вероятностных мер многогранников. Плоский случай.

Рассмотрим множество уровня А(ф) = (г: Ф(г) < ф} . В силу (1)

А(ф) = (г : атг + Ь1 < ф}, а значит Р(ф) = А(ф)). Таким образом задача вычисления Р (ф) сводится к нахождению вероятности попадания случайного вектора % в многоугольник, заданный системой линейных неравенств. Для начала определим вершины многоугольника. Составим из исходной системы

+ Ь =ф

атг + ь, =ф . ^удем искать все пересечения прямых,

подсистемы вида

входящих в исходные ограничения. Максимально имеем

2 т! т(т -1)

Ст =-=- подсистем при т > 3 . Для подсистем, одно из

(т - 2)!2! 2

неравенств которой является линейной комбинацией второго, решений не

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

порядке нахождения. Подставляя все найденные значения в исходную систему неравенств, найдем точки, удовлетворяющие всем ограничениям. Эти точки и будут вершинами, задающими многоугольник. Максимальное их число составит т.

Для удобства применения предлагаемых ниже алгоритмов нахождения вероятностных мер многоугольников полезно указать вершины ребер многоугольника. Достаточно будет найти все ребра многоугольника и последовательно составлять треугольники из ребра и заданного центра. Под ребрами в данной процедуре будем понимать наборы из координат двух вершин.

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

столбца матрицы АВТ вектор с, будем получать столбцы с двумя нулевыми элементами, номера которых и будут соответствовать парам вершин, образующих ребро.

Теперь рассмотрим алгоритм вычисления оценок (ф) и Р- (ф) для найденных фигур.

1. Если начало координат находится внутри многоугольника, то переходим к шагу 3, иначе к шагу 2.

Рис. 2. Иллюстрация к п. 2 алгоритма

2. Вычисляем вероятности попадания в

две фигуры G1 и G1 и А(ф),

см. рис. 2.

Для этого из центра проводим векторы через

все вершины многоугольника и центр

ОпА1,...,ОпАт . Вычисляем углы между всеми

парами

векторов

как

агеоо8

Х1Х2 + У1 У2

7 2 . 2 I 2 . 2

Х1 + У1\1 Х2 + У 2

и находим два

вектора, угол между которыми будет максимальным. Находим вершины, соответствующие этим векторам. Обозначим их АтЫ и Атах . Проведем

прямую

У - Уа.

х - х

А

У А - У А хА - хА

Атах Ат.п Атах Ат.п

-, разделяющую все ребра многоугольника

на два множества. В первое войдут все ребра, вершины которых расположены от центра до прямой, т.е. если выполнено

Уа, - УАт

х А - х А

Ак Ат

У А - УА. хА

А - х А

тах " *т1п Атах Ат1п

< 0, для каждой точки ребра. Во второе

множество - все остальные. Первая фигура будет образована ребрами из первого множества плюс ребра ОпАтах и ОпАтЛп , вторая - ребрами из второго множества, а так же ОпАтах и ОпАпйп . Применяя к данным

многоугольникам шаги 3-9 найдем оценки Г + и Г- для каждой фигуры. Вычитая из оценки вероятности попадания в большую фигуру оценку

вероятности попадания в меньшую, получим искомое значение. Переходим к шагу 3.

3. Проверяем многоугольник на выпуклость. Для этого нужно подставить координаты вершин в исходную систему. Если многоугольник выпуклый, т.е. при подстановке все неравенства выполняются, переходим к шагу 4.

4. Нумеруем ребра многоугольника в порядке их нахождения gm и переходим к шагу 5.

5. Разделяем многоугольник на треугольники, образованные центром Оп и ребрами ОпЛ1, ОпЛ>,ОпЛт_х,.., ОпАт, см. рис. 3. Вероятность попадания в каждый треугольник будем искать, деля их на более мелкие. Переходим к шагу 6.

Рис. 2. Способ разбиения многоугольника на треугольники

11

6. Лучами, выходящими из начала координат Оп , делим каждый

треугольник на к более мелких, к определяется исходя из заданной точности. Для этого находим величину угла ААхОп\ как

ZÁ1OnÁ2 = arceos

f 'QnAl2 + |QnA|2 -| Á Á2I2 л

21OnAf |OnÁ2|2

, и делим его на к углов,

. / Á O Á2 ^ 1

величины которых одинаковы и равны Ау =-1 n 2 . Найдем величину h

к

высоты треугольника Á1OnÁ2, проведенной из точки On. Переходим к шагу 7.

7. Рассмотрим каждый треугольник OnÁ1B1,OnB1B2,..,OnBkÁ2. Для упрощения доказательства сходимости алгоритма, а так же для вычисления априорной оценки точности вычисления функции F (ф) модифицируем разбиение треугольников на более мелкие. Если высота, проведенная из точки On находится внутри треугольника Á1OnÁ2, будем

разбивать треугольник BlOnBl+1 , в котором лежит высота на 2 треугольника. Соответственно необходимо переобозначить вершины B¡+. = Bh, Bl+2 = Bl+.Bk+. = Bk, где Bh - вершина, найденная при высоте.

Для треугольника OnÁ1Bl найдем величину углаAÁ1OnBl = у1 = 1Ау. Здесь

же найдем большую и меньшую стороны для каждого рассматриваемого треугольника

R = On\-Á-, ro = OnÁi, r = r ..

l n 1 sin(n-yl -ZOÁÁ) n l l-1

Если Я1 < г1 то необходимо переобозначить эти величины. Переходим к

шагу 8.

8. Оценка вероятности попадания в каждый треугольник есть вероятность попадания в круговой сектор соответствующего радиуса. Найдем оценки снизу и сверху для всех треугольников, используя Я1 и г1:

' - *?Л

1 - е

Ау

С г

- г1

2п

1 - е

Ау

2п

Переходим к шагу 9.

9. Суммируя все Рк+ и ^ , получаем оценки вероятности попадания в

каждый треугольник ОпА1 А2, ОпА2 Л3...ОпЛт-1 Ат:

р+ = V Р++, Р" = У Р

т / 1 к ' т / 1 к

Переходим к шагу 10.

10. Суммируя все (ф) и Рт (ф) по т, получим оценки вероятности

попадания в многоугольник сверху и снизу.

По построению и с учетом дополнительного разбиения в п.7 алгоритма можно заметить, что Я1 = г1+1 либо г = Я1+1, следовательно, можно сразу найти гарантированную границу погрешности для каждого большого треугольника

Л _Г> ¿2

Р +- р -

т т

-ОЛ -ОА

2

е 2 - е 2

Ау

2п

2

2

к

к

Суммируя погрешности для каждого треугольника, найдем гарантированную границу погрешности для всей фигуры

р р

^ -оА ~опА Л е 2 — е 2

Ау

/

2п

Из последнего очевидно, что при Ау^ 0 погрешность £ = р — р ^ 0 , следовательно алгоритм сходится для любой заданной точности. Так же это позволяет априорно вычислить к, т.е. число разбиений каждого треугольника, что значительно уменьшает время работы программы.

5. Оценка вероятностных мер многогранников. Трёхмерный случай.

Рассмотри трехмерный случай. Для решения вспомогательной задачи вычисления Р(ф) = е А(ф)) , где А(ф) = (г : ат г + Ъ1 < ф} , предварительно определяем геометрические параметры многогранника А(ф). С этой целью все грани А(ф) нумеруются и для каждой грани ищутся координаты всех прилегающих к ней вершин.

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

Способ разделения грани на треугольники.

Пусть имеется некоторая грань, у которой N вершин г1,..., .

- 1 N

1. Находим центр тяжести грани (выпуклого многоугольника) г = ^ ■

2. Составляем массив треугольников (массив наборов по 3 точки), беря последовательно каждые 2 соседние точки грани и центр тяжести.

Таким образом грань будет разделена на N треугольников. Далее будем делить эти треугольники на более мелкие, с учетом заданной точности.

Способ разделения треугольника на более мелкие.

1. Находим середины для каждой из сторон треугольника г*

яг _ + Х1+1

2

2. Соединив все новые точки и вершины исходного треугольника, получим 4 более мелких треугольника.

3. Повторим эти шаги до тех пор, пока сторона нового треугольника превосходит заданную.

Так, в зависимости от требуемой точности, разделим исходный треугольник на более мелкие и далее будем рассматривать каждый отдельно.

Составим из центра многогранника и полученных выше треугольников пирамиды. Нижняя оценка вероятности попадания в данную пирамиду будет равна вероятности попадания в сектор шара, образованный меньшим из ребер пирамиды.

Вероятность попадания в шар радиуса Я находим по формуле [6]

'2 -У

Р((X,У,г) с Вк) = 2Ф'(Я) -1 - Я.1-е 2,

к'

п

г 2

* 1 Г« —

где Ф (Я) = Г I е 2 йг - функция распределения нормального закона N (0,1). ^/2пJ-'x,

Вероятность попадания в часть шара, вырезаемую треугольной пирамидой со стороной равной радиусу шара Я

р((х ,у, г) с б) = ^ р((х ,у, г) с вк),

4п к

где & Б - телесный угол, вырезаемый треугольной пирамидой.

Треугольник с координатами вершин г1, г2, г3 виден из начала координат под

(г1, г2, г3)

телесным углом &Б = 2arctg------------— , где (г1, г2,г3) -

Г1Г2Г3 + (Г1 ■ Г2)г, + (г2 ■ Гз)Г1 + (Гз ■ Г1)Г2

смешанное произведение векторов, (г- ■ г})- скалярное произведение.

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

Заметим, что если точка пересечения высоты, проведенной из центра многогранника к каждому из треугольников, будет лежать внутри этого треугольника, то за меньший радиус нужно принимать величину найденной высоты. Проверить это можно, если найти все двугранные углы у основания пирамиды. Если все они острые, то высота, проведенная к основанию пирамиды будет лежать внутри.

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

6. Описание программного комплекса«ра8р18.ехе».

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

Для начала необходимо выбрать, в каком пространстве будут производиться расчеты - в 2-х или 3-х мерном. Выбор осуществляется переходом на соответствующую вкладку.

Система ограничений задается вводом линейных коэффициентов в систему неравенств. Неравенства заданы в форме (1). Программа позволяет добавлять нужное количество неравенств. Предельное количество ограничений - 50.

Далее задается точность оценки вероятности в пределах от 0,1 до 0,0000001.

Общая форма разделена на 2 рабочие области. В левой программа решает

вспомогательную задачу, т.е. вероятность попадания гауссовского вектора в

17

заданную область, а в правой вычисляется квантильный критерий для нее же. В первой рабочей области применены алгоритмы, описанные в пунктах 4 и 5, а во второй алгоритмы, описанные в пунктах 3-5. В каждой рабочей области есть кнопка «Вычислить».

После ввода всех данных, нажимаем кнопку «Вычислить».

В окне «Результаты вычислений» на каждой рабочей области отображаются верхняя и нижняя оценки вероятности и квантили, а также среднее величин.

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

7. Примеры расчётов.

Для плоского случая были выполнены расчеты вероятности попадания в многоугольники А и В с координатами вершин А1(-1;1), А2(1;2), А3(2;0), А4(1;-1), А5(0;-1) и В1(-1,1), В2(1,1), В3(1,-1), В4(-1,-1), см. рис. 4.

Таблица 1. Расчет вероятности попадания двумерного гауссовского вектора в многоугольник А.

Точность Количество Верхняя Нижняя Среднее

оценивания разбиений оценка оценка оценок

каждого вероятности вероятности

треугольника

0.01 48 0.4711 0 . 4611 0.4662

0.005 97 0.4686 0.4636 0.4661

0.001 478 0.4666 0.4656 0.4661

Таблица 2. Расчет вероятности попадания двумерного гауссовского вектора в

многоугольник В.

Точность Количество Верхняя Нижняя Среднее

оценивания разбиений оценка оценка оценок

каждого вероятности вероятности

треугольника

0.01 41 0.9132 0.9033 0.9082

0.005 83 0.9065 0.9015 0.90399

0.001 416 0.9012 0.9002 0.9007

Рис. 4. Иллюстрация к таблицам 1 и 2. Заметим, что для последнего многоугольника вероятность вычисляется

аналитически как 2(Ф0(а))2, где Ф0(а) = Ф*(а) -1/2 - функция Лапласа. Для а=1 это значение составляет 0.46605. Сравнивая это значение с найденным можно сделать вывод о точности работы программного комплекса.

8. Заключение.

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

Библиографический список

1.Малышев В.В., Кибзун А.И. Анализ и синтез высокоточного управления летательными аппаратами. М.: Машиностроение, 1987. — 302 с.

2.Кузьмин В.П., Ярошевский В.А. Оценка предельных отклонений фазовых координат динамической системы при случайных возмущениях. М.: Физматлит, 1995 — 296 с.

3.Кибзун А.И, Кан Ю.С. Задачи стохастического программирования с вероятностными критериями. М.: Физматлит, 2009. — 372 с.

4.Кан Ю.С., Травин А.А. О приближенном вычислении квантильного критерия // Автоматика и телемеханика. 2013. №6. С. 57-65.

5.Наумов А.В., Иванов С.В. Исследование задачи стохастического линейного программирования с квантильным критерием // Автоматика и телемеханика. 2011. №2. С. 142-158.

6.Вентцель Е.С. Теория вероятностей. М.: Высшая школа, 1999.— 576 с.