Построение оптимальных траекторий в многомерных пространствах на основе физических моделей Текст научной статьи по специальности «Математика»

X МОДЕЛИРОВАНИЕ СИСТЕМ И ПРОЦЕССОВ

УДК. 519.95

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

С. Д. Субочев,

канд. техн. наук, вед. программист

Санкт-Петербургский государственный университет аэрокосмического приборостроения

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

We derive an expression for the radius of a small arc of the optimal trajectory in two physical models: the motion of a point with unit velocity and the propagation of a light ray. The principle of minimal propagation time in three-dimensional physical space is generalised to high-dimensional space and is used for the elaboration of the corresponding digital method of fire adjustment in contstructing the optimal high-dimensional trajectories. The optimality of the obtained trajectories is justified by a formal mathematical gradient method.

Введение

Пусть заданы некоторая положительная функция потерь или затрат от п переменных ^г) = F(x1 , х2,..., хп) > 0, а также начальная и конечная точки г^ = [х1Ь, х2Ь,..., хт]т и

т

= [х1е, х2е, ..., хте] соответственно.

Оптимальной траекторией является линия Ь, соединяющая эти две точки, такая, что криволинейный интеграл I вдоль этой линии от заданной функции минимален:

Ь: {1 = $F(r)dZ = шт|. (1)

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

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

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

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

Известны многомерные сферические координаты [4, 5]. Эти координаты задаются углами, которые являются обобщенными физическими аналогами углов широты и долготы.

В статье развивается ранее означенный и нами подход [6] с обобщением физических аналогий на многомерные пространства для поставленной оптимизационной задачи.

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

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

этой точки Дг при разложении в ряд Тейлора с точностью до 2-го порядка малости по времени Д^ записывается так:

Аг = [Аж1, Ах2, Ахп]Т = УAt + WAt2/2, (2)

где

V/ Г ^ ' 1Т

= г =Аі , X2,..., Хп] ,

W// V /г гг

= г =Аі, *2,...,

// тТ

(3)

(4)

— векторы первых и вторых производных координат по времени соответственно (и где приняты

обозначения х, = и х" = ! х).

1 dt 1 dt2

Вектор V по своему физическому смыслу является вектором скорости, а вектор W — вектором ускорения в многомерном пространстве. Наложим ограничение, что модуль вектора скорости всегда имеет единичное значение

(5)

Дифференцируя (5) по времени, получаем, что скалярное произведение векторов скорости и ускорения тождественно равно нулю:

d^

VI - (У, ^ - 0.

(6)

Условие (6) означает ортогональность векторов скорости и ускорения, при этом можно положить, что и в многомерном пространстве при постоянной скорости движения полное ускорение является только нормальным, аналогично как и в трехмерном.

Направим касательную ось 0т параллельно вектору скорости V на середине дуги, а нормальную ось 0п — по вектору ускорения W (рис. 1), распо-

■ Рис. 1. Малая дуга оптимальной траектории

30 X"

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

їье8 = [-АТ °]Т, гепа = [т °]Т.

(7)

Тогда отрезок прямой, соединяющий по оси 0т эти точки с координатами (7), является хордой дуги. Так как любое малое перемещение Дг с точностью до 2-го порядка малости раскладывается по двум ортогональным базисным векторам V и W, малую дугу траектории можно считать плоской, т. е. лежащей в плоскости рисунка 0тп. Радиус кривизны траектории Rк и модуль нормального ускорения связаны общеизвестным соотношением, которое при единичной скорости записывается так:

Н = |VII2 Л-1 - Л-1. (8)

При этом векторы скорости и ускорения в проекциях на координатные оси

У = [т, п]Т = [1, 0]Т;

W = [т', П]Т =[0, Л-1 ]Т,

(9)

(10)

откуда уравнения траектории относительно середины дуги с точностью до 2-го порядка малости при -Дт < t < Дт равны

Пт

^2Лк1/2-П0,

где высота параболического сегмента дуги П0 =Ат2Лк-1/2.

(11)

(12)

Вместо функции от т переменных F(x1, х2,..., хт) рассмотрим соответствующее ей сечение F(т, п) в плоскости 0тп. При фиксированных концах (7) элементарная дуга окружности полностью определяется радиусом кривизны. Определим оптимальный радиус кривизны этой малой дуги такой, что криволинейный интеграл

Дт

1(Лк) = / F[nт(т/Лк)] т(т/Л^т, (13)

-Дт

зависящий от параметра Лк, будет минимален, т. е.

1(Лкр1) = ш1ш {1(Лк)}. (14)

С точностью до 2-го порядка малости

Пт(т / Лл ) =т2Л^/2-П0 (15)

— уравнение траектории (в котором произведена замена переменной t = т),

т(т / Лк^т = . 1 +

*Пт

dт

dт =

.2 Л

1 + -

2Лк

dт (16)

— дифференциал длины дуги.

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

*Ьев = [0, -П0І , ^ = [Ат, 0] ,

(17)

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

Дт

1(ЛК) = / F[nт(т/ Лк)] т(т/ Л^т; (18)

ЦЛ^1) = ш1п{/(Лк)}.

(19)

В этой новой задаче, как и прежде, зафиксирована конечная точка. Для новой начальной точки = [0,-П0]Т зафиксирована координата т = 0, а координата п0, рассчитываемая как (12), является варьируемой по радиусу кривизны. При этом векторы скорости и ускорения, записанные в проекциях на координатные оси как V = [1,0]т и W = [0, Л-1]Т, имеют постоянные направления, независимые от радиуса кривизны. Для старой начальной точки гЬег = [-Дт, 0]т направления векторов скорости и ускорения варьируются в зависимости от радиуса кривизны, что существенно усложняет задачу оптимизации. Этим и объясняется переход ко второй постановке задачи, решение которой, строго говоря, является квазиоптимальным, так как новая начальная точка не зафиксирована по координате п.

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

дF дF

F(0 + т, 0 + п) = F + — т + —п +

Эт Эп

. д2¥ т2 +

Эт2 2 дтдп^ дп2 2 '

(20)

где численные значения функции и ее соответствующих частных производных берутся в нулевой точке в начале координат 0тп (т = 0, п = 0).

Положим, что Лк > 1 и запишем минимизируемый по параметру Лк интеграл в виде суммы двух интегралов

1(Лк) ~ 10(Лк) + 12(Лк),

(21)

которые берутся соответственно от нулевого и 2-го членов ряда Тейлора (20):

Ат

ї0( Л к) = / Рт(т / Л^т,

(22)

Ат

12(Лк) = [^ п(т / Лк)т(т / Лк)!т. (23)

0 Эп

Подставляя в формулы для интегралов (22) и (23) уравнения траектории (15) и дифференциала дуги (16), производя интегрирование по т и затем дифференцирование по Лк, из условия

Л.

dЛ„

= 0

(24)

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

ДГ -К)-1 ¥,

(25)

д¥

где Fn =— — производная по нормальному на-

1 Эп

правлению.

При решении уравнения (24) был учтен 3-й порядок малости Дт3 в интегралах 10 и 12, берущихся от нулевого и 2-го членов ряда Тейлора (20).

Остальные члены ряда Тейлора, имеющие более высокие порядки малости, отбрасываются.

Итак, чтобы найти Лкр\ надо задать координаты г и вектор единичной скорости V (или единичный вектор касательного направления т - V).

Если модуль некоторого вектора скорости не нормирован к единице: ||У^|| Ф1, то следует перейти к единичной скорости (или единичному вектору направления т)

У = Ун к -т.

(26)

Далее необходимо найти вектор нормальной составляющей градиента функции

Ур = УР -(УР,т)т,

(27)

рассчитать производную по направлению нормали

¥ = УР

п п

(28)

и, подставляя (28) в (25), определить оптимальный радиус кривизны

Д

°Р^

-1

¥ =||УРп|| 1 ¥.

V П/ || П ||

При этом единичный вектор нормали

П=|1УРЛ-1 УРп.

(29)

(30)

Движение по оси 0т является заданным, или вынужденным, а по оси 0п — варьируемым, выбранным из условия минимума (1). Результат имеет очевидный смысл: при увеличении радиуса кривизны слагаемое 10(Лк) = FДt (21) растет, так как увеличивается длина дуги Дt, а слагаемое 12(Лк) =

~ - Fп' Дt3/Лк убывает, так как дуга выгибается в сторону меньших значений функции, и при условии (29) достигается минимум.

Обобщение модели распространения света по оптимальной траектории на многомерное пространство

Определим радиус кривизны малой дуги оптимальной траектории на основе известного факта, что свет в оптической среде с переменной оптической плотностью распространяется из начальной точки в конечную по некоторой траектории, такой, что время распространения — минимально. Допустим, что свет распространяется по некоторой кривой Ь, приходя из некоторой начальной точки в конечную (рис. 2). Разобьем кривую Ь на п

п

малых отрезков Д1; так, что ^Д1; = Ь. Обозначая

1=1

модуль скорости света на малом отрезке ||^|| и учитывая, что время распространения света на отрезке Д^ = ||Уг||-1 Д1, запишем суммарное время распространения Т по линии Ь как предельный переход к соответствующему криволинейному интегралу:

■ Рис. 2. Элементарные участки траектории луча света

П^~ А 1

Т=*£ Й!=Ь||у'(,)|(31)

где ЦУсСОЦ — модуль скорости света в оптической среде как функция положения в зависимости от текущей длины I на траектории Ь. Запишем известное соотношение в виде

М| = FC-1(Z)|| ^||, (32)

где Fc(l) — абсолютный показатель преломления оптической среды; Ц^Ц — модуль скорости света в вакууме. Тогда, подставляя (32) в (31), запишем принцип минимума времени распространения света

Ь :ш1шН ^с<1>!11 = / ^с(1)!1, (33)

[ Ьу J Ь

где Ьу — некоторая варьируемая траектория; Ь — траектория луча света, вдоль которой «автоматически» достигается минимум интеграла от функции Fc(l) — переменного показателя преломления оптической среды. Из курса оптики известна соответствующая закону Снеллиуса [7] простая иллюстрация преломления луча света при прохождении через границу двух слоев с разными показателями преломления F0 и F1 (рис. 3, дополненный необходимыми построениями). Луч света при переходе границы Г уменьшает скорость со значения У0 до значения У1 (так как F1 > F0); пгр — нормаль к граничной поверхности, разделяющей слои с различной оптической плотностью (градиент УР направлен по нормали пгр); а0 — угол падения; а1 —

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

угол преломления; Ф0В0 — длина фронта падающей волны, а Ф1В1 — преломленной; точки М0 и М1 — середины этих фронтов соответственно. Прямые углы отмечены точками.

Правый крайний падающий луч за время Дt проходит малый путь

В0В1 = У^. (34)

Левый крайний преломленный луч проходит путь

(35)

Ф0Ф1 = У1Дt.

Очевидно, что через точки М0 и М1 можно провести дугу окружности касательно к векторам скоростей. Тогда радиус кривизны этой дуги

Лк = Ц0М0 = Ц1М1 = Да ДЬ01,

(36)

где точка Ц0 образуется пересечением продолженных прямых линий фронтов Ф0В0 и Ф1В1; Да — угол между прямыми Ц0 М0 и Ц0 М1; ДЬ01 — длина дуги, приходящей из точки М0 в точку М1.

С точностью до 1-го порядка малости заменим длину дуги ДЬ01 на путь крайнего правого падающего луча

ДЬ01 ~ДБ = У0Д t ~ У^01Д t.

(37)

Прямая Ц1В1 параллельна прямой Ц0-В0, тогда с точностью до 1-го порядка малости (см. треугольник Ф1В1Ф2)

где

Да = ДУ Дt / Ф1В1 ~ДУ Дt / Ф0В0, Ф0В0 = ДБ ^ а0 = ДЬ01 ctg а0;

У У ДУ = У0-У1 = -1^-

0

0

(38)

(39)

(40)

где приращение показателя преломления вдоль линии ДЬ

01

ДF = F1 - F0. Расписывая (38) далее, имеем

Да^

= VвДtF0 2 эт а0

ДF

ДЬ01 сое а0

(41)

(42)

Отношение ДF/ДЬ01 в пределе при ДЬ01 — 0 представляет собой производную по направлению Ь (по касательной оси 0т), которая равна проекции градиента на это направление:

ДF

дF

— = 11Ш

Эт —>0 Д Ь

= (УР,^ = ||У^| сов «0, (43)

'01

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

пНИ «0.

(44)

Тогда, подставляя (44) и (43) в (42), а (42) в (36), после несложных преобразований окончательно имеем

F.

(45)

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

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

Алгоритм поиска оптимальной траектории методом пристрелки

Входными величинами на каждый 1-й шаг являются местоположение текущей точки Г = [х1;, х2;, хП1 ]Т и заданное направление движения в виде вектора единичной скорости или в виде единичного вектора направления V; -т; = [х!;, х2;,..., х'п; ] . Кроме того, должны быть рассчитаны радиус кривизны Лк; и вектор нормали п; по формулам (27)-

(30).

Далее рассчитываются соответствующие выходные величины ,-го шага

Г;+1 = Г; + Л; Дп; + т; Дт^

(46)

где Дт; — малое перемещение по касательной (которое задается одинаковым для всех шагов);

Дп; = Дт^ 0;

— малое перемещение по нормали.

т;+1 = в1ш 0;п; + Сов 0;тI,

где

= агсэт

^ДтЛ

Лк

(47)

(48)

(49)

— угол раствора ,-й дуги траектории.

Таким образом, если в заданной начальной точке г0 - rbeg = [х1Ь, х2Ь,..., хпЬ ] задать и некоторое начальное направление т0, то текущая точка, перемещаясь по шагам, будет двигаться по определенной траектории, которая в общем случае пройдет мимо заданной конечной точки Геш(1 = [х1е, х2е,..., хпе ] .

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

Г; - Гк

< Лсф ^11Г

;+1

- Г

beg

> Л,

сф

(50)

2

где

Л = г _____________г

‘'сф eпd beg *

(51)

На этой сфере должны заканчиваться все пробные (или пристрелочные) траектории. При выполнении условия (51), решая уравнение

Лс2ф -|- Гbeg||2 = 0 (52)

найдем экстраполированную точку траектории Г;э(^ф), которая лежит на этой сфере. Траектория экстраполируется рядом Тейлора до 5-го порядка малости

гэ(0 ~г + Х

^тг т

^гтт!’

(53)

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

dt2

dt2 Лк так и для нечетных:

п

лк

(54)

= а3т

dt , dt3

т 6_г т

ЛТ dt5 Д4

л,3 п2’ 1^5 л4 , (55)

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

Далее уравнение (52) решается методом последовательных приближений

^+1 = ^ + [Лс2ф - Бкв(^ )]/Бкв(Ь ), (56)

где у - номер итерации;

Бкв(^) = ||Гэ(^) -Гbeg||2 (57)

— квадратический путь экстраполируемой точки;

Бкв(t^) = ||Гэ(^) - ГЬе^|гэ ) (58)

— производная квадратического пути;

drэ Д dmr ^m-1

Гэ (^; ) - - = Х----------------

1 dt (т-1)!

(59)

— производная экстраполируемого радиуса-вектора.

Определив конечную точку некоторой пристрелочной траектории, лежащую на сфере Гпр -Гэ(£сф), где ^ф — время на последней итерации, при котором экстраполируемая точка достигает сферы, рассчитаем направляющие косинусы прямой линии, соединяющей начало траектории и точку у-й пристрелки Гпру

е = Л 1(г — г )

пр/ хф' пр/ л■beg/’

(60)

при этом направляющие косинусы прямой линии, соединяющей начало и конец требуемой оптимальной траектории:

■'тр

Лсф (reпd rbeg ).

(61)

Этим векторам направлений епр/ и етр/ некоторым образом соответствуют векторы начальных касательных направлений тпр/ и ттр пристрелочной траектории Тпр/ и требуемой оптимальной траектории Ттр (рис. 4, где траектория Ттр условно наклонена к нам, а Тпр/ — от нас). Вектор касательного направления оптимальной траектории ттр/ пока не известен. Следовательно, надо в вектор 0тп вносить такие поправки, чтобы вектор епр/ совместился с вектором етр.

Это проделаем таким образом. Рассчитаем поправку на плоскости П, натянутой на два вектора е„„: и е„.

пр/

тр

Аепр/ епр/ етр.

(62)

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

пр/

I (тпр/, епр/) .

(63)

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

|-/тр

| етр) ||

пр/

(64)

Ненормированный (/ + 1)-й вектор пристрелки сформируем как

-пр/+1

= Стр = Спр/ + Аепр/ .

(65)

Сп

Рис. 4. Соотношение векторов в методе пристрелки

т

Тогда с учетом подобия треугольников, образованных парами векторов ттр, тпру и стр, спру, нормированный у + 1-й вектор пристрелки сформируем так:

тпр/+1 = ттр = тпр/ + (ттр , етр ) Аепр/ .

(66)

Условием окончания итераций пристрелки (66) можно выбрать уменьшение угловой ошибки на сфере ниже минимального заданного предела

||Депр/|| <8ш1ш (67)

либо задавать заведомо большее количество циклов пристрелки, чем необходимо для выполнения условия (67). Чтобы не накапливались погрешности, единичный вектор направления, рассчитываемый по формулам (47) и (49), следует нормировать на каждом шаге вновь по формуле (26). Заметим, что, строго говоря, точка начала координат 0тп, относительно которой рассчитывается радиус кривизны, нам неизвестна, а фактически известна точка начала дуги с координатами гЬед = [0, -п0]т. Как показало промежуточное моделирование (результаты которого здесь не приводятся), процесс итерационного поиска точки начала координат от точки начала дуги через радиус кривизны не всегда приводит к необходимому уменьшению интеграла. Поэтому в формулах (29) и (45) следует рассчитывать постоянное значение функции и ее производной по нормали относительно доступной точки начала дуги по приближенным формулам

2(0,0) « 2(0,-^), 2(0,0) « 2(0,-п0). (68)

Вычисление криволинейного интеграла вдоль оптимальной траектории

Этот интеграл требуется не только для проверки различных траекторий на оптимальность, но также и для оценивания затрат или потерь наблюдателем. При выводе формулы (29) при разложении в ряд Тейлора были учтены варьируемые составляющие - в интегралах 10(Лк) и 12(Лк), которые оказывают наибольшее влияние на минимизацию интеграла вдоль малой дуги. Однако при большей протяженности траектории следует учесть и составляющие 11(Лк) и 13(Лк), которые вносят хотя и менее варьируемый — более постоянный, но зато существенный вклад в минимизируемый интеграл. Зафиксировав найденный оптимальный радиус (45), определим длину дуги как Дt , = Лк ;0; , где 0; — угол раствора дуги. Далее, опуская индекс I, возьмем четыре учтенные до 3-го порядка малости составляющие интеграла вдоль дуги окружности:

Дt

10 = | Fdt = FДt; (69)

0

11 = Tт(t)dt = ^^(1 - сов0);

дт дт

Эт

(70)

At

I* эр э 2 о

12 = [т-п,.^t = — (ЛКДtсов 0-Лк в1ш 0); (71) 0 дп дп

Э2.

13 = [ Э 2 Tт(t)dt = Э_р.^ дt-в1ш20 |. (72)

3 J Эт2 о Л а„2 I « 1 ' '

4 Эт2

Составляющая 10 имеет 1-й порядок малости. Несложно показать, что составляющие 11; 12 и 13 имеют 2-й и 3-й порядки малости соответственно:

Э2 Д#2

Эт 2

; 12 «-

Э2 Д£_ Эп 3Лк

, 18 «-

Э22 ДТ Эт2 6

На последнем шаге траектории, при котором пересекается сфера, в формулы (69)-(72) следует подставлять Д# = #сф.

Итоговый криволинейный интеграл вдоль оптимальной траектории

^еш! 13 3

I «2 2 1п,+2I (#сф),

;=0 п=0 п=0

(73)

где ;еш! — номер предпоследней точки траектории, лежащей еще внутри сферы, ;еш! определяется из условия (50). Второе слагаемое в (73) рассчитывается для интегралов при уменьшенном времени движения #сф от предпоследней точки г;еш! до заданной конечной геш!.

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

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

где

2(г) = (г2, Сг2) + 2сошл,

1,567 0,501 0,432 0,324 0,179 0,501 1,932 0,179 0,218 0,091 С = 0,432 0,179 2,147 0,107 0,215 0,324 0,218 0,107 2,254 0,137 0,179 0,091 0,215 0,137 2,702

— матрица формы 4-й степени аргументов;

(74)

г2 -

[(х1 х1ц) , (х2 х2ц) , ..., (х5 х5ц) ]

— вектор квадратических отклонений от центра осевой симметрии функции;

т

-м = \ 'V' -V* -V* -V* "V* —

хц _ |_х1ц, х2ц, х3ц, х4ц, х5ц J

= [8,562 8,321 8,678 8,324 8,879]т

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

Начальная и конечная точки задавались соответственно так:

гЬе!5 =[0,000 0,567 0,268 0,345 0,191 ]т,

геш! =[1,621 1,235 1,339 1,976 1,752^.

В табл. 1 приводятся относительные минимальные значения криволинейного интеграла 1мо вдоль траекторий, найденные методом пристрелки при увеличении элементарного шага Дт с 2-11 до

2-1 при количестве циклов пристрелки 48.

Количество итераций для определения времени пересечения конечной сферы, отсчитываемого от предпоследней точки траектории #сф, было выбрано равным 16.

При Дт = 2-11 получается минимальное значение 1м.а = 44 284,808, которое и принимается за абсолютный минимум (при еще меньших значениях Дт уже сказываются погрешности вычислений, которые искажают истинную картину). Относительные превышения абсолютного минимума

8!м.0 = 1м1а(1м.о - 1м.а) ‘100%.

Отношение относительных превышений минимума при уменьшении шага Дт [в единицах «количество раз»] к = [м.оС^п+О]-1 81м.о(Дтп ).

Так, например, число 2,07 в 6-й строке табл. 1 следует понимать так: при увеличении элементарного шага с 2-7 = 1/128 до 2-6 = 1/64 относительное превышение увеличивается в к = 0,352/0,170 = = 2,07 раз. Из 4-й колонки видно, что чем меньше элементарный шаг, тем точнее траектория приближается к минимальному значению интеграла, что и должно быть при правильной работе алгоритма. При увеличении элементарного шага в два раза относительные превышения увеличиваются также примерно в два раза. То есть превышение абсолютного минимума почти линейно зависит от элементарного шага, из чего следует, что порядки малости учитываемых и отбрасываемых величин во всех предшествующих выводах вполне согласованы.

Здесь автор не рассматривает неоднозначные решения задачи оптимальных траекторий методом

■ Таблица 1

Аг 1м.о % к, раз

2-11 44284,808 - -

2-10 44284,834 0,011 -

2-9 44294,887 0,034 3,09

2-8 44314,981 0,079 2,32

2-7 44355,135 0,170 2,15

2-6 44435,579 0,352 2,07

2-5 44596,732 0,716 2,03

2-4 44920,129 1,446 2,02

2-3 45554,171 2,878 1,99

2-2 46788,159 5,665 1,97

2-1 49118,901 10,930 1,93

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

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

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

В качестве тестовой была задана функция (74). Траектория задавалась полиномами 8-й степени таким образом, что координаты Х2, Х3, Х4, Х5 являются функциями координаты х1 — собственно

аргумента Х}(х1) = к;0 + к;1х1 + \2х1 +... + к^х^, где

, = 2, 3, 4, 5.

Начальная точка

гЬе!5 = [0,000; 0,567; 0,268; 0,345; 0,191],

а конечная точка

Геш! = 0,400000; 0,748993; 0,472549;

0,563411; 0,418735.

При таком задании траектории 28 коэффициентов: к22-к28, к32-к38, к42-к48, к52-к58 — свободно варьируются, а 4 коэффициента (к21, к31, к41, к51) являются зависимыми от всех предыдущих коэффициентов и рассчитываются соответствующим образом так, что траектория всегда приходит в конечную точку. Криволинейный интеграл, рассчитанный вдоль траектории по несложным, но громоздким формулам, которые здесь не приводятся, минимизируется по тем же 28 коэффициентам градиентным методом. Прежде всего, проверим соответствие решений, получаемых методом пристрелки и градиентным методом непосредственной минимизации криволинейного интеграла. Для вышеприведенных данных вначале была построена траектория методом пристрелки и рассчитано минимальное значение криволинейного интеграла: 1ма = 10 711, 179 548 291...,которое и принимается за абсолютное или «истинное».

Далее в качестве начального приближения для градиентного метода было задано четыре девятки отсчетов по координатам

{Х2(х1;), Х3(х1,), Х4(х1,), Х5(хИ)}, ; = 0, 1, ...,8.

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

х1; = 0,05;, ; = 0,1,...,8.

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

1км = 10 711,179 583 802...

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

Так, например, число 3,32 в 5-й колонке рассчитано так:

8!отн = 1к.м ~ 1м.а • 100 • 10-8(%) =

1м.а

= 3,32 •10-8(%).

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

Незначительное относительное отклонение 8Іотн = 3,32• 10-8(%) объясняется погрешностями вычислений. В табл. 3 приводятся относительные расхождения по координатам ХДхи), где / = 1, 2, 3, 4; в сечении гиперплоскостями х1і = 0,05і, где і = = 1, 2, ..., 7, квазиоптимальной траектории (найденной градиентным методом) и оптимальной (построенной методом пристрелки). Начальные и конечные точки обеих траекторий не варьируются и всегда совпадают (і = 0 и і = 8). Для большей наглядности и критичности этих оценок из всех координат Х/(х1і) вычтены постоянная и линейная составляющие, т. е. сравнивались относительные разности координат дуговых сегментов траекторий

X- (х1і) = Х/ктаг(Хіі) - Х/°Р*Х) • 100%,

/отН' и/ V \

ХІ°рї(Х1і)

где Х/^(Х1і) и Х/°рі(Хіі) — координаты дуговых сегментов квазиоптимальной и оптимальной траекторий соответственно.

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

Далее в качестве начального приближения для градиентного метода была задана прямая, соединяющая начальную и конечную точки траектории. Градиентный метод за 13 итераций от начального значения интеграла I = 10 711,803155 583... сошелся к квазиминимальному значению 1к м = = 10 711,183 207 083...

В табл. 4 отражены эти итерации.

Во 2-й колонке записано итерационное уменьшение минимизируемого интеграла I, в 3-й колонке — абсолютное отклонение от квазиминималь-ного значения 5!км, в 4-й колонке — абсолютное отклонение от абсолютного минимального значения 81м.а.

■ Таблица 3

х1і % % ЗX5aIн, %

0,05 -0,561 -0,053 0,440 0,675

0,10 -0,615 -0,081 0,422 0,674

0,15 -0,674 -0,110 0,402 0,671

0,20 -0,736 -0,142 0,380 0,668

0,25 -0,803 -0,176 0,356 0,665

0,30 -0,875 -0,214 0,329 0,660

3,50 -0,952 -0,254 0,300 0,655

■ Таблица 2

№ итерации і 2 3 4 5

ЗІотн • ^ % 1,45 2,88 3,30 3,31 3,32

Таблица 4

j I (дробная часть) 51к.м (дробная часть) 1м.а (дробная часть)

0 ,803155583 ,619949 ,623607

1 ,212249541 ,029042 ,032701

2 ,185817233 ,002610 ,006269

3 ,184263480 ,001056 ,004715

4 ,183971188 ,000764 ,004423

5 ,183875095 ,000668 ,004327

6 ,183835496 ,000628 ,004287

7 ,183773571 ,000566 ,004225

8 ,183653609 ,000447 ,004105

9 ,183631619 ,000425 ,004083

10 ,183607433 ,000400 ,004059

11 ,183408103 ,000201 ,003860

12 ,183267060 ,000060 ,003719

13 ,183207083 ,000000 ,003659

Целая десятичная часть всех интегралов равна 10 711, а целые десятичные части величин 5!к.м и 8!м.а на всех итерациях равны нулю.

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

Литература

1. Колесников Н. Е., Макаров К. А. Построение оптимальной траектории методом изохрон // Фундаментальные проблемы естествознания и техники: Сб. тр. Междунар. конгресса. СПб.: Изд-во СПбГУ, 2000. Т. 1. № 1. С.168-169.

2. Балк М. Б., Болтянский В. Г. Геометрия масс. М.: Наука, 1987. 160 с.

3. Калиткин Н. Н. Численные методы. М.: Наука, 1978. 512 с.

4. Peter F., Swaszek, John B. Thomas. Multidimensional Spherical Coordinates. Quantization // IEEE. Transaction of Information theory. July 1983. Vol. It-29. N 4. Р.570-576.

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

Возможные применения

Полученные результаты можно рекомендовать, например, для оптимизации управления систем в многомерном обобщенном пространстве состояний или параметров [8], где необходимо перевести систему из заданного начального состояния в заданное конечное с минимальными затратами, потерями и т. д. При этом системы могут быть не только технические но и, например, экономические и другие. Оптимальная траектория вышеописанным методом пристрелки получается в виде множества координат точек , = 0, 1, 2, ... ;еш!, которые разработчик или исследователь при необходимости может, например, аппроксимировать удобными для дальнейшего использования функциями. Кроме того, должна быть задана (или известна) функция потерь или затрат. Если эта функция не задана, то ее следует либо вывести аналитически, либо рассчитать численно, либо получить экспериментально. Далее найденное оптимальное движение системы в многомерном пространстве состояний следует пересчитать в программу управления (если это требуется). Однако все эти задачи уже выходят за рамки статьи и здесь не рассматриваются.

5. Субочев С. Д. Применение многомерных сферических координат для численного интегрирования и некоторых других задач // Информационно-управ-ляющие системы. 2005. № 3(16). С. 15-22.

6. Субочев С. Д. Построение оптимальных многомерных траекторий на основе физических моделей // Исследование, разработка и применение высоких технологий в промышленности: Сб. тр. 1-й Междунар. научно-практической конф. СПб.: Изд-во Политехнического университета, 2005. Т. 1. С. 85-86.

7. Яворский Б. М., Детлаф А. А. Справочник по физике. М.: Наука, 1964. 847 с.

8. Гроздовский Г. Л., Иванов Ю. А., Токарев В. В. Механика космического полета. Проблемы оптимизации. М.: Наука, 1975. 536 с.