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

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

201

>> sdelay([],[],[],'term') % закончить компиляцию

Следует привести к единому масштабу все временные парамет-

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

передачи), они задаются, как правило, в секундах. Не стоит забывать,

что при линеаризации звено чистого запаздывания (Transport Delay)

моделируется в Simulink приближенно разложением в ряд Паде. При

нулевом порядке аппроксимации Паде, который выставлен по умол-

чанию, задержка проявляется на экране приборов (Scope), но у линеа-

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

нице и задержка отсутствует, независимо от того, какая величина ее

была установлена.

Используя разложение Паде, MATLAB воспроизводит задержку

в виде провала (Undershoot) или нескольких колебаний в начале пере-

ходного процесса. Во избежание этого можно использовать график с

экрана Scope, открыв в нем кнопкой Parameters окно с панелью Data

history, установив флажок в опции Save data to workspace и формат

Array вывода данных в переменную ScopeData. Такая возможность

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

графика, построенного командой plot(ScopeData(:,1), ScopeData(:,2)),

grid, можно отнести неравномерный шаг по времени и отсутствие ав-

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

командой S=stepinfo(ScopeData(:,2), ScopeData(:,1)). Равномерный шаг

квантования по времени Fixed-step задается на панели Solver меню

Simulation-Configuration Parameters (по умолчанию стоит переменный

Variable-step). В режиме auto шаг выравнивается с параметром Sample

time в блоках To Workspace.

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

пространство в виде массивов tout, yout используются элементы Out,

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

данных. Построить график можно командой plot(tout, yout), получить

показатели качества командой S=stepinfo(yout, tout), при каждом за-

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

С помощью функции linearize() может выполняться линеариза-

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

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

мости подключать к схеме другие источники или приемники, либо

размыкать схему в реальности.

Рассмотрим систему с именем mdl (рисунок 3.13). Блоки Display

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

данном случае коэффициентов a1, a2, a3, блоки Constant и Scope для

202

линеаризации необязательны. У сумматоров отображаем имя коман-

дами контекстного меню Format-Show Name и Format-Flip Name.

Рисунок 3.13

Задаем значения переменных и вызываем модель.

>> clc; clear; a1=1; a2=2; a3=3; % задаем параметры

>> open mdl % открыть S-модель, сделать необходимые установки

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

чалом на выходе сумматора Sum и окончанием на выходе блока w4,

задавая точки входа и выхода io, точку рабочего режима OpPoint.

Функция io=linio('blockname',portnum,type,openloop) формирует

контрольные точки системы с указанием имени блока 'blockname', на

выходе которого задается точка, номера порта portnum, типа type точ-

ки – вход In или выход Out, статуса цепи в контрольной точке

openloop – размыкать On или нет (по умолчанию) систему. Получение

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

дится функциями getlinio() и setlinio(). Места размещения контроль-

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

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

весного режима, которую отыскивает функция findop(). Первый вари-

ант обращения opt=findop('model',ops) требует предварительного оп-

ределения спецификации функцией ops=operspec('model'), второй ва-

риант opt=findop('model', tsnapshot) подразумевает указание момента в

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

равновесное состояние Trim Condition (проще всего взять конечное

время моделирования). При нескольких точках равновесия у системы

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

вой линеаризированной модели. Даже если команда завершена точкой

с запятой, выводится сообщение Operating Point Search Report.

>> io(1)=linio('mdl/Sum',1,'in'); % задаем точку входа

>> io(2)=linio('mdl/w4',1,'out','on'); % задаем точку выхода

203

>> setlinio('mdl',io); % ставим контрольные точки

>> os=operspec('mdl'); % спецификация рабочей точки

>> OpPoint=findop('mdl',os); % расчет по первому варианту

Operating Point Search Report:

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

Operating Report for the Model mdl.

(Time-Varying Components Evaluated at time t=0)

Operating point specifications were successfully met.

States:

----------

(1.) mdl/w1

x: 0.875 dx: 2.22e-016 (0)

(2.) mdl/w2

x: 0.375 dx: -2.22e-016 (0)

(3.) mdl/w3

x: 0.125 dx: 0 (0)

Inputs: None

----------

Outputs: None

----------

>> sys=linearize('mdl',OpPoint,io);

>> w_openloop=tf(sys)

Transfer function from input "Sum" to output "w4":

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

s^3 + 6 s^2 + 12 s + 7

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

системе, это можно проверить вручную. Заметим, что для правильно-

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

блок w4 с единичным коэффициентом усиления, чтобы отделиться от

точки x1, иначе программа размыкала в контрольной точке 2 и внут-

ренний контур обратной связи, начинающийся также в точке х1. Блок

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

При втором варианте применения функции findop() сообщение

Operating Point Search Report не выводится.

>> clear; a1=1; a2=2; a3=3;

>> open mdl

>> op=findop('mdl',20);

>> io(1)=linio('mdl/Sum',1,'in');

>> io(2)=linio('mdl/w4',1,'out','on');

>> setlinio('mdl',io);

>> sys=linearize('mdl',op,io);

204

>> w=tf(sys)

Transfer function from input "Sum" to output "w4":

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

s^3 + 6 s^2 + 12 s + 7

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

цию idss, которая описывает модель в специальной форме простран-

ства состояний, где e(t) – шумы, К – матрица коэффициентов при них.

x(t+Ts) = A x(t) + B u(t) + K e(t); x(0)=x0

y(t) = C x(t) + D u(t) + e(t)

Для непрерывных систем Ts = 0. Все свойства можно получить

командой get(idss), а их описание idprops idss, сюда относятся:

A, B, C, D, K, X0 – матрицы пространства состояний,

SSParameterization – индикатор метода параметризации матриц:

'free' (по умолчанию) – все параметры матриц A, B, C выбира-

ются свободно, вид D, K, X0 зависит от свойств 'nk', 'Distur-

banceModel' и 'InitialState'.

'canonical' или 'сan' – A, B и C приводятся к наблюдаемой кано-

нической форме в зависимости от параметра 'CanonicalIndices' (обыч-

но auto), вид D, K, и X0 зависит от 'nk', 'DisturbanceModel', 'InitialState'.

'structured' – вид описания определяется структурой матриц As,

Bs, Cs, Ds, Ks и X0s.

Структура этих матриц может задаваться командами

SET(model,'As',An,'Bs',Bn,'Cs',Cn,'Ds',Dn,'Ks',Kn,'X0s',X0n) или укоро-

ченной SETSTRUC(model,An,Bn,Cn,Dn,Kn,X0n), для idss-модели оп-

ределяется только в случае 'structured'. Значение NaN у какого-либо

элемента этих матриц означает его настройку по результатам оцени-

вания, например, As = [1 NaN; 0 NaN], A = [1 2; 0 3] означает, что эле-

менты a

и a

будут подбираться при номинальных/начальных зна-

чениях 2 и 3.

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

структурной схемы на выходе (сигнал е) и на инвертирующем входе

(сигнал yos) главного сумматора Sum включены блоки To Workspace

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

которой этот блок подключен (рисунок 3.14), в обоих установлены

шаг квантования Sample time 0.1 с и формат сохранения данных Save

format типа массив Array. На вкладке Simulation-Configuration Parame-

ters-Solver заданы время моделирования Stop time 20 c, шаг постоян-

ный Fixed-step величиной step size 0.1 c (одинаковый с блоками To

205

Workspace) и метод решения дифференциального уравнения ode4

(Runge-Kutta). На этот раз источник сигнала обязателен, в качестве его

использован блок Step с указанием нулевого времени скачка Step time.

Рисунок 3.14

Создаем структурную заготовку idm, наследующую все сведе-

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

сток, заполняя матрицы As, Bs,Cs, Ds значениями NaN (для после-

дующей оценки), и оставляя матрицы Ks и X0s нулевыми (шумы от-

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

командой sim, формируем командой iddata(выход, вход, шаг) специ-

альный массив данных da для идентификации c входным е, выходным

yos сигналами и шагом Sample time, установленным в блоках To

Workspace и на вкладке Solver.

>> idm=idss(sys); % структурный каркас модели

>> idm.As=ones(size(sys.A))*NaN; % все параметры оцениваются

>> idm.Bs=ones(size(sys.B))*NaN;

>> idm.Cs=ones(size(sys.C))*NaN;

>> idm.Ds=ones(size(sys.D))*NaN;

>> idm.Ks=zeros(size(sys.B)); % шумы отсутствуют

>> idm.x0s=zeros(size(sys.B)); % начальные условия нулевые

>> sim mdl % моделируем

>> da=iddata(yos,e,0.1); % формируем данные для идентификации

>> m=pem(da, idm); % идентификация моделью с заданной структурой

>> m1=tf(m);

>> wi=m1(1) % выводим главную передаточную функцию

Transfer function from input "u1" to output "y1":

-2.636e-005 s^2 + 1.419e-005 s + 0.9999

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

s^3 + 5.998 s^2 + 12.05 s + 6.999

Input groups:

Name Channels

Measured 1

206

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

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

построение переходной характеристики обеих моделей (рисунок 3.15).

Характеристика системы wi отображается круглыми маркерами

через 10 отсчетов.

>> t=tout(1:10:length(tout)); step(w,’k’,tout); hold on;

>> step(wi,'k.',t); legend('w','wi','location','best');

>> hold off

0 5 10 15 20

0.05

0.1

0.15

From: Sum To: w 4

Step Response

Amplitude

w i

Рисунок 3.15

Хотя с помощью выходных данных e, yos блоков To Workspace

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

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

том состоянии, о чем не следует забывать (рисунок 3.16). У разомкну-

той системы входной сигнал e(t) имел бы вид единичного скачка.

>> plot(tout,e,'--',tout,yos), title 'Динамика системы'

>> legend('e','yos','location','best'), xlabel 'Время, с'

Рисунок 3.16

207

Если параметр Decimation в блоке To Workspace отличался от

единицы, например, замер производился через n шагов, то для постро-

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

менить шаг временного массива. Из-за интерполяции разреженность

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

>> t=tout(1:n:length(tout));

>> plot(t,e,'.',tout,yos), title 'Динамика системы'

На основе линеаризации фрагментов схемы с указанием кон-

трольных точек функционируют блоки-графоанализаторы, включен-

ные в библиотеку Simulink Control Design-Linear Analysis Plots: Bode

Plot, Linear Step Response Plot, Pole-Zero Plot, Gain and Phase Margin

Plot (Боде, Николс, Найквист, таблица), Nichols Plot, Singular Value

Plot. Все они позволяют установить на графике проектные ограниче-

ния по соответствующим параметрам и выдать сообщение при их вы-

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

лизатора к схеме (рисунок 3.17) рассмотрим на примере одного блока

Linear Step Response Plot, остальные применяются аналогично.

Рисунок 3.17

При вызове свойств блока открывается окно Block Parameters с

панелями линеаризации Linearizations, установки ограничений

Bounds, сохранения результатов Logging и сигнализации выхода за

рамки ограничений Assertion. Нажатием кнопки [+] справа от таблицы

Linearization inputs/outputs открывается дополнительное окошко Model

signal, после чего нужно щелкнуть мышью в схеме линию, на которой

устанавливается вход в схему, и кнопкой [«] передать сведения в таб-

лицу, затем аналогично задать точку выхода и в колонке Configuration

выбрать Output (рисунок 3.18). При необходимости размыкаем цепь в

соответствующей точке, установив флажок в опции Open Loop. Кноп-

ка [×] удаляет ошибочно выбранную точку. После изменений обяза-

тельно нужно щелкнуть OK или Apply. Строим график ошибки по

скачку задания, вход – сигнал r(t), выход – сигнал e(t).

208

Рисунок 3.18

На вкладке Logging можно выбрать сохранение данных в рабо-

чее пространство Save data to workspace и имя этой структуры (по

умолчанию sys), которая включает SS модель структурной схемы.

>> sys % структура сохраненных данных

sys =

time: 0

values: [1x1 ss]

blockName: 'linear_xy/Linear Step Response Plot'

>> sys.values % описание SS модели

a =

linear_xy/In linear_xy/In

linear_xy/In 0 1

linear_xy/In -2 -1

b =

Step

linear_xy/In 0

linear_xy/In 1

c =

linear_xy/In linear_xy/In

linear_xy/Su -2 -1

d =

Step

linear_xy/Su 1

Continuous-time model.

>> [a,b,c,d]=ssdata(sys.values) % извлечение матриц

a =

0 1

-2 -1

b =

c =

-2 -1

d =

209

Кнопка Show plot выводит график после моделирования, опция

Show plot on block open обеспечивает вызов окна с графиком Linear

Step Response Plot (рисунок 3.19) вместо окна Block Parameters при

щелчке на изображении блока (в окне графика есть кнопка вызова ок-

на параметров Open Block Parameters Dialog).

Рисунок 3.19

У блоков имеются особенности: в годографе Найквиста отсутст-

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

фике step (рисунок 3.19) не вывелась табличка для Peak Response, на

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

Метод фазовой плоскости используется для исследования сис-

тем первого и второго порядка. Плоскость с двумя координатами x

, характеризующими поведение системы, и траектория движения

изображающей точки на ней называются фазовыми. Совокупность фа-

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

вым портретом. При описании системы в пространстве состояний ме-

тодом фазовых переменных за координату x

обычно принимают вы-

ход y(t), за координату x

– скорость его изменения или dx

/dt.

Значение выходной переменной в i-й момент времени определя-

ется по абсциссе изображающей точки. Развернув фазовую траекто-

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

реходном процессе траектория сходится к особой точке, при неустой-

чивом – удаляется от нее, автоколебаниям соответствует граница ус-

тойчивости. Амплитуда автоколебаний равна расстоянию по оси абс-

210

цисс от центра фазовой траектории до ее границы, а частота автоколе-

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

отклонения по оси ординат к амплитуде автоколебаний.

Исследуем методом фазовой плоскости поведение системы (см.

рисунок 3.17). На первый (верхний) вход прибора XY Graph из биб-

лиотеки Sinks подают сигнал x

, откладываемый по оси абсцисс и за-

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

Время симуляции (окно Simulation Stop Time на панели инструментов)

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

Переходная характеристика начинается в нуле (рисунок 3.20),

заканчивается при значении 0.5, процесс сходящийся.

а б

Рисунок 3.20

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

ных начальных условиях, у устойчивых автономных систем фазовая

траектория сходится к нулю. Зафиксировав значения x1, x2 с помо-

щью блоков To Workspace, получим наглядный график (рисунок 3.21).

>> plot(x1,x2),grid, title 'Фазовая траектория';

>> xlabel 'x1', ylabel 'x2'

0 0.2 0.4 0.6 0.8

-0.2

0.2

0.4

0.6

Фазовая траектория

Рисунок 3.21