ТПРК2-Курсова

Ця стаття відноситься до групи довідкових статей для студентів ФІН.

Ця стаття не закінчена. Будь ласка, допишіть її.

Курсова робота складається з декількох пунктів. Усі вони виконуються у середовищі Матлаб.
Важливо також розуміти те, що ви робите під час курсової. Для цього слід читати матеріали, що їх дає Руслан Костянтинович. Нижче поданий практичний підхід до розв'язку поставлених задач.

Пошук матричної експоненти

Дано: матриця A
Задача: знайти матрицю e^At

Перелік дій (копіпейст з консолі Матлаба, версія r2010а):

1. Ввести матрицю А - рядки розділяються крапками з комою, наприклад

>>A=[-9 7 1 5 6; 9 8 7 -3 7; -1 -6 0 2 -10; 5 -4 -4 9 -9; -9 4 6 1 -8]

запише в пам'ять матрицю: $A=\begin{pmatrix} -9 & 7 & 1 & 5 & 6 \\ 9 & 8 & 7 & -3 & 7 \\ -1 & -6 & 0 & 2 & -10 \\ 5 & -4 & -4 & 9 & -9 \\ -9 & 4 & 6 & 1 & -8 \end{pmatrix}$

2. Знайти власні числа та власні вектори А:

>>[V, D] = eig(A);

Перша матриця (V) буде складатися з власних векторів матриці А, а друга (D) буде мати на діагоналях її власні числа.

3. Обернути матрицю з власних векторів:

>>V_inv = inv(V);

Нагадаю також, що наявність крапки з комою після команди означає, що на екран не буде нічого виводитись (якщо, звісно, ви не припустилися синтаксичної помилки).

4. Оскільки матрична експонента - це, фактично, елемент розкладу в ряд Тейлора по часу, то треба ввести відповідну змінну "час" (t):

>> t = sym('t');

5. Наступний крок - знайти матрицю e^λt, що буде містити експоненту від власних чисел матриці А на своїх діагоналях (і нулі на інших позиціях)

>> e_lambda = [exp(D(1,1)*t) 0 0 0 0; 0 exp(D(2,2)*t) 0 0 0; 0 0 exp(D(3,3)*t) 0 0; 0 0 0 exp(D(4,4)*t) 0; 0 0 0 0 exp(D(5,5)*t)];

6. Тепер є всі необхідні матриці для підрахунку матричної експоненти за формулою V ^{− 1}e^λtV

>> A_exp = V*e_lambda*V_inv;

Не лякайтеся виводу!

Вирішення однорідної системи

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

>> initial_state = [1;0;0;0;0];
>> solution = A_exp*initial_state;

Перевірка на стійкість

Для цього достатньо розглянути власні числа матриці А:

>>[V, D] = eig(A);

З цих двох матриць цікавість в даному питанні несе матриця D. Якщо на її діагоналі присутній хоча б один додатний елемент, то вона система не є асимптотично-стійкою. Очевидно, якщо всі числа на діагоналі від'ємні, то система буде асимптотично-стійкою.

Полюси

Полюсами системи є ті самі власні числа матриці А. В якості полюсів стійкості перерахуйте від'ємні діагональні елементи матриці D (якщо такі є). В якості полюсів нестійкості перерахуйте додатні діагональні елементи матриці D (знову ж таки, якщо такі є).

Перевірка на керованість

1. Для цього завдання на додачу до матриці А йдуть матриці B та С. Система диференціальних рівнянь, що її повторюють на лекціях кожного разу, заслуговує на увагу:

$\dot{x(t)}=Ax(t)+Bu(t) \dot{y(t)}=Cx(t)$ Як бачимо, тут не вистачає деякої матриці D (за логікою); Матлаб вважає її наявність потрібною, щоб побудувати відповідну модель (простору станів). Після вводу матриць (див. початок статті) можна побудувати відповідну модель простору станів:

>> System = ss(A, B, C, 0);

2. Матриця керованості даної моделі - специфічне утворення, що про нього можна взнати на лекціях. Практично отримати таку матрицю тривіально:

>> Controllability_matrix = ctrb(System);

3. Для того, щоб отримати інформацію про керованість системи, потрібно порахувати ранг її матриці керованості:

>> m = rank(Controllability_matrix);

Якщо це значення дорівнює рангу матриці А, то система є повністю керованою. Відповідні полюси керованості визначаються власними числами матриці А аналогічно до полюсів стійкості (тобто "нестійкі" власні числа відповідають полюсам некерованості). Також потрібно зауважити, що канонічна форма керованості збігається з початковим представленням.

Якщо система не є повністю керованою

На щастя, багатьох студентів оминала доля отримувати подібні системи. Тому раджу відшукати помилку у введених даних. Якщо така відсутня, то проконсультуйтеся, будь-ласка, з викладачем, та поділіться враженнями у цій статті :)

Перевірка на стабілізованість

Якщо ваша система є повністю керованою, то вона автоматично є стабілізованою.

В іншому випадку

Див. попередній пункт, якщо система не є повністю керованою

Перевірка на відновлюваність

Для виконання цього кроку потрібні ті самі маніпуляції, що й під час обрахунків для керованості. Єдина різниця - те, що властивість відновлюваності задається матрицею відновлюваності. Особливості побудови розповідаються на лекціях; з практичної точки зору потрібні наступні команди (з урахуванням того, що ми вже ввели матриці та побудували з них систему):

>> Observability_matrix = obsv(System);
>> m = rank(Observability_matrix);

Наступні міркування про повну відновлюваність системи, її полюси відновлюваності та канонічну форму відновлюваності аналогічні відповідним міркуванням про повну керованість, полюси керованості та канонічну форму керованості.

Якщо система не є повністю відновлюваною

Див. той пункт, де система не є повністю керованою.

Імпульсна характеристика системи

Отримати її легко:

>> impulse(system)

Приклад виводу:

Головне - інтерпретація. Вона залежить від стійкості системи. За наявності полюсів нестійкості обов'язково знайдуться графіки, що будуть нескінченно зростати (тобто, прямувати до плюс- або мінус-нескінченності). Тобто, якісь із вхідних даних не можна буде відслідкувати за будь-який період часу, що пройшов після їх введення.

Пошук матричної передавальної функції

Для виконання цього кроку необхідно спершу знайти резольвенту матриці А.

1. Вводимо додаткову змінну s та шукаємо резольвенту за формулою (sI − A) ^{− 1}. Матриця А, як правило, є квадратною й відомою (і кількість її рядків дорівнює, скажімо, dim); це знадобиться для побудови відповідної діагональної матриці (функція eye в Матлабі).

>> s=sym('s');
>> resolvent = inv(s*eye(dim)-A);

2. Матрична передавальна функція знаходиться за формулою H(s) = C(sI − A) ^{− 1}B.

>> H_s = C*resolvent*B;

Знову не лякайтеся виводу!

Замикання системи диференційних рівнянь зі сталими коефіцієнтами за допомогою зворотного зв'язку

1. Першим кроком в даному випадку є власноручне (!) задання полюсів системи. Тобто в данному випадку "конструювання" системи йде в зворотньому напрямку. Для завдання буває потрібним задати полюси три рази, кожен раз вдесятеро менші:

>> poles_1 = [-1, -2, -3, -4, -5];
>> poles_2 = [-10, -20, -30, -40, -50];
>> poles_3 = [-100, -200, -300, -400, -500];

2. Система була побудована раніше (а дані до неї введені ще раніше). Матриця зворотнього зв'язку будується, як завжди, дуже просто:

>> F_1 = place(A, B, poles_1);
>> F_2 = place(A, B, poles_2);
>> F_3 = place(A, B, poles_3);

В даному прикладі є три матриці - по одній на кожен набір полюсів.

3. За допомогою нової матриці можна сконструювати нову, замкнену систему, повне формулювання якої зустрічається в лекціях. В Матлабі ці системи будуються просто:

>> A_prime_1 = A-(B*F_1);
>> system_1 = ss(A_prime_1, B, C, 0);
>> A_prime_2 = A-(B*F_2);
>> system_2 = ss(A_prime_2, B, C, 0);
>> A_prime_3 = A-(B*F_3);
>> system_3 = ss(A_prime_3, B, C, 0);

4. Наступним кроком є порівняння імпульсних характеристик отриманих систем:

>> impulse(system_1)
>> impulse(system_2)
>> impulse(system_3)

Замкнені системи будуть мати асимптотичну стійкість (тобто на всіх графіках траєкторія прямує до нуля). Швидкість "занулення" траєкторій зростає разом із зменшенням полюсів (саме зменшення, а не зменшення модуля!).

Пошук оптимального лінійного детермінованого регулятора

1. Спершу потрібно (як і у випадку з полюсами) довільно задати параметри. Тут це будуть вагові матриці:

>> R_1 = [1 0 0 0 0; 0 2 0 0 0; 0 0 3 0 0; 0 0 0 4 0; 0 0 0 0 5];
>> R_2 = [6 0 0 ; 0 7 0; 0 0 8];

2. Система була побудована раніше, а для пошуку матриці оптимального зворотнього зв'язку та розв'язку відповідного системі матричного рівняння Ріккатті застосовується наступна команда:

>> [F, P, e] = lqr(system, R_1, R_2);

3. Будується відповідна система (з оптимальним зворотним зв'язком) та береться її імпульсна характеристика:

>> A_prime_opt = A — (B*F);
>> system_opt = ss(A_prime_opt, B, C, 0);
>> impulse(system_opt)

Аналіз показує, що побудована система є асимптотично-стійкою

Чому швидкість затухання траєкторії системи з оптимальним зворотнім зв'язком є меншою від тієї, що отримується для системи зі зворотнім зв'язком з дуже низькими полюсами?

Відповідь на це питання викладача досить нетривіальна. Виявляється, що швидкість затухання не є єдиним критерієм побудови системи. Можливо, реальні затрати на досягнення тих величин, що відповідають дуже низьким полюсам, є непомірно великими. Таким чином, завдяки цьому запитанню студент пригадує, що на цьому предметі моделюється поведінка динамічної системи, а отже враховуються реальні характеристики.

Чому по імпульсній характеристиці системи зі зворотним зв'язком неможливо визначити, чи є дана система оптимальною?

Відповідь на це питання у визначеннях імпульсної характеристики та оптимального зворотного зв'язку.

Імпульсна характеристика системи визначає реакцію i-ї компоненти вихідного вектора на імпульс, прикладений до j-ї компоненти вхідного вектора, причому по всіх інших компонентах вхід був нульовим.

Оптимальний зворотній зв'язок може бути знайдений після оптимізації квадратичного критерія (лекційний матеріал).

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