Можно получить более подробный график импульсной характеристики, если выполнить функцию cra
с одним аргументом zdan (Рис. 2. 4):
Необходимо отметить, что на графиках по оси абсцисс откладываются промежутки времени τ
= t
i
– t
i
-1
, а по оси ординат значения корреляционных функций для входного u
2
и выходного у2
сигналов; значения взаимокорреляционой функции и импульсной характеристики. Из полученных характеристик следует, что с увеличением τ
наблюдается резкий спад корреляционной зависимости входного сигнала, что свидетельствует о слабой взаимосвязи между сечениями процесса, соответствующими произвольным моментам времени (процесс более близок к белому шуму, а автокорреляционная функция к дельта-функции). Выходная величина наоборот более плавно изменяет свои состояния от одного момента времени к другому и, следовательно, взаимосвязь между предыдущим и последующим значениями выходного сигнала более тесная, чем у входного.
Для получения частотных характеристик экспериментальных данных воспользуемся функциями оценивания частотных характеристик. Функция spa
возвращает частотные характеристики одномерного объекта
и оценки спектральной плотности его сигналов для обобщенной линейной модели объекта:
· М – ширина временного окна (по умолчанию М = min(30, length(z) /10), где length(z) – число строк матрицы z);
· w – вектор частот для расчета частотных характеристик (по умолчанию [1: 128]/128*pi/T);
· maxsize – параметр, определяющий максимальный размер матриц, создаваемых в процессе вычислений.
· z_spe – матрица спектральных плотностей входного и выходного сигналов.
Построим диаграмму Боде (АЧХ, ФЧХ), используя функции spa и bodeplot и данные, полученные при изучении теплового объекта, и содержащиеся в файле dryer2
Результаты моделирования (АЧХ построена в логарифмическом масштабе) без доверительного коридора представлены на рис. 2. 5.
Полученные зависимости подтверждают высокочастотную составляющую значений входного и выходного сигналов. Границы изменения частот на графиках установлены по умолчанию.
Для получения частотных характеристик вместе с доверительным коридором шириной в три среднеквадратических отклонения в пакете System Identification Toolbox MATLAB имеется следующие возможности:
· вычисление g – оценки АФХ и ФЧХ в частотном формате и phiv – оценки спектральной плотности шума с помощью команды
где 'sd' – указывает на сплошную линию доверительного коридора (по умолчанию эта линия штриховая); 3 – величина доверительного коридора в три среднеквадратических отклонения; 'fill' – способ заливки доверительного коридора (желтым цветом).
Построим АЧХ, ФЧХ, используя функции spa, bodeplot, logspace и данные, полученные в файле Project24 с соответствующим доверительным коридором:
Результаты моделирования представлены на рис. 2. 6.
Результаты моделирования представлены на рис. 2. 7.
Полученный график оценки спектральной плотности шума с доверительным коридором показывает наличие равномерного распределения мощности сигнала по частотному спектру с последующим спадом мощности на частоте выше 1, 1 рад/с.
Далее необходимо выполнить параметрическое оценивание ТОУ.
2
. 8. Параметрическое оценивание данных
Параметрическое оценивание экспериментальных данных проводится с целью определения параметров модели заданной структуры путем минимизации выбранного критерия качества модели (чаще всего – среднего квадрата рассогласования выходов объекта и его постулируемой модели).
Для проведения параметрического оценивания
массив экспериментальных данных необходимо разделить условно на две части (не обязательно равные)
>> zdanv=zdan(1:500);
>> zdane=zdan(501:1000).
Первая часть массива данных будет использоваться для параметрического оценивания и построения модели ТОУ. Вторая часть необходима будет для верификации (проверки качества) модели, определения адекватности полученной модели и определения погрешностей идентификации. Необходимо отметить, что параметрическая идентификация в пакете S
ystem Identification Toolbox MATLAB
выполняется в дискретном виде и полученные модели, являются дискретными.
В пакете System Identification Toolbox рассмотрены различные виды моделей ТОУ, которые с различной степенью достоверности описывают объект автоматизации. Для выбора наиболее приемлемой структуры и вида моделей при параметрическом оценивании экспериментальных данных в пакете System Identification Toolbox MATLAB имеются специальные функции
· параметрического оценивания,
· задания структуры модели,
· изменения и уточнения структуры модели и выбора структуры модели.
Функция оценивания ar
оценивает параметры модели авторегресии:
A
(
z
)
y
(
t
) =
e
(
t
)
,
где: A
(
z
) =
1 + a
1
z
– 1
+ a
2
z
– 2
+...+ a
na
z
–
na
, т.е. коэффициенты полинома A
(
z
)
, при моделировании скалярных временных последовательностей. Функция имеет синтаксис:
th = ar(y,n)
Или другое написание, позволяющее изменять параметры моделирования:
[th,refl]=ar(y,n,approach,win,maxsize,T)
где аргументы:
y – вектор-столбец данных, содержащий N элементов;
n – порядок модели (число оцениваемых коэффициентов);
approach – аргумент (строковая переменная) определяет метод оценивания:
• 'ls' – метод наименьших квадратов;
• 'yw' – метод Юла-Уокера;
• 'burg' – метод Бэрга (комбинация метода наименьших квадратов с минимизацией гармонического среднего);
• 'gl' – метод с использованием гармонического среднего;
Если любое из данных значений заканчивается нулем (например, burg0), то вычисление сопровождается оцениванием корреляционных функций.
Аргумент win (строковая переменная) используется в случае отсутствия части данных:
• win =‘now’ – используются только имеющиеся данные (используется по умолчанию, за исключением случая approach = ‘yw’);
• window = 'prw’ – отсутствующие начальные данные заменяются нулями, так что суммирование начинается с нулевого момента времени;
• window = 'pow’ – последующие отсутствующие данные заменяются нулями, так что суммирование расширяется до момента времени N + n;
• window = ‘ppw’ – и начальные, и последующие отсутствующие данные заменяются нулями (используется в алгоритме Юла-Уокера);
Арумент maxsize определяет максимальную размерность задачи; Т – интервал дискретизации.
Возвращаемые величины:
• th – полученная модель авторегрессии в тета-формате (внутреннем матричном формате представления параметрических моделей пакета System Identification);
• relf – информация о коэффициентах и функции потерь;
Для использования функция параметрического оценивания ar
необходимо из массива экспериментальных данных, записанных в файле dan
выделить выходную переменную у
с помощью команды
>> y=dan.y,
что равносильно команде
>> y=get(dan,'y')
>> th =ar(y,4)
Discrete-time IDPOLY model: A(q)y(t) = e(t)
A(q) = 1 - 1.348 q^-1 + 0.6695 q^-2 - 0.2531 q^-3 - 0.04431 q^-4
Estimated using AR ('fb'/'now') from data set y
Loss function 0.0140559 and FPE 0.0141588
Sampling interval: 1
Полная информация о модели авторегрессии th может быть получена с помощью команды:
>> present(th)
Discrete-time IDPOLY model: A(q)y(t) = e(t)
A(q) = 1 - 1.348 (+-0.03022) q^-1 + 0.6695 (+-0.05018) q^-2
- 0.2531 (+-0.05018) q^-3 - 0.04431 (+-0.03023) q^-4
Estimated using AR ('fb'/'now') from data set y
Loss function 0.0140559 and FPE 0.0141588
Sampling interval: 1
Created: 17-Dec-2009 10:51:00
Last modified: 17-Dec-2009 10:51:00
В информации приведены сведения о том, что модель является дискретной и для оценивания ее параметров используется прямой-обратный метод (разновидность метода наименьших квадратов), на что указывает строковая переменная 'fb' (используется по умолчанию); для построения модели используются только имеющиеся данные у
, на что указывает строковая переменная 'now' (используется по умолчанию); определены: функция потерь Loss function, как остаточная сумма квадратов ошибки и так называемый теоретический информационный критерий Акейке (Akaike's Information Theoretic Criterion – AIC) FPE; интервал дискретизации Sampling interval.
Следующая функция arx
оценивает параметры модели AR и ARX: параметры модели ARX, представленной зависимостью:
A
(z
)y(t
) = B
(z
) u
(t
) + e
(t
) ,
или в развернутом виде:
y(t
) + а1
y(t-1
) + …+ аna
y(t-n
) = b1
u(t) + b2
u(t - 1) + …+ bnb
u(t - m) + e(t).
Здесь и ниже e
(t
) – дискретный белый шум.
B(z) = b1
+ b2
z-1
+ …+ bbn
z-nb + 1
Функция имеет следующий синтаксис:
dar = arx(z,nn).
Или другое написание, позволяющее изменять параметры моделирования:
dar = arx(z,nn,maxsize,T),
где z – экспериментальные данные;
nn – задаваемые параметры модели (аргумент nn содержит три параметра: na – порядок ( число коэффициентов) полинома A
(
z
)
; nb – порядок полинома B
(
z
)
; nk – величина задержки;
maxsize - максимальная размерность задачи;
Т – интервал дискретизации.
Естественно задаться вопросом: Какую степень полинома выбрать? Известно, что с увеличением порядка полиномов улучшается степень адекватности модели реальному объекту. Однако при этом получаются громоздкие выражения, и увеличивается время моделирования. Поэтому для нахождения оптимального порядка полиномов можно воспользоваться функциями выбора структуры модели:
Функция arxstruc
вычисляет функции потерь для ряда различных конкурирующих ARX моделей с одним выходом:
v = arxstruc(ze,zv,NN),
или v = arxstruc(ze,zv,NN, maxsize);
где: ze,zv – соответственно, матрицы экспериментальных данных для оценивания и верификации моделей;
NN – матрица задания конкурирующих структур со строками вида nn = [na nb nk];
maxsize - максимальная размерность задачи.
Возвращаемая величина v – матрица, верхние элементы каждого столбца которой (кроме последнего) являются значениями функции потерь для ARX моделей, структура которых отображается последующими элементами столбцов (т. е. каждый столбец соответствует одной модели). Первый элемент последнего столбца – это число значений экспериментальных данных для верификации моделей.
Функция selstruc
осуществляет выбор наилучшей структуры модели из ряда возможных вариантов
[nn,vmod]=selstruc(v)
[nn,vmod]=selstruc(v,с),
где: v – матрица, возвращаемая функцией arxstruc;
с – строковая переменная, определяющая вывод графика или критерий отбора наилучшей структуры:
при с = ‘plot’ выводится график зависимости функции потерь от числа оцениваемых коэффициентов модели
при с = ‘log’ выводится график логарифма функции потерь;
при с = ‘aic’ график не выводится, но возвращается структура, минимизирующая теоретический информационный критерий Акейке (AIC).
при с = ‘mdl’ возвращается структура, обеспечивающая минимум критерия Риссанена минимальной длины описания;
при с
равном некоторому численному значению а,
выбирается структура, которая минимизирует значение функции потерь vmod
= v(1 + a(d/N)), где N – объем выборки экспериментальных данных, используемых для оценивания; d – число оцениваемых коэффициентов модели; v – значение функции потерь.
Возвращаемые величины: nn – выбранная структура; vmod – значение соответствующего критерия.
Например, для данных dryer2 можно задать пределы изменения порядка модели:
>> NN=struc(1:10,1:10,1);
Вычислить функции потерь:
>> v=arxstruc(zdane,zdanv,NN);
И выбрать наилучшую структуру порядков полиномов:
>> [nn,vmod]=selstruc(v,'plot'),
где 'plot' – строковая переменная, определяющая вывод графика зависимости функции потерь от числа оцениваемых коэффициентов модели (рис. 2. 8).
Рис. 2. 8. Окно выбора структуры модели
В появившемся окне столбики указывают на величину функции потерь. При подведении курсора к соответствующему столбику, в правом поле окна отразятся значения порядков полиномов na, nb, nk. В поле графика появятся рекомендации по выбору цвета столбика. Воспользуемся рекомендацией, указанной в поле графика и выберем столбик, окрашенный
красным цветом для оптимального значения порядков полиномов и нажмем кнопку Select.
Взамен строковой переменной 'plot' возможны варианты:
• 'log' – выводится график логарифма функции потерь;
• 'aic' – график не выводится, но возвращается структура, минимизирующая так называемый теоретический информационный критерий Акейке (Akaike's Information Theoretic Criterion – AIC) FPE:
,
где v
– значение функции потерь, d
– число оцениваемых коэффициентов модели, N
– объем экспериментальных данных, используемых для оценивания;
• ‘mdl’ – возвращается структура, обеспечивающая минимум так называемого критерия Риссанена минимальной длины описания (Rissanen’s Minimum Description Lngth – MDL):
;
• при строковой переменной, равной некоторому численному значению a
, выбирается структура, которая минимизирует:
;
• vmod – значение соответствующего критерия.
Выбор наилучшей структуры порядков полиномов можно осуществить и с помощью более простой команды:
>> nn=selstruc(v,0)
MATLAB возвращает:
nn =
8 2 1
С учетом выбранной структуры модели определим вид модели ARX, выполнив функцию arx
:
>> darx=arx(zdanv,nn)
Возвращается матрица из 100 столбцов и 4 строк с значениями различных критериев: vmod =
vmod =
Columns 1 through 8
0.0107 0.0078 0.0080 0.0078 0.0079 0.0079 0.0079 0.0079
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 9 through 16
0.0079 0.0079 0.0084 0.0072 0.0073 0.0073 0.0072 0.0072
1.0000 1.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
9.0000 10.0000 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 17 through 24
0.0073 0.0073 0.0073 0.0073 0.0079 0.0070 0.0072 0.0072
2.0000 2.0000 2.0000 2.0000 3.0000 3.0000 3.0000 3.0000
7.0000 8.0000 9.0000 10.0000 1.0000 2.0000 3.0000 4.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 25 through 32
0.0072 0.0072 0.0072 0.0072 0.0073 0.0073 0.0079 0.0071
3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 4.0000 4.0000
5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 1.0000 2.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 33 through 40
0.0072 0.0072 0.0072 0.0072 0.0073 0.0073 0.0073 0.0073
4.0000 4.0000 4.0000 4.0000 4.0000 4.0000 4.0000 4.0000
3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 41 through 48
0.0080 0.0071 0.0071 0.0072 0.0071 0.0072 0.0073 0.0074
5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 49 through 56
0.0074 0.0074 0.0080 0.0070 0.0071 0.0071 0.0071 0.0071
5.0000 5.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000
9.0000 10.0000 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 57 through 64
0.0073 0.0073 0.0073 0.0074 0.0080 0.0070 0.0071 0.0071
6.0000 6.0000 6.0000 6.0000 7.0000 7.0000 7.0000 7.0000
7.0000 8.0000 9.0000 10.0000 1.0000 2.0000 3.0000 4.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 65 through 72
0.0071 0.0071 0.0073 0.0074 0.0074 0.0074 0.0080 0.0070
7.0000 7.0000 7.0000 7.0000 7.0000 7.0000 8.0000 8.0000
5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 1.0000 2.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 73 through 80
0.0071 0.0071 0.0071 0.0071 0.0073 0.0074 0.0074 0.0074
8.0000 8.0000 8.0000 8.0000 8.0000 8.0000 8.0000 8.0000
3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 81 through 88
0.0080 0.0070 0.0071 0.0071 0.0071 0.0071 0.0073 0.0074
9.0000 9.0000 9.0000 9.0000 9.0000 9.0000 9.0000 9.0000
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 89 through 96
0.0074 0.0074 0.0080 0.0070 0.0071 0.0071 0.0071 0.0072
9.0000 9.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000
9.0000 10.0000 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 97 through 100
0.0073 0.0074 0.0074 0.0074
10.0000 10.0000 10.0000 10.0000
7.0000 8.0000 9.0000 10.0000
1.0000 1.0000 1.0000 1.0000
Возвращается дискретная модель, представленная в тета - формате (внутренним видом матричных моделей).
Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t)
A(q) = 1 - 1.01 q^-1 + 0.3552 q^-2 - 0.03471 q^-3 - 0.1432 q^-4
+ 0.1302 q^-5 - 0.0128 q^-6 - 0.08582 q^-7 + 0.06296 q^-8
B(q) = 0.1367 q^-1 + 0.07335 q^-2
Estimated using ARX from data set zdanv
Loss function 0.00666153 and FPE 0.00693343
Sampling interval: 0.08
Функция armax
оценивает параметры ARMAX модели:
>> darmax = armax(zdanv,[2 2 2 1])
Аргументы функции:
zdanv – вектор экспериментальных данных; [na nb nc nk] – степени полиномов и величина задержки.
Возвращается дискретная модель, представленная в тета – формате:
Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + C(q)e(t)
A(q) = 1 - 0.8733 q^-1 + 0.1567 q^-2
B(q) = 0.1331 q^-1 + 0.1028 q^-2
C(q) = 1 + 0.1854 q^-1 - 0.01339 q^-2
Estimated using ARMAX from data set zdanv
Loss function 0.00787129 and FPE 0.00806524
Sampling interval: 0.08
Функция oe
оценивает параметры ОЕ модели:
>> zoe=oe(zdanv,[2 2 1])
Возвращается дискретная модель, представленная в тета – формате:
Discrete-time IDPOLY model: y(t) = [B(q)/F(q)]u(t) + e(t)
B(q) = 0.1478 q^-1 + 0.1052 q^-2
F(q) = 1 - 0.8219 q^-1 + 0.1102 q^-2
Estimated using OE from data set zdanv
Loss function 0.020577 and FPE 0.0209102
Sampling interval: 0.08
Функция bj
оценивает параметры модели Бокса-Дженкинса:
>> zbj=bj(zdanv,[2 2 2 2 1])
Возвращается дискретная модель, представленная в тета – формате:
Discrete-time IDPOLY model: y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)
B(q) = 0.1334 q^-1 + 0.101 q^-2
C(q) = 1 - 0.1222 q^-1 - 0.1405 q^-2
D(q) = 1 - 1.148 q^-1 + 0.3494 q^-2
F(q) = 1 - 0.8958 q^-1 + 0.1813 q^-2
Estimated using BJ from data set zdanv
Loss function 0.00699912 and FPE 0.0072286
Sampling interval: 0.08
Функция n
4
sid
используется для оценивания параметров моделей переменных состояния в канонической форме при произвольном числе входов и выходов:
[zn4s,AO] = n4sid(z,order,ny,auxord),
где: z – матрица экспериментальных данных
order – задает порядок модели. Если данный аргумент вводится как вектор – строка, то предварительные расчеты выполняются по моделям всех заданных порядков (по умолчанию от первого по десятый), с выводом графика, позволяющего выбрать оптимальный порядок. Если order = ‘best’(по умолчанию), то выбирается модель наилучшего порядка;
ny – количество выходов (по умолчанию ny = 1);
auxord – дополнительный порядок, используемый алгоритмом оценивания. Данный порядок должен быть больше, чем порядок, задаваемый параметром order (по умолчанию дополнительный порядок равен (1.2*order+3)). Если данный аргумент вводится как вектор – строка, то выбирается модель наилучшего порядка.
Для рассматриваемого примера project24 имеем:
>>zn4s=n4sid(zdanv,[1:10],[1:10]),
где в первых квадратных скобках задается интервал порядков модели order, и предварительные расчеты выполняются по моделям для всех заданных порядков от 1 до 10 с выводом графика, позволяющего выбрать оптимальный порядок. После этого необходимо в командной строке MATLAB набрать этот порядок и продолжить вычисление коэффициентов модели, нажав клавишу enter (рис. 2. 9). Во вторых квадратных скобках задается так называемый дополнительный порядок, используемый алгоритмом оценивания (по умолчанию дополнительный порядок равен (1.2*order+3)). При этом выбирается оптимальный порядок без вывода соответствующего графика.
Результатом выполнения команды является вывод результатов оценивания:
Warning: Input arguments must be scalar.
> In n4sid>transf at 1044
In n4sid at 134
Select model order:('Return' gives default)
При нажатии In n4sid>transf at 1027, In n4sid at 134 (синий цвет: имя модели, построенной в тета - формате) появится окно редактора М-файла программы.
При нажатии Enter появится Order chosen to 3 (Закажите выбранный 3)
State-space model: x(t + Ts) = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
Рис.
2. 9. График для выбора оптимального порядка модели
A =
x1 x2 x3
x1 0.77993 0.54376 -0.013931
x2 -0.12996 0.073872 0.19929
x3 -0.1067 0.016064 -0.47392
B =
расход газа
x1 -0.022796
x2 -0.048006
x3 -0.034112
C =
x1 x2 x3
температура -4.6742 -0.54702 0.0027609
D =
расход газа
температура 0
K =
температура
x1 -0.20139
x2 0.075372
x3 0.02383
x(0) =
x1 0.10598
x2 0.033411
x3 -0.0020287
Estimated using N4SID from data set zdanv
Loss function 0.00753932 and FPE 0.00791011
Sampling interval: 0.08
Функция pem
оценивает параметры обобщенной многомерной линейной модели:
>> zpem=pem(zdanv)
State-space model: x(t+Ts) = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3 x4 x5
x1 0.79771 0.52932 0.017441 -0.043818 -0.055481
x2 -0.089871 -0.056516 0.34751 0.4433 -0.14175
x3 0.12179 -0.31736 -0.29564 -0.43427 0.29463
x4 0.096387 -0.0086843 -0.063652 0.64711 -0.21098
x5 -0.012304 0.025405 -0.18701 0.69982 1.0514
B =
расход газа
x1 -0.021221
x2 -0.053799
x3 -0.040284
x4 -0.00059208
x5 -0.032983
C =
x1 x2 x3 x4 x5
температура -4.675 -0.54794 0.010925 0.068154 -0.1202
D =
расход газа
температура 0
K =
температура
x1 -0.22248
x2 0.037523
x3 0.024536
x4 0.11512
x5 -0.068061
x(0) =
x1 0.10985
x2 0.0039575
x3 0.071289
x4 -0.15615
x5 -0.15783
Estimated using PEM from data set zdanv
Loss function 0.00664935 and FPE 0.00720346
Sampling interval: 0.08
2
. 9. Функции преобразования моделей
Для дальнейшего использования, полученных моделей при анализе и синтезе систем автоматизации технологических процессов в пакете System Identification Toolbox MATLAB имеются специальные функции, позволяющие выполнить преобразование этих моделей из тета-формата (внутреннего вида матричных моделей, являющегося дискретным) в другие виды и в частности из дискретной в непрерывную модель в виде передаточной функции.
Функция th2a
rx
преобразует модель тета-формата в ARX модель:
Функция имеет синтаксис:
>> [A,B]=th2arx(darx)
где darx - модель тета-формата
A =
Columns 1 through 8
1.0000 -1.0105 0.3552 -0.0347 -0.1432 0.1302 -0.0128 -0.0858
Column 9
0.0630
B =
0 0.1367 0.0733
Функция th2ff
вычисляет частотные характеристики и соответствующие стандартные отклонения по модели в тета-формате.
В качестве аргумента функции может выступать любая из рассмотренных ранее моделей, например darx:
>> [g,phiv]=th2ff(darx)
IDFRD model g.
Contains Frequency Response Data for 1 output and 1 input
and SpectrumData for disturbances at 1 output
at 129 frequency points, ranging from 0.1 rad/s to 39.27 rad/s.
Output Channels: температура
Input Channels: расход газа
Sampling time: 0.08
Estimated from data set zdanv using ARX.
IDFRD model phiv.
Contains SpectrumData for 1 signal
at 126 frequency points, ranging from 0.1 rad/s to 39.27 rad/s.
Output Channels: температура
Sampling time: 0.08
Estimated from data set zdanv using ARX.
Функция th
2
poly
преобразует матрицу модели тета-формата в матрицы обобщенной (многомерной) линейной модели:
>> [A,B,C,D,K,lan,T]=th2poly(zpem)
A = 1.0000 -2.1441 1.5148 -0.3280 0.1302 -0.1081
B = 0 0.1322 -0.0698 -0.0946 0.0197 0.0668
C = 1.0000 -1.1083 -0.0108 0.1446 0.3438 -0.0371
D =1
K = 1
lan = 0.0069
T = 0.0800
Здесь параметр lan определяет интенсивность шума наблюдений.
Функция th2ss
преобразует тета-модель в модель для переменных состояния. В качестве аргумента функции может выступать любая из рассмотренных ранее моделей, например darmax:
>> [A,B,C,D,K,x0]=th2ss(darmax)
A =
0.8733 1.0000
-0.1567 0
B =
0.1331
0.1028
C =
1 0
D =
0
K =
1.0587
-0.1701
x0 =
0
0
Функция th2tf
преобразует модель тета-формата многомерного объекта в вектор передаточных функций, связанных с выбранным входом:
>> [num,den]=th2tf(zn4s)
num = 0 0.1327 0.1566 0.0575
den = 1.0000 -0.3799 -0.2810 0.0749
Команда tf
служит для представления передаточной функции в виде отношения:
>> zzn4s=tf(num,den,0.08)
Возвращает:
Transfer function:
0.1327 z^2 + 0.1566 z + 0.0575
------------------------------------
z^3 - 0.3799 z^2 - 0.281 z + 0.07493
Sampling time: 0.08
Функция thd
2
thc
преобразует дискретную модель в непрерывную
Например: преобразуем дискретную модель тета-формата zn4s (модель переменных состояния в канонической форме при произвольном числе входов и выходов) в непрерывную модель и представим ее в виде передаточной функции. Для этого необходимо сначала выполнить функцию thd
2
thc
преобразующую дискретную модель в непрерывную, затем выполнить функцию th2tf
преобразующую модель тета-формата многомерного объекта в вектор передаточных функций, связанных с выбранным входом, а затем команду tf
для представления передаточной функции в виде отношения:
>> sn4s=thd2thc(zn4s);
>> [num,den]=th2tf(sn4s);
>> sysn4s=tf(num,den)
Transfer function:
Transfer function:
-0.891 s^2 + 77.33 s + 746.9
---------------------------------
s^3 + 32.39 s^2 + 308.9 s + 891.7
Для обратного преобразования непрерывной модели в дискретную модель существует функция thc2thd.
Функция th2zp
рассчитывает нули, полюса и статические коэффициенты передачи (коэффициенты усиления) модели тета - формата zn4s многомерного объекта:
>> [zepo,k]=th2zp(zn4s)
zepo =
1.0000 61.0000 21.0000 81.0000
-0.5898 + 0.2921i 0.2754 + 0.1282i -0.4946 0.6660
-0.5898 - 0.2921i 0.3531 0.6365 0.0651
Inf Inf 0.2380 0.2121
k =
1.0000
0.8376
0.0648
Информацию о нулях и полюсах модели zn4s
можно получить с помощью команды
>> [zero,polus]=getzp(zepo)
zero =
-0.5898 + 0.2921i
-0.5898 - 0.2921i
polus =
-0.4946
0.6365
0.2380
С помощью команды zpplot
можно построить график нулей и полюсов модели zn4s
>>zpplot(zpform(zepo))
На рис. 2. 10. представлен график нулей (обозначены кружком) и полюсов (обозначены крестиком) модели zn4s, который получен с помощью команды zpplot.
Рис.
2
.10. Графики нулей и полюсов модели zn4s
Данные графика показывают, что модель является устойчивой: полюса модели находятся внутри окружности с радиусом, равным 1 и проходящей через точку с координатами (–1; j0).
2. 10
. Проверка адекватности модели
Одним из важных этапов идентификации объектов автоматизации является проверка качества модели по выбранному критерию близости выхода модели и объекта, т.е проверка ее адекватности. В пакете System Identification Toolbox MATLAB в качестве такого критерия принята оценка адекватности модели fit, которая рассчитывается по формуле: fit = norm (yh – y)/
, где norm – норма вектора; yh и y – выходы модели и объекта соответственно; N – количество элементов массива данных.
Для проверки адекватности полученных ранее моделей воспользуемся функцией:
>> compare(zdane,zn4s,zpem,zoe,zbj,darx,darmax).
где: zdane – выход объекта;
zn4s,zpem,zoe,zbj,darx,darmax – выходы моделей zn4s,zpem,zoe,zbj,darx,darmax
Рис.
2
. 11. Графики выходов объекта и моделей.
Результатом выполнения команды является вывод графика выходов объекта и построенных моделей (Рис. 2. 11). На графике цветными линиями представлены выходы полученных моделей и значения критерия адекватности, выраженного в процентах. Наилучшие показатели имеют модели darx, zn4s и zpem.
Для проверки адекватности модели zn4s воспользуемся функцией:
>>compare(zdane,zn4s)
Результат выполнения команды является вывод графика объекта на рис. 2. 12.
Рис.
2
. 12. Графики выходов объекта и модели
zn4s.
В пакете System Identification Toolbox MATLAB имеется возможность прогнозировать ошибку моделирования при заданном входном воздействии u
(
t
)
и известной выходной координате объекта y
(
t
)
. Оценивание производится методом прогноза ошибки Preictive Error Method, сокращенно PEM, который заключается в следующем. Пусть модель исследуемого объекта имеет вид так называемой обобщенной линейной модели
y
(t
) = W
(z
) u
(t
) + v
(t
),
где W
(z
) – дискретная передаточная функция любой из ранее рассмотренных моделей. При этом шум v
(t
) может быть представлен как
v
(t
) = H
(z
) e
(t
),
где e
(z
) – дискретный белый шум, который собственно и характеризует ошибку модели; H
(z
) – некоторый полином от z
, приводящий дискретный белый шум к реальным помехам при измерении выходных параметров объекта.
Из данных выражений следует, что
e
(t
) = H
-1
(z
) [y
(t
) – W
(z
) u
(t
)].
Функция resid вычисляет остаточную ошибку e для заданной модели, а также r – матрицу значений автокорреляционной функции процесса e
(t
) и значения взаимокорреляционой функции между остаточными ошибками e
(t
) и выходами объекта автоматизации y
(t
) вместе с соответствующими 99 %-ми доверительными коридорами. Кроме указанных значений выводятся графики данных функций. В качестве примера сравним остаточные ошибки и соответствующие корреляционные функции для полученных моделей darx и zbj, имеющих максимальную и минимальную оценки адекватности с помощью команд:
>> [e,r]=resid(zdan,darx)
>> [e1,r1]=resid(zdan,zbj)
Приведенные графики (рис. 2. 13 и 2 14) характеризуют равномерное распределение остаточных ошибок во всем диапазоне изменения интервалов времени τ
. Причем значения остаточных ошибок
для модели darx практически в два раза больше, чем для модели zbj. Для вывода графиков необходимо выполнить команду resid(r).
Рис. 2. 13. График автокорреляционной и взаимокорреляционной функций для модели
zbj
Рис.
2
. 14. График автокорреляционной и взаимокорреляционной функций для модели
darx
После выполнения функции:
[e,r]=resid(zdan,darx)
MATLAB возвращает:
Time domain data set with 1097 samples.
Sampling interval: 0.08
Outputs Unit (if specified)
e@температура гр.С 100
Inputs Unit (if specified)
u1
r =
1.0e+003 *
Columns 1 through 8
0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000
0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000
0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000
0.0002 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000
Columns 9 through 16
-0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000
-0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000
-0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000
-0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000
Columns 17 through 24
-0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000
-0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000
0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000
0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000
Columns 25 through 27
0.0000 1.0970 0.0010
-0.0000 0 0
0.0000 0 0
-0.0000 0 0
После выполнения команды >> resid(r) выводится график автокорреляционной и взаимокорреляционной функций для модели.
Таким образом, в ходе оценки адекватности различных моделей объекта автоматизации технологического процесса тепловой обработки материалов определены модели darx, zn4s и zpem, значения критерия адекватности которых максимальны и, следовательно, могут быть использованы в дальнейшем при анализе и синтезе систем автоматизации.
2. 11
. Анализ модели технического объекта управления
Для анализа модели ТОУ возьмем модель zn4s, имеющую один из наилучших показателей адекватности.
• zzn4s – дискретная модель в виде передаточной функции
0.1327 z^2 + 0.1566 z + 0.0575
------------------------------------
z^3 - 0.3799 z^2 - 0.281 z + 0.07493
• sysn4s – непрерывная модель в виде передаточной функции
-0.891 s^2 + 77.33 s + 746.9
---------------------------------
s^3 + 32.39 s^2 + 308.9 s + 891.7
Приведенные виды являются одной и той же моделью, записанной в разных формах и форматах. Проанализируем динамические характеристики модели. Построим переходную характеристику ТОУ для дискретной и непрерывной моделей и определим основные показатели переходного процесса. Для этого можно воспользоваться функцией step
. Функция step рассчитывает и строит реакцию модели на единичную ступенчатую функцию, т. е. возвращает переходную функцию системы:
step(sys)
step(sys, t)
step(sys1,sys2,….,sysN, t)
step(sys1,’PlotStyle1’,….,sysN, ’PlotStyleN’)
[y,t,x] = step(sys)
Д
ля моделей, заданных в пространстве состояний, начальные условия принимаются нулевыми. Аргументы функции следующие:
- sys,sys1,sys2,….,sysN – имена моделей для которых строятся переходные функции;
- t – аргумент, задающий момент окончания моделирования – либо в форме t = Tfinal (в секундах), либо в форме t = 0:dt:Tfinal.
Для дискретных моделей значение dt должно равняться интервалу дискретизации, для непрерывных моделей – быть достаточно малым, чтобы учесть наиболее быстрые изменения переходного процесса;
- ’PlotStyle1’,….,’PlotStyleN’ – строковые переменные, задающие стили (типы линий) при выводе нескольких графиков одновременно.
Возвращаемые величины:
- графики переходных процессов;
- y, x, t – соответственно, векторы, содержащие значения переходного процесса, переменных состояния и моментов времени (при возвращении данных величин график переходного процесса не отображается).
Выполним построение переходной характеристики ТОУ, представленной дискретной zzn4s инепрерывной sysn4s моделями и определим основные показатели переходного процесса, используя функцию step:
>>step(zzn4s,sysn4s)
После выполнения команды step MATLAB возвращает графики переходного процесса (Рис. 2. 15). Нажатие левой клавиши мыши в любом месте на графике переходного процесса приводит к появлению всплывающей информационной подсказки о величине текущего численного значения переходного процесса и моменте времени.
Нажатие правой клавиши в любом месте на графике переходного процесса приводит к появлению всплывающего меню редакции окна всплывающей информационной подсказки.
Рис.
2
. 15. Графики переходных процессов модели
z
zn4s и
sy
sn4s
На графиках переходных процессов ступенчатой линией представлен переходной процесс дискретной модели, а сплошной линией – непрерывной модели. Кроме того, в поле графика указаны основные характеристики переходного процесса:
• время регулирования (Setting time) – 0,769 с для обоих моделей;
• установившееся значение выходной координаты – 0,838 для обеих моделей.
Для построения импульсной характеристики моделей необходимо воспользоваться командой:
>>impulse(zzn4s,sysn4s).
После выполнения команды impulse MATLAB возвращает графики (Рис. 2. 16).
Основными характеристиками модели ТОУ при подаче на вход единичного импульсного воздействия являются:
• пиковая амплитуда (Peak amplitude) составляет для дискретной модели 0,207 а для непрерывной – 2,79.
• время регулирования составляет для дискретной модели 0,922 и для непрерывной модели – 0,863 с.
Для определения статического коэффициента усиления модели ТОУ можно использовать команду dcgain
:
>> k=dcgain(sysn4s)
После выполнения команды получим: k = 0.8376.
Рис.
2
. 16. Графики импульсной характеристики
Для определения частотной характеристики моделей используем команду bode:
Рис.2. 17. Частотные характеристики моделей
Выполним построение частотной характеристики ТОУ, представленной дискретной zzn4s и непрерывной sysn4s моделями (Рис. 2. 17).
Н
а графиках частотных характеристик указаны значения запасов устойчивости по амплитуде (Gain Margin), которые для дискретной модели составляет 29,7 dB, а для непрерывной модели – бесконечность.
Значения запасов устойчивости можно определить также и в режиме командной строки MATLAB с помощью команд:
>> [Gm,Pm,Wcg,Wcp]=margin(sysn4s) – для непрерывной модели:
MATLAB возвращает:
Gm =
26.5077
Pm =
Inf
Wcg =
48.5667
Wcp =
NaN
>> [Gm1,Pm1,Wcg1,Wcp1]=margin(zzn4s) – для дискретной модели:
MATLAB возвращает:
Gm1 =
9.0385
Pm1 =
Inf
Wcg1 =
21.0461
Wcp1 =
NaN
где Gm – запас устойчивости по амплитуде в натуральных величинах на частоте Wcg, Pm – запас устойчивости по фазе на частоте Wcp.
Для определения запасов устойчивости в логарифмическом масштабе необходимо выполнить следующие операции:
>> Gmlog=20*log10(Gm1) – для дискретной модели:
Gmlog =
19.1219
>> Gmlog=20*log10(Gm) – для непрерывной модели:
Gmlog =
28.4675
Как видно, определение запасов устойчивости последним способом позволяет значительно точнее вычислять эти значения, чем на графиках частотных характеристик. Анализ частотных характеристик показывает, что модели zzn4s и sysn4s являются устойчивыми с соответствующими запасами устойчивости по амплитуде. Запас устойчивости по фазе равен бесконечности.
Этот вывод подтверждается так же комплексной амплитудно-фазовой характеристикой АФХ (называется диаграммой Найквиста, Рис. 2. 18), так как годограф АФХ не пресек
ает точку комплексной плоскости с координатами –1, j0.
Рис.
2
. 18. Годограф АФХ для непрерывной и дискретной моделей
Для построения АФХ необходимо воспользоваться командой:
>>nyquist(zzn4s,sysn4s),
Определить устойчивость моделей можно с помощью карты нулей и полюсов по расположению нулей моделей относительно окружности с единичным радиусом на комплексной плоскости, как это было показано на рис. 2. 10. Построить карту нулей и полюсов моделей можно так же с помощью команды pzmap(zzn4s,sysn4s), либо – pzmap(zn4s,sn4s).
Построим график изменения e
(t
) и определим основные статистические характеристики помехи с помощь команды plot (e) (Рис. 2. 19).
Для получения статистических характеристик необходимо в строке меню графика в позиции Tools
выбрать опцию Data
statistics
. Результатом выполнения команды явится окно, в котором будут указаны основные статистические характеристики случайного процесса изменения во времени e
(t
),(Рис. 2. 20), к которым относятся:
• min и max – минимальное и максимальное значения помехи.
Для нашего случая – -0,2373 и 0,2086 соответственно;
• mean – арифметическое среднее значение (0,001403);
• median – медиана процесса (0,003994);
• std – среднеквадратическое отклонение (0,0805);
• range – диапазон изменения помехи от минимального до максимального значения (1.12).Во всех случаях размерность аддитивной помехи такая же, как и выходная величина объекта автоматизации – о
С.
Рис. 5. 19. График аддитивной помехи
e
(
t
)
Рис. 5. 20. Статистические характеристики
e
(
t
)
Полученные статистические характеристики помехи могут быть полезны в дальнейшем при синтезе системы автоматического регулирования температуры теплового объекта автоматизации.
Для решения задач анализа и синтеза систем управления важно знать ответ на другой не менее важный вопрос, чем полученные временные, частотные и статистические характеристики: обладает ли объект свойством управляемости в смысле возможности его перевода из заданной начальной точки (или области) в заданную конечную точку (или область). Решение проблемы управляемости основано на анализе уравнений переменных состояния вида:
,
,
где A
,
B
,
C
,
D
– матрицы соответствующих размеров, v
(t
) – коррелированный белый шум наблюдений. Возможна и другая (так называемая обновленная или каноническая) форма представления данной модели:
,
,
где К – некоторая матрица (вектор столбец), е(t) – дискретный белый шум (скаляр),
и формулируется следующим образом: объект называется вполне управляемым, если выбором управляющего воздействия u(t) на интервале времени [t
0
,
tk
] можно перевести его из любого начального состояния
y
(
t
0
) в произвольное заранее заданное конечное состояние
y
(
tk
).
Критерием управляемости линейных стационарных объектов является условие: для того чтобы объект был вполне управляем, необходимо и достаточно, чтобы ранг матрицы управляемости
MU
= (B AB A2
B
… An-1
B)
равнялся размерности вектора состояний n
rang MU
=
n
.
В пакете Control System Toolbox имеется функция ctrb
, формирующая матрицу управляемости в пространстве состояний. Для того, чтобы воспользоваться этой функцией необходимо вычислить матрицы A
, B
, C
, D
с помощью команды:
>>[A,B,C,D]=ssdata(sn4s)
A =
-0.8930 16.3384 4.0253
-4.7215 -22.0535 -3.5128
-1.0484 -2.5116 -9.4429
B =
0.3680
-1.5178
-0.3597
C =
-4.6742 -0.5470 0.0028
D =
0
Следует обратить внимание, что для расчета матриц используется непрерывная модель, так как дискретная модель имеет другие значения, а в критерии управляемости используются матрицы линейных непрерывных стационарных объектов.
Вычислим матрицу управляемости:
>> Mu=ctrb(A,B)
Mu =
0.3680 -26.5754 590.3514
-1.5178 32.9991 -626.2378
-0.3597 6.8234 -119.4511
Определим ранг матрицы управляемости:
>> n=rank(Mu)
n =
3.
Таким образом, для исследуемой модели объекта размерность вектора состояний, определяемая размером матриц A
и B
равна трем и ранг матрицы управляемости MU
также равен трем, что позволяет сделать вывод о том, что объект автоматизации является вполне управляемым, т.е. для него имеется такое управляющее воздействие u(t)
, которое способно перевести на интервале времени [t0
,
tk
] объект из любого начального состояния y
(
t
0
)
в произвольное заранее заданное конечное состояние y
(
tk
)
.
При синтезе оптимальных систем с обратной связью сами управления получаются как функции от фазовых координат. В общем случае фазовые координаты являются абстрактными величинами и не могут быть исследованы. Поддается измерению (наблюдению) вектор y
=
(y
1
, …, yk
)T
, который обычно называют выходным вектором или выходной переменной, а его координаты – выходными величинами. Выходная переменная функционально связана с фазовыми координатами, и для реализации управления с обратной связью необходимо определить фазовые координаты по измеренным значениям выходной переменной. В связи с этим возникает проблема наблюдаемости, заключающаяся в установлении возможности состояния определения состояния объекта (фазового вектора) по измеренным значениям выходной переменной на некотором интервале.
Решение проблемы наблюдаемости основано на анализе уравнений переменных состояния вида
или Y
(
p
) =
W
(
p
)*
U
(
p
)
и формулируется следующим образом: объект называется вполне наблюдаемым, если по реакции y(t
1
) на выходе объекта
, на интервале
времени
[t
0
, t
1
] при заданном управляющем воздействии u(t) можно определить начальное состояние вектора переменных состояния x(t), являющихся фазовыми координатами объекта.
Критерием наблюдаемости линейных стационарных объектов является
условие: для того, чтобы объект был вполне наблюдаемым, необходимо и достаточно, чтобы ранг матрицы наблюдаемости
М
Y
= (CT
AT
CT
(AT
)2
CT
… (
AT
)
n
-1C
)
равнялся размерности вектора состояния
n =
rang M
Y
.
Определим матрицу наблюдаемости и ее ранг с помощью функций пакета Control System Toolbox:
>> My=obsv(A,C)
My =
1.0e+003 *
-0.0047 -0.0005 0.0000
0.0068 -0.0643 -0.0169
0.3154 1.5712 0.4129
>> n=rank(My)
n =3
Таким образом, для исследуемой модели объекта размерность вектора состояний, определяемая размером матриц A
и С
равна трем и ранг матрицы наблюдаемости M
Y
также равен трем, что позволяет сделать вывод о том, что объект автоматизации является вполне наблюдаемым, т.е. для него всегда можно определить по значениям выходной величины y
(
t
)
вектор переменных состояния, необходимый для синтеза системы управления.
2. 12
. Основные результаты идентификации технического объекта автоматизации
Идентификация распылительной сушилки проводилась с целью получения модели объекта, необходимой для синтеза системы автоматизации и получения основных характеристик объекта автоматизации.
В результате проведенного эксперимента был получен массив данных, состоящий из 1097 значений входного параметра распылительной сушилки – расхода газа в м3
/час и 1097 значений выходного параметра – температуры отходящих газов в градусах Цельсия, измеренных через временные промежутки 0, 08 с.
В ходе идентификации были получены следующие результаты:
1. Обработаны и преобразованы данные в единый файл, содержащий необходимую информацию о входных и выходных параметрах объекта, их значениях и размерностях измерения. Получены графические зависимости изменения температуры отходящих газов от расхода горючего газа на входе распылительной сушилки.
2. Проведено непараметрическое оценивание исходных данных для определения статистических характеристик массивов исходных данных.
3. В результате параметри
ческого оценивания экспериментальных данных, проведенного с целью определения параметров модели заданной структуры путем минимизации выбранного критерия качества модели, были получены различные структуры и виды моделей распылительной сушилки:
– модель авторегрессии;
– модель авторегрессии с дополнительным входом;
– модель авторегрессии скользящего среднего;
– модель «вход-выход»;
– модель Бокса-Дженкинса;
– модель для переменных состояния.
4. Проверка адекватности моделей показала, что наилучшей степенью адекватности (55.28%) обладает модель для BJ. Получены значения автокорреляционной функции ошибок процесса и значения взаимокорреляционой функции между остаточными ошибками и выходами объекта автоматизации вместе с соответствующими 99 %-ми доверительными коридорами.
5. Проведенное преобразование моделей позволило получить вид передаточных функций распылительной сушилки в дискретном и непрерывном видах, необходимых для дальнейшего анализа и синтеза системы автоматизации:
0.1327 z^2 + 0.1566 z + 0.0575
W(z) = ----------------------------------------------
z^3 - 0.3799 z^2 - 0.281 z + 0.07493
-0.891 s^2 + 77.33 s + 746.9
W(s) = ----------------------------------------
s^3 + 32.39 s^2 + 308.9 s + 891.7
Проведенный анализ модели распылительной сушилки позволил определить основные статические и динамические характеристики объекта автоматизации:
– время регулирования – 0,863 с;
– запас устойчивости по амплитуде – 29,7 дБ;
– запас устойчивости по фазе
– бесконечность.
7. Анализ управляемости и наблюдаемости объекта автоматизации показал, что распылительная сушилка является вполне управляемой и наблюдаемой. На неё можно подавать управляющие воздействия для перевода её из одного начального состояния в произвольное заранее заданное конечное состояние и для этого заранее заданного управляющего воздействия можно определить (измерить) начальное состояние вектора переменных состояния.
ГЛАВА 3. ПОСТРОЕНИЕ СИСТЕМЫ УПРАВЛЕНИЯ ТЕХНОЛОГИЧЕСКИМ
ПРОЦЕССОМ
3.1.
Задание структуры системы автоматического управления, проверка системы управления на устойчивость
Функция rltool
открывает графический интерфейс, позволяющий проектировать корректирующее звено в замкнутой одномерной системе управления методом корневого годографа. Функция имеет следующий синтаксис:
rltool(sys,comp,LocationFlag,FeedbackSign)
где: sys – имя модели одномерного объекта;
comp – имя корректирующего звена (компенсатора);
LocationFlag – переменная, задающая позицию компенсатора в системе: 1 – в прямом тракте системы, 2 – в цепи обратной связи;
FeedbackSign – тип обратной связи: -1 – отрицательная обратная связь, 1 - положительная обратная связь.
Но удобнее работать с окном интерфейса SISO Design for System FeedbackConfig.
Выполнение функции rltool
без аргументов приводит к появлению
основного окна интерфейса SISO Design for System FeedbackConfig, реализующего метод корневого годографа (Рис. 3. 13). Кнопка +/- позволяет установить отрицательную (-) обратную связь. Кнопка FS позволяет выбрать структуру системы и позицию компенсатора в системе.
Выберем структуру с расположением компенсатора «С» в прямом тракте системы (Рис. 3. 13) и отрицательную обратную связь. Далее необходимо задать модели звеньев структурной схемы: «F», «C», «G», «H». Для этого следует сделать следующее:
В меню окна интерфейса SISO Design for System FeedbackConfig необходимо выполнить команды: File
Import
.
В открывшемся окне (Рис. 3. 14) загрузки модели Import
System
Data
выберем модель sysn4s. Переключатель Import from указывает, из какой области загружается модель. Модель sysn4s находится в рабочей области MATLAB, т. е. в Workspace. В поле System Data окна Import System Data (Рис. 3. 14) обозначена структурная схема замкнутой системы. В ней «F», «G», «H» звенья модели которые можно загружать. Звено, обозначенное буквой «С», это то компенсирующее динамическое звено структуру и параметры которого нужно определить.
Рис
. 3. 13.
Окна
интерфейса
SISO Design for System FeedbackConfig
.
Рис.
3
. 14. Окно загрузки модели
Import
System
Data
Далее необходимо выполнить загрузку модели технического объекта управления sysn4s в звено «G» нажатием кнопки со стрелкой, указывающей на звено «G». Модели звеньев «F» и «H» выберем по умолчанию (это пропорциональные звенья с единичным коэффициентом передачи). Сохраним, полученную модель под именем mysys1. Подтвердим сохранение нажатием кнопки ОК. Окно загрузки при этом закроется, а основное окно интерфейса приобретет вид, показанный на рис. 3. 15.
Рис
. 6. 15.
Основное
окно
интерфейса
SISO Design for System mysys1
В графической части окна на комплексной плоскости нулей и полюсов отображены полюсы и нули замкнутой системы, а также линии их перемещения при изменении коэффициента передачи компенсатора от заданного значения до бесконечности. Система имеет три полюса и два нуля (это подтверждается видом аналитического выражения передаточной функции ТОУ, которую можно просмотреть, если щелкнуть ЛК на блоке «G» структурной схемы замкнутой системы и в открывшемся окне System
Data
в поле Plant
Model
:
sy
sn4s
нажать кнопку Show Transfer Function). Передаточная функция имеет в числителе полином второй степени, а в знаменателе полином третьей степени.
Из расположения полюсов (крестики) на комплексной плоскости следует, что замкнутая система достаточно устойчива, так как все три полюса находятся в левой полуплоскости и достаточно далеко от границы
устойчивости. В этом еще можно убедиться, просмотрев график переходного процесса замкнутой системы, если в меню Analysis
выполнить команду Response to Step Command. Данный выбор приведет к открытию окна интерактивного обозревателя LTI
-
Viewer
for
SISO
Design
Tool
(Рис. 3. 16). Как видно из графика сходящегося переходного процесса (кривая r to y) время переходного процесса достаточно мало (с данным корректирующим звеном пропорционального типа с коэффициентом пропорциональности равным единице). Система устойчива.
Однако следует отметить, что при явной устойчивости системы наблюдается некоторое перерегулирование переходного процесса. Следовательно, можно попытаться скорректировать переходный процесс, сделав его апериодическим, т. е улучшить динамические свойства системы. Сделать это можно путем подбора передаточной функции компенсирующего звена «С».
Рис. 3. 16. Графики переходных процессов в системе
расположенных над графическим окном.
В поле Current Compensator окна SISO Design for System mysys4 отразится текущая передаточная функция компенсатора. Необходимо также помнить, что после выполнения меню Analysis в произведенной сессии дальнейшие изменения в структуре системы не будут отражены в результатах повторного выполнения меню Analysis. Поэтому для дальнейшего анализы при коррекции системы необходимо загружать новое окно интерфейса, выполнив повторно в режиме командной строки функцию rltool
.
Далее необходимо выполнить анализ построенной замкнутой системы управления с целью определения ее параметров и, сравнив их с заданными в техническом задании параметрами, сделать вывод о необходимости корректировки системы или убедиться в отсутствии такой необходимости. Просмотреть все характеристики можно выполнив в меню Analysis окна SISO Design for System mysys4 команды: Response to Step Command; Rejection of Step Disturbance; Closed-Loop Bode; Compensator Bode; Open-Loop Nyquist. После выполнения команд появится окно обозревателя LTIViewer (Рис. 3. 17)
Рис. 3
. 17
. О
кно обозревателя
LTIViewer
При выполнении в меню Tools команды Draw Simulink Diagram (изобразить диаграмму Simulink) можно перейти к моделированию функциональной схемы в среде Simulink (Рис. 3. 18).
Рис. 3. 21. Переход в среду
Simulink
Таким образом, в пункте 3. 1. 3 мы освоили алгоритм построения структурной схемы замкнутой системы управления. Оценили устойчивость системы. Полученные навыки используем для формирования и оптимизации системы управления ТОУ (выбранного варианта объекта управления). Рассмотрим процесс построения и оптимизации системы управления сушилки клинкера (модель сушилки уже получена – это модель sysn4s).
3. 2. Построение структуры системы автоматического регулирования
установки обжига клинкера
Необходимым условием надежной устойчивой работы автоматизированной системы регулирования является правильный выбор типа регулятора и его настроек, гарантирующий требуемое качество регулирования. В зависимости от свойств объектов управления, определяемых его передаточной функцией и параметрами, и предполагаемого вида переходного процесса
выбирается тип и настройка линейных регуляторов.
Согласно исходных данных переходный процесс должен быть апериодическим с малым временем регулирования и малым перерегулированием.
На основании заданных значений передаточных функций датчика, усилительно-преобразовательного устройства, исполнительного механизма (справочные данные) и построенной модели объекта регулирования sysn4s выполним в SIMULINK построение замкнутой системы автоматического регулирования обжига клинкера.
Предварительный вариант системы автоматического регулирования уже получен. Система оптимизирована по характеру переходного процесса и представлена в среде Simulink (Рис. 6. 22). Необходимо скорректировать полученную Simulink-модель системы, включив в нее недостающие элементы: модель датчика, модель усилительно-преобразовательного устройства и модель исполнительного механизма.
Структурно - функциональная блок-схема системы автоматического регулирования представлена на рис. 3. 22.
р
3 22.
Структурно - функциональная блок-схема системы автоматического регулирования
ЗС
– задающий сигнал; Р
– регулятор; УПУ
– усилительно - преобразовательное устройство; ИМ
– исполнительный механизм; ОУ
– объект управления; ДОС
– датчик обратной связи.
В соответствии со структурно-функциональной блок-схемой (Рис. 3. 20) системы автоматического регулирования выполним коррекцию топологии Simulink-модели системы (Рис. 3 21, дополнив ее блоками, имеющими передаточные функции в соответствие со справочными данными: Wдос
= 0.4, Wупу
=15(0.22 + 1); Wим
= 0.19(0.37 + 1) и включим в качестве задающего сигнала единичный скачек (блок Step, Рис. 3 23.)
Выполним команду Start
simulation
в окне модели asd
99*.
В окне Output
осциллографа будем наблюдать переходный процесс (Рис. 3. 25). Регистрация параметров переходной характеристики показывает, что имеющиеся показатели качества не удовлетворяют заданным: Большое время переходного процесса, появилось перерегулирование
Заданные показатели качества и запасы устойчивости:
Время регулирования ≤0.2 с
Статическая ошибка ≤0,05
Перерегулирование ≤1 %
Время нарастания ≤0.1 с
Устойчивость по амплитуде ≥10 дБ
Устойчивость по фазе от 30 до 80 градусов.
3. 23.
Simulink
-
модель системы
Рис.
3. 24
. График переходного процесса
Исходя из выше изложенных рекомендаций и учитывая, что вид переходной характеристики должен соответствовать апериодическому процессу, выполним процедуру оптимизации построенной системы управления. Для оптимизации параметров регулятора воспользуемся пакетом прикладных программ для построения систем управления Nonlinear Control Design (NCD) Blockset, реализующий метод динамической оптимизации.
ГЛАВА 4 ОПТИМИЗАЦИЯ ПАРАМЕТРОВ МОДЕЛИРУЕМОЙ СИСТЕМЫ
Для оптимизации параметров регулятора воспользуемся пакетом прикладных программ для построения систем управления Nonlinear Control Design (NCD) Blockset, который реализует метод динамической оптимизации. Этот инструмент, строго говоря, представляющий собой набор блоков, разработанных для использования с Simulink, автоматически настраивает параметры моделируемых систем, основываясь на определённых пользователем ограничениях на их временные характеристики. Сеанс в среде Simulink с использованием возможностей и блоков NCD Blockset состоит из ряда стадий:
· Создание модели системы из стандартных блоков в среде Simulink.
· Соединение входа блока NCD Outport с теми точками системы, на сигналы которых накладываются ограничения. Этими сигналами могут быть, например выходы системы, их среднеквадратические отклонения и т.д.
· Задание в режиме командной строки MATLAB начальных значений параметров, подлежащих оптимизации.
· Раскрытие блоков двойным щелчком мыши на пиктограмме NCD Outport
· Изменение конфигураций и размеров областей ограничений для сигналов с помощью мыши.
· Задание интервалов дискретизации в меню блока NCD Outport (один или два процента от длительности процесса моделирования) и указание идентификаторов параметров системы, подлежащих оптимизации.
· Задание параметров системы и указание их номинальных значений.
· Сохранение сформированных ограничений в виде файла с помощью команды Save (позднее они могут быть загружены с помощью команды load).
· Процесс оптимизации системы инициализируется нажатием кнопки Start.
Преобразуем Simulink-модель системы, включив в нее дополнительно пропорциональное звено (П-регулятор) с коэффициентом пропорциональности kp
(Рис. 4.1, Gain1). Для этого в окне Function
Parameters
редактора компоненты Gain1 выставим kp
.
(Рис. 4. 2). Подключим к выходу системы блок Signal Constraint из библиотеки Simulink,
Рис.
4
. 1. Преобразованная
Simulink-модель системы управления
Рис. 4. 2. Окно редактора пропорционального звена
содержащийся в разделе Simulink Response Optimization
. В данной операции контролируемым сигналом является реакция системы на единичный скачек, т. е. ее переходная функция. Оптимизируемым параметром является коэффициент kp
.
На переходную функцию накладываются ограничения: максимальное перерегулирование – не более 5%; время нарастания – не более 3 с; длительность переходного процесса - не более 6 с.
Для выполнения процедуры оптимизации наберем в командной строке MATLAB начальное значение настраиваемого параметра kp = 2 и введем его. Далее двойным щелчком мыши откроем рабочее окно блока Signal Constraint (Рис. 4. 3).
Рис.
4 3.
Рабочее окно блока
Signal
Constraint
.
В графической части окна показаны границы контролируемого сигнала, установленные по умолчанию. Для изменения границ в соответствии с заданными значениями используется указатель мыши, позволяющий перемещать линии в вертикальном и горизонтальном направлении. Точную установку линий ограничения можно выполнить, выделяя требуемую линию двойным щелчком левой клавиши мыши. При этом откроется окно редактора Edit
Constraint
(Рис. 4. 4), где можно установить диапазон длины и уровня выделенной линии.
Рис
. 4.4.
Окно
редактора
Edit Constraint.
Следующий этап состоит в объявлении переменных, подлежащих оптимизации. Выбор команды меню Optimization ››Tuned Parameters приведет к открытию диалогового окна задания настраиваемых параметров Tuned Parameters (Рис. 4. 5). В котором, после нажатия кнопки Add
,
появится диалоговое окно Add
Parameters
,
в нижнем поле
Рис.
4
. 5. Окна задания настраиваемых параметров
Tuned
Parameters
которого необходимо набрать имя коэффициента пропорциональности kp
,
подтвердив операцию нажатием кнопки ОК (Рис. 4. 6).
Рис.
4
. 6. Диалоговое окно
Add
Parameters
Появится график переходного процесса, подлежащего корректировке. Теперь необходимо запустить процесс оптимизации, нажав кнопку Start optimization. По окончании процесса оптимизации появится семейство графиков переходного процесса (Рис. 4. 7), в которых отражена динамика оптимизации при различных значениях коэффициента пропорциональности П – регулятора. Совокупность графиков содержит финальный график оптимального переходного процесса. График вписывается в установленные уровни.
Появится также окно выходной информации MATLAB (Рис. 4. 8), где содержится информация о процессе оптимизации и значение kp, соответствующее найденной оптимальной величине параметра П – регулятора. Характер оптимизированного переходного процесса можно также просмотреть на экране осциллографа (Рис. 4. 9).
Рис
. 4. 7.
Диалоговое
окно
Signal Constraint
max Directional First-order
Iter S-count f(x) constraint Step-size derivative optimality Procedure
0 1 0 2538
1 6 0 312.9 2.15 0 1 infeasible
2 9 0 1079 1.02 0 1 infeasible
3 12 0 202.2 0.831 0 1 infeasible
4 15 0 160.7 0.198 0 1 infeasible
5 18 0 203.8 0.199 0 1 Hessian modified; infeasible
6 21 0 168.6 0.203 0 1 Hessian modified twice; infeasible
7 24 0 202.8 0.202 0 1 Hessian modified; infeasible
8 27 0 163.5 0.2 0 1 Hessian modified twice; infeasible
9 30 0 203.5 0.2 0 1 Hessian modified; infeasible
10 33 0 166.7 0.202 0 1 Hessian modified twice; infeasible
11 36 0 203 0.201 0 1 Hessian modified; infeasible
12 39 0 164.6 0.2 0 1 Hessian modified twice; infeasible
13 42 0 203.3 0.201 0 1 Hessian modified; infeasible
14 45 0 166 0.201 0 1 Hessian modified twice; infeasible
15 48 0 203.1 0.201 0 1 Hessian modified; infeasible
16 51 0 165.1 0.201 0 1 Hessian modified twice; infeasible
17 54 0 203.2 0.201 0 1 Hessian modified; infeasible
18 57 0 165.7 0.201 0 1 Hessian modified twice; infeasible
19 60 0 203.2 0.201 0 1 Hessian modified; infeasible
20 63 0 165.3 0.201 0 1 Hessian modified twice; infeasible
21 66 0 203.2 0.201 0 1 Hessian modified; infeasible
22 69 0 165.5 0.201 0 1 Hessian modified twice; infeasible
23 72 0 203.2 0.201 0 1 Hessian modified; infeasible
24 75 0 165.4 0.201 0 1 Hessian modified twice; infeasible
25 78 0 203.2 0.201 0 1 Hessian modified; infeasible
26 81 0 165.5 0.201 0 1 Hessian modified twice; infeasible
27 84 0 203.2 0.201 0 1 Hessian modified; infeasible
28 87 0 165.4 0.201 0 1 Hessian modified twice; infeasible
29 90 0 203.2 0.201 0 1 Hessian modified; infeasible
30 93 0 165.5 0.201 0 1 Hessian modified twice; infeasible
31 96 0 203.2 0.201 0 1 Hessian modified; infeasible
32 99 0 165.4 0.201 0 1 Hessian modified twice; infeasible
33 102 0 203.2 0.201 0 1 Hessian modified; infeasible
34 105 0 165.5 0.201 0 1 Hessian modified twice; infeasible
35 108 0 203.2 0.201 0 1 Hessian modified; infeasible
36 111 0 165.5 0.201 0 1 Hessian modified twice; infeasible
37 114 0 203.2 0.201 0 1 Hessian modified; infeasible
38 117 0 165.5 0.201 0 1 Hessian modified twice; infeasible
39 120 0 203.2 0.201 0 1 Hessian modified; infeasible
40 123 0 165.5 0.201 0 1 Hessian modified twice; infeasible
41 126 0 203.2 0.201 0 1 Hessian modified; infeasible
42 129 0 165.5 0.201 0 1 Hessian modified twice; infeasible
43 132 0 203.2 0.201 0 1 Hessian modified; infeasible
44 135 0 165.5 0.201 0 1 Hessian modified twice; infeasible
45 138 0 203.2 0.201 0 1 Hessian modified; infeasible
46 141 0 165.5 0.201 0 1 Hessian modified twice; infeasible
47 144 0 203.2 0.201 0 1 Hessian modified; infeasible
48 147 0 165.5 0.201 0 1 Hessian modified twice; infeasible
49 150 0 203.2 0.201 0 1 Hessian modified; infeasible
50 153 0 165.5 0.201 0 1 Hessian modified twice; infeasible
51 156 0 203.2 0.201 0 1 Hessian modified; infeasible
52 159 0 165.5 0.201 0 1 Hessian modified twice; infeasible
53 162 0 203.2 0.201 0 1 Hessian modified; infeasible
54 165 0 165.5 0.201 0 1 Hessian modified twice; infeasible
55 168 0 203.2 0.201 0 1 Hessian modified; infeasible
56 171 0 165.5 0.201 0 1 Hessian modified twice; infeasible
57 174 0 203.2 0.201 0 1 Hessian modified; infeasible
58 177 0 165.5 0.201 0 1 Hessian modified twice; infeasible
59 180 0 203.2 0.201 0 1 Hessian modified; infeasible
60 183 0 165.5 0.201 0 1 Hessian modified twice; infeasible
61 186 0 203.2 0.201 0 1 Hessian modified; infeasible
62 189 0 165.5 0.201 0 1 Hessian modified twice; infeasible
63 192 0 203.2 0.201 0 1 Hessian modified; infeasible
64 195 0 165.5 0.201 0 1 Hessian modified twice; infeasible
65 198 0 203.2 0.201 0 1 Hessian modified; infeasible
66 201 0 165.5 0.201 0 1 Hessian modified twice; infeasible
67 204 0 203.2 0.201 0 1 Hessian modified; infeasible
68 207 0 165.5 0.201 0 1 Hessian modified twice; infeasible
69 210 0 203.2 0.201 0 1 Hessian modified; infeasible
70 213 0 165.5 0.201 0 1 Hessian modified twice; infeasible
71 216 0 203.2 0.201 0 1 Hessian modified; infeasible
72 219 0 165.5 0.201 0 1 Hessian modified twice; infeasible
73 222 0 203.2 0.201 0 1 Hessian modified; infeasible
74 225 0 165.5 0.201 0 1 Hessian modified twice; infeasible
75 228 0 203.2 0.201 0 1 Hessian modified; infeasible
76 231 0 165.5 0.201 0 1 Hessian modified twice; infeasible
77 234 0 203.2 0.201 0 1 Hessian modified; infeasible
78 237 0 165.5 0.201 0 1 Hessian modified twice; infeasible
79 240 0 203.2 0.201 0 1 Hessian modified; infeasible
80 243 0 165.5 0.201 0 1 Hessian modified twice; infeasible
81 246 0 203.2 0.201 0 1 Hessian modified; infeasible
82 249 0 165.5 0.201 0 1 Hessian modified twice; infeasible
83 252 0 203.2 0.201 0 1 Hessian modified; infeasible
84 255 0 165.5 0.201 0 1 Hessian modified twice; infeasible
85 258 0 203.2 0.201 0 1 Hessian modified; infeasible
86 261 0 165.5 0.201 0 1 Hessian modified twice; infeasible
87 264 0 203.2 0.201 0 1 Hessian modified; infeasible
88 267 0 165.5 0.201 0 1 Hessian modified twice; infeasible
89 270 0 203.2 0.201 0 1 Hessian modified; infeasible
90 273 0 165.5 0.201 0 1 Hessian modified twice; infeasible
91 276 0 203.2 0.201 0 1 Hessian modified; infeasible
92 279 0 165.5 0.201 0 1 Hessian modified twice; infeasible
93 282 0 203.2 0.201 0 1 Hessian modified; infeasible
94 285 0 165.5 0.201 0 1 Hessian modified twice; infeasible
95 288 0 203.2 0.201 0 1 Hessian modified; infeasible
96 291 0 165.5 0.201 0 1 Hessian modified twice; infeasible
97 294 0 203.2 0.201 0 1 Hessian modified; infeasible
98 297 0 165.5 0.201 0 1 Hessian modified twice; infeasible
99 300 0 203.2 0.201 0 1 Hessian modified; infeasible
100 303 0 165.5 0.201 0 1 Hessian modified twice; infeasible
Maximum number of iterations exceeded.
Restart or go to Optimization Options to increase the maximum of iterations.
kp =
0.0437
Рис.
4. 8.
окно выходной информации
Рис.
4. 9
. Характер оптимизированного переходного процесса
Исходя из приоритета характеристик переходного процесса в нашем случае наилучший будет: kp = 0,0437.
ГЛАВА 5. АНАЛИЗ КАЧЕСТВА СИСТЕМЫ УПРАВЛЕНИЯ
Необходимо выполнить анализ построенной системы управления и дать оценку ее качества по основным показателям. Анализ снятой переходной характеристики системы после выполнения оптимизации показывает, что новые показатели качества переходного процесса:
Время регулирования составляет 6 с.
Установившееся значение – 0,18
Время нарастания – 3
Статическая ошибка – 0
Перерегулирование - 0 %
удовлетворяют заданным показателям.
Заключение
В данном курсовом проекте проведена идентификация объекта автоматического регулирования.
Проведена проверка на наблюдаемость и управляемость объекта управления.
На основе анализа переходных характеристик объекта управления был выбран наиболее подходящий для данного переходного процесса П – регулятор.
Проведена оптимизация настроечных параметров этого регулятора.
В результате введения в систему П - регулятора были получены следующие параметры системы:
- Время переходного процесса 11 с.;
- Время нарастания – 10 с.
- Перерегулирование – 0%;
- Статическая ошибка – нет;
- Запас по фазе – 70 градусов;
Учитывая полученные значения и принятые допущения параметров системы можно утверждать, что выполнены все поставленные в задании на курсовую работу требовани
.