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

На правах рукописи

Травин Андрей Александрович, ^

АЛГОРИТМЫ ОЦЕНКИ КВАНТИЛЬНОГО КРИТЕРИЯ С ЗАДАННОЙ ТОЧНОСТЬЮ В ЗАДАЧАХ СТОХАСТИЧЕСКОГО ПРОГРАММИРОВАНИЯ С КУСОЧНО-ЛИНЕЙНЫМИ И КВАДРАТИЧНЫМИ ФУНКЦИЯМИ ПОТЕРЬ

Специальность 05.13.01 Системный анализ, управление и обработка информации (авиационная и ракетно-космическая техника)

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

'и го

Москва, 2015

005562764

Работа выполнена на кафедре "Теория вероятностей" Московского авиационного института (национального исследовательского университета, МАИ).

Научный руководитель: доктор физико-математических наук, профессор

Кан Юрий Сергеевич

Официальные оппоненты: доктор физико-математических наук, профессор,

главный научный сотрудник НИИ механики МГУ им. М.В. Ломоносова Морозов Виктор Михайлович

доктор технических наук, старший научный сотрудник, руководитель программы аэрокосмических исследований ФГУП ЦАГИ имени профессора Н.Е. Жуковского Филатьев Александр Сергеевич

Ведущая организация: ФГАОУ ВПО «Национальный исследовательский

ядерный университет «МИФИ»

Защита состоится "23" октября 2015 г. в 10 ч. 00 мин. на заседании Диссертационного совета Д 212.125.04 Московского авиационного института по адресу: 125993, Москва, А-80, ГСП-3, Волоколамское ш., 4, аудитория 301 ГАК.

С диссертацией можно ознакомиться в библиотеке Московского авиационного института (национального исследовательского университета) и на сайте МАИ по ссылке: http://goo.gl/bX9BZ3.

Отзыв на автореферат, заверенный гербовой печатью организации, просьба направлять по указанному адресу в двух экземплярах. Автореферат разослан " о1 /" 2015 г.

Ученый секретарь диссертационного совета Д 212.125.04, кандидат физико-математических наук

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Актуальность темы.

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

Ярким примером современных сложных управляемых динамических систем являются ракетно-космические и летательные аппараты (ЛА), такие как пассажирские самолеты, орбитальные спутники и ракеты-носители. Эти Л А функционируют в условиях возмущений, связанных с отклонениями параметров атмосферы, аэродинамическими параметрами, разбросом тяги двигательных установок, ошибками навигационных систем, погрешностями систем управления. Случайные возмущения приводят к тому, что показатель свойств, количественно характеризующий достижение заданной цели функционирования ЛА, так же является случайным.

Для решения задач анализа точности, оптимизации, управления и оценивания состояния в динамических системах в настоящее время известен ряд подходов.

Задачи анализа движения орбитального корабля при спуске в атмосфере и при автоматической посадке рассмотрены в работах В.В. Малышева, А.И. Кибзуна, Ю.П. Семенова, В.А. Тимченко, Ю.Г. Сихарулидзе, A.A. Лебедева, задачи динамики и баллистики ракет в работах C.B. Беневольского, Н.Ф. Герасюта, Л.Н. Лысенко, Н.М. Иванова, задачи оптимизации, управления и оценивания состояния ЛА различных типов - Ю.П. Доброленским, В.М. Кейном, H.H. Красовским, A.C. Филатьевым, В.М. Морозовым. В этих работах для задач исследования движения ЛА обычно используются детерминированные, минимаксные или среднеквадратические критерии, а нелинейные модели движения во многих случаях линеаризуются.

В работах В.В. Малышева и А.И. Кибзуна предложен иной подход, основанный на использовании вероятностных критериев качества и наиболее адекватных, как правило, нелинейных моделей движения управляемых объектов. К вероятностным критериям относятся вероятностный и квантильиый функционал (в конечномерном случае -функции). Эти критерии позволяют оценить точность системы управления с учетом ограничений по надежности.

Математическая теория конечномерных оптимизационных задач с функциями вероятности и квантили в роли критериев оптимизации является одним из актуальных направлений стохастического программирования. Задачам стохастического программирования посвящены работы В. Купера, Р. Леппа, К.Марти, Б.Т. Поляка, А. Прекопы, Э. Райка, Р. Ветса, Ю.М. Ермольева, Ю.С. Кана, А.И. Кибзуна, A.B. Наумова, Д.Б. Юдина и многих других. Прикладная значимость этой теории обусловлена тем, что ука-

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

Практическая важность этих задач подтверждается также тем, что к их исследованию приложили руку известные российские специалисты в области математической теории управления: Афанасьев В.Н., Колмановский В.В., Бахшиян БД., Назиров P.P., Зубов В.И., Кац И.Я., Красовский H.H., Кузьмин В.П., Ярошевский В.А. Среди зарубежных исследователей можно выделить К. Борелла, А. Прекопу, Леппа Р., Тамм Э.,

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

Функция вероятности в задачах анализа систем управления техническими объектами, как правило, имеет вид кратного интеграла

где Ф(х) - исходная целевая функция, а /(х) - плотность вероятности вектора £ случайных параметров. В прикладных задачах параметр <р обычно моделирует допустимую точность системы управления. Таким образом, функцию вероятности можно иитерпри-тировать как характеристику "алгоритмической надежности"системы управления.

Важной характеристикой качества системы управления в условиях стохастической неопределенности является функция квантили (квантильный критерий)

где а € (О,1) - заданная доверительная вероятность. Ее можно интерпритировать как гарантированную с заданной вероятностью (надежностью) точность системы управления.

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

Юби Э.

ф(х)<<р

xa = min{(£ : Pv> а},

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

Для подсчета вероятности в инженерной практике используется метод статистических испытаний (Монте-Карло), обладающий простотой, универсальностью и позволяющий исследовать широкий класс моделей. Однако для проведения расчетов с использованием этого метода требуются большие вычислительные затраты, поскольку метод Моите-Карло основан па проведении многократных статистических испытаний исследуемой модели и оценке вероятности как частоты успешных испытаний. Поэтому, к сожалению, существует большой разрыв между достижениями теории и практическим применением этих достижений, во всяком случае - к очень важной области, относящейся к исследованию задач динамики полета. Численные методы для решения задач анализа точности, возникающих при исследовании движения авиационных и космических ЛА, в которых в качестве критерия непосредственно использовались вероятность и квантиль рассматривались Малышевым В.В., К.А Карпом, Ермаковым С.М., Михайловым Г.А., Вентцель Е.С., Володиным В.Д., Кибзуиом А.И.

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

Tifc+l = fk{Uk),

где к - номер шага (итерации). При этом на каждом шаге приходится для некоторого точностного функционала Ф/с(и, £) подбирать параметр ip^ таким образом, чтобы

f(<M"fc,f) <¥>*) >а.

Это необходимо для того, чтобы множество {i Е й° : ^к(ик,0 < фк} было доверительным. В работе A.B. Наумова и C.B. Иванова эта задача обозначена, но не решена ввиду ее сложности. На решение этой проблемы и направлена настоящая диссертация. Рассматриваемая задача является по сути задачей анализа систем управления и в дальнейшем зависимость от и (или щ) в тексте опускается.

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

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

Для достижения цели предполагается:

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

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

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

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

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

Достоверность результатов. Достоверность результатов обеспечивается:

1. Строгостью постановок и доказательств утверждений.

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

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

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

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

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

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

в области техники — рассмотрена задача вероятностного анализа рассеивания точек падения фрагментов летательного аппарата для оценки района поиска фрагментов.

Апробация работы.

Результаты работы докладывались па следующих научных конференциях: Международная конференция «Системный анализ, управление и навигация» (Евпатория, Украина, 2010, 2011 г.г.), 40-я международная молодежная научная конференция «Га-гаринские чтения» (Москва, 2014 г.), Международная конференция «Системный анализ, управление и навигация» (Анапа, 2014 г.), научный семинар кафедры "Теория вероятпостей"МАИ под руководством профессора А.И. Кибзупа, научный семинар в Центральном аэрогидродииамическом институте имени профессора Н. Е. Жуковского, научный семинар кафедры "Прикладная математика"НИЯУ МИФИ.

Публикации.

Основные результаты диссертации опубликованы в 3 статьях в журналах, входящих в перечень ВАК [1-3], а также в трудах и тезисах международных научных конференций [4-7]. Получено 1 свидетельство о государственной регистрации программ для ЭВМ [8].

Структура и объем диссертации.

Диссертация состоит из введения, трех глав, заключения и списка литературы (103 источника). Объем диссертации включает 99 машинописных страниц, включая 23 рисунка, 26 таблиц и 2 приложения.

СОДЕРЖАНИЕ РАБОТЫ '

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

В первой главе приводится обзор результатов в области вероятностного анализа. Приводится описание нового метода построения двусторонних оценок для кваитилыюго критерия. Суть его в следующем.

Пусть случайный вектор размерности от, Ф(е) : Rm Я1 - измеримая по Борелю функция. Тогда т/ = Ф(£) является случайной величиной с функцией распределения

(1) F{x) = P(V <х)= Р(Ф(Й < х).

Квантильный критерий (функция квантили) определяется выражением

(2) х„ = min{:c : F{x) > а},

где а е (0,1) — доверительная вероятность.

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

Рассмотрим пп и 17+ — последовательности случайных величин с (функциями распределения F~(x) и F+(x) соответственно, сходящиеся по распределению к 77 , причем < H*) < -ВД Vx.

Требуется оценить функцию квантили ха с заданной точностью е, используя последовательности F£(x) , F~(x). Предполагается, что значения функций F+(x) и F~(х) вычислить намного проще, чем F(x).

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

Теорем а 8. Пусть для некоторых хк,ук 6 Л1 справедливы неравенства

(3) F~{xk) > tt, F+Ы < а.

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

(4) Ук < ха < хк.

Более того, если функция F(x) непрерывна на некотором отрезке [а, Ь] и ха £ (а, 6) - единственный корень уравнения F{x) — а, то для каждого к можно подобрать хк,ук, удовлетворяющие (4) и такие, что

(5) хк -ук~) 0 при к ->■ 00.

Под индексом к в нижеописанных процедурах будем поиимать шаг алгоритма, т.е. число шагов, за которые достигается конкретная точность оценивания квантили £ = 1хк~2/kl- Индекс п далее характеризует точность оценивания функций F+(x), F~{x). Например для описанной ниже процедуры нахождения вероятности попадания двумерного гауссовского вектора в эллипс, п есть число разбиений эллипса на секторы. Чем больше п, тем лучше точность оценки. В этом конкретном примере найдена априорная зависимость точности оценивания - iv(x)| от числа разбиений.

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

1. Выбираем точность s и полагаем п = 1,к = 1, где п — счетчик, отвечающий за точность F+ - F~ оценки функции распределения F, а к - шаг алгоритма. На первом шаге полагаем аг — a,bi = b, т.е. определяем границы отрезка, в котором предположительно находится искомая величина.

о гт г l 1 2ак + Ьк 2 Ьк + ак

2. Делим отрезок [ак,Ьк\ на три части и находим —^—- = ск и —^—- = dk.

Вычисляем значение функций F+,F~ в точках ck,dk и переходим к шагу 3.

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

3.1 Если F*(ck) < а < F~(dk), то значению ак присваиваем значение ск, а значению Ък значение dk и увеличиваем к на единицу. Переходим к шагу 4.

3.2 Вели Р~(ск) > а, то значение ак не меняем, а значению Ьк присваиваем значение ск и увеличиваем к на единицу. Переходим к шагу 4.

3.3 Если < а и < " < -Р+(4), то значение Ьк не меняем, а значению ак присваиваем значение ск и увеличиваем к на единицу. Переходим к шагу 4.

3.4 Если 4) > а и Р~(ск) < а < Р£{ск), то значение ак не меняем, а значению Ьк присваиваем значение 4 и увеличиваем к на единицу. Переходим к шагу 4.

3.5 Если 4) < а, то значение Ьк не меняем, а значению ак присваиваем значение 4 и увеличиваем к на единицу. Переходим к шагу 4.

3.6 Если Р~[ск) < а < Р+(ск) и F~i.dk) < а < 4), то увеличиваем п на единицу и переходим к шагу 3.

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

При этом в силу предложенной процедуры на каждом шаге алгоритма неравенство (4) остается выполненным, и хм - ук+1 <Цхк- ук), откуда и следует (5).

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

Чк + Ьк

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

Из описания алгоритма видно, что ключевым моментом в успешной реализации метода является конструирование последовательностей и Далее по тексту диссертации эти последовательности строятся для кусочно-линейной и квадратичной функции потерь, зависящей от двумерного гауссовского вектора, а также для кусочно-линейной функции потерь, зависящей от трехмерного гауссовского вектора. Во всех перечисленных случаях приводятся алгоритмы вычисления значений функций ^(х) и Р~(х), которые используются в изложенном выше алгоритме. Случай квадратичной функции потерь, исследуется в первой главе, кусочно-линейные функции рассматриваются во второй.

Приводится описание процедуры вычисления квантилей нормы двумерного гауссовского вектора, в основе которой вышеописанный алгоритм.

Двумерный случайный вектор £ = (£ь&)Т распределён по нормальному закону с нулевым математическим ожиданием и невырожденной ковариационной матрицей К. Для заданного а е (0,1) определяется квантиль хи уровня а распределения нормы ¡¡£[|. В частном случае а = 1/2 такая задача возникает при оценке кругового вероятностного отклонения точки падения космического аппарата и рассмотрена в главе 3.

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

дится к вычислению интеграла

(6) F(R) = J/ е"2^ 72'dxdy.

X2W<TP

Заметим, что данный интеграл при 7 1 не вычисляется аналитически. При вычислении данной величины на компьютере, используя стандартные методы, для некоторых исходных данных можно получить F(R) > 1, что исключает возможность оценки квантильного критерия. Например, используя математический пакет Maple, при R = 5 и 7 = 0,0001 получаем F(R) = 1,000000357.

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

Z

Рис. 1: Оценка вероятности попадания в эллипс

В качестве аргументов для оценок выступают вычисленные с помощью п.2

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

1. Часть эллипса, расположенная в рассматриваемом квадранте, делится на п секторов с одинаковыми углами = как показано на рис. 1. Переходим к шагу 2.

2. Для оценки вероятности попадания в сектор эллипса рассматриваются два круговых сектора со сторонами Л„ и г„ = (см. рис. 1), между которыми расположен эллиптический сектор, гг = Я0 = Л. Вычисляем угол <рп, а так же большую и меньшую стороны круговых секторов.

л „ Я

д/72 cos2(Vn) + sin^n)

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

3. Так как вероятность попадания в круговой сектор с углом Д<р вычисляется как

т дчч-(.-.-*)£

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

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

4. Суммируем оценки,полученные на шаге 3.3. по всем эллиптическим секторам, тем самым находим оценки для а.к,Ьк,ск,(1к.

(8) Г=4^ 1-е"^ = 1-

где Щ и Г{ - больший и меньший радиусы для каждого рассматриваемого эллиптического сектора, п — число разбиений эллипса на секторы.

Таким образом Р~ < < При увеличении п можно добиться нужной

р+ р-

точности оценивания и оценить вероятность попадания в эллипс .Р(Я.) =---.

С учетом специфического деления эллипса на секторы, становится возможным вычисление гарантирующей границы погрешности оценок. Л е м м а 5. Справедливо соотношение

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

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

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

(9) Ф(0 = шах{аГе + Ь;},

1 = 1,ТП

где 0+ — детерминированный п-мерный вектор, Ь; — детерминированная константа. Предполагается, что параметры функции (9) таковы, что она достигает своего минимума в некоторой точке га , и множество {г : Ф(г) ^ ¡р} является ограниченным выпуклым многогранником в Я" для любого > Ф(^о). Предполагается также, что

тезп{2 : Ф(2) = Ф(г0)} = О,

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

Если обозначить 77 = <&(£)■ то вероятностный критерии, выражением

Р(р) = РШ) ^ Ч>\

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

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

<ра = тш{^> : Р[у>) > а}.

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

Далее приводится алгоритм оценки вероятностной меры многоугольника. Рассматривается множество уровня А(<р) = {г : Ф(г) < у). В силу (9) А[(р) = {г : ¿[г + к < уз}, а значит Р{<р) — Р(£ € А{<р)). Таким образом задача вычисления Р(<р) сводится к нахождению вероятности попадания случайного вектора £ в многоугольник, заданный системой линейных неравенств. Предварительно определяются геометрические параметры многоугольника А (<р), находятся все вершины и ребра (наборы из двух соседних вершин).

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

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

2. Вычисляем вероятности попадания в две фигуры б) и С] и А(<р) , см. рис. 2. Применяя к данным многоугольникам шаги 3-8 найдем оценки Р+ и Я" для каждой фигуры. Вычитая из оценки вероятности попадания в большую фигуру оценку вероятности попадания в меньшую, получим искомое значение. Переходим к шагу 3.

3. Нумеруем ребра многоугольника в порядке их нахождения дй-.,дт и переходим к шагу 4.

4. Разделяем многоугольник на треугольники, образованные центром 0„ и ребрами ОпАиОпА2,ОпАп-и-,ОпАт, см. рис. 3. Вероятность попадания в каждый треугольник будем искать, деля их на более мелкие. Если высота, проведенная из точки 0„,

находится внутри треугольника А10„А2 , будем разбивать этот треугольник па 2 найденной высотой Н. Соответственно необходимо переобозначить вершины и увеличить их количество на единицу Ак+2 = Ак+\ ,..., Ат+1 = Ат , Ак+1 = А,, ,где Ак — вершина, найденная при высоте. Шаги 5-8 необходимо применить к каждому треугольнику. Переходим к шагу 5.

Рис. 3: Оценка вероятности попадания в многоугольник

5. Лучами, выходящими из начала координат 0„ , делим каждый треугольник на п более мелких. Для этого находим величину угла АА\ОпА2 как

= агссоэ I^.Г + ^Г-И.У

и делим его на п углов, величины которых одинаковы и равны Д7 = ¿А\°чАгш Переходим к шагу 6.

6. Рассмотрим каждый треугольник ОпАхВх, ОпВхВ2,.., ОпВкА2. Для треугольника 0„А\В1 найдем величину угла /.А\ОпВ1 = 71 = (Д7- Здесь же найдем большую и меньшую стороны для каждого рассматриваемого треугольника

(10) Я, = 0„Л, ^ОпАгА2 = =

52П(7Г — 7; — АО„АхА2)

Если < Г| то необходимо переобозначить эти величины. Переходим к шагу 7.

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

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

8. Суммируя все и Я", получаем оценки вероятности попадания в каждый треугольник 0„А2А3...0„Ат_:Ат:

(12) ^Е*?.*1™ = £*■»■

к к

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

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

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

Л е м м а 6. Справедливо соотношение

( -|А7

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

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

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

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

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

N

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

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

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

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

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

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

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

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

Рис. 4: Способ разделения граней

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

Р((Х,У, Л) С Вк) = 2Ф*(Я) - 1 - Д^е-т,

где Ф*(Я) = /^оо®-^"^ — функция распределения нормального закона N(0,1).

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

Р((х,у,г) сс) = ^рцх.у.г) с до.

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

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

г. о 4 (п,г2,г3)

Пд = 2 агсЬд-

ПГ2г3 + (п • г2)г3 + (г2 • Г3)Г1 + (гз • Г1 )г2' где (г],г2,гз) — смешанное произведение векторов, • т^-)— скалярное произведение.

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

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

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

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

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

При обосновании требований к параметрам объектов наземной инфраструктуры информационных систем слежения за полетом ЛА одним из важнейших параметров, используемых для учета влияния рассеивания на результаты пуска является КВО, характеризующее степень рассеивания точек падения на поверхность Земли. КВО является мерой кучности пусков при круговом рассеивании ЛА и представляет собой радиус круга, вероятность попадания в который равна 0,5 при условии совмещения центра нормального закона распределения ошибок пуска с центром круга.

В соответствии с терминологией, представленной во второй главе, КВО представляет собой квантильиый критерий качества надежности 0,5 для функции потерь, равной величине случайного отклонения точки падения фрагмента от центра нормального закона распределения. При этом КВО обычно бывает известным для типовых траекторий, например для траекторий максимальной дальности. Поэтому предлагается иайти зависимость КВО от других параметров полета, используя опубликованные данные характеристик ЛА.

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

По заданным значениям параметров разброса вектора скорости ЛА в начале ПУТ Ко и КВО для траектории максимальной дальности известна ктах требуется определить зависимость КВО точки падения фрагмента ЛА от полной сферической дальности полета и угла наклона траектории ЛА в начале ПУТ.

В качестве исходных данных для модельных расчетов использовались характеристики Л А Трайдент И, заимствованные из работы Л. Гронлунда и Д. Райта. Матрица Ко в расчетах принята единичной, что соответствует сферической модели рассеивания по скорости отделения фрагмента в конце АУТ.

Величина ктях принята равной 100 ед. Это приводит к определению искомой зависимости КВО в процентах от ктах . Например, при КВО=145 истинное значение искомого КВО рассчитывается по формуле КВО = 1.45 х ктах, где ктлх выражено в метрах.

Основные допущения следующие. Земля предполагается сферической. Вращение Земли не учитывается. Влияние атмосферы на активном участке траектории (АУТ)

не учитывается. Влияние атмосферы на пассивном участке траектории (ПУТ) учитывается в соответствии с моделью торможения для фрагмента конусовидной формы со скругленной вершиной. Граница плотных слоев атмосферы принимается равной 90 км. Случайные возмущения, приводящие к рассеиванию точек падения фрагментов ЛА, моделируются нормальным распределением разброса вектора скорости ЛА в начале ПУТ с нулевым математическим ожиданием и ковариационной матрицей а1Ко > где матрица Ко задана, а а1 - скалярный параметр. Распределение указанных случайных возмущений полагается одинаковым для всех допустимых траекторий ЛА. Величина КВО ктях для траектории максимальной дальности известна.

Далее описывается алгоритм оценки кругового вероятностного отклонения.

Для удобства обозначений принято: нижний индекс «0» у некоторого кинематического параметра означает, что этот параметр берется на начало ПУТ, а нижний индекс «а» говорит о том, что это параметр относится к моменту входа в плотные слои атмосферы.

В расчетах используются следующие вспомогательные системы координат (СК):

ха, Уа% 2о — произвольная инерциальная, в которой задаются данные на начало ПУТ. В расчетах использована абсолютная геоцентрическая СК (АГСК).

Ха,Уа<га — инерциальная, связана с точкой Оа пересечения невозмущенной траектории фрагмента ЛА с границей плотных слоев атмосферы. Ось Оау„ направлена по внешней нормали к границе атмосферы, ось Оаха ортогональна оси Оауа и направлена в плоскости невозмущенной орбиты в направлении движения фрагмента ЛА. Ось 0Я2„ дополняет систему до правой.

В, 5 — плоская СК на поверхности Земли, связанная с точкой падения фрагмента ЛА при невозмущенном движении. 5 — отклонение по дальности в плоскости невозмущенной орбиты, В — отклонение по боку.

Некоторые обозначения.

V = — вектор скорости.

г = (х,у,г)Т — вектор положения.

ва = а,гс1д(ууа/и1а)— угол входа в атмосферу.

Ь — Ь (|г>„|, — вп)— путь по Земле при торможении фрагмента в атмосфере. Ьу = дЬ/,ЭК1, Ь0 = дЬ/д {-0а) .

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

С = дц/дУа — матрица баллистических производных размера 5x3.

А = — матрица частных производных размера 2x5, да = (ха, га, Уьа,уга)т,

£> = (Ъв)т.

Разработанный алгоритм включает в себя следующие три этапа.

Этап 1. Подготовка исходных данных для вероятностного анализа в результате моделирования дискретных пучков траекторий ЛА на ПУТ, соответствующих различным задаваемым с шагом 500 км значениям сферической дальности полета и различным значениям угла бросания в начале ПУТ. Моделирование проводится с помощью разработанной для этой цели программы. На этом же этапе используется программа для

определения значений функции Ь (|и„|, -0а) для смоделированного набора значений скорости и угла входа в атмосферу.

Этап 2. Определение параметров Атщ, АП1ах для смоделированных на этапе 1 траекторий. Это наименьшее и наибольшее собственные значения (2х2)-матрицы АСКп(угАт.

Этап 3. Для каждой смоделированной траектории искомое КВО к определяется по формуле:

к = ст\/Атах/(7), где / (7) — КВО для двумерного нормального закона с плотностью

р(х'у)=2^вфН(12+ 7)}'

а 7 = . Для нахождения / (7) используется процедура, описанная в первой главе и реализованная как функция программы.

Параметр а определяется по известному КВО ктах для траектории максимальной дальности. Для этого сначала выполняются первые два этапа алгоритма. Затем вычисляем

^_ ктах

\/Атах/ (7)

В этой же главе описывается алгоритм определения матрицы С баллистических производных. Для его реализации необходимо рассчитать эту матрицу в СК ха,уа, га. Далее подробно описано получение баллистических производным путем дифференцирования в АГСК по Ц/о и у2о векторного интеграла площадей, условий входа в атмосферу и интеграла энергии. После этого получим 5 линейных уравнений относительно шести неизвестных производных. Первую из них, т.е. предлагается определить, используя численные метод, основанный на методе наименьших квадратов. Тогда для остальных пяти производных получаем систему из пяти линейных уравнений, решаемую аналитически.

Далее следует описание моделирования траекторий на пассивном участке, исходными данными для которого являются требуемая сферическая дальность полета 1Сф и определяемые ниже параметры конца активного участка траектории - дальность и высота полета: 1К к И.к. Требуется определить в векторной форме радиус вектор и скорость полета на момент начала баллистического полета, а также координаты точек падения в целевой системе координат, соответствующие выбранным начальным условиям (НУ). Величина и наклон вектора скорости на момент окончания АУТ варьируются с помощью датчика случайных чисел в определенных пределах.

В сферической системе координат дальность баллистического полета Iе определяется параметрами конца АУТ, выступающими в качестве начальных условий баллистического полета (НУ БП): 1б = где вк — наклон вектора скорости Ук в конце АУТ (в точке К); Лд- — высота конечной точки АУТ. Если требуется определить значения перечисленных параметров для случая полета на заданную дальность (что требуется в настоящем исследовании), то возникает краевая задача, имеющая аналитическое решение (исчерпывающе подробно эта задача исследована Д.А. Погореловым).

Для исследования влияния вариаций начальных условий БП на рассеивание точек падения параметры скорости и положения центра масс ЛА в точке окончания полета на АУТ пересчитываются из плоской сферической системы координат (СК) в пространственную СК. Для удобства решения рассматриваемой задачи на этапе расчета координат точек падения целесообразно использовать «абсолютную геоцентрическую систему координат».

Затем производится расчет координат точки падения и скорости полета в ней, соответствующих НУ У/с и г/с (как в случае номинального, так и для возмущенного полета на ПУТ).

Описывается процесс атмосферного торможения, который существенно зависит от аэродинамических свойств фрагмента ЛА. В настоящей работе в качестве фрагмента рассматривается осесимметричное твердое тело конической формы со скругленной вершиной.

В итоге были произведены расчеты КВО в зависимости от угла бросания с шагом 1000 км по дальности полета. Результаты расчета для дальности 8000 км представлены в таблице 1.

Можно отметить два обстоятельства. Во-первых, видно, что для «настильных» траекторий, характеризующихся малыми углами бросания, рассеивание фрагментов возрастает в разы и может оказаться неприемлемым. Это возрастание носит «взрывной» характер, о чем свидетельствует также табл. 1. Сравнивая эти результаты с результатами, полученными с использованием предложенных моделей, но без учета аэродинамического торможения, можно сделать вывод о том, что отмеченный «взрывной» эффект является следствием именно аэродинамического торможения.

Таблица 1: Расчет для дальности 8000 километров

ек у/^тах 7 /(7) к.

15 15.88 6867.0 11453 0.0017 0.676 58.1

25 25.63 6683.6 4967 0.008 0.681 25.4

35 35.39 6754.2 4386 0.01 0.682 22.4

45 45.18 7096.3 4558 0.01 0.682 23.3

55 55.01 7804.3 5001 0.011 0.683 25.6

Во-вторых, можно заметить существование критического угла бросания, выше ко-трого функция КВО практически не меняется во всем диапазоне допустимых дальностей. В рассмотренном примере это критическое значение равно примерно 25 градусам. При углах бросания ниже этого критического значения рассеивание фрагментов резко возрастает. Более точное определение критического значения угла бросания требует использования более точных моделей движения ЛА и более точных моделей внешних факторов (форма Земли, свойства атмосферы).

ОСНОВНЫЕ РЕЗУЛЬТАТЫ, ВЫНОСИМЫЕ НА ЗАЩИТУ

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

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

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

4. Решена прикладная задача вероятностного анализа рассеивания точек падения фрагментов летательных аппаратов.-Использован алгоритм п. 2.

Публикации в журналах из перечня ВАК

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

2. Кан Ю.С., Травин A.A. Программный комплекс вероятностного анализа систем с кусочно-линейной структурой // Электронный журнал "Труды МАИ 2014, № 72.

3. Гончаренко В.И., Кан Ю.С., Травин A.A. Математическое и программное обеспечение анализа рассеивания точек падения фрагментов летательных аппаратов // Электронный журнал "Труды МАИ 2012, № 61.

Публикации в других изданиях

4. Кан Ю.С., Травин A.A. Применение метода линеаризации к решению задачи анализа движения КА. / 15-я международная научная конференция "Системный анализ, управление и навигация 2010, С. 152-154.

5. Кан Ю.С., Травин A.A. Вероятность поражения защищенной цели. / 16-я международная научная конференция "Системный анализ, управление и навигация 2011, С. 152-153.

6. Травин A.A. Вычисление квантильного и вероятностного критерия для систем с кусочно-линейной структурой / Научные труды Международной молодежной научной конференции "XL Гагаринские чтения". М.: МАТИ. 2014, С. 182-183.

7. Кан Ю.С., Травин A.A. Алгоритмы вероятностного анализа систем с кусочно-линейной структурой / 19-я международная научная конференция "Системный анализ, управление и навигация 2014, С. 120-122.

Программы, зарегистрированные в реестре программ для ЭВМ

8. Кан Ю.С., Травин A.A. Вычисление квантильного критерия нормы двумерного гауссовского вектора с заданной точностью. Свидетельство о государственной регистрации программ для ЭВМ № 2013616158 от 27.06.2013 г.

Множительный центр МАИ (НИУ) Заказ от/?. 0& 2015" г. Тираж ЮО экз.