Бороденко В.А. Исследование систем управления в среде MATLAB

Подождите немного. Документ загружается.

251

>> w=tf([10 20 30],[1 2 20 20 2 1]);

>> allmargin(w)

ans =

GainMargin: [0.1225 0.8500 1.8775]

GMFrequency: [0.4744 1.0000 4.2160]

PhaseMargin: 1.6240

PMFrequency: 1.0677

DelayMargin: 0.0265

DMFrequency: 1.0677

Stable: 1

Обычная формулировка критерия Найквиста не позволяет сде-

лать заключение о неустойчивости в замкнутом состоянии системы с

передаточной функцией w=tf(4,[1 2 3 4 5]) (рисунок 4.28), так как

АФЧХ разомкнутой системы на первый взгляд не охватывает точку с

координатами (-1, j0).

Рисунок 4.28

Здесь можно использовать для проверки правило штриховки –

при движении по кривой в направлении увеличения частоты (по

стрелке) наносят штриховку справа от кривой; система устойчива в

замкнутом состоянии, если критическая точка не попадает в заштри-

хованную область. Поскольку на графике кривая изменяет направле-

ние движения, выворачиваясь наизнанку, штриховка из внутренней

области выливается наружу и захватывает критическую точку, систе-

ма будет неустойчивой.

Отметим попутно, что для системы W(s)=4/(s

+2s

+3s

+4s+5)

функция margin() дает бесконечное значение запаса устойчивости Gm

и Pm = 147.1°, хотя одновременно выводится противоречащее этим

запасам предупреждение о неустойчивости системы в замкнутом со-

252

стоянии The closed-loop system is unstable, полученное не из диаграм-

мы Боде, а необъявленным способом.

>> w=tf(4,[1 2 3 4 5])

>> [gm,pm]=margin(w)

Warning: The closed-loop system is unstable.

gm =

Inf

pm =

147.1087

Неясно, почему для системы с внутренней задержкой (рисунок

4.29) делается вывод о невозможности оценить устойчивость Closed

Loop Stable? Not known, хотя АФЧХ не охватывает точку с координа-

тами (-1, j0), запас по амплитуде равен 2.29 дБ и по фазе 180°.

Nyquist Diagram

Real Axis

Imaginary Axis

-1 -0.5 0 0.5 1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

System: sys

Gain Margin (dB): 2.29

At frequency (rad/sec): 0.665

Closed Loop Stable? Not know n

Рисунок 4.29

Особенный случай показан на рисунке 4.30 для системы с пере-

даточной функцией w=tf(1,[1 3 2 6 4 5 2]). Запасы устойчивости про-

грамма версии 7.01 находит и система кажется устойчивой (даже при

использовании правила штриховки), хотя заведомо известно, что в

замкнутом состоянии эта система будет неустойчивой – об этом гово-

рит и вывод Closed Loop Stable? No. При данных условиях думается,

что этот вывод программа делает не с помощью критерия Найквиста,

а иным способом, вероятнее всего самостоятельно замкнув систему и

вычислив ее полюсы. В версии же 7.11 маркеры и оценки устойчиво-

сти системы не выводятся вообще без каких-либо сообщений об

ошибке, но так же программа поступает и в том случае, когда запасы

устойчивости равны 100 % – что же прикажете думать?

253

>> w=tf(1,[1 3 2 6 4 5 2]); % разомкнутая система

>> [num,den]=tfdata(w,'v');

>> roots(num+den) % полюсы замкнутой системы

ans =

-2.9137

0.5109 + 1.1043i

0.5109 - 1.1043i

-0.2378 + 1.0212i

-0.2378 - 1.0212i

-0.6325

а б

Рисунок 4.30

По корням видно, что замкнутая система неустойчива (есть

комплексная сопряженная пара полюсов 0.5109 ± j1.1043 с положи-

тельной вещественной частью), однако из АФЧХ разомкнутой систе-

мы это отнюдь не следует (рисунок 4.30, а), годограф не охватывает

точку с координатами (-1, j0). Кривая Михайлова (рисунок 4.30, б) в

этом отношении предпочтительнее – она правильно указывает на не-

устойчивость замкнутой системы, т. к. при порядке системы n = 6 го-

дограф проходит лишь два квадранта комплексной плоскости, уходя

далее в бесконечность.

Оценим устойчивость этой же системы после замыкания, по-

строив график с помощью функции margin() – заведомо неустойчивая

система характеризуется бесконечными запасами устойчивости по

амплитуде и фазе (рисунок 4.31).

>> w=tf(1,[1 3 2 6 4 5 2]); % разомкнутая система

>> margin(w)

В то же время диаграмма Боде, построенная командой bode(),

при выборе опции Characteristics-Minimum Stability Margins никаких

запасов устойчивости вообще не отображает, поскольку нет частот

пересечения ЛАЧХ с осью абсцисс и ЛФЧХ с линией -180°. Но, как

254

уже говорилось, такое поведение присуще программе и при стопро-

центных запасах устойчивости системы, например, sys=tf(1,[1 2 3]).

-250

-200

-150

-100

-50

Magnitude (dB)

-2

-1

-180

-135

-90

-45

Phase (deg)

Bode Diagram

Gm = Inf dB (at Inf rad/sec) , Pm = Inf

Frequency (rad/sec)

Рисунок 4.31

Посвятим несколько слов оценке устойчивости системы с пере-

даточной функцией W(s)=s/(s^3+2s^2+3s). Хотя по графикам Bode,

Nyquist, Nichols она сверхустойчива после замыкания, никаких выво-

дов по устойчивости графики не делают и показ запасов устойчивости

(опции All/Minimum Stability Margins) не осуществляют. График

функции margin() дает запасы Gm = Inf, Pm = Inf, однако в численном

выводе этой же функции указывается, что замкнутая система неус-

тойчива, т. е. в разных формах представления данных функция проти-

воречит сама себе. Фактически замкнутая система устойчива, по-

скольку переходный процесс устанавливается со значением 0.25, но

нормальных выводов от MATLAB мы не получаем.

Сложно пользоваться значениями запасов устойчивости, если

вывод программы противоречит представляемым характеристикам и

способам определения показателей, изложенным в учебной литерату-

ре, а для систем со 100 % запасами устойчивости значения запасов во-

обще не отображаются. Поскольку сама программа MATLAB не поль-

зуется для оценки критерием Найквиста или диаграммой Боде, а про-

веряет устойчивость иным способом, напрашивается вывод, что для

SS, TF и ZPK моделей критерий Найквиста имеет ограниченное при-

менение лишь для расчета запасов устойчивости, и необходимо вни-

мательно следить за общим заключением об устойчивости системы,

255

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

же время для FRD-модели это единственный способ оценить устойчи-

вость системы после замыкания.

Для оценки устойчивости систем регулирования может исполь-

зоваться и теорема Ляпунова об асимптотической устойчивости. Для

того, чтобы нулевое решение автономной линейной системы





было асимптотически устойчиво, необходимо и достаточно, чтобы для

произвольной положительно определенной матрицы Q существовала

положительно определенная матрица Р, удовлетворяющая уравнению

P + PA = -Q. Часто матрицу Р находят, взяв в качестве Q единич-

ную матрицу.

Пусть матрица собственных коэффициентов MIMO-системы

равна



















, выберем в качестве Q единичную матрицу и вос-

пользуемся штатной функцией X = LYAP(A,Q), решающей уравнение

Ляпунова A*X + X*A' + Q = 0.

>> A=[-1 1; -1 -1]; P=lyap(A, eye(2))

P =

0.5000 0

0 0.5000

>> eig(P)

ans =

0.5000

Квадратичная форма V или соответствующая ей матрица Р явля-

ются положительно определенными в том случае, если собственные

значения λ матрицы Р, которые для симметричной матрицы действи-

тельны, все положительны.

Полученная матрица















500

050

является положительно опре-

деленной, т. к. оба собственных значения 0.5 положительны, следова-

тельно, рассматриваемая линейная система асимптотически устойчива

по Ляпунову.

Функция hsvplot(sys) строит гистограмму сингулярных величин

системы в соответствии с энергетическим вкладом переменных со-

стояния, разделяя их на устойчивые (синие) – три моды, и неустойчи-

вые (красные) – две моды (рисунок 4.32). Ее свойства выводятся опе-

ратором hsvoptions, задаются функцией setoptions().

>> h=hsvplot(tf(4,[1 2 3 4 5]));

>> setoptions(h,'Yscale','log','Offset',0.3)

256

Рисунок 4.32

С помощью MATLAB можно осуществить декомпозицию сис-

темы на две параллельно соединенные устойчивую и неустойчивую

части W(s)=W

(s)+W

(s) посредством функции stabsep(sys). В первую

подсистему объединяются моды с отрицательной действительной ча-

стью, во вторую – моды с нулевой и положительной действительной

частью, обусловливающие нестабильность системы.

>> w=tf(1,[1 1 6 6 25 25]);

>> [w1,w2]=stabsep(w) % декомпозиция системы

Transfer function: % устойчивая часть w1

0.01875 s^2 + 0.0875 s + 0.1938

-------------------------------

s^3 + 3 s^2 + 7 s + 5

Transfer function: % неустойчивая часть w2

-0.01875 s + 0.00625

--------------------

s^2 - 2 s + 5

>> eig([w w1 w2]) %распределение полюсов по подсистемам w, w1, w2

ans =

1.0000 + 2.0000i -1.0000 + 2.0000i 1.0000 + 2.0000i

1.0000 - 2.0000i -1.0000 - 2.0000i 1.0000 - 2.0000i

-1.0000 + 2.0000i -1.0000

-1.0000 - 2.0000i

-1.0000

>> ww=w1+w2 % обратная композиция системы для проверки точности

Transfer function:

257

3.469e-018 s^4 + 1.388e-017 s^3 - 2.498e-016 s^2 + 2.776e-016 s + 1

-------------------------------------------------------------------

s^5 + s^4 + 6 s^3 + 6 s^2 + 25 s + 25

избавляясь от членов полинома числителя со значением менее 1e-15

>> [num,den]=tfdata(ww,’v’); nn=num.*(abs(num)>1.e-15);

>> ww=tf(nn,den)

Transfer function:

-------------------------------------

s^5 + s^4 + 6 s^3 + 6 s^2 + 25 s + 25

Оператор stabsepoptions выводит свойства stabsep(), в числе ко-

торых имеется возможность представления неустойчивыми полюсов,

находящихся от мнимой оси на расстоянии меньшем, чем отступ off-

set (по умолчанию 1е-08). Например, система, находящаяся практиче-

ски на границе устойчивости (имеется два мнимых корня 0 ± 1.4142i),

за счет погрешности вычислений признана практически устойчивой.

>> q=tf([1 2 3],[1 2 3 4 2]);

>> [q1,q2]=stabsep(q)

Transfer function: % устойчивая часть q1

s^2 + 2 s + 3

-----------------------------

s^4 + 2 s^3 + 3 s^2 + 4 s + 2

Transfer function: % неустойчивая часть q2

Зададим допустимое расстояние до мнимой оси 0.001, теперь

декомпозиция на устойчивую и неустойчивую части производится.

>> [q1,q2]=stabsep(q,'Offset',0.001)

Transfer function: % устойчивая часть q1

0.4444 s + 1.111

----------------

s^2 + 2 s + 1

Transfer function: % неустойчивая часть q2 с мнимыми полюсами

-0.4444 s + 0.7778

----------------------

s^2 + 1.166e-015 s + 2

258

4.4 Оценка качества процесса регулирования

Анализ качества регулирования в MATLAB для обычных ли-

нейных систем производят с помощью корневых и прямых оценок.

Если, вызвав окно графического анализатора ltiview(sys), вы-

брать в контекстном меню графиков тип графика Plot Types-Pole/Zero,

для одномерной системы будет выведена карта размещения нулей

Zero (кружки) и полюсов Pole (крестики) системы на комплексной

плоскости (рисунок 4.33). Аналогичные действия выполняет опция

I/O Pole/Zero или функции pzmap(sys), pzplot(sys).

>> sys=tf([1 2 3],[1 2 3 4]); pzplot(sys)

Рисунок 4.33

Из основных показателей качества выводится только значение

перерегулирования Overshoot, которое дается и для нулей, хотя и не

имеет в этом случае смысла; оно отображается для каждого полюса

Pole отдельно, а не для системы в целом. Время регулирования нужно

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

зоне 5 % t

рег

≈ 3/0.175 = 17.14 с.

Для быстрой оценки основных показателей качества предлага-

ется функция pmetrics(sys, delta), где delta – задаваемая зона Δ контро-

ля времени регулирования, по умолчанию равная 0.05. Анализируют-

ся только полюсы с отрицательной действительной частью в порядке

удаления от мнимой оси.

function pmetrics(w, delta)

% Root metrics of step function

259

if nargin==1 % если зона delta не задана, будет 0.05

delta = 0.05;

end

p = esort(pole(w)); idn=0;

for id=1:length(p)

if real(p(id))<0 & imag(p(id))<=0

idn = idn + 1; ep(idn) = p(id);

end

ep = ep'; tr = -log(1/delta)./real(ep);

ov = 100*exp(pi.*real(ep)./abs(imag(ep)));

Poles______________SettlingTime_____Overshoot = [ep tr ov]

end

Выводимый результат

>> v=[0 1 -0.1-1i -0.1+1i -0.5+1i -0.5-1i -0.2 -10];

>> w=zpk([],v,1); pmetrics(w)

Poles_______________SettlingTime_____Overshoot =

-0.1000 + 1.0000i 29.9573 73.0403

-0.2000 14.9787 0

-0.5000 + 1.0000i 5.9915 20.7880

-10.0000 0.2996 0

Корневой годограф одномерной системы вызывается командой

rlocus(sys), rlocus(sys, K) или rlocus(sys, [k1 k2 ...]), если задан вектор

значений К или необходим просмотр годографа для заданных значе-

ний коэффициента k1 k2 и т. д., например

>> rlocus(syso,0:0.1:10),grid

Обратите внимание на невозможность задания вектора k край-

ними значениями в фигурных скобках, как в частотных функциях.

Команды [R,K]=rlocus(sys), R=rlocus(sys) или R=rlocus(sys,K)

выводят вектор R комплексных полюсов системы на заданных значе-

ниях коэффициента k. Отдельно полюса системы могут быть получе-

ны командой P=rlocus(sys,0), нули – командой Z=rlocus(sys,Inf). Вари-

антом с расширенными возможностями является функция rlocusplot().

Масштабная сетка содержит окружности равных значений соб-

ственных частот (частот сопряжения) Natural Frequencies с центром в

начале координат и лучи равных значений коэффициента демпфиро-

вания Damping, у которых значение 1 совпадает с отрицательной дей-

ствительной полуосью, а значение 0 – с мнимой осью (рисунок 4.34).

Ни то, ни другое не говорит прямо об основных показателях качества,

например, перерегулирование связано с коэффициентом демпфирова-

260

ния ξ зависимостью )1/exp(



 , т. е. требует пересчета.

Рисунок 4.34

Для каждой точки годографа можно получить с помощью мыши

текущее значение коэффициента передачи Gain, значение корня Pole,

значения коэффициента демпфирования Damping и частоты сопряже-

ния Frequency (рад/с), значение перерегулирования Overshoot (%).

Время регулирования придется определять самостоятельно по форму-

ле k/|α|, где α – это действительная часть корня, k = 3 для зоны Δ = 5 %

и k = 4 для Δ = 2 %. Все показатели определяются по действительной a

и мнимой b частям, выполним, например, расчет для правого полюса

>> a=10.1; b=-18.7;

>> tr=3/abs(a) % время регулирования для Δ = 5 %, с

tr =

0.2970

>> sigma=exp(pi*a/abs(b)) % перерегулирование Overshoot, о. е.

sigma =

5.4564

>> wn=sqrt(a*a+b*b) % частота сопряжения Frequency ω

, рад/с

wn =

21.2532

>> ksi=-a/wn % коэффициент демпфирования Damping ξ

ksi =

-0.4752

У неустойчивой системы действительная часть хотя бы одного

полюса положительна, а значение коэффициента демпфирования это-