СИНТЕЗ ИТЕРАЦИОННЫХ АЛГОРИТМОВРЕШЕНИЯ КРАЕВЫХ ЗАДАЧ И НЕЛИНЕЙНЫХ УРАВНЕНИЙ
Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений
9
УДК 629.191
В. И. МИРОНОВ, Ю. В. МИРОНОВ, Р. М. ЮСУПОВ
СИНТЕЗ ИТЕРАЦИОННЫХ АЛГОРИТМОВ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ И НЕЛИНЕЙНЫХ УРАВНЕНИЙ
Рассматривается метод синтеза итерационных алгоритмов решения нелинейных краевых задач и уравнений, основанный на использовании приближенных решений упрощенных модельных задач. Приведен пример задачи расчета импульсной программы полета космического аппарата.
Ключевые слова: итерационный алгоритм, краевая задача, нелинейные уравнения.
Введение. При рассмотрении многих задач, возникающих в различных областях науки и техники, приходится применять методы решения соответствующих краевых задач и нелинейных уравнений (в частности, при синтезе высокоэффективных алгоритмов управления подвижными объектами, а также при идентификации и оценивании их динамического состояния). Это обусловливает необходимость совершенствования существующих методов решения краевых задач и нелинейных уравнений и разработки новых.
К настоящему времени разработано большое число различных универсальных методов решения таких задач, они приведены в литературе, в частности, в работах [1—6]. Вместе с тем в различных областях знаний получены приближенные решения многих упрощенных модельных задач, имеющие ограниченное применение. Так, например, применительно к динамике полета космических аппаратов получено множество аналитических и упрощенных численных алгоритмов решения задач маневрирования в модельных гравитационных полях: однородном, линеаризованном, однородном центральном, квазиньютоновском и ньютоновском. Обычно эти алгоритмы используются для определения начального приближения при решении более сложных задач. Однако они также могут быть использованы в качестве базовых элементов при создании алгоритмов численного решения усложненных задач. Некоторые вопросы синтеза алгоритмов такого рода рассматривались в работах [7].
В настоящей работе рассматривается метод решения краевых задач и нелинейных уравнений — метод приближенного корректирующего оператора (ПКО), который позволяет использовать возможные упрощенные алгоритмы приближенного расчета в схеме численного поиска точного решения полной задачи. Такой подход расширяет конструктивный базис синтеза быстродействующих алгоритмов решения указанных краевых задач и нелинейных уравнений.
Метод приближенного корректирующего оператора. Пусть поведение объекта описывается векторным дифференциальным уравнением
x& = ϕ(x,q,t), x(t0 ) = x0 , t ∈[t0 ,T ] ,
где x — п-мерный вектор фазового состояния; q — т-мерный вектор управляющих пара-
метров. Требуется найти вектор q , который переводит управляемый объект из исходного со-
стояния x(t0 ) = x0 на требуемую траекторию движения, обеспечивающую в заданный момент времени Т достижение заданных граничных условий pT = p[x(Т )] . Будем считать, что вектор pT имеет размерность m и m ≤ n . Предполагается, что функции ϕ(x,q,t) и p[x(Т )] определены в некоторых заданных областях изменения аргументов, непрерывны и обладают необходимой степенью гладкости, так что обеспечивается единственность решения рассматриваемой задачи.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
10 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
Предположим, что известен точный нелинейный оператор связи А между вектором начального состояния объекта x0 , вектором управляющих параметров q и состоянием объекта в конечный момент времени p[x(Т )]. Оператор А задается интегрированием соответствую-
щей системы дифференциальных уравнений. Тогда требуемое значение вектора q , обеспечи-
вающее достижение заданной терминальной точки pT , будет удовлетворять уравнению
pT = А(x0 , q,T ) .
(1)
Допустим далее, что известен оператор A1 , устанавливающий приближенную связь ме-
жду векторами x0 , pT и q
pT = А1(x0 ,q,T ) .
(2)
Будем также считать, что для приближенного оператора A1 известен обратный оператор
A1−1 по вектору параметров управления q . В практической ситуации оператор A1 может
формироваться путем приближенного учета физических факторов, определяющих движение, либо формального упрощения точного оператора интегрирования A , либо комбинацией этих приемов.
Используя приближенный оператор A1 , запишем уравнение (1) в следующей эквивалентной форме:
А1(q) = pT − ∆A(q) ; ∆А(q) = А(q) − A1(q) .
(3)
Подвергнем левую и правую части равенства (3) операторному преобразованию A1−1 :
q = А1−1[pT − ∆A(q)].
Для решения этого нелинейного операторного уравнения применим метод последовательных приближений. В результате получим итерационный процесс
q(k+1) = А1−1[pT − ∆A(q(k) )] , k = 0,1, 2,...
(4)
Значение вектора q в первом приближении находим из уравнения (4) при ∆A(q(0) ) ≡ 0 ,
так что q(1) = А1−1[pT ] . Согласно (4), на каждом шаге решения задачи необходимо вычислять значение разност-
ного оператора ∆A(q(k) ) . Получим более эффективный алгоритм, в котором исключается не-
обходимость определения этого оператора. Рассмотрим несколько итераций. Согласно общей схеме расчета, величина q в первом
приближении находится из условия pT = А1(q(1) ) , так что q(1) = А1−1[pT ] .
На втором шаге значение q(2) находится из условия
(5)
pT = А1(q(2) ) + А(q(1) ) − A1(q(1) ) . Учитывая (5), получаем уравнение
2pT = А(q(2) ) + А(q(1) ) , отсюда
q(2) = А1−1[2pT − A(q(1) )] = q(2) = А1−1[pT − ∆(1) ]; ∆(1) = A(q(1) ) − pT .
На третьем шаге величина q(3) находится из условия
(6)
pT = А1(q(3) ) + А(q(2) ) − A1(q(2) ) , однако из выражения (6) следует, что
(7)
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений
11
А(q(2) ) = 2pT − А(q(1) ) . Тогда после подстановки (8) в (7) находим:
pT = А1(q(3) ) + ∆(1) + ∆(2) ; ∆(2) = A(q(2) ) − pT , откуда следует, что
q(3) = А1−1[pT − ∆(1) − ∆(2) ] . Продолжив этот процесс, приходим к следующей вычислительной схеме:
(8) (9)
∑q( k +1)
= А1−1 ⎡⎢pT
−
k
∆(i)
⎤ ⎥
;
k = 1, 2,... ,
⎣⎢ i=1 ⎦⎥
∆(i) = A(q(i) ) − pT = p(q(i) , x0 ,T ) − pT ; p(q(i) , x0 ,T ) ≡ A(q(i) ) .
Выражения (10) удобно представить в следующем виде:
(10)
∑q(k +1)
= M ⎡⎢pT
−
k
∆(i
)
⎤ ⎥
,
k = 1, 2,3,...,
⎣⎢ i=1 ⎦⎥
(11)
где M[•] ≡ A1−1[•] — приближенный корректирующий оператор, определяющий алгоритм
приближенного решения краевой задачи (1).
Условия сходимости данного метода устанавливаются на основе принципа сжимающих
отображений [1, 2]. Можно показать, что вычислительный процесс сходится со скоростью
геометрической прогрессии. Очевидно, что чем ближе приближенные операторы A1 и A1−1 к точным операторам A и A−1 , тем выше скорость сходимости.
В целом, метод ПКО отличается достаточно высокой экономичностью с вычислительной точки зрения, поскольку значение неизвестного вектора q полностью уточняется на каж-
дом итерационном шаге путем однократного интегрирования дифференциальных уравнений
движения объекта.
Рассмотренный выше метод ПКО естественным образом распространяется на задачи
решения нелинейных уравнений. Согласно этому методу, для решения уравнения
λ(q) = 0
(12)
применяется соотношение
∑q k +1
=
M
⎡ ⎢−
k
λ(qi ,T )⎥⎤ ;
⎣⎢ i=0
⎥⎦
k
= 1, 2, 3, ... ,
(13)
где M [•] — ПКО, определяющий алгоритм решения приближенного уравнения λ% (q) = 0 .
Особенности и варианты применения метода ПКО. Выбор приближенных моделей в методе ПКО может быть осуществлен различными способами с учетом специфики исходных зависимостей и условий решаемой задачи. Для этого могут применяться как формальные приемы упрощения исходных моделей, так и методы их аппроксимации. Во всех случаях необходимо стремиться к тому, чтобы приближенный алгоритм решения задачи был сравнительно простым и обеспечивалась достаточно быстрая сходимость итерационного процесса.
Важной особенностью метода ПКО является то обстоятельство, что на каждой итерации
значения pT (q) и λ(q) вычисляются один раз. Этим обеспечивается высокая экономичность
вычислений. Ниже это будет показано на примере. Значение оператора M [•] может уточ-
няться в ходе вычислительного процесса по результатам каждой итерации или через несколько итераций. В этом случае можно говорить о комбинированных вариантах использования метода ПКО.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
12 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
С общих методологических позиций очевидно, что в качестве оператора M [•] могут
рассматриваться и операторы различных известных методов численного решения нелинейных уравнений, таких как метод Ньютона и др.
При линейном представлении λ% (q) вида
λ% (q) = Aq + b
из (13) следует
qk+1 = qk − Mλ(qk ) ; M = A−1 , что по форме напоминает модифицированный (упрощенный) метод Ньютона и в частном
случае при
M
=
⎡ ⎢⎣
∂λ ∂q
⎤ −1 ⎥ ⎦q=q0
совпадает с ним.
Оператор M может уточняться в ходе итераций по некоторому правилу
M k = M k (qk ,qk−1, ..., qk−s ) ,
тогда соотношение (13) при s = 1 принимает вид
qk+1 = qk − M k (qk )λ(qk ) .
(14)
В частном случае при
Mk
=
⎡ ⎢ ⎣
∂λ ∂q
⎤ −1 ⎥ ⎦q=qk
из (14) следует алгоритм
метода Ньютона.
В целом, метод ПКО выражает достаточно общую идеологию конструирования алго-
ритмов решения краевых задач и нелинейных уравнений. При использовании более простых
корректирующих операторов M может наблюдаться замедление сходимости вычислений в
окрестности решения. В этом случае целесообразно предусмотреть специальные меры для
сокращения числа итераций. Для этого можно применить параметрическую модификацию
метода ПКО либо перейти в ходе итераций к комбинированному использованию метода ПКО
с другими методами, например, с одним из методов секущих [1].
Параметрическая модификация метода ПКО может быть представлена в виде
∑q( k +1)
= M ⎢⎡pT
− αk
k
∆(i
)
⎤ ⎥
,
k = 1, 2,3,...
⎢⎣ i=1 ⎦⎥
Рациональным выбором параметров αk можно добиться ускорения сходимости вычис-
лительного процесса. Возможные способы определения этих параметров аналогичны рас-
смотренным в [6].
Для сокращения числа итераций можно использовать известные методы ускорения схо-
димости вычислительных процессов, основанные на линейной интерполяции и использова-
нии асимптотических свойств линейно сходящихся последовательностей [1, 2, 5]. Однако эти
условия могут нарушаться на одной из стадий вычислительного процесса, что может привес-
ти к его расходимости. При комбинировании метода ПКО с другими более предпочтитель-
ным является использование методов, основанных на конечно-разностной аппроксимации
матрицы Якоби
∂λ ∂q
в методе Ньютона (14), характерной для методов секущих [1, 4]. Извест-
но, что такие методы обеспечивают сверхлинейную сходимость и имеют порядок, по крайней
мере,
1+ 2
5
≈ 1, 618 .
Линейная
интерполяция
λ(q)
в
Rm
может
привести
к
необходимости
использования ряда других методов секущих [1, 3].
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений
13
Проведенный анализ и накопленный опыт решения прикладных задач позволяют рекомендовать для комбинированного использования метода ПКО следующий способ конечно-
разностной
аппроксимации
оператора
M
=
⎡ ∂λ ⎤−1
⎢ ⎣
∂q
⎥ ⎦
на
каждом
шаге,
начиная
с
некоторой
ите-
рации, который применительно к решению уравнения (12) определяется выражением
M
=
⎡ ∂λ ⎤−1 ⎣⎢∂q ⎥⎦
≈ QΛ−1 ,
Q = [∆q1, ∆q2 , ..., ∆qm ] ; Λ = [λ1, λ2 , ..., λm ] ; ∆qi = qi+1 − qi .
Значения ∆qi и λi , необходимые для проведения вычислений, определяются по резуль-
татам первых т итераций, согласно основной процедуре метода ПКО (13). При переходе к
новой итерации проводится циклический сдвиг элементов матриц Q и Λ , при котором эле-
менты ∆q1 и λ1 исключаются, все остальные элементы сдвигаются влево на одну позицию, а
на место сдвинутых элементов ∆qm и λm ставятся элементы ∆qm+1 и λm+1 . Для ускорения
вычислений можно применить пошаговую аппроксимацию обратных матриц Λ−1 на основе метода Шульца быстрых обращений [1, 4].
В заключение отметим, что рассмотренные способы обеспечения сходимости позволяют не только ускорить процессы решения задач при выбранном приближенном корректирующем операторе, но и расширить область практического использования простых операторов.
Пример. Рассмотрим особенности применения метода ПКО на примере решения задачи расчета импульсной программы полета космического аппарата из некоторого исходного со-
стояния, определяемого значениями его фазовых координат x0 , y0 , z0 ,Vx0 ,Vy0 ,Vz0 на момент
t0 , в требуемое конечное состояние xT , yT , zT за заданное время T . Будем считать, что полет происходит в нормальном гравитационном поле Земли. C уче-
том [8] уравнения движения в абсолютной геоцентрической экваториальной системе отсчета представим в виде
x& = Vx ; y& = Vy ; z& = Vz ;
V&x = −ax; V&y = −ay; V&z = (2bc − a)z;
a = b[α00 + c(d − 1)] ; b = R0r−3 ; c = 1, 5α20 R02r−2 ; d = 5z2r−2 ;
J20 = −0, 001 082 627 ; α00 = 62 564 951 м2/c2 ;
α20 = −67 889, 273 м2/c2 ; R0 = 6371 км ; r = (x2 + y2 + z2 )1/ 2 .
В качестве ПКО воспользуемся приближенным аналитическим решением этой задачи, соответствующим движению в однородном центральном гравитационном поле. В этом случае динамика объекта описывается упрощенными уравнениями:
r& = V; V& = −ω2r ,
где ω — угловая скорость орбитального движения спутника в момент t0 ;
ω = ω(x0 ) = π0r0−3 ; r0 = x02 + y02 + z02 ; π0 = 398 600, 44 км3 / с2 .
Для такой модели задача импульсного полета имеет следующее решение:
δV
= ω rТ
− cos ωTr0 − ω−1 sin ωTV0 sin ωT
;
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
14 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
( )δV =
δVT δV
1/ 2
;
α
=
δV −1δV
,
где δV — вектор требуемого импульса скорости; δV — его модуль; α — вектор направляющих косинусов импульса δV .
Совокупность приведенных соотношений определяет значение ПКО M[•] . Решение ис-
ходной задачи производится в соответствии с методом ПКО:
k
∑q(k+1) = M [rT − ∆r(i) (T )] ; q = [δV, δV , α]T ; k = 1, 2,... i=1
На каждой итерации производится вычисление невязок ∆r(i) (T ) путем однократного
интегрирования приведенной выше полной системы дифференциальных уравнений движения объекта в нормальном гравитационном поле. При этом каждый раз изменяются начальные
условия интегрирования по значениям элементов вектора скорости V0k+1 = V0 + δVk . Приведем некоторые результаты численных расчетов. Пусть в исходном состоянии кос-
мический аппарат находится на экваториальной круговой орбите с высотой h = 900 км и имеет следующие начальные значения фазовых координат:
x0 = 7243 км; y0 = 634 км; z0 = 0 ;
Vx0 = −0, 645 км/с; Vy0 = 7,376 км/с; Vz0 = 0 . Требуется определить импульсное управление, обеспечивающее попадание космического аппарата за время T = 500 с в заданную точку xT = 5762 км; yT = 4596 км; zT = 0 . Результаты расчета импульсной программы управления по итерациям приведены в таблице. Здесь δρ(T ) — значение модуля координатного промаха относительно терминаль-
ной точки.
Номер итерации
0 1 2 3 4
δρ(T ) , км
165,446 22,415 2,806 0,348 0,043
αx
–0,4317 –0,6534 –0,6304 –0,6333 –0,6330
αy
0,9020 0,7570 0,7763 0,7739 0,7742
δVx , м/с
–515 –848 –801 –807 –806
δVy , м/с
1,075 0,982 0,986 0,986 0,986
δV , м/с
1,192 1,297 1,270 1,274 1,273
Данные таблицы свидетельствуют о хорошей сходимости рассмотренной реализации метода ПКО.
Работа выполнена при поддержке РФФИ (проект № 09-08-00259).
СПИСОК ЛИТЕРАТУРЫ
1. Вержбицкий В. М. Основы численных методов. М.: Высш. школа, 2002. 840 с. 2. Красносельский М. А. и др. Приближенное решение операторных уравнений. М.: Наука, 1969. 455 с. 3. Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими
неизвестными. М.: Мир, 1975. 560 с. 4. Островский А. М. Решение уравнений и систем уравнений. М.: ИЛ, 1963. 383 с. 5. Трауб Д. Ф., Вожьняковский Х. Общая теория оптимальных алгоритмов. М.: Мир, 1983. 382 с. 6. Федоренко Р. П. Приближенное решение задач оптимального управления. М.: Наука, 1978. 488 с.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
Навигационное обеспечение посадки воздушных судов с применением ГЛОНАСС-технологий 15
7. Миронов В. И. Конструктивный метод решения краевых задач управляемого движения // Алгоритмы и программы исследования систем управления. Вып. 6. Л.: ВИКИ им. А. Ф. Можайского, 1980. С. 70—74.
8. Космические траекторные измерения / Под ред. П. А. Агаджанова, В. Е. Дулевича, А. А. Коростелева. М.: Сов. радио, 1969. 504 с.
Вячеслав Иванович Миронов Юрий Вячеславович Миронов
Рафаэль Мидхатович Юсупов
Сведения об авторах — д-р техн. наук, профессор; Санкт-Петербургский институт информати-
ки и автоматизации РАН; E-mail: mironuv@yandex.ru — д-р техн. наук; Санкт-Петербургский институт информатики и автома-
тизации РАН; ст. научный сотрудник; E-mail: mironuv@yandex.ru — д-р техн. наук, профессор; Санкт-Петербургский институт информатики и автоматизации РАН; член-корреспондент РАН; E-mail: spiiran@iias.spb.su
Рекомендована институтом
Поступила в редакцию 28.05.09 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
9
УДК 629.191
В. И. МИРОНОВ, Ю. В. МИРОНОВ, Р. М. ЮСУПОВ
СИНТЕЗ ИТЕРАЦИОННЫХ АЛГОРИТМОВ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ И НЕЛИНЕЙНЫХ УРАВНЕНИЙ
Рассматривается метод синтеза итерационных алгоритмов решения нелинейных краевых задач и уравнений, основанный на использовании приближенных решений упрощенных модельных задач. Приведен пример задачи расчета импульсной программы полета космического аппарата.
Ключевые слова: итерационный алгоритм, краевая задача, нелинейные уравнения.
Введение. При рассмотрении многих задач, возникающих в различных областях науки и техники, приходится применять методы решения соответствующих краевых задач и нелинейных уравнений (в частности, при синтезе высокоэффективных алгоритмов управления подвижными объектами, а также при идентификации и оценивании их динамического состояния). Это обусловливает необходимость совершенствования существующих методов решения краевых задач и нелинейных уравнений и разработки новых.
К настоящему времени разработано большое число различных универсальных методов решения таких задач, они приведены в литературе, в частности, в работах [1—6]. Вместе с тем в различных областях знаний получены приближенные решения многих упрощенных модельных задач, имеющие ограниченное применение. Так, например, применительно к динамике полета космических аппаратов получено множество аналитических и упрощенных численных алгоритмов решения задач маневрирования в модельных гравитационных полях: однородном, линеаризованном, однородном центральном, квазиньютоновском и ньютоновском. Обычно эти алгоритмы используются для определения начального приближения при решении более сложных задач. Однако они также могут быть использованы в качестве базовых элементов при создании алгоритмов численного решения усложненных задач. Некоторые вопросы синтеза алгоритмов такого рода рассматривались в работах [7].
В настоящей работе рассматривается метод решения краевых задач и нелинейных уравнений — метод приближенного корректирующего оператора (ПКО), который позволяет использовать возможные упрощенные алгоритмы приближенного расчета в схеме численного поиска точного решения полной задачи. Такой подход расширяет конструктивный базис синтеза быстродействующих алгоритмов решения указанных краевых задач и нелинейных уравнений.
Метод приближенного корректирующего оператора. Пусть поведение объекта описывается векторным дифференциальным уравнением
x& = ϕ(x,q,t), x(t0 ) = x0 , t ∈[t0 ,T ] ,
где x — п-мерный вектор фазового состояния; q — т-мерный вектор управляющих пара-
метров. Требуется найти вектор q , который переводит управляемый объект из исходного со-
стояния x(t0 ) = x0 на требуемую траекторию движения, обеспечивающую в заданный момент времени Т достижение заданных граничных условий pT = p[x(Т )] . Будем считать, что вектор pT имеет размерность m и m ≤ n . Предполагается, что функции ϕ(x,q,t) и p[x(Т )] определены в некоторых заданных областях изменения аргументов, непрерывны и обладают необходимой степенью гладкости, так что обеспечивается единственность решения рассматриваемой задачи.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
10 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
Предположим, что известен точный нелинейный оператор связи А между вектором начального состояния объекта x0 , вектором управляющих параметров q и состоянием объекта в конечный момент времени p[x(Т )]. Оператор А задается интегрированием соответствую-
щей системы дифференциальных уравнений. Тогда требуемое значение вектора q , обеспечи-
вающее достижение заданной терминальной точки pT , будет удовлетворять уравнению
pT = А(x0 , q,T ) .
(1)
Допустим далее, что известен оператор A1 , устанавливающий приближенную связь ме-
жду векторами x0 , pT и q
pT = А1(x0 ,q,T ) .
(2)
Будем также считать, что для приближенного оператора A1 известен обратный оператор
A1−1 по вектору параметров управления q . В практической ситуации оператор A1 может
формироваться путем приближенного учета физических факторов, определяющих движение, либо формального упрощения точного оператора интегрирования A , либо комбинацией этих приемов.
Используя приближенный оператор A1 , запишем уравнение (1) в следующей эквивалентной форме:
А1(q) = pT − ∆A(q) ; ∆А(q) = А(q) − A1(q) .
(3)
Подвергнем левую и правую части равенства (3) операторному преобразованию A1−1 :
q = А1−1[pT − ∆A(q)].
Для решения этого нелинейного операторного уравнения применим метод последовательных приближений. В результате получим итерационный процесс
q(k+1) = А1−1[pT − ∆A(q(k) )] , k = 0,1, 2,...
(4)
Значение вектора q в первом приближении находим из уравнения (4) при ∆A(q(0) ) ≡ 0 ,
так что q(1) = А1−1[pT ] . Согласно (4), на каждом шаге решения задачи необходимо вычислять значение разност-
ного оператора ∆A(q(k) ) . Получим более эффективный алгоритм, в котором исключается не-
обходимость определения этого оператора. Рассмотрим несколько итераций. Согласно общей схеме расчета, величина q в первом
приближении находится из условия pT = А1(q(1) ) , так что q(1) = А1−1[pT ] .
На втором шаге значение q(2) находится из условия
(5)
pT = А1(q(2) ) + А(q(1) ) − A1(q(1) ) . Учитывая (5), получаем уравнение
2pT = А(q(2) ) + А(q(1) ) , отсюда
q(2) = А1−1[2pT − A(q(1) )] = q(2) = А1−1[pT − ∆(1) ]; ∆(1) = A(q(1) ) − pT .
На третьем шаге величина q(3) находится из условия
(6)
pT = А1(q(3) ) + А(q(2) ) − A1(q(2) ) , однако из выражения (6) следует, что
(7)
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений
11
А(q(2) ) = 2pT − А(q(1) ) . Тогда после подстановки (8) в (7) находим:
pT = А1(q(3) ) + ∆(1) + ∆(2) ; ∆(2) = A(q(2) ) − pT , откуда следует, что
q(3) = А1−1[pT − ∆(1) − ∆(2) ] . Продолжив этот процесс, приходим к следующей вычислительной схеме:
(8) (9)
∑q( k +1)
= А1−1 ⎡⎢pT
−
k
∆(i)
⎤ ⎥
;
k = 1, 2,... ,
⎣⎢ i=1 ⎦⎥
∆(i) = A(q(i) ) − pT = p(q(i) , x0 ,T ) − pT ; p(q(i) , x0 ,T ) ≡ A(q(i) ) .
Выражения (10) удобно представить в следующем виде:
(10)
∑q(k +1)
= M ⎡⎢pT
−
k
∆(i
)
⎤ ⎥
,
k = 1, 2,3,...,
⎣⎢ i=1 ⎦⎥
(11)
где M[•] ≡ A1−1[•] — приближенный корректирующий оператор, определяющий алгоритм
приближенного решения краевой задачи (1).
Условия сходимости данного метода устанавливаются на основе принципа сжимающих
отображений [1, 2]. Можно показать, что вычислительный процесс сходится со скоростью
геометрической прогрессии. Очевидно, что чем ближе приближенные операторы A1 и A1−1 к точным операторам A и A−1 , тем выше скорость сходимости.
В целом, метод ПКО отличается достаточно высокой экономичностью с вычислительной точки зрения, поскольку значение неизвестного вектора q полностью уточняется на каж-
дом итерационном шаге путем однократного интегрирования дифференциальных уравнений
движения объекта.
Рассмотренный выше метод ПКО естественным образом распространяется на задачи
решения нелинейных уравнений. Согласно этому методу, для решения уравнения
λ(q) = 0
(12)
применяется соотношение
∑q k +1
=
M
⎡ ⎢−
k
λ(qi ,T )⎥⎤ ;
⎣⎢ i=0
⎥⎦
k
= 1, 2, 3, ... ,
(13)
где M [•] — ПКО, определяющий алгоритм решения приближенного уравнения λ% (q) = 0 .
Особенности и варианты применения метода ПКО. Выбор приближенных моделей в методе ПКО может быть осуществлен различными способами с учетом специфики исходных зависимостей и условий решаемой задачи. Для этого могут применяться как формальные приемы упрощения исходных моделей, так и методы их аппроксимации. Во всех случаях необходимо стремиться к тому, чтобы приближенный алгоритм решения задачи был сравнительно простым и обеспечивалась достаточно быстрая сходимость итерационного процесса.
Важной особенностью метода ПКО является то обстоятельство, что на каждой итерации
значения pT (q) и λ(q) вычисляются один раз. Этим обеспечивается высокая экономичность
вычислений. Ниже это будет показано на примере. Значение оператора M [•] может уточ-
няться в ходе вычислительного процесса по результатам каждой итерации или через несколько итераций. В этом случае можно говорить о комбинированных вариантах использования метода ПКО.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
12 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
С общих методологических позиций очевидно, что в качестве оператора M [•] могут
рассматриваться и операторы различных известных методов численного решения нелинейных уравнений, таких как метод Ньютона и др.
При линейном представлении λ% (q) вида
λ% (q) = Aq + b
из (13) следует
qk+1 = qk − Mλ(qk ) ; M = A−1 , что по форме напоминает модифицированный (упрощенный) метод Ньютона и в частном
случае при
M
=
⎡ ⎢⎣
∂λ ∂q
⎤ −1 ⎥ ⎦q=q0
совпадает с ним.
Оператор M может уточняться в ходе итераций по некоторому правилу
M k = M k (qk ,qk−1, ..., qk−s ) ,
тогда соотношение (13) при s = 1 принимает вид
qk+1 = qk − M k (qk )λ(qk ) .
(14)
В частном случае при
Mk
=
⎡ ⎢ ⎣
∂λ ∂q
⎤ −1 ⎥ ⎦q=qk
из (14) следует алгоритм
метода Ньютона.
В целом, метод ПКО выражает достаточно общую идеологию конструирования алго-
ритмов решения краевых задач и нелинейных уравнений. При использовании более простых
корректирующих операторов M может наблюдаться замедление сходимости вычислений в
окрестности решения. В этом случае целесообразно предусмотреть специальные меры для
сокращения числа итераций. Для этого можно применить параметрическую модификацию
метода ПКО либо перейти в ходе итераций к комбинированному использованию метода ПКО
с другими методами, например, с одним из методов секущих [1].
Параметрическая модификация метода ПКО может быть представлена в виде
∑q( k +1)
= M ⎢⎡pT
− αk
k
∆(i
)
⎤ ⎥
,
k = 1, 2,3,...
⎢⎣ i=1 ⎦⎥
Рациональным выбором параметров αk можно добиться ускорения сходимости вычис-
лительного процесса. Возможные способы определения этих параметров аналогичны рас-
смотренным в [6].
Для сокращения числа итераций можно использовать известные методы ускорения схо-
димости вычислительных процессов, основанные на линейной интерполяции и использова-
нии асимптотических свойств линейно сходящихся последовательностей [1, 2, 5]. Однако эти
условия могут нарушаться на одной из стадий вычислительного процесса, что может привес-
ти к его расходимости. При комбинировании метода ПКО с другими более предпочтитель-
ным является использование методов, основанных на конечно-разностной аппроксимации
матрицы Якоби
∂λ ∂q
в методе Ньютона (14), характерной для методов секущих [1, 4]. Извест-
но, что такие методы обеспечивают сверхлинейную сходимость и имеют порядок, по крайней
мере,
1+ 2
5
≈ 1, 618 .
Линейная
интерполяция
λ(q)
в
Rm
может
привести
к
необходимости
использования ряда других методов секущих [1, 3].
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений
13
Проведенный анализ и накопленный опыт решения прикладных задач позволяют рекомендовать для комбинированного использования метода ПКО следующий способ конечно-
разностной
аппроксимации
оператора
M
=
⎡ ∂λ ⎤−1
⎢ ⎣
∂q
⎥ ⎦
на
каждом
шаге,
начиная
с
некоторой
ите-
рации, который применительно к решению уравнения (12) определяется выражением
M
=
⎡ ∂λ ⎤−1 ⎣⎢∂q ⎥⎦
≈ QΛ−1 ,
Q = [∆q1, ∆q2 , ..., ∆qm ] ; Λ = [λ1, λ2 , ..., λm ] ; ∆qi = qi+1 − qi .
Значения ∆qi и λi , необходимые для проведения вычислений, определяются по резуль-
татам первых т итераций, согласно основной процедуре метода ПКО (13). При переходе к
новой итерации проводится циклический сдвиг элементов матриц Q и Λ , при котором эле-
менты ∆q1 и λ1 исключаются, все остальные элементы сдвигаются влево на одну позицию, а
на место сдвинутых элементов ∆qm и λm ставятся элементы ∆qm+1 и λm+1 . Для ускорения
вычислений можно применить пошаговую аппроксимацию обратных матриц Λ−1 на основе метода Шульца быстрых обращений [1, 4].
В заключение отметим, что рассмотренные способы обеспечения сходимости позволяют не только ускорить процессы решения задач при выбранном приближенном корректирующем операторе, но и расширить область практического использования простых операторов.
Пример. Рассмотрим особенности применения метода ПКО на примере решения задачи расчета импульсной программы полета космического аппарата из некоторого исходного со-
стояния, определяемого значениями его фазовых координат x0 , y0 , z0 ,Vx0 ,Vy0 ,Vz0 на момент
t0 , в требуемое конечное состояние xT , yT , zT за заданное время T . Будем считать, что полет происходит в нормальном гравитационном поле Земли. C уче-
том [8] уравнения движения в абсолютной геоцентрической экваториальной системе отсчета представим в виде
x& = Vx ; y& = Vy ; z& = Vz ;
V&x = −ax; V&y = −ay; V&z = (2bc − a)z;
a = b[α00 + c(d − 1)] ; b = R0r−3 ; c = 1, 5α20 R02r−2 ; d = 5z2r−2 ;
J20 = −0, 001 082 627 ; α00 = 62 564 951 м2/c2 ;
α20 = −67 889, 273 м2/c2 ; R0 = 6371 км ; r = (x2 + y2 + z2 )1/ 2 .
В качестве ПКО воспользуемся приближенным аналитическим решением этой задачи, соответствующим движению в однородном центральном гравитационном поле. В этом случае динамика объекта описывается упрощенными уравнениями:
r& = V; V& = −ω2r ,
где ω — угловая скорость орбитального движения спутника в момент t0 ;
ω = ω(x0 ) = π0r0−3 ; r0 = x02 + y02 + z02 ; π0 = 398 600, 44 км3 / с2 .
Для такой модели задача импульсного полета имеет следующее решение:
δV
= ω rТ
− cos ωTr0 − ω−1 sin ωTV0 sin ωT
;
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
14 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
( )δV =
δVT δV
1/ 2
;
α
=
δV −1δV
,
где δV — вектор требуемого импульса скорости; δV — его модуль; α — вектор направляющих косинусов импульса δV .
Совокупность приведенных соотношений определяет значение ПКО M[•] . Решение ис-
ходной задачи производится в соответствии с методом ПКО:
k
∑q(k+1) = M [rT − ∆r(i) (T )] ; q = [δV, δV , α]T ; k = 1, 2,... i=1
На каждой итерации производится вычисление невязок ∆r(i) (T ) путем однократного
интегрирования приведенной выше полной системы дифференциальных уравнений движения объекта в нормальном гравитационном поле. При этом каждый раз изменяются начальные
условия интегрирования по значениям элементов вектора скорости V0k+1 = V0 + δVk . Приведем некоторые результаты численных расчетов. Пусть в исходном состоянии кос-
мический аппарат находится на экваториальной круговой орбите с высотой h = 900 км и имеет следующие начальные значения фазовых координат:
x0 = 7243 км; y0 = 634 км; z0 = 0 ;
Vx0 = −0, 645 км/с; Vy0 = 7,376 км/с; Vz0 = 0 . Требуется определить импульсное управление, обеспечивающее попадание космического аппарата за время T = 500 с в заданную точку xT = 5762 км; yT = 4596 км; zT = 0 . Результаты расчета импульсной программы управления по итерациям приведены в таблице. Здесь δρ(T ) — значение модуля координатного промаха относительно терминаль-
ной точки.
Номер итерации
0 1 2 3 4
δρ(T ) , км
165,446 22,415 2,806 0,348 0,043
αx
–0,4317 –0,6534 –0,6304 –0,6333 –0,6330
αy
0,9020 0,7570 0,7763 0,7739 0,7742
δVx , м/с
–515 –848 –801 –807 –806
δVy , м/с
1,075 0,982 0,986 0,986 0,986
δV , м/с
1,192 1,297 1,270 1,274 1,273
Данные таблицы свидетельствуют о хорошей сходимости рассмотренной реализации метода ПКО.
Работа выполнена при поддержке РФФИ (проект № 09-08-00259).
СПИСОК ЛИТЕРАТУРЫ
1. Вержбицкий В. М. Основы численных методов. М.: Высш. школа, 2002. 840 с. 2. Красносельский М. А. и др. Приближенное решение операторных уравнений. М.: Наука, 1969. 455 с. 3. Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими
неизвестными. М.: Мир, 1975. 560 с. 4. Островский А. М. Решение уравнений и систем уравнений. М.: ИЛ, 1963. 383 с. 5. Трауб Д. Ф., Вожьняковский Х. Общая теория оптимальных алгоритмов. М.: Мир, 1983. 382 с. 6. Федоренко Р. П. Приближенное решение задач оптимального управления. М.: Наука, 1978. 488 с.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1
Навигационное обеспечение посадки воздушных судов с применением ГЛОНАСС-технологий 15
7. Миронов В. И. Конструктивный метод решения краевых задач управляемого движения // Алгоритмы и программы исследования систем управления. Вып. 6. Л.: ВИКИ им. А. Ф. Можайского, 1980. С. 70—74.
8. Космические траекторные измерения / Под ред. П. А. Агаджанова, В. Е. Дулевича, А. А. Коростелева. М.: Сов. радио, 1969. 504 с.
Вячеслав Иванович Миронов Юрий Вячеславович Миронов
Рафаэль Мидхатович Юсупов
Сведения об авторах — д-р техн. наук, профессор; Санкт-Петербургский институт информати-
ки и автоматизации РАН; E-mail: mironuv@yandex.ru — д-р техн. наук; Санкт-Петербургский институт информатики и автома-
тизации РАН; ст. научный сотрудник; E-mail: mironuv@yandex.ru — д-р техн. наук, профессор; Санкт-Петербургский институт информатики и автоматизации РАН; член-корреспондент РАН; E-mail: spiiran@iias.spb.su
Рекомендована институтом
Поступила в редакцию 28.05.09 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 1