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

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

211

4 Исследование систем управления

4.1 Моделирование временных характеристик

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

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

(Step Response), импульсную функцию (Impulse Response), реакцию на

произвольное воздействие (Linear Simulation Results), реакцию на на-

чальные условия (Response to Initial Conditions) и аналитическое раз-

ложение на простые дроби (Residue). Расчет временных характеристик

выполняется в виде таблицы и с построением графиков.

Команда ctrlpref вызывает окно Control System Toolbox Prefe-

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

единицы измерения (Units), формат подписей и осей (Style), зону

ошибки Δ для времени регулирования и способ оценки времени на-

растания (Options), параметры программы проектирования регулятора

для одномерных объектов (закладка SISO Tool).

Переходную характеристику строит при обращении без левого

аргумента функция step(sys), stepplot(sys) или графический анализатор

ltiview(имя_системы) по умолчанию Вы можете установить отличи-

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

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

фике. Например, система 1 – сплошная красная (red) линия, система 2

– штриховая зеленая: step(sys1,'r',sys2,'g--'). Опция контекстного меню

Normalize вписывает все характеристики в графическое окно незави-

симо от масштаба для их наилучшего обзора. Следует иметь в виду,

что в качестве последнего значения характеристики Final Value на са-

мом деле всегда указывается значение коэффициента усиления k

уст

Команда [y, t] = step(sys) выводит столбцами значения реакции и

времени, команда [y,t,x] = step(sys) – реакцию, время и переменные

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

дать конечным значением tf и вектором Т, например, step(sys, tf) или

step(sys, T), несколькими значениями step(sys,[t1 t2 ...]).

Устаревшие команды step(num,den,t), step(a,b,c,d,iu,t) допускают

обращение по отдельным элементам TF и SS моделей. Индекс iu в

MIMO системах указывает номер входа, для которого строятся харак-

теристики. Можно выбрать отдельную характеристику, например, от

первого входа ко второму выходу step(a,b(:,1),c(2,:),d(2,1)). Имеются

отличия для устаревших форм в порядке перечисления табулируемых

параметров [y,x,t]=step(a,b,c,d) и в форме вывода времени, как векто-

ра-строки, а не столбца.

212

Если схема собрана в Simulink, переходную характеристику

можно увидеть на экране осциллографа Scope, подключив ко входу

системы блок Constant и задав в нем необходимую величину (1 для

переходной характеристики, произвольное значение – для кривой раз-

гона). Но лучше использовать для этого блок Step, в свойствах кото-

рого можно установить не только начальную и конечную величину

сигнала (Initial-Final Value), но и момент его формирования (Step

Time). Для пуска симуляции служит кнопка Start Simulation на панели

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

цесса. При несоответствии масштаба Scope иногда кажется, что реак-

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

масштаб графика опцией AutoScale или выбором пределов в Axes

Properties контекстного меню экрана.

Использовав в контекстном меню графика LTI Viewer команду

Plot Types-Impulse, получим импульсную (весовую) характеристику

системы, обозначаемую обычно w(t) или g(t). Если вызывать эту и

другие характеристики не через графопостроитель ltiview(), а отдель-

ной функцией, например, impulse(sys), выводится сразу рисунок (Fig-

ure). Форма [y,t]=impulse(sys) выводит таблицу значений y, t, форма

[y,t,x] = impulse(sys) добавляет переменные состояния для много-

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

значения tf или в виде вектора Т, например, impulse(sys, tf) или im-

pulse(sys, T). Большими возможностями обладает форма impulseplot().

Все, сказанное выше относительно анализа нескольких систем

сразу или устаревших форм обращения для функции step(), относится

и к функции impulse(), например, особенности табулирования

[y,x,t]=impulse(num,den).

Вычисление реакции на начальные условия функцией initial()

производится в MATLAB только при описании системы в пространст-

ве состояний для вектора начальных условий x(0), который, например,

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

величины. Так, реакция многомерной системы sys на начальные усло-

вия x1(0)=1; x2(0)=0; x3(0)=0; x4(0)=0 с нанесением сетки определяет-

ся командами

>>x0=[1 0 0 0]; initial(sys,x0), grid

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

х(0) задаются внутри системы, а не на ее выходах, т. е. без учета мат-

рицы с. Хотя в MATLAB переменные состояния обычно назначаются

от входа и первой в традиционном понимании переменной у(0) здесь

будет соответствовать n-ая, т. е. последний элемент вектора х(0), при

213

балансировке матриц или построении модифицированных канониче-

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

менных. Поэтому при задании начальных значений переменных со-

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

переменным состояния MATLAB.

При формировании структурной схемы модели в Simulink для

блоков State-Space, LTI system, интеграторов или звеньев чистого за-

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

менных на выходе. Кроме того, в библиотеку Simulink Extras-

Additional Linear включены специальные блоки с заданием начальных

условий на входах и выходах (State-Space with initial outputs, Transfer

Fcn with initial outputs, Zero-Pole with initial outputs), с заданием на-

чальных значений переменных состояния (Transfer Fcn with initial

states, Zero-Pole with initial states).

Использование в контекстном меню графика LTI Viewer коман-

ды Plot Types-Initial Condition вызывает окно Linear Simulation Tool с

открытой вкладкой Initial States, где можно выставить начальное зна-

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

вектор начальных условий x(0) из Workspace кнопкой Import State

Vector. Процесс моделирования запускается кнопкой Simulate.

Наконец, в MATLAB с помощью функции lsim() возможен рас-

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

отсчетами (матрица u) в течение интервала времени (вектор t). Число

столбцов матрицы u должно быть равно числу входов системы, вектор

t формируется с равномерным шагом. Построим график реакции сис-

темы sys с единственным входом на гармоническое воздействие sin(t)

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

>> t = 0:0.01:5; u = sin(t); lsim(sys,u,t)

Реакция одновременно на начальные условия и вынуждающее

воздействие требует формы обращения lsim(sys, u, t, x0) и вычисляет-

ся только в пространстве состояний. Наличие левой части выражения

отменяет построение графика: при обращении y = lsim(sys,u,t,x0) вы-

числяется матрица y, размер которой определяется числом отсчетов и

выходов системы; при обращении [y,t,x] = lsim(sys,u,t,x0) вычисляется

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

отсчета времени. Напомним, что график может быть получен с помо-

щью функций lsimplot(), ltiview('lsim',sys,u,t,x0) и непосредственно из

LTI Viewer через контекстное меню Plot Types-Linear Simulation.

На графике отображается и входной сигнал u(t), если в контек-

стном меню графика LTI Viewer выбрана опция Show Input.

214

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

кривой между точками отсчета (один для всех систем) в виде строко-

вой переменной 'zoh' – фиксация нулевого порядка или 'foh' – фикса-

ция первого порядка (линейная интерполяция), например, командой

lsim(sys1,sys2,u,t,'zoh'). По умолчанию метод интерполяции выбира-

ется автоматически по наилучшему сглаживанию кривой.

Сравним на графике изменение во времени переменных состоя-

ния х

… х

и выхода y системы sys, описанной SS-моделью (рисунок

4.1). Обозначения y, x1, x2, x3 нанесены на готовый график.

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

>> t=0:0.1:30; u=t.^0; % воздействие в виде единичного скачка

>> [y,t,x] = lsim(sys,u,t); % нулевые начальные условия

>> plot(t,y,t,x),grid; title('Переходный процесс');

>> xlabel('Время, с')

Рисунок 4.1

В соответствии с нумерацией переменных программа строит

первый график синим цветом, второй – зеленым, третий – красным и

т. д. Форма и численные значения графика y(t) будут соотноситься с

графиками переменных состояния так, как это определяет матрица с.

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

и амплитуда реакции в два раза выше согласно коэффициенту с(1) = 2.

Очевидно, что при с = [1 1 0] будет выполняться суммирование двух

составляющих y(t) = x

(t) + x

(t).

Если дополнительные параметры не указаны при обращении из

командной строки, открывается пустое окно Linear Simulation Results

и окно ввода данных Linear Simulation Tool (рисунок 4.2).

На вкладке Input Signals в области Timing зададим время окон-

чания процесса End Time 25 с и интервал отсчета 0.1 с (время начала

215

установлено в нуль). На вкладке System Inputs число каналов соответ-

ствует числу входов. Инструмент Linear Simulation Tool позволяет

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

лый шум) или как объект, заданный в Workspace или MS Excel, в виде

mat, ASCII или текстового файла. Выберем проектирование сигнала

Design Signal, откроется окно Signal Designer с четырьмя типовыми

функциями Sine Wave, Square Wave, Step Function, White Noise, оста-

новимся на периодическом сигнале прямоугольной формы Square

Wave. Задаваемые параметры: частота повторения Frequency (Hz) 0.1,

амплитуда 1, длительность моделирования Duration (sec) 25 c совпада-

ет с заданной раньше. Ввод данных осуществляется кнопкой Insert,

они отображаются на вкладке System Inputs. Количество отсчетов

равно 251.

Рисунок 4.2

Запустив моделирование кнопкой Simulate, получим для систе-

мы с ПФ W(s) = 10/(s

+ 2s

+ 3s + 4) график (рисунок 4.3). При необ-

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

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

ют вкладку Initial States.

216

Linear Simulation Results

Time (sec)

Amplitude

0 5 10 15 20 25

-5

Рисунок 4.3

Вместо автоматического выбора способа интерполяции по наи-

лучшему виду кривой (Automatic по умолчанию) может быть выбран

вручную метод либо с фиксацией нулевого порядка Zero Order Hold,

либо с линейной интерполяцией First Order Hold. Импортировать го-

товые массивы для времени и сигнала позволяют кнопки Import time,

Import Signal (с выбором вектора-строки или вектора-столбца). Ввод

данных возможен для нескольких каналов (графиков) сразу.

Сформировать в рабочем пространстве массивы временных t и

сигнальных u отсчетов позволяет функция [u, t]=gensig(Type, Tau), где

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

ная Type может принимать значения 'sin', 'square', 'pulse' соответствен-

но для сигналов типа синусоида, меандр, импульс. Расширенное об-

ращение [u, t]=gensig(Type, Tau, Tf, Ts) дает возможность задать об-

щую длительность сигнала Tf и шаг отсчета Ts.

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

нала sin(), cos(), качающейся частоты chirp(), прямоугольного square(),

линейно изменяющегося (треугольника или пилы) sawtooth(), обычно

требующие указания периода.

Вычисление свободной составляющей вызывает трудности при

традиционном задании начальных значений рассогласования у(0),

скорости y'(0) и ускорения y''(0), т. е. показателей, фиксируемых на

выходе системы. Для использования функции lsim() в подобном слу-

чае следует преобразовать вектор y(0) в вектор х(0), разделив его на

элемент матрицы с, соответствующий переменной состояния, назна-

ченной на выходе системы (остальные элементы матрицы с должны

быть нулевыми). Пусть заданы начальные условия y(0)=-1.5, y'(0)=75.

С учетом особенностей MATLAB формируем вектор y(0) в обратном

порядке и преобразуем его в вектор х(0), кривая реакции системы дей-

ствительно начинается со значения y(0) = -1.5.

217

>> num=5; den=[0.048 0.4 10]; % задаем ПФ системы

>> y0=[75 -1.5]; % вектор начальных условий

>> [a,b,c,d]=tf2ss(num,den); % переход к канонической форме

>> sys=ss(a,b,c,d); % модель в пространстве состояний

>> x0=y0/c(2)

x0 =

0.7200 -0.0144

>> t=0:0.01:2; u=t.^0; % воздействие в виде единичного скачка

>> lsim(sys,u,t,x0)

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

ставлении (с единственной единицей в матрице с), тогда преобразова-

ния начальных значений не требуется. Чтобы избежать указанных за-

труднений, разработаны функции dtoy() и xsim(), предназначенные

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

процесса по передаточной функции одномерной системы. Функция

dtoy(den,x0) по характеристическому полиному и начальным значени-

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

ной составляющей, либо полином N

(s) числителя ее изображения.

function[num0,y0,t0]=dtoy(den,x0)

% Расчет свободной составляющей по знаменателю ПФ

[m,n]=size(x0);

if m>n % перевести вектор-столбец в вектор-строку

x0=x0';n=m;

end

for i=1:n % расчет числителя N

(s) изображения Y

св

(s)

num(i)=sum(x0(1:i).*fliplr(den(1:i)));

end

w=tf(num,den); % формирование изображения Y

св

(s)

if nargout==0 % если только построить график

impulse(w), grid;

title('Свободная составляющая');

ylabel('Yсв(t)'); xlabel('Время')

else

num0=num; [y0,t0]=impulse(w); % иначе выводим таблицу

end

Результат использования созданной функции (рисунок 4.4).

>> den=[0.048 0.4 10]; x0=[-1.5 75]; [n0]=dtoy(den,x0)

n0 = % вывод полинома N

(s) числителя изображения Y

св

(s)

-0.0720 3.0000

>> [n0,y0,t0]=dtoy(den,x0) % вывод N

(s) и таблицы отклика

>> dtoy(den,x0) % вывод только графика

218

Рисунок 4.4

Для расчета произвольного переходного процесса, включа-

ющего вынужденную и свободную составляющие, одномерной систе-

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

xsim(), являющаяся аналогом штатной функции lsim().

function[y,t]=xsim(w,u,t,x0)

% Расчет произвольной реакции системы

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

[m,n]=size(x0);

if m>n % перевести вектор-столбец в вектор-строку

x0=x0';n=m;

end

for i=1:n % расчет числителя изображения Yсв(s)

n0(i)=sum(x0(1:i).*fliplr(den(1:i)));

end

w0=tf(n0,den);

y0=impulse(w0,t);

yt=lsim(w,u,t);

y=y0+yt;

if nargout==0 % если нет левых (выходных) аргументов

plot(t,y,'linewidth',2); hold on;

plot(t,yt,'--k',t,y0,':k',t,u,'k');

title('График выходной величины'); xlabel('Время, с');

legend('Y(t)','Yв(t)','Yo(t)','u');

hold off

end

Таблица значений выходной величины y(t), полученных на за-

данном массиве t, выводится при обращении y=xsim(w,u,t,x0), в этом

219

случае график не строится. Пример расчета переходного процесса для

сложного вынуждающего воздействия и заданных на выходе системы

начальных условий с представлением графиков (рисунок 4.5).

>> w=tf(5,[0.048 0.4 10]); % система

>> x0=[-1.5 75]; % начальные условия

>> t=0:0.001:0.5; % интервал времени

>> u=5+100*t.*(t<0.1)+10*(t>=0.1); % управляющий сигнал

>> xsim(w,u,t,x0)

Рисунок 4.5

Реакцию системы в аналитическом виде (оригинал) по ее изо-

бражению Лапласа можно найти с помощью функции residue(), задав

полиномы b(m), a(n) числителя и знаменателя рациональной дроби

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

residue() выводит значения вычетов в полюсах (коэффициентов чис-

лителей), полюсов и свободной части разложения на простые дроби.

Вычеты (элементы вектора r) печатаются в том же порядке, в каком

выводятся полюса (вектор p) – первый сверху вычет соответствует

первому полюсу, т. е. корню характеристического уравнения, и т. п.

Пользуясь таблицей соответствий оригиналов и их изображений по

Лапласу, можно составить оригинал для действительного полюса λ с

вычетом α

)exp( t









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





– действи-

220

тельная и мнимая части вычета,





– аналогичные части полюса)

   

.)sin()cos(exp ttt2

















Плюс перед  справа записывается при несовпадении знака мни-

мой части у соответствующих друг другу вычета и полюса.

Коэффициенты для кратных полюсов выводятся в порядке воз-

растания кратности, т. е. сначала коэффициент и полюс в степени 1,

затем в степени 2 и так до максимальной степени корня. В качестве

примера найдем коэффициенты разложения для функции 5/(s+1)^3.

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

r =

p =

-1.0000

k =

[]

Напомним, что оригинал подобного разложения должен быть

записан по таблице соответствия в виде 0∙e

-t

+ 0∙t∙e

-t

+ 5/2∙t

∙e

-t

Обращение [num, den]=residue(r, p, k) позволяет восстановить по

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

правильность вычисления коэффициентов числителей простых дро-

бей, если они выполнялись не в среде MATLAB. Отсутствующее зна-

чение k должно вводиться при этом в виде пустой матрицы. Заметим,

что элемент вектора k будет ненулевым только при равенстве степе-

ней полиномов числителя и знаменателя изображения.

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

строить временную характеристику, например, организовав цикл типа

FOR … END в интервале времени от 0 до t

уст

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

ва точек (время t

уст

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

из столбца р по формуле 10/min|Re(p)|). Длинную формулу можно

продолжить на другой строке, прервав знаком троеточия. Внимание!

Вместо матричного умножения используется поэлементное перемно-

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

Найдем оригинал передаточной функции с помощью разложе-

ния на простые дроби и обратного преобразования Лапласа. Коэффи-

циенты числителей и корни знаменателей разложения на простые