Основи статистичної динаміки систем управління
Характеристика динамічних систем. Методи первинного аналізу кореляційних сигналів. Обчислення логарифмічних частотних характеристик. Визначення спектральної щільності випадкового процесу. Cтатистичні характеристики сигналів на виході лінійної системи.
| Рубрика | Производство и технологии |
| Вид | методичка |
| Язык | украинский |
| Дата добавления | 10.03.2016 |
| Размер файла | 2,3 M |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
end;
20. Додайте в скрипт команди для побудови графіків обох спектрів на одному полі та перевірки значень дисперсій:
plot(w, S, w, Sk);
fprintf('Dr = %g\n', Dr)
fprintf('Dr(S): %g\n', trapz(w,S)/pi )
fprintf('Dr(Sk): %g\n', trapz(w,Sk)/pi )
Завантажте скрипт та поясніть отримані результати.
21. Повторіть розрахунки для кута зустрічі з хвилею 135 градусів та поясніть отримані результати.
Використовуючи наведену вище інформацію виконати наступне (варіанти для розрахунків отримати у викладача):
- побудувати і аналізувати спектри морського хвилювання;
- вивчити методи апроксимації спектрів хвилювання раціональними функціями;
- змоделювати випадкові процеси за допомогою формуючих фільтрів, що збуджуються білим шумом.
5. Оптимізація лінійної системи при випадкових збуреннях. Задача оптимальної фільтрації
Мета роботи: ознайомитись з методами та особливостями синтезу оптимальних лінійних фільтрів Вінера. Навчитись моделювати фільтри в SIMULINK-MATLAB та оцінити результати фільтрації.
Оптимальна фільтрація.
Розглянемо просте завдання оптимізації лінійної системи при випадкових збуреннях - задачу оптимальної фільтрації. При вимірах корисний сигнал спотворюється перешкодою , причому і корисний сигнал, і перешкода - стаціонарні центровані випадкові сигнали. Передбачається, що перешкода і корисний сигнал статистично незалежні, тобто, ніяк не зв'язані, породжуються різними джерелами. Спектральна щільність сигналу і перешкоди відома.
Завдання полягає в тому, щоб знайти передавальну функцію фільтру (безперервної системи), який виділяє сигнал з його суміші з перешкодою щонайкраще (будує його найкращу оцінку ), рис. 5.1.
При заданих спектральних щільностях і можна побудувати відповідні їм формуючі фільтри і (тут і далі ми перейдемо до змінної ):
, , (5.1)
Тоді схему можна зобразити в наступному вигляді, рис. 5.2.
Тут і - незалежні білі шуми одиничної інтенсивності. Для ясності додамо на блок-схему формування помилки :
Таким чином, потрібно знайти передавальну функцію фільтру , який забезпечує мінімальну дисперсію помилки . Пригадаємо, що знання спектральної щільності помилки дозволяє знайти її середній квадрат, рівний дисперсії для центрованих процесів:
, (5.2)
В останній рівності від кутової частоти ми перейшли до комплексної змінної. Якщо перешкоди немає, спектральна щільність помилки обчислюється по передавальній функції від входу до виходу :
,
де
, (5.3)
Аналогічно, якщо є перешкода, але немає сигналу, спектральна щільність помилки обчислюється по передавальній функції від входу до виходу :
, де , (5.4)
Якщо ж діють обидва сигнали, але вони статистично незалежні, спектральна щільність помилки дорівнює сумі приведених вище "окремих" спектрів:
, (5.5)
Тут і далі для скорочення запису у функцій змінної не враховано аргумент, а зірочка (верхній індекс) позначає заміну на . Можна згрупувати доданки трохи по-іншому:
, (5.6)
Виділивши в цьому виразі невідому передавальну функцію і "все інше", запишемо його в загальному вигляді:
, (5.7)
де, і - деякі функції, причому і (в даному випадку
,
і ). Потрібно знайти стійку передавальну функцію , яка забезпечує мінімум функціонала:
, (5.8)
Цю задачу вперше вирішив американський математик Норберт Вінер в 1940-х роках. Ми коротко розглянемо сучасний варіант рішення задачі Вінера, яке запропонували Х. Боде і К. Шенон.
Функція є деяким спектром, для якого можна побудувати формуючий фільтр , такий що всі його нулі (корені чисельника) і полюси (корені знаменника) мають негативні дійсні частини:
, (5.9)
Така операція називається спектральною факторизацією (розкладанням на співмножники). Тоді вираження для спектру помилки можна представити у вигляді:
, (5.10)
де
.
Ми тільки що виконали виділення повного квадрата в . Ця назва пов'язана з тим, що вираз при підстановці
є квадратом модуля частотної характеристики . В даному випадку .
Щоб отримати стійкий фільтр, виконаємо сепарацію (виділимо стійку і нестійку частину виразу):
, (5.11)
Якщо вираз - не строго правильна функція (степінь чисельника більше або дорівнює степені знаменника), існує безліч різних сепарацій. Серед них ми виберемо правильну сепарацію, при якій , - строго правильна функція (степінь її чисельника менше степені знаменника). Тоді вираз можна представити у вигляді:
де,
,
Використовуючи теорію функцій комплексних змінних, можна довести, що:
Таким чином, алгоритм побудови оптимального фільтру включає дві операції: спектральну факторизацію і сепарацію.
Для обчислення дисперсії помилки при використанні оптимального фільтру потрібно для отриманої функції знайти спектральну щільність і обчислити інтеграл. У середовищі MATLAB існує і інший спосіб: факторизувати отриману спектральну щільність помилки:
а потім обчислити (за допомогою функції norm) норму стійкої передавальної функції . Ця норма є СКВВ виходу системи при дії на її вхід одиничного білого шуму, тобто .
Основна частина команд вводиться в командному вікні середовища Matlab. Команди, які необхідно використовувати в інших вікнах, позначені іконками відповідних програм.
|
Етап виконання завдання |
Команди та ілюстрації |
|
|
1. Створіть новий m-файл (скрипт) з командами, які очищують робочу область та задають формуючі фільтри та спектральні щільності корисного сигналу та перешкоди (білий шум): clear all; Fx = tf(1, [1 1]); Fn = tf(0.1); Sx = Fx*Fx'; Sn = Fn*Fn'; |
||
|
2. Скопіюйте в робочу папку файли wiener.m та h2opt.m, додайте в скрипт команду [С,errOpt] = wiener( Sx, Sn ) та завантажте скрипт клавішею F5. Запишіть у звіт передаточну функцію оптимального фільтру та СКВВ помилки фільтрації (значення errOpt). |
||
|
3. Додайте в скрипт команди для підготовки даних до моделювання [nX,dX] = tfdata(Fx, 'v'); [nN,dN] = tfdata(Fn, 'v'); [nС,dС] = tfdata(С, 'v'); b = bandwidth(Fx*С); % полоса пропускання Fx(s)*W(s) tau = 2*pi/100/b; % інтервал кореляції для джерела шуму |
||
|
4. Побудуйте Simulink-модель задачі фільтрації: Для того щоб вивести на нижній осцилограф два сигнали, використовуйте мультиплексор (блок Mux з групи Signal Routing). Збережіть модель з ім'ям lab5sim.mdl. |
||
|
5. В параметрах джерел білого шуму встановіть одиничну інтенсивність (Noise Power), інтервал кореляції tau (Sample Time) та разное початкове значення послідовності випадкових чисел (Seed). |
||
|
6. В параметрах блока Noise model (блок Transfer Fcn з групи Continuous) задайте чисельник nN та знаменник dN (з робочої області). Для моделі корисного сигналу (Signal model) задайте чисельник nX та знаменник dX, а для фильтру (Filter) - nW та dW. |
||
|
7. В параметрах першого осцилографа (позначеного X+N) зніміть флажок, що обмежує кількість запам'ятовуваних даних (Limit data points), та задайте вивід результатів в масив xn. Аналогічно встановіть вивід другогоосцилографа (Error) в масив err, а третього (X, Xe) - в масив xxe. |
||
|
8. Встановіть час моделювання 50 с (меню Simulation - Simulation parameters - Stop Time). Додайте в скрипт команди для запуска моделювання та розрахунку СКВВ помилки за результатами моделювання: sim('lab5sim') errSim = std(err(:,2)) Порівняйте цей результат з отриманим раніше теоретичним значенням. |
||
|
9. Додайте команди для побудови графіка суміші "сигнал+шум", а також вихідного та відновленого сигналів (на окремому полі). figure(1) plot ( xn(:,1), xn(:,2) ); title('Signal + noise'); figure(2) title('Original and reconstructed signals'); plot(xxe(:,1),xxe(:,2),xxe(:,1),xxe(:,3)); legend('Original', 'Reconstructed'); Завантажте скрипт та скопіюйте графіки у звіт. |
print -dmeta |
|
|
10. Повторіть розрахунки для та . Додайте у ззвіт графіки теоретичного СКВВ помилки та СКВВ, отриманого при моделюванні. Зробіть висновки про вплив інтенсивності шуму вимірювань. |
||
|
11. В циклі побудуйте серію оптимальних регуляторів при зміні дисперсії перешкоди від 0,01 до 1 та побудуйте графік. Для цього в скрипт необхідно додати такий код: varN = [0.01:0.02:0.1 0.2:0.1:1]; % масив дисперсій перешкоди sigmaE = []; % пустий масив СКВВ помилки for i=1:length(varN) % цикл по всім варіантам Sn = tf(varN(i)); % спектр шума [C,errOpt] = wiener ( Sx, Sn ); % синтез регулятора C % показати регулятор sigmaE(i) = errOpt; % запам'ятати СКВВ помилки end; figure(3); % графік в новому вікні plot (varN, sigmaE ); Зробіть висновки щодо графіка. |
||
|
12. Додайте в масив varN великі (від 2 до 100) значення дисперсії перешкоди varN = [0.01:0.02:0.1 0.2:0.1:1 2:10 20:10:100]; і повторіть розрахунки. Поглянете, до якого значення прагне СКВВ помилки в оптимальній системи при збільшенні дисперсії перешкоди. Який при цьому виходить регулятор? Чи вирішує система свою задачу? Спробуйте обґрунтувати ці результати теоретично. |
6. Обчислення показників якості систем управління
Мета роботи: ознайомитись з програмною реалізацією типових процедур проектування для обчислення показників якості систем управління, використовуючи програмне середовище MATLAB.
Застосування спектральних алгоритмів ідентифікації та синтезу потребує значного обсягу математичних операцій з поліномами та матрицями. Для виконання таких розрахунків на ЕОМ можна скористатись програмним забезпеченням, розробленим спеціально для розв'язання подібних задач, або застосувати універсальні програмні системи для математичних обчислень.
1. Обчислення показників якості систем
Показник якості являє собою кількісну оцінку якості системи. В системах управління та вимірювальних системах найважливішою стороною якості є точність, яка головним чином характеризує успішність виконання системою покладених на неї задач.
Обчислення інтегральної квадратичної помилки
Нагадаємо, що якість багатовимірної динамічної системи при детермінованих сигналах оцінюють інтегральною квадратичною помилкою. Обчислення інтегральної квадратичної помилки може бути зведено до обчислювання інтегралів від дробово-раціональних функцій комплексної змінної вигляду .
(6.1)
де
, ,
причому a(s) має нулі лише в лівій півплощині змінної s =jщ.
Обчислити інтеграли такого вигляду (дисперсійні інтеграли) можна за допомогою спеціальних таблиць або програм для ЕОМ .
Програма-функція, що подана далі, написана мовою системи MATLAB і дозволяє обчислювати значення інтеграла (6.1) на основі спеціальних таблиць.
ss = coloss(b, а).
Вхідні змінні: а - поліном знаменника дробово-раціональної функції (підінтегрального виразу) ; b - поліном чисельника цієї функції. Вихідна змінна: ss - обчислене значення інтеграла.
Текст функції coloss(b,a) такий:
function ss = coloss (b,a)
ss = 0.0;
na = length(a); % Визначення степеня полінома знаменника
nb=length(b); b0=zeros(1,(na-nb-1)); b=[b0,b];
for k=1:1:na
if (a(k) <= 0.0) % Перевірка умови стійкості системи
ss = NaN; disp('Система нестійка'); return; end;
alfa = a(k)/a(k+1); beta = b(k)/a(k+1);
ss = ss + beta*beta/alfa; k2 = k + 2; if (k2 > na) break; end;
for i=k2:2:(na-1)
a(i) = a(i) - alfa*a(i+1); b(i) = b(i) - beta*a(i+1); end; end
ss = ss/2.0;% Значення інтегралу
Приклад 1. Знайти інтегральну квадратичну помилку одновимірної системи, що працює при детермінованих сигналах, якщо зображення сигналу помилки цієї системи має вигляд:
Значення інтегральної квадратичної помилки системи обчислюють за допомогою таких операторів:
b = [3 1 12 3 9 1]; а = [1 3 5 12 6 9 1];
disp(coloss(b,a));
Обчислення середньої квадратичної помилки
Як відомо, показником якості системи при випадкових впливах є середня квадратична помилка. З порівняння формули видно, що інтегральні показники якості систем при детермінованих і випадкових сигналах за своєю формою відрізняються лише коефіцієнтом перед символом інтеграла. В такому разі середню квадратичну помилку системи можна обчислити за допомогою функції coloss(b,a), помножуючи результат на 2р.
Приклад 2. Обчислити середню квадратичну помилку одновимірної динамічної системи, що працює при випадкових сигналах, якщо спектральна щільність сигналу помилки має вигляд:
Для використання функції coloss(b,a) останній вираз треба подати у вигляді добутку комплексно спряжених множників:
маємо:
b(s) = s + 2; a(s) = (s + l)(s + 3) = s2 + 4s + 3 .
В такому разі середню квадратичну помилку системи за допомогою функції coloss(b,a) можна обчислити так:
b = [1 2];а=[1 4 3]; s = coloss(b,a); disp(s*2*pi);
1.8326 % Значення інтегралу.
2. Факторизація спектральної щільності
Факторизацією спектральної щільності (за Вінером) називають подання спектральної щільності як парної функції частоти у вигляді комплексно спряжених співмножників:
Sxx = S(s)S(-s), (6.2)
Нулі і полюси співмножника S(-s) лежать в правій півплощині (ППП) комплексної змінної, а нулі і полюси S(s) - у лівій півплощині (ЛПП). В більшості розрахунків спектральну щільність стаціонарних випадкових процесів апроксимують дробово-раціональною функцією:
Операція факторизації цієї дробово-раціональної функції полягає у визначенні поліномів B(s), B(-s), A(s), A(-s) і включає знаходження її нулів і полюсів, сортування знайдених коренів, симетрично розміщених в лівій півплощині і правій півплощині комплексної змінної та відтворення відповідних поліномів чисельника та знаменника.
Такі операції в системі MATLAB можна виконати за допомогою функції FCWN, звернення до якої має вигляд:
function [Gpl,Got,tol]=FCWN(w)
%Вхід: SISO об'єкт w типа tf; вихід: два SISO об'єкта Gpl и Got типу zpk,
%пов'язаних з w співвідношенням:
%w(s)=Gpl(s)*Got(-s);
%Gpl - стійка частина w; Got - нестійка частина
%Пошук нулів z та полюсів p, коефіцієнта підсилення k для вихідного SISO об'єкта
[z,p,k]=zpkdata(w,'v');
%Перевірка додатної визначеності вихідного полінома чисельника
zz=abs(real(z));
if any(zz)<1e-06
tol=1;
else
tol=0;
end
% Перевірка додатної визначеності вихідного полінома знаменника
zz=abs(real(p));
if any(zz)<1e-06
tol=2;
end
%сортування на стійкі sp та нестійкі np полюси
j=0;np1=0;
for i=1:length(p)
if real(p(i))<0
j=j+1;sp(j)=p(i);
else
np1=np1+1;np(np1)=p(i);
end
end
% сортування на стійкі sz та нестійкі nz нулі
j=0;nz1=0;
for i=1:length(z)
if real(z(i))<0
j=j+1;sz(j)=z(i);
else
nz1=nz1+1;nz(nz1)=z(i);
end
end
k=sqrt(abs(k));
%Формування результату факторизації
%Визначення знака коефіцієнта k нестійкої частини
if rem(nz1,2)==0
if rem(np1,2)==0
k1=k;
else
k1=-k;
end
else
if rem(np1,2)==0
k1=-k;
else
k1=k;
end
Gpl=zpk(sz,sp,k);Got=zpk(nz,np,k1);
Приклад 3. Виконати факторизацію спектральної щільності, яка подана дробово-раціональною функцією такого вигляду:
Після факторизації отримуємо дробово-раціональну функцію, у якої всі нулі та полюси розташовані в ЛПП:
3.Сепарація функції комплексної змінної
Сепарацією називають представлення функції комплексної змінної у вигляді цілої частини і двох складових: в першій складові й особливості лежать тільки в ЛПП, а у другої - тільки в ППП:
N(s)=N0(s)+N+(s)+N- (s),
де N(s) - функція комплексної змінної s=jщ, яка підлягає сепарації; N0(s) - ціла частина, яка є поліномом аргументу s; N+(s) - правильний дріб аргументу s з особливостями тільки в ЛПП; N- (s) - правильний дріб аргументу s з особливостями тільки в ППП.
Вказані математичні операції в системі MATLAB можна виконати за допомогою функції [r,p,k] = residue(b,a).
Ця функція обчислює залишки, полюси та поліном цілої частини ДРФ, яка подана як відношення двох поліномів b(s) та a(s) за спадаючим степенем s:
де r1...rn - вектор-стовпець лишків; p1,...pn - вектор-стовпець полюсів; k- вектор-рядок коефіцієнтів цілої частини дробово-раціональної функції.
Для полюсів (коренів знаменника) кратності т розкладення на прості дроби включає додатково члени:
Функція [b,a] = residue(r,p,k) з двома вхідними та трьома вихідними параметрами обчислює дробово-раціональну функцію у вигляді відношення двох поліномів b та а.
На базі функції residue написана функція sep. Звернення до неї має вигляд:
[bs,as,ks] = sep(b,a),
Вхідні змінні:
b - коефіцієнти полінома чисельника; а - коефіцієнти полінома знаменника.
Вихідні змінні: bs, as, ks - коефіцієнти полінома чисельника та знаменника сепарованої дробово-раціональної функції, що має особливості в ЛПП, та поліном її цілої частини. Текст цієї функції має вигляд:
function [bs,as,ks] = sep(b,a)
[r,p,k] = residue(b,a); m=length(p); j=0;
for i=1:m
if(real(p(i))<=0) j=j+1; rs(j)=r(i); ps(j)=p(i); end; end
[bs,as,ks]=residue(rs,ps,k);
Приклад 4. Виконати сепарацію дробово-раціональної функції:
яка мас три кратних полюси при s = -1. Для виконання сепарації треба ввести поліноми чисельника та знаменника і скористуватись функцією sер(b,а):
b=[2 6 5 -5]; а=[1 2 0 -2 -1];
[bs,as,ks]=sep(b,а);
У результаті виконання операторів отримуємо поліноми чисельника bs та знаменника as передаточної функції, що має полюси тільки в ЛПП:
4. Обчислення логарифмічних частотних характеристик
Для обчислення логарифмічних частотних характеристик доцільно використати функцію freqs, яка виконує розрахунок частотної характеристики динамічної системи за відомою передаточною функцією, яка задана векторами коефіцієнтів чисельника b та знаменника а.
Функція [H,w] = freqs(b,a,n) автоматично вибирає n значень відліків на частотах w, для яких обчислюється комплексна частотна характеристика. Діапазон частот визначається автоматично. Якщо n не вказано, то за умовчанням n = 200.
Щоб сформувати дані для побудови графіків амплітудно-частотної та фазо-частотної характеристик необхідно обчислити модуль та фазу вектора Н за виразами:
mag = abs(H); phase = angle(H).
Для перетворення одиниць вимірювання в герци, градуси, децибели використовують співвідношення:
f=w/(2*pi); phase = phase* 180/рі; mag = 20*logl0(mag).
Функція freqs(b,a,w) без вихідних аргументів обчислює частотну характеристику системи та виводить в графічне вікно графіки амплітудно-частотної та фазочастотної характеристик.
Приклад 5. Обчислити комплексну частотну характеристику системи з передавальною функцією
Побудувати графіки амплітудно-частотної та фазочастотної характеристик. Для вирішення задачі можна застосувати таку програму:
b=[5]; a=[0.01 0.05 1]; [H,w]=freqs(b,a);
f=w/(2*pi); mag=20*log10(abs(H)); phase=angle(H)*180/pi;
subplot(2,1,1), semilogx(f,mag,'k'), grid;
title('Амплітудно-частотна характеристика','Fontsize', 12);
ylabel(' Амплітуда, Дб','Fontsize', 12);
subplot(2,1,2), semilogx(f,phase,'k'),grid;
title('Фазо-частотна характеристика','Fontsize', 12);
xlabel('Частота, Гц','Fontsize',12);
ylabel('Фаза, rpaд','Fontsize', 12);
Обчислювати логарифмічні частотні характеристики для lti-моделей динамічних систем в tf-, zpk- та ss-формах зручно за допомогою функцій і команд bode в декількох варіантах використання.
Команда bode(sys) будує на екрані графіки логарифмічних частотних характеристик для lti-моделі sys. Модель може бути безперервною або дискретною, одновимірною або багатовимірною. Для багатовимірної моделі ця команда будує множину логарифмічних частотних характеристик для кожного каналу системи від входу до виходу. Кількість точок графіка та діапазон частот визначаються автоматично.
Команда bode(sys, w) будує графіки логарифмічних частотних характеристик у заданому діапазоні частот (в рад/с), що заданий вектором w = {wmin, wmax}. Для побудови частотних характеристик на заданих частотах створюється вектор значень частот w. Щоб сформувати масив рівновіддалених у логарифмічному масштабі частот w слід застосувати команду logspase.
Команда bode(sysl, sys2,....sysN, w) дозволяє будувати частотні характеристики для декількох моделей на одному графіку. Вектор частот w = {wmin, wmax} не є обов'язковим. Всі моделі повинні мати однакову кількість входів і виходів і можуть бути як безперервними, так і дискретними. Можна також задати колір, тип лінії і маркера для кожної системи, наприклад, так bode(sysl,'r',sys2,'y - ',sys3,'gx').
Для створення більш інформативних графіків доцільно використовувати функції: [mag,phase,w]= bode(sys,...) та [mag,phase]= bode(sys,w). Функції повертають вектор частоти w та масиви амплітудних mag та фазових phase частотних характеристик. Частота обчислюється в радіанах в секунду, фаза - в градусах. Графік при цьому не виводиться. Вихідні аргументи mag, phase є тривимірними масивами, остання розмірність яких - частота. Для подання амплітудної характеристики в децибелах відповідний масив можна перетворити таким чином: magdb = 20*logl0(mag).
Приклад 6. Визначити амплітудно-частотні та фазочастотні характеристики систем з передавальними функціями:
Амплітудно-частотну характеристику подати у децибелах. Вивести на графіку заголовок, позначення осей та пояснення. Один з варіантів вирішення задачі такий:
b=[4];
a=[0.01 0.02 1];
G=tf(b,a);
bl=[0.5];al=[0.04 0.05 1];
G1=tf(bl,al);
[mag,phase,w]=bode(G);
[mag1,phase1]=bode(G1,w);
mag=20*log10(squeeze(mag));ph=squeeze(phase);
mag1=20*log10(squeeze(mag1));
ph1 = squeeze(phase1);
subplot(2,1,1),hn=semilogx(w,mag,'-k',w,mag1,'-ko'), grid on;
set(hn,'markersize',4); legend('G(s)','G1(s)');
title('Амлітудно-частотна характеристика','fontsize', 12);
ylabel('Амплітуда, Дб','fontsize', 12); subplot(2,1,2),hn=semilogx(w,ph,'-k',w,ph1,'-ko'), grid on;
set(hn,'markersize',4);
legend('G(s)','G1(s)')
title('Фазо-частотна характеристика', 'fontsize', 12);
ylabel('Фаза, rpaд','fontsize', 12);
ylabel('Частота, рад/c','fontsize', 12);
5. Обчислення перехідної функції
Для lti-моделей динамічних систем в tf-, zpk- та ss-формах перехідну функцію можна обчислити за допомогою групи функцій і команд step.
Команда step(sys) будує на екрані графік перехідної функції для lti-моделі sys. Модель може бути безперервною або дискретною, одновимірною або багатовимірною. Для багатовимірної моделі ця команда будує множину перехідних функцій для кожного каналу системи. Тривалість процесу визначається автоматично і встановлюється такою, щоб визначити його основні особливості.
Команда step(sys, t) будує графіки перехідної функції в заданому діапазоні часу. Діапазон може бути заданий моментом закінчення t = tf (в секундах) або у вигляді вектора t = 0:ts:tf. Діапазон часу слід вибирати таким, щоб відобразити перехідний процес з необхідною точністю.
Команди step(sysl, sys2, ... sysN), step(sysK, sys2, ... sysN, t) дозволяють будувати перехідні функції для декількох моделей на одному графіку. Всі моделі повинні мати однакову кількість входів і виходів і можуть бути як безперервними, так і дискретними.
Команди step(sysl,sl,...sysN,sN), step(sysl,sl,...sysN,sN,t) дозволяють будувати перехідні функції для декількох моделей на одному графіку і позначити їх різними стилями. Рядкові змінні sl,...,sN можуть мати до трьох символів, які позначають тип лінії, тип маркера та їх колір.
Функції [у, t, х]= step(sys) та [у, t, х]= step(sys,t) обчислюють перехідні функції для векторів виходів у, векторів часу t, значення змінних стану х. Графіки на екран не виводяться.
6.Обчислення імпульсної перехідної функції
Для lti-моделей динамічних систем в tf-, zpk- та ss-формах імпульсну перехідну функцію можна обчислити за допомогою групи функцій і команд impulse. Властивості команд і функцій impulse подібні властивостям команд і функцій step.
Команда impulse(sys) будує на екрані графік імпульсної перехідної функції для lti-моделі sys. Модель може бути безперервною або дискретною, одномірною або багатовимірною. Для багатовимірної моделі ця команда будує множину імпульсних перехідних функцій для кожного каналу системи. Тривалість процесу визначається автоматично і встановлюється такою, щоб визначити основні особливості перехідних процесів.
Команда impulse(sys, t) будує графіки перехідної функції в заданому діапазоні часу. Діапазон може бути заданий моментом закінчення t = tf (в секундах) або у вигляді вектора t = 0:ts:tf. Діапазон часу слід вибирати таким, щоб відобразити перехідний процес з необхідною точністю.
Команди impulse(sys 1 ,sys2,...sysN) та impulse(sysl,sys2,... sysN, t) дозволяють будувати імпульсні перехідні функції для декількох моделей на одному графіку. Всі моделі повинні мати однакову кількість входів і виходів.
Команди impulse(sysl,sl,...sysN,sN) та impulse(sysl, si,... sysN, sN, t) дозволяють будувати перехідні функції для декількох моделей на одному графіку і позначити їх різними стилями.
7. Функції перетворення передаточних функцій.
Функція tsys = ctranspose(sys) використовують для побудови спряженої моделі.
Функція isys=inv(sys) виконує інверсію моделі і дає можливість отримати обернену передавальну функцію або матрицю передавальних функцій системи.
Під час перетворення передавальних функцій систем за допомогою вказаних вище операцій результати, як правило, ускладнюються, тому що цими операціями не передбачається скорочень в передавальних функціях. Що б позбавитись таких ускладнень, у виразах доцільно використовувати функцію sysr=minreal(sys,tol), яка виконує скорочення однакових нулів та полюсів передавальної функції. При скороченні функція виконує пошук однакових у межах допуску tol нулів та полюсів. Якщо допуск не вказано, то за умовчанням його значення tol=sqrt(eps), де eps - машинна точність.
8. Математичні операції з поліномами
Поліноми від комплексної змінної s = jщ використовують для опису моделей динамічних систем у формі передавальних функцій, а також динамічних характеристик випадкових та детермінованих сигналів в них.
В системі MATLAB поліноми задаються у вигляді вектора-рядка, елементами якого є коефіцієнти полінома.
Коефіцієнти вводяться в квадратних дужках і відділяються пропусками або комами. Вводити треба всі коефіцієнти, включно з нульовими. Для доступу до окремих коефіцієнтів поліномів використовують круглі дужки. В системі MATLAB передбачено виконання декількох операцій над поліномами, деякі з яких доцільно застосувати для реалізації описаних спектральних алгоритмів.
Множення поліномів здійснює функція с = conv(a,b), яка повертає вектор-рядок с коефіцієнтів полінома, який є добутком заданих поліномів а та b.
Ділення поліномів виконує функція [q,r] = deconv(b,a). Вона повертає коефіцієнти полінома q, що є часткою від ділення поліномів а і b, а також залишок у векторі r так, що b = conv(a,q) + r.
Корені полінома обчислює функція r = roots(p). Вона для заданого полінома р повертає вектор-стовпчик r, елементами якого є корені полінома р.
Побудувати поліном за його коренями дозволяє функція р = роlу(r). Для заданого вектора-стовпця r коренів деякого полінома функція обчислює вектор-рядок р коефіцієнтів цього полінома.
Обчислити значення полінома за заданим значенням його аргументу можна за допомогою функції у = polyval(p, х). Аргументами функції є заданий вектор коефіцієнтів полінома р та задане значення аргументу х.
7. Оптимальний стохастичний спостерігач (фільтр Калмана)
Мета роботи: навчитись відновлювати вектор стану системи при стохастичних збуреннях за допомогою стохастичного спостерігача (фільтр Калмана), використовуючи програмне середовище MATLAB.
1-й етап. Припустимо, що рівняння реальної системи має вигляд:
х = Ах + Ви + щ1, (7.1)
y = Сх + щ2, (7.2)
Через щ1, позначається шум. який збурює стан, а через щ2 - шум спостережень чи вимірювань.
Таким чином, на першому етапі нам необхідно задати матриці [A,B,C,D] системи в просторі станів, а також матриці коваріацій шумів стану та вимірювань.
Розглянемо випадок синтезу оптимального стохастичного спостерігача для літального апарата Стан системи збурюється турбулентністю атмосфери, яка є кольоровим шумом. Стандартна задача синтезу оптимального стохастичного спостерігача припускає, що на систему діє білий шум. Для приведення поставленої задачі до стандартного виду необхідно використовувати формуючий фільтр (фільтр, на вхід якого поступає білий шум з. а на виході маємо кольоровий шум щ1, спектральні щільності якого описують турбулентність атмосфери) в структурі об'єкта, тобто розширити об'єкт. Тобто, для синтезу ми будемо використовувати об'єкт з розширеним простором стану, що описується матрицями [Aех,Bех,Cех,Dех].
2-й етап. Допустимо, то спостерігач повного порядку виду:
відповідає системі, що описується рівняннями (7.1), (7.2).
Тоді помилка відновлення визначається виразом:
, (7.3)
Задача оптимального вибору величини К називається задачею оптимального спостереження (задачею побудови оптимального спостерігача), тобто нам необхідно вибрати такі коефіцієнти спостерігача, щоб мінімізувати вираз (7.3).
Матрицю К можна отримати, вирішивши алгебраїчне рівняння Ріккаті для спостерігача:
де Q - рішення рівняння Ріккаті; V1, V2 - інтенсивності білого шуму.
Тоді матриця коефіцієнтів підсилення оптимального спостерігача буде дорівнювати:
Схема з'єднання об'єкта зі спостерігачем наведена на рис. 7.1.
Таким чином, на другому етапі ми здійснюємо синтез фільтра Калмана.
В результаті синтезу отримаємо спостерігач, що описується четвіркою матриць в просторі стану, на виходе якого маємо відновлений вектор стану системи (рис. 7.1).
В пакеті MATLAB можна вирішити задачу побудови оптимального стохастичного спостерігача, виконавши наступну послідовність дій:
1. Задати четвірку матриць об'єкта:
2. В реальних системах виміряти стан точно не можна через шуми вимірювань, крім того на систему завжди діє збурення. Задамо матрицю збурень G і створимо модель збуреного об'єкта в просторі стану.
3. Стан системи, на яку діють збурення можна відновити за допомогою оптимального стохастичного спостерігача (фільтра Калмана). Для цього необхідно задати матрицю коваріацій шумів вимірювань Qn та матрицю коваріацій шумів стану Rn. Задачу синтезу в MATLAB вирішуємо за допомогою оператора kalman.
4. Коли стан системи відновлено, можна застосувати закони синтезу оптимального детермінованого регулятора.
V =[1.5 1.5 2.298 0.8 0.5];|
Q=diag(V);
V1=[1 1.17];
R =diag(VI);
N=zeros(5,2);
N(3.2)=-0.624;
[K,S,E]=lqr(A-L*C,B,Q,R,N)
5. Для отримання характеристик системи записуємо матриці стану замкненої системи. За допомогою операторів impulse та step одержуємо імпульсну та перехідну характеристики замкненої системи
figure (1)
impulse(A-L*C-B*K,B,C,D)
figure(2)
step(A-L*C-B*K,B,C,D)
6. Якість синтезованої системи оцінюємо за допомогою Н2 норми:
Ac=A-L*C-B*K;
Необхідно за допомогою стохастичного оптимального спостерігача відновити стан системи та аналітично побудувати оптимальний детермінований регулятор.
8. Структурна ідентифікція динамічного об'єкта
Мета роботи: вивчити методи алгоритми структурної ідентифікації динаміки об'єкту. Розробити програмну реалізацію алгоритму структурної ідентифікації.
Структурна ідентифікація дозволяє визначити як порядки і кількість диференційних рівнянь, які описують динаміку об'єкта, що ідентифікується, так і коефіцієнти цих рівнянь, крім того новітні методи структурної ідентифікації в частотній області можуть бути використані для оцінювання динамічних властивостей збурень, дія яких проявляється в реальних умовах експлуатації.
Новітні алгоритми структурної ідентифікації динаміки об'єкту і діючого на нього узагальненого збурення володіють рядом переваг, серед яких слід відзначити точність; забезпечення побудови моделей динаміки в режимі нормальної експлуатації; зручність, простота та інформативність; наукова обґрунтованість; можливість використання достатньо простої апаратури нечутливої до зовнішніх впливів.
Застосування їх у практиці істотно скорочує строки та обсяги експериментів, що мають за мету одержання певних обсягів інформації.
Запропонований відомим українським вченим Блохіним Л.М. метод не тільки дає можливість отримати найкращі за мінімумом дисперсії похибки оцінки динамічних характеристик досліджуваного об'єкту за даними вимірювання вхідного та вихідного сигналів в реальних експлуатаційних умовах, а також визначити спектральну щільність неконтрольованого збурення.
У досліджуваному випадку об'єкт ідентифікації (рис. 1.1) описується звичайним диференційним рівнянням вигляду:
, (8.1)
де ш - зовнішні стохастичні збурення, які діють на об'єкт управління;
u - сигнал керування; x - вихідний сигнал;
Р і М - операторні поліноми.
При використанні обраного методу вважається, що:
- вхідний сигнал u і вихідний сигнал х представляють собою стаціонарні випадкові процеси з відомими спектральними та взаємними спектральними щільностями Suu(s), Sхх(s) та Sux(s), які були визначені за результатами вимірювань сигналів u і х;
- збурення ш не корельовано з вхідним сигналом u.
Задача структурної ідентифікації полягає в тому, щоб за статистичними даними “вхід - вихід” отримати найкращі за точністю оцінки динамічних характеристик об'єкту управління та збурення на класі дробово-раціональних функцій.
Згідно обраного алгоритму шукана оптимальна матриця оцінок визначається з виразу:
, (8.2)
де Ф11 - передаточна функція об'єкту; Ф12 - передаточна функція фільтру, який формує динамічні характеристики збурення, приведеного до виходу системи із заданою спектральною щільністю; R0 - ваговий коефцієнт; D - результат факторизації блочної матриці Syy:
; (8.3)
Т0+Т+ - стійка частина результату сепарації блочної дробово-раціональної матриці:
, (8.4)
де “*” - знак Ермітового спряження матриць.
Матриця взаємних спектральних щільностей Syx
(8.5)
де S?x - результат факторизації виразу Sx?S?x , який входить до рівняння зв'язку:
. (1.6)
Для аналізу якості отриманих моделей використовується функціонал, який представляє собою дисперсію похибки ідентифікації вигляду:
. (8.7)
де - спектральна щільність помилки ідентифікації.
Значення може бути розраховане з використанням теореми Вінера-Хінчина:
. (8.8)
Тоді функціонал якості (1.7) з урахуванням (1.8) набуває вигляду:
, (8.9)
де
; ; .
Для реалізації даного методу необхідно виконати наступні етапи:
1. спектральні та взаємні спектральні щільності векторів вхідного та вихідного сигналів Suu(s), Sхх(s) та Sux(s) обчислені за результатами ідентифікаційного експерименту, підставити до рівняння зв'язку і визначити ;
2. скласти матриці S'yy і S'yx;
3. отримати вирази для D і D* після виконання факторизації матриці ;
4. визначити матрицю
;
5. виконати сепарацію матриці Т і знайти вираз для ;
6. визначити блокову матрицю ;
7. після лівостороннього видалення полюсів Ф11 отримати шукані оцінки операторних поліномів М та Р;
8. визначити спектральну щільність узагальненого збурення S'шш.
Для одновимірної системи обчислювальні операції виконуються з дробово-раціональними функціями від змінної s=jщ.
На наступному етапі виконується аналіз якості ідентифікації відповідно до функціоналу 8.9 та за необхідності відбувається спрощення (редукування) отриманих моделей.
Спектральна щільність вхідного сигналу
Спектральна щільність вихідного сигналу
де в=уш/уu - співвідношення шум/сигнал
Взаємна спектральна щільність сигналів
Підставляючи вихідні дані в рівняння зв'язку маємо (етап 1):
Розкладаючи на комплексно спряжені множники отримуємо наступні спектральні щільності:
та .
Знайдемо матриці (етап 2):
і .
Знаходимо матриці D, D*, а також D-1, D*-1 (етап 3). Спочатку знаходимо:
.
Отриману матрицю розкладаємо на дві ермітово-спряжені:
, .
Відповідні обернені матриці:
, .
Використовуючи знайдені значення обчислюємо матрицю Т (етап 4)
.
Відповідно до алгоритму необхідно виконати сепарацію отриманої матриці. Проте отримана матриця не містить полюсів(коренів знаменника) у ППП, тому сепарація матриці не потрібна. Отже:
.
Визначаємо блочну матрицю Ф (етап 6):
З отриманої блочної матриці визначається передавальна функція динамічного об'єкта :
кореляційний частотний логарифмічний спектральний
передаточна функція фільтру, який формує динамічні характеристики збурення, приведеного до виходу системи:
.
Визначаємо спектральну щільність збурення:
.
Література
1. Блохин Л. Н. Оптимальные системы стабилизации. - К.:Техніка, 1982. - 144 c.
2. Блохін Л. М., Буриченко М. Ю. Статистична динаміка систем управління: Підручник. - К.: НАУ, 2003. - 208 с.
3. Венгеров А. А., Щаренский В. А. Прикладные вопросы оптимальной линейной фильтрации. М., Энергоиздат, 1982. 190с.
4. Кузовков Н. Т. Модальное управление и наблюдающие устройства. М., Машиностроение, 1976. 184с.
5. Солодовников В. В. Основы теории и элементов систем автоматического регулирования. М., Машиностроение, 1985. 510с.
6. Черных И. В. SIMULINK: среда создания инженерных приложений / под общ. ред. к.т.н. В. Г. Потемкина. - М.: ДИАЛОГ-МИФИ, 2003. - 496 с.
7. MATLAB в инженерных и научных расчетах / [Дащенко А. Ф., Кириллов В. Х., Коломиец Л. В., Оробей В. Ф.]. - О.: Астропринт, 2003. - 210с.
Размещено на Allbest.ru
Подобные документы
Властивості та функціональне призначення елементів системи автоматичного керування. Принцип дії, функціональна схема, рівняння динаміки. Синтез коректувального пристрою методом логарифмічних частотних характеристик. Граничний коефіцієнт підсилення.
курсовая работа [2,9 M], добавлен 22.09.2013Опис принципової схеми та принципу дії гідравлічного слідкуючого приводу. Складання рівнянь динаміки системи автоматичного керування та їх лінеаризація. Створення структурної схеми даної системи та аналіз її стійкості. Побудова частотних характеристик.
курсовая работа [252,1 K], добавлен 31.07.2013Визначення передаточних функцій, статичних та динамічних характеристик об’єкта регулювання. Структурна схема одноконтурної системи автоматичного регулювання. Особливості аналізу стійкості, кореляції. Годограф Михайлова. Оцінка чутливості системи.
курсовая работа [1,1 M], добавлен 10.01.2015Застосування неруйнівного контролю для визначення показників якості матеріалів без порушення їх властивостей та функціонування. Класифікація сигналів та методів дефектоскопії. Аналіз придатності виробів на підставі норм бракування та умов експлуатації.
курсовая работа [283,3 K], добавлен 11.09.2014Проблеми забезпечення необхідних властивостей лінійних автоматичних систем. Застосовування спеціальних пристроїв, для корегування динамічних властивостей системи таким чином, щоб забезпечувалася необхідна якість її функціонування. Методи їх підключення.
контрольная работа [605,5 K], добавлен 23.02.2011Мета впровадження автоматичних систем управління у виробництво. Елементи робочого процесу в парокотельній установці. Вибір структури моделі об'єкта регулювання та розрахунок її параметрів. Розрахунок параметрів настроювання автоматичних регуляторів.
курсовая работа [986,6 K], добавлен 06.10.2014Характеристика основних положень термодинаміки. Аналіз термодинамічних процесів ідеального газу. Поняття, структура та призначення теплового насосу. Принцип розрахунку теплообмінних апаратів. Методи термодинамічного аналізу енерго-технологічних систем.
учебное пособие [2,5 M], добавлен 28.11.2010Опис видів котлів-утилізаторів і характеристика автоматичної системи регуляції температури перегрітої пари на виході з котла-утилізатора КУ-80. Розрахунок метрологічних характеристик вимірювальних каналів АСР. Структурна схема функцій і надійності АСР.
дипломная работа [1,4 M], добавлен 31.03.2011Основи управління якістю та її забезпечення в лабораторіях. Виникнення систем управління якістю. Поняття якості результатів діяльності для лабораторії. Розробка системи управління якістю випробувальної лабораторії. Проведення сертифікаційних випробувань.
дипломная работа [4,0 M], добавлен 15.12.2011Автоматизація роботи підприємств по виготовленню бетонних ростворів, автоматичне управління технологічним процесом. Теоретичні основи технологічного процесу в окремих технологічних апаратах і машинах. Розроблення системи автоматичного керування.
курсовая работа [2,5 M], добавлен 26.09.2009
