Приклад розрахунку робастного контролера (H-infinity control)

Що таке робастний контролер і навіщо ускладнювати собі життя? Чому нас не влаштовує стандартний, всіма впізнаваний, ПІД-регулятор?

Відповідь криється в самій назві, з англ. «robustness» — The quality of being strong and not easy to break (Властивість бути сильним і складно зламати). У випадку з контролером це означає, що він повинен бути «жорстким», непіддатливою до змін об'єкта управління. Наприклад: мат. моделі DC мотора є 3 основних параметри: опір і індуктивність обмотки, і постійні Кт Ке, які рівні між собою. Для розрахунку класичного ПІД регулятора, дивляться в даташит, беруть ті 3 параметри і розраховують коефіцієнти ПІД, начебто все просто, що ще потрібно. Але мотор — це реальна система, в якій ці 3 коефіцієнта не постійні, наприклад в слідстві високочастотної динаміки, яку складно описати чи потрібно високий порядок системи. , Наприклад: Rдаташит=1 Ом, а насправді R знаходитися в інтервалі [0.9,1.1] Ом. Так у показники якості у випадку з ПІ регулятором можуть виходити за задані, а робастний контролер враховує невизначеності і здатний утримати показники якості замкнутої системи в потрібному інтервалі.

На цьому етапі виникає логічні питання: А як знайти цей інтервал? Знаходиться він за допомогою параметричної ідентифікації моделі. На Хабре нещодавно описали метод МНК («Параметрична ідентифікація лінійної динамічної системи»), проте дає одне значення идентифицируемого параметра, як в даташіте. Для знаходження інтервалу значень ми використовували «Sparse Semidefinite Programming of Relaxation Поліноміальні Optimization Problems» додаток для матлаба. Якщо буде цікаво, можу написати окрему статтю як користуватись даним пакетом та як провести ідентифікацію.

Сподіваюся, зараз стало хоч не багато зрозуміліше навіщо робастное управління потрібно.

Не буду особливо вдаватися в теорію, тому що я її не дуже зрозумів; я покажу етапи, які необхідно пройти, щоб отримати контролер. Кому цікаво, можна почитати «Essentials Of Robust Control, Kemin Zhou, John C. Doyle, Prentice Hall» або документацію Matlab.

Ми будемо користуватися наступною схемою системи управління:

image
Рис 1: Блок схема системи управління

Думаю, тут все зрозуміло (di – обурення).

Постановка завдання:
Знайти: Gc(s) і Gf(s) які задовольняють всі задані умови.

Дано: image– об'єкт управління з невизначеностями, які задані проміжки часу: image
Для подальших розрахунків будемо використовувати номінальні значення, а невизначеності врахуємо далі.
Kn=(16+9)/2=12.5,p1n=0.8,p2n=2.5, відповідно отримав номінальний об'єкт управління:
image
image
Типи збурень:
  • image— ступінчасте збурення підсилювача з амплітудою D_a0
  • image— синусоїдальна обурення об'єкта управління з амплітудою a_p і частотою w_p
  • image— синусоїдальна обурення сенсора з амплітудою a_s і частотою w_s
Або іншими словами це характеристики шуму елементів системи.

Робочі характеристики (performance specifications)
  1. Коефіцієнт посилення системи управління Kd=1
  2. Стале помилка при вхідному впливі Ramp R0=1: image

  3. Стале помилка в присутності da: image
  4. Стале помилка в присутності dp: image
  5. Стале помилка в присутності ds: image
  6. Перерегулювання image
  7. Час регулювання image
  8. Час наростання image


Рішення
Дуже детально розписувати не буду, тільки основні моменти.

Одним з ключових кроків H∞ методі визначення вхідних і вихідних вагових функцій. Ці вагові функції використовуються для нормування входу і виходу і відображення тимчасових і частотних залежностей вхідних збурень і робочих характеристик вихідних змінних (помилок) [1]. Чесно кажучи, мені це практично ні про що не говорить, особливо якщо вперше стикаєшся з даним методом. У двох словах, вагові функції використовуються для завдання потрібний властивостей системи управління. Наприклад, у зворотному зв'язку варто сенсор, який має шум, зазвичай високочастотний, так от, вагова функція буде свого роду кордоном, яку контролер не повинен перетинати, що б відфільтрувати шум сенсора.

Нижче будемо виводити ці вагові функції на основі робочих характеристик.

  1. imageТут все просто
  2. В цьому пункті нам потрібно визначити скільки полюсів на нулі необхідно мати для контролера (позначимо µ) щоб задовольнити Умова 2. Для цього скористаємося таблицею: image
    Рис 2: Помилки щодо положення, швидкості та прискорення

    System type або астатізм (p) – позначає кількість полюсів в нулі, об'єкта управління, в нашому випадку p=1 (система з астатизмом першого порядку). Так як наш об'єкт управління і є система з астатизмом 1 порядку, у такому разі полюсів в нулі у контролера бути не повинно. Скористаємося формулою:
    μ+p=h
    Де h – порядок вхідного сигналу, для ramp h=1;
    μ=h-p=1-1=0
    Тепер скористаємося теоремою кінцевих значень, щоб знайти e_r∞
    image
    Де e_r – помилка стеження,
    image
    yr – дійсний вихід (див. Рис 1), уd – бажаний вихід
    image
    Рис 3: Визначення помилки спостереження

    В результаті отримаємо таке:
    image
    Це формула усталеної помилки, невідома в даному рівнянні S (Sensitivity function)
    image
    Де image– sensitivity function, L(s) – loop function. Без поняття як вони переводяться на російську, залишу англійські назви. Також complementary sensitivity function image(як видно з формул, в склад S і T входять Gc – розраховується контролер, відповідно кордону для S і T ми знайдемо з помилок і робочих характеристик, вагові функції визначимо з S і T, а матлаб з вагових функцій знайде бажаний контролер.)
    В двох словах про S і T [1].
    • Sensitivity function, S(s), описує вихід y(s) як функцію обурення da і dp, також пов'язує помилку стеження і вхідний вплив (для низьких частот) image
    • Complementary sensitivity function, T(s), пов'язує вихід системи з вхідним впливом, а також показує на скільки шум сенсора ds впливає на вихід системи (для високих частот).
      image
    image
    Рис 4: Діаграма Боде для S,T і L

    З графіка видно що S послаблюємо низькочастотні збудження, у той час як Т послаблює високочастотні обурення.
  3. Скористаємося теоремою кінцевих значень для e_da
    image
    image
    Так як у нас вийшло два нерівності знайдемо умову задовольняє обох:
    image
    Це умова говорить, де Sensitivity function повинна перетинати 0 dB axis діаграми Боде.
    image
    Рис 5: Боде для S

  4. Dp у нас гармонійне обурення, низькочастотне. Створимо маску для S в районі частот wp
    image
    Це значення показує що S повинен знаходитись нижче -32 dB для частот wp, що б отфильровать обурення
    image
    Рис. 6: Маска для S

  5. Ds також гармонійне обурення високих частот, тут свою роль буде виконувати T
    Виконаємо теж саме:
    image
    image
    Рис 7: Маска для Т

    Комплексний порядок вагових функцій визначається з умови маски і частоти перетину з 0 віссю. У нашому випадки від wp до w приблизно одна декада, а так як у нас є -32 dB то S повинен бути мінімум 2 порядку. Теж саме стосується і Т.
    У підсумку схематично обмеження мають наступний вигляд для S і Т відповідно:
    imageimage
    Рис 8: Всі маски для S і T

  6. Для перекладу часових характеристик скористаємося графіками
    image
    Рис 9: Графік залежності коефіцієнту демпфування від переригулирования

    Знаючи переригулирование знайдемо коефіцієнт демфирования для 10% ->
    ξ=0.59
    Знаючи коефіцієнт демфирования знайдемо (резонансне) максимально припустиме значення S і T
    image
    Рис 10: Графік залежності S_p0 і T_p0 від коефіцієнту демфирования

    S_p0=1.35
    T_p0=1.05

  7. З часу регулювання і настроювання знайдемо наскільки швидка повинна бути система управління
    image
    Далі знайдемо натуральну частоту для S
    image
    image
    Натуральна частота для T знаходимо з діаграми Боде. За умовою 5 в частоті 40 rad/s T повинен знаходитись нижче -46 dB, а значить при нахилі в -40 дб/дек натуральна частота повинна знаходитися нижче 4 rad/s. Ладу Боде, підбираємо оптимальне значення, я отримав що:
    image
    image
    Рис 11: Боде T-функції

    Після цього у нас є всі дані для побудови S T, які потім трансформуємо у вагових функцій. S T мають наступну форму:
    image
    Зазвичай для побудови вагових функцій використовують Баттерворта коефіцієнти
    image
    Вагові функції мають вигляд:
    image
    imageвсе просто, і нас є все необхідне, що б підставити в формулу. imageпотрібно ще кілька розрахунків.
    Так як робочі характеристики дають нам граничні умови в рамках яких наш контролер буде виконувати всі умови. Вагові функції об'єднують в собі всі умови і потім використовуються як права і ліва межа в яких перебувати контролер задовольняє всі умови.
    image
    • imageцей параметр називають generalized DC gain, він відображає поведінки для низьких частот (s=0)

    • image, w1 вибираємо приблизно біля частоти збурення або на одну декаду нижче (поставимо в 0,0005 rad/s)

    • • LMI solver не сприймає нулі функції(zero in origine), тому s замінимо на нуль біля початку координат (s+0.0005)


В результаті отримаємо:
image
Generalized plant
Hinf метод або метод мінімізації нескінченної норми Hinfinity відноситься до загальної формулюванні проблеми управління яка грунтується на схемі подання системи управління зі зворотним зв'язком:

image
Рис 12: Узагальнений схема системи управління

Перейдемо до пункту розрахунку контролера і що для цього потрібно. Вивчення алгоритмів розрахунку нам не пояснювали, говорили: «робіть ось так і все вийде», але принцип цілком логічний. Контролер виходить в ході вирішення оптимізаційної задачі:
image
image– замкнута передатна функція W до Z.
Зараз нам потрібно скласти generalized plant (пунктирний прямокутник рис нижче). Gpn ми вже визначили, це номінальний об'єкт управління. Gc — контролер, який ми отримаємо в результаті. W1 = Ws(s), W2 = max (WT(s),Wu(s)) – це вагові функції визначені раніше. Wu(s) – це вагова функція невизначеності, давайте визначимо її.
image
Рис 13: Розкрита схеми управління

Wu(s)
Припустимо що маємо мультиплікативні невизначеності в об'єкті управління, можна зобразити так:
image
Рис 14: Мультиплікативну невизначеність

І так для знаходження Wu скористаємося матлаб. Нам потрібно побудувати Bode все можливих відхилень від номінального об'єкта управління, і потім побудувати передатну функцію яка опише всі ці невизначеності:
image
Зробимо приблизно за 4 проходи для кожного параметра і побудуємо bode. В результаті отримаємо таке:
image
Рис 15: Диграмма Боде невизначеностей

Wu буде лежати вище цих ліній. У матлаб є tool який дозволяємо мишкою вказати точки і по цим точкам будує передатну функцію.

Кодmf = ginput(20);
magg = vpck(mf(:,2),mf(:,1));
Wa = fitmag(magg);
[A,B,C,D]=unpck(Wa);
[Z,P,K]=ss2zp(A,B,C,D);
Wu = zpk(Z,P,K)


Після введення точок старіє крива і пропонується ввести ступінь передатної функції, введемо ступінь 2.
Ось що вийшло:
image
Тепер визначимо W2, для цього побудуємо Wt і Wu:
З графіка видно що Wt більше Wu значить W2 = Wt.
image
Рис 16: Визначення W2

Далі нам потрібно в симулинке побудувати generalized plant як зроблено нижче:
image
Рис 17: Блок схема Generalized plant в simulink

І зберегти під ім'ям, наприклад g_plant.mdl
Один важливий момент:
image— not proper tf, якщо ми її залишимо так то нам видасть помилку. З цього замінюємо на imageі потім додамо два нуля до виходу z2 з допомогою “sderiv“.
Реалізація в матабp = 2; % частота нулів в районі частоти Т функції
[Am,Bm,Cm,Dm] = linmod('g_plant');
M = ltisys(Am,Bm,Cm,Dm);
M = sderiv(M,2,[1/p 1]);
M = sderiv(M,2,[1/p 1]);
[gopt,Gcmod] = hinflmi(M,[1 1],0,0.01,[0,0,0]);
[Ac,Bc,Cc,Dc] = ltiss(Gcmod);
[Gc_z,Gc_p,Gc_k] = ss2zp(Ac,Bc,Cc,Dc);
Gc_op = zpk(Gc_z,Gc_p,Gc_k)

Після виконання цього коду отримуємо контролер:
image
В принципі можна так і залишити, але зазвичай видаляють низько — та високочастотні нулі і полюси. Таким чином ми зменшуємо порядок контролера. І отримуємо наступний контролер:
image
отримаємо ось такий Hichols chart:
image
Рис 18: Hichols chart розімкнутої системи з отриманим контролером

І Step response:
image
Рис 19: Перехідна характеристика замкненої системи з контролером

А тепер саме солодке. Вийшов наш контролер робастным чи ні. Для цього експерименту просто потрібно змінити наш об'єкт управління (коефіцієнти к, р1, р2) і подивитися step response та потрібними характеристиками, в нашому випадку це перерегулювання, час регулювання для 5 % і час наростання [2].
imageimageimage
Рис 20: Тимчасові характеристики для різних параметрів об'єкт управління

Побудувавши 20 різних перехідних характеристик, я виділив максимальні значення для кожної часової характеристики:
• Максимальне переригулирование – 7.89 %
Макс час наростання – 2.94 сек
Макс час регелирования 5% — 5.21 сек
І о диво, характеристики там де потрібно не тільки для номінального об'єкта, але і для інтервалу параметрів.
А зараз порівняємо з классичиским ПІД контролером, і подивимося коштувала гра свічок чи ні.

ПІД расчитал pidtool для номінального об'єкта управління (див вище):

image
Рис 21: Pidtool

Отримаємо такий контролер:
image

Тепер H-infinity vs PID:

image
Рис 22: H-infinity vs PID

Видно що ПІД не справляється з такими невизначеностями і ПХ виходить за задані обмеження, в той час як робастні контролер «жорстко» тримають систему в заданих інтервалах перерегулювання, часу наростання і регулювання.

Що б не подовжувати статтю і не втомлювати читача, опущу перевірку характеристики 2-5 (помилки), скажу, що в разі робастного контролера всі помилки знаходяться нижче заданих, також був проведений тест з іншими параметрами об'єкта:

image

Помилки перебували нижче заданих, що означає, що даний контролер повністю справляється з поставленим завданням. У той час як ПІД не впорався лише з пунктом 4 (помилка dp).

На цьому все за розрахунком контролера. Критикуйте, запитуйте.

Посилання на файл матлаб: drive.google.com/open?id=0B2bwDQSqccAZTWFhN3kxcy1SY0E

Список літератури

1. Guidelines for the Selection of Weighting Functions for H-Infinity Control
2. it.mathworks.com/help/robust/examples/robustness-of-servo-controller-for-dc-motor.html
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.