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

Труды МАИ. Выпуск № 87 www.mai.ru/science/trudy/_

УДК 51-71

О решении задачи корректирования скалярного терминального состояния летательного аппарата при произвольном распределении

мультипликативного возмущения.

Игнатов А.Н.

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

e-mail: alexei. ignatov1@gmail. com.

Аннотация

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

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

Введение

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

Среди работ по коррекции летательных аппаратов с вероятностным и

УаЯ-критерием выделим обширную монографию [5]. В ней для поиска оптимального

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

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

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

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

Пусть - величина ошибки по какому-либо (одному) параметру траектории

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

корректирующего воздействия, ^ - параметр, характеризующий влияние корректирующего воздействия на величину ошибки, гх > 0 . Тогда г2 - величина ошибки после проведения коррекции - определяется по формуле

г2 = ^ + (1 + X).

Будем считать начальное состояние ^ случайным, = Х0 . Предположим, что случайные величины Х0 и х1 одинаково распределены, независимы, центрированы и имеют распределения с плотностями /(х) и / (х) соответственно. Величина ^ считается заданной. Также предположим, что до и после проведения коррекции состояния системы ^ и ^ могут быть точно измерены. Предложенная модель, с одной стороны, является частным случаем модели корректирования летательного аппарата, рассмотренной в [5], так как рассматривается только одна коррекция. А, с другой стороны, является ее обобщением, так как в [5] исследовалось лишь нормальное распределение помех.

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

Вместо задачи

zi) = argmax \<ф),

и1( z1)

где под записью щ (z) понимается, что управление щ ищется как функция от состояния системы z1, и дополнительно предполагается, что управление щ (z1) ищется в классе непрерывных функций, которая решалась в [5], будем решать задачу

uv(sa,si,...,%,s^+i) = arg max z2\<q>\

u1( s0's1'"-'sN 'sN +1)

где под записью щ(s0,s,-•,^,s^+1) понимается, что управление щ ищется как функция от сегмента (интервала, полуинтервала) s

So =(-«, z1), S1 =[z1, z 2), S2 = [z 2, z3),..., Sn =[zN, zN +1), Sn+1 =[-z\+x), в который попадает состояние z , где N + 2 - число сегментов разбиения, то есть управление представляет собой кусочно-постоянную функцию, а zN+1 = -z1 . Таким образом, управление щ (s0, s,-•, ^, s^+1) представляется в виде

щ1(s0, s1,., sN, sN+1) '

щ0, щ1,

щ1 , z1 G sN+1,

(1)

где и\ - некоторые действительные числа, I = 0, N +1.

Выбор сегментов

Компоненты гг , составляющие сегменты разбиения могут быть вычислены, например, по формуле

z1 ^ so-

z1 е s1

г'+1 = £ _ ^ = £ (1 _ И \ (2)

N V N)

где 1 = 1, N. При выборе гг по формуле (2) сегменты «,..., ^ имеют равную длину

н = г»1 _ ✓ = г1 Г1 _ 211_ г1 Г1 _ 1 = _2 £1

V

N

V

N ) N

Вычислим вероятность попадания состояния £ в сегменты бесконечной длины Р({Хо е «,} + {Х е ^+1}) =

= Р(Хо е «о) + Р(Хо е «N,1) _Р({Хо е «о,Хо е «N,1}). (3)

Поскольку сегменты « и по построению выбраны непересекающимися, то последнее слагаемое в (3) равно нулю, следовательно,

Р({Хо е «о},{Хо е «N,1}) = Р(_® < Хо < £Х),

, Р(_21 < Хо < +Х) = 1 _ 21) , 21), (4)

где

х

х)= | /{)(г)йг

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

г1 =1пДг : 1 _^о(_г), ^(т) >а,г< о}. (5)

Таким образом, сегменты при объединении составляют все множество действительных чисел, при этом если N^да, то длина сегментов «,«2,...,^ конечной длины стремится

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

Сведение исходной двухшаговой задачи к набору одношаговых задач Согласно формуле полной вероятности [8] получаем

def N+1

Pp(ui(% Si, — , SN, SN+i)) = P(| Z2 < () = Z2 < Z1 e S) :

i=0

N+1

= ep(| z1 + ^So,S1,.,sn,sn+1)(1 + X l)\<<P,z1 G st) ^ max . (6)

i=0 u1(s0's1'-"'sN ,sN+1 )

Учитывая представление управления щ(s0,s,-•,^,s^+i) в форме (1), задача (6) принимает вид

N +1

p,(ul(so,sn,sn+1)) = ZP(| z1 + (1 + Xl)|<(p,z1 G s) ^ max . (7)

i=0 u^u^^u^1

Поскольку управления щ0,щ1,...,u1N+1 друг с другом не связаны совместными ограничениями, а вероятности, стоящие под знаком суммы в (7), по определению неотрицательны, то для поиска оптимальной стратегии ~(s0,s1,.,sN,sN+1) , которая

представима в виде

Up(S0, S1, — , SN , SN+1) <

Z1 e S0' zi e Si,

~zi e SN+1,

нужно решить задачи

~ = arg max P(| Zi + ^1^1(1 + Xi) |< p, Zi e st), (8)

u

при 1 = о, N, 1 , а оптимальное значение критерия в классе кусочно-постоянных управлений равно

N ,1

РД~Д % SN, «N,1)) = ЕР(1 £1+ (1 + Х1) < £1 е ).

Детерминированный эквивалент

Рассмотрим подробнее 1 -е слагаемое функции Рр(и1(«о,«1,...,sN,«у+1)). При и[ =о

имеем

Р = Р(|£ , !1и1 (1, Х)|<р, £1 е ) =

= Р(| £1 |< £1 е ) = Р(| Хо |< р, Хо е )= | /о(УоЖ (9)

[_р,р]п«г-

При и{ > о получаем

Р1 (и[) = Р(| £ , гхи\ (1, Х) |< р, £ е «) = Р(_р < £ , (1, Х) < р, £ е «) =

= Р(_р < Хо , Х<и[ (1, Х1) < р, Хо е «г) =

( р_Уо х

= Р

р_ Х

о _ 1 < Х1 _ 1,Хо е I = /(Уо)

Хи

1 1

I /1( У^У1

_р_ Уо

йуоУ (10)

При и{ < о получаем

Р1 (и1) = Р(| £ , гхи\ (1 + Х) |< р, £ е «) = Р(_р < £ , гхи\(1 + Х) < р, £ е «) =

= Р(_р < Хо , Хй(1, X) < р, Хо е «г) =

о

£

V Х1и1

^- <Р-У0 1

= Р

Р-^ - 1< X <-РР^0 - 1, X е

V ¿1 и1

Щ

\/0( У0)

1 1

I /(У^

Р-у0 -1

^0

(11)

Таким образом, функция Рр(и^ )) представима в виде

N+1

Рр (и1(% 5N, 5N+l)) = ЕРг (и1)'

1=0

где

def

Рг (Щ) =

и

Ре', Р (и1), и Р- (и1), и

= 0, >0, <0,

где Р0', Р'+ (и[), Р1(и') определяются по формулам (9), (10), (11) соответственно. Для

решения задачи (8) при некотором фиксированном 1, 1 = 0, N +1 необходимо решить

задачи

и+ =агвзир Р+ (и1), Р+ = зир Р+ (и1),

и[ >0

и' >0

(12)

и- = аг§ эир Р- (и1), Р-' =эир Р- (и1).

(13)

Таким образом, решение задачи (8) выражается через решение вспомогательных задач (12), (13) следующим образом

ир = <

и-, Р- =тах{Р-, Р0', Р+ }, 0, Р- = тах{Р-, Р, Р+ }, и+, Р+ = тах{Р-, Р0', Р+ }.

Однако на практике аналитическое вычисление значения функции Р (и') в некоторой

5

V чи1

У

<

и <0

и. <0

фиксированной точке и1 , то есть вычисление повторного интеграла, может быть

9

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

Решение задачи поиска оптимального управления при помощи дискретизации вероятностной меры Проведем некоторые преобразования функции Pp(ui(s0,si,.,sN,sN+i)). Поскольку

вероятность произведения двух событий A и B меньше вероятности одного из них

P(AB) < P(B), P( AB) < P( A),

то имеют место следующие неравенства

P(| Zi + i^O + Xi) |< p,Zi e s0) < P(Zi e s0), (14)

P(| zi + ?iuf+1(1 + Xi) |< p,zi e stf+i) < P(zi e s„+1). (15)

Учитывая (5), (7), (14), (15), получаем

N+i

Pp(Ui(So,S,...,Sn,Sn+i)) = £P(| Zi + ii«i(i + Xi) |< p,Zi e s) <

i=0

N

< EP(| Zi + t U (i + Xi) |< p, Zi e st) + av i=i

Воспользовавшись методом средних прямоугольников [9] аппроксимируем i -е

слагаемое в (7), 1 = 1, N

Р(| 71 + гхи\ (1 + X,) |< р, 71 е 5,) = Р(| X + (1 + X,) |< р, X е 5,) ;

¥0

V + 7'

Р

„1 1 „1+1

г + г

2

+¿Щ (1 + X)

<р

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

РаГ°Х<5^ 51,., SN, SN+1)) =

N {71 + '

Ж

1=1

Р

„1 I „1+1

г + г

+ 1хи\ (1 + X)

<р

являющуюся аппроксимацией (7), и поставим задачу максимизации

рpappr°X(ul(5o,5l,...,5N,5N+l)) ^ тах.

и1

(16)

Для поиска оптимальных стратегий в (16) дискретизируем случайную величину X

Пусть х\ , к = 1, К - реализации случайной величины X , сгенерированные согласно плотности или функции распределения случайной величины X и

упорядоченные по возрастанию. Определим меру этих точек как р\ =1/К, к = 1, К • Составим случайную величину X со значениями х\ и вероятностной мерой,

сосредоточенной в этих точках Р{Х = х\} = рк.

Таблица 1. Таблица распределения случайной величины X •

2

2

2

%!2 К 1

р р1 Р2 К\ Р\

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

заданные числа. Подобные ограничения на управление можно встретить, например, в [11]. Таким образом, используя вместо непрерывной случайной величины X ее дискретный аналог ~, получаем задачи

К = агБ тах Л/0

иНи1о«"иир]

^+1'+ ^ ( р

I +I

+ (1 + X)

(17)

' = 1,N , и приближенное значение Р^(щ(я0,^,•..,^,^+1)) функционала вероятности

Р(и1( % ^ N, SN +1))

N

р

'=1 V

2

I +1

+ (1 +

(18)

Оценим сверху максимальное по модулю значение реализаций случайной величины

гхи\ (1 + X)

I (1 + ~1) 1=1^ II и1 II1 + ~ |<| I тах{| иь„ |, | и^ |} 11 + ~ <

<| | тах{| | иир |}(1 + тах{| % |. | % 1 |}).

2

2

2

Зафиксируем некоторое I, I = 1, N и решим задачу (17), сведя ее к следующей задаче смешанного целочисленного линейного программирования

1 K

Kl k=l M,e[MlowMup]

tu(l + Xk) < (<p - 0.5(Z + z1+l))Sk+ (l -sKJ)Z,k = l,...,Ki, (20)

tu (l + xk) > (-y - 0.5(z1 + z1+l ))^k,, - (1 - ^k,j )Z, k = l,..., Ki, (21)

8Kl e{0,l>,k = l,...,Kl, (22)

где Z =| tl | max{| UioW |,| u^ |}(l + max{| xj |. | xK |}).

При фиксированном допустимом управлении щ значения переменных 5к^г показывают выполнение неравенств

Z1 + Zl+l ~ т1' W+l

^ I ^ 1 у — т -г N

-y--< tyu[(l + x) < <р - -

2

для реализаций : 1 - если неравенства выполняются, 0 - если неравенства не

выполняются. Отметим, что задача (19)-(22) схожа по структуре с решаемой в [12]

задачей смешанного целочисленного линейного программирования.

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

определяемой соотношениями (19)-(22). В [5] было показано, что у оптимального

управления имеется зона нечувствительности к некоторым значениям состояния ^ .

Предложим подобную структуру управления и в рассматриваемой нами задаче поиска

оптимального кусочно-постоянного управления.

Зафиксируем управление щ = 0 . Получим задачу целочисленного линейного

13

программирования

- 1 ^

Р, = — ^ тах, К1 к=1

0 < (р-0.5(+ гм))Зк_ , + (1 ,,к = 1,...,К:,

0 > (-р -0.5(+ г'41))^,, -(1 -8к,,)?,к = 1,...,К:,

^, е {0,1},к = 1,.,К,.

(24)

(25)

(26)

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

V- [ 0, Р = Р,

и = •, Л _

р [и*, Р > Р,,

(27)

где и* - решение задачи (19)-(22). Таким образом, приближенным значением Рр(ир(я0,^,.,^,^!)) оптимального значения критерия Рр(ир(з0,,^+1)) является

N

Рр(~р(э0, , 3Х+1)) ~ Рр(ир( 30, , зх+1)) = ^//0

где управление

N ^ г+1^

2 + 2

,=1

Р ,

ир (30 , , . , , +1 )

uр, е ^

ир, г е ^

uр, е ^

е ,

к, е +1

(28)

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

<

находятся при помощи формулы (27). Отметим, что управление для сегментов ^ и в (28) доопределяется при помощи управления на сегментах ^ и ^ соответственно.

Пример

Пусть X ~ Ы(0,0.82) , X ~ Ч0,0.52) , а также ^ =1, иир = -и1оте = 10 и р = 1.15 . Зафиксируем ах = 0.000177 . Решив уравнение (5), получаем г1 = -3 . Проанализируем зависимость оптимального значения критерия Рр(йр^0, ^,..., ^, ^+1)) от числа сегментов разбиения конечной длины и количества реализаций К0 и К случайных величин X и X соответственно. Также подставим стратегию (28) в (7), чтобы оценить качество полученной аппроксимации.

Таблица 2. Зависимость оптимального значения критерия А(и„(я0,^,...,^,^+1))

от К, N.

*1 N Рр(г1р(^0, . ", SN, ^N+1)) Píp(йv(s0, S1,., SN, SN+1))

5000 10 0,98594 0,96854

20 0,98658 0,96939

10000 10 0,98593 0,96841

20 0,98647 0,96939

50 0,98661 0,96961

75 0,98628 0,9826

150 0,98627 0,98272

15000 10 0,98634 0,96844

20 0,98671 0,96939

Как следует из таблицы 2, увеличение количества реализаций к для больших к не сильно влияет на приближенное значение критерия. Поэтому при N > 20 рассматривался только случай к = 10000. Увеличение количества сегментов разбиения влияет на точность получаемого решения: чем больше сегментов, тем решение, получаемое при помощи предложенной процедуры, точнее. Отметим, что помимо значения критерия, как точного, так и приближенного, интересен вид самого получаемого управления.

Рисунок 1. Вид оптимального управления при N = 150, к = 10000.

Как следует из рисунка 1, получаемое управление имеет зону нечувствительности, которое практически совпадает с зоной нечувствительности для точного решения, полученного в классе позиционных стратегий [5]. Для дальнейшего увеличения точности получаемого решения необходимо еще больше увеличить число N, что однако потребует большего числа времени и мощности компьютера. Результаты в таблице 2 получены с помощью пакета ILOG CPLEX [13] на персональном компьютере (Intel Core i5 4690, 3.5 GHz, 8 GB DDR3 RAM).

Заключение

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

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

1. Братусь А.С., Черноусько Ф.Л. Численное решение задач оптимальной коррекции при случайных возмущениях. // Журнал вычислительной математики и математической физики. 1974. Т. 14. № 1. С. 68-78.

2. Ярошевский В.А., Парышева Г.В. Оптимальное распределение

корректирующих импульсов при однопараметрической коррекции. // Космические исследования. 1966. Т. IV. Вып. 1. С. 3-16.

3. Малышев В.В., Красильщиков М.Н., Бобронников В.Т., Нестеренко О.П., Федоров А.В. Спутниковые системы мониторинга. - М.: Изд-во МАИ, 2000. - 568 с.

4. Седов Л.И. Механика в СССР за 50 лет. Том 1. Общая и прикладная механика. -М.: Наука, 1968. - 416 с.

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

6. Кибзун А.И., Хромова О.М. О коррекции положения стохастической системы по квантильному критерию // Труды МАИ, 2014, №72: http://www.mai.ru/science/trudy/published.php?ID=47323

7. Азанов В.М., Кан Ю.С. Оптимизация коррекции околокруговой орбиты искусственного спутника Земли по вероятностному критерию // Труды Института Системного Анализа РАН. 2015. Т. 65. №2. С. 18-26.

8. Гмурман В.Е. Теория вероятностей и математическая статистика. - М.: Высшая школа, 2003. - 479с.

9. Самарский А.А., Гулин А.В. Численные методы: Учеб, пособие для вузов. - М.: Наука, 1989. - 432 с.

10. Rice S.O. Mathematical Analysis of Random Noise. // Bell System Technical Journal. 1945. Vo. 24. №1. P. 46-156.

11. Богуславский И.А. Методы навигации и управления по неполной

статистической информации. - М.: Машиностроение, 1970. - 256 с.

12. Наумов А.В., Иванов С.В. Задача распределения инвестиций в развитие отраслей наземного космического комплекса // Труды МАИ, 2012, №50: http://www.mai.ru/science/trudy/published.php?ID=28677

13. IBM ILOG CPLEX V12.1. User's Manual for CPLEX. - International Business Machines Corporation, 2009. - 952 p.