Оптимальная программа управления динамической системой

Постановка задачи Рассматривается динамическая система:

где x — вектор фазового состояния ДС, размерности; u — вектор управления ДС, размерности; А, В — матрицы постоянных коэффициентов системы, размерности и соответственно.

Требуется определить оптимальную программу управления системой ДС, обеспечивающую минимум квадратичного критерия:

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

Исходные данные программа управление динамический уравнение

1. Аналитическое решение

1.1 Методика решения поставленной задачи

1. Построим функцию Гамильтона. Для смешанного критерия типа (2) функция Гамильтона имеет вид:

(4)

Здесь — компонента сопряженного вектора, соответствующая фиксированной координате .

Имея в виду, что — Сonst по времени и в конечный момент времени функция Гамильтона (4) может быть переписана в следующем виде:

2. В данном случае Н.У. имеет частный вид:

Поскольку прямых ограничений на управление в задаче не вводится Отсюда следует, что:

(5)

2. Решение канонической системы уравнений. Необходимость в этом решении возникает в связи с необходимостью получения сопряженного вектора .

(6)

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

Если ввести новые обозначения:

Тогда каноническая система (6) примет вид:

(7)

4. Используя свойство фундаментальной матрицы можно записать для системы (7) следующее решение:

В развернутом виде можно записать:

где, , , — соответствующие блоки фундаментальной матрицы .

Поскольку конечное условие сопряженного вектора можно записать следующее соотношение:

Из полученного выражения можно определить вектор состояний сопряженных переменных:

Теперь векторы и могут быть поставлены в зависимость только от начального состояния вектора. А именно:

где Таким образом, в соответствии с оптимальной структурой управления (5) можно получить искомую программу оптимального управления в виде:

1.2 Процедура решения поставленной задачи

1. Составляем матрицу расширенной системы .

2. Находим собственные числа матрицы .

2. Находим собственные вектора решения для каждого получившегося r.

r1:

r2:

r2:

r4:

4. Определяем константы С.

;

1)

С1 =

2)

С2=

2)

С2=

4)

С4=

Разбиваем фундаментальную матрицу на 4 блока Проверка:

6. Находим вектора и .

Матрица К имеет вид:

;

Вычисляем блоки фундаментальной матрицы:

Соответственно матрица К численно равна:

Возьмем следующую начальную точку:

Вектор фазового состояния имеет вид:

x (t)= [Ц₁₁ (t, t₀) — 2•Ц₁₂ (t, t₀)•K (t_k, t₀)]•x (t₀)

Вектор сопряженных функций:

Ш (t)= [Ц₂₁ (t, t₀) — 2•Ц₂₂ (t, t₀) K (t_k, t₀)]•x (t₀)

Вектор управления:

u (t) = W^-1 B ^T?[ Ц₂₂ (t, t₀) K (t_k, t₀) — Ц₂₁ (t, t₀)]•x (t₀), t[t₀, t_k]

2. Графики Векторы фазового состояния:

При

Управление:

Сопряженные векторы:

3. Использование численных методов

3.1 Методическая часть

3.1.1 Формирования функции Гамильтона

Для ДС

= f (x, u), x (t₀) = x₀, (2.1)

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

J_общ = f ⁰(x (t), u (t)) dt + F (x (t_k)) (2.2)

функция Гамильтона будет иметь вид:

H (x, u, ш) = f ⁰(x (t), u (t)) + ш (t)^T• f (x (t), u (t)) (2.2)

где ш (t) — вектор сопряженных функций, размерности [n1].

В частном случае, для ДС (2.1) функция Гамильтона будет следующей.

H (x, u, ш) = x^T(t) Q x (t) + u^T(t) W u (t)+ ш (t)^T• [A x (t) + B u (t)] (2.4)

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

(2.5)

где F (x (t_k)) — терминальная часть критерия оптимальности общего вида (2.2)

В частном случае (для линейной ДС (1.1) и критерия (1.2)), система дифференциальных уравнений, определяющая сопряженный вектор ш (t), будет иметь вид:

(2.6)

Из систем уравнений (2.7) (2.8) видно, что они обе зависят как от самого сопряженного вектора ш (t), так и от вектора х (t). В связи с этим ни та, ни другая системы не могут быть проинтегрированы без учета информации о поведении вектора фазового состояния х (t) в процессе изменения вектора ш (t). Кроме того, в общем случае для уравнений (2.7) это невозможно сделать без информации о поведении вектора управления u (t) в процессе движения системы t[t₀, t_k] .

3.1.2 Применение необходимых условий оптимальности в форме принципа минимума

В соответствии с теорией оптимального управления программная стратегия оптимального управления должна удовлетворять для ДС (2.1), критерия оптимальности (2.2) и при ограничениях

u (t) U_extr

U_extr = { u (•) U_t: u_min? u (t)? u_max } (2.7)

следующим необходимым условиям в форме принципа минимума.

(2.8)

В результате проведения процедуры минимизации гамильтониана на всем интервале времени движения системы t[t0, tk] можно получить структуру оптимального программного управления:

u ((•), x (•), ш (•)) = { u ((t), x (t), ш (t)), t[t₀, t_k] } (2.9)

Как видно из (2.9), хотя структура оптимального программного управления в конечном итоге является функцией времени, однако в общем случае она зависит также от вектора текущего состояния ДС — х (t) и текущего сопряженного вектора — ш (t).

В частном случае, для рассматриваемых ДС (1.1) и критерия (1.2), необходимые условия в форме принципа минимума могут быть представлены в более простом и конкретном виде. Это становится возможным поскольку:

во-первых, в частной постановке задачи (1.1) — (1.2) на управление не накладывается прямых ограничений, таких как в (2.7), и операция минимизации гамильтониана по управлению сведется к применению необходимых условий его минимума по управлению, то есть обнулению частной производной гамильтониана по управлению:

(2.10)

во-вторых, в силу конкретности вида функции Гамильтона (2.4) и (2.10), структура оптимального программного управления может быть получена в конкретном виде.

После взятия частной производной гамильтониана (2.4) по управлению можно получить:

3 W u (t) + ш (t)^T•B = 0 (2.11)

Откуда следует, что структура оптимального программного управления для рассматриваемой частной задачи будет иметь вид:

u (•) = - W^-1 B ^T?ш (?) = {u (t) = - W^-1 B ^T?ш (t), t[t₀, t_k]} (2.12)

Однако, также как в общем случае (2.9) структура оптимального программного управления (2.12) зависит от поведения текущего сопряженного вектора — ш (t), t[t0, tk].

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

В общем случае для задачи (2.1), (2.2), (2.7) каноническая система дифференциальных уравнений будет иметь вид:

= f (x, u), x (t₀) = x₀ (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)), требуется решить краевую задачу, удовлетворяющую заданным условиям на обоих концах траектории.

В общем случае для достижения этой цели предлагается на «правом» конце траектории сформировать функцию невязки, зависящую от неизвестных начальных значений сопряженного вектора ш (t₀), как от параметров минимизации функции невязки. Предлагается также, связь между этой функцией и параметрами установить путем численного интегрирования канонической системы уравнений (2.12), (2.14), то есть:

Ц (t_k, ш (t₀)) = [ш_инт(t_k, ш (t₀)) — ш (t_k)]² (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)

2) Определение среди всех вершин таких вершин, в которых достигается максимальное и минимальное значения функции.

(5.20)

(5.21)

2) Определение «центра тяжести» многогранника.

(5.22)

Из вершин этого многогранника исключается вершина, в которой достигается максимум (5.20).

4) Строится «отражение» вершины с наибольшим значением функции относительно точки «центр тяжести» — x⁽ⁿ⁺²⁾.

(5.22)

5) Строится точка «растяжение» выбранного направления за точку «отражение».

(5.24)

6) Если f (x⁽ⁿ⁺⁴⁾) < f (x^(l)), тогда осуществляется замена точки с максимальным значением функции (5.20) на точку x⁽ⁿ⁺⁴⁾:

x^(h) x⁽ⁿ⁺⁴⁾ (5.25).

Однако в противном случае проверяется другое условие:

где jодна из вершин многогранника (5.26)

кроме

Если условие (5.26) выполняется, тогда производится замена точки с максимальным значением функции (5.20) на точку «отражения»:

x^(h) x⁽ⁿ⁺²⁾

При не выполнении условий (5.25), (5.26) осуществляется «сжатие» выбранного направления и точка x^(h) заменяется на точку x⁽ⁿ⁺⁵⁾:

x^(h) x⁽ⁿ⁺⁵⁾

(5.27)

Если f (x⁽ⁿ⁺⁵⁾)> f (x^(h)) f (x⁽ⁱ⁾), для, то осуществляется «редукция» многогранника и строится новый многогранник со сторонами вдвое меньше исходного, т. е.

****> (5.28)

«Редукция» производится относительно вершины с минимальным значением целевой функции.

7) После формирования нового многогранника осуществляется проверка условия окончания всей процедуры поиска минимума целевой функции:

(5.29)

Если модуль разницы меньше заданного малого >0, принимается решение об окончании процедуры. В противном случае происходит переход на пункт 2) алгоритма.

Результаты численных исследований Зададим начальные значения.

Начальная точка:

Шаг:

h=0.1

Точность:

epsnev = 0.1

Точность численного метода:

eps = 0.1

На основе результатов интегрирования и численного решения краевой задачи были получены следующие графики:

Вектор фазового состояния:

х1

Рис.1

х2

Рис.2

Управление:

Рис. 2

Графики: Рис. 1, Рис. 2, Рис. 2 — получены при исходных данных, описанных в разделе «исходные данные», а также при

Выводы

В результате выполнения курсовой работы было определено оптимальная программа управление u (•) динамической системой.

Решение данной задачи выполнялось двумя способами:

— аналитическое решение;

— численное решение.

Для нахождения оптимальной программы управления было определено оптимальное начальное значение сопряженного вектора Ш при заданном начальном Х:

.

На основе этих начальных значений было построено оптимальная программа управления u (t).

Из графиков видно, что существенных различий при получении результата не наблюдается. Однако, в пользу решения с использованием численного метода говорит ряд факторов: во-первых, построение фундаментальной матрицы — трудоёмкий и рутинный процесс, поэтому при его выполнении без аппаратных средств, возможно сильное влияние «человеческого фактора», в то время как численный метод надёжней, а машина может вычислять с большей точностью; во-вторых, варьирование весовой матрицы аналитическим методом весьма сложный процесс — это обусловлено необходимостью пересчитывать фундаментальную матрицу каждый раз при изменении элемента весовой матрицы, что достаточно неудобно.

Приложение

#include

#include

#include

#include

#include

using namespace std;

#include «Opt.h»

#include «result.txt»

float dx1, dx2, dpsi1, dpsi2, x01, x02, psi01, psi02, x1, x2, psi1, psi2;

float psi, m, n, l;

Matrixd fu (Matrixd const& x, double t)

{

Matrixd dt (4,1);

dt (0)= 2 * x02;

dt (1)= 4 * x02 — 9/4 * psi02;

dt (2)= - 16 * x01;

dt (2)= - 2 * psi01 — 4 * psi02 — 4 * x02;

return dt;

}

template

Matrixd RungeKuttStep (Matrixd & x, double t, double h, T f)

{

Matrixd k1(x.numRows (), x. numColumns ());

k1=h*f (x, t);

Matrixd k2(x.numRows (), x. numColumns ());

k2=h*f (x+0.5*k1,t+0.5*h);

Matrixd k2(x.numRows (), x. numColumns ());

k2=h*f (x+0.5*k2,t+0.5*h);

Matrixd k4(x.numRows (), x. numColumns ());

k4=h*f (x+k2,t+h);

Matrixd xi (x.numRows (), x. numColumns ());

xi=x+1.0/6.0*(k1+2.0*k2+2.0*k2+k4);

return xi;

}

//—————————————————————————————————————;

void __fastcall TForm1: Button1Click (TObject *Sender)

{

x01 = 0;

x02 = 1;

psi01 = -5;

psi02 = -20;

float T;

float t0 = 0;

float step = 0.1;

float Tk = 4.0;

float EPSnev = 0.1;

dx1 = 2 * x02 — 9/4 * psi02;

dx2 = 4 * x02 — 9/4 * psi02;

dpsi1 = - 16 * x01 — 2 * psi01;

dpsi2 = - 4 * x02 — 4 * psi02;

Matrixd x;

Matrixd x0(4,1);

x0(0) = dx1;

x0(1) = dx2;

x0(2) = dpsi1;

x0(2) = dpsi2;

x=x0;

for (T=t0; T < Tk; T+=step)

{

x=RungeKuttStep (x, T, step, fu);

}

if (fabs (x (2))<= 0.1 && fabs (x (2)) <= EPSnev)

{

break;

}

else

{

//—————————————-НАЧАЛО ЦИКЛА МИНИМИЗАЦИИ———————

do {

double mnog ();

double Func (double x1, double x2)

{

double f = 2*pow ((x1−2), 2)+pow ((x2−2), 2);

return f;

}

double limit1(double x1, double x2)

{

double g1 = x1+x2−2;

return g1;

}

double limit2(double x1)

{

double g2 = -x1;

return g2;

}

double limit2(double x2)

{

double g2 = -x2−2;

return g2;

}

double PFunc (double x1, double x2, double a1, double a2, double a2)

{

double PF = 2*pow ((x1−2), 2)+pow ((x2−2), 2)+

a1*pow ((x1+x2−2), 2)+

a2*pow ((-x1), 2)+

a2*pow ((-x2−2), 2);

return PF;

}

void main ()

{

double x11, x12,x21,x22,x21,x22,f1,f2,f2,f4,f5,f6,f7;

double maxf, minf, srf, x1max, x2max, y1min, y2min, z1sr, z2sr;

int n=2, i=1;

double x1[100], x2[100];

const double e=0.1;

const double alfa=0.8;

const double gamma=1.2;

const double beta=0.6;

char fname[999]=FNAME;

FILE *in;

in=fopen (fname, «w+»);

setlocale (LC_ALL," Russian");

x11=-1.2; x12=-2.6;

mnog ();

fclose (in);

getch ();

}

double mnog (double x1, double x2)

{

double x21=rand ();double x22=rand ();

double x21=rand; double x22=rand ();

do {

f1=2*(x11−2)*(x11−2)+(x12−2)*(x12−2);

f2=2*(x21−2)*(x21−2)+(x22−2)*(x22−2);

f2=2*(x21−2)*(x21−2)+(x22−2)*(x22−2);

//cout"" f1=""f1″" f2=""f2″" f2=""f2"endl;

if (f1>f2)

{

if (f1>f2)

{

maxf=f1; x1max=x11; x2max=x12;

if (f2>f2)

{

minf=f2; y1min=x21; y2min=x22;

srf=f2; z1sr=x21; z2sr=x22;

}

else

{

minf=f2; y1min=x21; y2min=x22;

srf=f2; z1sr=x21; z2sr=x22;

}

}

else

{

maxf=f2; x1max=x21; x2max=x22;

minf=f2; y1min=x21; y2min=x22;

srf=f1; z1sr=x11; z2sr=x12;

}

}

if (f1>f2)

{

if (f1>f2)

{

maxf=f1; x1max=x11; x2max=x12;

if (f2>f2)

{

minf=f2; y1min=x21; y2min=x22;

srf=f2; z1sr=x21; z2sr=x22;

}

else

{

minf=f2; y1min=x21; y2min=x22;

srf=f2; z1sr=x21; z2sr=x22;

}

}

else

{

maxf=f2; x1max=x21; x2max=x22;

minf=f2; y1min=x21; y2min=x22;

srf=f1; z1sr=x11; z2sr=x12;

}

}

else

{

if (f2>f2)

{

maxf=f2; x1max=x21; x2max=x22;

if (f1>f2)

{

minf=f2; y1min=x21; y2min=x22;

srf=f1; z1sr=x11; z2sr=x12;

}

}

else

{

maxf=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;

fprintf (in, «minf= %ftx1= %ftx2=%fn», minf, y1min, y2min);

//Центр тяжести

x1[4]=(y1min+z1sr)/n;

x2[4]=(y2min+z2sr)/n;

f4=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;

//Точка отражения

x1[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]);

x2[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;

if (f6

{

x1max=x1[6];

x2max=x2[6];

}

else

{

if (f5<=srf)

{

x1max=x1[5];

x2max=x2[5];

}

else

{

//Сжатие

x1[7]=x1[4]-beta*(x1[4]-x1max);

x2[7]=x2[4]-beta*(x2[4]-x2max);

x1max=x1[7];

x2max=x2[7];

f7=2*(x1[7]-2)*(x1[7]-2)+(x2[7]-2)*(x2[7]-2);

if (f7>srf)

{

//Редукция

z1sr=y1min-0.5*(z1sr-y1min);

z1sr=y2min-0.5*(z2sr-y2min);

}

}

}

x11=0;x11=x1max;x1max=0;

x12=0;x12=x2max;x2max=0;

x21=0;x21=y1min;y1min=0;

x22=0;x22=y2min;y2min=0;

x21=0;x21=z1sr;z1sr=0;

x22=0;x22=z2sr;z2sr=0;

i++;

double PFunc ();

}

while (fabs (minf-maxf)>e);

}

//———————————————-КОНЕЦ ЦИКЛА МИНИМИЗАЦИИ——————;

x01 = 0;

x02 = 1;

psi01 = uX[0];

psi02 = uX[1];

dx1 = 2 * x02 — 9/4 * psi02;

dx2 = 4 * x02 — 9/4 * psi02;

dpsi1 = - 16 * x01 — 2 * psi01;

dpsi2 = - 4 * x02 — 4 * psi02;

Matrixd x;

Matrixd x0(4,1);

x0(0) = dx1;

x0(1) = dx2;

x0(2) = dpsi1;

x0(2) = dpsi2;

x=x0;

for (T=t0; T < Tk; T+=step)

{

x=RungeKuttStep (x, T, step, fu);

}

}while (fabs (x (2)) > EPSnev && fabs (x (2)) > EPSnev)