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

МАТЕМАТИКА

УДК 519.6

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

Е. А. Андреева, И. С. Мазурова

APPLICATION OF FAST AUTOMATIC DIFFERENTIATION FOR OPTIMIZATION OF SYSTEMS OF INTEGRO-DIFFERENTIAL EQUATIONS

E. A. Andreeva, I. S. Mazurova

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

In this paper, the authorbi solve the optimal control problem in the predator-prey model described by a system of Volterra integro-differential equations. The maximum principle for the initial continuous problem is formulated considering the defined constraints and functional form. The initial continuous problem is reduced to the discrete optimal control problem. For the discrete optimal control problem formulas for calculating the gradient of the target function using the method of direct differentiation and the method of fast automatic differentiation are received. The numerical method is developed to construct an approximate optimal solution on the basis of fast automatic differentiation method. It is shown that the obtained numerical results of the approximate optimal control correspond to the maximum principle for the initial continuous problem with a prescribed accuracy.

Ключевые слова: оптимальное управление, интегро-дифференциальные уравнения, модель хищник-жертва, метод быстрого автоматического дифференцирования.

Keywords: optimal control, integro-differential equations, prey-predator model, fast automatic differentiation method.

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

В ВЦ РАН под руководством Ю. Г. Евтушенко [1] разработана методология быстрого автоматического дифференцирования (БАД-методология), позволяющая с единых позиций определять градиенты для явно и неявно определенных функций и для вычислительных процессов, которые являются результатом дискретной аппроксимации непрерывных систем, описываемых дифференциальными и интегро-дифферен-циальными уравнениями. Этот подход был применен для построения оптимального решения систем, описываемых интегро-дифференциальными уравнениями типа Вольтера.

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

f m Л

X (t) = x (t) ег -Z ед (О IV i=i )

-Z ълУj(t) X(t) - Ri)- ui(t x

j=i

i = 1, m,

f n Л

-aj -ТслУ1(t) +

l=1

y , (t) = y, (t)

d,i X(t) - Ri'

V l=1

m 1

+y,(t)ZrJi J Fji(t-т)X(0-R)dT-vj(t),

1=1 t-

j = 1, n,

при заданных начальных условиях

X(0) = x, X (в) = (pt (в\ в e [-r, 0], y, (0) = y0 i = 1, m, j = 1, n.

(1)

(2)

(3)

jv s j

Здесь m, n - количество классов хищников и жертв соответственно, xi (t), i = 1, m - численность

популяции жертв i-ого класса, у, (t), j = 1, n - численность популяции хищников j-ого класса в момент времени t e[0,^], Ri - численность популяции

жертв i-ого класса, недоступная хищникам,

Е. А. рАндреева, И. С. Мазурова

47

МАТЕМАТИКА

e,, a j, a,,, bji, cji, dji, у j, - действительные

положительные коэффициенты, характеризующие взаимодействие популяций хищников и жертв. Функции распределения Fji (t — т) характеризуют состояния системы в прошлом, описывая влияние доступной для хищника жертвы на характерном отрезке времени, связанном с ростом популяции жертв. Функциями управления являются функции u, (t) - скорость отлова популяции жертв i-ого класса, Vj (t) - скорость

отлова популяции хищников j-ого класса, удовлетво -ряющие ограничениям

0 ^ U (t) ^ U max , 0 ^ Vj (t) ^ Vj

г = 1, m 3 = 1, «,

(4)

где Uj max, Vj max - максимальные скорости отлова популяции жертв и хищников соответственно, далее

U = (U1,...,um), v = (V1,...,vn).

Математическая модель взаимодействия популяций с произвольным конечным числом хищников и жертв (1) - (3) является обобщением модели взаимодействия двух популяций, построенной в работе [3].

Оптимальное управление отловом строится из условия оптимизации заданного функционала T

J (u, V) = f fo (t, x, y, u, v)dt +

о (5)

+Ф( x(T), y(T)).

В работе [2] показано, что оптимальное управление U (t), V (t), t G [0, T] в задаче (1) - (5) удовлетворяет принципу максимума:

max

0^i^Umax,0^Vj^max |_

-Aofo(t, x (t), y(t X а,щ)— £ рг (t )Щ — £ (t )щ,

j=1 j=1

m ___

— ЛЛ (t, x (t), У(tX u (tX v (t)) — £ Рг(t )ui(t) — £ rj(t )Vj (tX

г=1

j=1

(6)

сопряженные функции

Рг (t), rj (t), г = 1, m, j = 1, n

являются решением системы интегро-дифферен-циальных уравнений:

( m

P(t) = —p(t) Ie,—£ ax(t)

V i=1

m n

+£ Pl (t)aUX, (t) + Рг (t)£ Ь3гУ3 (t) ы j=1 (7)

—£ rj(t) yj(t) dj, -

j=1

n l+I

—£ r3 (t) y3 (t у f F3, (t—t j=1 m,

j=1 t

n

r3(t) = r3 (t )j3 + r3(t )£ c3,y,(t) +

i=1

m

+£ p i(t )bb (xi(t)—R)+

i=1

n m

+£ r(t )y(t )cj-— r(t )£ di(x(t)—R)—

—r3 (t)£yi f fj, (t—t)(x, (T—R )dT 3 =1 n,

i=1 t—r

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

СФ

pt (T) = —Л — (x(T), y(T)), i = 1,

Cx,

Pi (t) = 0 при t > T,

m,

/ ч СФ

p (T) = — A0 —(x(T), y(T)), 3 = 1, n,

dy3

rj (t ) = 0 при t > T.

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

Дискретная задача оптимального управления состоит в минимизации функции:

(8)

(9)

I = A t £ Fx( x*, y*, uk , v*) +

к = 0

+ F0 (* , yq )- inf,

при заданных рекуррентных соотношениях:

(10)

P+1 = x* +Atx*

f m

e, —£

V l=1

a,lxl —

(11)

—At£ bj,yk (^— R,)—Atu*, j=1, m,

j=1

yj+1 = yj + Atyj ■

•(—aj—£ c3,yik+£ dfl (^—Ri)) ■

(12)

l=1 l=1

m к—1

+At2yj£yj £ FXs • (xf — Ri) — Atvj, j = 1,/

i=1 s=*—v

начальных условиях

x = P (ti Xl = —^,0, yj =^; i = 1, m, j = 1, n

(13)

48 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 2

МАТЕМАТИКА

и ограничениях на управление

0 < ик < и ,0 < vk < v. ,

(14)

1 = 1, m, j = 1, n.

Для рекуррентных соотношений (11) - (12) введем следующие обозначения:

xl = F,(l, X11, Y1l,U1l, V1l), (15)

yl = F2( l, X21, Y21, U21, V21), (16)

, Y ,U ,V - заданные наборы векторов

X , yl , ul , vl соответственно, i = 1,2 , l принимает значения от 1 до q, в общем виде будем писать l е D.

Будем говорить, что соотношения (15)-(16) задают q-шаговый процесс. Формулы (15)-(16) определя-l l 7

ют вектора x , у на l -ом шаге.

Для каждого j е D введем индексные наборы

Q1 j , R1 j , S1 j W1 j , Q 2 j , R 2 j , S 2 j W 2 j , содержащие индексы всех векторов x l , y l иul ,Vl l е D принадлежащих соответственно наборам X 1 j , Y 1 j и U 1 j , V 1 j в (15):

Q1 j ={i е D : x1 е Xj},

S1 j ={i е D : U1 е Uj},

R1 j = {i е D : у‘ е Yj},

W1 j = {i е D : : v1 е Vj},

X 2 j , Y 2 j и U 2 j , V 2 j в

Q2j = {i е D : x1 е Xj},

S2j = {i е D : U1 е Uj},

R 2 j = {1 е D : у‘ е Yj},

W2j = {i е D : v1 е Vj}.

На j-шаге множества Q j, S1j, R1 j, W^- опре

и

деляют индексы всех «входных» векторов xi , и i ,

у i , v i , содержащихся во множествах X J , UJ ,

Y j , Vj соответственно, необходимых для вычисления правой части (15), аналогично на j-шаге множества Q2 j , S2 j, R2 j, W2 j задают множества

всех индексов «входных» векторов xi , u. , у 1 , v., необходимых для вычисления правой части (16). Введем сопряженные индексные наборы:

Q1 j = {1 е D : xj е X1},

S1 j = {i е D : uj е U1},

R1 j = {i е D : у1 е YJ},

W1 j = {i е D : v1 е V j },

Q2j = {i е D : yJ е Y1},

S2j = {i е D : vj е V1},

R2j = {i е D : у1 е YJ},

W2j = {i е D : v1 е Vj}.

На j-шаге множество Q1j определяет номера

всех тех шагов 1 е D , для которых xJ е X1 .

Между элементами введенных множеств есть связь. Если г е Q\ j, s е Sj j , то j е QXi, j е Sxs.

С помощью введенных множеств можно записать:

X1 j = {x1 : i е Qx j}, U1 j = {u1 : i е Sx. }.

Из приведенных определений следует, что если

xq е X1 j, уq е Y 1 j, то имеет место следующие функциональные зависимости:

xq = F1 (q,..., xj,...),

у4 = F2(q,..., у1,...). Следовательно, наборы

Q1 j , R1 j , S1 j W1 j , Q 2 j, R 2 j , S 2 j W2 j могут быть названы наборами индексов входных переменных x, у и U, v соответственно, в то

время как Q1 J', R1 j- S1 j W1 J, Q 2 j, R 2 j,

S 2 j W 2 j называются наборами индексов выходных векторов x, у и U, v соответственно.

Определим матрицы производных процесса по управлениюЛ

NUf = dxjT / диг, N2j = ду]Т / 5иг

, (17)

MUj = dxjT / dv, M2j = ду^ / dv1

где

Nm = ^ui-FU (j, X], Yj ,Uj ,Vj) + = Z N1,^1^(j',Xj,Yj,Uj,Vj) +

кеЙ1 j

+ Z N2,Fl (j, Xj, Yj ,Uj ,Vj),

,еа ,

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 2 49

МАТЕМАТИКА

м 2 j=S2 jFlt (j, XJ, YJ, UJ, VJ) + + Z м 2 lkFTlvt (j, XJ, YJ, UJ, V1) +

k g R

2 j

+

Z м ulFlk (j, X1, Y1, U1, V1),

k gQ.

2j

1, i g S 2

1, i g S i : I -ь * <„> 2 ,■

i e S1 j7 Я7 = | o, i e S2,,

мij = Z MmFl, (j,X1, Y1,U1, Vj) +

k gQ1 j

+ Z MMFT (J, X1, Yj,U j, V1),

k gRi

«2, = Z N (j, Xj, Yj, Uj, V‘)

k gR

+

2 j

+ Z Ni,k-FLk (j,Xj,Y‘,Uj, Vj).

8x'+u 8x'

M

1i+1 j 8vJ 8vJ

8

dvJ

8uJ

(1 + At(e - ax') - Atax1 -

-AtbTy) + -8y- (-AtbT x' ).

8yi+1

M2i+1 j = ~dv j = ~At5i+1 j

dy_

dvJ

(1 + At(-a - cTy1 + dT (x1 - R)) - AtcT y

Ti

i-1 8v'

+Arf Z Fi-‘(xs - R)) + T7(Aty’dT)

q-1 oxs

+At2 Z —-ysyTFs1,

Lu8v]

s=i+1 1

N.

8y+1 8У

2i+11 8uJ 8uJ

+dT ( x1 - R)) + AtylcT +

i -1

(1 + At (-a- cTy

8xi

+AtfT Z Fi-s (xs - R)) + — (Aty'dT)

s=i-„ 8uJ

q-1 8xs +At2 Z —yyTFs =+1 8и]У Г

я J1, i + 1 = J

1+11 lo, i + 1 * j'

Градиент целевой функции по u l, l = 0, q - 1 для процесса (15) - (16) запишется следующим образом:

d Q ( и, v )

du

= I ,,i (x, У, u, v) +

+Z N uk1 xk(x, y, u, v) +

k=1 q

+Z N 2 ikIy.(x, y, u, v).

k=1

(18)

Градиент целевой функции по v l , l = 0, q - 1 для процесса (15) - (16) запишется следующим образом:

d Q (u, v )

dvl

= I,(x, y, u, v) +

k gQ2 j

Для процесса (15) - (16) формула (17) запишется следующим образом:

.. 8x1+1 8x',.. i

n1'+1J = ~8йг = -AtS'+1J+ ~й7(1+At (e - ax) -

- Atax1 - AtbTy1) + -8y- (-AtbTx1),

+ Z M 2lkIy k (x, y, u , v) +

k = 1 q

(19)

+Z м ш1^(x, y, u, v).

k=1

Основное вычислительное время при расчете градиента по формулам (18) - (19) затрачивается на нахождение матриц N 1ij , N 2 ij , M 1ij , M 2 ij , для

этого необходимо решить 2q(n + т) раз систему 2 г \2

уравнений с q (n + т ) неизвестными.

Введем множители Лагранжа pl, rl, тогда функция Лагранжа запишется в следующем виде:

L( x, y, u, v, p, r ) = I (x, y, u, v) +

+Z [ FT (l, X1, Yl ,Ul ,Vl) - xlT ]pl +

l=1

+Z [ FT (l, X1, Yl, Ul ,Vl) - ylT Y.

l=1

Множители Лагранжа определяются из решения системы уравнений:

L(x, y, u, v, p, r) = у (x, y, u, v) +

+ Z FlxT(l,Xl,Yl,Ul,Vl)pl +

IgQu

+

Z F2/(l,X1,Yl,Ul,Vl)rl -p1 = 0,

,gQ2'

Lyy (x, y, u, v, p, r) = Iy/ (x, y, u, v) +

+ Z V (l, Xl, Y1 ,Ul ,Vl)pl +

IgRu

+

Z F2/(l,X1,Yl,Ul,Vl)rl -r1 = 0.

1gR?i

l

50 Вестник Кемеровского государственного университета 2014 № 4 Т. 2

МАТЕМАТИКА

Тогда множители Лагранжа p\, r'k, i = 1, да, j = 1, n, k = 0, q для процесса (15) - (16) запишутся следующим образом:

к

Pk =Jxk (x, У, u, v) + pk+1 +Mpf+1

Xj

n n

к+1 V'' I ,.k , A.V''„k+1 .k

j=1

j=1

да > да

-Z ailxl -A'ZPf+lalixl -

l=1 l=1

k+v n

A2 Z Z rs+1 Vsy-Fsrk i

s = k+1 j=1

(20)

rj = ^k (x,y,u,v) + fj + - Atrj + -Д'.+ Хзд

l =1

да n m

A'ZPk+4(xl -Rl)-A'Zrf+1yfcij + Atrj+1 £dji(xf -Ri) +

l=1 l=1 l=1

да к-1

,2.к+1^ X"' d-sx.s

+ A'2 rf+1 Z^jl Z Fji- S (x/ - Rl), j = 1, n.

l=1 s=k-r

(21)

Градиент целевой функции по

uk, i = 1, да, k = 0, q - 1

для процесса (15) - (16) запишется следующим образом:

dQ(u,v)

duk = Lu) ( X, У , U , V, P , r ) =

= I uk (x, У, U, v) + Ap+1.

(22)

Градиент целевой функции по

Vj , j = 1, n, k = 0, q - 1

для процесса (15) - (16) запишется следующим образом:

d Q (u, v)

------k— = Lv*(x> У>u >v> P,r) =

dvj

= I vk (x, У, u, v) + A 'r

(23)

k +1

Основное вычислительное время при расчете градиента по формулам (22) - (23) тратиться на нахождение векторов p l и rl, l = 0, q , для этого необходимо решить 2 системы уравнений с q (n + да) неизвестными.

Выражения для вычисления градиента «прямым» и «обратным» методами математически эквивалентны, но с вычислительной точки зрения имеется существенно! отличие. Показано, что для расчета градиента по формулам «прямого» метода дифференцирования необходимо решить 2q(n + да) раз систе-2 2

му уравнений с q (n + да) неизвестными. Для

расчета градиента по формулам БАД необходимо решить 2 системы уравнений с q(n + да) неизвестными.

Для нахождения градиента по формулам прямого метода дифференцирования и метода БАД необходимо вычислить:

Для прямого метода дифференцирования: Для метода БАД:

q 2 да 2 значений dx jT / 5ui , q 2 даn значений dy j / д ui, q 2 даn значений д xj / д vi, q 2n 2 значений dyjT / дvi, dl (u ) qда значений , du‘ qn - dI (u) qn значений , dv1 qда значений p 1 , qn значений r i , dI (u) qда значений , du‘ qn й dI (u) qn значений , dvi

22 всего q (n + да ) + q (n + да ) значений всего 2q(n + да) значений

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 2 51

МАТЕМАТИКА

Таким образом, применение прямого метода вычисление градиента целевой функции в q (n + m) раз более трудоемко, чем расчеты по формулам БАД.

На основе метода быстрого автоматического дифференцирования разработана программа, позволяющая построить оптимальное решение задачи (10) -(14). Ниже на рисунках 1 - 4 представлены графики зависимости оптимального управления ux (t), Vj (t) и численности популяции Xj (t) yj (t) при различных значениях запаздывания r при следующих параметрах системы: T=20; m = 1, n = 1,

X10 = 6, У10 = 5, R = 0, ^ = 0.75, cn = 0, dy = 0, ax = 0.75, уц = 0.1, ац = 0.0375, b11 = 0.0375, 0 < ut < 0.1, i = 1,2,

At = 0.01,

A = 7, B = 13,5, M = N = 1000,

I(u, v) = M ■ max( - xq,0) +

+N ■ max ( - yq ,0 ),

точность метода e = 0,0000001.

Рис. 2. График зависимости У1 (t) при различных значениях запаздывания r

ол -и10% -

■ 4il(t)opt_i-0.1

П ПЛ ul O.lopt. i=1.1

П ГЫ

п m u 1 (t) op t_i=1.5

п

с 5 10 15 20 t.

Рис. 3. График зависимости ul (t) при различных значениях запаздывания r

52 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 2

МАТЕМАТИКА

■----'У 1 (t)opt_i=0.1

_ V1 (t) о-p t_l=1.1

1-----v 1 (t) op t_i=1.5

20 t.

Рис. 4. График зависимости vx (t) при различных значениях запаздывания r

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

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

в q (n + m) раз более трудоемко, чем расчеты по формулам быстрого автоматического дифференцирования, где q (n + m ) - размерность управления. На

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

Литература

1. Евтушенко Ю. Г. Оптимизация и быстрое автоматическое дифференцирование. М.: Научное издание ВЦ РАН, 2013. 144 с.

2. Andreeva E. A., Mazurova I. S. Optimal control of predator-prey model with distributed delay // Mathematical Modelling and Geometry. 2013. Vol. 1. № 3. P. 38 - 48. Режим доступа: http://mmg.tversu.ru

3. Beretta E., Capassa V., Rinaldi F. Global stability results for a generalized Lotka_Volterra system with distributed delays. Journal of Math. Biology. 1988. № 26.

Информация об авторах:

Андреева Елена Аркадьевна - доктор физико-математических наук, профессор, заведующая кафедрой компьютерной безопасности и математических методов управления математического факультета Тверского государственного университета, [email protected].

Elena А. Andreeva - Dcotor of Physics and Mathematics, Professor, Head of the Department of Computer Security and Mathematical Management Methods, Tver State University.

Мазурова Ирина Сергеевна - соискатель кафедры компьютерной безопасности и математических методов управления математического факультета Тверского государственного университета, IrinaSMazurova@gmail. com.

Irina S. Mazurova - post-graduate student at the Department of Computer Security and Mathematical Management Methods, Tver State University.

(Научный руководитель - Е. А. Андреева)

Статья поступила в редколлегию 22.07.2014 г.

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 2 53