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

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

241

4.8722.10908.11928.86.1

234







ssss

и строим график (рисунок 4.18) последовательностью команд

>> k=tf(-[1.6 8.928 11.908 10.22 87.4],[1.04 1.67 2.2 1.8]);

>> p=nyquistoptions;

>> p.Title.String='D-разбиение';

>> nyquist(k,p)

Оставляем на графике кривые как для положительных, так и для

отрицательных частот.

Рисунок 4.18

Значения k, соответствующие устойчивости системы, лежат в

области D0 на отрезке действительной оси (Real Axis) от отметки 47.1

до плюс бесконечности, что требует проверки. Подставив значение

координаты точки на этом отрезке k = 60 в характеристическое урав-

нение, убеждаемся, что это действительно область устойчивости, так

как все корни уравнения имеют отрицательную действительную

часть. Аналогично проверяются области D1, D2 при значениях k = -60

и k = 0.

>> k=60;

>> D0=[1.6 8.928+k*1.04 11.908+k*1.67 10.22+k*2.2 87.4+k*1.8]

D0 =

1.6000 71.3280 112.1080 142.2200 195.4000

>> roots(D0)

ans =

-42.9970

-1.5063

242

-0.0384 + 1.3726i

-0.0384 - 1.3726i

Рассмотрим пример D-разбиения по двум параметрам.

Для построения методом D-разбиения области устойчивости в

пространстве параметров Т и k системы четвертого порядка с характе-

ристическим уравнением D(s)=Ts

+ (2T + 1)s

+ (2T+2)s

+ 2s + k = 0

воспользуемся достаточным условием Δ

= a

- a

) - a

≥ 0

критерия Гурвица и сформируем по нему функцию Z(T, k), подставив

значения параметров. При вычислениях заменяем Т на х и k на y, что-

бы не объявлять символьные переменные, поскольку переменные x, y

и z для двумерных функций MATLAB принимает по умолчанию. За-

дав пределы изменения постоянной времени 0 ≤ Т ≤ 2 и коэффициента

усиления 0 ≤ k ≤ 4 с шагом 0.1, построили с помощью функции con-

tourf() график с залитой зеленым цветом областью устойчивости D(0)

(рисунок 4.19).

D-разбиение по двум параметрам

Постоянная времени Т

Коэффициент k

0 0.5 1 1.5 2

Рисунок 4.19

>> [x,y]=meshgrid(0:0.1:2,0:0.1:4); % сетка значений

>> z=2*((2*x+1).*(2*x+2)-x*2)-(2*x+1).^2.*y; % функция

% построить контурный график и залить область, выделенную на уровне 0-0

>> [c,h]=contourf(x,y,z,[0 0]);

>> set(h,'showtext','on'); % включить печать значений уровней

>> title('D-разбиение по двум параметрам'); grid;

>> xlabel('Постоянная времени Т'); ylabel('Коэффициент k')

График не учитывает требований, вытекающих из необходимых

условий устойчивости a

> 0, а именно: Т > 0 – вытекает из требования

= T > 0, и k > 0 – обусловлено требованием a

= k > 0. Эти условия

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

Если необходима оценка устойчивости по Найквисту, то следует

задавать для анализа передаточную функцию разомкнутой (Open-

Loop) системы, а на графике оценивать положение особой точки с ко-

243

ординатами (-1, j0) – АФЧХ неустойчивой в замкнутом состоянии

(Closed-Loop) системы охватывает эту точку. АФЧХ выводится опци-

ей Plot Types-Nyquist в контекстном меню графоанализатора ltiview()

или отдельной командой nyquist(), nyquistplot(). Для правильного

представления углов желательно использовать команду axis equal.

График более понятен, если убрать кривую для отрицательных частот.

При выборе в контекстном меню графика опций Characteristics-

Stability Margins в точках, отмеченных маркерами, указываются запа-

сы устойчивости по амплитуде Gain Margin в децибелах, по фазе Phase

Margin в градусах, по запаздыванию Delay Margin в секундах и соот-

ветствующие им частоты в рад/с, а также дается общая оценка устой-

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

ной связью. Значение Inf говорит о том, что по данному параметру

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

АФЧХ первого рода) или нулевые значения – об отсутствии запасов.

Обратите внимание на то, что в СНГ максимальным значением запаса

по фазе (100 %) считается не бесконечность как в MATLAB, а 180°.

После перехода от All Stability Margins к Minimum Stability Margins

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

менений не произойдет.

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

зомкнутой системы можно получить одним из следующих способов:

- вычесть соответственно из коэффициентов полинома знамена-

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

ентов полинома числителя;

- разомкнуть систему, созданную в среде Simulink: реально, уб-

рав (Cut) главную обратную связь и подключив к точке разрыва выход

Out1, либо виртуально, описанным путем, после чего передать пара-

метры линеаризированной модели в рабочее пространство MATLAB.

Первый способ может применяться только в случае охвата всей

системы главной единичной отрицательной обратной связью, второй

более универсален.

АФЧХ разомкнутой системы без нулей w=tf(3,[1 2 3 2]) прохо-

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

оси, столько квадрантов, каков порядок системы (рисунок 4.20). Она

не охватывает точку с координатами (-1, j0), следовательно, система в

замкнутом состоянии будет устойчивой. Об этом гласит и заключение

по каждому из запасов Closed Loop Stable? Yes.

В СНГ запас устойчивости принято считать равным расстоянию

от критической точки до точки пересечения АФЧХ с отрицательной

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

= 1 – 0.747= 0.253. На

244

диаграмме Найквиста запас по модулю Gm (дБ) вычисляется по от-

ношению единицы к отсекаемому на указанной полуоси отрезку

>> Gm=20*log10(1/0.747)

Gm =

2.5336

Рисунок 4.20

Запас устойчивости Dm по чистой задержке в секундах можно

вычислить, используя значения запаса по фазе Pm = 17.1° и соответст-

вующей ему частоты FPm = 1.56 рад/с (с пересчетом в градусы)

>> Dm=17.1/1.56/180*pi

Dm =

0.1913

Особенность АФЧХ нейтральной системы w=tf(5,[1 2 3 0]),

имеющей в разомкнутом состоянии нулевой полюс (рисунок 4.21) –

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

Рисунок 4.21

245

Распространенным инструментом оценки устойчивости систем

является диаграмма Боде или ЛЧХ, вызов которой производится из

командной строки обращением bode(sys), bode(num,den), bode(a,b,c,d),

bode(sys1,sys2,...), с указанием крайних частот bode(sys,{wmin,

wmax}) или всего диапазона частот bode(sys, logspace(-2,2,100)),

bode(sys, fr), из LTI Viewer опцией меню Plot Types-Bode. Дополни-

тельными возможностями обладает функция bodeplot(). Выведем ЛЧХ

для той же системы (рисунок 4.22) с отмеченными запасами устойчи-

вости – опция Characteristics-Minimum Stability Margins.

Bode Diagram

Frequency (rad/sec)

-150

-100

-50

System: w

Gain Margin (dB): 2.5

At frequency (rad/sec): 1.73

Closed Loop Stable? Yes

Magnitude (dB)

-2

-1

-270

-180

-90

System: w

Phase Margin (deg): 17.1

Delay Margin (sec): 0.192

At frequency (rad/sec): 1.56

Closed Loop Stable? Yes

Phase (deg)

Рисунок 4.22

Диаграмма Николса является полной аналогией годографа

Найквиста и предназначена для тех же целей. Она имеет три отличия:

- если годограф Найквиста строится в координатах действитель-

ная часть-мнимая часть на комплексной плоскости, то годограф Ни-

колса строится в дополняющих их координатах фаза (градусы)-

амплитуда (децибелы);

- положение характеристики оценивается относительно точки с

координатами (0 дБ, -180°);

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

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

Для работы с экспериментальными данными диаграмма Никол-

са удобнее годографа Найквиста, так как обычно в опытах произво-

дится замер отношения амплитуд и разности фаз гармонических сиг-

налов на входе и выходе объекта. Построим оба графика для одной и

той же разомкнутой системы (рисунок 4.23).

246

Рисунок 4.23

Оценки запасов устойчивости полностью совпадают. Все про-

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

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

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

На диаграмме Николса (рисунок 4.24) неустойчивой в замкну-

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

тойчивости по амплитуде отрицателен, т. е. отсутствует.

Рисунок 4.24

247

Замкнутая система устойчива, если годограф Николса разомк-

нутой системы проходит ниже точки с координатами (0 дБ, -180°). За-

пас устойчивости по амплитуде (Gain Margin) равен расстоянию от

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

фазе (Phase Margin) – расстоянию от указанной точки до кривой по

горизонтали вправо. В этой же точке рассчитывается запас по чистой

задержке (Delay Margin).

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

функция margin(), отображающая минимальные из возможных запасы

устойчивости (в одну сторону), и allmargin(), осуществляющая вывод

всех значений запасов устойчивости в виде структуры. При обраще-

нии [Gm, Pm, GmF, PmF] = margin() данная функция возвращает наи-

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

Gm и фазе Pm (в градусах) и соответствующие им частоты GmF, PmF,

а при отсутствии левой части margin() строит диаграмму Боде с вы-

численными значениями запасов в заголовке.

>> w=tf(3,[1 2 3 2]); % берем рассмотренную ранее систему

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

gm =

1.3338

pm =

17.1340

gmf =

1.7323

pmf =

1.5599

>> S=allmargin(w)

S =

GainMargin: 1.3338 % запасы устойчивости по амплитуде

GMFrequency: 1.7323 % частоты пересечения линии -180° (рад/с)

PhaseMargin: 17.1340 % запасы по фазе (град)

PMFrequency: 1.5599 % частоты среза (рад/с)

DelayMargin: 0.1917 % запасы по задержке (с)

DMFrequency: 1.5599 % частоты среза (рад/с)

Stable: 1 % система после замыкания устойчива

(иначе 0 или NaN для FRD систем)

В функциях allmargin() и margin() версии 7.11 Gm – это относи-

тельная величина К=1/А

-180°

=1/0.747=1.3338, которая вычисляется по

отношению единицы к отсекаемому АФЧХ на отрицательной дейст-

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

функцией margin(), отображается другое значение Gm, равное 20lg(K).

В margin() версии 7.01 Gm = 20lg(K)=2.5336 и измеряется в дБ.

248

В форме [gm,pm,gmf,pmf] = margin(mag,phase,fr) или S = allmar-

gin(mag,phase,fr) функции находят запасы устойчивости и соответст-

венные частоты по готовым массивам значений амплитуды mag и фа-

зы phase на диапазоне частот fr. Для повышения точности расчетов в

вычислении промежуточных точек применяется интерполяция.

Функция loopmargin() рассчитывает запасы устойчивости от-

дельно по всем контурам обратной связи для MIMO систем с квадрат-

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

в форме cm = loopmargin(sys) для SISO систем.

>> sys=ss([0 1 0;0 0 1;-4 -3 -2],[0 0 1;0 0 2]',[1 0 0; 0 6 0],0);

>> [cm,dm,mm]=loopmargin(sys)

cm = % классические запасы устойчивости

2x1 struct array with fields:

GainMargin

GMFrequency

PhaseMargin

PMFrequency

DelayMargin

DMFrequency

Stable

dm = % круговые запасы устойчивости

2x1 struct array with fields:

GainMargin

PhaseMargin

Frequency

mm = % допустимые изменения коэффициента усиления и фазы

GainMargin: [0.5945 1.6820]

PhaseMargin: [-28.5339 28.5339]

Frequency: 3.9288

>> cm.GainMargin

ans =

26.0001

ans =

Inf

Исследуем ту же систему w=tf(3,[1 2 3 2]) при введенной транс-

портной задержке 2 с, существенно превышающей показанный на

графике nyquist допустимый запас Delay Margin 0.192 с.

>> set(w,'outputdelay',2);

>> nyquist(w)

Параметры системы изменились настолько (рисунок 4.25), что

она стала неустойчивой, запасы по фазе и модулю отрицательные. Ха-

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

АФЧХ спиралью вокруг начала координат.

249

Рисунок 4.25

MATLAB позволяет просмотреть сразу несколько графиков и

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

системы. Например (рисунок 4.26), при значениях коэффициента уси-

ления разомкнутой системы k = 0.5, 1.0 и 2.0 замкнутая система будет

соответственно устойчивой, на границе устойчивости и неустойчивой

(объединяем значения вариантов коэффициента полинома числителя

фигурными скобками).

Рисунок 4.26

>> sys=tf({0.5;1;2},[1 2 3 4 1])

Transfer function from input to output...

0.5

#1: -----------------------------

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

250

#2: -----------------------------

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

#3: -----------------------------

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

>> nyquist(sys)

Иногда приходится менять масштаб графика, чтобы проверить

устойчивость системы. Так, при автоматическом выводе годографа

Найквиста для системы с ПФ w=tf([10 20 30],[1 2 20 20 2 1]) участок с

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

а). После изменения масштаба с помощью опции меню Zoom In или

опции контекстного меню Properties-Limits на пределы -10 … 0.5 по

оси Х и -0.3…0.3 по оси У становится очевидным, что система устой-

чива, так как кривая обходит критическую точку (рисунок 4.27, б).

Другим способом детального исследования части характеристики яв-

ляется отказ от автоматического задания диапазона частот и само-

стоятельный его выбор, в этом случае кривая бывает обычно более

сглаженной по сравнению с первым способом.

а б

Рисунок 4.27

Это пример системы с двумя видами запаса устойчивости по

амплитуде – на увеличение 5.47 дБ и на уменьшение коэффициента

усиления -1.41 дБ. Он показывает, что отрицательное значение запаса

устойчивости по амплитуде в MATLAB еще не является признаком

неустойчивости системы – это может быть просто запас на уменьше-

ние коэффициента усиления, что обязательно следует анализировать.

Более того, функция allmargin() выводит под названием GainMargin

все три пересечения АФЧХ с линией -180° на частотах 0.4744, 1.000 и

4.216 рад/с, даже то (на частоте 0.4744), которое запас устойчивости

по амплитуде не характеризует и не отвечает смыслу GainMargin.