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

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

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

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

Итерационная формула для нахождения коррекций &u^J (t) искомого управления u^J+1(t) имеет вид

.  , _    .    . дН

AuJ(t) = kJG 1 д j ф,    kJ = const > 0, где kJ - длина шага, G - постоянная матрица влияния (зачастую единичная), Н = H(t, х, u, ^) = ^т f(t, х, и) - функция Гамильтона.

Достоинством градиентного метода является очень большая область сходимости, поэтому даже при очень плохих начальных приближениях u⁰(t) искомый минимум, как правило, находится без затруднений. С другой стороны, итерационный процесс может привести к соседнему минимуму, и тогда встает вопрос о выборе приближения u⁰(t) более близкого к оптимальному u * (t). Более того, метод обычно медленно сходится вблизи искомого минимума, поэтому к вопросу о выборе начального приближения все же стоит подходить основательно.

Extr-H-метод. Extr-H-метод представляет собой модификацию метода градиента, призванную улучшить его плохую сходимость в окрестности минимума J_min за счет того, что коррекция ku¹ (t) выбирается не пропорционально градиенту, а вычисляется непосредственно из необходимого условия максимума функции Гамильтона H .

Здесь коррекции вычисляются по формуле

Auj(t) = uJ*(t) — uJ(t), или hui(t) = kJ[uJ*(t)

—

uJ(t)],     kJ < 1,

если ku^J(t) принимает слишком большие значения (u^J (t) - аналитически найденное управление).

Определенная трудность в применении Extr-H-метода состоит в том, что часто не представляется возможным выразить аналитически управление uJ (t) и уже на этом этапе приходится приступать к итерационному процессу. Именно этим можно объяснить тот факт, что метод используется относительно редко.

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

Поправка Ди¹ (t) здесь вычисляется как

s'

Метод сопряженных градиентов устраняет проблему медленной сходимости вблизи вычисляемого минимума J, не требуя аналитически выражать u¹ (t) как в Extr-H-методе. Поэтому он является наиболее оптимальным из рассмотренных.

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

Идея метода состоит в поиске такой коррекции △u⁷(t), которая минимизировала бы главную линейную часть 6J = (grad h(x⁰,T),y(ty) приращения функционала J (здесь y(t) есть поправка по траектории фазовой координаты).

Упрощенное регулирование смешанной популяции. Рассмотрим хорошо известную модель Т. Мальтуса, описывающую динамику изолированной популяции:

x(t) = ax(t), x(0) = x0 > 0, где x(t) - численность популяции в момент времени t, a - показатель изменения ее численности, связанный с уровнями рождаемости и смертности (если a > 0, то численность растет, иначе – убывает). Очевидно, что эта модель приводит к экспоненциальному росту населения и не применима в долгой перспективе. Однако в качестве жесткой модели ее можно использовать в расчетах на недолговременный период.

Теперь допустим, что у некоторых членов этой популяции может появиться признак, который выделит их в отдельную популяцию, динамика которой отлична от динамики первой. Тогда, заменяя х на х1, модель примет вид x1(t') = (a- p)x1(t),    х1(0) = х10 > 0,

X2(t) = fixi(t),     х2(0) = х2о > 0, где x2(t) - численность малой популяции в момент времени t, /3 - показатель изменения численности этой популяции аналогичный а.

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

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

J = -(х1(Т) + х₂(ТУ) ^ min.

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

X1(t) = (а-/)X1(t) + Yiu(t), х1(0) = хю > 0, X2(t) = 3X1(t) - Y2u2(t),    х2(0) = х2о > 0, где Y1, Y2 — некоторые коэффициенты, вообще говоря, зависящие от текущей численности популяции. Но эти зависимости мы выяснять не будем.

Найдем для начала аналитическое решение задачи. Функция Гамильтона имеет вид:

Н(1,х,и,-ф) = ^1((а — /)х1 + Y1u) + ^₂(/х₁ — y₂u²).

Тогда оптимальное управление u*(t) вычисляется как

u*(t) =

Y1^1 2Y2^2

Сопряженная система имеет вид:

^1 = —(a — /W1

3^2,    ^1(Т) = 1,

^2(Т) = 1.

^1= 0,

Решая ее, получим

^1⁽t) = ^ _ _a (/ — ae^ ^a)(T °),    ^2(0 = 1.

Таким образом:

u*(t) = k(p - ae^(p-a)(T-t)),    к =     ^-^ '

Теперь перейдем к расчетам. Установим такие значения констант: a = —4 • 10^-3, p = -10^-3,

xi0 = 10 0 0, x20 = 40 0,    -i = -2 = 1, T = 1.

Подставляя эти значения в выражение для оптимального управления, получим

4e3*io-3(i-t) — i u*(t) =-------------------.

Этому управлению отвечает значение функционала/(1) = 1396,29.

Заметим сразу, что на исследуемом интервале [0;1] оптимальное управление u*(t) меняется очень слабо - функция монотонно убывает от u*(0) = 0.502 до u*(T) = 0.5. Основываясь на этом, зададим довольно «плохое» начальное приближение u⁰(t) = 10.

Для численного интегрирования будем использовать хорошо известную схему Хойна для уравнений типа у = /(x, у)

y;+i = У; + h * k1 + k2,     ki = /(xi,yi),     k2 = f(xi+i,yi + h * ki), имеющую второй порядок аппроксимации. Однако, основываясь на постановке задачи, имеет смысл запоминать номер итерации L*, на которой уже не происходит изменение целой части функционала J. Для вычисления интегралов будем использовать формулу трапеций:

b

[            (b — a) z 4

J /(x)dx « ^-A/(a) + /(b))

a с оценкой остаточного члена

Л(/) =

max|/"(x)|(b 12a) . [a,b]             12

Экспериментально выяснено, что достаточно хорошим значением длины шага является k = 0.5. Но, забегая вперед, стоит упомянуть, что при этом ее значении уже с помощью обычного метода градиента, требуемая точность (при h = 0.01 – до четвертого знака после запятой) достигается уже на пятой итерации (а целая часть функционала перестает изменяться на четвертой). Поэтому возьмем длину шага k = 0.4.

Для применения метода Шатровского было выбрано допустимое множество U = [0.5; 2.5] и допустимое значение функционала /доп = 1397, причем из-за особенностей итерационного процесса функционал оценивался по абсолютному значению.

Результаты вычислений представлены в таблице (L - количество итераций, J^L -значение функционала J, u^L(t) - значение найденного управления в момент времени t).

Полученные результаты полностью подтверждают теоретические рассуждения – ExtrH-метод показал наилучшие результаты, метод сопряженных градиентов показал более быструю сходимость, нежели обычный метод градиента. Единственным явным несоответствием с аналитическим решением является то, что табличная функция u^L(t) возрастает на выбранном интервале исследования, в отличие от убывающего управления u'(t).