Алгоритм реализации схемы Беллмана

ОТЧЕТ

По лабораторной работе №7

по курсу Теория оптимального управления

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

ОГУ 010200.6006.02 О

Руководитель

_____________ Иванова Ю.П.

“___ ”______________2011 г.

Исполнитель

студент гр. 08 ПриМ

_______________Кияева Е. А.

“___ ” _____________2011 г.

Оренбург 2011

Содержание

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

Алгоритм реализации схемы Беллмана….…………….………….……7

Практическая часть……………………………………………………….8

Приложение А – Код программы…….……………………………….11

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

Рассмотрим модель очистки воды от загрязнения органическими отходами.

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

(1)

(2)

(3)

(4)

(5)

(6)

Здесь (мг. /ед.врем. л) – количество отходов, удаленных из единицы объема в единицу времени;

(мг. /л) – концентрация органических отходов в момент времени t;

( мг. /л.) – дефицит кислорода в воде в момент времени t;

а (ед.врем.) – коэффициент, характеризующий стоимость работ в единицу времени;

(1/ ед.врем.) – коэффициент отбора кислорода за единицу времени;

(1/ед.врем.) – коэффициент реаэрации в единицу времени;

– начальное значение в момент дефицита кислорода и концентрации отходов соответственно;

В уравнении (2) производная (мг. /ед.врем. л) характеризует скорость разложения отходов в единицу времени.

В уравнении (3) (мг. /ед.врем. л) характеризует динамику дефицита кислорода в единицу времени.

В неравенстве (мг. /ед.врем. л) – максимальное количество отходов, удаленных из единицы объема в единицу времени.

Целью управления является минимизация суммарных затрат на проведение очистных работ, уменьшение концентрации отходов и дефицита кислорода в течении рассматриваемого периода времени Т.

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

Применим для решения задачи (1) – (6) схему Беллмана. Для определенности будем считать, что в задаче и заданы параллельные фазовые ограничения.

Алгоритм реализации схемы Беллмана

Шаг	Действие
	Задать разбиения: а) отрезка времени : б) – ограничения по в) – ограничения по г) – ограничения по u ;
Первый этап (по убывающему индексу)
	Вычислить значения функции Беллмана в точках терминального множества , где
	Начало цикла по ;
	Положить – число проколотых точек. (Точка называется проколотой, если она не принадлежит множеству достижимости) Организовать цикл по точкам множества ;
	Вычислить координаты точки
	Положить – число недопустимых управлений в точке ;
	Организовать цикл по допустимым управлениям ;
	Вычислить управления ;
	Вычислить точку ;
	Проверить выполнение фазового ограничения (принадлежит ли точка множеству ):
	Найти ближайшую к непроколотую точку сетки ; если все ближайшие точки сетки (т.е. точки, находящиеся на расстоянии меньше чем ) являются проколотыми, то и переход к шагу 14;
	Вычислить ;
	Найти ;
	Положить ;
	Перейти к следующему шагу цикла по n, т.е. положить и проверить условие окончания цикла: ;
	Проверяем, не является ли точка проколотой:
	Перейти к следующей итерации цикла по i,j, т.е. а) положить и проверить условие б) положить и проверить условие
	Проверить, все ли точки сетки являются непроколотыми:
	Перейти к следующей итерации цикла по , т.е. полагаем и проверить условие окончания цикла:
	Получить функцию Беллмана и − синтез управления;
Второй этап
	Найти ближайшую к непроколотую точку сетки . Если все ближайшие точки сетки (т.е. точки, находящиеся на расстоянии меньше, чем ) являются проколотыми, то “нет решений” и идти к шагу 32, иначе − к шагу 21;
	Вычислить ;
Третий этап (по возрастающему индексу) Определение оптимальной траектории и программного управления
	Положить ;
	Организовать цикл по полагаем ;
	Найти ближайшую к непроколотую точку сетки . Если все ближайшие точки сетки (т.е. точки, находящиеся на расстоянии меньше, чем ) являются проколотыми, то “нет решений” и надо перейти к шагу 32, иначе − к шагу 25;
	Определить соответствующее ей управление – синтез и полагаем ;
	Вычислить точку ;
	Перейти к следующему шагу цикла по , т.е. положить и проверить условие
	Получить приближение оптимальной траектории и программного управления, вычислить искомое значение функционала I;
	Проверить точность полученного решения:
	Для получения более точного решения увеличиваем число точек разбиения множеств . Перейти к шагу 0;
	− решение задачи. Конец алгоритма.

Практическая часть

1) Дискретная аппроксимация.

Разобьем отрезок [0,T] на N частей точками и ,приняв их в качестве узловых, интеграл (1) заменим квадратурной формулой прямоугольников. Дифференциальные уравнения (2), (3) - разностным уравнением с помощью схемы Эйлера. В результате придем к дискретной задаче оптимального управления вида:

2) Результаты выполнения программы

Концентрация органических отходов и дефицит кислорода в воде с течением времени уменьшаются. Количество отходов, удаленных из единицы объема в единицу времени, остается постоянным.

Таким образом, при увеличении продолжительности очистки воды количество отходов и дефицит кислорода в воде будут уменьшаться.

Приложение А

(обязательное)

Код программы

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm1 *Form1;

struct point

{bool f;double x; double y;};

double dt, dx1,dx2,du,eps,J,J_,k1=0.0004321,k2=0.008281,a,T;

int N,m1,m2,m3;

typedef double **massiv;

typedef double *masv;

struct sov

{massiv u,V;};

typedef sov *m;

double Bk(double x, double y, double u)

{

return (a*u+x+y)*dt;

}

double f(double x1,double x2, double u,int i)

{

switch (i)

{

case 0: return x1-(k1*x1+u)*dt;

case 1: return x2+(k1*x1-k2*x2)*dt;

}

}

double I(massiv x, masv u)

{

double s=0;

for (int i=0;i<N;i++)

s+=(a*u[i]+x[0][i]+x[1][i])*dt;

return s;

}

bool D(point x)

{

if ((x.x>=1)&&(x.x<=10))

if ((x.y>=-10)&&(x.y<=10)) return true;

else return false;

else return false;

}

int round (double a)

{

if (a-int (a)>0.5) return int(a)+1;

if (a-int(a)<=0.5) return int (a);

}

point prok(point x)

{

point a;

double d;

if (dx1>dx2) d=dx1;

else d=dx2;

a.f=false;

int i0,j0;

a.x=round((x.x-1)/dx1);

a.y=round((x.y+10)/dx2);

if ((1+a.x*dx1>=1)&&(1+a.x*dx1<=10)&&(-10+a.y*dx2>=-10)&&(-10+a.y*dx2<=10)) a.f=true;

else

if (((x.x+d>=1)&&(x.x<1))||((x.x-d>10)&&(x.x>10))||((x.y+d>=-10)&&(x.y<-10))||((x.y-d<=10)&&(x.y>10)))

{

a.f=true;

if (x.x<1) a.x=0;

if (x.x>10) a.x=m1;

if (x.y<-10) a.y=0;

if (x.y>10) a.y=m2;

}

return a;

}

double max(double a, double b)

{

if (a>b) return a;

else return b;

}

//---------------------------------------------------------------------------

void __fastcall TForm1::Button1Click(TObject *Sender)

{

int i0,j0,i_0,j_0;

eps=StrToFloat(Edit6->Text);

N=StrToInt(Edit4->Text);

m1=StrToInt(Edit1->Text);

m2=StrToInt(Edit2->Text);

m3=StrToInt(Edit3->Text);

a=StrToFloat(Edit7->Text);

T=StrToFloat(Edit12->Text);

bool flag=false;

int notx,notu;

point xk,x1;

double up;

mas U;

mass X;

m x;

do

{

dt=T/N;

dx1=9.0/m1;

dx2=20.0/m2;

du=1.0/m3;

x=new sov [N+1];

X=new double* [2];

U=new double [N+2];

for (int i=0;i<=N;i++)

{

x[i].u=new double *[m1+1];

x[i].V=new double *[m1+1];

for (int j=0;j<=m1;j++)

{

x[i].u[j]=new double [m2+1];

x[i].V[j]=new double [m2+1];

}

}

for (int i=0;i<2;i++)

X[i]=new double [N+2];

for (int i=0;i<=m1;i++)

for (int j=0;j<=m2;j++)

x[N].V[i][j]=0;

for (int k=N-1;k>=0;k--)

{

notx=0;

for (int i=0;i<=m1;i++)

for (int j=0;j<=m2;j++)

{

x1.x=1+i*dx1;

x1.y=-10+j*dx2;

notu=0;

for (int n=0;n<=m3;n++)

{

up=n*du;

xk.x=f(x1.x,x1.y,up,0);

xk.y=f(x1.x,x1.y,up,1);

if (!D(xk)) notu+=1;

else

{

if (prok(xk).f)

{

i0=prok(xk).x;

j0=prok(xk).y;

double fn=Bk(x1.x,x1.y,up)+x[k+1].V[i0][j0];

if ((fn<x[k].V[i][j]) ||(n==0))

{

x[k].V[i][j]=fn;

x[k].u[i][j]=up;

}

}

else notu++;

}

}

if (notu==m3) notx+=1;

}

}

point b;

b.x=X[0][0]=5;

b.y=X[1][0]=2;

int i_,j_;

if (prok(b).f)

{

i_=prok(b).x;

j_=prok(b).y;

}

b.x=X[0][0]=5;

b.y=X[1][0]=2;

for (int k=0;k<N;k++)

{

if (prok(b).f)

{

i0=prok(b).x;

j0=prok(b).y;

U[k]=x[k].u[i0][j0];

b.x=X[0][k+1]=f(X[0][k],X[1][k],U[k],0);

b.y=X[1][k+1]=f(X[0][k],X[1][k],U[k],1);

}

}

J=I(X,U);

flag=true;

if ((fabs(J-x[0].V[i_][j_])<=eps)||((J-J_)<=eps)) flag=false;

else

{

J_=J;

m1=2*m1;

m2=2*m2;

m3=2*m3;

N=2*N;

}

}

while (flag);

Edit5->Text=FloatToStr(J);

StringGrid1->Cells[0][0]="dt";

StringGrid1->Cells[1][0]="x1(t)";

StringGrid1->Cells[2][0]="x2(t)";

StringGrid1->Cells[3][0]="u(t)";

for (int i=0;i<=N;i++)

{

StringGrid1->Cells[0][i+1]=FloatToStrF(i*dt,ffGeneral,5,2);

StringGrid1->Cells[1][i+1]=FloatToStrF(X[0][i],ffGeneral,5,2);

StringGrid1->Cells[2][i+1]=FloatToStrF(X[1][i],ffGeneral,5,2);

if (i%10==0)

{

Chart1->Series[0]->AddXY(StrToFloat(i*dt),StrToFloat(X[0][i]),"",clRed);

Chart1->Series[1]->AddXY(StrToFloat(i*dt),StrToFloat(X[1][i]),"",clBlue);

}

if (i!=N)

{

StringGrid1->Cells[3][i+1]=FloatToStrF(U[i],ffGeneral,5,2);

if (i%10==0)

{

Chart2->Series[0]->AddXY(StrToFloat(i*dt),StrToFloat(U[i]),"",clRed);

}

}

}

delete []X;

delete [] x;

delete []U;

}