Оптимальная программа управления динамической системой
Постановка задачи
Рассматривается динамическая система:
,
где
x - вектор фазового состояния ДС, размерности 
; u - вектор управления ДС, размерности 
; А,В - матрицы постоянных коэффициентов системы,
размерности 
и 
соответственно.
Требуется
определить оптимальную программу управления 
системой
ДС, обеспечивающую минимум квадратичного критерия:
где
- матрицы весовых коэффициентов, определяющие
значимости соответствующих переменных 
.
Исходные
данные
программа управление динамический уравнение
1. Аналитическое решение
.1 Методика решения поставленной задачи
. Построим функцию Гамильтона. Для смешанного критерия типа (2) функция
Гамильтона имеет вид:
(4)
Здесь
- компонента сопряженного вектора, соответствующая
фиксированной координате 
.
Имея
в виду, что - Сonst по времени и в конечный момент времени 
функция Гамильтона (4) может быть переписана в
следующем виде:
. В данном случае Н.У. имеет частный вид:
Поскольку
прямых ограничений на управление в задаче не вводится
Отсюда
следует, что:
(5)
2.
Решение канонической системы уравнений. Необходимость в этом решении возникает
в связи с необходимостью получения сопряженного вектора 
.
(6)
Эта
система замкнутая (она зависит сама от себя). Решение полученной краевой задачи
для системы линейных дифференциальных уравнений может быть получено
аналитически с помощью переходной (фундаментальной) матрицы 
.
Если
ввести новые обозначения:
Тогда
каноническая система (6) примет вид:
(7)
.
Используя свойство фундаментальной матрицы можно записать для системы (7)
следующее решение:
В
развернутом виде можно записать:
где
, 
, 
, 
-
соответствующие блоки фундаментальной матрицы 
.
Поскольку
конечное условие сопряженного вектора
можно
записать следующее соотношение:
Из
полученного выражения можно определить вектор состояний сопряженных переменных:
Теперь
векторы 
и 
могут
быть поставлены в зависимость только от начального состояния вектора 
. А именно:
где
Таким
образом, в соответствии с оптимальной структурой управления (5) можно получить
искомую программу оптимального управления в виде:
.2
Процедура решения поставленной задачи
.
Составляем матрицу расширенной системы 
.
.
Находим собственные числа матрицы .
.
Находим собственные вектора решения для каждого получившегося r.
r1:
r2:
:
:
.
Определяем константы С.
;
)
С1
=
)
С2=
)
С2=
)
С4=
Разбиваем
фундаментальную матрицу на 4 блока
Проверка:
. Находим
вектора и .
Матрица
К имеет вид:
;
Вычисляем
блоки фундаментальной матрицы:
Соответственно
матрица К численно равна:
Возьмем следующую начальную точку:
Вектор фазового состояния имеет вид:
x(t)= [Φ11 (t, t0) - 2∙Φ12 (t, t0)∙K(tk, t0)]∙x(t0)
Вектор сопряженных функций:
Ψ(t)= [Φ21 (t, t0) - 2∙Φ22 (t, t0) K(tk, t0)]∙x(t0)
Вектор управления:
u(t) = W-1 B T∙[ Φ22 (t, t0) K(tk, t0) - 
Φ21 (t, t0)]∙x(t0), tÎ[t0, tk]
2.
Графики
Векторы фазового состояния:
При
Управление:
Сопряженные векторы:
3. Использование численных методов
3.1 Методическая часть
.1.1 Формирования функции Гамильтона
Для ДС
= f(x, u) , x(t0) = x0,
(2.1)
и
критерия оптимальности
Jобщ = 
f 0(x(t), u(t)) dt + F (x(tk))
(2.2)
функция
Гамильтона будет иметь вид:
где
ψ(t)
- вектор сопряженных функций, размерности
[n´1].
В
частном случае, для ДС (2.1) функция Гамильтона будет следующей.
H(x, u, ψ) = xT(t) Q x(t) + uT(t)
W u(t)+ ψ(t)T∙ [A x(t) + B u(t)] (2.4)
Соответственно,
в общем случае вектор сопряженных функций определяется как результат
интегрирования системы дифференциальных уравнений:
(2.5)
где
F(x(tk)) - терминальная часть критерия оптимальности общего
вида (2.2)
В
частном случае (для линейной ДС (1.1) и критерия (1.2) ), система
дифференциальных уравнений, определяющая сопряженный вектор ψ(t), будет иметь вид:
(2.6)
Из
систем уравнений (2.7) (2.8) видно, что они обе зависят как от самого
сопряженного вектора ψ(t)
, так и от вектора х(t). В связи с этим
ни та, ни другая системы не могут быть проинтегрированы без учета информации о
поведении вектора фазового состояния х(t) в процессе изменения вектора ψ(t). Кроме того, в общем случае для уравнений (2.7) это
невозможно сделать без информации о поведении вектора управления u(t) в
процессе движения системы tÎ[t0, tk]
.
3.1.2 Применение необходимых условий оптимальности в форме
принципа минимума
В соответствии с теорией оптимального управления программная стратегия
оптимального управления должна удовлетворять для ДС (2.1), критерия
оптимальности (2.2) и при ограничениях
u(t) Î Uextr
Uextr = { u(∙)Î Ut: umin ≤ u(t) ≤ umax } (2.7)
следующим необходимым условиям в форме принципа минимума [2].
(2.8)
В
результате проведения процедуры минимизации гамильтониана на всем интервале
времени движения системы tÎ[t0, tk] можно получить
структуру оптимального программного управления:
u((∙), x(∙),
ψ(∙)) = { u((t), x(t), ψ(t)), tÎ[t0, tk]
} (2.9)
Как
видно из (2.9), хотя структура оптимального программного управления в конечном
итоге является функцией времени, однако в общем случае она зависит также от
вектора текущего состояния ДС - х(t) и текущего сопряженного вектора - ψ(t).
В
частном случае, для рассматриваемых ДС (1.1) и критерия (1.2), необходимые
условия в форме принципа минимума могут быть представлены в более простом и
конкретном виде. Это становится возможным поскольку:
во-первых,
в частной постановке задачи (1.1) - (1.2) на управление не накладывается прямых
ограничений, таких как в (2.7), и операция минимизации гамильтониана по
управлению сведется к применению необходимых условий его минимума по
управлению, то есть обнулению частной производной гамильтониана по управлению:
(2.10)
во-вторых,
в силу конкретности вида функции Гамильтона (2.4) и (2.10), структура
оптимального программного управления может быть получена в конкретном виде.
После
взятия частной производной гамильтониана (2.4) по управлению можно получить:
W
u(t) + ψ(t)T∙B
= 0 (2.11)
Откуда
следует, что структура оптимального программного управления для рассматриваемой
частной задачи будет иметь вид:
(∙)
= - W-1 B T∙ψ(∙)
= {u(t) = - W-1 B T∙ψ(t),
tÎ[t0, tk]} (2.12)
Однако,
также как в общем случае (2.9) структура оптимального программного управления
(2.12) зависит от поведения текущего сопряженного вектора - ψ(t),tÎ[t0, tk].
.1.3 Определение оптимального программного управления путем
решения краевой задачи для канонической системы дифференциальных уравнений
В общем случае для задачи (2.1), (2.2), (2.7) каноническая система
дифференциальных уравнений будет иметь вид:
= f(x, u)
, x(t0) = x0
(2.12)
(2.14)
Общая
размерность системы (2.12), (2.14) будет [3n ´ 1].
Принципиально
важной особенностью этой системы является то, что первая система (2.12)
«привязана» к вектору начального состояния x(t0) (на «левом» конце траектории),
а вторая система (2.14), которую необходимо интегрировать совместно с первой,
задана в момент окончания движения ψ(tk) (на «правом» конце траектории).
Учитывая
эту особенность системы (2.12), (2.14), для определения траектории оптимального
движения системы (2.1) - x(∙) и соответствующего сопряженного вектора ψ(∙),
необходимых для окончательного
определения оптимальной программы управления (см. (2.9)), требуется решить
краевую задачу, удовлетворяющую заданным условиям на обоих концах траектории.
В
общем случае для достижения этой цели предлагается на «правом» конце траектории
сформировать функцию невязки, зависящую от неизвестных начальных значений сопряженного
вектора ψ(t0), как от
параметров минимизации функции невязки. Предлагается также, связь между этой
функцией и параметрами установить путем численного интегрирования канонической
системы уравнений (2.12), (2.14), то есть:
Φ(tk , ψ(t0)) = [ψинт(tk
,
ψ(t0)) - ψ(tk)]2 (2.15)
где
Φ(tk , ψ(t0)) - функция
невязки, вычисляемая в конечный момент времени tk ; 
-
сопряженный вектор, получаемый как результат численного интегрирования системы
(2.12), (2.14) при начальных условиях х(t0), ψ(t0) ; ψ (tk) - требуемое
конечное значение сопряженного вектора, вычисляемое согласно (2.14) по формуле:
. (2.16)
Минимизация
невязки (2.15) должна осуществляться итеративно с помощью одного из методов
математического программирования нулевого порядка (при выполнении данной
курсовой работы был использован метод случайного поиска с направляющим конусом
- см. Приложение 1). Очевидно, что для каждой новой итерации необходимо заново
интегрировать уравнения (2.12), (2.14) с новыми значениями 
.
(2.17)
3.2 Блок-схема алгоритма численного решения краевой задачи
Метод деформируемого многогранника
Суть метода заключается в использовании многогранника с n+1 -ой вершиной в n-мерном пространстве области
определения целевой функции для определения в его окрестности либо минимума
функции, либо направления уменьшения функции. Затем, исходя из информации о
значениях функции в вершинах многогранника, построение нового многогранника,
включающего или приближающегося к минимуму функции, путем условной «деформации»
исходного с целью получения новой информации о нахождении минимума функции.
Алгоритм метода
1) Формирование в n-мерном
пространстве многогранника путем задания (n+1)-ой его вершин. В частности, путем сдвига начала координат
на заданную величину h по
всем n координатам.
(5.29)
) Определение среди всех вершин таких вершин, в которых достигается
максимальное и минимальное значения функции.
(5.20)
(5.21)
) Определение «центра тяжести» многогранника.
(5.22)
Из вершин этого многогранника исключается вершина, в которой достигается
максимум (5.20).
) Строится «отражение» вершины с наибольшим значением функции
относительно точки «центр тяжести» - x(n+2).
(5.22)
) Строится точка «растяжение» выбранного направления за точку
«отражение».
(5.24)
) Если f(x(n+4)) < f(x(l)) , тогда
осуществляется замена точки с максимальным значением функции (5.20) на точку x(n+4):
(h) Þ x(n+4) (5.25).
Однако в противном случае проверяется другое условие:
, где j-
одна из вершин многогранника (5.26)
кроме
Если условие (5.26) выполняется, тогда производится замена точки с
максимальным значением функции (5.20) на точку «отражения»:
x(h) Þ x(n+2)
При не выполнении условий (5.25), (5.26) осуществляется «сжатие»
выбранного направления и точка x(h) заменяется на точку x(n+5):
x(h) Þ x(n+5)
, 
(5.27)
Если
f(x(n+5))> f(x(h)) f(x(i)), для 
, то
осуществляется «редукция» многогранника и строится новый многогранник со
сторонами вдвое меньше исходного , т.е.
****>
(5.28)
«Редукция» производится относительно вершины с минимальным значением
целевой функции.
) После формирования нового многогранника осуществляется проверка
условия окончания всей процедуры поиска минимума целевой функции:
(5.29)
Результаты численных исследований
Зададим начальные значения.
Начальная точка:
Шаг:
h=0.1
Точность:
epsnev = 0.1
Точность
численного метода:
eps = 0.1
На основе результатов интегрирования и численного решения краевой задачи
были получены следующие графики:
Вектор фазового состояния:
х1
Рис.1
х2
Рис.2
Управление:
Рис. 2
Графики:
Рис. 1, Рис. 2, Рис. 2 - получены при исходных данных, описанных в разделе
«исходные данные», а также при 
Выводы
В результате выполнения курсовой работы было определено оптимальная
программа управление u(∙)
динамической системой.
Решение данной задачи выполнялось двумя способами:
аналитическое решение;
численное решение.
Для нахождения оптимальной программы управления было определено
оптимальное начальное значение сопряженного вектора Ψ
при заданном начальном
Х:
, 
.
На
основе этих начальных значений было построено оптимальная программа управления u(t).
Из
графиков видно, что существенных различий при получении результата не
наблюдается. Однако, в пользу решения с использованием численного метода
говорит ряд факторов: во-первых, построение фундаментальной матрицы -
трудоёмкий и рутинный процесс, поэтому при его выполнении без аппаратных
средств, возможно сильное влияние «человеческого фактора», в то время как
численный метод надёжней, а машина может вычислять с большей точностью;
во-вторых, варьирование весовой матрицы аналитическим методом весьма сложный
процесс - это обусловлено необходимостью пересчитывать фундаментальную матрицу
каждый раз при изменении элемента весовой матрицы, что достаточно неудобно.
Приложение
#include <iostream>
#include <math.h>
#include <conio.h>
#include <iomanip>
#include <locale.h>namespace std;
#include "Opt.h"
#include "result.txt"dx1, dx2, dpsi1, dpsi2, x01, x02,
psi01, psi02, x1, x2, psi1, psi2;psi, m, n, l;
Matrixd fu (Matrixd const& x,double t)
{dt(4,1);(0)= 2 * x02;(1)= 4 * x02 - 9/4 * psi02;
dt(2)= - 16 * x01;(2)= - 2 * psi01 - 4 * psi02 - 4 * x02;dt;
}<typename T>RungeKuttStep(Matrixd & x,double
t,double h,T f)
{k1(x.numRows(),x.numColumns());=h*f(x,t);k2(x.numRows(),x.numColumns());=h*f(x+0.5*k1,t+0.5*h);k2(x.numRows(),x.numColumns());=h*f(x+0.5*k2,t+0.5*h);k4(x.numRows(),x.numColumns());=h*f(x+k2,t+h);xi(x.numRows(),x.numColumns());
xi=x+1.0/6.0*(k1+2.0*k2+2.0*k2+k4);xi;
}
//---------------------------------------------------------------------------__fastcall
TForm1::Button1Click(TObject *Sender)
{= 0;= 1;= -5;= -20;T;t0 = 0;step = 0.1;Tk = 4.0;EPSnev =
0.1;= 2 * x02 - 9/4 * psi02;= 4 * x02 - 9/4 * psi02;= - 16 * x01 - 2 * psi01;=
- 4 * x02 - 4 * psi02;
Matrixd x;x0(4,1);(0) = dx1;
x0(1) = dx2;(2) = dpsi1;(2) = dpsi2;=x0;(T=t0; T < Tk;
T+=step)
{=RungeKuttStep(x,T,step,fu);
}(fabs(x(2))<= 0.1 && fabs(x(2)) <= EPSnev)
{;
}
{
//---------------------------НАЧАЛО ЦИКЛА МИНИМИЗАЦИИ--------------
do {mnog();Func(double x1, double x2)
{f = 2*pow((x1-2),2)+pow((x2-2),2);f;
}limit1(double x1,double x2)
{g1 = x1+x2-2;g1;
}limit2(double x1)
{g2 = -x1;
return g2;
}limit2(double x2)
{g2 = -x2-2;g2;
}PFunc(double x1,double x2,double a1,double a2,double a2)
{PF =
2*pow((x1-2),2)+pow((x2-2),2)+*pow((x1+x2-2),2)+*pow((-x1),2)+*pow((-x2-2),2);PF;
}main()
{x11,x12,x21,x22,x21,x22,f1,f2,f2,f4,f5,f6,f7;maxf,minf,srf,x1max,x2max,y1min,y2min,z1sr,z2sr;n=2,
i=1;x1[100], x2[100];double e=0.1;double alfa=0.8;double gamma=1.2;double
beta=0.6;fname[999]=FNAME;*in;=fopen(fname,
"w+");(LC_ALL,"Russian");=-1.2; x12=-2.6;();(in);();
}mnog(double x1, double x2)
{x21=rand();double x22=rand();x21=rand; double
x22=rand();{=2*(x11-2)*(x11-2)+(x12-2)*(x12-2);=2*(x21-2)*(x21-2)+(x22-2)*(x22-2);=2*(x21-2)*(x21-2)+(x22-2)*(x22-2);
{(f1>f2)
{=f1; x1max=x11; x2max=x12;
if(f2>f2)
{=f2; y1min=x21; y2min=x22;
srf=f2; z1sr=x21; z2sr=x22;
}
{=f2; y1min=x21; y2min=x22;
srf=f2; z1sr=x21; z2sr=x22;
}
}
{=f2; x1max=x21; x2max=x22;
minf=f2; y1min=x21; y2min=x22;
srf=f1; z1sr=x11; z2sr=x12;
}
}(f1>f2)
{(f1>f2)
{=f1; x1max=x11; x2max=x12;
if(f2>f2)
{=f2; y1min=x21; y2min=x22;
srf=f2; z1sr=x21; z2sr=x22;
}
{=f2; y1min=x21; y2min=x22;
srf=f2; z1sr=x21; z2sr=x22;
}
}
{=f2; x1max=x21; x2max=x22;
minf=f2; y1min=x21; y2min=x22;
srf=f1; z1sr=x11; z2sr=x12;
}
}
{(f2>f2)
{=f2; x1max=x21; x2max=x22;
if(f1>f2)
{=f2; y1min=x21; y2min=x22;
srf=f1; z1sr=x11; z2sr=x12;
}
}
{=f2; x1max=x21; x2max=x22;
minf=f1; y1min=x11; y2min=x12;
srf=f2; z1sr=x21; z2sr=x22;
}
}
/*cout«"maxf="«maxf«" x1="«x1max«"
x2="«x2max«endl
«"minf="«minf«" x1="«y1min«"
x2="«y2min«endl«"srf="
«srf«" x1="«z1sr«" x2="«z2sr«endl;*/
cout«"minf="«minf«" x1="«y1min«"
x2="«y2min«endl;(in, "minf= %f\tx1=
%f\tx2=%f\n",minf,y1min,y2min);
//Центр тяжести[4]=(y1min+z1sr)/n;[4]=(y2min+z2sr)/n;=2*(x1[4]-2)*(x1[4]-2)+(x2[4]-2)*(x2[4]-2);
//cout«"x1(n+2)="«x1[4]«"
x2(n+2)="«x2[4]«endl;
//Точка отражения[5]=x1[4]+alfa*(x1[4]-x1max);
x2[5]=x2[4]+alfa*(x2[4]-x2max);
f5=2*(x1[5]-2)*(x1[5]-2)+(x2[5]-2)*(x2[5]-2);
//cout«"x1(n+2)="«x1[5]«" x2(n+2)="«x2[5]«endl;
//Точка растяжения
x1[6]=x1[4]+gamma*(x1[5]-x1[4]);[6]=x2[4]+gamma*(x2[5]-x2[4]);
f6=2*(x1[6]-2)*(x1[6]-2)+(x2[6]-2)*(x2[6]-2);
//cout«"x1(n+4)="«x1[6]«"
x2(n+4)="«x2[6]«endl;(f6<maxf)
{max=x1[6];max=x2[6];
}
{(f5<=srf)
{max=x1[5];max=x2[5];
}
{
//Сжатие[7]=x1[4]-beta*(x1[4]-x1max);
x2[7]=x2[4]-beta*(x2[4]-x2max);
x1max=x1[7];max=x2[7];=2*(x1[7]-2)*(x1[7]-2)+(x2[7]-2)*(x2[7]-2);(f7>srf)
{
//Редукцияsr=y1min-0.5*(z1sr-y1min);sr=y2min-0.5*(z2sr-y2min);
}
}
}=0;x11=x1max;x1max=0;=0;x12=x2max;x2max=0;=0;x21=y1min;y1min=0;=0;x22=y2min;y2min=0;
double PFunc();
}(fabs(minf-maxf)>e);
}
//-------------------------------КОНЕЦ ЦИКЛА МИНИМИЗАЦИИ-------------=
0;= 1;
psi01 = uX[0];= uX[1];= 2 * x02 - 9/4 * psi02;= 4 * x02 - 9/4
* psi02;= - 16 * x01 - 2 * psi01;= - 4 * x02 - 4 * psi02;x;x0(4,1);(0) = dx1;
x0(1) = dx2;(2) = dpsi1;(2) = dpsi2;=x0;(T=t0; T < Tk;
T+=step)
{=RungeKuttStep(x,T,step,fu);
}
}while (fabs(x(2)) > EPSnev && fabs(x(2)) >
EPSnev)