ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И ПРОЦЕССЫ УПРАВЛЕНИЯ
N 1, 2004 Электронный журнал, рег. N П23275 от 07.03.97
http: //www. neva. ru/journal e-mail: diff@osipenko. stu.neva. ru
Теория обыкновенных дифференциальных уравнений
МЕТОД КУСОЧНО-ЛИНЕЙНОЙ АППРОКСИМАЦИИ ДЛЯ РЕШЕНИЯ ЗАДАЧ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ
В.И.БОЛДЫРЕВ
Россия, 630090, Новосибирск, пр. ак. Коптюга, NN 4, Институт математики им. С.Л.Соболева СО РАН,
Кибернетика, e-mail: vibold@math.nsc.ru
Аннотация
В работе получены следующие основные результаты:
1. Разработан симплексный метод для решения задачи минимизации псевдовыпуклой функции на выпуклом компактном множестве. Доказана сходимость получаемой алгоритмом минимизирующей последовательности.
2. Разработан метод последовательных приближений для задачи минимизации псевдовыпуклого функционала, определенного на множестве конечных состояний линейной системы управления. Дано обоснование предлагаемого алгоритма.
3. Дан способ построения практически "реализуемых"управлений для решения линейных задач с линейными терминальными ограничениями и без них.
4. Разработан алгоритм для решения краевых нелинейных задач управления.
Введение
Математическое моделирование многих динамических процессов, возникающих на практике (промышленное производство, экономика, экология, химия, биология, движение летательных аппаратов и т. д.) является в настоящее время основным инструментом получения знаний о их поведении при различных способах воздействия. Одна из главных целей моделирования — поиск такого управляющего воздействия, при котором достигается в некотором смысле "максимальный эффект". Например, минимальные затраты ресурса (времени) на производство единицы продукции или перевод управляемого объекта из начального состояния в заданное конечное.
Наиболее удобным и распространенным средством описания динамических процессов являются дифференциальные уравнения. Возникающие при этом задачи, как правило, хорошо известны в теории оптимального управления. Однако, подавляющее большинство из них не имеют простого (аналитического) решения и требуют разработки численных методов.
Актуальной проблеме разработки эффективных методов численного решения линейных и нелинейных задач оптимального управления посвящена настоящая работа.
Линейные задачи оптимального управления привлекают внимание исследователей по двум основным причинам:
— возможностью простого описания и получения конструктивных эффективных методов (алгоритмов) решения;
— возможностью использования линеаризации нелинейных задач и, следовательно, управлением нелинейными системами с помощью линеаризованных моделей.
Как правило, для линейных задач удается доказать локальную и глобальную сходимость предлагаемых методов решения, исследовать также скорость сходимости и получить ее оценки. Мы ограничимся рассмотрением детерминированных задач, в которых управляемая система описывается обыкновенными дифференциальными уравнениями без запаздываний в фазовых координатах и управлении. Исследование свойств решений любых задач преследует конечную и, по-видимому, главную цель — найти решение задачи в явном виде, но, если этого сделать невозможно, предложить численный метод его поиска. Методы решения опираются на теоретические факты и, естественно, используют теоремы о виде оптимального управ-
ления, о конечности числа переключений, о существовании оптимального управления, о единственности и телесности областей достижимости линейных систем. Следует подчеркнуть, что телесность области достижимости зависит только от области управления и параметров системы, а не от вида функционала задачи, фазовых и терминальных ограничений.
Известно, что принцип максимума сводит задачу оптимального управления к решению двухточечной краевой задачи для системы дифференциальных уравнений. Характерным для задач оптимального управления является то, что аналитическое решение задачи удается получить лишь в редких случаях. В связи с этим большую роль играют численные методы построения оптимального управления. Потребности практики и бурный прогресс вычислительной техники стимулировали разработку вычислительных методов оптимального управления.
Работы в области создания численных методов оптимального управления интенсивно ведутся с середины прошлого века, с момента появления принципа максимума Л. С. Понтрягина [68] и динамического программирования Р. Беллмана [10]. Трудности при решении задач оптимального управления вызваны необходимостью решать краевую задачу, большой размерностью систем, наличием ограничений на управление и фазовые координаты, многоэкстремальностью. Это привело к большому разнообразию вычислительных методов их преодоления.
Методы решения задач оптимального управления можно условно разбить на:
1) градиентные методы в пространстве управлений, основанные на формуле приращения для первой вариации функционала (к ним можно отнести некоторые методы из [18], [61], [75], [92], [95]);
2) методы, связанные с варьированием и перебором траекторий в пространстве фазовых состояний (см. [7], [60], [83], [97], [116]);
3) методы решения двухточечной краевой задачи, полученной из принципа максимума (к ним можно отнести некоторые методы, описание которых можно найти в работах [18], [20], [39], [62], [78]);
4) методы, использующие дискретизацию задачи с последующим применением линейного и нелинейного программирования (см., например, [38],
[84]);
5) методы последовательных приближений, основанные на процедурах линеаризации и вариации управлений (см., например, [1], [5], [21]-[25], [26],
[45], [46], [48], [51], [52], [53], [56], [58], [59], [73], [74], [76], [80], [81], [84], [85]);
6) методы штрафных функционалов или модифицированного функционала Лагранжа (см., например, [23], [34], [38], [58], [67], [84], [105], [112]).
Первые работы по численным методам решения задач оптимального управления сделаны Д. Е. Охоцимским и Т. М. Энеевым [65], Л. И. Шатров-ским [87], А. Брайсоном и В. Денхемом [17]. В последствии было опубликовано значительное количество статей и книг на эту тему. Можно выделить несколько различных направлений в разработке численных методов, существенно отличающихся друг от друга. Прежде всего, следует упомянуть прямые методы, основанные на спуске в пространстве управлений. Ряд исследований связан с непрямыми методами, в которых с помощью принципа максимума Л. С. Понтрягина исходная задача редуцировалась к краевой. Цикл работ по численным методам оптимального управления был выполнен Н. Н. Моисеевым [61] и его учениками И. А. Крыловым, Ф. Л. Черноусько, И. А. Вателем и др. В этих работах был разработан подход, основанный на вариациях в пространстве состояний.
Другое направление, развиваемое Р. П. Федоренко [84], базируется на использовании идей метода линеаризации. Методы условного градиента и проекции градиента были перенесены В. Ф. Демьяновым и А. М. Рубиновым [37] на задачи оптимального управления. Указанный Н. Н. Моисеевым [61] в 1971 году подход, основанный на методах нелинейного программирования, нашел широкое развитие в работах Э. Полака [66], Ю. Н. Ермольева, В. П. Гуленко, Т. И. Царенко [40], А. И. Пропоя [69], Д. Табака и Б. Куо [77], Ю. Г. Евтушенко [38] и др. Обзор различных численных методов оптимального управления содержится в [18], [37], [52], [58], [62], [84], [85], [86].
Рассмотрим методы численного решения задачи оптимального управления, не касаясь редукции ее к краевой задаче и методов решения краевой задачи.
Пусть динамический процесс описывается системой нелинейных обыкновенных дифференциальных уравнений
х = /(х,и,1), I Е Т =[¿0,^1], ж(£о) = х0. (1)
Здесь х — п-мерный вектор фазовых координат; и — й-мерный вектор управления, принадлежащий классу кусочно-непрерывных функций, и(Ь) Е и, где и С — выпуклое компактное множество.
Требуется найти допустимое управление и(Ь), минимизирующее функ-
ционал
Fo(u) = ^o(x(ti)) (2)
при дополнительных условиях
^(м) = ^г(х(^)) = 0 (< 0), г = 1,т (т < п). (3)
Предполагается, что вектор-функция ](х,м,£) непрерывна по £ и непрерывно дифференцируема по х и м; функции ^¿(х), г = 0,т, — непрерывно дифференцируемы по х.
Для каждого функционала ^¿(м), г = 0,т, введем гамильтониан
Н(^,х,м,г) = (х, м, ¿)) и сопряженную систему
^* = -Д(х,М)У, = (4)
Тогда градиенты функционалов F(u), i = 0,m, вычисляются по формулам
F(u) = -K(^,x,u,t).
Метод последовательных приближений для решения задач оптимального управления со свободным правым концом (m = 0) использует принцип максимума Л. С. Понтрягина. Основу метода составляют различные типы линеаризации динамической системы вместе с соответствующими процедурами варьирования управлений [5], [26], [52], [56], [58], [73], [74], [76], [84].
Пусть для допустимого управления u(t) вычислена траектория x(t) и решение сопряженной системы (4) для i = 0. Введем функцию
AvH[t] = AvH(^(t),x(t),u(t),t) = H(^(t),x(t),v,t) -H(^(t),x(t),u(t),t)
и пусть функция U = u(t), t E T, такая, что u(t) в каждый момент t Е T является решением следующей задачи
Au{t)H[t] = max AvH[t], t E T. (5)
veU
Пусть для заданного управления uk (t) вычислена траектория xk (t), решение сопряженной системы (4) для i = 0, управление uk(t) из условия максимума (5). Тогда последовательные приближения в соответствии с
методом [52] можно записать так: мк+1(г) = ик(г), к = 1,.... Но, к сожалению, для этого простого варианта метода нельзя гарантировать сходимость. Различные модификации этого метода, позволяющие улучшить его сходимость, изложены в [58].
Построим функцию
Айкн[г] = ,хк,ик,г) -п(фк,хк,ик,г).
Найдем Тк = тик из условия максимума
Лик^н[Тк] = т^х^к (¿)н[г]. (6)
Последовательные приближения в соответствии с методом [26] строятся по формулам
или
где
\ик(г), г е т - т£к, т£к = {г е т : Аик(4)Н[г] > £кАи^н^]}
Тк = [тк - ^к(тк - го), Тк + £к(г1 - Тк)],
£к е [0,1]
отыскивается из условия убывания функционала ^0(м).
Из условия максимума (6) следует, что если АикН[тк] =0, то допустимое управление ик (г) удовлетворяет принципу максимума. Методы типа [52],[58], также как и итерационный процесс (7), нуждаются в дополнении, связанном с тем, что принцип максимума может вырождаться. Поэтому в [26] итерационный процесс (7) дополняется вычислительной процедурой, основанной на необходимом условии оптимальности особого управления (необходимом условии оптимальности второго порядка). Описанию модификаций метода улучшения допустимого управления второго порядка для задач оптимального управления посвящена работа [8]. Такие методы являются достаточно громоздскими, поэтому в вычислительной практике они не получили широкого применения. Более предпочтительным в этом отношении является метод [5], полученный на основе линеаризации управляемой системы без общепринятой коррекции линейного члена в разложении правой части.
Наиболее эффективным в рамках принципа максимума является метод улучшения управления [73]. Этот метод дает следующую процедуру варьирования управления и(£):
ф I <д; (8)
где
V= arg шах[^*/(ж(£),г>,£)], 9=
А^/[¿] = /(ж(*),и(*),*) - /(ж(*),и(*),*).
Параметр д > 0 подбирается так, чтобы обеспечить уменьшение функционала ^(и). Используемая в данном методе линеаризация однозначно ориентирована на получение принципа максимума. При такой линеаризации формируется и выводится в остаток слагаемое Аи}/ж[£]Аж (для возмущенного управления = и(£) + Аи(£)). Такая операция снижает качество аппроксимации Аи}^ (и) и сказывается на эффективности соответствующих методов. Это обстоятельство учтено в методе [5]. Процедура варьирования в этом методе выглядит следующим образом:
и(^,р) = 1 9(^) <Р (9)
[ > Р.
Для каждого р > 0 находится решение сопряженной системы
^ = -/г(ж(*),и(*,^,р),*)*^, ^(¿1) = -^0Х(ж(^х)). (10)
В результате формируются управление ир(£) = и(£,^р(£), р) и параметр р > 0 подбирается так, чтобы обеспечить уменьшение функционала ^(ир) < ^ (и).
Разница между двумя приведенными способами варьирования вполне наглядна. В (7) сопряженная траектория зафиксирована относитель-
но управления и(£) и варьирование производится на базе заранее известной функции переключения 9(£, и)). В (9), (10) происходит совместное варьирование управления и сопряженной вектор-функции. Управление ир(£) вырабатывается на основе функции переключения 9(£,^р(¿)) с возмущенным решением = сопряженной системы. Процедура варьирования (9), (10) для каждого р требует, в отличие от (7), дополнительного интегрирования сопряженной системы.
Численное решение задач оптимального управления с терминальными ограничениями существенно усложняется по сравнению с решением задач со свободным правым концом. Это связано с необходимостью учитывать дополнительные ограничения, причем итерационные методы зависят от типа ограничений. В задачах с ограничениями-неравенствами [74] методы последовательных приближений конструируются в классе допустимых управлений с сохранением всех функциональных ограничений на каждой итерации. В задачах с ограничениями-равенствами методы работают в классе доступных управлений, когда функциональные ограничения не выполнены.
Рассмотрим вначале задачу оптимального управления (1)-(3) при ограничениях типа неравенства. Развиваемый В. А. Срочко в [74] подход связан с линеаризацией функционалов задачи в рамках процедур игольчатого варьирования. Пусть и(£), £ Е Т — допустимое управление, х(£) = х(£,и) — соответствующая ему фазовая траектория. Выделим набор индексов активных ограничений вместе с индексом целевого функционала
Как известно, для решения задачи об улучшении управления и(£) используется принцип максимума. В традиционной формулировке принципа максимума обычно присутствуют множители Лагранжа, вопрос об отыскании которых остается открытым. Поэтому в монографии В. А. Срочко [74] другая форма представления необходимых условий (принципа максимума, дифференциального принципа максимума) более приемлема для реализации. Вводится характеристика допустимого управления
Здесь V — множество допустимых управлений. Очевидно, что 6(и) > 0, и управление и(£) удовлетворяет принципу максимума тогда и только тогда, когда 6(и) = 0. Аналогичным образом в случае дифференциального
I = {0} и {г = 1, т : ^(и) = 0}. Будем использовать обозначения
Д [£, и] = Н(фг(£), х(£), V(£),£) - Н(фг(£), х(£), и(£), £), КМ= К №(£),х(£),и(£),£).
т
принципа максимума вводится величина:
n(u) = maxmin / (HU[t,u], v(t) — u(t))dt.
veV iel J T
Тогда равенство n(u) = 0 эквивалентно выполнению дифференциального принципа максимума для управления u(t).
Для варьирования управления u(t) выбирается некоторое управление v Е V и вводится функция варьирования x(t), t Е T, из пространства ), принимающая только два значения {0,1}. Функция варьирования нормируется с помощью параметра а Е [0,1]: fT x(t)dt = a(ti — to). Обозначим через X^; семейство функций варьирования. Совокупность варьированных управлений строится в виде
Uv,x(t) = u(t) + x(t)(v(t) — u(t)), t Е T.
Для того, чтобы при варьировании управления u Е V обеспечить уменьшение функционалов Vi(u), i Е I, необходимо при малых а выбрать параметры варьирования v Е V, x Е X^; из условия минимизации соответствующих вариаций ^Vi(u,v,x), i Е I, т. е. решить задачу
min У x(t)Av(t)Hi[t,u]dt ^ max, v Е V, x Е Xl (11)
T
Метод последовательных допустимых приближений для решения задачи (1)-(3) построен на основе процедуры игольчатого улучшения (метод игольчатой линеаризации). Пусть uk(t), t Е T, — допустимое управление; xk(t), t Е T, — соответствующая фазовая траектория; £k — положительное число. Выделим индексное множество
I*(е*) = {0} и {г = 1,т : Уг(и*) > -е*}.
Для г Е I* (е*) обозначим
А^ нг* [г] = А« (г), ж* (г), и* (г), г),
где (г) — решение сопряженной системы
V = (г),и*(г),г), = -^(ж*(¿1)).
Решение задачи (11) находится в результате декомпозиции ее на две подзадачи. Решается первая вспомогательная задача (она получается при
а = 1 (x(t) = 1)):
ielk (£fc)
min I Av(t)Hik[t]dt ^ max, v e V. (12)
T
;r.ki
Пусть uk (t), t e T, — решение этой задачи;
ök = min I Auk(t)Hlk [t]dt. ielk (£k) J T
Пусть
gk (t) = Aük(t)Hik [t] Ik (sk).
Решается вторая вспомогательная задача с параметром а e [0,1] (при v(t) = Uk (t)):
min i x(t)gk(t)dt ^ max, x e X1. (13)
ielk(£k) J T
Пусть Ха(t), t e T, — решение этой задачи;
öt = min i xt(t)gf (t)dt.
ielk(£k) J T
Тогда семейство варьируемых управлений имеет вид
ut (t) = uk (t) + xt (t)(uk (t) - uk (t)), t e T.
Переход к очередному приближению осуществляется по правилу: а) если ök < sk, то sk+1 = sk/2, uk+1(t) = uk (t);
б) если ök > Sk, то sk+i = Sk, uk+1(t) = utk (t). При этом шаг ak e (0,1] выбирается из условия уменьшения функционала при сохранении ограничений:
F°(uk+1) < F°(uk), Fi(uk+1) < 0, i = 17m.
Доказывается сходимость метода по невязке принципа максимума.
Замечание. Если в задаче (1)-(3) функциональное ограничение неактивно на управлении u e V (V(u) < 0, i = 1,m), то решение вспомогательных задач (12), (13) можно представить в явном виде. В этом случае I = {0}, и вспомогательное управление определяется условием максимума функции Понтрягина (соответствующее целевому функционалу):
u(t) = arg maxH°(^°(t), x(t), v, t), t e T.
veU
Оптимальная функция варьирования имеет вид
0, £о(£) < м«; Х«(£) = ^ 1, £о(*) > м«;
0 V 1, ро(^) = М«,
где м« — множители Лагранжа, обеспечивающие интегральное условие
/ Х«(£)^£ = «(¿1 — £0). В общем случае, когда I = {0}, вспомогательные
т
задачи требуют численного решения. Для минимаксных задач построен специализированный метод игольчатого улучшения.
Метод локального варьирования, связанный с дифференциальным принципом максимума, описывается подобным образом. Обозначим
Нк [£] = нж (¿),хк (¿),ик (£),£),
лк
и первая вспомогательная задача принимает вид:
Пусть ик(£), £ Е Т, — ее решение;
ш1п / (ник— ик(£))^£ —> шах, V Е V.
¿е4 (ей) У т
(14)
Обозначим
Пк = ш1п / (ник[¿],ик(£) — ик(£))^£.
т
к
(¿) = (ник [£],ик (£) — ик (£)); X«2 = {х Е ЫТ) : 0 < х(£) < а, £ Е Т} множество функций варьирования. Вторая вспомогательная задача имеет вид:
шт / — шах, х Е Х(
«е/й(ей ) У т
Пусть Х«(£), £ Е Т, — ее решение;
(15)
«е/й (ей)
т
Образуем семейство варьируемых управлений
ика (t) = uk (t) + ха (t)(uk (t) - uk (t)), t G T.
Переход к очередному приближению осуществляется следующим образом:
а) если щ < £к, то £к+1 = £к/2, ик+^£) = ик(£);
б) если Пк > £к, то £к+1 = £к, ик+1(£) = икак (£).
Шаг ак Е (0,1] выбирается из условия уменьшения функционала при сохранении ограничений. Доказывается сходимость метода по невязке дифференциального принципа максимума. Вспомогательные интегральные задачи (14), (15) решаются в рамках метода параметризации с помощью процедуры нелокального спуска для семейства квадратичных функционалов.
Замечание. Если на управлении и Е V все ограничения неактивны (I = {0}), то решения вспомогательных задач имеют вид
u(t) = arg max(H0[t,u],v), veU
Xa(t) = a, t G T.
В результате, получаем стандартную процедуру слабого улучшения, соответствующую методу условного градиента:
ua(t) = u(t) + a(u(t) - u(t)), t G T.
Рассмотрим теперь задачу оптимального управления (1)-(3) при ограничениях типа равенства. В этом случае, цель итерации — уменьшение функционала Лагранжа и функционала невязки типа максимума по ограничениям задачи.
Введем функционал, характеризующий невязку выполнения ограничений-равенств в задаче (1)-(3):
F (u) = max |Fi(u)|.
1 <i<m
Пусть Л = (Ао, А1,..., Am) — вектор множителей Лагранжа в задаче (1)-(3) с условием А0 G {0,1}. Введем общую сопряженную вектор-функцию
^(t,A) = Y A#(t)
¿=о
и соответствующую функцию Понтрягина
т
Н(А, х, м, ¿) = ^ АгНг(^\ х, м, ¿).
¿=0
Понятно, что А) является решением системы
т
¿=0
Причем
Н(А,^,х,м,£) = (х,м,£)).
Обозначим
Д^НМ] = Н(А,^(£,А),х(£),^),;£) -Н(А,^(£,А),х(£),м(£),;£).
Отметим, что функция Н обслуживает функционал Лагранжа задачи (1)-(3)
m
L(u,A) = ^ Аг Fi(u).
¿=o
Говорят, что управление u(t) удовлетворяет принципу максимума (ПМ) в задаче (1)-(3), если найдется вектор А = 0, А0 > 0 такой, что выполняется условие
AvH[t, А] < 0, v Е U, t Е T,
и удовлетворяет дифференциальному принципу максимума (ДПМ), если выполняется условие
(HM[t,A], v - u(t)) < 0, v E U, t E T.
Для того, чтобы выбрать параметры варьирования v E V, х E X из условия минимизации соответствующих вариаций необходимо решить задачу
J x(t)Av(t)H0[t]dt ^ max, v E V, х E X;
T
J x(t)Av(i)H'[t]dt = aFi(u), i = 1,m; (16)
T
J x(t)dt = a(ti - to). T
Здесь
X = {х Е ¿те(Т): х(г) Е {0,1}}.
Проводится декомпозиция задачи (16). Полагая а = 1 (х(г) = 1), получаем первую вспомогательную задачу на поиск управления "и(£):
J Av(t)H°(^°(t), x(t), u(t), t)dt ^ max, v Е V;
Tr . . _ (17)
Av(t)Hi(^i(t), x(t), u(t), t)dt = Fi(u), i = 1, m.
T
Пусть u(t), t Е T, — ее решение, Л = (Л°,..., Лт), Л° Е {0,1}, — соответствующий вектор множителей Лагранжа. Согласно принципу максимума для задачи (17)
u(t) = arg minH(A,^(t,A),x(t),v,t), t Е T.
veU
Образуем вспомогательные функции
g (t,A) = Aö(t)H[t,A],
gi(t) = Au(t)Hi[t], i = 1, m. Для нахождения функции варьирования x(t) решается вторая вспомо-
гательная задача
J x(t)g(t A)dt ^ max, x Е X;
T
у х(гМг)^ = г = 1,т; (18)
т
У = а(^1 - ¿о).
т
Пусть х« (¿) — решение задачи (18). Образуем а-параметрическое семейство управлений варьирования
и«(г) = и(г) + х«(г)(и(г) - и(г)), г е т.
Вспомогательные задачи (17), (18) позволяют построить семейство управлений Е V, а Е [0,1], со свойством локального (для малых а > 0) улучшения по двум функционалам ^(и), Ь(и, А):
^(и«) (и) ( ^(и) > 0),
£(иа,А) <Ь(м,Л) ( ¿(и) > 0).
Метод последовательного улучшения доступных управлений на основе процедуры слабого варьирования строится подобным образом. Слабое варьирование управления и(£), £ £ Т, проводится с помощью функций варьирования х(£), £ £ Т, из множества
= {х £ МТ) : 0 < х(£) < а, £ £ Т}.
Здесь а £ [0,1] — параметр варьирования. Задача в вариациях формулируется точно также, как в случае игольчатого варьирования:
J x(t)(HU[t], v(t) - u(t))dt ^ max, v G V, x G X2;
T
J x(t)(HU(t),v(t) - u(t))dt = aFi(u), i = 1,m.
T
Проводится декомпозиция этой задачи на две подзадачи. Полагая x(t) = 1, получаем первую вспомогательную задачу на поиск управления v(t):
/ (H00[t], v(t) - u(t))dt ^ max, v G V;
T
(19)
J (HU(t),v(t) - u(t))dt = Fi(u), i = 1,m.
T
Пусть u(t), t G T, — ее решение, Л = (Ло,..., Äm), Ло = 0 V 1, — соответствующий вектор множителей. Согласно принципу максимума
u(t) = arg min(HM[t, Л], v), t G T.
vgU
Образуем вспомогательные функции
w(t,A) = (Hu[t,A],ü(t) - u(t)),
w,(t) = (HU[t],u(t) - u(t)), i = 1,m.
Для нахождения функции варьирования x(t) решается вторая вспомогательная задача
J x(t)w(t, A)dt ^ max, x G X^2;
7 _ (20)
x(t)wj(t)dt = aFi(u), i = 1,m.
T
Поскольку функция х(£) = а является допустимой, то решение задачи (20) очевидно: х«(£) = а, £ Е Т. Тогда семейство управлений варьирования имеет вид
иа(£) = и(£) + а(и(£) - и(£)), £ Е Т.
Шаг а Е [0,1] выбирается из условия улучшения по двум функционалам ^(и) и Л).
Из приведенных результатов видно, что для варьирования управления конструируются однопараметрическое семейство управлений и формулируются специальные вспомогательные задачи. Подходящее значение параметра подбирается с помощью линейного поиска. Однако, возможности линейного поиска ограничены и приводят к медленной сходимости минимизирующей последовательности. Поэтому в [12], [33], [41], [81] используюется многопараметрическое семейство управлений. Применение многопараметрического (пространственного) поиска существенно расширяет возможности варьирования: кривая спуска (наискорейшего спуска) в фазовом пространстве отслеживается с помощью кусочно-линейной аппроксимации (не с помощью отрезков, а с помощью последовательности симплексов). При этом вместо второй вспомогательной задачи конструируется специальная задача для нахождения параметров. Кроме того, отметим, что такой способ отслеживания кривых спуска может использоваться при минимизации негладких функций конечного состояния (типа максимума вогнутых функций).
Следует подчеркнуть, что приведенный список работ далеко не полный. Однако, уже сам факт существования множества методов решения задач оптимального управления говорит о том, что нельзя однозначно установить превосходство одного алгоритма над другим. По поводу выбора алгоритма можно лишь сказать, что необходимо стараться минимизировать следующий приближенный показатель:
(степень обусловленности) х (время счета одной итерации)
(скорость сходимости)
Если обозначить скорость сходимости через г, то
г = 1 для алгоритмов сходящихся линейно,
г = 2 для алгоритмов, сходящихся квадратично, и
1 < г < 2 для алгоритмов, сходящихся сверхлинейно. Степень обусловленности — это коэффициент, который должен отражать
отклонения формы области поиска от сферической, т. е. быть малым для почти сферической выпуклой области и большим, если область имеет неудачную форму. Фактическое его значение зависит и от конкретного алгоритма. Естественно, показатель обусловленности выбирается на основе опыта. Поэтому актуальна проблема разработки новых численных методов решения задач оптимального управления, максимально учитывающих специфику решаемых задач и обладающих малой вычислительной трудоемкостью. Решению этой проблемы и посвящена работа.
Научная новизна определяется следующими результатами, полученными автором. На основе линейной интерполяции и кусочно-линейной аппроксимации функции на выпуклом компактном множестве разработаны:
— симплексный метод решения задачи минимизации выпуклой функции на выпуклом компактном множестве;
— приближенный симплексный метод и его модификация;
— метод последовательных приближений для задачи минимизации выпуклого функционала, определенного на множестве конечных состояний линейной системы управления;
— применение симплексного метода для решения линейных задач управления с линейными терминальными ограничениями;
— способ построения практически "реализуемых"управлений для линейных систем с детерминированными возмущениями и терминальными ограничениями;
— метод кусочно-линейной аппроксимации для решения краевых нелинейных задач, возникающих в результате применения принципа максимума Л.С.Понтрягина.
Результаты опубликованы в работах [11], [12], [13], [15], [14].
Работа состоит из трех разделов и введения.
В разделе 1 для задачи нелинейного программирования — задачи минимизации псевдовыпуклой функции на выпуклом компактном множестве — описывается симплексный алгоритм (точный) для нахождения решения (алгоритм 1.1). Алгоритм 1.1 генерирует последовательность симплексов, на каждом из которых в качестве вспомогательной задачи находится минимум (точный) целевой функции с помощью алгоритма 1.2. Описываются два варианта приближенного симплексного алгоритма (алгоритм 1.3), в котором отыскание локального минимума заменяется определением достаточно грубого его приближения. Приводится доказательство сходимости
симплексного алгоритма. Для задачи минимизации псевдовыпуклого функционала, определенного на множестве конечных состояний линейной системы управления, описывается алгоритм для нахождения оптимального управления в виде выпуклой комбинации экстремальных управлений вспомогательных линейных задачах Майера.
В разделе 2 для линейной системы управления с линейными терминальными ограничениями решается задача минимизации выпуклого функционала на множестве конечных состояний. Для эффективного решения вспомогательных экстремальных задач, получаемых в результате применения метода модифицированной функции Лагранжа и метода параметризации целевой функции, используется симплексный алгоритм. Дается способ построения практически "реализуемых"управлений, с помощью которых линейная система переводится из любой начальной точки в некоторую конечную точку при дополнительных ограничениях и без них. Эти управления формируются в виде линейных комбинаций управлений, вычисляемых до начала процесса управления.
В разделе 3 для решения краевых нелинейных задач, возникающих в результате применения принципа максимума Л. С. Понтрягина, описывается алгоритм, являющийся комбинацией метода неподвижной точки и квазиньютоновского метода. Использование кусочно-линейной аппроксимации функции на триангуляциях позволяет расширить область применимости алгоритма для нахождения решения краевых задач.
В каждом разделе приведены результаты счета, оформленные в виде таблиц. Программы, по которым производились расчеты, можно найти на интернетсайте http://viboldirev.narod.ru
1 Симплексный алгоритм для решения экстремальных задач
В данном разделе рассматривается задача минимизации псевдовыпуклой функции на выпуклом компактном множестве. Предлагается симплексный алгоритм для нахождения решения. Приводится доказательство его сходимости. Дается описание приближенного симплексного алгоритма. Рассматривается задача оптимального управления со свободным правым концом. Критерием качества является псевдовыпуклый функционал, определенный на множестве конечных состояний линейной системы управления. Предлагается алгоритм для нахождения решения. Приводятся результаты счета. Основные результаты данного раздела опубликованы в работах [12], [13].
Задача минимизации квадратичной формы на выпуклом компактном множестве Q в п-мерном евклидовом пространстве была решена Гильбертом [104], который предложил находить минимум квадратичной формы с помощью вычисления локальных минимумов на подходящим образом выбранной последовательности линейных сегментов во множестве Q. Барр [96] обобщил итеративную процедуру Гильберта так, что минимизация производится на последовательности выпуклых многогранников (метод выпуклых оболочек, см. так же [41], [81], [33]). В [110], как альтернатива процедуре Барра, был предложен симплексный алгоритм, в котором на каждом шаге минимизация выполняется на выпуклой оболочке п + 1-ой точки (т. е. на п-симплексе). Алгоритмы, в которых аппроксимация решений осуществляется с помощью симплексов, в дальнейшем будем называть симплексными.
В данном разделе решается задача минимизации псевдовыпуклой функции ^(ж), определенной на выпуклом компактном множестве Q евклидова пространства [12]. Поиск минимума функции обычно в этом случае осуществляется фактически с помощью опорных к множеству Q гиперплоскостей (внешнее представление множества). Любая же точка выпуклого множества является линейной комбинацией не более п + 1 крайней точки этого множества (внутреннее представление). Поэтому целесообразно использовать именно это внутреннее представление множества при поиске минимума функции конечного состояния. В то же время, в результате решения вспомогательной задачи линейного программирования получается крайняя точка множества, тем самым, полученное решение вспомогательной задачи является "мостом"между внешним и внутренним представлениями множе-
ства.
Предлагаемый в [12] алгоритм решения задачи минимизации псевдовыпуклой функции на выпуклом компактном множестве строится на основе методов [106]. В работе [106] для нахождения максимума псевдовогнутой функции (см. [108], стр. 140-144]), заданной на многограннике предлагается алгоритм, который находит точный максимум целевой функции за конечное число шагов. Этот алгоритм использует внутреннее представление допустимого множества, разлагая его на симплексы переменной размерности. Используя линейную программу, алгоритм вырабатывает последовательность симплексов, которые содержат локальные относительно симплексов максимумы. Получаемая последовательность локальных максимумов строго возрастающая, и наибольший из них является глобальным максимумом на допустимом множестве. Нахождение локальных максимумов на этих симплексах сводится к максимизации целевой функции на аффинных многообразиях (см. [71], стр. 19-29).
Предлагаемый алгоритм вырабатывает с помощью решения вспомогательной задачи линейного программирования симплексы в допустимом множестве. На этих симплексах осуществляется минимизация псевдовыпуклой целевой функции. При этом предполагается, что функция конечного состояния должна допускать решение в замкнутой форме уравнений, вытекающих из необходимых условий оптимальности. В идейном плане способ нахождения локального относительно симплекса минимума такой же, как и в [106]. В отличие от [106] процесс нахождения глобального на множестве достижимости минимума будет уже бесконечным. Однако, использование симплексов в процессе минимизации делает вычислительную процедуру достаточно эффективной.
Одним из способов решения задачи оптимального управления с помощью принципа максимума Л. С. Понтрягина [68] является сведение ее решения к решению краевой задачи (см., например, [70], [61], [66], [84]). Для решения краевой задачи, возникающей в результате применения принципа максимума, предложены различные приближенные методы (методы прогонки, метод Ньютона, квазиньютоновские методы, методы градиентного типа — проекции градиента, условного градиента [61], [66], [84], [32], [56], [41], [81]). При нахождении минимума функционала на множестве конечных состояний линейной системы управления решение исходной задачи сводится к решению серии вспомогательных линейных задач оптимального управления (задач Майера), к нахождению минимума целевой функции и к
определению допустимого управления, переводящего систему в найденную точку минимума. В случае решения задачи минимизации нормы конечного состояния управляемой системы решение исходной задачи сводится к вычислению минимума квадратичной формы на компактном выпуклом множестве Q в n-мерном евклидовом пространстве.
В данном разделе рассматривается задача минимизации псевдовыпуклого функционала ^(ж), определенного на множестве Q конечных состояний линейной системы управления. Решение исходной задачи сводится:
1) к решению задачи нелинейного программирования ^(ж) ^ min (к
xeQ
нахождению вектора ж Е Q), а затем:
2) к нахождению управления, приводящего траекторию системы в точку ж.
Приводимый в этом разделе алгоритм решения задачи минимизации функционала на множестве конечных состояний линейной системы управления строится на основе методов [33], [106]. В работе [33] описывается алгоритм для решения задачи минимизации нормы конечного состояния управляемой системы. Оптимальное управление в исходной задаче находится в виде выпуклой комбинации оптимальных управлений вспомогательных задач Майера, поиск минимума квадратичной функции (нормы конечного состояния) на множестве достижимости состоит в последовательном решении задач квадратичного программирования — задач по определению точки выпуклой оболочки (симплекса), ближайшей к началу координат (метод выпуклых оболочек, см. также [41], [81]). Задача нахождения минимума квадратичной функции на симплексе решается с помощью конечно-шаговой процедуры, на каждой итерации которой определяется безусловный минимум квадратичной функции.
В работе [106] для нахождения максимума псевдовогнутой функции (см. [108], стр. 140-144]), заданной на многограннике предлагается алгоритм, который находит точный максимум целевой функции за конечное число шагов. Этот алгоритм использует внутреннее представление допустимого множества. Предлагаемый в [12] алгоритм вырабатывает с помощью решения вспомогательной задачи Майера симплексы во множестве достижимости. На этих симплексах, в отличие от [33], осуществляется минимизация псевдовыпуклой целевой функции.
1.1. Основные понятия и определения
Основным инструментом в нашем построении является кусочно-линейная аппроксимация (РЬ-аппроксимация) отображения / : Ип ^ Кп. Эта аппроксимация определяется через значения f в вершинах п-мерных симплексов.
Пусть Р = {х°,...,хп} — множество (аффинный базис). Аффинная оболочка множества Р является линейным многообразием М = {х : х =
п
Ра, а^ = 1}. Если точка х Е М, то она может быть представлена в виде
г=°
п
х = Ра, ^^а^ = 1. (1.1)
г=°
Выпуклая оболочка аффинно независимых точек множества Р представляет собой п-мерный симплекс (п-симплекс): а = {х : х = Ра, а^ > 0, г =
_ п
0,п, а» = 1}.
г=°
Величина diam а = шах^х1 — х2||то : хг Е а, г = 1,2} называется диаметром симплекса. Говорят, что отображение Ьрх = Ах + с является линейной аппроксимацией отображения / : Ип ^ Ип, если /(хг) = Ахг + с, г = 0, п. Здесь А — действительная невырожденная п х п-матрица, с — действительный п-вектор.
Пусть точка х лежит в некотором симплексе а (не обязательно единственном) , и поэтому имеет барицентрическое представление относительно вершин симплекса а:
пп
х = а^хг, ^^а^ = 1, а > 0, г = 0,п. (1.2)
г=° г=0
Значение Ьр в точке х Е а дается соотношением
п
Ьр х = ^ аг/(хг), (1.3)
г=°
где а^ удовлетворяет (1.2). Если х лежит в пересечении двух или более симплексов, то х лежит в относительной внутренности некоторой ]-грани т. Для любого ] Е {0,... ,п} ]-грань т в а есть выпуклая оболочка некоторых (^ + 1) вершин симплекса т: т = [хг°,... ]. Только ненулевым а^ в формуле (1.2) соответствуют вершины т, которые являются общими для
любого симплекса т, содержащего x. Таким образом, определение отображения Lp зависит только от x и не зависит от выбора частного симплекса, содержащего x. Отображение Lp явно обладает интерполяционным свойством.
Пусть AX = (x0 - xn,...,xn-1 - xn), AF = (f (x0 - xn),...,f (xn-1 - xn)).
n-1
Из альтернативного представления формул (1.1) и (1.3) x = xn + ^ a*(x* -
¿=0
n-1
xn) и LPx = f (xn) + a*(f(x*) - f (xn)) следует, что
¿=0
Lpx = AF(AX)-1x + (f (xn) - AF(AX)-1xn). (1.4)
Триангуляция пространства Rn есть набор n-симплексов T = {а* : i = 1, 2,...}, объединение которых покрывает пространство Rn, элементы которого пересекаются попарно по одной (самое большее) общей j-грани, причем пересечение этого объединения с любым ограниченным множеством состоит из конечного числа симплексов.
Величина mesh T = sup {diam а : а G T} называется мелкостью триангуляции T. Приведем в качестве примера триангуляцию Фройденталя [107], [94]. Эта триангуляция пространства Rn определяется точкой x0 и имеет мелкость, пропорциональную заданному положительному числу 5. Будем обозначать ее Tn[x0,5] и определять следующим образом. Вершины триангуляции являются узловыми точками множества x0 + 5Zn, где Z обозначает множество целых чисел. Тогда n-симплекс а может быть определен вершиной z0 G Zn и перестановкой п множества {1,... , n}, а именно: а = [x0,..., xn], x0 = x0 + 5z0 и x* = xi-1 + 5en(i), i = 1, n. Здесь e1,..., en — стандартный ортонормированный базис. Диаметр любого симплекса а в Tn[x0, 5] есть 5, так что мелкостью триангуляции Фройденталя является также 5.
Кусочно-линейное отображение Ft : Rn ^ Rn на триангуляции определяется значениями отображения f в вершинах симплексов триангуляции T. Внутри симплексов всего пространства Rn отображение Ft определяется интерполяцией. Значит,
n
Ftx = ^af (x*), x G а G T, (1.5)
¿=0
где а* определяется в силу (1.2). Аналогично формуле (1.4) имеем:
Ftx = AF (AX )-1x + (f (xn) - AF (AX )-1xn). (1.6)
Известно, что идея метода секущих для нахождения решения системы нелинейных алгебраических уравнений основано на замене нелинейного отображения f линейным интерполирующим отображением. Нуль линейного отображения принимается в качестве очередного приближения к решению нелинейной системы, т.е. на (k + 1)-й итерации xk+1 является решением линейной системы Lpx = 0. Используя представления (1.3), (1.4) и применяя стандартную нумерацию точек в зависимости от номера итерации, получим следующие два вида метода секущих
n
xk+1 = ^ агхк-\ k = 0,1,..., (1.7)
¿=0
где коэффициенты разложения а находятся в результате решения линей-
nn
ной системы af (xk-i) = 0, ^ а = 1;
¿=0 ¿=0
xk+1 = Xk - AXk(AFk)-1f (xk), k = 0,1,... . (1.8)
Методы кусочно-линейной аппроксимации (PL-методы) успешно используются для решения системы нелинейных алгебраических уравнений f (x) = 0 или для нахождения неподвижных точек x таких, что f (x) = x. Показано, что определенные реализации этих методов имеют квадратичную скорость сходимости, когда они подходят достаточно близко к решению, и в этих случаях (без использования производных отображения) вычислительные затраты сравнимы с затратами метода Ньютона [114], [117].
Можно показать, что если f имеет якобион f'(x), который непрерывен по Липшицу в окрестности симплекса а, то существуют константы K1 и K2 такие, что
||f (x) - Ftx|| < K1(diam а)2 и ||f (x) - AF(AX)-1|| < K2(diam а),
где x E а и нормы являются sup-нормами. Для того, чтобы гарантировать локальное существование решений уравнения, необходимо предположить регулярность отображения f, т.е. матрица f'(x) должна иметь полный ранг. Аналогичное условие для Ft состоит в том, что AF имеет полный ранг.
При работе PL-метода генерируется последовательность симплексов ак, диаметры которых £k = diam ak стремятся к нулю, и на этих симплексах находятся решения xk линейной системы Ftx = 0 (см. (1.5), (1.6)), т.е.
n n _
системы af (xik) = 0, ^ a^ = 1, a^ > 0, i = 0,n (здесь xik — вершины
¿=0 ¿=0
симплекса ). Подчеркнем дополнительно, что симплексы принадлежат некоторой триангуляции. В гомотопном методе, в методе непрерывной деформации эти n-симплексы являются n-гранями (n + 1)-симплексов. Обычно, такие грани называют полностью размеченными или полными. В (n + 1)-симплексе имеются две полностью размеченные n-грани. Этот (n + 1)-симплекс, содержащий полностью размеченную грань, и сама n-грань называются поперечными. Таким образом, PL-алгоритмы могут рассматриваться как непрерывная деформация Fk (кусочно-линейной аппроксимации отображения f на симплексах a"k триангуляции) к F= f и как отслеживание пути из нулей xk аппроксимации Fk при k ^ то. При определенных условиях на отображение f можно показать, что путь xk сходится к нулю отображения f (см., например, [113], [114], [117]).
Рассмотрим теперь применение метода линейной аппроксимации и метода кусочно-линейной аппроксимации для решения задач оптимизации. Пусть задана псевдовыпуклая функция Rn ^ R1. Функция D С Rn ^ R1 называется псевдовыпуклой на выпуклом множестве D0 С D (см. [108], стр. 140-144]), если для любых G D0 имеет место
— z) > 0 ^ ^(ж) > ^(z).
Очевидно, что задача безусловной минимизации ^(ж) —> min сводится к
' xgR"
решению нелинейной системы алгебраических уравнений g(x) = ^'(ж) = 0. Решение этих уравнений можно осуществить указанными выше методами секущих и PL-методами.
Рассмотрим задачу нелинейного программирования
^(ж) ^ min, (1.9)
xgQ
где ^(ж) — псевдовыпуклая функция, а Q — замкнутое выпуклое компактное множество конечномерного пространства V. Пусть V* — сопряженное пространство. Определим отображение z : V* ^ V следующим образом:
z(g) = arg min(g, ж). (1.10)
xgQ
В алгоритмах оптимизации именно это отображение служит фундаментом для построения минимизирующей последовательности. При решении задачи нелинейного программирования вида (1.9) V = Rn, g = ^'(y), g G V*, y G Q. Для задачи оптимального управления пространство V является пространством фазовых координат, а элемент g G V* определяется краевым условием ^(i1) = —g сопряженной системы (в частности, для задачи минимизации функции ^(ж) конечного состояния g = ^'(y), y G Q).
1.2. Алгоритм решения задачи нелинейного программирования
Рассмотрим задачу нелинейного программирования
^(ж) ^ min, (1.11)
xgQ
где ^(ж) — псевдовыпуклая функция, а Q — выпуклое компактное множество. Ниже приводится алгоритм решения задачи (1.11), являющийся обобщением метода условного градиента. Суть метода условного градиента заключается в том, что на каждом шаге проводится линеаризация функции ^(ж) в точке yk, решается вспомогательная линейная задача минимизации
(^'(yk), ж) ^ min,
xgQ
а минимизатор zk G dQ в этой задаче определяет направление движения для получения следующего приближения. Таким образом,
(^'(yk),zk) < (^'(yk),ж) для всех ж G Q, zk G dQ, yk+i = + 0Ä(zk - ), о < < 1, где — решение задачи одномерной минимизации
+ 0(zk - yk)) ^ mm .
Более эффективным методом для решения задачи (1.11) является симплексный метод (метод выпуклых оболочек [96], [33], [41]). Он может быть интерпретирован как обобщение метода условного градиента (решается вспомогательная задача минимизации целевой функции не на отрезке, а на симплексе). Пусть в результате вычислений на k-й итерации получены крайние точки z0,...,zq, q < n, множества Q. Множество (выпуклая оболочка q + 1 точки)
( q q
ak = < ж G Q : ж = ^^ о^, > 0, i = 0, q, ^^ а^ = 1
I ¿=о ¿=о
есть q-мерный симплекс. На каждом шаге симплексного метода проводится линеаризация функции ^(ж) в точке yk, решается вспомогательная линейная задача минимизации (^(yk),ж) ^ min, а минимизатор zk "замещает"в
xgQ
симплексе a"k вершину zго. Эта точка является противолежащей вершиной
по отношению к грани, в которой находится точка yk. В дальнейшем базис, определяющий такую грань минимальной размерности, будем обозначать Sk. Полученный таким образом очередной симплекс ak+1 и симплекс a"k имеют общую грань, в которой находится точка yk. Такие симплексы называются смежными, а каждый из них — поперечным. Следующее приближение yk+1 получается в результате решения задачи минимизации вида ^(ж) ^ min . Таким образом,
xefffc+1
(^'(yk),zk) < (^'(yk),ж) для всех ж Е Q, ^(yk+1) < ^(yk).
Обозначим P = {z0,..., zq} — аффинный базис, Д(Р) = [z0,..., zq] — выпуклая оболочка P (замкнутый симплекс или просто симплекс), ri a — внутренность симплекса a относительно его афинной оболочки.
Алгоритм 1.1 (решение задачи нелинейного программирования):
0. Задать z Е Q. Вычислить g0 = ^'(z). Найти крайнюю точку z0 Е Q,
решив задачу (g0, ж) ^ min. Положить y0 = z0, P1 = {z0}, a1 = [z0], k =1.
xeQ
1. Вычислить gk = ^'(yk-1). Найти крайнюю точку zk Е Q, решив задачу (gk, ж) ^ min.
xеQ
2. a) Если (gk, zk — yk-1) = 0, процесс прекращается;
b) Если (gk, zk — yk—1) < 0, то дополнить базис Pk вектором zk: Pk = Pk U {zk}, ak = Д(Pk). Перейти к следующему шагу.
3. Найти решение yk задачи нелинейного программирования ^(ж) ^ min, приняв в качестве начальной точки yk—1. Положить Pk+1 = Sk, a^+1 =
xеafc
Tk, k = k + 1. Перейти к шагу 1.
Приведенный алгоритм является прямым методом в том смысле, что получаемая им алгоритмическая последовательность определяется на допустимом множестве. Он вырабатывает строго убывающую последовательность {^(yk)} значений функции ^(ж). В то же время алгоритмом генерируется последовательность смежных симплексов {ak}, для каждого из которых yk является локальным минимизатором функции ^(ж), т.е.
^(yk) = min ^(ж). (1.12)
Вершинами симплекса ak являются крайние точки zk множества Q, получаемые с помощью решения вспомогательной задачи математического
программирования (^'(yfc),x) ^ min. Симплекс определяется аффинным
xgQ
базисом Pk = {z0,... , zq}, 0 < q < n. Но так как локальный минимум yk может лежать на грани размерности p < n — 1 симплекса (в симплексе ak меньшей размерности, чем размерность ), то некоторые весовые коэффициенты (по меньшей мере, 1) в барицентрическом представлении yk будут равны 0. Пусть
yk = Pk w = z0w0 + ... + zqwq.
Тогда можно записать
= {zj,j G J(w)}, Tk = A(Sk),
где J(w) = {j : Wj = 0, j = 0,q}. Очевидно, что Sk С Pk, Tk С ak.
Пусть у(т) — кривая спуска, проходящая через точки yk. Тогда можно сказать, что эта кривая отслеживается последовательностью смежных симплексов, а сама последовательность симплексов образует кусочно-линейное многообразие. Таким образом, кривая спуска содержится в этом многообразии и аппроксимируется ломанной, полигонами которой являются отрезки [yk—1,yk]. Такой способ задания кривых является уникальным. Известно, что, если минимум ^(ж) на Q достигается в точке X, то условие (^'(X), ж — X) > 0 для всех ж G Q является необходимым условием минимума (и достаточным в выпуклом случае). Эквивалентно (^'(yk—1),z — yk—1) < 0 для некоторого z G Q означает, что yk—1 не является минимизатором на Q (см. шаг 2. b) алгоритма 1.1). Поэтому величина £(ж) = sup(^'(x),x — z)
zgQ
может быть неотрицательной и равной нулю в точках, в которых выполняется условие экстремума. Следовательно, она может служить мерой близости к минимуму. Справедливо следующее утверждение (см. теорему 3.1 из [106]).
Лемма 1.1. Пусть ^ является псевдовыпуклой, и пусть z — решение задачи (^'(y), ж) ^ min. Тогда y — минимизирует ^ на Q в том и только
xgQ
том случае, когда имеет место равенство
(p'(y),z — y) = 0.
Лемма 1.2 ([43], стр. 48) (необходимое условие глобального минимума функции на множестве). Пусть ^(х) — выпуклая и непрерывно дифференцируемая функция; Q — выпуклое компактное множество. Точка у Е ^
является решением задачи (1.11) тогда и только тогда, когда
(^(y),z - y) > 0 z G Q.
1.3. Минимизация функции на симплексе
В предыдущем пункте для нахождения минимизатора функции ^(ж) на выпуклом множестве Q необходимо было решать (шаг 3 алгоритма 1.1) вспомогательную задачу
^(ж) ^ min . (1.13)
Здесь а — q-симплекс, а С Rn (0 < q < n).
Пусть начальная точка z = Pw, Р = {z0,...,zq}, а = Д(Р), M = {ж : ж = Ра,е*а = 1} — q-многообразие, где e = (1,..., 1)*. Решение задачи (1.13) сводится к решению конечного числа (не превосходящего n) задач минимизации функций ^(ж) на многообразиях Mp уменьшающейся размерности p < n — 1.
Алгоритм 1.2 (решение задачи минимизации на симплексе):
1. Найти решение ж' задачи минимизации
^(ж) ^ min, (1.14)
жем
ж' = Pw'.
a) Если > 0 для всех г = 0, д, то ж' С п а является минимизатором ж функции ( на а. Полагаем у = ж', $ = Р, процесс прекращается.
b) Если по меньшей мере одно эд' < 0, то ж' С п а. Получили ((ж') < ((ук-1). Перейти к следующему шагу.
2. Определить точку жр пересечения отрезка ж'] с границей симплек-
$ + ... + zqwp.
са а. Пусть жр = Pwp = z0wp + ... + zqwp. Получим жр = z + Л (ж' — z),
0 < Л < 1, wp = Aw' + (1 — A)wj > 0, i = 0,q, причем для некоторого j wp = 0, A =( max
j ' V т/ \ W-j
v «gJ(w) j
3. Удаляем из Р вершины zj, ^ С J(^р), которые не являются носителями жр, что приводит к уменьшению аффинного базиса. Полученный базис Рр порождает подсимплекс ар С а.
a) Если ар имеет нулевую размерность, у = жр является его единственной точкой и, таким образом, относительным локальным минимизатором ^ на ар. Полагаем у = жр, Б = Рр, процесс прекращается.
b) Если ар имеет размерность > 0, полагаем г = жр, Р = Рр, М = Мр, а = ар, переходим к шагу 1.
В дальнейшем будут полезны следующие леммы (см. [106], лемма 3.4).
Лемма 1.3. Пусть ^ псевдовыпукла и пусть ^(ж') < ^(г). Тогда любое жр = Аж' + (1 — А)г, 0 < А < 1, удовлетворяет неравенству ^(жр) <
).
Доказательство. Допустим противное, что ^(жр) > ^(г). Отсюда в силу псевдовыпуклости следует: (^'(г),жр — г) > 0. Преобразуя выражение для жр, получим ж' — г = -(жр — г), 1 < - < Домножим левую и правую части равенства скалярно на ^'(г):
(^'(г),ж' — г) = -^'(г), жр — г) > 0. А
Из этого неравенства в силу псевдовыпуклости функции ^ следует, что <р'(ж') > ^(г). Получили противоречие, которое доказывает лемму.
Лемма 1.4 [106]. Пусть у и г такие, что ^(ж) > ^(у) для всех ж Е [у, г]. Тогда (^'(у),г — у) > 0.
Решение задачи (1.13) условной минимизации является конечно-шаговой последовательной процедурой нахождения критических точек функционала ^(ж) на многообразиях Мр (уменьшающейся размерности), содержащих грани симплекса . Отметим, что удобнее осуществлять минимизацию на многообразиях в пространстве барицентрических координат, так как используется внутреннее представление множества (. Если найденная на шаге 1 алгоритма 1.2 критическая точка ж' Е п а, то цель итераций при переходе от минимизатора ук—1 Е п С ( к минимизатору ук Е п ак С ( достигнута, а именно ук = ж', ^(ук) < ^(ук—1). Если же точка ж' Е г! а, то находится жр — точка пересечения отрезка [ж', ук—1] с границей а С М. Дальнейшее удаление из базиса Р ненесущих вершин не влияет на положение жр в а, так как в этом случае изменяется только барицентрическая система координат (и> заменяется на и>р). Точка же жр находится в п ар С а, ар С а, то есть в подсимплексе, соответствующем подбазису Рр С Р.
По лемме 3.1 жр удовлетворяет неравенству ^(жр) < ^(ук—1) и, таким образом, новый минимизатор ж' на соответствующем подмногообразии Мр
будет удовлетворять неравенству ^(ж') < ^(жр) < ^(yk_1). Если ж' G ri ар, то итерация завершилась успешно. В противном случае, ж' G ri ар процедура повторяется при а = ар, Р = Pp, M = Mp.
Заметим, что наиболее трудоемким в алгоритме 1.2 является шаг 1 — нахождение критической точки функции ^(ж) на многообразии. Очевидно, что min^(ж) = min^(Ра). В правой части последнего равенства берется
xGM а
минимум по а при условии е*а = 1. Можно легко избавиться от ограничения е*а = 1. Пусть
.,0 _,q _,q-1 _,q\ о f „
^q-1
D = (z0 - zq,..., zq-1 - zq), в = (ао,..., аq-l)*.
n-1
Тогда ап = 1 - ^ вг, а задача (1.14) преобразуется к задаче без огра-
г=о
ничений ^(De + zq) ^ min. Поэтому для нахождения критической точки функции ^(ж) нужно решить относительно в систему уравнений.
DV'(De + zq) = 0. (1.15)
В силу сказанного, ясно что для большей эффективности алгоритма необходимо чтобы система (1.15) допускала решение в замкнутой форме. Это имеет место, например, если ^(ж) — положительно определенная квадратичная функция.
1.4. Сходимость симплексного алгоритма
Утверждения о сходимости алгоритма 1.1 приводятся в следующей теореме.
Теорема 1.1 (сходимость алгоритма 1.1). Пусть Q — выпуклое компактное множество в Rn, ^(ж) — псевдовыпуклая функция на
Rn, {yk},
k = 0,1..., — последовательность, генерируемая алгоритмом 1.1. Тогда:
1) если последовательность {yk} конечна (т.е. {уо,у1,...,у1 , у1+1}), то у1 является решением задачи (1.11);
2) если последовательность {yk} бесконечна, то каждая ее точка сгущения является решением задачи (1.11), т.е. ^(yk) монотонно убывает и lim (^'(yk-1),yk-1 - zk) = 0.
k^to
Доказательство. Доказательство сходимости минимизирующей последовательности проводится по схеме доказательства теоремы 3 [109] с учетом специфики симплексного алгоритма.
Предположим, что алгоритм генерирует только конечную последовательность у0,... , у1, у1+1. Тогда ((у1+1) = ((у1). Действительно, на (1 + 1)-й итерации на шаге 1 алгоритма будет получена крайняя точка г1+1, то есть
(('(у1 ),г - ¿1+1) > 0 г Е д. (1.16)
Пусть у1+1 Е I = [у1, г1+1] такое, что
((у1+1) < ((г) г Е I.
Значит,
^(у1+1) < ((у1). (1.17)
В силу же пункта 3 алгоритма
((у1+1) < ((г) г Е а/+1.
В то же время, у1+1 Е [у1, г1+1] С 07+1. Поэтому ((у1+1) < ((у1+1). Учитывая (1.17) и предположение ((у1+1) = ((у1), получаем ((у1+1) = ((у1) < ((г) для всех г Е I. Но в силу леммы 4 это означает, что
((/(у1 ),^+1 - у1) > 0. (1.18)
Из (1.16) и (1.18) следует, что
(('(у1 ),г - у1) > 0 г Е д
и, следовательно, в силу леммы 2, у1 — решение задачи (1.11).
Предположим теперь, что алгоритм генерирует бесконечную последовательность {ук}, и пусть у будет точкой сгущения этой последовательности. Тогда существует неограниченное подмножество К целых чисел такое, что подпоследовательность {ук}к сходится к у.
Пусть ук Е = [ук-1,гк] такое, что ((ук) < ((г) для всех г Е . Но множество д компактно по предположению, поэтому существует неограниченное подмножество К С К такое, что подпоследовательности }к, {ук}к, {ук+1}к сходятся к некоторым точкам г, у, у соответственно. Тот факт, что эти последовательности сходятся, обеспечивается тем, что последовательность {ук}к сходится к у. В силу непрерывности отображения ('(•) и скалярного произведения справедливы следующие соотношения:
(('(у), г) < (('(у), г) г Е д, (1.19)
((у) < ((г) г Е I* = [у,4 (1.20)
У Е Г, (1.21)
< ^(У). (1.22)
Предположим теперь, что ^(у) < ^(у) — 5 для некоторого 5 > 0. Тогда существует I такое, что ^(ук+1) < ^(ук) — 5/2 для всех к > I в К, то есть монотонно убывающая последовательность {^(ук)} не ограничена снизу. Но это невозможно, поскольку ^ — непрерывное отображение, а Q компактно. Отсюда следует, что не существует 5 > 0 такое, что ^(у) < ^(У) — 5, то есть
> ^(У). (1.23)
Формулы (1.20)—(1.23) означают, что ^(у) = ^(у) = ^(у) < для всех г Е I*. В силу же леммы 4 имеем
— у) > 0.
Но тогда из формулы (1.19) следует, что (^'(у), г — у) > 0 для всех г Е Q, и поэтому у — решение задачи (1.16).
Теорема доказана.
Отметим, что теорема 1.1 свидетельствует о том, что последовательность локальных на а"к минимизаторов сходится, в силу выпуклости (псевдовыпуклости) функции к глобальному на Q минимизатору этой функции (см. лемму 2).
1.5. О приближенном решении вспомогательных задач
Симплексный алгоритм, описанный в пункте 1.2, требует на каждой итерации решения вспомогательной задачи — минимизации функции ^(ж) на симплексе. Но точное решение этой задачи само может быть предельной точкой итерационного процесса. Следовательно, симплексный метод нуждается в модификации шага 3.
Решение вспомогательных задач на шаге 3 симплексного метода может быть связано в общем случае с большими вычислительными затратами. В частном же случае, используя конкретный вид функции ^(ж) (например, ^(ж) = ||ж||), задачу (1.13) на шаге 3 можно решать приближенно, а именно, можно решать задачу (1.14) лишь при р = д. Это же означает, что вовсе не обязательно решать задачу (1.13) до конца, а можно ограничиться определением отделяющей опорной гиперплоскости. Далее до конца пункта будем полагать ^(ж) = ||ж||.
Пусть
Y = L(P) = (z0,... , zq) — матрица, определяющая симплекс;
M = r(P) = {ж G Rq—1 : x = Yw, e*w = 1} — аффинная оболочка P, т.е. q-мерное многообразие. Пусть z' = Y7 G а С M;
J(Y) = {j : Yj = 0, j = S = P\{zj'}, где j' : Yj' = 0, т.е. j' G J(Y); t = A(S) — (q — 1)-мерная грань симплекса а; Z = L(S) — матрица, определяющая грань т.
Пусть p = q, ^(x) = ||ж||. В этом случае критическая точка ж' = 0. Поэтому для барицентрического представления ж' необходимо решить систему линейных алгебраических уравнений вида Ya = 0, e*a = 1. Однако, для вычисления очередной крайней точки допустимого множества Q можно использовать не градиент функции ^(ж) в точке минимума (найденного на предыдущей итерации), а его некоторое приближение. В качестве такого приближения можно взять, например, вектор нормали к гиперплоскости, содержащей общую грань смежных симплексов. Напомним, что именно в этой грани находится точка минимума, определенная на предыдущей итерации. Пусть гиперплоскость описывается уравнением (ж, w) — 1 = 0. Так как вершины грани лежат в этой плоскости, то для определения нормали получаем систему алгебраических уравнений Z*w — e = 0. Приведем упрощенный вариант приближенного симплексного алгоритма.
Алгоритм 1.3 (приближенный симплексный алгоритм):
0. Задать Z G Q. Вычислить g0 = z/||z||. Найти крайнюю точку z0 G Q,
решив задачу ^°,ж) ^ min. Положить v0 = z0, P1 = {z0}, а1 = A(P1),
xgQ
p =1, в0 = 1; А = 0, i = 1,n; k = 1.
1.
a) если p < n, вычислить gk = vk—1/||vk—1||;
b) если p = n найти решение wk—1 G Rn сопряженной системы линейных алгебраических уравнений:
Zk*_ 1w — e = 0.
Положить gk = wk—1/|wk—1|;
c) определить крайнюю точку zk G Q, решив задачу (gk, ж) ^ min. Если
XgQ
p < n, положить p = p +1.
2. Положить £к = (дк, — -ик Если £к = 0, то процесс прекращается; в противном случае дополнить базис Рк вектором , а именно:
Рк = Рк и {¿к},
П = ¿(Рк),
д = |Рк | (число элементов Рк),
^к = Д(Рк),
Мк = Г(Рк). Перейти к следующему шагу.
3. Найти критическую точку функции ^(ж) на многообразии Мк, т.е. решить систему линейных алгебраических уравнений относительно а:
Ук а = 0, е*а = 1.
Если аг > 0 для всех i = 0, q, то ж' = Ykа £ ri ak. В этом случае
a) если q < n, то положить vk = ж', вг = аг, i = 0,q, Sk = Pk, Zk = L(Sk), k = k + 1, перейти к шагу 1;
b) если q = n, процесс прекращается, так как произошел "захват"гло-бального минимизатора симплексом; если же какое-либо а^ < 0, перейти к следующему шагу.
4. Определить точку z' пересечения отрезка [x',vk-1] с границей симплекса ak. Пусть z' = Yky и z' = ж' + A(vk—1 — ж'), 0 < А < 1; y« = аг + А(вг — аг) > 0, i = 0,q, причем для некоторого j' имеем Yj' = 0; А = max .
5. Удалить из базиса Pk точку zj', j' £ J(y).
Положить Sk = Pk\{zj"}, Zk = L(Sk), Pk+1 = Sk, ak+i = Tk, vk = z', в = Y, k = k + 1; перейти к шагу 1.
Сходимость последовательности {^(vk)} значений функции к минимуму гарантируется теоремой 1.1, так как выполнены все необходимые для нее условия. Основное требование теоремы строгого убывания значений ^(vk) обеспечивается шагами 3 и 4 приведенного алгоритма (см. также [12], лемма 2]).
Приведем еще одну версию приближенного симплексного алгоритма, основанную на линейной интерполяции отображения z : Rn ^ Rn, а именно z(g) = arg min(g,x). В силу шага 1 c) алгоритма 1.3 zk = z(gk), k =
xGQ
0,1,.... Поэтому, получив на шаге 4 алгоритма 1.3 коэффициенты разло-
жения 7«, г Е J(7), на шаге 1 Ь) можно не решать сопряженную систему, а положить и>к—1 = ^ 7«^. Остальные операции алгоритма 1.3 остаются
«е7 (7)
без изменений.
1.6. Постановка задачи оптимального управления
Рассмотрим следующую задачу оптимального управления. На траекториях управляемой системы
Ж = A(t)x + B(t)u, t е T = [t0,ti], x(t0) = ж0, u(t) e U (1.24)
минимизировать функционал
F (u) = ^(x(ti)). (1.25)
Здесь ж — n-вектор фазового пространства V; u — s-вектор управления; A(t), B (t) — n x n- и n x s-матрицы соответственно.
Будем считать, что матрицы A(t), B(t) определены при t е T, вещественны и непрерывны. Предположим, что функция ^(ж) непрерывно дифференцируема по ж и псевдовыпукла на Rn, U С Rs — выпуклое компактное множество. Предположим, что система (1.24) полностью управляема. Класс U допустимых управлений определим как множество кусочно-непрерывных на T вектор-функций u(t) со значениями в U. Обозначим через Q множество достижимости системы (1.24) в момент t = t1.
Согласно [55] (стр. 178, теорема 1.1) Q — выпуклый компакт. Обозначим через ж(^и) решение системы (1.24) при фиксированном допустимом управлении u(t). Пусть u(t) — допустимое управление, удовлетворяющее условию максимума
(V>(t,u),A(tb(t,u) + B(t)v) ^ max, t е T. (1.26)
veU
Здесь ^(t,u) — решение сопряженной системы
^ = —A*(t)^, ^(ti) = —р'(ж(иЛ)). (1.27)
Заметим, что управление u(t) является оптимальным управлением в
задаче минимизации линейного функционала (д,ж(^^)) ^ min на тра-
veU
екториях линейной системы (1.24) при g = , u)) и при фиксирован-
ном управлении u(t). Отсюда следует, что конец траектории ж(^й) (вектор
ж(^,й)) является решением задачи минимизации
— ж(^,и)) ^ min .
Очевидно, что z = ж(^,й) — крайняя точка множества Q. В силу принципа максимума нахождение оптимального решения {и^),ж(^} сводится к
решению серии вспомогательных линейных задач оптимального управления (задач Майера) вида
(g, ж(£1, u)) ^ min,
ugU
причем (1.27) принимает вид
0 = —A*(t)0, 0(tO = —g. (1.28)
Таким образом, решение исходной задачи (1.24), (1.25) сводится:
1) к решению задачи нелинейного программирования ^(ж) ^ min (к
xgQ
нахождению вектора ж G Q), затем
2) к нахождению управления u(t), приводящего траекторию системы (1.24) в точку ж, т.е. ж(^1,и) = ж.
1.7. Симплексный алгоритм решения задачи оптимального
управления
В этом пункте будет дано формальное обоснование симплексного алгоритма для решения задачи (1.24), (1.25). Подробное его описание для решения задачи нелинейного программирования, эквивалентной исходной задаче, приведено в пункте1.2.
Рассмотрим задачу минимизации выпуклого функционала F(u) = ^(ж(^1,и)), u G U, на множестве конечных состояний системы (1.24). В сокращенной форме эту задачу можно представить в виде
^(ж) ^ min, (1.29)
xgQ
где ^(ж) — выпуклая и непрерывно дифференцируемая функция, Q — выпуклое компактное множество.
Введем функцию
H(0^,v,t) = (0,А^)ж + B(t)v), t G T, v G U. Определим экстремальное управление с помощью соотношения
(0(t,u), A(tk(t,u) + B(t)u(t)) = max(0(t,u), A(tk(t,u) + B(t)v), tGT,
vgU
где 0(t,u) — решение сопряженной системы:
0 = —A*(t)0, 0(t1,u) = —(ж(t1,u)).
Или, иначе
(^(t,u),B(t)u(t)) = max(^(t,u),B(t)v), t G T. (1.30)
vGU
Заметим, что управление u(t) одновременно является оптимальным управлением в задаче минимизации линейного функционала
(('(x(t^u)),x(ti,v)) ^ min
vgU
на траекториях линейной системы (1.24), где u(t) — фиксированное управление.
Одним из методов решения задачи (1.29) является метод условного градиента. Пусть на k-й итерации уже построены: uk (t) — допустимое управление, xk(t) = x(t, uk), uk(t) — экстремальное управление (1.30). Идея метода состоит в том, что для уменьшения целевого функционала используется решение задачи одномерной минимизации:
((yk + 0(zk - yk)) ^ min , (1.31)
где yk = x(t1,uk), zk = x(t1,uk). Пусть (0 — решение задачи (1.31), yk+1 = yk + 0(zk — yk). Тогда управление, приводящее траекторию системы (1.24) в точку yk+1, имеет вид
uk+1(t) = uk (t) + 0(uk (t) — uk (t)).
Более эффективным методом для решения задачи (1.24) является симплексный метод (метод выпуклых оболочек [96], [33], [41]). В работах [96], [33], [41] симплексный метод рассматривается для квадратичной функции ((ж) вида ((ж) = ||ж — ж'||2, где ж' — фиксированный вектор.
Идея симплексного метода заключается в том, что для уменьшения целевого функционала ((ж) используется решение задачи многопараметрической минимизации — задачи минимизации целевой функции на симплексе (выпуклой оболочке n + 1 точек множества достижимости). Пусть в результате вычислений на k-й итерации получены крайние точки множества достижимости z0k = x(t1,u0k),...,znk = x(t1,unk). Так как любая точка выпуклого множества Q является выпуклой комбинацией не более чем n + 1 крайних точек из этого множества, в качестве параметров можно взять коэффициенты барицентрического представления точек этого множества. Поэтому введем вектор параметров а = (a0,...,aq)*. Обозначим
Yk = (z0k,...,znk). Пусть
ak = < z G Rn : z = Ykа, а > 0, i = 0, q, ^^ а^ = 1 l ¿=0
— q-мерный симплекс, ^k С Q. Тогда уменьшение целевой функции ^(ж) на множестве Q может быть достигнуто в результате решения задачи минимизации вида
^(z(a)) = min ^(z). (1.32)
В результате решения задачи (1.32) будет найдена точка yk+1 = z(а) = Ykа, а управление, приводящее траекторию системы (1.24) в точку yk+1 будет иметь вид
q
uk+1(t) = ^ аг uik (t).
¿=0
Таким образом, рассматриваемый метод работает в классе экстремальных управлений, удовлетворяющих условию максимума функции Понтрягина.
Пусть
ДиН(t) = Н(^(t,u),x(t,u),ü(t),t) - Н(^(t,u)x(t,u),u(t),t).
В силу (1.25) ДиH(t) >0, t G T. Обозначим
ti
5(u) = J ДиН (t)dt.
to
Если £(u) = 0, то ДиН(t) = 0 почти для всех t G T, то есть u(t) удовлетворяет принципу максимума Понтрягина. Поэтому величину £(u) можно рассматривать как критерий выполнения (нарушения) необходимых условий оптимальности.
Обозначим: U(t) — оптимальное управление; u(t) = x(t,u); Xk(t) = x(t,uk); (t) = ^(t,uk). В силу выпуклости ^(ж) имеем:
0 < F(uk) - F(и) = ^(x(ti,uk)) - ^(x(tbw)) <
< (^/(x(t1,uk)),x(t1,uk) - x(t1,u)) < < (^(x(tbuk)),x(t1,uk) - x(t1,uk)) = 4.
Однако, ^k можно записать в иной форме:
* -/н,*,, — „ + — ж.„+
to
tl
+ (/, B (t)uk — B (t)uk )]dt = J , B (t)uk — B (t)uk )dt.
to
Из неравенств 0 < F(uk) — F(u) < и из условия lim = 0 следует, что
k^to
последовательность {uk(t)} минимизирующая. То, что ^k ^ 0 при k то, утверждается в теореме пункта 1.4. Таким образом, основываясь на (1.32), можно построить симплексный алгоритм для решения задачи оптимального управления (1.24), (1.25). Подробное описание алгоритма для решения экстремальных задач дано в пунктах 1.2 и 1.3. Поэтому остановимся здесь на специфике алгоритма для решения задач оптимального управления.
Для работы алгоритма в начале определяется приближение z конечного состояния системы с помощью задания нулевого приближения управления u(t), t G T, т.е. вычисляется z = ж(^,й). После этого находятся крайние точки zk множества достижимости Q как результат действия отображения z : V* ^ V, т.е. zk = z(gk), gk G V*. Здесь gk определяет краевое условие ^(t1) = —gk для сопряженной системы. Это наиболее трудоемкий по объему вычислений шаг в алгоритме. Для задачи оптимального управления он состоит из следующих действий (определяющих отображение z(g)):
1) вычисление решения ^ (t) сопряженной системы (1.27) при gk = ('(yk—1), а именно ^(t1) = —gk, и нахождение при этом экстремального управления uk (t) из условия максимума (1.26);
2) вычисление решения системы (1.24) при u = uk(t) и нахождение крайней точки zk = z(gk) = ). На следующем шаге алгоритма находится локальный относительно симплекса ak минимизатор yk G ak. Этот поиск можно осуществить [12] с помощью нахождения критических точек функционала ((ж) на многообразиях, содержащих грани симплекса ak.
Итак, пусть на некоторой (к — 1)-й итерации получены: yk—1, Pk — аффинный базис, ak = A(Pk) — симплекс. Тогда на очередной к-й итерации алгоритма выполняются следующие операции.
1. Вычисляется gk = ('(yk—1), крайняя точка zk G Q в результате решения задачи (gk, ж) ^ min.
xgQ
2. Если (gk, zk-yk 1) = 0, то процесс прекращается; в противном случае базис Рk дополняется вектором zk: Рk = Рk U {zk}, ak = Д(Рk).
3. Отыскивается решение yk задачи нелинейного программирования ^(ж) ^ min, принимая в качестве начальной точки yk-1; процесс повторя-
xgafc
ется с п. 1.
1.8. Численные примеры
В этом пункте приводятся в качестве иллюстрации работы алгоритма 1.1 результаты численных расчетов двух тестовых примеров. Примеры решались с точностью £ = 10-4 (равенство нулю скалярного произведения проверялось с точностью до £, см. шаг 2 алгоритма 1.1). За единицу трудоемкости работы алгоритма 1.1 принимается "итерация", равная объему вычислений при интегрировании сопряженной и прямой системы, т. е. объему вычислений, требуемых для получения очередного приближения (очередной крайней точки множества достижимости).
Обозначения:
k — номер итерации;
gk — краевое условие для сопряженной системы (см. (1.28));
zk = ж^,t1) — крайняя точка множества достижимости;
Mgk) = (gk ,zk).
Пример 1.1.
ж 1 = ж2,
= жз,
ж3 = и;
ж(0) = (1, 0,0)*, T =[0, 3], U = [-1,1],
F =
= ^(ж) = 1+(ж,ж) ,
где (ж, ж) = ж2 + ж2 + ж3. Начальное приближение мо(^) = 0, t G T.
Результаты счета:
F = 0.006352,
g= (0.0816, -0.1191,0.0638)*,
х = (0.0413, -0.0603,0.0323)*, и(г) = {а : г е (Г"1, Г], г = 1,7},
г 0 1 2 3 4 5 6 7
гг 0.0 0.7561 0.7743 0.7863 2.2614 2.2988 2.3269 3.0
аг — —1.0 —0.7828 0.1396 1.0 0.0775 —0.1396 —1.0
Ход итераций иллюстрируется таблицей 1.1. Таблица 1.1
к 9к г к к)
1 0.5 0 0 —3.5 —4.5 —3 —1.75 0.371429
2 0.466939 —0.323265 —0.21551 —1.29737 —0.97866 0.753047 —0.451717 0.181074
3 0.296572 —0.48208 0.278262 1.57109 1.18975 0.500291 0.0315992 0.0457256
4 0.233283 —0.29256 0.137586 —1.36322 —1.35445 —0.491784 0.0105798 0.0158861
5 0.148127 —0.19017 0.0495741 0.407119 0.580816 0.957537 —0.00267961 0.0142739
6 0.126143 —0.178173 0.0838308 0.118448 0.110583 0.264196 0.0173863 0.0118349
7 0.0975113 —0.155281 0.109809 —0.698147 —1.10755 —0.869661 0.00840837 0.00894662
Продолжение Таблицы 1.1
8 0.0898542 -0.139606 0.0852702 0.21512 -0.0357704 -0.126697 0.0135197 0.0079771
9 0.0943983 -0.130487 0.0722009 -0.984196 -1.08686 -0.530831 0.0105886 0.00746415
10 0.0887548 -0.12816 0.0699338 -0.282145 -0.378542 -0.145752 0.013279 0.00706527
11 0.0843498 -0.126153 0.0680933 0.356225 0.219152 0.15536 0.0129797 0.00661472
12 0.084435 -0.121868 0.0628975 0.0142296 -0.0526816 0.0812663 0.0127332 0.00649501
13 0.0811125 -0.12022 0.0666701 0.0408247 -0.0918303 -0.0258179 0.01263 0.00636339
14 0.0811994 -0.119569 0.0638795 0.151056 0.0433215 0.0853687 0.012539 0.0063519
15 0.0816549 -0.119104 0.0638224 -0.0207134 -0.116907 0.00495144 0.0125487 0.0063519
Пример 1.2. (См. [106].)
ж 1 = ж2 + и, ж2 = -и,
ж 3 ж1;
ж(0) = (0.5,0, 0)*,
Т =[0,1.5],
и = [-1,1],
1
Р = = (с,х) + (ж*Аж) 2 , где с = (0,0,1)*,
0.38 0.39 0.35 А = V= 0.39 0.42 0.35 , 0.35 0.35 0.34
0.2 0.1 0.3 V = 0.3 0.4 0.3 . 0.5 0.5 0.4
Начальное приближение:
и0(г) = 0, г е Т.
Результаты счета: Р = 1.226594;
Ж = (0.6114, 0.6423, 1.5644)*; ж = (-0.0197, 0.7887, 0.4679)*;
Ж(г) = {а : г е (гг-1,гг], г = !~3},
где
г 0 1 2 3
гг 0.0 0.3187 1.4631 1.5
аг — 1.0 —1.0 1.0
Ход итераций иллюстрируется таблицей 1.2.
Таблица 1.2
к С* гк Мс * *)
1 0.610845 0.617595 1.58047 -0.0626435 0.926013 0.429424 1.21233 1.22856
2 0.610208 0.644187 1.56084 -0.0197915 0.788765 0.467918 1.22638 1.22659
3 0.611419 0.642348 1.56445 -0.0251575 0.803247 0.464054 1.22657 1.22659
2 Решение линейной задачи оптимального управления с терминальными ограничениями
Для линейной системы управления с терминальными ограничениями рассматривается задача минимизации выпуклой функции конечного состояния. Задача оптимального управления редуцируется к задаче нелинейного программирования — минимизировать выпуклую функцию на множестве при дополнительных ограничениях. Используя метод модифицированной функции Лагранжа, решение этой задачи сводится к решению последовательности задач минимизации выпуклой функции на выпуклом компактном множестве. Для эффективного решения полученных вспомогательных экстремальных задач предлагается использовать симплексный алгоритм.
В данном разделе описывается способ построения финитных управлений для линейных систем. Приводится конкретный вид управления для случая линейных систем с детерминированными возмущениями и систем с терминальными ограничениями. Получаемые финитные управления являются линейными комбинациями управлений, переводящих за фиксированное время единичные орты Яп в 0 для линейных нестационарных систем. Такие финитные управления удобны в практических приложениях, так как они могут быть вычислены до начала процесса управления. Это свойство управлений является весьма важным при реализации быстродействующих процессов. В качестве иллюстрации приводятся расчеты тестовых примеров. Данный раздел написан на основе работ [11], [13].
Создание эффективных вычислительных методов теории оптимального управления является одним из приоритетных направлений в теории оптимальных процессов. Эти методы разработаны и показали свою эффективность в задачах со свободным правым концом (см., например, [58], [66], [26]). Для задач оптимального управления с краевыми условиями появляются дополнительные затруднения.
Отметим основные подходы к численному решению задач оптимального управления для систем, описываемых обыкновенными дифференциальными уравнениями:
1) методы последовательных приближений, основанные на процедурах линеаризации и варьирования управлений [58], [84], [76], [74];
2) методы штрафных функционалов или модифицированного функционала Лагранжа [58], [84], [38], [23];
3) дискретизация задачи с последующим применением методов линейного и нелинейного программирования [84], [38].
Основная трудность при разработке методов первого типа состоит в построении множеств варьирования управления. Рассматриваемый ниже метод примыкает ко второй группе и использует идеи математического программирования применительно к задачам оптимального управления. При использовании методов выпуклого программирования решение исходной задачи обычно определяется как предел минимизирующей последовательности решений более простых (вспомогательных) экстремальных задач. В данном разделе предлагается алгоритм, рассчитанный на использование для решения вспомогательных задач быстросходящегося метода (симплексного алгоритма). Подробное описание симплексного алгоритма приводится в пунктах 1.2 и 1.3 раздела 1 для решения задачи нелинейного программирования, а в пунктах 1.6 и 1.7 раздела 1 — для решения задач оптимального управления со свободным правым концом.
2.1. Постановка задачи
Рассмотрим управляемую систему
х = А(г)х + В(г)и, г е Т =[£о,£1], х(го) = хо , и(г) е и. (2.1)
Здесь х — фазовый п-вектор; и — й-вектор управления; и — выпуклое компактное множество, А(г), В (г) — п х п- и п х й- матрицы соответственно. Будем считать, что матрицы А(г), В (г) определены при г е Т, вещественны и непрерывны.
Предположим, что система (2.1) полностью управляема. Класс и допустимых управлений определим как множество кусочно-непрерывных на Т вектор-функций и(г) со значениями Ь е и. Обозначим через х(г,и) решение системы (2.1) при фиксированном допустимом управлении и(г).
Задача. Минимизировать функционал
^о(и) = ^о(х(г1, и)), и е И, (2.2)
на траекториях системы (2.1) при дополнительном условии
х(г1, и) е Я = {х : Сх + А = 0}, (2.3)
где г(х) = Сх + А е Ят — заданная вектор-функция терминального ограничения.
Предположим, что функция ^0(ж) выпуклая и непрерывно дифференцируемая, а для ограничений (2.3) выполняются условия регулярности: компоненты вектор-функции г(ж) имеют линейно-независимые градиенты.
Обозначим через Q множество достижимости системы (2.1) в момент t = t1: Q = {x(t1,u) : u G U}. Как известно, при сделанных предположениях Q — выпуклый компакт. Перепишем исходную задачу (2.1)-(2.3) в сокращенной форме:
^0(ж) ^ min, Сж + d = 0, ж G Q, (2.4)
и сделаем дополнительные предположения, связанные с корректным переходом к двойственной задаче. Каждому вектору z = (z1,..., zm)* сопоставим множество R(z) = {ж : Сж + d = z} и предположим, что:
1) int Q П R(z) = O;
2) Q П R(z) = O в некоторой окрестности точки z0 = 0.
Искомое решение и задачи (2.4) может быть определено как предел последовательности минимумов модифицированной функции Лагранжа [58], то есть решение исходной задачи (2.1)-(2.3) может быть определено как предел минимизирующей последовательности решений задач оптимального управления со свободным правым концом. Применение метода модифицированной функции Лагранжа приводит к необходимости решать задачи минимизации выпуклой функции на выпуклом компактном множестве. Для этой цели в данной работе предлагается использовать симплексный алгоритм.
Основу симплексного алгоритма для решения задач оптимального управления составляют: свойство суперпозиции линейных систем; сим-плексность; замещения; сходимость последовательности локальных минимумов вспомогательных задач к глобальному. В результате работы симплексного алгоритма оптимальное решение и задачи (2.4) получается в виде выпуклой комбинации точек выпуклого компактного множества (множества достижимости). Оптимальное управление, приводящее систему в эту точку, в силу линейности системы (2.1), также получается в виде линейной комбинации некоторых вспомогательных (полученных на предыдущих итерациях вычислительного процесса) управлений. Поэтому, если оптимальное управление содержит участки особого управления, то такие участки в приближенном кусочно-постоянном управлении получаются автоматически, что является достоинством предлагаемого подхода. Приближенный симплексный алгоритм, приведенный в пункте 1.4, в сочетании с
теоремой об отделимости выпуклых множеств в конечномерном пространстве может служить основой создания алгоритмов для решения многих задач оптимального управления (в том числе и с подвижными концами).
2.2. Модифицированная функция Лагранжа
Рассмотрим в качестве исходной задачу
^o(x) ^ min, Cx + d = 0, x G Q. (2.5)
Введем функцию Лагранжа L(x,A) = ^0(x) + (A,Cx + d) на Q x Rm. Известно (см., например, [10, стр. 111]), что если функция Лагранжа L, соответствующая задаче (2.5), имеет на Q x Rm седловую точку (x,A), то x является решением задачи (2.5). Однако, алгоритм, осуществляющий непосредственный поиск седловой точки функции Лагранжа сходится весьма медленно. Это объясняется тем, что радиус кривизны поверхности |^0(x(A)), r(x(A))} может оказаться слишком большим. Здесь x(A) = arg min L(x,A) — точка множества Q, в которой заданная функ-
xgQ
ция L достигает своего минимума на Q, а множители Лагранжа A задают параметризацию этой поверхности.
Определим
Li(A) = inf L(x, A) и
xgQ
L2(x) = sup L(x, A).
Agßm
Очевидно, что L1(A) < L2(x) при всех x G Q, A G Rm.
Лемма 2.1 ([35], стр. 112). Точка (x,A) G Q x Rm является седловой для функции L тогда и только тогда, когда L1(A) = L2(x).
Введем в рассмотрение модифицированную функцию Лагранжа
M(x, A,a) = <p0(x) + (A, Cx + d) + ^||Cx + d||2 на Q x Rm
2
(a G R1, a > 0 — коэффициент штрафа).
Задача
M(x,A,a) ^ min (2.6)
xgQ
неэквивалентна (2.5), но аппроксимирует ее с тем большей точностью, чем больше a. Коэффициент штрафа a в отличие от метода штрафов уже не
обязательно должен быть достаточно велик. Слагаемое а||г(х)||2 имеет целью лишь обеспечить выпуклость упомянутой поверхности с не очень большим радиусом кривизны.
Аналогично тому, как это делается с обычной функцией Лагранжа, определим
Mi (Л, а) = inf М (х,А,а), (2.7)
xgQ
а также
М2(х,а) = sup M(х,А,а).
Лбйт
Используя модифицированную функцию Лагранжа M, тогда можно записать
M2(x, а) ^ min
xgQ
и в качестве двойственной — рассмотреть задачу
Mi (А, а) ^ max . (2.8)
Лей™
При этом имеет место неравенство [10, стр. 116]
sup М1(А,а) < inf М2(х,а).
Известно достаточно большое количество алгоритмов (см., например, [38], [35]), базирующихся на решении двойственной задачи (2.8). Предложенные алгоритмы отличаются друг от друга способом решения вспомогательной задачи минимизации (2.7), а также способом пересчета двойственных переменных А и параметра а. Алгоритм данного раздела в основе содержит алгоритм [72], но отличается от него и от остальных алгоритмов тем, что для решения вспомогательной задачи (2.7) используется симплексный алгоритм [12]. Подробное описание его приводится в разделе 1.
Отметим, что ради устойчивости вычислительного процесса в алгоритме используется масштабирование в = 1, а вместо функции M(x, А, а) — функция
M(x, А, в) = вМ(x, А, а), в = 1.
а
Очевидно, что функции М(х,А,в) и М(х,А,а) имеют одни и те же седло-вые точки. Приведем алгоритм отыскания седловой точки модифицированной функции Лагранжа M(x, А, в).
Пусть заданы константы (определяемые опытным путем) С1, П; 7 £ (0,1), 5 £ (0,1) и известны на (1 + 1)-м шаге итерационного процесса х1, Л1, в/. Определим решение х/+1 задачи
M(ж, Л1 ,ß) ^ min .
xgQ
(2.9)
Для отыскания решения х/+1 задачи (2.9) используем симплексный алгоритм. Обозначим через ¿/ число итераций симплексного алгоритма, необходимых для определения х/+1. Эта величина используется в качестве характеристики трудоемкости решения вспомогательной задачи (2.9) симплексным алгоритмом.
Полагаем
Л/+1 =
Л1+ß/ г(ж/+1), г/ >п, Л1+ß/ г(ж/+1), г/ < п, Л1, г/ < п,
и и
Определяем
ß/+1 =
^ß/,
ß/,
|г(ж/+1)||2 < y ^г(ж/ )||2, |г(ж/+1)||2 >y |г(ж/ )||2.
г/ < п,
(2.10)
(2.11)
Для полученных значений ß/+1 вновь решается задача (2.9). Таким
образом общая схема алгоритма имеет следующий вид:
1. Задать жо, Ло, ß0, г'. Положить l = 0.
2. Определить решение ж/+1 задачи (2.9) симплексным алгоритмом.
3. Если |г(ж/)|| < г', вычисления прекращаются. Если ß/ < c1, вычисления прекращаются, так как система ограничений несовместна.
4. Вычислить Л/+1 по формуле (2.10).
5. Вычислить ß/+1 по формуле (2.11).
6. Положить l = l + 1 и перейти к пункту 2.
Заметим, что при поиске решения двойственной задачи (2.8) необходимо вычислять значения функции (2.7). Однако такое вычисление само по себе требует решения экстремальной задачи (2.9). Для ее решения предлагается использовать симплексный алгоритм.
2.3. Метод параметризации целевой функции
В предыдущем пункте исходная задача оптимального управления
F0(u) ^ min при дополнительном условии ж(^, u) G R или в иной форме
ие U
^о(ж) ^ min редуцировалась к серии задач без дополнительных ограничений ^(ж(^1,м)) = M(ж(^1 , и)Л,а) ^ min при фиксированных Л, а или
ugU
^(ж) = M(ж, Л, а) ^ min при фиксированных Л, а. В данном пункте рас-
xgQ
сматривается сведение исходной задачи к другой серии задач без дополнительных ограничений. Введем функцию
N (ж,т) = [^о(ж) — т]2 + ||Сж + d||2. (2.12)
Пусть при фиксированном т(т < ^(ж)) решается т-задача:
^(ж) = N(ж,т) ^ min . (2.13)
xgQ
Очевидно, что при f = ^(ж) имеем N(ж(г)) = 0. Пусть ж(т) — решение задачи (2.13) при фиксированном т.
Лемма 2.2 ([38], стр. 185). Пусть т1 < т2 и ж(т1), ж(т2) — решения т-задачи. Тогда ^0(ж(т1)) < ^0(ж(т2)).
Из леммы следует, что значения целевой функции ^0(ж(т)) является монотонно возрастающей функцией переменной т. Предположим, что известна некоторая нижняя оценка т0 оптимального значения ^0(ж(т)), т.е. т0 < ^0(ж). Такую оценку можно получить, например, численно решив вспомогательную задачу
^0(ж) ^ min (2.14)
xgQ
или F0(u) ^ min. Метод параметризации целевой функции состоит в после-
ugU
довательном нахождении точек ж(т&) и увеличении параметра таким образом, что стремится к ^0(ж). Ключевым моментом метода является способ изменения параметра т от итерации к итерации. Морисоном Д. предложен следующий вариант.
Пусть на k-й итерации известно некоторое значение параметра < ^0(ж). Из решения вспомогательной задачи (2.13) найдем точку ж(т^). Если N(ж(т&),т&) = 0, то расчеты заканчиваются. В противном случае положим
т^+1 = т^ + y/N (ж(т* ),т*), (2.15)
и итеративный процесс продолжается. Доказательство сходимости итеративного процесса приводится в [38].
Отметим, что в методе параметризации целевой функции величина Р = %0(Х)—т] является коэффициентом штрафа. Метод параметризации поэтому, по существу, отличается от метода внешних штрафных функций
только тем, что в нем автоматически определена политика изменения вспомогательного параметра т, в то время как в методе внешних штрафных
функций пользователь должен специально определять правило изменения р.
Понятно, что эффективность метода в целом связана с возможностью "быстрого"решения т-задачи (2.13). В связи с этим, для решения вспомогательных задач (2.13), (2.14) предлагается использовать симплексный алгоритм, описанный в разделе 1.
2.4. Построение финитных управлений для линейных систем без
ограничений
Рассмотрим сначала построение управлений в линейных управляемых системах в простейшем случае. Пусть задана система линейных обыкновенных дифференциальных уравнений вида:
х = А(г)х + В(г)и + /(г), г е [го,г1], (2.16)
где х — фазовый п-вектор; и — й-вектор управления, А(г), В (г) — п х пи п х й-матрицы соответственно. Будем считать, что матрицы А(г), В (г) и вектор ](г) определены при г е [го,г1], вещественны и непрерывны.
Пусть заданы вектора хо, х1 — начальное и конечное состояние системы (2.16). Предполагается, что матрицы А(г) и В (г) образуют полностью управляемую пару {А(г),В(г)}. Допустимыми управлениями будем считать кусочно-непрерывные функции вида
и(г) = N (г)А, (2.17)
где А е Яп — постоянный вектор, N (г) — й х п-кусочно-непрерывная ограниченная на [го,г1] матрица.
Пусть х(г) — решение системы (2.16) при управлении (2.17). Требуется найти вектор-функцию и = и(г), чтобы решение системы (2.16), начинающееся при г = го в точке х(го) = хо, попадало при г = г1 в точку х(г1) = х1.
Задача нахождения оптимального управления, переводящего линейную систему в некоторое конечное состояние, исследована достаточно полно. Однако, численная реализация предложенных алгоритмов решения этих задач сопряжена в силу их итерационности со значительными затратами машинного времени. Это является одним из препятствий применения таких алгоритмов для быстропротекающих процессов. Предлагаемый способ
построения финитных управлений позволяет преодолеть, в некоторой степени, эту трудность. Задачу нахождения управлений вида (2.17) условно можно разбить на две: нахождение матрицы N(г), затем вектора Л. Столбцами матрицы N(г) являются управления, переводящие за время ¿1 — г0 единичные орты Яп в начало координат. Такие управления могут быть вычислены до начала процесса управления. Параметр же Л вычисляется достаточно просто.
Непрерывные управления вида (2.17) были введены в 1960 году Р. Калманом [42] при рассмотрении достаточных условий управляемости линейной системы, а также для задачи минимизации энергии управления (см., например, [6]).
Кусочно-непрерывные управления вида (2.17) были получены в работе [2] при рассмотрении задачи оптимального быстродействия в результате решения трансцендентных уравнений. При ограничении управлений (2.17) по амплитуде в [3] описывается многогранная область начальных условий для системы управления, из которой возможен перевод системы в начало координат за фиксированное время. Кусками поверхности этого многогранника являются (п — 1)-симплексы, вершины которых лежат на осях координат. Здесь приводится простая и общая схема получения финитных управлений.
Предлагаемый способ построения финитного управления базируется на свойстве суперпозиции линейных систем. Возьмем в качестве ]-го столбца матрицы N(г) управление V7 (г), ] = 1,п, переводящее линейную систему
// = А(г)у + В (гу, у(*о) = е (2.18)
за время Т = г1 — г0 из точки е7 в 0
N7(г) = N(г)е7 = V7(г), ; = Т"П.
Здесь е7 — единичный орт, ^ (г) — ]-й столбец матрицы N (г). Так как {А(г),В(г)} полностью управляемая пара, то такие кусочно-непрерывные управления существуют [42], [44].
Пусть Ф(г) — фундаментальная матрица решений системы
// = А(г)/,
причем Ф(г0) = Е — единичная матрица. В силу указанного выше выбора
матрицы N(t) имеет место соотношение:
ti
Ф^У + J Ф(^)Ф-1(£)В(t)Nj(t)dt = 0,
to
отсюда
ti
ej + J Ф-1ВДВ(t)Nj (t)dt = 0
to
или в матричной форме:
ti
E + ^ Ф-1(^В(t)N(t)dt = 0. (2.19)
to
Отметим, что матрица N(t) определяется указанным способом неоднозначно. Эта неоднозначность может быть устранена дополнительными требованиями
F (Nj) ^ min , где F — некоторый функционал.
Определим теперь параметр Л для управления (2.17). Подставим управление (2.17) в систему (2.16), запишем ее решение в форме Коши, при t = t1, положив
ti
У1 = / Ф(^)Ф-1(0/(t)dt = 0.
to
В результате получим:
ti
ж1 = Ф(t1)ж0 + Ф(^)У Ф-1(^В (t)N (t)dU + y1.
to
Но в силу соотношения (2.19)
ti
J Ф-1(^В(t)N(t)dt = -E,
to
поэтому
ж1 - У1 = Ф(tl)(жо - Л)
отсюда
Л = хо - Ф-1(£1)(Ж1 - ух).
Таким образом, справедлива следующая
Теорема 2.1. Пусть система (2.16) полностью управляема. Тогда финитное управление, переводящее систему (2.16) из точки х0 в точку хх за фиксированное время ¿х — ¿0, имеет вид:
и(*) = N(¿)(хо + Ф—х(^х)(ух - хх)), ¿1
Ух = / Ф(^х)Ф-1(^ = 0, (2.20)
¿о
причем ] -м столбцом матрицы N(¿) является управление, переводящее систему (2.18) за время ¿х - ¿0 из точки е-7 в начало координат (здесь е-7 — ^-й единичный орт).
Приведем конкретный вид финитных управлений (2.20) для случая, когда матрица N(¿) удовлетворяет некоторым дополнительным требованиям.
Пусть качество управления вспомогательной системы (2.18) оценивается функционалом:
¿1
Р(V) = / /о(и(т))^г.
¿0
Задача состоит в нахождении управлений V-7 (£), которые переводят систему (2.18) из состояния у(£0) = е-7 в состояние у (¿а) = 0 и минимизируют критерий качества:
F ) = min F ), j = 1,n.
Положим
H(^,y,v) = (Ay(t) + Bv(t),^(t) - fo(v(t)),
где ^(t) — удовлетворяет дифференциальному уравнению
^ = -A*(t)^. (2.21)
Пусть v(t) — оптимальное управление, а y(t) — соответствующая оптимальная траектория. Тогда согласно принципу максимума Л. С. Понтрягина [68] найдется непрерывная вектор-функция ^(t) ф 0, являющаяся решением уравнения (2.21) и такая, что выполняется условие максимума
H(^^(t),y(t),V(t)) = (A(t)y(t),^(t)) + maX[(ß(t)v,#)) - fo(v)]. (2.22)
1. Предположим, что столбцы матрицы N (г) должны удовлетворять требованиям:
2J Nj*(t)Nj(t)dt ^ min, j = 1,n. (2.23)
to
Это означает, что необходимо построить оптимальное управление для системы (2.18), переводящее ее из точки y = ej в точку y = 0, за время T = ti — t0 и доставляющее функционалу
ti
2
F (vj) = i vj*(t)vj (t)dt (2.23')
¿о
наименьшее возможное значение. В силу условия максимума (2.22) оптимальное управление V7 имеет вид:
V' (г) = в *(г)^' (г), (2.24)
где (г) — решение сопряженной системы (2.21) с некоторым начальным условием 0(г0) = с7, подлежащим определению, (г) = Ф—1*(г)с7. Пусть Д(г) = Ф—1(г)В(г). Тогда оптимальное управление запишется в виде:
V (г) = я*(г)с'. (2.25)
Подставим управление (2.25) в систему (2.18). Так как это управление переводит систему из точки е- в 0 за время г1 — г0 имеем:
Ф(г1)е7 = у Ф(г1)л(г)л*(г)с7'^ = 0.
¿0
Отсюда, обозначив
¿1
я(т) = у я(г)я*(г)^г, ¿0
получим
с7' = — £(Т )е7.
Таким образом, оптимальное управление (2.24) примет вид:
v
j (t) = —R*(t)D—1(T )ej.
t
i
t
i
Но так как по построению является j-м столбцом матрицы NV(t), получим
NV(t) = -R*(t)D-1(T), а соответствующее финитное управление (2.20) имеет вид:
u(t) = —R(t)D—1 (T)[xo + Ф-1(^1)(у1 - xi)]. (2.26)
Отметим, что финитное управление (2.26) является оптимальным управлением для системы (2.16), переводящим ее из точки x0 в точку x1 и минимизирующее функционал
ti
F(u) = 2 J u*(t)u(t)dt.
to
Однако, это свойство теряется (как видно из последующих пунктов 2-4), если на управление вспомогательной системы (2.18) наложены дополнительные ограничения по амплитуде.
2. Пусть столбцы матрицы N(t) удовлетворяют требованию (2.23), причем
N(t)l< j, i = М; j = 1^, (2.27)
где /j — константы. В этом случае необходимо найти удовлетворяющее ограничению (2.27) управление для системы (2.18), переводящее ее из точки y = ej в y = 0, за время t1 — t0 и доставляющее минимум функционалу (2.23'). В силу (2.22) оптимальное управление в данном случае имеет вид:
(t) = ^(B*(t)^j(t),/j), j = 1~n, ^ G Rs
для почти всех t G [t0,t1]. Здесь функция Rs x Rs ^ Rs означает следующее:
/ j\ I zi , zi G ( 1 г , 1 г);
^(z,/) = 4 , . , ч
[ / i sign (Zi) ,
/г > 0, i = 17s, z, / G Rs.
В данном случае начальное условие сопряженной системы (2.21) ^(t0) = cj, j = 1,n, может быть найдено по сравнению с предыдущим пунктом уже алгоритмически [6]. Тогда финитное управление (2.20) имеет вид:
n
u(t) = ^ Aj<^(B*(t)^j(t),/j) (2.28)
j=1
для почти всех Ь £ [Ь0 , ¿1],
Л = жо + Ф-1(Ь1)(у1 - ж1).
Приведенные в пунктах 1, 2 финитные управления (2.26) и (2.28) являются непрерывными ограниченными функциями. В следующих двух пунктах приводятся финитные управления, являющиеся кусочно-постоянными функциями.
3. Предположим, что столбцы матрицы N(Ь) удовлетворяют требованиям
1 а
JE IN(t)|dt ^ min, j = 1,n,
—>
\ - ■ г / \
, г=1
¿0
а также выполняются ограничения (2.27). Поэтому управления ^ (Ь), минимизирующие функционал
1 а
з) = I VА из
F(vj) = / Е Iv1 (t)|dt, j = 1,n,
to i=1
и переводящие систему (2.18) из y = ej в 0 за время t1 — to, имеют вид:
(t) = ^o(B*(t)^j(t),1j), j = 1,n, ^o G Rs для почти всех t G [to ,t1]. Функция Rs x Rs —> Rs означает следующее
/Л ) 0, Z G (—1,1); li Sign Zi ,
1г > 0, г = 1,5, г £ Яа, 1 £ Яа.
В данном случае управление (2.20) имеет вид:
п
и(Ь) = £ Лз ^о(В (Ь),1з)
з=1
для почти всех Ь £ [Ь0 , Ь1],
Л = жо + Ф-1(Ь1)(у1 - ж1).
4. Предположим, что столбцы матрицы N(Ь) должны быть "минимальны по амплитуде". Смысл этих слов таков: найдутся числа > 0, ^ = 1,п,
при которых заданное время T = ti — to является оптимальным по быстродействию при переводе системы (2.18) управлениями вида (t) = ^(t) из точки y = ej в 0 (см., например, [88], [41], [81]). При этом считается, что управления (t) стеснены ограничениями по амплитуде вида
IV(t)i < , i = м; j = 17n
где l > 0 E Rs — const. Подобная задача является обратной в смысле Н. Н. Красовского по отношению к задаче быстродействия с теми же ограничениями. Оптимальные управления v(t) и числа ^ могут быть найдены алгоритмом из [88], [41], [81].
В силу принципа максимума Л. С. Понтрягина оптимальное управление (t) в данном случае имеет вид
(t) = ^ sign (B*(t)^j(t)), j = 17П.
Тогда финитное управление (2.20) имеет вид:
n
u(t) = ^ Ajsign (B*(t)^j(t)) j=i
для почти всех t E [t0 , ti],
A = xo + Ф-i(ti)(yi — xi).
Таким образом, для построения финитных управлений достаточно найти управления до начала процесса регулирования, запомнить их точки переключения (если они существуют), а затем во время процесса управления сформировать финитные управления. В следующем пункте рассмотрено построение финитных управлений для линейных систем с различными терминальными ограничениями.
2.5. Построение финитных управлений для линейных систем с
ограничениями
Рассмотрим теперь построение финитных управлений для линейных систем с терминальными ограничениями.
1. Пусть задана система линейных дифференциальных уравнений вида (2.16) с начальными условиями x(t0) = x0. Рассмотрим множество
Q = {x : c**x — b, = 0, i = 17m}. (2.29)
Здесь С — n-мерные вектора, а b — вещественные константы. Требуется построить управление u = u(t), переводящее систему (2.16) из точки ж0 в точку ж1 множества Q, ближайшую к началу координат.
Запишем (2.29) в матричной форме
Q = {ж : С*ж - b = 0}, (2.30)
где С = (с1,... ,cm) — n х m-матрица, b - m-вектор. Пусть матрица С * имеет ранг r. Найдем наилучшее приближенное решение ж (по методу наименьших квадратов) системы:
С*ж - b = 0, (2.31)
то есть
||С*ж - bll = min ||С*ж - b||. (2.32) Это решение определяется [31] по формуле
ж = (С *)+b = (C+)*b, (2.33)
где С + — псевдообратная матрица для матрицы С.
Покажем, что вектор ж имеет наименьшую длину по сравнению с другими векторами, для которых реализуется минимум (2.32). Действительно, обозначим совокупность всех векторов С*ж, ж G Rn, образующих r-мерное подпространство пространства Rm, через S1. Соотношение (2.32) означает, что С*ж представляет собой ортогональную проекцию вектора b на подпространстве S1. В то же время совокупность всех векторов ж G Rn, удовлетворяющих условию С*ж = 0 образует подпространство R2 в Rn. Это подпространство является ортогональным дополнением к r-мерному подпространству R1 из Rn.
Пусть ж (ж = ж) — любой другой вектор из Rn, для которого реализуется минимум (2.32). Тогда
С *ж = С *ж.
Следовательно,
С*(ж - ж) = 0.
То есть
ж - ж G R2.
Но так как
гр I Гр _ Гр Гр - Гр _[_ ( Гр _ Гр\
fA/ _\_ >Aj >AJ >AJ - >AJ | \ >AJ >AJ I ^
то по теореме Пифагора находим:
ЦхЦ2 = уху2 + ух - Х||2 > уху2. (2.34)
Таким образом в силу (2.34) точка х = Х множества Q является ближайшей к началу координат.
Теперь можно найти финитное управление, переводящее систему (2.16) из точки х0 в точку множества Q, ближайшую к началу координат. Подставляя в (2.20) вместо х1 точку х, определяемую соотношением (2.33) получим искомое управление. Таким образом, справедлива
Теорема 2.2. Пусть система (2.16) полностью управляема, а ее решения стеснены дополнительно линейными ограничениями С*х(^1) — Ь = 0. Тогда финитное управление, переводящее систему (2.16) из точки х0 в ближайшую к началу координат точку х множества (2.30), имеет вид:
и(Ь) = N (*)[хо + Ф—1(^1)(у1 — х)],
(2.35)
У1 = / Ф(^1)Ф—1(^)/
¿0
причем столбцами матрицы N(£) являются управления, переводящие систему (2.18) за время — Ь0 из точек е-7 в начало координат (здесь е-7, ] = 1,п, — единичные орты).
Конкретный вид матрицы N(£) может быть определен также, как в предыдущем пункте.
2. Пусть задана управляемая система вида (2.16) с начальными условиями х(Ь0) = х0 и пусть задана некоторая кривая х = £(Ь), Ь £ [Ь0 , ¿1]. Требуется построить финитное управление, доставляющее наименьшее значение величине Цх(Ь1) — £(Ь1)у при условиях сг*х(Ь1) — Ь = 0, % = 1,т, где С — п-мерные вектора, а Ь — вещественные константы. Считаем, что точка £ = £(Ь1) не принадлежит множеству Q.
Как и в предыдущем пункте найдем вначале точку х, в которую нужно перевести систему (2.16). Ясно, что эта точка является ортогональной проекцией точки £, на множество Q. Следовательно, нужно решить следующую задачу: минимизировать функцию
р(х) = ух — £ у (2.36)
при ограничениях (2.30). Пусть
х = у + £. (2.37)
Сделаем замену переменных (2.37) в формулах (2.36) и (2.31). В результате, для нахождения проекции точки £ на множество Q достаточно решить относительно у (по методу наименьших квадратов) систему
С*у + С*£ - Ь = 0. (2.38)
Решение у системы (2.38) имеет вид:
У = -(С *)+(С * £ - Ь).
Используя свойства псевдообратной матрицы (Р*)+ = (Р+)*, (РР+)* = получим окончательно:
у = -СС+£ + С+Ь.
Осуществляя обратную замену переменных, получим
х = £ - СС+£ + (С +)*Ь. (2.39)
В данном случае финитное управление определяется так же, как для задачи пункта 1 по формуле (2.35) только с той разницей, что точка х находится по формуле (2.39).
3. Рассмотрим более общий случай. Пусть задана управляемая система (2.16) с начальными условиями х(£0) = х0. Требуется построить финитное управление, доставляющее наименьшее значение выпуклой дифференцируемой функции конечного состояния д(х(£х)) при линейных ограничениях
С*х(£х) - Ь = 0. (2.40)
Так как система (2.16) полностью управляема, то финитное управление может переводить точку х0 в любое конечное положение хх. Подставляя в (2.20) вместо хх точку х, доставляющую минимум функции д(х) на множестве Q, найдем требуемое финитное управление.
Для нахождения оптимальной точки х может быть применен метод проекции градиента. Идея метода состоит в следующем [43]. Для вычисления очередного приближения к оптимальной точке направление йк = хк - а д'(хк), к = 0,1, 2,..., проектируется на множество Q ив качестве направления спуска берется
) - хк.
Здесь д'(хк) — градиент функции д(х) в точке хк, ак > 0 — шаг спуска, р(й*) — проекция направления йк на множество Q. Так как ограничения (2.40) линейны, метод проекции градиента может быть записан в виде:
хк+х = (Е - СС +)(хк - акд'(хк)), х0 = (С+)*Ь, к = 0,1, 2,... . (2.41)
Таким образом, финитное управление в данном случае определяется по формуле (2.35), причем точка х является предельной точкой последовательности (2.41).
2.6. Численные примеры
Приведем в качестве иллюстрации работы симплексного алгоритма результаты численных расчетов известных по литературе тестовых примеров. Так как симплексный алгоритм ориентирован на линейную систему управления, то приводимая нелинейная система
£ = I(£) + Ви, £(*о) = хо , г £ Т, и £ и,
традиционным способом линеаризуется около опорной траектории:
хр = Ар(£)хр + Ви + ср(г), хр(го) = х0 , г £ Т, и £ и,
Ар(г) = | (£ р-1(г)),
с.(г) = I(£р-1(г)) -1(£р-1(г))£р-1(г),
£р (г) = £(г,и*р (г)),
и*р(г) — управление, доставляющее минимум функционалу
Ыхр(гО).
Предварительно исходная задача (если необходимо) приводится к виду (2.1)-(2.3). Критериями окончания процесса линеаризации являются
1Ы№)) - ^о(хр(^1))| < 4 Ы^-1^)) - Ы№))1 < 4
где , е'2 — заданные заранее величины. Ниже в примерах полагалось ^ = = 0.002. Отметим, что линеаризация около опорной траектории требует достаточно хорошего начального приближения для линейных систем.
Обозначения:
р — номер линеаризованной задачи;
I — номер вспомогательной задачи (2.9) метода модифицированной функции Лагранжа при решении р-й линеаризованной задачи;
— число задач (2.9) при решении р-й линеаризованной задачи;
к — номер итерации;
к
— краевое условие для сопряженной системы;
к
гк — крайняя точка множества достижимости;
к
ук — точка минимума на симплексе а"к;
^(ук) — значение функции ^(ж) в точке минимума ук;
) — значение опорной функции (дк, ¿к). Пример 2.1 [76].
ж 1 = ж2 + и, ж2 = —и, Жз = ж1 /2; ж(0) = (0.5,0, 0)*, = [0,1.5],
и =[-1,1];
ад = жз(.(1.5)),
Г1 = ж1, г2 = ж2 (см. (2.2), (2.3)).
Возьмем в качестве минимизируемой функции (ж) = 2(ж,ж), ж = (ж1 ,ж2 ,ж3)*. Начальное приближение и0(г) = 0, г Е Т.
При реализации метода модифицированной функции Лагранжа использовались следующие величины (см. формулы (2.10), (2.11)):
П = 3п, 7 = 0.5, 6 = 0.1, е' = 10-4, С1 = 10—9.
Результаты счета: р =1, 11 = 4, к = 15, <ж0 = 0.044921; ж = (1.0 х 10—4, 1.0 • 10—4, 3.0)*, ж = (0.0, 0.0, 0.2997)*, ЗД = {а : г Е (¿г-1,^г], г = 1,3}, где
г 0 1 2 3
г 0.0 0.5001 0.9997 1.5
а — 1.0 — 1.0 • 10—4 —0.9996
Ход итераций иллюстрируется таблицей 2.1.
В приведенном примере 1 отрезок [г1, г2] является участком особого управления. Он был получен автоматически. Однако, это скорее "счаст-ливая"случайность, чем закономерный факт. Алгоритмы, использующие построение последовательности крайних точек (построение последовательности экстремальных управлений для линейных задач Майера), могут плохо сходиться. Дело в том, что оптимальное управление в исходной задаче
может оказаться особым. Поэтому при вычислении крайних точек (экстремальных управлений) необходимо учитывать [26] возможность появления (обнаружения) особых участков управления. Если же этого не делать, алгоритм может попросту расходиться из-за неоднозначности определения управления. Такая ситуация уже возникает в линейных задачах с терминальными ограничениями.
Таблица 2.1
р, 1, к дк г к ук к к)
1 1 1 0.5 0.0 0.1875 -0.1164 0.68590000 0.0891 0.3779 0.1359 0.168 -0.0415 -0.0136000
1 1 2 0.2558 0.2718 0.168 0.875 -1.5 0.4687 0.434 0.0454 0.1827 -0.1051 -0.0099
1 1 3 0.368 0.09090000 0.1827 -0.1094000 0.2506 0.1975 0.434 0.0454 0.1827 0.01860000 -0.0099
1 2 4 0.5 0.0 0.0188 -0.1249000 0.5187 0.1461 0.2296 0.2245 0.1696 -0.0597 0.04670000
1 2 5 0.2026 0.2469 0.017 0.875 -1.5 0.4626 0.1149 0.0346 0.222 -0.1852 0.0046
1 2 6 0.07640000 0.0381 0.0222 0.069 -0.3808 0.3477 0.0543 0.00720000 0.2442 -0.0015 0.0019
1 2 7 0.0097 0.0079 0.0244 -0.0776 0.0646 0.25520000 0.0543 0.00720000 0.2442 0.0060 0.0019
1 3 8 0.50540000 7.0 • 10-4 0.0019 -0.125 0.499 0.1816 0.1947 0.2437 0.18460000 -0.0625 0.04960000
1 3 9 0.1971 0.2469 0.00180000 0.875 -1.5 0.4667 0.0998 0.04960000 0.2457 -0.197 0.0067
Продолжение Таблицы 2.1
1 3 10 0.1012 0.05080000 0.0025 0.1222 —0.4944 0.41400000 6.00000000 2.0 • 10—4 0.2977 —0.0117 4.0 • 10—4
1 3 11 0.0011 9.00000000 0.0030 —0.084 0.09530000 0.3004 6.00000000 2.0 • 10—4 0.2977 9.00000000 4.0 • 10—4
1 4 12 0.5006 1.0 • 10—4 2.0 • 10—4 —0.125 0.4998 0.17020000 0.195 0.2439 0.1791 —0.0625 0.0489
1 4 13 0.1953 0.2442 2.0 • 10—4 0.875 — 1.5 0.47940000 0.1 0.05 0.2398 —0.1954000 0.0063
1 4 14 0.10020000 0.05010000 2.0 • 10—4 0.12480000 —0.4995 0.4252 0.0 0.0 0.2977 —0.0124000 0.0
1 4 15 1.0 • 10—4 1.0 • 10—4 3.00000000 —0.0503000 —0.0465 0.3355 0.0 0.0 0.2977 1.0 • 10—4 0.0
Пример 2.2 [26].
х 1 = х2,
¿2 = — й1п(х1) + М, х(0) = (5, 0)*; М1 = [0, 5],
и = [—1,1].
Возьмем в качестве минимизируемой функции ^о(х) =
(хх, х), х —
(¿1, ¿2)*.
Начальное приближение:
( 1, г е [0,1.0],
и0(г) = I -1, г е (1.0,4.5],
[ 1, г е (4.5,5.0].
Результаты счета:
р =1, к = 6, = 3.453019, д = (0.9247, -0.3806)*, ж = (3.1907, -1.3201)*, ЗД = {аг : г е (гг-1,гг], г = 175} где:
г 0 1 2 3 4 5
гг 0.0 0.9519 0.9698 4.5610 4.5845 5.0000
аг — 1.0000 -0.6079 -1.0000 0.6079 1.0000
Ход итераций иллюстрируется таблицей 2.2: Таблица 2.2
Р,к дк хк ук Мдк )
1 1 0.9305 -0.3662 3.16220000 -1.3847 3.16220000 -1.3847 3.4496 3.4521
1 2 0.916 -0.4011 3.2485 -1.181 3.1783 -1.3467 3.44940000 3.4518
1 3 0.92080000 -0.3901 3.21920000 -1.248 3.186 -1.3276000 3.451 3.45150000
1 4 0.9231 -0.3847000 3.20530000 -1.2805 3.1897 -1.3182 3.4512 3.45130000
1 5 0.9242 -0.3819 3.19860000 -1.2963 3.1915 -1.3136 3.4512 3.4512
1 6 0.92470000 -0.3806 3.19540000 -1.3039 3.1915 -1.3136 3.4512 3.4512
3 Решение нелинейной задачи оптимального управления
Для решения краевой задачи, возникающей в результате применения принципа максимума Л. С. Понтрягина, предлагается алгоритм, являющийся комбинацией метода неподвижной точки и квазиньютоновского метода. Первый из них обладает широкой областью сходимости и применим для достаточно широкого класса функций, а второй имеет локальную сверхлинейную сходимость. Приводятся результаты счета. Результаты данного раздела опубликованы в работах [14], [15].
Известно, что принцип максимума Л. С. Понтрягина [68] позволяет решать достаточно широкий класс задач теории оптимальных процессов. Однако, численная реализация такого подхода может быть затруднена отсутствием эффективного способа нахождения начальных условий для решения сопряженной системы, с помощью которой формулируются необходимые условия оптимальности. Отметим следующие два направления для построения приближенных методов решения задач оптимального управления:
1) одним из них является сведение решения задачи оптимального управления к решению краевой задачи (см., например, [61], [66], [68], [70],
[84]);
2) другое направление составляют методы последовательных приближений, связанных с построением минимизирующей последовательности управлений (в функциональных пространствах) [26], [41], [53], [58], [61], [66], [73], [76], [84].
Предлагаемый метод примыкает к первой группе. Решение краевых задач, возникающих в результате применения принципа максимума, сводится либо к решению системы нелинейных алгебраических уравнений, либо к максимизации вогнутой функции в конечномерном пространстве. Для решения алгебраических уравнений в этом случае можно воспользоваться, например, методами прогонки, методом Ньютона или квазиньютоновскими методами [61], [66], [84], а в случае максимизации функции — методами градиентного типа (проекции градиента, условного градиента) [28], [29], [41]. Однако упомянутые выше методы в некоторых случаях либо сходятся плохо, либо вообще не сходятся. Квазиньютоновские методы могут оказаться предпочтительнее в силу их значительно более быстрой сходимости. При этом они требуют вычислений функции не больше, чем градиентные методы. На практике же проблема заключается в том, что обычно неизвестно
достигнута ли уже область сходимости (например, методом Ньютона) или нет. Даже более того, неизвестно, будет ли оптимальная точка иметь такую область.
Отметим особо алгоритмы [33], [41], которые получают аппроксимацию решений использованием симплексов. Эти алгоритмы позволяют находить оптимальные управления и в том случае, когда множество достижимости не строго выпукло (например, имеет форму многогранника). В этих случаях методы градиентного типа сходятся очень медленно, а применение симплексных алгоритмов может существенно ускорять итерационный процесс и фактически обеспечивать вычислительную устойчивость. Это объясняется тем, что в вычислительных целях непрерывная система управления раздискречивается. Поэтому алгоритмы применяются фактически к дискретным системам управления, а для дискретных систем управления множества достижимости могут оказаться многогранниками (как это имеет место, например, в случае задачи минимизации топлива).
В данном разделе предлагается [14], [15] для решения уравнений, возникающих в результате применения принципа максимума, использовать комбинированный метод: метод неподвижной точки для непрерывных отображений [79] и квазиньютоновский метод [19], [99], [103]. Параметрический метод приближенного вычисления неподвижных точек непрерывных отображений обладает достаточно широкой областью сходимости ввиду того, что использует для их нахождения кусочно-линейную аппроксимацию функции на триангуляциях. При этом за "глобальную"сходимость метода приходится расплачиваться большим количеством вычислений значений функции. Поэтому, оборвав этот процесс с некоторой точностью, можно продолжить решение квазиньютоновским методом, который имеет при хорошем начальном приближении локальную сверхлинейную сходимость. Таким образом, можно построить алгоритм, который обладает как свойством "глобальной"сходимости метода неподвижной точки (симплексного метода), так и свойством локальной сверхлинейной сходимости квазиньютоновских методов. Главное, однако, состоит в том, что свойство "глобаль-ной"сходимости симплексных алгоритмов остается в силе в любом случае, даже без какого-либо предположения о гладкости функции. Симплексный метод часто работает лучше других методов, когда начальное приближение плохое или когда отображение имеет некоторые патологические свойства.
3.1. Постановка задачи
Рассмотрим следующую задачу оптимального управления: на траекториях управляемой системы
ж = /(ж и г Е Т = [г0, г1] , (3 1)
ж(г0) = ж0, и(г) е и (
минимизировать функционал
Ф0(и) = ^0(ж(г1)) (3.2)
при дополнительных условиях
Фг(и) = ((ж(г1)) = 0, г = 1,?. (3.3)
Предположим, что функции (¿(ж), г = 0,?, дифференцируемы на Яр, вектор-функция ](ж, и, г) и матрица частных производных /Ж(ж,и,г) непрерывна по своим аргументам на Яр х и х Т, множество и С Яг компактно. Допустимыми управлениями будем считать кусочно-непрерывные на Т вектор-функции и(г) со значениями в и.
Пусть и(г) — допустимое управление, а ж(г) — соответствующая ему траектория. Положим Н(ж,^,и, г) = (ж,и,г), где функция ^(г) удовлетворяет дифференциальному уравнению
^ = —/х*(ж,и,г)^ (3.4)
с краевым условием
^(г1) = —(0х(ж(г1)) — ^ Лг(гх(ж(г1)) (3.5)
¿=1
для некоторого параметра Л = (Л1,... , ЛЧ).
Говорят, что управление и(г) удовлетворяет принципу максимума, если найдется вектор Л = 0 такой, что выполняется условие максимума
н (ж(г),^(г),м(г),г) = тах н (ж(г),^(г),и,г). (3.6)
иеи
Пусть и(г) — оптимальное управление, а ж(г) — соответствующая оптимальная траектория. В этом случае принцип максимума [68] утверждает, что найдется непрерывная вектор-функция ^(г) ф 0 такая, что ^(г) является решением уравнения (3.4) с краевым условием (3.5) для некоторого
Л = Л и оптимального процесса {&(£),£(£)}, причем выполняется условие максимума (3.6), т.е.
Н(£(£),ф(г),й(г),£) = тахН(¿(¿), ф(£), и, ¿).
Пусть и(£) можно однозначно определить из условия максимума (3.6) как функцию ж(£),ф(£) : и(£) = 7(ж(£),ф(£)). Подставляя управление и(£) в системы (3.1) и (3.4), получим систему 2р уравнений
Xх = / (х,7 ф = -/Лх,7 (х,ф),0ф (3.7)
решения которой определяются заданием ж(£о), ф(£о) и Л. Поэтому, решая систему (3.7) с краевым условием
ж(£о) = х0, ф(^) = )) - ^Х(х(^))Л, ( = ((1 ,...,(), (3.8)
получим решение исходной задачи оптимального управления. Таким образом, получена краевая задача для системы обыкновенных дифференциальных уравнений с достаточно сложной структурой правой части, множество решений которой содержит решение исходной задачи.
Известно что всякий метод решения нелинейной краевой задачи сводит, по существу, ее решение к решению некоторой нелинейной системы уравнений относительно некоторых параметров. Далее эта система решается каким-либо методом решения нелинейных систем. Выбор вспомогательных переменных (параметризация задачи) часто определяет успешность ее решения. Приведем два способа определения вектора неизвестных параметров.
1) Положим £ = ф(£о). В качестве неизвестных переменных примем (£, Л). Интегрируя систему (3.7) "слева направо"от ¿о до ¿1 при начальном условии (хо,£), определим , £) и , £), а тем самым будет установлена функциональная зависимость (р + д)-мерного вектора г = (ф (¿1 , £),
,£))) от (р + д)-мерного вектора д = (£,Л). Таким образом, в этом случае решение краевой задачи (3.7),(3.8) сводится к решению системы нелинейных уравнений
ф(¿1,£) + ,£)) + ,£))Л = 0, р(ж(*1 ,£)) = 0 (3.9)
относительно неизвестных переменных £, Л.
2) Положим £ = ж(^1). В качестве неизвестных переменных примем (£, Л). Найдем ф(^1) из условия трансверсальности (3.8). Интегрируя систему (3.7) "справа налево"в обратном времени, при начальном условии
(£, , £, Л)), определим х(£0 , £, Л). В результате, будет установлена функциональная зависимость (р + q)-мерного вектора г = (х(£0 , £,Л),^(£)) от (р + q)-мерного вектора д = (£,Л). В этом случае система нелинейных уравнений имеет вид
относительно неизвестных переменных £, Л. Таким образом, формально решение краевой задачи сводится к решению системы р + q нелинейных уравнений г(д) = 0. Отметим здесь лишь, что иногда предпочтительнее выбирать второй способ параметризации задачи: в этом случае решение системы нелинейных уравнений (3.10) может оказаться более устойчивым, чем решение системы (3.9).
Началом конструктивной аппроксимации неподвижных точек непрерывных отображений послужил алгоритм Скарфа [115]. Алгоритмы, которые вырабатывают более точную аппроксимацию, основывались на гомо-топиях [100], [101]. Использовать гомотопии в симплексных алгоритмах для получения аппроксимации неподвижных точек было предложено Ивзом. Симплексные алгоритмы, использующие векторную разметку, включают в себя кусочно-линейную аппроксимацию функции на триангуляциях [101], [102], [93], [94]. Более подробные и общие обзоры работ по этой тематике можно найти, например, в [79], [93], [94], [107]. Там же приведена обширная библиография по данным вопросам.
Рассмотрим систему нелинейных алгебраических уравнений
где г : Яп ^ Яп — непрерывная функция. Будем предполагать, что существует решение д системы (3.11).
Одним из способов решения уравнения (3.11) является сведение решения исходной задачи к решению более простой задачи, вводя в рассмотрение целое семейство отображений /(т, д), зависящих от некоторого параметра т. При этом для некоторого конкретного значения т (например, при т = £0) получается система гх(д) = 0, имеющая известное решение д0, а при т = 0 — исходное уравнение г(д) = 0, решение которого необходимо найти. Такой переход можно осуществить, например, с помощью линейной
х0 - х(*0 ,£,Л) = 0, <р(£) = 0
(3.10)
3.2. Метод неподвижной точки
г (д) = 0,
(3.11)
гомотопии [64], стр. 226]:
1(т,д) = 0(д - до) + (1 - 0)*(д), (3.12)
где 0 = т/^о, 0 < 0 < 1.
Решение до уравнения 1(^о,д) = 0 известно, а уравнение 1(0, д) = 0 нужно решить. Рассмотрим уравнение
1(т,д) = 0, т £ [0,^о]. (3.13)
Предположим, что уравнение (3.13) для каждого т £ [0,£о] имеет решение д = д(т), непрерывно зависящее от т. Геометрически это означает что точка д описывает пространственную кривую
С = {д £ Яп : д = д(т), т £ [0Л]} ,
одним концом которой служит некоторая заданная точка до, а другим — точка, являющаяся решением д = д(0) уравнения 1(0, д) = *(д) = 0. Поэтому необходимо уметь приближенно отслеживать кривую С.
Гомотопный путь С можно отслеживать либо с помощью численного решения определенной задачи Коши (алгоритмы продолжения), которой удовлетворяет кривая С, либо с помощью кусочно-линейной аппроксимации (симплексами) вспомогательного отображения 1(т, д) (симплексные алгоритмы). Подробные обзоры работ по методам продолжения и симплексным алгоритмам можно найти, например, в [79], [93], [94], [107].
Пусть у = (т, д), р = [уо, у1,..., уп] есть п-мерный симплекс (п-симплекс), а в (у) = Вд + с — аффинное отображение, совпадающее со значениями функции 1(у) в вершинах симплекса р, то есть
Вдг + с = 1(уг) , г = 0, п .
Пусть Ь = пч . Матрица Ь называется матрицей меток, а
Здесь п х п-матрица В и вектор с определяются из последних равенств. 1 1 ... 1
1(уо) ...... 1(уп)
сама функция 1(у) — векторной разметкой (векторным помечиванием).
Легко можно показать, что если точки 1(уо),...,1(уп) аффинно независимы (находятся в общем положении), то нуль аффинной функции в (у)
п
может быть представлен в виде [64], стр. 190-192] д = ^ а^дг, причем
__г=о
коэффициенты а^, г = 0,п определяются в результате решения системы
Ьа = ео . (3.14)
Здесь е0 = (1,0,..., 0)* — единичный орт. Так как предполагается, что кривая С пересекает п-мерную грань (п + 1)-симплекса, то ноль должен принадлежать п-симплексу, и в этом случае будет а > 0. Поэтому говорят, что п-симплекс полностью размечен в соответствии с отображением /, если существует решение а > 0 системы (3.14). Таким образом, можно сказать, что симплексные алгоритмы являются, по существу, методами приближенного отслеживания кривой С посредством (п + 1)-мерных симплексов, определенных на триангуляциях и имеющих две полностью раз-меченне п-мерные грани. Одна полностью размеченная грань считается входом для кривой С, а другая ее выходом. Так как обе грани полностью размечены, то соответствующие им матрицы меток невырождены. Поэтому переход от входной грани к выходной осуществляется (как в линейном программировании) введением в базис на некоторое место нового вектор-столбца. В результате, получаем базис, соответствующий выходной грани. Симплексный алгоритм вырабатывает конечную или бесконечную последовательность полностью размеченных п-мерных граней и соответствующих (п + 1)-симплексов.
Ниже описывается алгоритм нахождения решения системы (3.11), основанный на триангуляциях [79], стр. 84-95] с непрерывным уменьшением размера сетки (уточняемые триангуляции). Эти уточняемые триангуляции становятся все более мелкими, когда гомотопный параметр т приближается к нулю. Описание указанных триангуляций приводится в работах [79],
[94].
Обозначим: а"к — (п + 1)-мерный симплекс на к-й итерации; Рк—1, Рк — его грани, причем уг = (тг,дг), г = 0,п, — вершины входной грани Рк—1, а уг = (т\ дг), г = 0,п, — вершины выходной грани рк, У = (у0,..., уп),
G = (д0,... , <Т ^ Ррк = * ... \ — матрица меток грани рк, Я =
1(у0) ... 1(уп)
(д0,...,дп) = Ь—, та = 2—^0, 5 = 0,1,....
Алгоритм 3.1 (гомотопный алгоритм):
Шаг 0. Задать начальную точку д0 Е Лп. Задать числа £0, £1, е2, (. Положить к = 0, т = 0.
Шаг 1. Начальный шаг: Положить т = т + 1;
а
) определить вектор п Е Яп+1: пГ = пГ = зГ 1, г = 1,п;
b) определить точку у-1, являющуюся опорной (центральной) при нахождении остальных точек у0,... ,уп симплекса а"к размерности п + 1, в соответствии с точкой пт;
c) вычислить остальные вершины у0,...,уп симплекса а"к размерности п + 1, т.е. вычислить вершины грани рк = [у0,..., уп];
d) вычислить матрицу разметки входной грани рк:
г ... 1
Рк /(у0) ... /(?/") '
Найти обратную матрицу 1
e) принять у-1 в качестве входного вектора у+, т.е. вектора который нужно ввести в базис к.
Шаг 2. Шаг линейного программирования: положить к = к + 1;
a) вычислить /(у+);
b) найти номер вектора, который должен быть замещен в базисе :
у- = у^О ;
c) удалить столбец (1, /(у-))* из базиса ЬРк-1 и ввести столбец (1, /(у+))*. Получаем новый базис к. В результате получаем также матрицу к1.
Шаг 3. Определение нового (п + 1)-симплекса а"к:
a) вычислить вершины у-1, у0,... , уп нового (п + 1)-симплекса ^к;
b) установить соответствие между вершинами симплексов о"к-1 и ^к. Найти точку у+, противолежащую грани рк.
Шаг 4. Переход на новый уровень: если (у+)к > (у+)к-1, то перейти на шаг 2, иначе — к следующему шагу.
Шаг 5. Вычисление очередной неподвижной точки:
a) вычислить неподвижную точку, лежащую в грани рк : ут = дт =
b) определить г(дт) и /(ут);
c) если ||/(ут)|| > £2, то перейти на шаг 2, иначе — к следующему шагу 6.
Шаг 6. Проверка точности решения: если ||г(дт)|| < £1 — остановиться, иначе — перейти к шагу 7.
Шаг 7. Рекурсивный цикл:
а) если
к (дт )||-||г (дт-1)
< ( и £2 > 10 8, положить £2 = 0.1 £2;
b) вычислить дто = max |gm — gm если дто > 0.5£m, то положить
i=1,n
^m+i = 0.5£m, иначе £m+1 = дт. Перейти к шагу 1.
В дальнейшем вычисления, требуемые для получения очередного приближения к решению уравнения (3.11) (очередной неподвижной точки) будут называться "этапом", а вычисления, необходимые для интегрирования прямой и сопряженной системы (3.7), (3.8), — "итерацией". Если на шаге 5 гомотопного алгоритма необходимая точность вычисления решения уравнения (3.11) не достигнута, то для ускорения итерационного процесса на этом шаге необходимо проделать дополнительно ньютоновские шаги. Подробное описание квазиньютоновского алгоритма дается в следующем пункте.
3.3. Квазиньютоновский метод
Рассмотрим систему нелинейных алгебраических уравнений
z(g) = 0 , (3.15)
где z : Rn ^ Rn, g e Rn.
Введем обозначения. Пусть zk = z (gk), Agk = gk — gk—1, Azk = zk — zk—1. Будем обозначать через z'(g) матрицу Якоби, вычисленную в точке g.
Если матрица Якоби z'(g) легко вычисляется, наиболее известным методом для решения системы (3.15) является метод Ньютона
gk+1 = gk — z'(gk)—1z(gk), k = 0,1.... (3.16)
Если же z'(g) недоступна для вычисления, то матрица z'(gk) (или z'(gk)—1) заменяется некоторой ее аппроксимацией.
Ниже будет рассматриваться квазиньютоновский метод в следующем виде
,k+1_ Jfe rwJ^ П. . _ П. ■ (Agk — DkAzk)w*k
w*k Azk
/+1 = sk — Dk z (gk), Dk+1 = D + (Agk .. DA* , k = 0,1....
Предложенный в этом пункте алгоритм является незначительной модификацией алгоритмов [103], [19]. Предположим, что на к-й итерации имеется множество линейно независимых векторов Д^ = {Д^, ] Е }, где = {^'и, ... — множество упорядоченных по убыванию целых чисел. Предположим так же, что вектор Дг^ может быть
заменен в матрице Д^к на вектор Д^к, так что в результате такой замены получится новый базис.
Обозначим
4 = 4\{Ы , Ак = {Д^ ,; е Jk}.
Тогда в качестве и>к для коррекции Д можно выбрать вектор, ортогональный ко всем Д^- , ] е 4к. Очевидно, что в этом случае вектор и>к имеет вид
= Д^к — Д5к ,
где Д^к — ортогональная проекция вектора Д^к на подпространство, натянутое на вектора Д^-, ] е 4. Пусть А+ — матрица, псевдообратная [31], стр. 32] для матрицы Ак. Легко можно показать, что
Д^к = АкА+Д^к .
Алгоритм 3.2 (квазиньютоновский алгоритм):
Шаг 0. Задать начальную точку д0, матрицу Д0, число £1. Вычислить г0. Положить к = 0.
Шаг 1. Составить множество 4 = {к—1,... , к—п}. Положить Д^к = Е (Е — единичная матрица).
Шаг 2. Если || гк || < £1, то остановиться.
Шаг 3. Найти Ддк = —Дкгк. Положить дк+1 = дк + Ддк. Вычислить
гк+1 и Дгк = гк+1 — гк.
Шаг 4. Определить номер ^ вектора, который необходимо исключить из базиса Д^к. Определить множество = 3к\{^кв}.
Шаг 5. Составить матрицу Ак = {Д^-, ] е 3'к}. Найти псевдообратную матрицу А+. Вычислить Д5к = А кА+Дгк, а также и>к = Дгк — Д5к.
Шаг 6. Вычислить матрицу Д+1 по формуле
Дк+1 = Д + (Ддк — Д Дгк)
и>*к Дгк
Шаг 7. Составить матрицу Д^к+1 = {Д^к, Д^-, ^ е 3'к} и соответствующее этой матрице множество 3к+1 = {к, 3'}.
Шаг 8. Положить к = к + 1 и перейти к шагу 2.
Алгоритм 3.2 сходится к д [19], если точка д0 достаточно близка к д, а норма ||Д0 — г'(д0)—1| достаточно мала.
На шаге 5 алгоритма 3.2 для вычисления псевдообратной матрицы и на шаге 4 при проверке на линейную зависимость векторов можно воспользоваться методом Гревилля [31], стр. 39-40] последовательного нахождения псевдообратной матрицы.
3.4. Общая схема алгоритма
В предыдущих пунктах были приведены алгоритмы решения нелинейных алгебраических уравнений: рестартовый гомотопный (раздел 2) и квазиньютоновский (раздел 3). В данном пункте приводится общая схема использования этих алгоритмов. Оба метода для получения приближения к решению уравнений вырабатывают последовательность смежных симплексов (два симплекса называются смежными, если они имеют общую грань). При этом гомотопный метод может быть применен для убывающей последовательности сеточных размеров, если использовать в качестве начальной точки на любом этапе аппроксимацию, полученную на предыдущем этапе. На тесную связь между этими двумя методами впервые, вероятно, указал Сейгал [113]. Он предложил чередовать шаги симплексного алгоритма с ньютоновскими шагами. Техника чередования PL-шагов (начальные буквы слов piecewise linear) и ньютоновских шагов для ускорения гомотопных алгоритмов обсуждалась в [114], [117].
В работе [15] предлагается более простой способ — комбинирование гомотопного и квазиньютоновского алгоритмов, а не чередование PL-шагов и квазиньютоновских шагов для глобализации метода Ньютона. А именно, после нахождения очередного приближения к решению системы уравнений гомотопным алгоритмом проделываются пробные ньютоновские шаги. Если они не успешны (нарушается монотонное убывание невязки отклонения левой части системы от нуля), вновь перезапускается гомотопный алгоритм, причем в качестве начальной точки принимается последнее приближенное решение системы уравнений.
Пусть задана начальная точка g0 е Rn и начальный сеточный размер ¿0. Предположим, что используется и перезапускается гомотопный алгоритм. Пусть на m-м этапе имеется приближенный нуль gm-1 функции z и сеточный размер ¿m-1. Для получения очередного приближения к решению уравнений проделываются вычисления в соответствии с шагами 1-5 гомотопного алгоритма. Предположим, что после многократного выполнения шагов 1-5 алгоритма 3.1 достигнута необходимая точность вычисления решения уравнения (3.13). В результате, будут найдены (см. шаг 5 алгоритма 3.1) gm, z(gm), а еще — Gm и Zm, соответствующие матрицам Ym и LPm. Пусть N(k,g) означает шаги 1-8 алгоритма 3.2 при gk = g. Тогда, используя матрицы Gm и Zm, определяем в качестве начального приближения Dk = AGmAZm1 и проделываем шаги квазиньютоновского алгоритма N(k, gm). Если алгоритм 3.2 не находит решения с необходимой точностью,
сеточный размер уменьшается и перезапускается вновь гомотопный алгоритм.
Алгоритм 3.3 (общий алгоритм).
Шаг 0. Задать начальную точку д0 е Лп. Задать числа ^0,£1,£2,^. Положить к = 0, т = 0.
Шаг 1. Проделывать шаги 1-5 гомотопного алгоритма (алгоритма 3.1) до тех пор, пока не будет достигнута необходимая точность приближения е2 для гомотопного алгоритма.
Шаг 2. Проверка точности решения: если ||г(дт)|| < £1 — остановиться, иначе — перейти к шагу 3.
Шаг 3. Проделывать шаги N (к, дт) квазиньютоновского алгоритма.
Шаг 4. Рекурсивный цикл:
а) если
|z (gm)||-||z (gm-1)
< Z и e2 > 10 8, положить e2 = 0.1e2;
^Г - g
m—1
; если > 0.5 £m, то положить
б) вычислить = max
i=1,n
¿m+1 = 0.5 im, иначе ¿m+1 = Mm• Перейти к шагу 1.
Легко видеть, что алгоритм 3.3 генерирует последовательность {gm} приближенных решений уравнения z(g) = 0. Последовательность {gm} ограничена и, следовательно, имеет по меньшей мере одну точку сгущения. Если же предположить, что z(g) имеет изолированные нули, то полученная алгоритмом 3.3 последовательность {gm} сходится при m ^ то к g.
3.5. Численные примеры
Для проверки эффективности предложенного алгоритма 3.3 были составлены программы и проведены расчеты задач, известных в литературе. Приведем здесь итоговые результаты проверки, которые демонстрируют типичное поведение алгоритма.
Примеры решались при следующих данных: точность решения задачи £1 = 10-4, точность решения задачи по методу неподвижной точки е2 =
0.01, Z = 0.1.
Обозначения: m — номер этапа, k — номер итерации, тmk — значения гомотопного параметра, gmk — краевое условие для сопряженной системы
(см. (3.7) — (3.10)), /mk = /(тmk, gmk) — значение функции разметки в методе неподвижной точки (см. (3.12)).
Пример 3.1 [26]:
¿1 = ¿2 , ¿2 = — sin x + u,
x(0) = (5,0)*, T = [0, 5], U = [—1,1], Фо = ¿2 + ¿2 .
Ограничения отсутствуют. Использовался второй вид параметризации (см. пункт 3.1). Начальные приближения: = 0.2, g0 = x(t1) = (3.1, —1)*. Ход итераций иллюстрируется таблицей 1.
Результаты счета: Ф0 = 11.90805, g = (3.178936, —1.342526)*. Оптимальное управление — релейное управление с двумя точками переключения t1, t2:
где t1 = 0.982443, t2 = 4.550369.
Пример 3.2 [76]:
¿1 = x2, ¿2 = — sin x1 + U, ¿3 = u,
¿(0) = (1,0, 0)*, T = [0, 2n], U = [0,1], Фо = ¿3, Ф1 = ¿1, Ф2 = ¿2.
Начальные приближения: = 0.2, g0 = ^(0) = (0,0, —1)*. Ход итераций иллюстрируется таблицей 2.
Результаты счета: Ф0 = 0.999637, g = (—0.997403,0, —1)*. Оптимальное управление — релейное управление с двумя точками переключения
где = 1.175244, £2 = 2.174881. Ниже в таблицах символом "*" помечаются итерации, получаемые квазиньютоновским алгоритмом.
u(t) =
Таблица 3.1. Результаты счета алгоритма 3.3.
т, к,ттк „тк У *(„ тк) Фо, II*(„тк)11, || /тк ||
1 9 0.1 3.2608 -1.29613 -0.164538 0.201425 12.3128 0.260086 0.00162065
1 * 10 3.11629 -1.37099 0.0918759 -0.143661 11.5909 0.170528
1 * 11 3.18719 -1.34041 -0.00215886 0.00999895 11.9549 0.0102294
1 * 12 3.18129 -1.34186 -0.000998527 0.00316346 11.9212 0.00331731
1 * 13 3.17897 -1.34252 -1.56394 • 10-5 5.02865 • 10-5 11.9082 5.26623 • 10-5
Таблица 3.2.
т, к,ттк „тк У г („тк) Фо, 1|г(„тк)||,
|| /тк ||
1 —0.9006 0.9264 0
37 —0.3151 0.3476 0.9895
0.1 — 1 0 0.01036
1 —0.9421 0.2905 0.6859
59 —0.06369 0.002756 0.2905
0.05 — 1 1.11 • 10—16 0.005619
1 —0.9608 0.05106 0.936
* 60 0.09059 0.07051 0.08706
— — 1 —2.22 • 10—16 —
2 —0.9805 0.01564 0.9778
69 0.05181 0.03983 0.04279
0.05 — 1 4.441 • 10—16 0.001046
2 —0.997 0.001503 0.9979
* 70 —0.0003432 —8.748 • 10—5 0.001506
— — 1 0 —
2 —0.9974 —0.0001966 0.9998
* 71 7.389 • 10—5 3.307 • 10—5 0.0001994
— — 1 0 —
2 —0.9974 — 1.005 • 10—6 0.9996
* 72 —3.527 • 10—7 —3.749 • 10—7 1.073 • 10—6
— — 1 0 —
Литература
[1] Аввакумов С. Н., Киселев Ю. Н., Орлов М. В. Методы решения задач оптимального управления на основе принципа максимума Понтряги-на // Труды Математического института им. В. А. Стеклова РАН. 1995. Т. 211. С. 3-31.
[2] Александров В. М. Квазиоптимальные процессы в автоматических системах. Известия АН СССР. Техническая кибернетика. Т. 13, № 3. 1975.
[3] Александров В. М. Приближенное решение задачи линейного быстродействия // Автоматика и телемеханика. 1998. № 12. С. 3-13.
[4] Александров В. М. Численный метод решения задачи линейного быстродействия // Журнал вычислительной математики и математической физики. 1998. Т. 38. № 6. С. 918-931.
[5] Антоник В. Г., Срочко В. А. Решение задач оптимального управления на основе методов линеаризации // Журнал вычислительной математики и математической физики. 1992. Т. 32. № 7. С. 979-991.
[6] Атанс М., Фалб П. Оптимальное управление. М.: Машиностроение, 1968.
[7] Баничук Н. В., Петров В. М., Черноусько Ф. Л. Численное решение вариационных и краевых задач методом локальных вариаций // Журнал вычислительной математики и математической физики. 1966. Т. 6. № 6. С. 947-961.
[8] Батурин В. А., Урбанович Д. Е. Методы улучшения второго порядка для задач оптимального управления // Известия РАН. Теория и системы управления. 1997. № 3. С. 99-103.
[9] Башков Е. А. Алгоритм решения задачи быстродействия по амплитуде и "мере" управления // Теория оптимальных процессов. Киев: Издательство Института кибернетики. 1974. С. 15-23.
[10] Белман Р. Динамическое программирование. М.: Издательство иностранная литература, 1960.
[11] Болдырев В. И. Построение финитных управлений для линейных систем с ограничениями // Новосибирск, 1992. 17 с. (Препринт / РАН. Сибирское отделение. Институт математики; № 8).
[12] Болдырев В. И. Минимизация псевдовыпуклого функционала на множестве конечных состояний линейной системы управления // Известия РАН. Теория и системы управления. 2002. № 3. С. 49-52 (V. I. Boldyrev. Minimization of a Pseudo-Convex Functional on a Set of Terminal States of a Linear Control System // Journal of Computer and Systems Sciences International. 2002. V. 41(3). P. 382-388).
[13] Болдырев В. И. Решение линейной задачи оптимального управления с терминальными ограничениями // Известия РАН. Теория и системы управления. 2003. № 5. С. 41-50 (V. I. Boldyrev. Solution of a Linear Optimal Control Problem With Terminal Constraints // Journal of Computer and Systems Sciences International. 2003. V. 42(3)).
[14] Болдырев В. И. Численное решение задачи линейного быстродействия // Фундаментальная и прикладная математика. 1999. Т. 5. № 3. С. 637-648.
[15] Болдырев В. И. Численное решение задач оптимального управления // Известия РАН. Теория и системы управления. 2000. № 3. С. 85-92. (англ. перевод: Boldyrev, V. I. Numerical solution of optimal control problems // Journal of Computer and Systems Sciences International. 2000. V. 39. № 3. P. 415-422.)
[16] Болтянский В. Г. Математические методы оптимального управления. М.: Наука, 1969.
[17] Брайсон А., Денхем В. Применение наискорейшего спуска к задачам оптимального управления // Ракетная техника и космонавтика. 1964. № 2.
[18] Брайсон А., Хо Ю-Ши. Прикладная теория оптимального управления. М.: Мир, 1972.
[19] Бурдаков О. П. Устойчивые варианты метода секущих для решения систем уравнений // Журнал вычислительной математики и математической физики. 1983. Т. 23. № 5. С. 1027-1040.
[20] Васильев О. В., Терлецкий В. А. Оптимальное управление краевой задачей // Труды Математического института им. В.А. Стеклова РАН. 1995. Т. 211. С. 121-130.
[21] Васильев О. В. Методы оптимизации в функциональных пространствах. Иркутск: Издательство Иркутского университета, 1979.
[22] Васильев О. В. К вопросу о численном решении задачи терминального управления // Труды Иркутского университета. Иркутск: Издательство Иркутского университета, 1968.
[23] Васильев Ф. П. Лекции по методам решения экстремальных задач. М.: Издательство Московского университета, 1974.
[24] Васильев О. В., Бельтюков Н. Б., Терлецкий В. А. Алгоритмы оптимизации динамических систем, основанные на принципе максимума // Вопросы кибернетики. М.: Наука, 1991. С. 17-38.
[25] Васильев О. В., Тятюшкин А. И. Опыт в решении задач оптимального управления на основе необходимых условий для оптимальности типа принципа максимума // Вопросы устойчивости и оптимизация динамических систем. Иркутск: Издательство Иркутского университета, 1983. С. 43-64.
[26] Васильев О. В., Тятюшкин А. И. Об одном методе решения задач оптимального управления основанного на принципе максимума // Журнал вычислительной математики и математической физики. 1981. Т. 21. № 6. С. 1376-1384.
[27] Васильев О. В., Тятюшкин А. И. К численному решению задач линейного быстродействия // Дифференциальные и интегральные уравнения. Иркутск: Издательство Иркутского университета. 1973. Вып. 2. С. 57-69.
[28] Габасов Р., Кириллова Ф. М. Построение последовательных приближений для некоторых задач оптимального управления // Автоматика и телемеханика. 1964. Т. 27. № 2. С. 5-17.
[29] Габасов Р., Кириллова Ф. М. Качественная теория оптимальных процессов. М.: Наука, 1971.
[30] Габасов Р., Кириллова Ф. М. Оптимизация линейных систем. Минск: Издательство Беларусского государственного университета, 1973. Т. 33. № 9. С. 31-62.
[31] Гантмахер Р. Теория матриц. М.: Наука, 1966.
[32] Гиндес В. Б. К задаче минимизации выпуклого функционала на множестве конечных состояний линейной системы управления // Журнал вычислительной математики и математической физики. 1966. Т. 6. № 6. С. 962-970.
[33] Гиндес В. Б. Один метод последовательных приближений для решения линейных задач оптимального управления // Журнал вычислительной математики и математической физики. 1970. Т. 10. № 1. С. 216-223.
[34] Грачев Н. И., Евтушенко Ю. Г. Библиотека программ для решения задач оптимального управления // Журнал вычислительной математики и математической физики. 1979. Т. 19. № 2. С. 367-387.
[35] Гроссман К., Каплан А. А. Нелинейное программирование на основе безусловной минимизации. Новосибирск: Наука, 1981.
[36] Демьянов В. Ф. К построению оптимальной программы в линейной системе // Автоматика и телемеханика. 1964. Т. 25. № 1. С. 3-11.
[37] Демьянов В. Ф., Рубинов А. М. Приближенные методы решения экстремальных задач. Л.: Издательство Ленинградского государственного университета, 1968.
[38] Евтушенко Ю. Г. Методы решения экстремальных задач и их применение в системах оптимизации. М.: Наука, 1982.
[39] Ермольев Ю. М., Гуленко В. П. О численных методах решения задач оптимального управления // Кибернетика. 1966. № 1. С. 72-78.
[40] Ермольев Ю. М., Гуленко В. П., Царенко Т. И. Конечно-разностный метод в задачах оптимального управления. Киев: Наукова думка, 1978.
[41] Жолудев А. И., Тятюшкин А. И., Эринчек Н. М. Численные методы оптимизации управляемых процессов // Известия АН СССР. Техническая кибернетика. 1989. № 4. С. 14-31.
[42] Калман Р. Е. Об общей теории управления I // Международный конгресс ИФАК. М.: Физматгиз, 1960.
[43] Карманов В. Г. Математическое программирование. М.: Наука, 1975.
[44] Квакернаак Х., Сиван Р. Линейные оптимальные системы управления. М.: Мир, 1977.
[45] Киселев Ю. Н. Оптимальное управление. М.: Издательство Московского государственного университета, 1988.
[46] Киселев Ю. Н. Быстро сходящиеся алгоритмы для линейного оптимального быстродействия // Кибернетика. 1990. Т. 62. № 6. С. 47-57.
[47] Киселев Ю. Н., Орлов М. В. Численные алгоритмы линейных быстродействий // Журнал вычислительной математики и математической физики. 1991. Т. 31. № 12. С. 1763-1771.
[48] Киселев Ю. Н. Построение точных решений для нелинейной задачи оптимального быстродействия специального вида // Фундаментальная и прикладная математика. 1997. Т. 3. Вып. 3. С. 847-868.
[49] Кирин Н. Е. К решению общей задачи линейного быстродействия // Автоматика и телемеханика. 1964. Т. 25. № 1. С. 16-22.
[50] Кирин Н. Е. Методы последовательных оценок в задачах оптимизации управляемых систем. Л.: издательство Ленинградского государственного университета, 1975.
[51] Кротов В. Ф., Гурман В. И. Методы и задачи оптимального управления. М.: Наука, 1973.
[52] Крылов И. А., Черноусько Ф. Л. О методе последовательных приближений для решения задач оптимального управления // Журнал вычислительной математики и математической физики. 1962. Т. 2. № 6. С. 1132-1139.
[53] Крылов И. А., Черноусько Ф. Л. Алгоритм метода последовательных приближений для задач оптимального управления // Журнал вычислительной математики и математической физики. 1972. Т. 12. № 1. С. 14-34.
[54] Левин А. Ю. Линейные оптимальные быстродействия и центрированные сечения // Вестник Ярославского университета. 1975. Вып. 12. С. 87-93.
[55] Ли Э. Б., Маркус Л. Основы теории оптимального управления. М.: Наука, 1972.
[56] Любушин А. А. Модификации и исследование сходимости метода последовательных приближений для задач оптимального управления // Журнал вычислительной математики и математической физики. 1979. Т. 19. № 6. С. 1414-1421.
[57] Любушин А. А. О применении модификации метода последовательных приближений для задач оптимального управления // Журнал вычислительной математики и математической физики. 1982. Т. 22. № 1. С. 30-35.
[58] Любушин А. А., Черноусько Ф. Л. Метод последовательных приближений для расчета оптимального управления // Известия АН СССР. Техническая кибернетика. 1983. № 2. С. 147-159.
[59] Милютин А. А., Илютович А. Е., Осмоловский Н. П., Чуканов С. В. Оптимальное управление в линейных системах. М.: Наука, 1993.
[60] Моисеев Н. Н. Численные методы теории оптимальных управлений, использующие вариации в пространстве состояний // Кибернетика. 1966. № 3. С. 1-29.
[61] Моисеев Н. Н. Численные методы в теории оптимальных систем. М.: Наука, 1971.
[62] Моисеев Н. Н. Элементы теории оптимальных систем. М.: Наука, 1975.
[63] Орлов М. В. Линейная задача быстродействия: численный алгоритм // Вестник Московского университета. Серия 15. Вычислительная математика и кибернетика. 1986. № 4. С. 41-46.
[64] Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем со многими неизвестными. М.: Мир, 1975.
[65] Охоцимский Д. Е., Энеев Т. М. Некоторые вариационные задачи, связанные с запуском искусственного спутника Земли // Успехи физических наук. 1957. Т. 63. № 1а. С. 36-51.
[66] Полак Э. Численные методы оптимизации. Единый подход. М.: Мир, 1974.
[67] Поляк Б. Т., Третьяков Н. В. Метод штрафных оценок для задач на условный экстремум // Журнал вычислительной математики и математической физики. 1973. Т. 13. № 1. С. 34-46.
[68] Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищенко Е. Ф. Математическая теория оптимальных процессов. М.: Физмат-гиз, 1969.
[69] Пропой А. И. Элементы теории оптимальных дискретных процессов. М.: Наука, 1973.
[70] Розоноэр Л. И. Принцип максимума Л. С. Понтрягина в теории оптимальных систем II // Автоматика и телемеханика. 1959. Т. 22. № 11. С. 1441-1457.
[71] Рокафеллар Р. Выпуклый анализ. М.: Мир, 1973.
[72] Скоков В. А. Алгоритм метода штрафных оценок для минимизации функций многих переменных при наличии ограничений // Программы и алгоритмы. ЦЭМИ АН СССР. 1973. Вып. 55.
[73] Срочко В. А. Вариационный принцип максимума и методы линеаризации в задачах оптимального управления. Иркутск: Издательство Иркутского университета, 1989.
[74] Срочко В. А. Итерационные методы решения задач оптимального управления. М.: Физматлит, 2000.
[75] Срочко В. А., Мамонова Н. В. Квазиградиентный метод решения задач оптимального управления // Известия высших учебных заведений. Математика. 1996. Т. 415. № 12. С. 84-91.
[76] Срочко В. А., Хамидулин Р. Г. Метод последовательных приближений в задачах оптимального управления с краевыми условиями // Журнал вычислительной математики и математической физики. 1986. Т. 26. № 4. С. 508-520.
[77] Табак Д., Куо Б. Оптимальное управление и математическое программирование. М.: Наука, 1975.
[78] Тихонов А. Н., Галкин В. Я., Заикин П. Н. О прямых методах решения задач оптимального управления // Журнал вычислительной математики и математической физики. 1967. Т. 7. № 2. С. 416-423.
[79] Тодд М. Дж. Вычисление неподвижных точек и приложения к экономике. М.: Наука, 1983.
[80] Тятюшкин А. И. Численное решение задач оптимального управления // Дифференциальные уравнения и численные методы. Новосибирск: Наука, 1986. С. 208-217.
[81] Тятюшкин А. И. Численные методы и программные средства оптимизации управляемых систем. Новосибирск: Наука, 1992.
[82] Тятюшкин А. И. Алгоритм поиска оптимального управления в задаче линейного быстродействия // Алгоритмы и программы решения задач линейной алгебры и математического программирования. Иркутск: Издательство Иркутского гос. университета. 1979. С. 115-128.
[83] Федоренко Р. П. К обоснованию метода вариаций в фазовом пространстве для численного решения задач оптимального управления // Журнал вычислительной математики и математической физики. 1969. Т. 9. № 6. С. 1396-1402.
[84] Федоренко Р. П. Приближенное решение задач оптимального управления. М.: Наука, 1978.
[85] Черноусько Ф. Л., Баничук Н. В. Вариационные задачи механики и управления. Численные методы. М.: Наука, 1973.
[86] Черноусько Ф. Л., Колмановский В. Б. Вычислительные и приближенные методы оптимального управления // Итоги науки и техники. Математический анализ. 1977. Т. 14. С. 101-166.
[87] Шатровский Л. И. Об одном численном методе решения задачи оптимального управления // Журнал вычислительной математики и математической физики. 1962. Т. 2. № 3. С. 488-490.
[88] Шевченко Г. В. Алгоритмы решения некоторых задач оптимального управления для линейных систем // Новосибирск, 1990. 29 с. (Препринт / РАН. Сибирское отделение Институт математики; № 11).
[89] Шевченко Г. В. Численный алгоритм решения линейной задачи оптимального быстродействия и его модификация // Новосибирск, 2002. 27 с. (Препринт / РАН. Сибирское отделение Институт математики; № 93).
[90] Шевченко Г. В. Численный алгоритм решения линейной задачи оптимального быстродействия // Журнал вычислительной математики и математической физики. 2002. Т. 42. № 8. С. 1184-1196.
[91] Шевченко Г. В. Линейная задача оптимального управления с выпуклым однородным функционалом // Фундаментальная и прикладная математика. 1999. Т. 5. № 3. С. 757-763.
[92] Энеев Т. М. О применении градиентного метода в задачах оптимального управления // Космические исследования. 1966. Т. 4. № 5. С. 651-669.
[93] Allgower E., Georg K. Simplicial and continuation methods for approximating Fixed Points and Solutions to systems of equations // SIAM Review. 1980. V. 22. ь 1. P. 28-85.
[94] Allgower E., Georg K. Numerical Continuation Methods. Introduction. New York, Berlin: Springer Verlag, 1990.
[95] Balakrishnan A. V. On a new computing technique in optimal control // SIAM J. Control. 1968. V. 6. № 2. P. 149-173.
[96] R. O. Barr. An efficient computational procedure for a generalized quadratic programming problem // SIAM J. Control. 1969. V. 7. № 3. P. 415-429.
[97] Bulirsch, R. Miele, A. Stoer, J. Well, K. H. Optimal Control. Calculus of variations, optimal control theory and numerical methods // ISNM. International Series of Numerical Mathematics. V. III. Basel: Birkhaeuser, 1993.
[98] Caetano M. A. L., Yoneyama T. New iterative method to solve optimal control problems with terminal constraints //J. Guad. Control Dyn. 1996. V. 19. № 1. P. 262-264.
[99] Dennis J.E., More Jr. & Jorge J. Quasi-Newton Methods, Motivation and Theory // SIAM Review. 1977. V. 19. b 1. P. 46-89.
100] Eaves B.C. Computing Kakutani Fixed Points // SIAM J. Appl. Math.
1971. V. 21. № 2. P. 236-244.
101] Eaves B. C. Homotopies for Computation of Fixed Points // Math. Prog.
1972. V. 3. № 1. P. 1-22.
102] Eaves B.C., Saigal R. Homotopies for Computation of Fixed Points on Unbounded Regions // Math. Prog. 1972. V. 3. № 2. P. 225-237.
103] Gay D.M., Schnabel R.B. Solving Systems of Nonlinear Equations by Broyden's Method with Projected Updates // Nonlinear Prog. 1978. № 3. P. 245-281.
104] E. G. Gilbert. An iterative procedure for computing the minimum of a quadratic form on a convex set // SIAM J. Control. 1966. V. 4. № 1. P. 61-80.
105] Hestenes M. R. Multiplier and gradient methods //J. Optimization Theory Appl. 1969. V. 4. № 5. P. 303-320.
106] Hohenbalken B. A finite algorithm to maximize certain pseudoconcave functions on polytopes // Math. Prog. 1975. V. 9. P. 189-206.
107] Laan G., Talman A. J. J. Simplicial Fixed Point Algorithms. Amsterdam: Mathematical Centre, 1980.
108] Mangasarian O. L. Nonlinear programming. McGraw-Hill, New York,
1969.
109] Meyer G. L. Accelerated Frank-Wolfe Algorithms // SIAM J. Control. 1974. V. 12. № 4. P. 655-663.
110] Pecsvaradi T., Kumpati S. Narendra. A new iterative procedure for the minimization of a quadratic form on a convex set // SIAM J. Control.
1970. V. 8. № 3. P. 396-402.
111] Polak E. An historical survey of computational methods in optimal control // SIAM Review. 1973. V. 15. № 2. P. 553-584.
112] Rockafellar R. T. Augmented Lagrange multiplier functions and duality in nonconvex programming // SIAM J. Control. 1974. V. 12. № 2. P. 268-285.
113] Saigal R. On The Convergence Rate of Algorithms for Solving Equations that are Based on Methods of Complementary Pivoting // Math. Oper. Res. 1977. № 2. P. 108-124.
114] Saigal R., Todd M.J. Efficient Acclleration Techniques for Fixed Point Algorithms // SIAM J. Num. An. 1978. № 15. P. 997-1007.
115] Scarf H. The Approximation of Fixed Points of a Continuous Mapping // SIAM J. Appl. Math. 1967. V. 15. № 5. P. 1328-1343.
116] Seywald H., Kumar R. R. Some recent developments in computational optimal control // IMA Vol. Math. Appl. 1997. V. 93. P. 203-233.
117] Todd M.J. Improving the Convergence of Fixed Point Algorithms // Complementarity and Fixed Point problems / Balinski M.L., Cottle R.W. Math. Prog. Study. 1978. № 7.
118] Yamashita, Y., Shima, M. Numerical computational method using generalized algorithm for optimal control problem with terminal constraints and free parameters //J. Nonlinear Anal., Theory Methods Appl. V. 30. № 4. 1997. P. 2285-2290.