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

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

221

дроби находим, определив числитель num и знаменатель den переда-

точной функции

>> v=[-.3 -1+i -1-i -2+0.5i -2-0.5i];

>> q=zpk([],v,3);

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

>> [r,p,k]=residue(num,den)

r =

-0.6056 + 0.5821i

-0.6056 - 0.5821i

0.2850 + 0.5390i

0.2850 - 0.5390i

0.6412

p =

-2.0000 + 0.5000i

-2.0000 - 0.5000i

-1.0000 + 1.0000i

-1.0000 - 1.0000i

-0.3000

k =

[]

Расчет значений оригинала в интервале 0-t

уст

(25 с) с шагом 0.1 с

и построение графика для двух пар комплексных сопряженных полю-

сов и одного действительного полюса

>> t=0:0.1:25;

>> y = .6412*exp(-.3*t)+...

2*exp(-2*t).*(-.6056*cos(.5*t)-.5821*sin(.5*t))+...

2*exp(-1*t).*(.285*cos(1*t)-.539*sin(1*t))

>> plot(t,y), grid, title 'Оригинал передаточной функции'

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

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

ниям должна совпадать с импульсной.

Разложение на простые дроби использует предлагаемая функция

w=resp(sys, tf), аргументами которой являются SISO объект sys и дли-

тельность моделирования tf, возвращает эта функция изображения

(передаточные функции) для каждой составляющей оригинала изо-

бражения sys и графики отдельных составляющих и их суммы.

function [w] = resp(sys, time)

% Построение составляющих оригинала ПФ SYS за время time

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

[r,p,k]=residue(num,den);

id=0;

for ik=1:length(r)

222

if imag(p(ik))>0 % комплексный корень

[num,den]=residue([r(ik) r(ik)'],[p(ik) p(ik)'],[]);

id=id+1; w(id)=tf(num,den);

end

if isreal(p(ik)) % действительный корень

id=id+1; w(id)=tf(r(ik), [1 -p(ik)]);

end

clf, hold on

impulse(sys,'.', time); % общий график

for ik=1:id

impulse(w(ik),'k',time); % отдельные составляющие

end

title 'Составляющие оригинала'

hold off

end

Выведем изображения отдельных составляющих и их графики

(рисунок 4.6) с помощью функции resp() за период 10 c, нумерация

изображений и графиков их оригиналов совпадает. Каждой паре ком-

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

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

ставлены функцией residue(), вывод изображений (функций передачи)

блокируется точкой с запятой.

Состав ляющие оригинала

Amplitude

0 2 4 6 8 10

-1.5

-1

-0.5

0.5

System: untitled3

Time (sec): 0.185

Amplitude: 0.607

System: untitled2

Time (sec): 1.73

Amplitude: -0.204

System: untitled1

Time (sec): 0.288

Amplitude: -0.768

System: sys

Time (sec): 2.86

Amplitude: 0.219

Рисунок 4.6

>> resp(sys,10)

Transfer function from input 1 to output:

-1.211 s - 3.004

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

s^2 + 4 s + 4.25

223

Transfer function from input 2 to output:

0.57 s - 0.508

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

s^2 + 2 s + 2

Transfer function from input 3 to output:

0.6412

-------

s + 0.3

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

ляющей, можно получить непосредственно в MATLAB, применяя об-

ратное преобразование Лапласа ilaplace(), и построить график функ-

цией ezplot(y,t), где у – исследуемая функция (рисунок 4.7).

>> syms s t

>> w=(0.57*s-0.508)/(s^2+2*s+2); t=0:0.01:10;

>> y=ilaplace(w)

y =

(57*(cos(t) - (539*sin(t))/285))/(100*exp(t))

>> ezplot(y,t)

0 2 4 6 8 10

-0.4

-0.2

0.2

0.4

0.6

(57 (cos(t) - (539 sin(t))/285))/(100 exp(t))

Рисунок 4.7

Таким образом, для каждой составляющей переходного процес-

са мы получили аналитическое выражение и возможность построения

отдельного графика.

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

математики, имеющийся в MATLAB. Символьные переменные долж-

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

взята в апострофы. Объявление производится оператором syms.

При преобразовании Лапласа выражение выводится с коэффи-

циентами, представленными в виде рациональных дробей. Построен-

ный функцией ezplot() график полностью совпадает с импульсной ха-

рактеристикой от этого же изображения. Еще один пример.

224

>> syms s t

>> q=1/(s^3+6*s^2+11*s+6); temp=ilaplace(q)

temp =

1/(2*exp(t)) - 1/exp(2*t) + 1/(2*exp(3*t))

Возвращаемый MATLAB версии 7.11 результат обратного пре-

образования Лапласа неудовлетворителен с точки зрения учебного

процесса (отрицательные полюсы представлены положительными

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

понятным для студентов был результат в версии 7.01, совпадающий с

описанием экспоненциальной функции в TF и ZPK моделях.

temp =

-exp(-2*t)+1/2*exp(-3*t)+1/2*exp(-t)

Функция dirac(t) в символической математике соответствует

единичному импульсу (дельта-функции), а heaviside(t) – единичному

скачку. Они функционируют как при прямом, так и при обратном

преобразовании Лапласа.

>> syms s t

>> q=heaviside(t)+2*exp(-3*t);

>> laplace(q)

ans =

1/s+2/(s+3)

>> w=(3*s^2+2*s+1)/(s^2+2*s+2);

>> temp=ilaplace(w) % обратное преобразование Лапласа

temp =

3*dirac(t) - (4*(cos(t) + sin(t)/4))/exp(t)

>> simplify(laplace(temp)) % прямое преобразование Лапласа

ans =

3 - (4*s + 5)/(s^2 + 2*s + 2)

Принцип построения TF модели нарушен и для использования

подобного результата в традиционном виде передаточной функции

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

>> ww=tf([4 5],[1 2 2]); sys=3-ww

Transfer function:

3 s^2 + 2 s + 1

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

s^2 + 2 s + 2

Реакцию на любое воздействие, имеющее преобразование Лап-

ласа, можно определить как произведение ПФ двух блоков, представ-

ляющих входное воздействие и систему. Например, на вход системы с

ПФ W(s) = (s+1)/(s^2+s+3) поступает сигнал 10e

-3t

, изображение кото-

225

рого равно R(s) = 10/(s+3). Строим временную характеристику, ис-

пользуя функцию impulse(), поскольку она не вносит дополнительных

полюсов или нулей в изображение реакции

>> w=tf([1,1],[1,1,3]); r=tf(10,[1,3]); rw=w*r; impulse(rw)

Вычислим теперь реакцию системы с ПФ W(s) = 1/(s+7) на про-

извольную функцию времени x(t) = 5e

-0.1t

+ 0.2t

аналитически

>> syms t s

>> xt=5*exp(-0.1*t)+0.2*t^2

xt =

5/exp(t/10) + t^2/5

>> xs=laplace(xt)

xs =

5/(s + 1/10) + 2/(5*s^3)

>> y=xs*1/(s+7)

y =

(5/(s + 1/10) + 2/(5*s^3))/(s + 7)

>> r=ilaplace(y)

r =

50/(69*exp(t/10)) - 85888/(118335*exp(7*t)) - (2*t)/245 + t^2/35 + 2/1715

Задав вектор отметок времени, можно построить график

>> t=0:0.1:30; yt=subs(r);

>> plot(t,yt),grid; title('Реакция')

Результат не очень удобен для использования, функции simplify

и numden не способствуют упрощению. Разовое применение функции

перевода в цифровой формат vpa() к результату

>> vpa(r)

ans =

0.72463768115942028985507246376812/exp(0.1*t) –

0.72580386191743778256644272615879/exp(7.0*t) – ...

выводит выражение, которое не приводится здесь целиком в связи с

нецелесообразностью полного воспроизведения – его в таком виде

сложно даже прочитать. Попробуем преобразовать выражение в

уменьшенный цифровой формат, который можно задать и отдельной

функцией digits(D), где D – число обрабатываемых разрядов. Однако,

поскольку влияние этой команды на последующие операции нежела-

тельно, лучше воспользоваться возможностью временного изменения

числа разрядов в самой функции vpa(). Двойное преобразование

функцией vpa() обеспечивает эффект уменьшения разрядной сетки

для всех чисел, входящих в преобразуемое выражение.

226

>> vpa(vpa(r,4),4)

ans =

0.7246/exp(0.1*t) - 0.7258/exp(7.0*t) - 0.008163*t + 0.02857*t^2 + 0.001166

Довольно часто выражение, возвращаемое обратным преобразо-

ванием Лапласа, неприемлемо для дальнейшего использования, на-

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

мительный оригинал

>> syms s

>> w=5/(s^3+2*s^2+3*s+4); temp=ilaplace(w)

temp =

5*sum(exp(r3*t)/(3*r3^2 + 4*r3 + 3), r3 in RootOf(s3^3 + 2*s3^2 + 3*s3 + 4, s3))

В версии 7.01 в этом случае частично помогало двойное приме-

нение функции перевода в цифровой формат vpa()

>> temp=vpa(vpa(temp),4)

temp =

1.094*exp(-1.651*t)-1.094*exp(-.1747*t)*cos(1.547*t)+1.044*exp(-.1747*t)*sin(1.547*t) +

+.2500*i*(2.087*exp(-.1747*t)*cos(1.547*t)+2.188*exp(-.1747*t)*sin(1.547*t)) +

+ .2500*i*(-2.087*exp(-.1747*t)*cos(1.547*t)-2.188*exp(-.1747*t)*sin(1.547*t))

хотя и полученное выражение усложнено, т. к. содержит семь состав-

ляющих переходного процесса, включая даже мнимые, в отличие от

трех составляющих, возвращаемых численной функцией residue() и

соответствующих оригиналу

1.0938e

-1.6506t

- 1.0938e

-0.1747t

cos(1.5469t) + 1.0436e

-0.1747t

sin(1.5469t).

>> [r,p,k]=residue(5,[1 2 3 4])

r =

1.0938

-0.5469 - 0.5218i

-0.5469 + 0.5218i

p =

-1.6506

-0.1747 + 1.5469i

-0.1747 - 1.5469i

k =

[]

В версии 7.11 результат (целиком не приводится) от преобразо-

вания не становится лучше

temp=

5.0*numeric::sum(exp(r3*t)/(3.0*r3^2 + 4.0*r3 + 3.0), r3 in

{((27^(1/2)*50^(1/2))/27 - 35/27)^(1/3) - 5/(9*((27^(1/2)*-

50^(1/2))/27 - 35/27)^(1/3)) - 2/3, 5/(18*((27^(1/2)*-

227

50^(1/2))/27 - 35/27)^(1/3)) - ... и так далее

Как прикажете это бесконечное выражение трактовать, в какой

форме записать в отчет по лабораторной работе или практическому

занятию? В результате использовать символьную математику MAT-

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

осторожно.

4.2 Моделирование частотных характеристик

Используя созданную модель системы, можно получить ее ос-

новные характеристики в частотной области – амплитудно-фазовую

характеристику (АФЧХ) или годограф Найквиста (Nyquist Diagram),

амплитудную (АЧХ) и фазовую (ФЧХ), действительную (ВЧХ) и

мнимую (МЧХ) характеристики, логарифмические характеристики

(ЛЧХ) или диаграмму Боде (Bode Diagram), реакцию на моногармони-

ческий сигнал.

Графики нетиповых функций обычно строят с помощью коман-

ды Plot, либо используя команды semilogx() и semilogy() – с логариф-

мическим масштабом только по оси x или только по оси y, loglog() – с

логарифмическим масштабом по обеим осям графика. Необходимые

подписи и оси на графике проще добавить в режиме редактирования

окна Figure через меню Insert, вид графиков корректируют после на-

жатия кнопки Edit Plot на панели инструментов фигуры.

Функция еvalfr() позволяет вычислить значение комплексного

коэффициента передачи (Frequency Response) системы W(jω) на одной

заданной частоте в виде действительной Real() и мнимой Imag() со-

ставляющих, например

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

>> fr=0+j*0.5; % перевод частоты 0.5 рад/с в комплексную

>> yf=evalfr(sys,fr) % расчет действительной и мнимой частей

yf =

0.2475 - 0.0972i

а его описание через амплитуду M(ω) и фазу P(ω) (в радианах) вектора

можно получить с помощью функций Abs() и Angle().

>> M=abs(yf)

M =

0.2659

>> Phi=angle(yf)

Phi =

-0.3743

228

Отклик для диапазона вещественных частот вычисляется функ-

цией freqresp(имя_системы, диапазон_частот). Диапазон частот можно

задавать в виде вектора fr = 0:0.1:10, в линейном fr = linspace(0,10,100)

или логарифмическом fr = logspace(-2,2,100) масштабе, причем во всех

случаях было задано 100 точек. В первых двух командах – для диапа-

зона частот от 0 до 10 рад/с, в последней команде для диапазона час-

тот от 10^-2 до 10^+2 рад/с. Если число точек не указано, по умолча-

нию функция linspace() использует 100 линейно распределенных, а

функция logspace() 50 логарифмически равномерно распределенных

отсчетов. При обращении y=freqresp(sys,fr) табулируется комплекс-

ный отклик в формате y(:, :, index)=Re ± Im×i.

Функция freqs(num,den) строит для выбираемых автоматически

200 логарифмически распределенных значений частот амплитудную

(в логарифмическом масштабе или lg

A(ω)) и фазную (в градусах)

характеристики по числителю и знаменателю передаточной функции

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

ра fr частот freqs(num,den,fr), количества N>2 частот freqs(num,den,N),

или табулирование комплексного отклика y = freqs(num, den, fr),

[y,fr]=freqs(num,den) в форме y = Re ±Im×i без вывода графика.

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

характеристики строятся выбором Plot Types-Bode в графопостроите-

ле ltiview() или командой bode(sys), bodeplot(sys) для диапазона

{wmin,wmax} или вектора fr частот bode(sys,fr), либо отдельных час-

тот bode(sys,[f1 f2 ...]). При этом ЛАЧХ (Magnitude) L(ω)=20lgA(ω)

измеряется в децибелах (dB), ЛФЧХ (Phase) откладывается в градусах

(degrees), а частота (Frequency) по оси абсцисс откладывается в рад/с и

в логарифмическом масштабе. Фазовый угол не будет выходить из

пределов +180° и -180°, если сбросить флажок в опции Unwrap phase

окна Properties-Options.

В форме с левым аргументом [mag, phi]=bode(sys,fr) или [mag,

phi, fr]=bode(sys) вместо графика табулируются значения амплитуды

mag(:,:,k) в дБ и фазы phi(:,:,k) в градусах на частотах fr(k) в рад/с. Как

у функций step(), impulse(), можно использовать устаревшие команды

bode(num,den,fr), bode(a,b,c,d,iu,fr), причем в форме с левым аргумен-

том mag и phi выводятся, как простые векторы-столбцы.

По команде damp(имя_системы) печатаются величины собст-

венных значений (Eigenvalue) или полюсов системы, соответствую-

щие им значения коэффициента демпфирования или затухания 

(Damping) и частот сопряжения. Обращением [frq, dmp, ev]=damp(sys)

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

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

229

Eigenvalue Damping Freq. (rad/s)

-1.75e-001 + 1.55e+000i 1.12e-001 1.56e+000

-1.75e-001 - 1.55e+000i 1.12e-001 1.56e+000

-1.65e+000 1.00e+000 1.65e+000

Если модель системы представлена в виде нормированных по-

линомов второго порядка

asas 

, частота сопряжения ω

может

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

, а коэф-

фициент демпфирования (затухания)  по формуле







 .

>> q=zpk(w)

Zero/pole/gain:

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

(s+1.651) (s^2 + 0.3494s + 2.423)

>> natural_freq=sqrt(2.423)

natural_freq =

1.5566

>> damping=0.3494/freq/2

damping =

0.1122

По команде dcgain(имя_системы) выводится значение коэффи-

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

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

ристик. При установке свойства Display Format ZPK-модели в состоя-

ние Frequency все перечисленные параметры (коэффициент DC Gain

0.25, частота сопряжения Natural Frequency 1.5566, коэффициент

демпфирования Damping 0.1122=0.2244/2) непосредственно наблюда-

ются в описании модели

>> set(q,'DisplayFormat','frequency'); q

Zero/pole/gain:

0.25

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

(1+s/1.651) (1 + 0.2244(s/1.557) + (s/1.557)^2)

Функция bandwidth(sys) вычисляет полосу пропускания частот-

ной характеристики SISO системы, fb = bandwidth(sys) возвращает

первую по возрастанию частоту в рад/с, на которой расхождение меж-

ду текущим и установившимся значением коэффициента усиления

достигает -3 дБ. Для FRD модели функция возвращает начальную час-

тоту приближения к коэффициенту усиления в установившемся ре-

230

жиме. Форма fb = bandwidth(sys,dbdrop) позволяет самостоятельно ус-

тановить образцовую величину расхождения dbdrop вместо стандарт-

ной -3 дБ или 70.79 % от начального значения.

Для SISO-систем штатные функции bodemag(sys) и sigma(sys),

sigmaplot() одинаково строят логарифмическую амплитудную харак-

теристику, различаясь заголовками Bode Diagram и Singular Values.

Для обеих можно задать диапазон частот wmin-wmax или частотный

вектор fr, например, sigma(sys,fr), sigma(sys,{wmin, wmax}). При об-

ращении sv=sigma(sys) или [sv,fr]=sigma(sys) выводятся в виде табли-

цы сингулярные значения σ

(дБ), служащие обобщенными пока-

зателями колебательности и робастности системы и представляющие

собой р = min(m,n) положительных квадратных корней из собствен-

ных значений матрицы W

)( WW

eig



где W – передаточная матрица размера n×m, i = 1...p, а σ

>σ

> ... >0

вычисляются в зависимости от частоты ω.

Для MIMO-систем число графиков Magnitude функции bodemag

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

sigma строит на одном графике характеристики, количество которых

равно наименьшему из числа входов и выходов. Это же число сингу-

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

для каждой частоты k и выбирается из массива командой sv(:, k).

Если применить команду plot(real(y),imag(y)) или просто plot(y),

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

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

y(jω). АФЧХ также выводится опцией Plot Types-Nyquist в контекст-

ном меню графического анализатора ltiview() или отдельной командой

nyquist(имя_системы), nyquistplot(имя). В командах nyquist(sys, fr) или

nyquist(sys, {wmin, wmax}) задается вектор частот fr или его гранич-

ные частоты wmin-wmax, вектор частот может быть задан непосредст-

венно nyquist(sys,0.01:0.1:10) или nyquist(sys,logspace(-2,2,100)). Об-

ратная АФЧХ, при необходимости, может быть построена командой

nyquist(inv(sys)). При обращении с левым аргументом табулируется

комплексный отклик отдельно [Re,Im]=nyquist(sys) или дополнитель-

но с вектором использованных частот [Re,Im,fr]=nyquist(sys) в форма-

те Re(:, :, k), Im(:, :, k), fr(k). Допустимы устаревшие формы обраще-

ния nyquist(num,den,fr) или nyquist(a,b,c,d,iu,fr), как у bode().

Рекомендуется отменить команду Negative Frequencies (отрица-

тельные частоты) в опции Show контекстного меню графика для вы-