Размещено на http://www.allbest.ru/

Постановка задачи

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

,

где 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].

(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(n+2)

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

x(h) x(n+5)

, (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 <iostream>

#include <math.h>

#include <conio.h>

#include <iomanip>

#include <locale.h>

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 <typename T>

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= %f\tx1= %f\tx2=%f\n",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<maxf)

{

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)

Размещено на Allbest.ru