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

кандидата физико-математических наук
Коннов, Александр Иванович
город
Москва
год
1999
специальность ВАК РФ
05.13.01
Диссертация по информатике, вычислительной технике и управлению на тему «Итеративные методы последовательного улучшения управления динамическими системами»

Автореферат диссертации по теме "Итеративные методы последовательного улучшения управления динамическими системами"

РОССИЙСКАЯ АКАДЕМИЯ НАУК

ИНСТИТУТ ПРОБЛЕМ УПРАВЛЕНИЯ *

им. В.А.Трапезникова ^ ^

РГ6 ОД ж

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

КОННОВ Александр Иванович

ИТЕРАТИВНЫЕ МЕТОДЫ ПОСЛЕДОВАТЕЛЬНОГО УЛУЧШЕНИЯ УПРАВЛЕНИЯ ДИНАМИЧЕСКИМИ СИСТЕМАМИ

Специальность: 05.13.01 — Управление в

технических системах

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

Москва - 1999

Работа выполнена в Институте проблем управления им. В.А.Трапезникова Российской Академии наук

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

профессор В.Ф.Кротов

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

а£4.С;6> 2000г. в /V час. на заседании диссертационного совета Д002.68.03 Института .проблем управления имени академика В.А.Трапезникова по адресу:-117806, Москва, Профсоюзная, д.65.

С диссертацией можно ознакомиться в библиотеке Института проблем управления РАН.

Автореферат разослан

Ученый секретарь диссертационного совета

кандидат технических наук С. А.Власов

Н.А.Бобылев

кандидат физико-математических наук К.Г.Григорьев-

Ведущая организация: Институт программных систем

Российской Академии наук.

Защита диссертации состоится

/>-/73, Ч с ^

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

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

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

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

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

Связь с планом. Исследования по теме дисертационной работы проводились в соответствии с плановой тематикой работ Института проблем управлени РАН в рамках тем, по которым аытор выступал ответ-ственнным исполнителем.

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

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

Публикации. По теме диссертации автором подготовлены две рабо-

ты. ;

Структура и содержание работы. Диссертация состоит из трех глав и списка литературы, включающего 45 наименований. Объем диссертации составляет 82 страницы текста, подготовленного в формате ЬаТеХ.'

СОДЕРЖАНИЕ РАБОТЫ. Глава 1 посвящена постановке задачи и состоянию вопроса. В ней рассмотрены наиболее распространеные ■ алгоритмы решения задач оптимального управления,

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

В случае непрерывних задач рассмотрим управляемую систему, состояние которой описывается системой уравнений вида

^ = • (1)'

с условиями

х(0) = X ¡С € Я" и<=иу

где функции /'(¿,х, и), г = 1,...,п определены при всех 1,х,и и дважды дифференцируемы по аргументам ( и I, х £ — заданный вектор, V С Ег — замкнутое множество, Еп — п-мерное евклидово пространство со скалярным. произведением (х,у). Компоненты ж(<) называются фазовыми координатами управляемой системы, а Еп — ее фазовым пространством. При этом траектория х(<) является кусочно-дифференцируемой функцией, а программа управления и(£) — кусочно-непрерывной. Множество допустимых процессов и) = (х(1),и(*)) с траекторией х(<) и программой управления и(1) обозначается через И. Программа управления такая что порожденный ею процесс XV = (х(*),и(<)) будет допустимым, называется допустимой программой управления.

На множестве Р допустимых процессов определен функционал цели

■ . т

/(«(<), «(<)) = I /°(<,х,«)^+^(х(Т))->шт (2)

о

где вектор-функция /°(<,х,и) и функция ^(х) определены при всех <, х, и и дважды дифференцируемы по аргументам I и х.

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

Рассматриваются дза подхода к решению этой задачи.

Во-первых, можно попытаться непосредственно отыскать такой элемент ui £ D, на котором минимизируемый функционал достигает своего минимального значения.

Во-вторых, можно попытаться отыскать такую последовательность элементов {ш,} € D, Кто I(w,) —> info I(uj) при s —* сю. Такая последовательность существует всегда, когда существует нижняя грань функционала I(w). Для корректного применения такого подхода, однако, одной только сходимости значений функционалов недостаточно. Для того, чтобы полученный результат можно было бы использовать в качестве конкретной рекомендации при решении задачи, часто требуется, чтобы последовательность элементов {у,} сама сходилась бы к какому-либо предельному элементу. Именно этот предельный элемент и будет рассматриваться в качестве приближенного решения задачи.

Поскольку такой элемент существует при достаточно общих препо-ложениях и поскольку само нахождение предельной.траектории еще не дает решения задачи, можно отказаться от требования существования •предельной траектории и ограничиться только сходимостью последовательности процессов. Итак, необходимо находить такую последовательность процессов, чт-;

/(u)J + 1) < I{w,) lim /(wt) = inf /(w)

1-4 OO IVgV

В случае дискретных задач оптимального управления'уравнения движения имеют вид

*(* + i) = /M0;v(0) « = о.....т-1, (з)

с условиями

х(0) = х х£Бп и £ U.

Координаты системы х € Еп, управление системы и € U С Ет. Время t S Л, где А = {0,1,2,... , Т} — подмножество целых чисел. Функции /' (t, х, и), i = 1,..., п определены при всех допустимых t, х, и и дважды дифференцируемы по х, х € Еп —заданный вектор, U С Ег — замкнутое множество, Е" — n-мерное евклидово пространство. Множество допустимых процессов w = (x(t),u(t)) с траекторией x(f) и программой . управления u(t) также обозначим через D.

Функционал цели для этой задачи имееет вид • Т-1

I(x(t),u(t))= £ f°(t,x,u)+F(x(T))-+min (4)

t=о

где вектор-функция /°(4,х,и) и функция Р(х) определены при всех 1, х, и и дважды дифференцируемы по х.

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

В 1.2 рассмотрены условия оптимальности и уравнения оптимального процесса.

При этом при описании методов, используемых при решении задач оптимального управления, будем следовать подходу, описанному в работах В.Ф.Кротова. Впедем необходимые обозначения.

Введем основные конструкции.

В случае непрерывной задачи рассмотрим класс П вещественных дифференцируемых функций : [0;У] х Еп -> II (функция Кро-

това) и соответствующие ей конструкции

Здесь вектор-функция /(¿,х,и) — правая часть уравнений движения (1), — подынтегральная часть минимизируемого функцио-

нала (2), а Р(х) — его терминальная часть. Через х) и обозначаются частные производные функции х) по I и по а: соответственно, ж), /(¿, х, и)) — евклидово произведение векторов х) и /(£,х,и) в пространстве Еп.

Для построенного таким образом функционала Ь(ги, <р) справедливо следующее утверждение: функционалы £(и/, ¡р) и 1(ги) тождественно равны для всякого допустимого процесса ю £ Б и функции х) £ П.

В случае дискретной задачи класс II функций представляет

собой отображения <р(1, х) : А х Еп —> К, где А — множество точек {0,1,2,..Т}. Функции Н(1,х,и), (7(х) и Ь{ы) при этом имеют вид

ад = Р{х) + <р{т,х)-<р{о,х)

(6)

т

о

С(х) = Р{х) + 1р(Т,х)-1р{ 0,х)

(8) (9)

■ Т-1

(Ю)

о

При этом вектор-функция f(t,x,u) представляет собой правую часть системы (3). Функционал L также на всяком допустимом процессе тождественно равен I.

Введем функции u(t,x), P(t,x) и функционал H(t,xl>,x,«), определяемые следующим равенством:

R(t,x,u(t,x)) = P(t,x) = maxR(t,x,u) (11)

и€£/

Н(г,ф,х,и) = (ф, f(t,x,u))~ f°{t,x,u) (12)

Через J обозначим, матрицу первых производных правой части уравнений движения

JtAt,x,u)- gxJ .

Иногда для краткости будем обозначать J(t, ro(i)i «о(0) через J(t), где (xo{i), uq(<)) € D — некоторый изначально заданный допустимый процесс.

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

Теорема 1 Пусть имеется функция *p{t,x) £ П и последовательность {ад,} С D, такая, что

1.

i i

J R{t,x,(t),u,(t))dt j ¡i{t)dt,

2.

т > -оо,

где = вир^кЕ" R(t¡x,гt), а т = ¡пГ 1££»С(1). Тогда эта последовательность минимизирует функцнонал 1.

Доказательство этого утверждения имеется в литературе.

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

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

В следующем параграфе приводится пример рассуждений для вывода уравнения Беллмана.

Однако само решение этого уравнения с использованием разностных схем требует колоссального объема памяти ЭВМ и растет эскпо-ненциально с увеличением размерности задачи — как где N —

число узлов разбиения по каждому измерению, а п — размерность задачи.

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

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

2. Методы, основанные на последовательном улучшении программы управления. Эти методы являются наиболее разработанными, но проблемы здесь связаны с учетом фазовых ограничений.

3'. Методы, основанные на последовательном приближении процессов к процессу, удовлетворяющему критерию оптимальности. Обычно при этом от предельного процесса требуется, чтобы он удовлетворял принципу максимума. При этом сами улучшаемые процессы могут не удовлетворять уравнениям движения.

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

1. интегрирование уравнений движения (1) с начальными условиями х(0) = х и программой управления ио(<) для получения соответствующей траектории а:о(<).

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

Наиболее широко распространенным является градиентный метод. Рассмотрим этот метод, используя прежние обозначения, использующие функцию ip(t, х). Будем искать новый процесс w в задаче улучшения, который бы был достаточно близок к улучшаемому процессу to, так что знак разности ДI = I(w) — I(wq) совпадал с линейной частью разности

Т

A I = AL = Gx(x{T))Ax{T) - J (Я, Ах + Я„Ди) dt

о

Уравнение для AL определяется видом функции ip(t,x). Потребуем, чтобы функция y>(ttx) удовлетворяла условиям

Rx(t,x0(t),u0{t)) =0, G,(x(T)) = 0

Эти уравнения содержат только первые производные функции <p(t, х) вдоль траектории xo(t). Выписав их, получим следующую систему уравнений для ф(1) = <px(t,x):

ф{1) + JT{t, r0(<), tio(<))#') ~ /?(*. *o(0. «о(«)) =0

ф(Т) + Fx(x(T)) = 0 ^

Для дискретных задач эти уравнения имеют вид

ф(1) - + 1) - /°(<, го(0,«о(«)) = 0

Ф{Т) + №(Т)) = 0

(14)

Это означает, что такие уравнения будут выполнены, если мы возьмем функцию <р{1,х) = ф(1)х, где вектор-функция ф{1,х) находится из решения системы (13) или (14). Будем называть такую функцию функцией локального улучшения программы управления. Затем

т

61{ги0)=5Цю0,<р) = ! Я«{1,хъ(*),«о(1))ЛиА (15)

о

при этом Ли(',£о(<). «о(*)) = Я„(«, ф(1),х0(1),-и0(1)).

Положим, что имеется заданная функция й"ио(*) и произвольный малый параметр е > 0, такой что

1. Правая часть равенства (15) положительна.

2. u{i;e) = uo(i) + e6u0(t) € U, t € [0;T]

3. Траектория.Го(<;е), определяемая программой управления u(i;e), уравнениями движения и начальными условиями, является допустимой, то есть w(e) = (x(t;e),u(t;e)) 6 D.

Тогда найдется е > 0, такое что

/(«>)< /(u'o), ,и>=ш(е)

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

1. Находим траекторию xo{t) из решения задачи Коши с заданной программой управления ио(<)

2. Находим вектор-функцию ф(1), решая линейную задачу Коши (13) и затем ftu(i,xo(t),^o(i)) . Отсюда определяем функцию tp(t,x) = y>(t)x

3. Строим вариацию программы управления öuo(t), при которой правая часть равенства (15) положительна.

4. При различных положительных е находим траекторию процесса xo(i;i), решая задачу Коши с программой управления и(<) =

+ eitio(i). Величина е подбирается из тех соображений, чтобы величина функционала I(w) уменьшалась.

Этот алгоритм относится к алгоритмам "челночного" типа.

Выражение (15) задает градиент функционала I(w) в пространстве программ управлений u(i). В силу этого, предложенный метод может рассматриваться как применение градиентного подхода к решению задачи.

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

Метод, рассматриваемый в данной работе, основан на методе глобального улучшения.

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

Опишем способ улучшения для непрерывных процессов. В случае дискретных процессов общий подход сохраняется, меняется только вид некоторых уравнений. Способ этот состоит в следующем: пусть задана допустимая программа управления «о(<), тогда

0. Из уравнения процесса (1) и начальных условий с программой и = ио(0 определяем траекторию хо(<) и процесс хид = (хо(*)> ио(')) € О

1. Строим функцию <р(1,х) £ П, такую, что

Д(4,хо(*).«о(0) = ттЩг,х,и0(г))

СЫТ)) = тахС(х) <1Й>

2. Определяем управление и({,х)'иэ условия (11) и из уравнений движения находим новый процесс и;.

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

Теорема 2. Для процессов ги и и>0, построенных описанным выше образом, справедливо неравенство 1(ю) < /(т^о). Если существует г £ [0,Т], такое что

Е(г,хо{т),щ{т))ф Р{т,х0{т)),

то 1{ги) < /(и/о).

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

Для реализации этой схемы улучшающая функция у>(1,х) задается в виде

<р(г, х) = <*,(<), х) + ¿(Дх, £(<)Дх> (17)

где Ах = х - хо(<), а- вектор-функцня ф(1) и матричная функция Е(4) подлежат определению. Здесь £(<) — матрица вторых производных

функции tf>\l,x). Первые необходимые условия выполнения неравенств (16)

dRjt.xc (t),u„(fl) = 0

дх

ООЫТ)) _ дх~ '

приводят «ас к сопряженной системе уравнений для ip(t).

Обозначим через A(t) матрицу квадратичной формы d2R((, ®o(i)i tio(t)), а через M(t) — матрицу квадратичной формы d2H(t, xo(t), i на процессе w0. Таким образом, Лу(<) = а F

à'mt,x o(t).Uo(t))

йг'бх'

Вторые необходимые условия выполнения уравнений (16) приводят к следующей системе дифференциальных неравенств:

." £*3Я>0 d*R = {Ax,A{t)Ax)

Непосредственным дифференцированием находим

Л(1) = ¿(í) + £(t) J(t) + + M(t)

Ъхдс вввх ^ ^К1 I

Для дискретных процессов гти уравнения имеют вид A(i) = -Щ) + fHWt + 1 )J(t) + M(i)

(18)

(19)

(20)

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

A«(t) = ¿¿(<) i,j = 1.....n

(21)

где 6(1) ~ {<?,(<)} — положительная вектор-функция.

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

элементами на главной диагонали, то есть

3« 'fix' ~~ T J ''J — 1' ••■>"

o(T)) _ i — 1 n * '

¿а'дя' — » — •!(•••»"

где a = {ffi} — положительный вектор. После этого уравнения (19) или (20), определяющие E(f), сводятся к решению линейного дифференциального уравнения с заданными значениями на правом конце.

' Таким образом, алгоритм выглядит следующим образом:

0. задается произвольная допустимая программа управления uo(i)

1. используя уравнения процесса, замкнутые программой и = щ(<), определяем траекторию xo(i). Для процесса wq =-(хо(0>uo(t)) ■вычисляется значение функционала /(«>о)-

2. используя уравнения сопряженную систему уравнений (13), определяем вектор-функция xp(t)

3. устанавливаем все функции J,- (<) i — 1,..., г» равными неотрицательной константе. Аналогично, устанавливаем равными неотрицательной константе а,-. Используя уравнения (19), определяем матрицу E(i) и функцию <p(t,x) в соответствии с (17)

4. для этой функции <p(t, г) находим управление й(i, г) в соответствии с уравнением (11)

5. находим траекторию х(<), используя уравнения движения с управлением u(t) — хI(i, x(i)) и получаем новый процесс«; = (x(t),u(t)).

После этого сравниваем значения функционалов 1(хи) и I(xv0) Если I(w) > I(wо), то увеличиваем значения 6 и а, после чего повторяем операции- 3-5. Способ подбора и улучшения векторов S и а основан в данном случае на интуитивных соображениях и на результатах • численных экспериментов.

Замечание 1. Для улучшения достаточно выполнения неравенств

J R{t,x{t),u0{tj)dt > j R{t,x0(t),u0{t))dt

G(x(T)) < G(x0(T)).

Последнее всегда справедливо, если справедливо (16). Но выполнение п.п. 1., 2., 3. алгоритма гарантирует только относительный минимум

функции R(t,x,u0(t)) в точке xq(<) при каждом t € (0;Т) и максимум функции G(x) в точке xq(T). Предполагается, что зазор между относительным минимумом (максимумом) этих функций и данным неравенством устраняется, при необходимости, указанным подбором настроечных параметров. Это предположение подтверждается опытом применения данного алгоритма.

Тем не менее, остается теоретический пробел в обосновании этого подхода: отсутствие теоремы существования решения уравнений (16) в виде (17).

Замечание 2. Для нахождения функции £(i) необходимо решать систему линейных дифференциальных уравнений размерности л(п+1)/2. То есть, размерность этой системы возрастает пропорционально квадрату размерности п. При большом п.затраты на решение этой системы могут стать неоправданно велики.. Это тем более неоправданно, что (21) и (22) не являются необходимыми условиями (16). Поэтому возможна и актуальна разработка более рациональных схем реализации (16), (17). Решение этих проблем и составляет основное содержание диссертации.

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

Для решения задач с нефиксированным временем

/ = F(T,x(T)) min .

вводится независимая переменная г 6 [0; 1] и управление р, а уравнения движения заменяются уравнениями

dt

Тт = ' •jjf = />/'(<>*.") ¿ = 1,...,"

после чего минимизируется функционал

/' = F(t(l),x(l)) min

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

дк(Т,х(Т)) = 0, к = 1.....д

В этом случае минимизируется функционал

1 = Р + ^2(Хкдк^Вд1) (23)

аг=1

Здесь А* — постоянные, к= 1,..В — заданная положительная константа, являющаяся параметром метода. После того, как улучшение достигнуто, величинам Аь присваиваются новые значения

\ь-*>1к+Вдк к= 1, —, д

Затем для новых значений А* функционал (23) снова минимизируется.

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

В главе 2 рассматривается описание алгоритма нового метода, его обоснование, возможные модификации и рекомендации по подбору па-раметрова численного метода. При этом функция функции <р^,х) ищется в виде у>(1,х) = я) + |(£(<)Дг, Дг).

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

1. Функция /°(*,г,и) дважды дифференцируема по г и найдутся такие вещественные К, М, что при достаточно больших по норме х, функцию /°(<, х, и) можно ограничить квадратичпой функцией7 а именно: х,и) < К\\х\\2 при \\х\\ >М

2. Функция /(¿, х, и) удовлетворяет следующему условию: найдутся такие вещественные К, М, что при достаточно больших по норме х, функцию /(*, х, и) можно ограничить линейной функцией, а именно: при всех < и и из области определения, /(<,£,«) < /<||а:|| при ||ж|| > М. При всех' ({, х, и) £ [0;Т] х Я" х V матрица =

ограничена по норме.

3. Функция Р{х) удовлетворяет следующему условию: найдутся та- ■ кие вещественные К,М, что при достаточно больших по норме х, функцию Р(я) можно ограничить квадратичной функцией, а именно: Р(х) < К||х||2 при ||х|| > М

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

■ Л(<,хо(<),ио(<)) = ттЯ(<,Е,«о(<))

С(х0(Т)) = т1хС(х) Е

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

Теорема 3. Если для задачи (1)-(2) (случай непрерывной задачи) выполнены дополнительные условия 1.-3., то для любого допустимого процесса хио'= (хо(<)> «о(*)) € О найдутся таки$ числа А, В и С, что, для Л1&боН^у.цкц1!и 1р{1, х) = (ф(1), Ах) + 1(т{1)(Ах, Ах) будут выполнены условия (24). При этом вектор-функция т/>(<) находится из решения системы (13), а функция сг(1) : [0; Т\'-Ь Я удовлетворяет неравенствам

сг(Т) + 2А < 0 • ¿¿■(i) - |or(i)|fl -i-C > 0

Из этой теоремы можно вывести следующее следствие:

Следствие. Функция E(t) может быть задана в вице Е(<) = (a(e"Y(T-t)_ 1) + Р)Е, где Е — единичная матрица, числа а,/? < 0, у > 0.

Теорема 4. Если для задачи (3)-(4) (случай дискретной задачи) выполнены дополнительные условия 1.-3., то для любого допустимого процесса щ = (zo(i)>uo(î)) 6 D найдутся такие числа А, В и С, что для любой функции <p(i, х) = (t/i(<), Ах) + |cr(i){Ar, Дг) будут выполнены условия (24). При этом вектор-функция 4>{t) находится из решения системы (13), а функция cr(t) : {0,1, . ■■ ,Т] R удовлетворяет неравенствам

v"

<r(T) + 2Л < 0 -i,r{t) + ar(t + 1)5 + С>0 ¿ = 0,1,...,Г-1

Доказательство этих утверждений имеется в приложении. В соответствии с этими теоремами предлагается искать функцию £(<) в виде £(£) = (Q(e7(T-t) — 1) + Р)Е, где а, /3 и у — вещественные числа; а,/3 < 0, у > О, Е — единичная матрица. В этом случае алгоритм поиска улучшения будет следующим:

0. задается произвольная допустимая программа управления u0^t)

1. используя уравнения процесса, с и = uo(t) определяется траектория xo(t). Для процесса wq = (xo(i), uo(0) вычисляется значение функционала f(wo)-

2. используя уравнения (13), определяется вектор-функция ф{1)

3. устанавливаем 7 > 0, а,Р < 0. Определяем функцию y(i,x) в виде <p(t, х) = (Щ, х) + Ах), <r(t) (а(е»<т-0- 1) +/?).

4. для этой функции <p(t, х) находим управление u(t, х) в соответствии с уравнением (11)

5. находим траекторию x(t), используя уравнения движения с упра-' влением u(t) = u(t,x(t)), находим новый процесс w и сравниваем

значения функционалов.

Если функционал не улучшился., то увеличиваем значение у и уменьшаем а и /?, после чего повторяем Операции 3-5. Согласно следствию faKne а, в, 7 существуют. При этом, вообще говоря, при подборе параметров Е{£) достаточно изменять только один из параметров.

В 2.2 рассматривается модификация рассматриваемого алгоритма, основанная на анализе уравнений, задающих условия для функции Е(<).

Рассмотрим его на примере решения непрерывной задачи'(1)-(2). В этом случае существует матрица перехода В, определяемая из условий

B(t) = -JT(t,zb(t),u0(i))B(t) (25)'.

В(0) = Е, Е — единичная матрица. [

Поскольку норма матрицы j|J|| ограничена, то матрица В существует и определена на всем отрезке [0;Т].

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

¿(0 + Е(0 J{t) + /T(i)ST(<) =; А(<) - M(t),

где A(t) — произвольная положительно определенная матрйца, имеет Ьид E(i) = B(t)S(t)BT (t), где матрица S(t) определяется из решения Уравнения

S(t) = B-l{t)(A(t)-M(t))(BT)-\t).

Если в качестве Л(<) возьмем такую положительно определенную матрицу, что A(t) — M(t) = kB(t)BT (i), к > 0, то уравнение для S(t) будет иметь вид 5(i) = кЕ, откуда получаем, что S{t) = 5(0) + ktE. Отсюда следует, что функция E(i) определяется из уравнения

Е(<) = B(<)5(0)ßT(<) + ktB{t)BT(t).

Значение S(0) определяется из граничных условий (19). Для того, чтобы матрица второй квадратичной формы G(x(T)) была отрицательно определена, необходимо, чтобы отрицательно определена была матрица

^ (F(x(T)) + B{T)S(0)B(T) + кТВ(Т)Вт(Т)), откуда получаем условие для 5(0): «

5(0) = В-'(Т) (с- ^l-F(r(T))) (В^-Г))"1 - кТЕ,

где С — произвольная отрицательно определенная матрица. Возьмем в качестве С такую отрицательно определенную, матрицу, что С- gf^F(x(T)) = 1В(Т)Вт(Т), КО. В этом случае 5(0) = (1-кТ)Е. Так получаем еще одно уравнение для определения функции Е[t).

E(t) = B(t)(l + k(t-T))BT(t) (26)

Подобный выбор функции £(i) позволяет избавиться от достаточно грубой оценки, когда получаются слишком большие величины £(<). Правда, за это надо платить вычислением матрицы B(t). В то же время такой выбор матрицы Е(<) более экономичен в вычислительном плане по сравнению с исходным методом, поскольку не требуется вычислять значение M(t) =

При вычислениях может оказаться удобнее вычислять матрицу В(<) при граййчных условиях В(Т) = Е (при этом ip(t) = В(*)^(Г)). В этом случае, повторяя все вышеприведенные рассуждения, получим, что функция £(<) будет иметь тот же самый вид

Е (t) = B(t)(l + Kt-T))BT(t)

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

вычислении функции V>(i) на шаге 2. алгоритма использовать не старую программу управления uo(i), а программу управления, выбранную из условия (11) максимизации R(t,x,v). То есть, при вычислении вектор-функции V'(i) решать систему уравнений

7 * =

| щ(1) = argrnaxH(t^(t),x0{t),u) У'1'

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

L(xo,UQ]tp) > L(xo,vii\<p) — из условия построения îli(t,x). Здесь й\(t) — программа управления, полученная на шаге 2. алгоритма из условия (27). В самом деле, поскольку L{x,u;<p) = G(z(T)) — /0Г R(t,x, u) dt, .то поскольку R(t, xo(t), uo(t)) < R(t, xo(t),îIi(t)), то и

L(xo, щ; <p) > Ь(х,щ]<р) — из условия построения <p[t,x), удовлетворяющей условиям (16).

L(x,u\\<p) > L(x,u;ip) — поскольку v. выбирается из условия (11). Отсюда следует, что L(xoiuo;y) > L(x,u-,tp). Следует заметить при этом, что процесс (xo(t),ui(t)) не будет допустимым, поскольку будут нарушены уравнения движения (1).

В 2.4 рассматриваются примеры использования описанных.выше алгоритмов на примере решения модельных задач.

В Примере 1. рассматривается следующая задача оптимального управления: на функции х(() и u(i) .наложены условия:

£ = и M < 1 £(0) = х(Т) = 0, требуется минимизировать функционал

/ = / (и2 - г2) dt -+ min.

J о ■

Решение этой задачи известно — при Т < ж имеется единственная экстремаль — пара функций x(i) = 0, u(t) = 0. Значение функционала при этом равно 0. При Т > 1г к ней добавляются экстремали, являющиеся сопряжением линейных функций и синусоид:

{±t при t < Т\

±к cos (< - Г/2) при n<t< г2 ±Т T t при t > 7*2

«(О = |

±1 при I < Т\

- Т/2) при п < < < т2 при < > Г2

Значение функционала при этом равно I = —51"?. Здесь п, тг, Л выбираются из условий гладкости, то есть х = ±1 при * = Г1,их = т1 при 1 = т3.

Кроме того, при Т = п?г, где п — целое,'добавляются также синусоиды, то есть

х(£) = ±вш{ и(/) = ±С08<

При этом значение функционала равно 0.

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

В Примере 2 рассматривалась задача о брахистохроне в следующей постановке: требуется минимизировать функционал

Ш1П

с условиями х = и; х(0) = хо", х(Т) = Х).

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

Т = 2;. х0 = 3; х, = 10

х(1) = {

В качестве начального приближения был выбран процесс с траекторией

3 при 0<<<1

3 + 7(< — 1) при 1 < < < 2

Пройедено сравнение решений этой задачи, полученных градиентным методом и методом, описанным в данной работе.

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

Поведение динамической системы описывается следующей системой уравнений:

т = —и

h = v

v = -g+ — (Vu - Ce-1,AVsgnv) m \ 1

iG[0;T] T = 70; «e[0;0,04] m(0) = 0; h(0) = 0; «(0) = 0

Здесь m — масса ракеты, h — вертикальная координата (высота), v — вертикальная скорость; V — постоянная, характеризующая величину реактивной тяги, д,С, 7 — заданные постоянные, связанные с силой тяготения, аэродинамическим сопротивлением и убыванием плотности атмосферы с высотой.

Задача решалась при значениях

д = 0,01 К = 2,0 С = 0,05 7 = 0,01.

Управление u(i) определяет режим расхода горючего.

Необходимо обеспечить максимальную высоту подъема за заданное время

Л (Г) ->• шах

при дополнительном ограничении на расход топлива т(Т) = 0,2. (Максимальная высота d заданное время при заданном расходе горючего).

Качественное решение &той задачи известно и выглядит оно следующим образом:

{0,04 при 0 < < < <i u'{t) при ¿1 < t < <2 0 при t2 < t < Т

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

Условие тп(Т) — 0, 2 обеспечивал использованием модифицированных множителей Лагранжа. .

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

Описанные в данной работе алгоритмы' обеспечивали достаточное улучшение на первых же итерациях, в зависимости от модификации метода — от 2 до 10. Кроме того, они были менее чувствительны к выбору начального процесса.

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

В примерах 4 и 5 рассмотрено применение описанного здесь алгоритма к решению дискретных задач.

В Примере 4 требуется минимизировать функционал

/ = — (6ui +4ti2 + "з) min при следующих условиях

(«1 + 2и2 + Зи3 < 5 2щ + иг + из < 4

и- {0,1}

Решение этой задачи известно: = иг = 1, иг = 0, 7 = —10.

В Примере 5 рассматривается сходная с предыдущей задача.

Алгоритм решения этих задач следующий: Задача заменяется эквивалентной ей многошаговой задачей оптимального управления. Приводится пример построения функции <p(t, х) при котором не происходит улучшения начального процесса и пример, при котором оптимальное управление достигаемся за одну итерацию.

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

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

Покажем, что увеличение абсолютных значений собственных чисел функции £(£) приведет к тому, что новая траектория будет "прижиматься" к исходной, В самом деле, рассмотрим выражение ^{Ах, Ах) — 2(Д/, Ах) для оценки "стремления" новой траектории отклониться

от исходной. Поскольку управление tl(i,x) выбирается нз условия (V* 4- £Ах,А/) —> шах то при достаточно больших Ах (таких, что ||т/>|| < ||£Ах||) программа управления будет в значительной степени определяться условием максимизации (ЕАх, А/). Поскольку E(i) отрицательно определена в каждый момент времен», то отсюда следует, что (Ах, А/) —> min. То есть, по мере отклонения от начальной траектории, на новую траекторию начинает действовать новый фактор, стремящийся приблизить ее к исходной траектории.

Эти соображения позволяют оценить и длину шага интегрирования на отрезке [0;Т]. При достаточно больших по модулю собственных числах £(<) процессы, получаемые с помощью этого алгоритма могут иметь вырожденные участки ("псевдоскользящие режимы"), объясняемые не структурой задачи, а параметрами численного метода.

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

Расчеты приведенных ниже задач производились на персональной ЭВМ IBM PC с процессором Pentiurn-166 и 32MB оперативной памяти. Для ввода параметров и вывода результатов использовался пакет MathCad.

Как правило, для расчетов использовалась равномерная сетка разбиения отрезка интегрирования с числом узлов от 1000 до 10000, что давало возможность контроля точности метода. Для приведенных ниже задач достаточным оказалось йспользования 1000 - 2000 узлов. В качестве метода интегрирования использовался метод Рунге-Кутта четвертого порядка.

В задаче о мягкой посадке рассатривалась следующая задача:

г

и ■ = -+--^

т -г г

Pcosd ( V2 ц

т

Р sin(? tii»

m г Р

(28)

т —--

с

Т -4 min

с фазовыми переменными г, и, v и m и управлением Р и в. Ра [0; Ртах],

= и

V =

6 е [0; 2тг}. Эта задача соответствует плоскому движению реактивного летательного аппарата в гравитационном поле. Р — величина тяги, в — угол, который составляет вектор тяги с радиус-вектором, направленным от центра планеты, ц = 7М — произведение гравитационной Постоянной на массу планеты, с — характеристика двигателя (скорость истечения газа), г — расстояние до центра планеты, и — нормальная составляющая скорости, v — ее тангенциальная составляющая, т — приведенная безразмерная масса.

r(0) -RL + h u(0) = 0 t)(0) = v0 m(0) = 1 .

r {T)~Rl u(0) = 0 v(0) = 0

Здесь Ri — радиус планеты, h — высота орбиты, vq —\/h/(Rl + h) — круговая скорость. Эта задача соответствует задаче мягкой посадки с круговой орбиты на поверхность планеты. В качестве планеты была выбрана Луна, Исходные данные ц — А,924 • 1012 м3/с2, Ртах = 4 м/с2, с = 3100м/с, Я/. = 1,737 • 10® м, h = 163 • 103м. Требуется минимизировать время спускания аппарата Т min

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

Задача решалась тремя способами — градиентным методом, методом, описанным в главе 1 (Вариант 1) н способом, предложенным в данной работе (Вариант 2).

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

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

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

В задаче об управлении летательным аппаратом в атмосфере рассматривалась следующая задача.' В качестве фазовых координат выбираются высота — у и проекции на плоскость горизонта —' х и z

соответственно, V — скорость аппарата, углы & и ф задают направление вектора скорости в полярных координатах. Зависимостью силой тяжести от высоты пренебрегаем, зависимость плотности от высоты считаем линейной. В безразмерных координатах уравнения движения выглядят следующим образом:

а = a(t) и /? = (3(t) — управляющие функции (углы атаки и скольжения) с ограничениями |а| < |ат1Х| и Щ < |/5юах|. Заданные функции имеют вид

На параметр наложено ограничение 0 < Ц\ < 1.

В расчетах принимались следующие начальные значения: го = 1, 5; 1?о = 0; ^о = 0; х0 = -350; уо = 10; го = 0 н следующие значения параметров: к = 9,2; атах = /?тах = 0,3; Л = 6,4; <т — 10; ро = 3,2; рц = 0,5; г = 1,14; а = Ь = 1; сх = 1,5; с0 = 0,035.

Решалась задача быстродействия. Требовалось за минимальное время достичь

Как и в предыдущем случе, метод, при котором необходимо вычислять матрицу £(<) из решения дифференциального уравнения, показывает более хорошие результаты. Но само решение этого уравнения требует больших вычислительных затрат.

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

i = kv eos $ eos ф y = kvsine z = kv eos в sin ф

cy — aa cx — bfí cx = c0 + ci(c* + cj)

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

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

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

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

4. Создан пакет программ для решения задач динамического программирован ия.

Основные результаты диссертации представлены к публикации в следующих работах

1. А.И.Коннов, В.Ф.Кротов О глобальных итеративных алгоритмах оптимизации процессов управления //Труды Института, М.: Институт проблем управления им. Трапезникова РАН 1999г., т. V, с. 44

2. А.И.Коннов, В.Ф.Кротов О глобальных методах последовательного улушения управляемых процессов // Автоматика и телемеханика 1999г., N10, с. 77

Зак.125.ТарД00.15ПУ.

Библиография Коннов, Александр Иванович, диссертация по теме Системный анализ, управление и обработка информации (по отраслям)

1. A.E. Bryson, Y.Ch.Ho Applied Optimal Control, Hemisphere Publishing Corporation, Washington, DC.

2. H.J.Kelley Gradient theory of optimal flight paths. MRS J.30(10) 1962.

3. H.J.Kelley, R.E.Kopp, H.G.Mayer Successive approximation techniques for trajectory optimization. Proceeding of Symphosium on Vehicle System Optimization. New York, 1961.

4. V.A.Kazakov, J.Somloi, D.J.Tannor Optimal Controlof quantum systems: I. A shortsighted optimization procedure, Chem. Phys., 172, 1983.

5. V.A.Kazakov, J.Somloi, Yu. Orlov Control of photochemical branching: new procedures for finding optimal pulses and global upper bounds, in Time-Dependent Quantum Molecular dynamics (J.Broeckhove and L.Lathouwers, eds.) PlenumPress, New York, 1992.

6. Vadim F.Krotov Global Methods in Optimal Control Theory. Marcel Dekker, Inc., New-York, 1996.

7. V.F.Krotov A technique of global bounds of optimal control theory, Control Cybern 17 (2-3) 1988.

8. H. M. Wagner Principles of Operations Research, Prentice Hall, Englewood Cliffs, N.J., 1969.

9. Алексеев В.M., Тихомиров В.M., Фомин C.B. Оптимальное управление. Москва, Наука 1979г.

10. Беллман Р. Динамическое программирование. Москва, Изд-во иностранной лит. 1960г.

11. Беллман Р., Дрейфус С. Прикладные задачи динамического программирования. Москва, Наука 1965г.

12. Бахвалов Н.С. Численные методы. Москва, Наука 1975г.

13. Болтянский В. Г. Математические методы оптимального управления. Москва, Наука 1969г.

14. Бутковский А.Г. Методы управления системами с распределенными параметрами. Москва, Наука 1975г.

15. Габасов Р., Кириллова Ф.М. Конструктивные методы оптимального управления. ИАН СССР Техническая кибернетика, N 2, 1983г.

16. Габасов Р. Кириллова Ф.М. Принципы максимума в теории оптимального управления. Минск, Наука и техника, 1974г.

17. Годунов С.К., Рябенький B.C. Разностные схемы. Москва, Наука 1977г.

18. Гурман В.И. Принцип расширения в задачах управления. Москва, Наука 1985г.

19. Гурман В.И., Батурин В.А., Данилина Е.В. Новый метод улучшения управляемых процессов, Наука, Новосибирск, 1983г.

20. Гурман В.И., Расина И.В. Практическое применение достаточных условий сильного относительного минимума. Autom. Remote Control, N10, 1979г.

21. Гурман В.И., Батурин В.А., Москаленко А.И. и др. Методы улучшения в вычислительном эксперименте. Новосибирск, Наука, 1988г.

22. Гурман В.И., Батурин В.А., Данилина Е.В. и др. Новые методы улучшения управляемых процессов. Новосибирск, Наука, 1987г.

23. В.А.Емеличев В.И.Комлик Метод построения последовательности планов для дискретных задач оптимизации, Наука, Москва, 1981г.

24. Еремин И.И., Астафьев H.H. Введение в теорию линейного и выпуклого программирования. Москва, Наука 1976 г.

25. Иослович И.О., Борщевский М.З. Некоторые задачи оптимизации стабилизации осесимметричного спутника. Космические исследования, вып.З 1966г.

26. Иоффе А.Д., Тихомиров В.М. Теория экстремальных задач. Москва, Наука, 1974г.

27. Казаков В. А., Кротов В.Ф. Оптимальное управление резонансным взаимодействием света и вещества. Автоматика и телемеханика, N4, 1983г.

28. Колмогоров А.Н., Фомин С.В. Элементы теории функций и функционального анализа. Москва, Наука 1976г.

29. Кротов В.Ф. Методы решения вариационных задач на основе достаточных условий абсолютного минимума. Ч 1-4. Автоматика и телемеханика, N12 т.23 1962г., N15 т.24 1963г., N7 т.25 1964г., N1 т.26 1965г.

30. Кротов В.Ф. Приближенный синтез оптимального управления. Автоматика и телемеханика N11 т25. 1964г.

31. Кротов В.Ф. Методы оптимизации процессов управления с минимаксным критерием. Ч 1-3. Автоматика и телемеханика N12 1973г., N1, N3 1974г.

32. Кротов В.Ф. Вычислительные алгоритмы решения и оптимизации управляемых систем уравнений. Ч I, И. ИАН СССР Техническая кибернетика, N5, N6 1975г.

33. Кротов В.Ф., Вукреев В.З., Гурман В.И. Новые методы вариационного исчисления в динамике полета. Москва, Машиностроение 1969г.

34. Кротов В.Ф., Гурман В. И. Методы и задачи оптимального управления Москва, Наука 1973г.

35. Кротов В.Ф., Фельдман И.Н. Итерационный метод решения задач оптимального управления. ИАН СССР Техническая кибернетика, N 2, 1983г.

36. ЛюбушинА.А., Черпоусъко Ф.Л. Метод последовательных приближений для расчета оптимального управления. ИАН СССР Техническая кибернетика, N 2, 1983г.

37. Марчук Г.И. Методы вычислительной математики Москва, Наука, 1977г.

38. Моисеев H.H. Численные методы в теории оптимальных систем. Москва, Наука 1971г.

39. Потрягин Л. С, Болтянский В.Г., Гамкрелидзе Р.В., Мищенко Е.В. Математическая теория оптимальных процессов. Москва, Физмат-гиз, 1976г.

40. Пропой А.И. Элементы теории оптимальных дискретных систем. Москва, Наука 1981г.

41. Стечкин С.Б., Субботин Ю.Н. Сплайны в вычислительной математике. Москва, Наука 1976г.

42. Табак Д., Куо B.C. Оптимальное управление и математическое программирование. Москва, Наука, 1975г.

43. Федоренко Р.П. Приближенное решение задач оптимального управления. Москва, Наука, 1978г.

44. Филиппов А.Ф. О некоторых вопросах теории оптимального регулирования. Вестник МГУ, N 2, 1959г.

45. Энеев Т.М. О применении градиентного метода в задачах теории оптимального управления. Космические исследования, т4. N 5, 1966г.