3.2 Знаходження оптимального закону керування

керування лінійний модель динамічний

Метод динамічного програмування

Принцип оптимальності. Оптимальне керування володіє тою властивістю, що, яким би не був початковий стан системи й початкове рішення (тобто оптимальне керування, отримане для інтервалу ), то наступне рішення (керування в інтервалі ) повинне бути оптимальним щодо стану, що виник у результаті першого рішення.

Принцип оптимальності стверджує, що вибір оптимального керування визначається тільки станом системи в даний момент часу і призначеним критерієм оптимальності, та не залежить від того, яким чином система прийшла в даний стан. Принцип оптимальності покладений в основу методу динамічного програмування.

Безперервний варіант динамічного програмування. Нехай потрібно мінімізувати функціонал:

 (12)

де - не фіксовано;

.

Одержуємо для цієї задачі необхідну умову оптимальності, заснована на принципі оптимальності. Припустимо, оптимальне керування й фазова траєкторія знайдені. Позначимо

Візьмемо на траєкторії проміжну крапку . Ділянка траєкторії від точки до точки відповідно до принципу оптимальності є оптимальною траєкторією.

Маємо .

Представимо інтеграл у вигляді суми двох інтегралів наступного вигляду:

З огляду на малість , перший інтеграл правої частини можна замінити:

де - величина високого порядку малості.

Тоді

Тепер припустимо, що уведена функція S безперервна й має безперервні часткові похідні по всіх аргументах, тобто існує

Зроблене припущення дозволяє функцію розкласти на ряд Тейлора:

(14)

Підставимо рівняння (14) у рівняння (13). Функція не залежить від керування в момент часу t, тому її можна винести за знак мінімуму. У результаті з рівняння (13) з урахуванням рівняння (14) і

(15)

Рівняння (15) являє собою функціональне рівняння Беллмана із граничною умовою . Щоб знайти керування, що надає мінімум вихідному функціоналові (1.58), прирівняємо до нуля частинну похідну по лівій частині (15). Тоді

(16)

Рівняння (16)разом з рівнянням

утворить систему функціональних рівнянь Беллмана, використовувану для визначення оптимального керування:

Розвязання

Для обєкту, математична модель якого описується диференціальним рівнянням:

знайти за методом динамічного програмування таке управління u(t), що заданий критерій якості I(x(t),u(t)) виду:

за допомогою методу динамічного програмування приймає найменше значення.

Перетворимо диференціальне рівняння 3 го порядку в систему диференційних рівнянь:

Критерій якості набуває вигляду:

,

де

- є підінтегральна функція, а функції мають зміст:

Обираємо функцію Беллмана у вигляді квадратичної форми:

Знаходимо частинні похідні:

Формуємо систему рівнянь Беллмана згідно з отриманими частковими похідними:

або

Знайдемо часткові похідні функції Беллмана:

Тоді перше рівняння системи запишеться:

після групування отримаємо:

Одержана тотожність буде виконуватись тільки в тому випадку, якщо всі коефіцієнти одночасно дорівнюватимуть нулю:

Знайдені рішення системи алгебраїчних рівнянь запишемо в таблиці, а ті розвязки, що нам підходять виділимо.

,

Отже одержимо функцію оптимального керування:

Модель зовнішніх збурень:

Виконавчий привід:

Передаточна функція прямих ланок:

Передаточна функція замкнутої системи:

Передаточна функція розімкнутої системи:

Передаточна за зовнішнім збуренням:

3.3 Програмний код

function dkr_tf

clc

clear

n=9;

b = n/(n+4);

d =30;

A_33 = 0.0143 ;

A_13 = 0.0496 ;

A_23 = 0.0675 ;

k_1 = ((d*A_13)/(2*b))

k_2 = ((d*A_23)/(2*b))

k_3 = ((d*A_33)/(b))

%u_opt = x*k1 + dx*k2 + ddx*k3

n= 9;

a_2 = 3*n+n/10;

a_1 = n/2+(3*n^2)/10;

a_0 = (n^2)/20;

Tvp=0.01;

kvp = 0.7 + 0.1*n;

Wvp = tf([kvp],[Tvp 1]);

T1 = 0.15;

Wd = tf([T1 0],[1]);

Wok = tf([d],[1 a_2 a_1 a_0])

%figure(11), step(Wok)

Wos = tf([k_3 k_2 k_1],[1])

w1 = series(Wvp,Wok)

Ws_c = feedback(w1, Wos) %loop system

Ws_u = series(w1, Wos) % break system

Ws_e = 1/(1+Ws_u) % loop system of error

Wf_c=(Wok/(1+Ws_u))*Wd

%figure(12),step(Ws_c)

m1= 0;

m2 =0;

m3 =0;

mn=0;

while (mn ~= 4)

mn = menu(Дослідження ДС,...

Переходные характеристики,...

Логарифмические характеристики,...

Інші,...

Вихід);

switch mn

case 1

while (m1~=5)

m1 = menu(Переходные характеристики,...

Перехідна характеристика замкнутої системи,...

Перехідна характеристика розімкнутої системи,...

Перехідна характеристика по помилці рег.,...

Перехідна характеристика по зовнішньому впливу.,...

Вихід);

switch m1

case 1

name=Перехідна характеристика замкнутої системи;

step(Ws_c, name)

case 2

name=Перехідна характеристика розімкнутої системи;

stepview(Ws_u, name)

case 3

name=Перехідна характеристика по помилці рег.;

stepview(Ws_e, name)

case 4

name=Перехідна характеристика по зовнішньому впливу;

stepview(Wf_c, name)

case 5

close;

end

end

case 2

while(m2~=4)

m2 = menu(Логарифмические характеристики,...

ЛАЧХ розімкненої системи,...

ЛАЧХ замкнутої системи,...

ЛАЧХ системи по зовн. впливу,...

Вихід);

switch m2

case 1

name = ЛАЧХ розімкненої системи;

get_bode(Ws_u, name )

case 2

name = ЛАЧХ замкнутої системи;

get_bode(Ws_c, name)

case 3

name = ЛАЧХ системи по зовн. впливу;

get_bode(Wf_c, name)

case 4

close;

end

end

case 3

while(m3~=4)

m3 = menu(Другие характеристики,...

Вплив синусоїдою ,...

Вплив імпульсами,...

Вплив синусоїдою по збуренню,...

Вихід);

switch m3

case 1

name = Вплив синусоїдою;

get_sin_response(Ws_c,name )

case 2

name = Вплив імпульсами;

get_pulse_response(Ws_c,name )

case 3

name = Вплив синусоїдою по збуренню;

get_sin_resp_f(Wf_c,name )

case 4

close;

end

end

case 4

close;

break;

end

end

function f = get_bode(Ws_u,name )

figure(3)

clf(reset);

[Gm,Pm,Wg,Wp] = margin(Ws_u);

subplot(3, 2, [1 2 3 4])

bode(Ws_u),grid on;

title(name)

subplot(3,2,5 )

description(1) = {[f Параметри ЛАЧХ]};

description(2) = { };

description(3) = {[fgain margin: m Gm]};

description(4) = {[fphase margin: m Pm]};

description(5) = {[fcrossover frequencies: m Wg]};

description(6) = {[fcrossover frequencies: m Wp ]};

axis([0 1 -1 0]);

axis off;

text(0,-1,description,VerticalAlignment,bottom);

subplot(3,2,6 )

textrow(4) = {[ ]};

textrow(5) = {[num2str(Gm)]};

textrow(6) = {[num2str(Pm)]};

textrow(7) = {[num2str(Wg)]};

textrow(8) = {[num2str(Wp)]};

axis([0 1 -1 0]);

axis off;

text(0,-1,textrow,VerticalAlignment,bottom);

function f = stepview(sys, name)

S = stepinfo(sys);

figure(2)

clf(reset);

subplot(3, 2, [1 2 3 4])

step(sys), grid on;

title(name);

[t_ss, os, n, t_r, y_inf] = get_step_info( sys ) ;

subplot(3,2,5 )

description(1) = {[fПараметри перехідного процесу для]};

description(2) = { };

description(3) = {[fУсталене значення m]};

description(4) = {[fЧас перехідного процесу m, с]};

description(5) = {[fЧас відпрацювання m]};

description(6) = {[fПеререгулювання m, %]};

description(7) = {[fКількість перебігів m, шт]};

axis([0 1 -1 0]);

axis off;

text(0,-1,description,VerticalAlignment,bottom);

subplot(3,2,6 )

textrow(4) = {[ ]};

textrow(5) = {[num2str(y_inf)]};

textrow(6) = {[num2str(t_ss)]};

textrow(7) = {[num2str(t_r)]};

textrow(8) = {[num2str(os)]};

textrow(9) = {[num2str(n)]};

axis([0 1 -1 0]);

axis off;

text(0,-1,textrow,VerticalAlignment,bottom);

function [t_ss, os, n, t_r, y_inf] = get_step_info( sys )

%get_step_info Compute step info

% self made step info function

t_ss=0

[y,t] = step(sys);

[val , ind] = max(y);

y_max = val;

% peak_time = t(ind)

l_y = length(y);

y_inf = y(l_y)

if (y_inf >0.005 || y_inf <-0.005 )

os = ((y_max-y_inf)/y_inf)*100;

for j = 1:l_y

ii = l_y - j;

if ii ==0

t_ss = 0

ii = 1

end

if (y(ii) > 1.02*y_inf) || (y(ii) < 0.98*y_inf)

t_ss = t(ii)

break;

end

end

for j = 1:l_y

if y(j) >= y_inf

t_r = t(j);

break;

end

end

n =0 ;

for j = 2:ii

if (y(j) < y_inf) && (y(j+1) > y_inf) || (y(j) > y_inf) && (y(j+1) < y_inf)

n= n+1;

end

end

if n==0

t_r = t_ss;

end

elseif (y_inf <=0.005) || (y_inf >= -0.005)

disp([kogda ustanovivshyesia znachenie blizko 0])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

os = y_max;

y_inf = 0;

for j = 1:l_y

ii = l_y - j;

if (y(ii) > -0.02*y_inf) || (y(ii) < 0.02*y_inf)

t_ss = t(ii);

break;

end

end

t_r = 0;

n =0 ;

for j = 2:ii

if (y(j) < y_inf) && (y(j+1) > y_inf) || (y(j) > y_inf) && (y(j+1) < y_inf)

n= n+1;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if n==0

t_r = t_ss;

end

end

function f = get_sin_response(sys,name )

figure(3)

clf(reset);

[z, t] = gensig(sin,0.65, 3, 0.01)

lsim(sys,z,t), grid on;

title(name);

function f = get_pulse_response(sys,name )

figure(3)

clf(reset);

[z, t] = gensig(pulse,0.35, 3 ), grid on;

lsim(sys,z,t)

title(name);

function f = get_sin_resp_f(sys,name )

figure(3)

clf(reset);

t=0:0.01:2;

z = 0.1*sin(16*t);

% [z, t] = gensig(pulse,0.35, 3 ), grid on;

lsim(sys,z,t),grid on;

title(name);