Метод динамического программирования

Метод динамического программирования был разработан Р. Беллманом. Он применим для решения задач оптимизации систем управления по различным критериям.

Пусть объект управления описывается системой вида

                                     .                                (1.34)

Критерий оптимальности имеет вид (1.7).

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

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

Введем в рассмотрение функциональное уравнение Беллмана

                        .                   (1.35)

В соответствии с изложенным, оптимальное управление остается оптимальным для каждого из рассматриваемых участков, поэтому (1.35) можно записать в виде

                 .            (1.36)

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

Рисунок 1.5 – Пояснение принципа оптимальности Беллмана

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

Пусть t – фиксированное значение времени, а Δ t – некоторый отрезок времени. Тогда представим уравнение Беллмана (1.35) в виде

         .    (1.37)

На основании принципа оптимальности динамического программирования (1.37) можно записать в виде

     . (1.38)

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

         .     (1.39)

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

                            .                       (1.40)

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

              .         (1.41)

                      .                 (1.42)

Далее фиксируется значение x и минимум в (1.41) вычисляется по тем значениям u, для которых точка , вычисляемая по (1.42), соответствует целевому положению. Таким образом по значению  можно приближенно определить значения  в некоторой области пространства состояний. Так как на интервале  управление , то найдя функцию , находим приближенно и .

Последующие шаги находятся аналогично.

Если сделать предположение о дифференцируемости функции , то можно записать уравнение Беллмана в непрерывной форме:

              .         (1.43)

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

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

                                      .                                 (1.44)

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

                                .                           (1.45)

В данной задаче объект (1.44) и функционал (1.45) не зависят от времени, поэтому производная .

Запишем уравнение (1.43), с учетом (1.44) и (1.45):

                     .                (1.46)

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

              .         (1.47)

Из уравнения (1.47) находим оптимальное управление:

                                     .                                (1.48)

Подставим выражение (1.48) в (1.46):

               .          (1.49)

Преобразуем (1.49):

                         .                    (1.50)

Решаем уравнение (1.50) относительно функции Беллмана , в результате чего получим:

                               .                          (1.51)

Из полученных двух решений асимптотическую устойчивость обеспечивает только одно:

                               .                          (1.52)

Подставляя (1.52) в (1.48), получаем окончательное выражение для оптимального управления:

                                   .                              (1.53)

1.6. Проектное задание 2

Для объекта управления, описываемого уравнением из табл. 1.2, используя уравнение Беллмана (1.43), синтезировать оптимальное по критерию, заданному в табл. 1.2, управление. Промоделировать полученную систему в Matlab.

Таблица 1.2 – Варианты заданий для метода динамического программирования

№	Уравнение объекта	Критерий
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

Продолжение таблицы 1.2

22
23
24
25
26
27
28
29
30

При моделировании вычислить значение критерия качества. Изменить полученный коэффициент усиления регулятора на +50 % и -50 % и убедится, что найденное управление доставляет минимум заданному критерию качества.

При моделировании использовать пример программы, приведенный ниже.

 

%---------------------------------------------

clc

clear all

close all

 

x0=[3;0];

 

[t,x]=ode45('bellman_optimal1',[0 6],x0);

 

figure(1); hold on; grid on; plot(t,x(:,1),'b');

figure(2); hold on; grid on; plot(t,x(:,2),'b');

%-------------------------------------------------

Вначале производится очистка экрана, переменных и закрытие графических окон.

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

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

В конце программы осуществляется построение переменной состояния и значений критерия качества.

ode-функции с именем bellman_optimal1, оформляемая в отдельном файле с тем же именем, представлена ниже.

%---------------------------------------------

function y=bellman_optimal1(t,x)

u=(1-sqrt(2))*x(1)*1.0;

y=[-x(1)+u; u^2+x(1)^2];

%-------------------------------------------------

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

Далее определяется оптимальное управление.

В конце функции задаются правые части уравнений объекта.

Моделирование необходимо провести в три этапа. На первом этапе управление умножается на коэффициент, равный 1.0 (см. вторую строку функции bellman_optimal1). При этом параметр, задающий цвет графика в операторе plot, равен 'b', что означает синий цвет.

После моделирования на экране появляются графики с изменением во времени переменной состояния и критерия качества.

Далее меняем коэффициент, на который умножается управление на 0.5 и 1.5, комментируем первые три строки в основной программе моделирования, изменяем цвет графиков и проводим моделирование снова.

В результате получаем графики, аналогичные графикам, представленным на рис. 1.6 и 1.7.

На рис. 1.6 и 1.7 синим цветом представлены результаты моделирования оптимальной системы, красным цветом моделирование при уменьшенном на 50 % коэффициенте усиления, а зеленым – при увеличенном на 50 % коэффициенте усиления регулятора. Видим, что изменение управления приводит к увеличению значения критерия качества.

Рисунок 1.6 – Переменная состояния системы (1.44), (1.53)

Рисунок 1.7 – Критерий качества (1.45) системы (1.44), (1.53)

 

Линейной программирование

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

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

Необходимо найти такие значения переменных x ₁, x ₂,..., x _n, при которых целевая функция

                          .                     (1.54)

имеет минимальное (максимальное) значение при ограничениях, заданных равенствами:

                       .                  (1.55)

где  – постоянные числа.

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

Пусть, например, задана целевая функция двух переменных

                                     .                                (1.56)

и совокупность неравенств

                                      .                                 (1.57)

Такая двумерная задача имеет геометрическую интерпретацию на плоскости. Для построения допустимой области проведем прямые x ₁ + x ₂ = 1, x ₁ = 0 и x ₂ = 0. В результате получим треугольник, показанный на рис. 1.8.

Проведем теперь прямую Q = 2, показанную на рис. 1.9. Видим, что она попадает в допустимую область. Если строить линии Q > 2, то они не будут попадать в требуемую область, поэтому можно не рассматривать эти линии и двигаться в сторону уменьшения целевой функции. Чтобы построить остальные линии, достаточно перенести прямую Q = 2 параллельно самой себе. Из рис. 1.9 видим, что в крайней точке x ₁ =0, x ₂ = 0 целевая функция Q достигает минимума, равного 0.

Рисунок 1.8 – Допустимая область

Рисунок 1.9 – Допустимая область

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

Рассмотрим симплекс-метод на примере. Пусть задана целевая функция

                    .               (1.58)

и условия

                                 .                            (1.59)

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

Это требование означает, что область, образованная условиями, ограниченна и не пуста, поэтому задача линейного программирования всегда имеет решение. Например в системе (1.59) базисными переменными являются x ₂, x ₅ и x ₄. Индексы базисных переменных определяют числа λ₁ = 2, λ₂ = 5, λ₃ = 4. Определение индексов базисных переменных является первым шагом симплеск-метода.

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

                            .                       (1.60)

При этом все переменные рассматриваются в положительной области.

Далее осуществляется построение симплекс-таблицы, которая для рассматриваемого примера имеет вид табл. 1.3.

Таблица 1.3 – Таблица симплекс метода

	1	3
2 5 3	3 -1 1	2 3 -4	5 9 1
	-5	35	-109

В первом столбце находятся индексы базовых переменных x ₂, x ₅, x ₄. В первой строке находятся индексы остальных переменных x ₁ и x ₃. В последнем столбце находятся коэффициенты правых частей уравнений системы (1.59). В последней строке находятся коэффициенты перед свободными переменными x ₁ и x ₃ в новой целевой функции (1.60) и свободный коэффициент, взятый со знаком минус. Наконец, в последней строке таблицы находятся коэффициенты перед свободными переменными x ₁ и x ₃ в уравнениях системы (1.59).

Если в уравнениях (1.59) положить свободные переменные x ₁ и x ₃, равными нулю, то получим точку, называемую вершиной: x ₁ = 0, x ₂ = 5, x ₃ = 0, x ₄ = 1, x ₅ = 9. Эта точка соответствует канонической форме задачи линейного программирования.

В общем случае симплекс-таблица имеет вид табл. 1.4.

Таблица 1.4 – Общий вид таблицы симплекс метода

	λm+1 … λn
λ1 … λm	a1λm+1 … a1λn         … amλm+1 … amλn	b 1 … b m
	pλm+1 … pλn	- Q 0

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

Если указанный критерий не выполняется, то нужно перейти из одного базиса в другой. Этот переход осуществляется следующими этапами.

1. Выбор разрешающего столбца. Среди элементов последней строки таблицы выбирается любой отрицательный p _j* < 0. Лучше всего выбирать наименьший коэффициент. Номер этого столбца j * называется разрешающим.

2. Выбор разрешающей строки. Если в выбранном столбце все коэффициенты a _ij* < 0, то минимума не существует. В противном случае необходимо для положительных коэффициентов вычислить отношения b _i/ a _ij*. Строка i, для которой указанное отношение минимально, называется разрешающей строкой i *. При этом элемент a _i*j* называют разрешающим элементом.

3. Замена базиса с помощью разрешающего элемента. Она осуществляется путем применения следующих правил:

                            .                       (1.61)

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

                                     .                                (1.62)

В соответствии с формулой (1.62) новый разрешающий элемент является обратной величиной к старому разрешающему элементу.

                            .                       (1.63)

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

                                   .                              (1.64)

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

                            .                       (1.65)

                                     .                                (1.66)

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

                    .               (1.67)

                          .                     (1.68)

Рассмотрим переход в новый базис на приведенном примере.

1. В табл. 1.3 в последней строке только один отрицательный коэффициент: −52. Этот коэффициент находится в столбце номер 1, поэтому j * = 1. Разрешающий столбец равен [3 −1 1]^T.

2. В разрешающем столбце два положительных элемента: a ₂₁ = 3 и a ₄₁ = 1. Вычисляем два числа: b ₂/ a ₂₁ = 5/3 и b ₄/ a ₄₁ = 1/1. Последнее число меньше, поэтому строка 4 является разрешающей: i * = 4. Разрешающий элемент выделен жирным шрифтом в табл. 1.3: a ₄₁ = 1.

3. Приступаем к созданию новой таблицы. В первом столбце на месте, где стояла цифра 4, ставим 1, а в первой строке вместо цифры 1 ставим 4. В результате получаем таблицу 1.5.

Таблица 1.5 – Индексы таблицы симплекс метода

	4	3
2 5 1	* * *	* * *	* * *
	*	*	*

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

                                    .                               (1.69)

По (1.62) вычисляем новый разрешающий элемент

                                  .                             (1.70)

4. Далее по формуле (1.63) вычисляем элементы нового разрешающего столбца. Берем разрешающий столбец и умножаем его на −1. Элементы полученного столбца заносятся в новую табл. 1.6 (разрешающий элемент уже находится в таблице).

5. В соответствии с формулой (1.64) в последней строке умножаем −5 на −1 и результат заносим в последнюю строку новой таблицы.

6. Производим вычисления по (1.65) и (1.66). Берем разрешающую строку [1  −4] и умножаем ее на 1. Полученная строка заносится в новую таблицу за исключением уже вычисленного разрешающего элемента. По формуле (1.66) берем 1 в последнем столбце и умножаем ее на 1. Тогда получаем табл. 1.6. Звездочками обозначены не вычисленные элементы.

Таблица 1.6 – Таблица с вычисленными разрешающими строкой и столбцом

	4	3
2 5 1	-3 1 1	* * -4	* * 1
	5	*	*

7. Далее вычисляем остальные элементы. По (1.67) получаем элементы , . Отметим, что в старых элементах номера i* = 4, j* = 1, а в новых элементах (с крышкой) i* = 1, j* = 1.

8. По (1.68) получаем .

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

                                    .                               (1.71)

В соответствии с (1.71): , .

Свободный элемент Q ₀ определяется в соответствии с выражением

                                  .                             (1.72)

Согласно последнему выражению получаем .

Результат преобразования показан в табл. 1.7.

Таблица 1.7 – Преобразованная таблица

	4	3
2 5 1	-3 1 1	14 -1 -4	2 10 1
	5	15	-104

Из последней таблицы видим, что все коэффициенты  больше нуля, поэтому она соответствует минимуму функционала. Точка минимума определяется равенствами: x ₄ = 0, x ₃ = 0, x ₂ = 2, x ₅ = 10, x ₁ = 1.

Для решения задачи линейного программирования в Matlab существует функция linprog. В общем случае решение задачи с использованием данного оператора имеет вид

                        .                   (1.73)

Задача (1.73) решается программой вида

x = linprog(f,A,b,Aeq,beq,lb,ub)

Если какие-то условия не заданы, то соответствующие матрицы и вектора задаются пустыми.

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

clc

clear all

close all

 

f=[48 13 2 17 3];

A=[];

b=[];

Aeq=[3 1 2 0 0; -1 0 3 0 1; 1 0 -4 1 0];

beq=[5;9;1];

lb=zeros(5,1);

X = LINPROG(f,A,b,Aeq,beq,lb)

В первых трех строках программы очищается экран, закрываются графические окна и очищаются все переменные.

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

В результате полечено следующее решение

X =

 

1.0000

2.0000

0.0000

0.0000

10.0000

Видим, что решение совпадает с выкладками, сделанными вручную.

 

1.8. Проектное задание 3

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

Таблица 1.8 – Преобразованная таблица

 

Проверить результат. Для этого вычислить функцию Q, при полученном оптимальном решении. В данном случае она равна 104. Далее изменяем полученное оптимальное решение следующим образом. Увеличим переменную x ₁ = 2, а переменную x ₂ = 0 оставим прежней. Подставим эти значения в (1.59). В результате получим:

–>

Очевидно, что переменная x ₃ = -0.5 отрицательная, значит нужно уменьшить переменную x ₁ =0.5, а переменную x ₂ = 0 оставим прежней. Снова подставим эти значения в (1.59). В результате получим:

–>

Решая систему, получаем: x ₁ =0.5, x ₂ = 0, x ₃ = 1.75, x ₄ = 7.5, x ₅ = 3.75.

Вычисляем функцию Q по (1.58): 166.25 > 104.