О ЛИНЕАРИЗАЦИИ МОДЕЛИ ВОЗМУЩЕННОГО ДВИЖЕНИЯ В ЗАДАЧЕ ВЕРОЯТНОСТНОГО АНАЛИЗА РАССЕИВАНИЯ БАЛЛИСТИЧЕСКИХ ТРАЕКТОРИЙ Текст научной статьи по специальности «Физика»

Труды МАИ. Выпуск № 99

http://trudymai.ru/

УДК 519.248:681.51

О линеаризации модели возмущенного движения в задаче вероятностного анализа рассеивания баллистических траекторий

Васильева С.Н.*, Кан Ю.С.**

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

* e-mail: sofia mai@mail. ru **e-mail: [email protected]

Аннотация

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

Ключевые слова: метод линеаризации, вероятностный анализ, круговое вероятное отклонение, рассеивание баллистических траекторий

1. Введение

В настоящей статье рассматривается задача вероятностного анализа рассеивания концов баллистических траекторий на поверхности Земли. В качестве характерисики рассеивания используется круговое вероятное отклонение (КВО). КВО определяется [1] как одна из характеристик точности срельбы и представляет собой длину радиуса круга, вероятность попадания в который равна 1/2 при условии совмещения центра нормального распределения ошибок с центром круга. Таким образом, КВО есть квантиль уровня 1/2 случайной величины, характеризующей расстояние между возмущенной и номинальной точками падения. Предполагается, что отклонение от расчетной траектории обусловлено только случайными возмущениями начальной скорости. Они считаются малыми по сравнению с модулем номинальной начальной скорости. Считается, что расчетная траектория является участком эллиптической Кеплеровой дуги. Выражения из теории Кеплера для расчета возмущенных траекторий нелинейно зависят от вектора начальной скорости. Поэтому КВО не удается найти аналитически. Для преодоления этой проблемы предлагается модифицировать метод линеаризации, предложенный в работе [2]. В частности, предлагается линеаризовать указанную выше нелинейную зависимость. Линеаризация осуществляется путем выделения линейной части

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

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

гауссовского вектора с заданной точностью. Метод основан на построении сходящихся последовательностей верхних и нижних оценок квантили.

Ранее задача по вычислению КВО рассматривалась в статье [5], в которой изучалась зависимость величины КВО от сферической дальности номинального полета при учете возмущений вектора скорости. В этой статье метод линеаризации фактически был использован, но его математическое обоснование не было представлено. В настоящей работе такое обоснование приводится впервые. В общей постановке исследуется вопрос о точности метода линеаризации в задачах оценки квантилей функции потерь, зависящей от малых случайных параметров. В разделе 6 доказано, что ошибка оценки квантильного критерия с использованием линеаризованной модели имеет тот же порядок, что и малый параметр. В общем нелинейном случае квантиль можно оценить методом Монте-Карло с использованием выборочных оценок, свойства которых изложены в [6], где также рассмотрены задачи стохастического программирования с квантильным критерием, частным случаем которого является КВО. Среди зарубежных публикаций, близких к данной тематике, можно указать, например, работы [7-18], в которых исследованы задачи стохастического программирования с вероятностными ограничениями. На модельном примере производится сравнение оценок КВО, полученных с применением метода линеаризации и метода Монте-Карло. Важность проблемы оценки адекватности математических моделей целевого применения летательных аппаратов отмечена в [19]. В разделе 8 приведены результаты расчетов, которые свидетельствуют о том, что ошибка в определении КВО с использованием метода

линеаризации по сравнению с методом Монте-Карло не превышает 1,5 % в широком диапазоне исходных данных.

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

Рассматривается возмущенное движение материальной точки в центральном гравитационном поле Земли. Начальные условия движения характеризуются детерминированными исходными координатами на поверхности Земли, заданными радиус-вектором я, и случайным вектором скорости V = Ун + Д V, где V - заданный детерминированный (номинальный) вектор скорости, ДV - вектор возмущения скорости. Компоненты вектора ДV независимы и одинаково распределены по нормальному закону с нулевым математическим ожиданием и дисперсией ц2 , причем значение параметра ц считается малым по сравнению со значением VI.

Подобная модель возмущений использовалась в работе [5] в постановке задачи вероятностного анализа возмущенного баллистического полета. Считается, что

меньше первой космической скорости и в отсутствии случайных возмущений (т.е. при ДV = 0 ) материальная точка в соответствии с кеплеровской теорией будет двигаться по эллиптической траектории, которая в некоторый момент времени Тк пересекает земную поверхность в точке падения яь . Из-за учета случайных возмущений происходит рассеивание точки падения, которую обозначим через я . Вектор я является случайным из-за случайности ДV. В качестве характеристики размеров области рассеивания выберем круговое вероятное отклонение (КВО) к, которое определяется условием:

р(1К - **)=12, (1)

где р - вероятность. Величина к подлежит оценке.

Отметим, что зависимость Як от приращения скорости АУ является нелинейной из-за нелинейности аналитических зависимостей теории кеплеровского движения. Это приводит к негауссовости закона распределения вектора Як и нетривиальности задачи вычисления к. Единственным работоспособным инструментом здесь представляется метод Монте-Карло, основанный на прямом статистическом моделировании возмущенных траекторий и оценке к по полученной выборке значений случайной величины - . Схема расчетов по

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

найдена относительно легко. Численный метод нахождения квантили заданного уровня для нормы двумерного гауссовского случайного вектора предложен в [4].

Подобный подход к оценке КВО был использован в [5], но без оценки его точности. Данная статья как раз и посвящена данному вопросу, т.е. оценке точности метода линеаризации по сравнению с методом Монте-Карло применительно к задаче вероятностного анализа рассеивания концов баллистических траекторий. Расчеты производятся как с учетом, так и без учета вращения Земли вокруг своей оси.

3. Используемые системы координат

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

Начальными параметрами являются координаты точки старта я в АГСК, модуль номинальной начальной скорости VI|, азимут 8 и угол бросания в. Конец

траектории як тоже лежит на поверхности Земли. Вектор V отличается от V тем, что в нем учитывается вращение Земли вокруг своей оси. Далее для удобства вместо

азимута используется угол а = ^-8. Значение величины VII и различные значения

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

При учете вращения Земли точка падения смещается против вращения Земли по сравнению с её положением на невращающейся Земле. Угловую скорость вращения Земли будем считать постоянной. Введем инерциальную систему координат Ох" у* 2*, которая будет совпадать с АГСК в начальный момент времени. Оси О2 и О2" совпадают. Эта система координат совпадает с АГСК на всём протяжении полета, если вращение Земли не учитывается. Если вращение Земли учитывается, то преобразование любого вектора X из исходной системы координат в систему координат Ох"у"2* осуществляется следующим образом:

X = л*х *, л* =

СОБ у Бт у 0 - Бту собу 0 0 0 1

2пТ

у =-

Т

вращ.

где Т = 86 400 [е.] - время полного обращения Земли вокруг своей оси, Т - время

полета (будет определено ниже).

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

г г г

определим направления векторов правой тройки гх , е2 , гъ , задающих её оси:

, [0 01]т X я ,_я

е1 =

Г г г

3 Я3 ' 2 31

[001] XЯ

где а х Ь - векторное произведение векторов а и Ь , Я — 6371 [км.] - радиус Земли.

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

Поскольку в начальный момент времени системы координат Охуг и Ох*у"г" совпадают, то и координаты векторов в этих системах тоже совпадают, поэтому из системы координат Охуг можно сразу перейти к системе Ох'у'2'. Но при нахождении координат точки падения с учетом вращения Земли необходимо из системы координат Ох 'у 2 ' перейти сначала в систему координат Ох'у" г", а уже из неё в систему координат Оху2 , проделав соответствующие обратные преобразования.

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

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

соб а Бт а 0 X " = А "XА "= - Бта соБа 0 0 0 1

матрицу А перехода из системы координат Ох у г" в систему Ох"у"г". Эта матрица будет однозначно определяться по радиус-вектору точки начала движения и вектору скорости, заданным в АГСК:

ез =

я

е2 =

КхУ

1МГ

Х" = АХ\А =

в\ ~е2Хез->

йг Й)'

Для определения расстояния от конца номинальной траектории до реальной точки падения материальной точки введем систему координат Оху2, ось 02 которой проходит через номинальную точку падения и центр Земли. Матрица перехода может быть определена по радиус-векторам точки приземления и точки старта следующими соотношениями:

ез =

е —

Якк я

Як х я X — лх, л

х Я|Г

е2 — е3 х е1,

(

(^2)

(е3)'

В этой системе координат при расчете КВО разность номинального и возмущенного векторов по третьей компоненте не существенна, поэтому будем рассматривать отклонения в плоскости Оху. Отклонения траекторий возмущенного движения по дальности АЬ будем откладывать по оси Оу, а отклонения по боку АБ по оси Ох .

4. Генерация исходных данных для расчетов В данном разделе приведены основные соотношения теории Кеплера,

позволяющие в плоском случае по заданной максимальной дальности полета Ь

определяется пучок из m навесных траектории одинаковом скорости, характеризующихся различными значениями в}, j = 1,..., m , и задаваемыми

значениями Ц < Lmx дальности полета, где Lmax - максимальная дальность.

Значение угловоИ дальности полета F в радианах связано с дальностью полета L следующим соотношением:

IT Lj Г"

Fj = R~3' j =1^

По номинальным значениям дальностей определим номинальные точки падения в системе координат Ox "y"Z":

R'Lj = (К sin Fj, 0, R3 cos Fj f.

Определим оптимальный угол бросания, под которым должна быть запущена материальная точка, чтобы угловая дальность была равна Lmx , а скорость в начальный момент времени была бы минимальной [20, (11.9)] (траектории минимальной скорости, одинаковой дальности):

3nax = ^arcrtg

f h ^

tg:

2 Sin Fmax J

P _ Lmax

5 max T-) R3

- iirii-\\R ii

где h = 11 11R™", - радиус-вектор точки падения, R - радиус-вектор начальной

IK, I

точки. Поскольку точки старта и падения лежат на поверхности Земли, то параметр к равен нулю. С учетом этого выражение для вычисления угла бросания можно представить в более простом виде:

, 1 Г ^}

Зпах = Т агссгё |

2 I 2

Для вычисления углов бросания для остальных т -1 дальностей полета понадобится значение модуля начальной скорости. При вычислении вектора скорости с учетом вращения Земли необходимо к номинальному начальному вектору скорости V = )т прибавить вектор скорости точки старта на поверхности Земли V, обусловленный вращением Земли, и лежащий в плоскости, касательной к поверхности Земли в этой точке, параллельно плоскости экватора:

V» = V + у3.

Абсолютное значение минимальной номинальной скорости определяется по формуле [20, (11.5)], которая после учета того, что точки старта и падения находятся на одном расстоянии от центра Земли, примет вид:

Г . _ . _ , V

Яз

*1п2£тахс1ё-^ + 2ссв2 £тах

V 2 У

Для вычисления углов в}, у = 1, т -1 необходимо решить квадратное уравнение,

полученное из выражения [20, (11.5)], которое связывает между собой номинальную скорость, угол бросания и угловую дальность полета:

|2 Я \\У\ I2 Я

ву-+1 = 0, у = 1, т-1.

2 по по

Решением данного квадратного уравнения является пара корней. Для навесных траекторий справедливы следующие соотношения:

=^ ^ыд _ 4+№ А

2 4

4

в. = arctg

ША+Я

2 2п

V

} = 1, ю -1.

Определим компоненты вектора скорости у учитывая, что его направление задано углом бросания в] и азимутом 5 . Для удобства перейдем в систему

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

скорости Земли в точке старта УЗ =

_ 24Ао

Т

, где А0г - расстояние от точки старта до

оси Ог в АГСК. Тогда вектор У" имеет вид:

У ' '=

соБа Бта 0 - Бт а соб а 0 0 0 1

ГУ' + УзЛ

у:

У

^ У ' + У3)соБа + у ' Бта ^

—(У + У31Бта + У2 соБа

У

у

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

уравнений:

с^в, =

(V

+ У ) сова + К Бта

V

(V ' + У3) Бт а = V' соб а;

(V;+V3 )2)2)2^

Данная система является разрешимой относительно компонент вектора скорости V' поэтому удается получить аналитические выражения для них:

V ' = соБв соа\\ VI ¡-V3,

V' = сОБв, БШа

2 ]

Уз' = Б1П в, „. „

Вектор номинальной скорости в начальный момент времени в АГСК определяется по следующей формуле:

V» =( А' )ТК

Также при учете вращения Земли, точка конца траектории полета смещается против вращения Земли на расстояние, пропорциональное времени полета. Для учета вращения Земли по найденным значениям скорости и угла бросания определим время полета. Для этого воспользуемся формулами [20, (9.59, 9.61)], поскольку рассматривается движение по эллиптической траектории между точками,

лежащими по разные стороны от её фокуса:

Т = Т (*-щ + 2г]).

Здесь

<

2

<

Р = агсБш

1 - 2к

У1 = е1 сО8 Р1 , ТР

= х/1 - 4к (1 - к) совв , а = , к = НА,

^ V ' 1 2 (1 - к) 2;т0

где е - эксцентриситет, а - большая полуось эллиптической орбиты. Приведенные формулы вытекают из [20, (9.59, 9.61)] с учетом того, что концы всех рассматриваемых траекторий находятся на поверхности Земли.

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

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

5. Расчет КВО методом Монте-Карло.

У.=Ун+АУ., А^ПЩО,/^2), где I = 1, п - номер моделируемой реализации, /3 - единичная матрица размера 3 х 3. Найдем матрицы Д. , которые являются матрицами А для каждой реализации вектора скорости V, ' = 1, п . Перейдем в систему координат, где движение при фиксированной реализации вектора скорости является плоским и происходит в плоскости Ох1г1. В силу этого факта вектор скорости имеет лишь две ненулевые

компоненты: V"=(К" 0 V")г. Угол бросания для каждой реализации вычислим по

V' — ^

известным компонентам скорости: дг= агс1ап^ , I = 1,п . Определим угловые

К

дальности ^, I = 1,п по формуле [20, (7.35)]:

F = 2 arctg

tg3

_ l

v cos2 e,R3 \\V\I2 J

Как и отмечалось ранее, плоское движение осуществляется в плоскости Ох"г", поэтому координаты точки конца номинальной траектории имеют следующий вид:

R = (R3 sinF, 0, R3 cosF) Осуществим обратный переход к системе координат Oxyz

RkH=A*TATK_

Rkl=A*T~^Rli = \n.

Получим статистическую выборку Лг=||Я - реализаций случайной величины А = ||Я - , являющейся расстоянием от номинальной точки падения до

фактического места касания поверхности Земли. Построим вариационный ряд А" <...<Апп реализаций случайной величины А . Выборочную оценку функции квантили уровня р можно определить по формуле из монографии [6, с.171]:

где [□] - целая часть числа.

Найденная оценка функции квантили и будет являться статистической оценкой для КВО, если положить р = 1.

6. Обоснование метода покомпонентной линеаризации функции потерь

Как отмечалось выше, КВО к является квантилью уровня 1/2 распределения случайной величины - |. В общем случае квантилью уровня р распределения случайной величины И называется действительное число

И р = т1п р}'

где Р =Р(И<^) - функция распределения случайной величины И . Свойства квантилей изложены в монографии [6]. Рассмотрим случайную величину - в системе координат Оху!, в которой отклонение траектории возмущенного движения по дальности АЬ будем откладывать по оси Оу, а отклонение по боку АБ по оси Ох . Поскольку отклонения вдоль оси О! незначительны, то опустим их и будем рассматривать отклонения концов возмущенных траекторий от номинальной точки падения в двумерном случае в плоскости Оху . Тогда справедливо следующее:

R - rA

fAB1 ,AL ,

л/AB2 +AL2 . С учетом этого вероятностное ограничение (1) может

быть представлено в следующем виде:

f fAB 1 1

р <к

v vAL J J

- 1

Очевидно, что

f f ABl 1 f f ABl 2 1

P < к = P

v vAL J J v vAL J J

Поэтому КВО может быть найдено из условия

к2 =

fAB 1

Для дальнейших выкладок потребуются следующие результаты.

Лемма 1. Пусть jj - случайные величины, Р = Р(£<^), Р1 = p(j<^), причем Pv>PJ,Vp, тогда: [j]^ >[£]^,Vpе(о, 1).

Доказательство. Введем обозначение [j]p=vp . Тогда P(^<^p) = Р^ >= p . Отсюда следует [£]p= min {з: p, > p}< min Pl > p} = [j]p, то есть [£]p <[j]p . Лемма доказана.

Лемма 2. Пусть X и 7 - случайные величины, причем Y неотрицательна, / > 0 - действительное число. Тогда для любого действительного числа R > о справедлива двухсторонняя оценка:

р(х<р)<р{x, где ек=р{у>я) (2) Доказательство. Оценим вероятность р(X + ¡У <р) . Учитывая, что У > 0 получим р(X + ¡У<р) = р(X<р-цУ)<р(X<р) . Тем самым обосновано второе неравенство в (2). Построим оценку снизу. р(X + ¡У <р) = рл1 + рл2, где

РК1 =р(X <р-¡У, У <Я), РК2 =р(X <р-¡У, У >Я).

Построим оценки для обоих слагаемых. Для первого слагаемого имеем следующую оценку:

р >р( X <р-лЯ, У < Я ) = р( X <р-лЯ, У < Я ) + р( X <р-лЯ, У > Я )-р( X <р-лЯ, У > Я ) = = р(X <р- ¡Я)-р(X <р-лЯ, У > Я)>р(X < (р-¡Я)-ед Ря1 <р( X + МУ <р)<р( X <р)

Для второго: 0<РЛ2 <р(У >Я) = ея.

В итоге получаем неравенство: р(X < р - ¡Я) < р(X + ¡У < р) < р(X <р) + ек. Лемма доказана.

Лемма 3. Пусть X и У - случайные величины, причем У отрицательна, ¡> 0 -действительное число. Тогда для любого действительного числа Я > 0 справедлива двухсторонняя оценка:

р^^^р^ + ^У <р)<Р(x<р + лЯ) + 5R где 8Я =р(У <-Я). Доказательство леммы 3 аналогично доказательству леммы 2.

Теорема. Пусть х, Y - случайные величины, л - детерминированный параметр. Тогда для УЯ > 0 з константы SR, sR : SR, 0 при R ^да , для любого p g (0,1) справедливо неравенство

[х 1-лЯ <[х+лY ] p <[х ] р+^+лЯ-

Доказательство. Введем обозначение:

Y = Y- + Y+, Y- = min {Y, 0}< 0, Y+ = max {Y, 0}> 0.

С учетом этого: Р (X + ¡Y < р) = Р (X < р - л (Y + Y+)).

Тогда справедлива следующая двухсторонняя оценка вероятности:

р( X <р- ¡Я )-s <р( X < р - ¡Y+ ) < р (х < р - л (Y- + Y+ )) < р (X < р - л Y- ) < р (X <р + лЯ ) + ^.

Рассмотрим левую часть неравенства р(X<р-¡Я)-s <р(X+лY<р). Согласно Лемме 1 отсюда следует, что [X + лY]р < [X + лЯ]р+е = [X]р+Е + лЯ . Аналогичные

рассуждения справедливы и для правой части неравенства: из

р(х+лY<р)<р(х<р+лЯ)+бя следует [X+лY]p >[х-лЯ]^ = [х]р_^ -лЯ.

Теорема доказана.

Применим теорему для нахождения оценки квантили квадрата нормы

отклонения. Для этого найдем разложение первого порядка вектора

(АВ 1

в ряд

^ Г AB 1 _ „ ^

Тейлора по возмущению вектора скорости: = DAV + r, где D - матрица Якоби

VAL У

(матрица баллистических производных), БАУ - линейная часть разложения, г -остаточный член в форме Лагранжа, АУ = /л^,^и N(0,^),

AV = AAV = = r = °(л2)• Поскольку матрица A является ортогональной, то вектор g= имеет нормальное распределение: n(0,/3). Возведем в квадрат

{ AB Л {AB Л 2 С другой

норму вектора lAL J = AB2 +AL2 . стороны

{ AB Л

,AL j

DAV + r\f =||DaV|f + ||r\f + 2(DAV,r). Здесь DAV имеет порядок малости л, r

- порядок /и2. Исходя из этого можно определить порядки слагаемых, входящих в квадрат нормы вектора отклонения. Первое слагаемое ||ВЛК||2 имеет порядок

малости u2 , 2(ВАК, г) - порядок малости ¡3 , ||г|Г - порядок малости u4 . В соответствии с этим справедливо следующее представление:

ГЛв)2

л2 (X + ¡Y + ц2Z), где X, Y, Z - случайные величины одного порядка.

Оценим квантиль уровня а случайной функции X + ¡У + ¡л22. Пусть для неё выполнено условие теоремы

Р([Х]Р <Х+1и¥ + 1и22<[х]р + в)>оу£>о. Сделаем замену Х = Х + /лУ, тогда согласно

теореме справедлива оценка квантили

Гх] -/л2 К <\х + /л22~] <Гх] +/Л2К.

_ _р—1 _ _р _ j р+ЕК1

По теореме оценим левую и правую части неравенства. Для левой части неравенства справедливо следующее [х] ^ = [Х + /лУ\р 5 <[х]р_з

Аналогичные выкладки справедливы и для правой части неравенства:

Гх] = \Х + /лУ] =\х]

2

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

Ир2 ¡1 / 4х + / Л¡1 /

Учитывая, что сумма и2Я + /Я имеет порядок малости /и , справедливо следующее:

[Х]^я^Я2 + ^(л)-[Х + Л + Л ^ ]р -[Х]2 -

Отметим, что ек и 5К зависят от величины Я , и при её стремлении к бесконечности стремятся к нулю. Если выполнено условие р([x\р< X <[x]^ + е)> 0, Уе> 0, то функция квантили [X]р непрерывна по р в точке

р , отсюда [X] 1 -^Я2 X], , [X^ +ея2 ЧX]р при Я

• Ю.

7. Схема расчетов по методу линеаризации

Метод линеаризации позволяет вместо исходной нелинейной задачи (1)

решать задачу с линейной моделью к = шт |р: р(|| БАУ\ -р)> 11. Поскольку левая и

правая части неравенства, стоящего под знаком вероятности неотрицательны, то

к = шт |р: р(||БАУ\|2 -р2 )> -11. Вектор, стоящий под знаком нормы имеет нормальное

распределение с нулевым математическим ожиданием и ковариационной матрицей К = л2ББТ . Рассматриваемая задача сводится к задаче оценки квантильного критерия уровня 1/2 нормы двумерного гауссовского вектора с известными параметрами распределения. В работе [4] был предложен метод, позволяющий

строить такую оценку параметра к с задаваемой точностью. Метод основан на построении верхней и нижней оценок квантили двумерного нормального вектора с задаваемой точностью. Задача о нахождении оценки квадрата нормы случайного вектора сводится к нахождению малой полуоси эллипса, вероятность попадания в который для реализаций двумерного стандартного нормального случайного вектора составляет 1/2. Отношение малой и большой полуосей эллипса известно и равно отношению меньшего собственного значения матрицы К к её большему собственному значению. Для нахождения малой полуоси используется метод, предложенный в [4], который является модификацией метода дихотомии, позволяющей находить оптимальное значение аргумента функции с задаваемой точностью.

Преимуществом линеаризованной модели является то, что величина к оказывается пропорциональна величине /. Этот факт позволяет не пересчитывать КВО для каждого значения / , а находить КВО путем умножения имеющихся результатов на отношение нового значения малого параметра к его старому значению. Данное обстоятельство отмечено в [5] и позволяет существенно упростить расчеты и избежать трудоемких вычислений по методу Монте-Карло.

8. Результаты расчетов

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

Таблица 1.

Ь, 103 км 6 7 8 9 10 11 12

в/ 59.45 54.15 48.72 43.07 37.06 30.25 18.02

Результаты, приведенные в таблице 1, получены согласно формулам раздела 4. При

проведении расчетов с учетом вращения Земли в нелинейном случае использовался метод Монте-Карло, описанный в разделе 5. Для азимута равного п/2 были получены следующие данные:

Таблица 2.

м / с 100 10 1 0.1 0.01

6 216.9228 21.6937 2.1697 0.2170 0.0217

7 229.7218 22.9628 2.2977 0.2298 0.0230

8 245.9681 24.6262 2.4639 0.2464 0.0246

9 268.1990 26.8010 2.6802 0.2680 0.0268

10 301.7416 30.2043 3.0214 0.3021 0.0302

11 364.0618 36.4017 3.6387 0.3638 0.0364

12 581.4578 58.2774 5.8193 0.5820 0.0582

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

Таблица 3.

м / с 100 10 1 0.1 0.01

6 213.8201 21.3820 2.1382 0.2138 0.0214

7 226.8283 22.6828 2.2683 0.2268 0.0227

8 243.6363 24.3636 2.4364 0.2436 0.0244

9 266.7932 26.6793 2.6679 0.2668 0.0267

10 301.7822 30.1782 3.0178 0.3018 0.0302

11 362.2104 36.2210 3.6221 0.3622 0.0362

12 577.6564 57.7656 5.7766 0.5777 0.0578

Относительные ошибки (в %):

м / с 100 10 1 0.1 0.01

Ь, 103

6 1.4510 1.4575 1.4732 1.4730 1.4731

7 1.2756 1.2341 1.2967 1.3066 1.3080

8 0.9571 1.0778 1.1314 1.1297 1.1303

9 0.5269 0.4561 0.4582 0.4557 0.4555

10 0.0134 0.0865 0.1176 0.1060 0.1048

11 0.5111 0.4987 0.4578 0.4353 0.4314

12 0.6581 0.8859 0.7403 0.7446 0.7451

Результаты расчетов для азимута, равного ж/ 3, с учетом вращения Земли по методу Монте-Карло:

Таблица 5.

м / с 100 10 1 0.1 0.01

Ь,103Км\^

6 215.2422 21.5459 2.1547 0.2155 0.0215

7 227.7667 22.7443 2.2749 0.2275 0.0227

8 244.8734 24.5013 2.4512 0.2451 0.0245

9 267.0662 26.7093 2.6716 0.2672 0.0267

10 301.2899 30.2095 3.0195 0.3020 0.0302

11 359.9150 36.0017 3.5969 0.3598 0.0360

12 572.8517 57.3469 5.7327 0.5732 0.0573

Расчеты по линеаризованной модели:

Таблица 6.

м / с 100 10 1 0.1 0.01

Ь, 103

6 213.8208 21.3821 2.1382 0.2138 0.0214

7 226.8290 22.6829 2.2683 0.2268 0.0227

8 243.6373 24.3637 2.4364 0.2436 0.0244

9 266.7944 26.6794 2.6679 0.2668 0.0267

10 301.7837 30.1784 3.0178 0.3018 0.0302

11 362.2127 36.2213 3.6221 0.3622 0.0362

12 577.6624 57.7662 5.7766 0.5777 0.0578

м / с 100 10 1 0.1 0.01

Ь, 103

6 0.6648 0.7660 0.7723 0.7628 0.7632

7 0.4134 0.2707 0.2897 0.2884 0.2872

8 0.5074 0.5646 0.6088 0.6146 0.6161

9 0.1019 0.1118 0.1372 0.1455 0.1475

10 0.1636 0.1031 0.0537 0.0659 0.0671

11 0.6343 0.6061 0.6961 0.6721 0.6735

12 0.8328 0.7260 0.7598 0.7773 0.7837

Далее приведены результаты расчетов в случае, когда Земля предполагается неподвижной и учет скорости вращения Земли не производится. Приведем таблицы с результатами расчетов при азимуте равном п/ 2. Результаты расчетов по методу Монте-Карло:

Таблица 8.

м / с 100 10 1 0.1 0.01

Ь, 103 к^^^

6 212.8787 21.2897 2.1291 0.2129 0.0213

7 226.4652 22.6351 2.2642 0.2264 0.0226

8 243.2685 24.3822 2.4396 0.2439 0.0244

9 266.9554 26.7276 2.6726 0.2672 0.0267

10 302.1888 30.2391 3.0253 0.3025 0.0302

11 361.0871 36.1604 3.6136 0.3613 0.0361

12 575.0542 57.4996 5.7568 0.5756 0.0576

Расчеты по линеаризованной модели:

Таблица 9.

м / с 100 10 1 0.1 0.01

Ь, 103

6 213.8201 21.3820 2.1382 0.2138 0.0214

7 226.8283 22.6828 2.2683 0.2268 0.0227

8 243.6363 24.3636 2.4364 0.2436 0.0244

9 266.7932 26.6793 2.6679 0.2668 0.0267

10 301.7822 30.1782 3.0178 0.3018 0.0302

11 362.2104 36.2210 3.6221 0.3622 0.0362

12 577.6564 57.7656 5.7766 0.5777 0.0578

м / с 100 10 1 0.1 0.01

Ь, 103

6 0.4403 0.4315 0.4267 0.4320 0.4325

7 0.1601 0.2104 0.1815 0.1728 0.1716

8 0.1510 0.0762 0.1346 0.1227 0.1241

9 0.0608 0.1811 0.1746 0.1639 0.1656

10 0.1347 0.2017 0.2468 0.2269 0.2249

11 0.3101 0.1673 0.2353 0.2519 0.2495

12 0.4505 0.4605 0.3420 0.3550 0.3563

Результаты расчетов по методу Монте-Карло для азимута, равного ж/ 3:

Таблица 11.

м / с 100 10 1 0.1 0.01

Ь,10з3км\^

6 214.2608 21.4251 2.1418 0.2142 0.0214

7 227.2456 22.7334 2.2728 0.2273 0.0227

8 243.7343 24.4381 2.4427 0.2442 0.0244

9 265.5993 26.6108 2.6612 0.2661 0.0266

10 302.5890 30.2947 3.0326 0.3033 0.0303

11 364.0885 36.4789 3.6488 0.3649 0.0365

12 581.7320 58.2252 5.8214 0.5821 0.0582

Расчеты по линеаризованной модели:

Таблица 12.

м / с 100 10 1 0.1 0.01

Ь,103км\^

6 213.8204 21.3820 2.1382 0.2138 0.0214

7 226.8286 22.6829 2.2683 0.2268 0.0227

8 243.6367 24.3637 2.4364 0.2436 0.0244

9 266.7937 26.6794 2.6679 0.2668 0.0267

10 301.7827 30.1783 3.0178 0.3018 0.0302

11 362.2112 36.2211 3.6221 0.3622 0.0362

12 577.6584 57.7658 5.7766 0.5777 0.0578

Относительная погрешность вычислений:

м / с 100 10 1 0.1 0.01

Ь, 103

6 0.2059 0.2014 0.1703 0.1808 0.1818

7 0.1839 0.2228 0.1972 0.2097 0.2110

8 0.0401 0.3057 0.2616 0.2485 0.2470

9 0.4477 0.2569 0.2530 0.2445 0.2445

10 0.2672 0.3858 0.4888 0.5098 0.5119

11 0.5183 0.7118 0.7373 0.7367 0.7366

12 0.7052 0.7952 0.7766 0.7744 0.7824

Как видно из таблиц, относительная погрешность метода линеаризации не превышает 1,5% как с учетом, так и без учета вращения Земли. Отметим, что в [3] вращение Земли не учитывалось.

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

В статье обосновано применение метода линеаризации к решению задач вероятностного анализа рассеивания концов баллистических траекторий. Относительная погрешность по сравнению с методом Монте-Карло, возникающая при использовании линеаризованной модели при подсчете КВО, достаточно мала и не превосходит 1,5% в широком диапазоне дальностей полета и азимутов. КВО, полученное таким способом, линейно зависит от среднеквадратического отклонения компонент скорости /. Этот факт позволяет не пересчитывать КВО для каждого значения / , а находить КВО путем умножения имеющихся результатов на отношение нового значения малого параметра к его старому значению. Такой подход может значительно сэкономить вычислительные ресурсы, необходимые для расчета КВО.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (РФФИ), проект №18-08-00595.

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

1. Сергеев И.Д., Яковлев В.Н., Соловцов Н.Е. и др. Военный энциклопедический словарь Ракетных войск стратегического назначения. - М.: Большая Российская энциклопедия, 1998. - 634 с.

2. Васильева С.Н., Кан Ю.С. Метод линеаризации для решения задачи квантильной оптимизации с функцией потерь, зависящей от вектора малых случайных параметров // Автоматика и телемеханика. 2017. № 7. С. 95 - 109.

3. Алферьев В.Л. Свойства матриц частных производных на Кеплеровой дуге // Двойные технологии. 2011. № 4(57). С. 14 - 21.

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

5. Гончаренко В.И., Кан Ю.С., Травин А.А. Математическое и программное обеспечение анализа рассеивания точек падения фрагментов летательных аппаратов // Труды МАИ. 2012. № 61. URL: http://trudymai.ru/published.php?ID=35615

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

7. Kogan A., Lejeune M.A. Threshold boolean form for joint probabilistic constraints with random technology matrix // Mathematical Programming, 2014, no. 147, pp. 391 - 427.

8. Kogan A., Lejeune M.A., Luedtke J. Erratum to: Threshold boolean form for joint probabilistic constraints with random technology matrix // Mathematical Programming, 2016, no. 155(1), pp. 617 - 620.

9. Guigues V., Juditsky A. Nemirovski A. Non-asymptotic condence bounds for the optimal value of a stochastic program // Optimization Methods & Software, 2017, no. 32(5), pp. 1033 - 1058.

10. Barrera J., Homem-de-Mello T., Moreno E., Pagnoncelli B.K., Canessa G. Chance-constrained problems and rare events: an importance sampling approach // Mathematical Programming. Ser. B, 2016, no. 157, pp. 153 - 189.

11. Guigues V. Henrion R. Joint dynamic probabilistic constraints with projected linear decision rules // Optimization Methods & Software, 2017, no. 32(5), pp. 1006 - 1032.

12. Bremer I., Henrion R., Möller A. Probabilistic constraints via SQP solver: Application to a renewable energy management problem // Computational Management Science, 2015, no. 12, pp. 435 - 459.

13. Luedtke J. A branch-and-cut decomposition algorithm for solving chanceconstrained mathematical programs with finite support // Mathematical Programming, 2014, no. 146(1-2), pp. 219 - 244.

14. Minoux M., Zorgati R. Convexity of gaussian chance constraints and of related probability maximization problems // Computational Statistics, 2016, no. 31(1), pp. 387 -408.

15. Wim van Ackooij. Eventual convexity of chance constrained feasible sets. Optimization. // A Journal of Mathematical Programming and Operations Research), 2015, no. 64(5), pp. 1263 - 1284.

16. Wim van Ackooij, R. Henrion. (Sub-) Gradient formulae for probability functions of random inequality systems under Gaussian distribution // SIAM Journal on Uncertainty Quantification, 2017, no. 5(1), pp. 63 - 87.

17. Wim van Ackooij, Sagastizabal C. Constrained bundle methods for upper inexact oracles with application to joint chance constrained energy problems // SIAM Journal on Optimization, 2014, no. 24(2), pp. 733 - 765.

18. Xie W., Ahmed S. On quantile cuts and their closure for chance constrained optimization problems // Mathematical Programming, 2017, ser. B, pp. 1 - 26.

19. Гончаренко В.И., Кобзарь А.А., Кучерявенко Д.С. Идентификация параметров движения летательных аппаратов на активном участке траектории с использованием дискретного вейвлет-преобразования // Труды МАИ. 2011. № 46. URL: http: //trudymai .ru/published.php?ID=25995

20. Погорелов Д.А. Теория кеплеровых движений летательных аппаратов. - М.: Физматгиз, 1961. - 106 с.