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

УДК 517.977 ББК 22.193

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

Батурин В. А.1

(ФГБУН Институт динамики систем и теории управления СО РАН, Иркутск) Черемных С. В.2

(Евразийский лингвистический институт в г. Иркутске (филиал ФГБОУВПО «Московский государственный лингвистический университет»), Иркутск)

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

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

1. Введение

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

1 Владимир Александрович Батурин, доктор физико-математических наук, профессор (rozen@icc.ru).

2 Светлана Викторовна Черемных, кандидат физико-математических наук (sv157@rambler.ru).

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

2. Задача оптимального управления и алгоритмы последовательных улучшений

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

Пусть задана управляемая система дифференциальных уравнений на отрезке [¿о, ¿1]: ёх Ж

где состояние х(0 е Я" - непрерывная и кусочно-дифференцируемая вектор-функция, а управление и(1) е и с Ят - кусочно-непрерывная вектор-функция.

Множество пар (х(0, и(1)), удовлетворяющих перечисленным условиям, называется множеством допустимых В. Предполагается, что В Ф 0.

На множестве В задан функционал ч

,и)

*о

Требуется минимизировать функционал I на множестве В, т.е. найти такую последовательность {(хДО, иДО)} с В, чтобы

выполнялось I(х,и) ^ I при 5 ^ да.

о

Поскольку в дальнейшем пойдет речь об алгоритмах второго порядка, то на функции £ £0, ^ накладывается условие, что они дважды дифференцируемы.

(1) = f (*, х, и), х(*о) = Хо ,

-1

I (х, и) = ^Х х(*г)) +1 f 0(*, х, и)ё*.

Поиск состояния х(0 и управления и(^), минимизирующих заданный функционал, осуществляется с помощью итерационной процедуры, в которой на каждом шаге решается задача улучшения: заданы начальные значения х1^) и требуется найти новые состояние Х\() и управление и11^) такие, что выполняется 1(х°, и111) < 1(х1, и1).

Для поставленной задачи разработана группа приближенных методов, основанных на принципе расширения. В работах [1, 3-8, 10, 14, 15, 16] описаны методы сильного и слабого улучшения для данной задачи и получены соответствующие алгоритмы. При выводе этих алгоритмов применен метод локализации, например, в [1, 10] конструируется положительно определенный функционал ^х1, и1, х, и), имеющий смысл штрафа за отклонение от опорной траектории, а затем задача улучшения исходного функционала сводится к задаче улучшения вспомогательного функционала

/„ (х, и) = аI(х, и) + (1 - а)3(х1, и1, х, и), где а е [0, 1] - скалярный параметр. В этих работах функционал ^х1, и1, х, и) задается в следующем виде: - в методе сильного улучшения:

(2) 3 (х1, и1, х, и) = —

| Ах'ЕАхЛ + [Ах(^)]' ЕАх(^)

- в методе слабого улучшения:

(ч

^Аи'ЕАиЖ + [Ах(^)] ' ЕАх(^)

(3) 3 (х1, и1, х, и) = 1

где Е - единичная матрица соответствующей размерности (тхт или пхп), Ах = х - х*((), Ли = и - и^(^). Введенный параметр а служит средством регулирования близости соседних процессов (х1, и1) и х, и11), и от правильного выбора значения этого параметра зависит успешность итерации. В реальных прикладных задачах нередко степень влияния разных компонент вектора управления на весь процесс существенно отличается, кроме того, могут на порядки отличаться и масштабы различных компонент управления и состояния. Поэтому имеет смысл в формулах (2) и (3) заменить единичные матрицы на диагональ-

ные, элементы которых выполняют роль весовых коэффициентов. Тогда функционал /(х7, и7, х, и) примет вид - в методе сильного улучшения:

I Ах' рьхйг + [Аг(^ )]'уАх^)

(4) 3ру (х1, и1, х, и) =1

2

V '0

- в методе слабого улучшения:

(5) 3(х1, и1, х, и) =1

1 1 > | Аи' рАийг + [Ах(^)] 'уАх(гх)

Л

где Р - диагональная матрица размерности п х п (или т х т) с положительными элементами, у - диагональная матрица размерности п х п с неотрицательными элементами.

Алгоритмы сильного и слабого улучшения второго порядка, описанные в работе [1], с учетом того, что функционал /(х7, и7, х, и) задается формулой (4) или (5), имеют вид:

Алгоритм сильного улучшения.

1. Задаем управление и7^), из системы (1) находим х7(().

2. Задаем значения параметров: а = а1, р = р1, у= у1.

3. При этих значениях параметров из системы

У=-на-она-на),

(6) ^ = "На + (1 -а)Р-оИах - Н^о-оН^о ,

у('х) = -а¥ х(х7 (',)), о('х) = -а¥х (х7 (',)) - (1 -а)у вычисляем У) и о(1).

Здесь У) - п-вектор, а(0 - п х п-симметрическая матрица, На (', х, р, и) = р'/ (', х, и) -а/0 (', х, и), На(',х,р) = тахНа(г,х,р,и),

р - п-вектор, все производные функции Н а вычисляются в точке х7^), у(0), производная На - в точке х7^), У), и7^)).

4. Из системы

= f{t,x,u{t,x,p)), x(t0) = x0,

dx dt

где

ù(t,x,p) = argmaxH"(t,x,p,u),

n a '

ugU

1,

p = f(t) + &(t)Ax, Ax = x - x (t), вычисляем x'Xt).

5. Находим

U11 (t) = u(t, x11 (0, wit) + <j(f)(xF (0 - x7 (0)).

6. Если Дх77, и71) > Iix7, и1), то изменяем значения параметров а, Д, уи повторяем процесс с пункта 3.

Алгоритм слабого улучшения.

1. Задаем управление u7(t), из системы (1) находим x7(t).

2. Задаем значения параметров: а= а, Д = Д, у= /.

3. При этих значениях параметров из системы

f=-на+(на + °ки )(на - (1 - а)дг на,

(7)

dg=-Hi-она - н>+(на+ощи ) х х (нна - (1 -а)рг\на + нао

ыщ

yit,) = -aFx(X (tj) , ait,) = -aFxx ( x1 (tj) - (1 -a)y вычисляем Щ) и a(t). В системе (7) производные функции Ha вычисляются в точке (t, xI(t), ip(t), uI(t)). 4. Из системы dx dt

где й = и1 + Mit, Ах) ,

■ = f(t,x,u), x(t0) = x0,

Au(t, Ax)=-(на - (i - а)?)-[на+(нао+на^,

Ax = x - x1 (t), вычисляем x'Xt). 5. Находим

u11 (t) = u1 (t) + Au(t, x11 - x1). 6. Если ¡(x11, u11) > I(x\ u1), то изменяем значения параметров а, Д, у и повторяем процесс с пункта 3.

3. Управление параметрами алгоритма

В вышеприведенных алгоритмах от выбранных значений а, Д, у зависит весь ход итерационного процесса: существование решения системы (6) или (7) на всем отрезке [to, ti], успешность итерации (т.е. удалось ли улучшить функционал), глубина улучшения функционала на итерации и, следовательно, необходимое количество итераций. Сформулируем задачу управления выбором параметров алгоритма: имеется элемент (xI(t), uI(t)), требуется найти такие значения параметров а, Д, у, чтобы найденные при этих значениях новые траектория и управление были наилучшими (в смысле наименьшего значения функционала).

Для решения сформулированной задачи рассмотрим задачу оптимального управления, решаемую с помощью алгоритма улучшения, как задачу с векторным параметром а, где a = (а, Д, Д, ..., Д^или m), у, У2, ..., у). После выполнения пунктов 1-5 алгоритма сильного или слабого улучшения получаем: решается система

& = f(t,x,u(t,x,a)), x(t0) = х0;

находятся X1 и U = u(t, X1, a); вычисляется значение функционала

ч

I{x" ,u") = F{x" (ij)) + J f (t,x11,u(t, x11,a))dt .

Введем обозначения:

/(t, x, u(l, x, a)) = g(t, x, a), /° (t, x, u{t, x, a)) = g° (t, x, a). Тогда получаем следующую задачу:

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

ti

J (a) = I (x, a) = F (x(t)) + j g0 (t, x, a)dt,

= gх, а), х(?о) = х0,

где х(У) - непрерывные и кусочно-дифференцируемые функции, удовлетворяющие на [¿о, ¿1] системе дифференциальных уравнений

ёх Ж

а - параметр-вектор размерности 2п + 1 (или т + п + 1).

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

В результате получаются следующие выражения для За 3^ и 3Г] при а = а, / = у= у1.

В алгоритме сильного улучшения:

Ч г -1-1

/оч ^=\н:^,хн,ф,ин){наии^,хП ,Ч, + а(х11 -х'Хи11)] X

XИи0,х11,уа+&а(х" -х1),и")ёГ,

Ч !

3Л =\н;а,хп,1у,ип).\наииа,хп,¥^а(хп -х'Хи11)] X (9) I 1 J

х И° (¿, х , ^ + (х" - х1), и г = 1,..., п;

4 г -1-1

Зг=\н;{1,х\у,,и11).\наии{1,х11,¥ + о{х11-х1),и11)\ х (Ю) ' * 1 J

И0 Ц, х11, + (х" - х1), и11 ,

] = 1, ...,п .

В алгоритме слабого улучшения:

За=-\Н'и (г, х" ,цг,иПу 1_Наии -(1-а1 )/3' X к 1

х[Нии (г,х1 ,¥а,и1) + $ ]• [Иаии - (1 - а 1]-1 х

(11) х\наа + (н^а+нах хх1 - х ^-[н: - (1 - а1]-1 х

х[ Ни (г, х ,Уа, и) + (Наи¥а а +

+Них (г, х ,Га ,и1))(х11 - х1)]}Л;

3Л = -\Н'и(г,х",) ■ \наии -(1-а1 )Р1 х к 1

X[н1 (г,У) - (1 - а1 )Ёг ]X[на -(1-а1)/З1

(12) х[на + (Н^а + На )(хи - х1 )]-[на - (1 - а1]-1 х х [Н1 (г, х1 , и1) + (Н^ь +

+Н0х(г,х1 ,и1 ))(хп - х1)]}сИ, I =

зУ] = -)н'и а, х11, тУ) ■ \на - а- «7)р1 ]-1 х к 1

хН0ш (г, х1 , и1) •[на - (1 - а1 ]-1 х

(13) х[на + (Н^а + на )(хл - х )]-[на - (1 - а1 ]-1 х

х" Н0(г, х1 , и1) + (Н^ +

+Н0ш(г,х1 ^,и1 ))(хп - х1)])а, ] = 1п. В формулах (8)-(13): н(г, х, р, и) = р' • /(г, х, и) - /0(г, х, и), На (г, х, р, и, а) = р' • /(г, х, и) - а/0 (г, х, и);

н°(г,х,р,и) = р' ■ /(г,х,и),

функция !//(/ ) находится из системы

где

H(t, х, р, а) = p'git, х, а) - g° (t, х, а),

a = (а1 ,р/ ,у/ ), i = 1,...,n(m), j = 1,...,n;

производные функции На, аргументы которых опущены, вычисляются в точке (t, x1(t), y(t), u1(t)); E, — квадратная матрица m x m, у которой элемент e,i = 1, а все остальные элементы равны нулю.

В формулах (8)-(10) функции y(t) и o(t) находятся из системы (6) при а= а, р = р1, у= у11, а в формулах (11)—(13) функции y(t) и o(t) находятся из системы (7). Функции у а (t), Oa(t), Ур (t), Ор (t), фу. (t), Оуj (t) находятся из систем,

полученных дифференцированием системы (6) (или системы (7)) по а, р, у, соответственно. Эти системы имеют следующий вид.

В методе сильного улучшения:

^ Уа = "На " Нфа " Оа (На " H^) " "о(Н0а + Н0уУа ),

d Оа ="[На,а + Ну " Р " О " О^а +

+ Н1х:шУа ] " [НЛа + НашшУа ]o " НЛ0а " 0 Уфу0 "

(14)

~'~Нухула] [" 'хуа 11 1хуула]0 ' 'ху0 а О а' 'уу' "О[Науа + ^ЛуууУа ]О " ОНауОа ,

Уа(Ч) = "F(х7 (tj) ,

Оа(^) = "Fxx (X7 (О) + У.

° (На" к) -ону

(15)

- он*w„ - На„„„уй О -

\хуур НхууУр0

Чу0р ~0pt НУу0 " ^УллуЛр0 0 " 0Нау0pi ■. 14

¥Л(к) = - а х1 (Г!)) = О, (г:) = 0 .

С гу, = -н Г - % (н а -н а) - %н аГ с % =-н аш н »-она

(16) ш

н 1ггУ%~ н 1Г°У, -

у,н ГГГ0 - оНТГу, % - оНггоУ>

сту^1) = -{\-а)Ё].

В методе слабого улучшения:

С Га = -Нх (г, х1 Га, и ) + [НХи (г, х' Га, ^ ) + %аКи ] х х [Н»ы - (1 - а )РГ • на - [Наш + оН»ы ] х х [наи - (1 - а Ж1 • [Нии (г, х1 Га, и1) + р] х х [Н а - (1 - а Ж1 • на + Ни + она ] х х[наи -(1 -а)Р] х • Ни(г,хГа,и1),

С %а=-Нхх х1 Га , и' ) - %аКх - НГ%а + (17) + [Нхи (г, х1 Га , и1) + ОаНа ] • [Щи - (1 - а)Р]- х х [Н»х + щ%] - [нхи + она ] • [наи - (1 - а ЖХ х

х [Нии (г, х1 Га и1) + Р] • [Щи - (1 - а х

х [нах + НГ%] + [Нахи + она ] • [Наии - (1 - а)Р]- х

х[Н их (г, х Га ,и1 ) + Н а ] ,

Га (к) = -Рх (х (г,)), %а (г,) = -Р„ (х1 (г,)) + У.

С ¥р=-н {г, х Г, и1) + [Н0и (г, х Г, и1) + Щи ] х х|Я;;;( -(1 -а)р\1 • наи -\на + он;и]• |я;;;, -(1 -а)р\1 х

х[Н(1,х,1//,и)-(\-а)Ё]-[Н -(1 -а)Р\-Н +

+ [Н1 + оЩи ] • [Я;и - (1 - а)РГ ■ К (V, X, у) _

С= -НхОх(г,х1 г,и1)-оЛна -Нах%Л +

+ [^ (г,х1 Гр„,и1) + орНа] •[наи - (1 - аЖ1 х х [н"х + я;;у| - [Я» + он;и ] • |я;;;( - (1 - а)р\1 х X [ЯИ°И (г, X, ^ у) - (1 - ] • [Н"и - (1 - аЖ' х (18) х [нах + наг%] + [Нахи +%на ] х

х [Н:и - (1 - а Ж1 • [НОх (г, х1 Г, и) + Н Гг%ь ]

Г (О = 0, (О = о.

С Гу=- Н0Х (г, х г, и1)+[Н0„ (г, х1 г, и1) + OyJ н;и ] х

х [ нат - (1 - а ж 1 • на - [ на + она ] • [на - (1 - а ж - х

х на (г, х1 Гу] и1) • [на - (1 - а)р]- • наа + +[на + она ] • [на - (1 - ажх • но (г, х1 г ,и1).

С о у, =-К (г, х Г и1) -оуНГх - НГоу] +

+ [Н°хи (г,х1 Г ,и1) + оу1 на] • [на - (1 - а)Ж]-1

х

(19)

х [на+наго] - [на+она] • нии - (1 - аж1 х

х на (г, х1 Гу] и1) • [на - (1 - а)р]- • [на+нг+

Г ]

+нхи +она ] • [на - (1 - а)рг1 • н (г, х г„и1

+Н Го, ]

^(0 = 0, сгу^1) = -{1-а)Ё].

В системах (14)—(19) выражение Нахху/уа понимается как матрица размерности п х п, элементами а- которой являются суммы

д¥к .

= У На

у ^ j да ,

k=1

аналогичным образом понимаются выражения H а Щц,уо, H а xwwYa, H ацт,Ца, матрица Et (Ej) - квадратная матрица nxn, у которой элемент e„ = 1 (ец = 1), а все остальные элементы равны нулю; все производные функции H а вычисляются в точке (t, x1, Ц); производные функции На, аргументы которых опущены, вычисляются в точке (t, x1, ц, U).

Чтобы, решая полученную задачу минимизации, не интегрировать заново систему (6) (или систему (7)) при каждом изменении значений параметров, можно вычислять функции Ц) и o(t) приблизительно, разложив их в ряды Тейлора по параметрам в окрестности точки ( ai, р, /):

V(t, а11, р11 у11 ) = W(t, а1, р1 у1 ) + ца (t, а1, р1 У ) х

n( m)

(20) х (а11 - а1 ) + £ Цр (t, а1 ,р1 у1 )(Pf -р ) +

i=1

n

+ (t, а1 ,р1 ,у1 )(У-у),

j=i j

o(t, а11 р11 у11 ) = o(t, а1 ,р1, У ) + oa (t, а1 ,р1 У ) х

n(m )

(21) х (а11 - а1 ) +^Ор (t, а1, р1 у1 )(р« -р ) +

i=1

n

+ ^o (t, а1 ,р1 у1 )(У-у).

j=i j

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

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

1. Задаем управление ииз системы (1) находим х1(().

2. Задаем значения параметров: а = а1, ¡ = (З1, у= у1.

3. Из системы (6) при а 1, З1, у1 вычисляем 1^(0 и 0(3).

4. Из системы

^ = /^,х,йЦ,х,р)), х(^) = х0,

где

й^,х,р) = а^тахН" (¿,х,р,и),

иШ

р = ^(г) + о(г)Ах, Ах = х - х7 (г), вычисляем Х1^).

5. Находим

г/7 (0 = X77 (О, И» + ^(ОСх77 (0 - X7 (0)).

6. Вычисляем За, 3(- и Зу. в точке ( а 1, 31, у1) по формулам (8)-

(10), при этом производные у и с по параметрам находятся из систем (14)-(16). Задаем

а11 = а1 -%3а, 37=33 -£Э(, у =У-$3У], ] = 1,2,..., и ,

где скалярный параметр. Переходим к решению одномерной задачи минимизации функционала 3( а11, З11, у11) по переменной £ . Функции и О) при а = а11, ¡ = 311, у= у11 вычисляем по формулам (20)—(21).

7. Находим и при о". /3". ;/". у/(/. с/1, /3". у"). 0(1(, а1,31, у11), а затем заново вычисляем Х1^) и и71^).

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

1. Задаем управление и^(0, из системы (1) находим х1(^).

2. Задаем значения параметров: а = а 1, ¡ = 31, у= у1.

3. Из системы (7) при а 1, 31, у1 вычисляем и О).

4. Из системы

<3х &

где

% = /Ц,х,й), хЦ0) = х0,

и = и1 + ААх),

Аи(Г, Ах) = ~(Наии - (1 - а)Р) 1 |Я,Г + (Я> + Я»)Аг],

Ах = х - х1 (г), вычисляем Х1^).

5. Находим

и11 (г) = и1 (г) + Аи (г, х11 - х1).

6. Вычисляем За, Зр- и Зу■ в точке (а1, р1, у1) по формулам

* ]

(11)—(13), при этом производные у и о по параметрам находятся

из систем (17)—(19).

Задаем

а 11 = а1 -40а, Р?=РЖ~43р, г = 1,2,..., т,

У? =у1 , ] = П где 4- скалярный параметр. Переходим к решению одномерной задачи минимизации функционала З( а11, р11, у11) по переменной 4 . Функции Г) и с() при а = а11, р = рш, у= у11 вычисляем по формулам (20)—(21).

7. Находим и при а". /?". ;/". у/(/. а". /?". ;/"). а(г, а11, р11, у11), а затем заново вычисляем Х1^) и и11^).

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

Алгоритм с Окигп = 0

1. Задается управление и7^), из системы (1) находится х7(^).

2. Задаются параметры оО = 0 и ( = (З1 > 0.

3. Из системы

Л Га = -Их Х Га, и' X

Л °а= -Ихх (^ Х Га , и' ) ИГх & Х Г, ^ ) -

(22) - И^, х1 Г, и1 )а а,

Га (О = - Ех (X (О),

(О = - Рхх (X&)); находятся ЦТоаЦ, 0) и ОиЦ, 0).

4. По формуле

ч

За (0,3' ) = -|Ии' (1, X' Га ,и' ) • (¡' )- • Ии (1, X' Га ,и'

Ч

вычисляется 0а(0, ¡7).

5. Задается а е (0, 1].

6. Находится значение а111 как решение задачи одномерной минимизации функционала 0 по параметру а. При этом функции ГЛ а) и о(1, а) находятся по приближенным формулам:

у/(г, а) = Га (г,0) • а, о(г, а) = о о (г, 0) • а .

7. Из системы ск Лг

где

.1 , к..г* Л„ „М а

и = и + Аи(1,Ах,а ,/?),

Аи(г, Ах, а11) = -(И° - (1 - а 11)¡' )-1 х

х [ ИО + (Иаи¥о(г, а11) + ИХ )Ах],

Ах = х - X , вычисляется Х7(0, а затем находится

и

' (г) = и' (г) + Аи(г, х11 - х1, а11, ().

Производные функции На вычисляются в точке

(Г, х1, у, а11), и1, а11).

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

Работоспособность и эффективность модифицированных алгоритмов были проверены на ряде тестовых примеров, в которых осуществлялось управление параметром а. Описания некоторых из этих примеров и результаты их решения исходными (далее они названы базовыми) и модифицированными алгоритмами можно найти в работах [2, 11, 12].

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

Пример 1.

ёх _ и 1

< Ж , 1(х, и) = ( (х2(0+и 2(г))Ж, г <е [0,1].

х(0)=1 о

Начальное управление и7^) = 0.

Результаты вычислений при а^ап = 1, шаге интегрирования Н = 0,001, точности вычисления функционала е = 0,001 представлены в таблице 1.

Таблица 1.

Базовый алгоритм Модифицированный алгоритм

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 1 1 0,7335367 0 1 1 1 0,7335367

Результаты вычислений при а^ап = 0,7, шаге интегрирования Н = 0,001, точности вычисления функционала е = 0,001 представлены в таблице 2.

Таблица 2.

Базовый алгоритм Модифицированный алгоритм

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 2 3 0,7 0,7 0,7 1 0,7460582 0,7360622 0,7355737 0 1 0,999958 1 0,7335368

Пример 2.

ёх,,

■ = и

I(х,и) = | (и(г) - х(г) +г)2ёг, г е [0,1].

ёг

о

Х ( 0 ) = 0

Х ( о ) = 1

Начальное управление и7^) = 0.

Результаты вычислений при а^ап = 1, шаге интегрирования Н = 0,001, точности вычисления функционала е = 0,001 представлены в таблице 3.

Таблица 3.

Базовый алгоритм Модифицированный алгоритм

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 2 1 1 0,3333333 0,8744572Е-16 0,3455601Е-29 0 1 2 1 1 0,3333333 0,8744572Е-16 0,3455601Е-29

Результаты вычислений при а^ап = 0,7, шаге интегрирования Н = 0,001, точности вычисления функционала е = 0,001 представлены в таблице 4.

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

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

Таблица 4._

Базовый алгоритм Модифицированный алгоритм

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 2 3 4 5 6 О О О О О О I 7 1 ООЗЗЗЗЗЗЗ 0,1068978 ОООЗ6З5779 ООО124116З 0^04239099 0^01447999 0^004946249 0 1 2 3 4 ОО85 ОО925 ОО9625 ОО98125 ООЗЗЗЗЗЗЗ ООО42ОО11З 0^01860381 ОО24О9146Е-О4 ОО8461745Е-О7

4. Задача управления эколого-экономической моделью

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

Рассмотрим модель замкнутой системы (отсутствует импорт и экспорт продуктов и ресурсов), описывающую совместную динамику экономики и природных ресурсов

v = Av + A& + Bz + R^f + Cv + p, p dt p dt И

df = -K(r - r*) - Dv - Ddv + z ,

dt p dt

где v(t) e Rn - вектор потока общественного продукта; z(t) e Rm - вектор интенсивности восстановления природных ресурсов; Z(t) e Rm - вектор мощностей восстановительных

отраслей; р(1) е Я" - вектор конечного непроизводственного потребления; г(1) е Ят - вектор, характеризующий состояние природных ресурсов; г - вектор естественного состояния ресурсов; А - матрица прямых затрат на производство продукта; Ар - матрица затрат на расширение производства; В - матрица затрат на восстановление ресурсов; Вр - матрица затрат на расширение мощности восстановления ресурсов; С - матрица, отражающая амортизационные затраты; К - матрица естественного восстановления и взаимного влияния ресурсов; Б, Бр - матрицы, отражающие расходы ресурсов, связанные с текущим производством и его расширением.

Очевидно, что величина производимого продукта у не может быть отрицательной и ограничена сверху некоторой величиной V (мощности производства), обусловленной объемом основных фондов и численностью промышленно-производственного персонала, т.е. выполняются ограничения: 0 < у< V; V> 0. Аналогично, интенсивность экзогенного (по отношению к природной среде) восстановления природных ресурсов не может превышать мощности восстановительных отраслей: 0 < г < Z; Z > 0. Кроме того, сильное отклонение состояния природных ресурсов от естественного может привести к необратимым для природы последствиям (если это загрязнения) или к невозможности продолжения производственной деятельности (если это запасы полезных ископаемых, сырья). Поэтому необходимо ограничить величину этого отклонения:

Я-тт < г < Ятах.

Ставится задача о поиске режима поведения эколого-экономической системы, соответствующего минимуму функционала полезности

и

I = 1

А ■ г (г) - г - А ■ р(и) ■

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

Сделаем предположения:

1. v(t) и r(t) - скалярные величины (т.е. однопродуктовая, одноресурсная модель);

2. v = V, z = Z (производство и восстановление идут на полную мощность);

3. C = 0 (амортизационные затраты включены в прямые). Поставим задачу оптимального управления следующим образом. Будем управлять интенсивностью изменения мощностей восстановительных отраслей Z(t) и конечным потреблением p, т.е. Z(t) = ui, p = иг. В результате сделанных предположений и введенных обозначений модель примет вид:

dz

— = U , dt 1

dL = [(1 - A)v - Bz - Bpu -ы2]/Ap ,

^f = - K (r - r*) - Dv - Dvdv + z . dt p dt

Добавим в функционал несколько слагаемых, имеющих значение штрафов за нарушение ограничений:

(min{0; z})4 - штраф за нарушение условия z > 0; (min{0; v})4 - штраф за нарушение условия v > 0; (v - V)2 - штраф за нарушение условия v = V; (min{0; r - Rmin})4 - штраф за нарушение условия r > Rmin. Кроме того, для регуляризации задачи добавим еще два слагаемых: si(ui)2 и S2(u2)2. Зададим отрезок времени [to, ti] и начальные значения z(to), v(to), r(to).

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

—1 = U , dt 1

dx -

= [Ax2 - Bxi - BpUi - u2]/Ap ,

£ = -K ( x3 - r*) - DX22 - Dpdr + xi,

X (0) = Xj0, x2 (0) = x°, x3 (0) = x° ,

I = | А-(хз ({) - г*)2 - и2 ■ е 5 + (шт{0; х1 }) + (шт{0; х2 }) +

о

+ (^2 -V)2 +(шт{0;(xз -Яшт})4 + sl(ul)2 + S2(u2)2 и(и) е Я2, и е [1о, 1,].

Данная задача решалась при следующих исходных данных:

■ начальный момент времени 10 = 0, конечный момент времени 11 = 1;

■ коэффициенты в уравнениях и функционале являются константами и имеют следующие значения:

А = 0,44; г* = 1,5; V = 1; В = 10; В = 0,005; Вр = 1,6; = 0,5; Бр = 0,001; Ар = 1,1; К = 0,1; А = 0,5; 5 = 0,02; ^ = 1; s2 = 1;

■ начальные значения вектора состояния х0 = 0,1, х0 = 1, х30 = 1.

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

1. Пусть планируется небольшой постоянный прирост мощностей восстановительной отрасли и постоянное непроизводственное потребление, т.е. начальное управление ы\(1) = 0,1, Ы21(1) = 0,1. Такое управление приведет к следующему состоянию эколого-экономической системы (опорный процесс):

- экзогенное восстановление природного ресурса неуклонно возрастает на рассматриваемом отрезке времени (от 0,1 до 0,2);

- производство общественного продукта резко уменьшается вплоть до отрицательных значений в конце периода;

- состояние природного ресурса улучшается (увеличение от 1 до 1,2), однако естественного уровня не достигает.

На рис. 1-5 представлены графики опорного и полученного в результате вычислений оптимального процессов.

Рис. 1. Первая компонента вектора управления

Рис. 2. Вторая компонента вектора управления

Рис. 3. Первая компонента вектора состояния

1,5 1,0 0,5 0,0 -0,5 -1,0

ч %

* ч ч

1111 н 1111111111 —; о о 1111 т О О 111111111111111111 1П ЧО о" о щи о 111111« 1 II 11111 III 00 _ о о «*

время

опорный процесс

оптимальным процесс

Рис. 4. Вторая компонента вектора состояния

Рис. 5. Третья компонента вектора состояния

Оптимальный процесс.

При заданных исходных условиях естественное восстановление природного ресурса протекает достаточно интенсивно, поэтому мощности восстановительной отрасли следует сворачивать (величина изменения мощности отрицательная и принимает значения на отрезке [0, 1] от -0,32 до -10-4). Непроизводственное потребление, наоборот, надо увеличивать (от 0,447 до 0,49). Такое управление приведет к следующему состоянию эколого-экономической системы:

- экзогенное восстановление ресурса прекратится, а в конце рассматриваемого периода начнется дополнительное потребление ресурса;

- производство общественного продукта будет оставаться приблизительно на одном уровне;

- состояние ресурса немного улучшится (увеличение с 1 до 1,05), хотя и не достигнет естественного уровня.

Итерационный процесс улучшения функционала отражен в следующих таблицах: результаты решения данной задачи алгоритмами слабого улучшения при а^ап = 1 - в таблице 5, при а1аг1 = 0,7 - в таблице 6; результаты решения данной задачи

алгоритмами сильного улучшения при а^ап = 1 - в таблице 7, при а^ан = 0,7 - в таблице 8.

Таблица 5.

Базовый алгоритм слабого улучшения Модифицированный алгоритм слабого улучшения

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 2 3 1 1 1 0,53700 -0,095142 -0,11308 -0,11309 0 1 2 3 1 0,9926 0,1799 0,53700 -0,095142 -0,11309 -0,11311

Таблица 6.

Базовый алгоритм слабого улучшения Модифицированный алгоритм слабого улучшения

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 2 3 4 5 6 О О О О О О | 7 7 7 7 7 7 1 0,53700 0,019446 -0,077500 -0,10304 -0,11021 -0,11225 -0,11284 0 1 2 3 0,85 0,99997 0,99997 0,53700 -0,039682 -0,11306 -0,11309

Таблица 7.

Базовый алгоритм Модифицированный алгоритм

сильного улучшения сильного улучшения

а Значение а Значение

1омер ерации еа ние арте ее г £ функционала и ри ец § & Й н т еа ние арте г 1 функционала

и нр со л п и нр За п

0 — 0,53700 0 — 0,53700

1 1 -0,095142 1 1 -0,095142

2 1 -0,11308 2 0,998 -0,11308

3 1 -0,11309 3 0,998 -0,11309

Таблица 8.

Базовый алгоритм Модифицированный алгоритм

сильного улучшения сильного улучшения

а Значение а Значение

и ри ец з я еа ние арте О 2 функционала и ри ец з я еа ние арте О 2 функционала

р о ет НЧ □ и чм нач амр За п р о ет НЧ □ и чм нач амр За п

0 — 0,53700 0 — 0,53700

1 0,7 -0,0013480 1 0,85 -0,048189

2 0,7 -0,088013 2 0,962 -0,11297

3 0,7 -0,10732 3 0,962 -0,11308

4 0,7 -0,11173

5 0,7 -0,11276

6 0,7 -0,11301

2. Добавим в функционал штраф за нарушение условия ¿(1) > 0, чтобы постараться не уменьшать мощности восстановительной отрасли. В этом случае функционал примет вид

I = | \ - (х3 (') - г*)2 - и2 ■ + (шп{0;х })4 +

о

+ (т1п{0;х2}) + (х2 -V)2 + (т1п{0;(хз -Л^}) + + ^ (щ )2 + ^ (и2 )2 +(т1п{0; щ })4

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

3. Вместо введения штрафа за нарушение условия ДУ) > 0 запланируем следующий режим изменения мощности восстановительной отрасли: Д0=0,01г. Тогда управляемая система и функционал примут вид

-х,

—1 = 0,01х ,

dt 1

dx -

= [АХ2 -ВХ1 -Вр • 0,01X1 -и1 ]/Ар ,

£ -К (хз - /) - ^ - ^ + Х1,

х (0) = X0, х (0) = х0, X (0) = х30, 1 г

I = | (х (t) - Г*)2 -щ • +(т1п{0;х })4 +(т1п{0;х})4 +

+ (х - V)2 +(т1п{0;(хз - ДтШ})4 + ^К)2

Управлять теперь будем только непроизводственным потреблением, т.е. р = ы\, и пусть начальное управление п\(1) = 0,1.

На рис. 6-9 представлены графики опорного и полученного в результате вычислений оптимального процессов.

0

Рис. 6. Управление

Рис. 7. Первая компонента вектора состояния

' опорный процесс оптимальный процесс

Рис. 8. Вторая компонента вектора состояния

опорный процесс оптимальный процесс

Рис. 9. Третья компонента вектора состояния

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

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

Итерационный процесс улучшения функционала отражен в следующих таблицах: результаты решения данной задачи алгоритмами слабого улучшения при а^ап = 1 - в таблице 9, при а^ап = 0,7 - в таблице 10; результаты решения данной задачи алгоритмами сильного улучшения при а^ап = 1 - в таблице 11, при а^ап = 0,7 - в таблице 12.

Таблица 9.

Базовый алгоритм Модифицированный алгоритм

слабого улучшения слабого улучшения

а Значение а Значение

Номер итерации Значение параметра функционала Номер итерации Значение параметра функционала

0 — 0,17018 0 — 0,17018

1 1 0,13655 1 0,9708 0,13654

2 1 0,13654 2 0,4757 0,13653

Таблица 10.

Базовый алгоритм Модифицированный алгоритм

слабого улучшения слабого улучшения

а Значение а Значение

Номер итерации Значение параметра функционала Номер итерации Значение параметра функционала

0 — 0,17018 0 — 0,17018

1 0,7 0,13814 1 0,993 0,13649

2 0,7 0,13656 2

3 0,7 0,13652 3

Таблица 11.

Базовый алгоритм Модифицированный алгоритм

сильного улучшения сильного улучшения

а Значение а Значение

Номер ггерации Значение параметра функционала Номер терации Значение параметра функционала

0 — 0,17018 0 — 0,17018

1 1 0,13655 1 1 0,13655

2 1 0,13654 2 1 0,13654

Таблица 12.

Базовый алгоритм сильного улучшения Модифицированный алгоритм сильного улучшения

Номер итерации Значение параметра а Значение функционала Номер итерации Значение параметра а Значение функционала

0 1 2 0,7 0,7 0,17018 0,13657 0,13652 0 1 2 0,7 0,7271 0,17018 0,13657 0,13640

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

Среди рассмотренных сценариев эколого-экономической системы наибольшую полезность (соответствующую наимень-

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

Литература

1. БАТУРИН В.А., УРБАНОВИЧ Д.Е. Приближенные методы оптимального управления, основанные на принципе расширения. - Новосибирск: Наука. Сиб. предприятие РАН, 1997. - 175 с.

2. БАТУРИН В.А., ЧЕРЕМНЫХ С В. Управление выбором параметров в алгоритмах слабого улучшения второго порядка для задач оптимального управления // Известия РАН. Теория и системы управления. - 2006. - №2. - С. 65-72.

3. ВАСИЛЬЕВ О.В., ТЯТЮШКИН А.И. Об одном методе решения задач оптимального управления, основанном на принципе максимума // Журн. вычисл. мат. и мат. физики. -1981. - №6. - С. 1376-1384.

4. ГУРМАН В.И., БАТУРИН В.А., РАСИНА ИВ. Приближенные методы оптимального управления. - Иркутск: Изд-во Иркут. ун-та, 1983. - 178 с.

5. ГУРМАН В.И., РАСИНА И.В. О практических приложениях достаточных условий сильного относительного минимума // Автоматика и телемеханика. - 1979. - №10. -С.12-18.

6. КРОТОВ В.Ф., ГУРМАН В.И. Методы и задачи оптимального управления. - М.: Наука, 1973. - 448 с.

7. КРОТОВ В.Ф., ФЕЛЬДМАН И.Н. Итерационный метод решения задач оптимального управления // Изв. АН СССР. Техн. кибернетика. - 1983. - №2. - С. 160-168.

8. ЛЮБУШИН А.А., ЧЕРНОУСЬКО Ф.Л. Метод последовательных приближений для расчета оптимального управления // Изв. АН СССР. Техн. кибернетика. - 1983. - №2. -С. 147-159.

9. Модели управления природными ресурсами / Под ред. В.И. Гурмана. - М.: Наука, Гл. ред. физ.-мат. лит-ры, 1981. -264 с.

10. Новые методы улучшения управляемых процессов /

B.И. Гурман [и др.]. - Новосибирск: Наука. Сиб. отд-ние, 1987. - 183 с.

11. ЧЕРЕМНЫХ С.В. Управление выбором параметров в алгоритмах сильного улучшения // Труды IX Международной Четаевской конференции «Аналитическая механика, устойчивость и управление движением» (Иркутск - оз. Байкал, 12-16 июня 2007 г.), Т. 3. - Иркутск, 2007. - С. 268-279.

12. ЧЕРЕМНЫХ С.В. Автоматизация выбора значений параметров в алгоритмах улучшения второго порядка // Труды XIV Байкальской Международной школы-семинара «Методы оптимизации и их приложения» (Иркутск - Северобай-кальск, 2-8 июля 2008 г.), Т. 2. - Иркутск, 2008. -

C. 211-221 (на диске).

13. Эколого-экономическая стратегия развития региона. Математическое моделирование и системный анализ на примере Байкальского региона / В.Е. Викулов, В.И. Гурман, Е.В. Данилина и др. - Новосибирск: Наука. Сиб. отд-ние, 1990. - 184 с.

14. BATURIN V.A. First-order improvement method for the problems of optimal control of logic-dynamic systems // Nonlinear Analysis. - 2006. - Vol. 64, No. 2. - P. 288-294.

15. JACOBSON D.H. New second-order and first-order algorithms for determining optimal control. A differential programming approach // J. Optimization Theory and Application. - 1968. -Vol. 2, No. 4. - P. 441-440.

16. ORTEGA J.M., RHEINBOLDT W.C. Iterative Solution of Nonlinear Equations in Several Variables. - Academic Press, NY, 1970. - 572 p.

THE SECOND ORDER METHODS FOR THE OPTIMAL CONTROL PROBLEMS

Vladimir Baturin, Institute of Systems Dynamics and Control Theory of the Siberian Branch of RAS, Irkutsk, Doctor of Science, professor (rozen@icc.ru).

Svetlana Cheremnykh, Irkutsk Euro-Asian Linguistic Institute of the Federal State Government-Funded Educational Organization of Higher Professional Education "Moscow State Linguistic University", Irkutsk, Cand.Sc. (sv157@rambler.ru).

Abstract: The paper considers the iterative improvement algorithms, the efficiency of which substantially depends on the chosen parameters values. The problem of control of these parameters is formulated and discussed. We designed the modified algorithms where parameters are automatically adjusted at each iteration. The original and modified algorithms are applied to solve the problem of optimal control for the ecology-economic system.

Keywords: optimal control, iterative improvement algorithm, strong improvement, weak improvement, adaptive parameters adjustment.

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

Поступила в редакцию 08.10.2015.

Опубликована 31.03.2016.