Главная Учебники - Разные Лекции (разные) - часть 35
|
Государственное образовательное учреждение высшего профессионального образования
“
Камчатский государственный технический университет
”
Факультет (филиал) ФИТ специальность (направление) УИ Кафедра Систем Управления Дисциплина: Автоматизация и проектирование средств и систем управления ПОЯСНИТЕЛЬНАЯ ЗАПИСКА
к курсовому проекту (работе) на тему:
«Автоматизированное проектирование системы управления технологическим процессом производства цемента
»
Студент: Рунго Олег Викторович Группа: 05-УИ шифр_____________________________________ Обозначение курсового проекта (работы)___________________________ Проект (работа) защищен(а) на оценку_____________________________ Руководитель проекта (работы)_____________________ Пюкке Г.А. подпись, дата инициалы и фамилия Члены комиссии________________________________________________ подпись, дата инициалы и фамилия _______________________________________________________________ подпись, дата инициалы и фамилия _________________________________________________________________ подпись, дата инициалы и фамилия ПЕТРОПАВЛОВСК-КАМЧАТСКИЙ 2009 г. Федеральное агентство по рыболовству
Государственное образовательное учреждение высшего профессионального образования
“
Камчатский государственный технический университет
”
Факультет (филиал) ФИТ специальность (направление) УИ Дисциплина: Автоматизация и проектирование средств и систем управления ЗАДАНИЕ НА КУРСОВОЙ ПРОЕКТ (РАБОТУ)
Студент шифр_______________группа 05-УИ 1 Тема «
Автоматизированное проектирование системы управления технологическим процессом производства цемента
»
2 Срок представления проекта (работы) к защите________________200_ г. 3 Исходные данные для разработки: Передаточная функция компоненты объекта автоматизации - 4 Содержание пояснительной записки: Титульный лист Задание Содержание Введение 1 Описание технологического процесса 2 Идентификация объекта автоматизации 3 Построение системы управления технологическим процессом 4 Оптимизация параметров моделируемой системы 5 Анализ качества системы управления Заключение Список использованных источников 5 Перечень графического материала: Отсутствует Руководитель проекта (работы)__________________ Пюкке Г.А. подпись, дата инициалы и фамилия Задание принял к исполнению___________________ подпись, дата инициалы и фамилия
1 Описание технологического процесса
1.1 Обоснование целесообразности и необходимости автоматизации технологического процесса…………………………………………………………………………..7 1.2 Описание технологического процесса и производственного оборудования…………………………………………………………………………………………………..8 1.3 Требования к системе автоматизации технологического процесса…………………………………………………………………………………………………………..12 2 Идентификация объекта автоматизации
2.1 Особенности построения моделей технологических объектов управления……………………………………………………………………………………………………..13 2.2 Виды моделей линейных стационарных динамических объектов………………………………………………………………………………………………………..16 2.3 Виды моделей пакета System Identification Toolbox…………………………………………………………………………………………………………..20 2.4 Основные операторы и функции пакета System Identification Toolbox…………………………………………………………………………………………………………..23 2.5 Пример использования пакета System Identification Toolbox для идентификации технологических объектов управления………………………………………………………………………………………….………….25 2.6 Обработка данных при построении модели объекта управления………………………………………………………………………………………………..…….29 2.7 Оценивание статистических и частотных характеристик исходных данных……………………………………………………………………………………………………………..33 2.8 Параметрическое оценивание данных……………………………………….……38 2.9 Функции преобразование моделей……………………………………………………42 2.10 Проверка адекватности модели………………………………………….………….49 2.11 Анализ модели технического объекта управления……………….…..56 2.12 Основные результаты идентификации технического объекта идентификации……………………………………………………………………………………………..71 3.1 Задание структуры системы автоматического управления, проверка системы управления на устойчивость…………………………………………………..74 3.2 Построение структуры системы автоматического регулирования установки обжига клинкера………………………………………………………………………76 4 Оптимизация параметров моделируемой системы
………………….80 5 Анализ качества системы управления
…………………………………….….85 Характерной особенностью современного этапа автоматизации производства состоит в том, что он опирается на революцию в вычислительной технике, на самое широкое использование микропроцессоров и контроллеров, а также на быстрое развитие робототехники, гибких производственных систем, интегрированных систем проектирования и управления, SCADA-систем разработки программного обеспечения. Целью автоматизации является снижение объёма ручного труда, обеспечение стабильности характеристик технологического процесса, обеспечение возможности наблюдения, анализа и управления параметрами технологического процесса человеком. Результатом этого процесса является получение автоматизированной системы. Автоматизация производства позволяет повысить качество и снизить себестоимость продукции. Автоматизированная система - это совокупность управляемого объекта и автоматизированных управляющих устройств, в которой часть функций управления выполняет человек. Автоматизированная система получает информацию от объекта управления, передаёт, преобразует и обрабатывает её, формирует управляющие команды и выполняет их на управляемом объекте. Человек определяет цели и критерии управления, корректирует их, если изменяются условия. Применение средств и систем автоматизации позволяет решать следующие задачи: · вести процесс с производительностью, максимально достижимой для данных производительных сил, автоматически учитывая непрерывные изменения технологических параметров, свойств исходных материалов, изменений в окружающей среде, ошибки операторов; · · автоматически управлять процессами в условиях вредных или опасных для человека. Решение поставленных задач предусматривает целый комплекс вопросов по проектированию и модернизации существующих и вновь разрабатываемых систем автоматизации технологических процессов и производств. Использование автоматизированной системы управления повышает надежность работы устройств и улучшает технико-экономические показатели за счет следующих усовершенствований: · реализации системы на базе средств, отвечающих современному уровню развития техники управления технологическими процессами; · реализации более сложных законов автоматического регулирования, точнее и полнее учитывающих специфику протекающих технологических процессов; · внедрения системы безопасного управления объектами повышенной опасности (розжиг, плановый и аварийный останов, опрессовка и т.д.); · создания комфортных условий работы для оперативного персонала, снижающих нагрузку оператора, облегчающих принятие им решений по управлению. Приближения решений к оптимальным, благодаря лучшему информационному обеспечению (представление данных в требуемом объеме, в удобном для восприятия виде в реальном времени); · повышения меры ответственности персонала за счет наличия в системе функций слежения и протоколирования действий персонала по управлению системой; · · выдачи технико-экономических показателей и объективной информации о технологическом процессе, которая может быть использована неоперативным инженерно-техническим и административным персоналом для решения производственных и организационно-экономических задач. В данном учебном пособии на конкретном примере одного из видов технологического процесса производства рассматривается методика анализа и синтеза системы автоматизации. Изложение материала базируется на использование возможностей современной интегрированной системы компьютерной математики MATLAB и её приложений. Рассмотренные в учебном пособии вопросы должны найти отражение в курсовом и дипломном проектировании по автоматизации технологических процессов и производств. Целью курсового проектирования по дисциплине "Автоматизация проектирования систем и средств управления" является закрепление знаний, выработка навыков проектирования систем с использованием элементов автоматизации проектных процедур, работы с технической литературой и данными Интернета: государственными и отраслевыми стандартами, каталогами заводов-изготовителей, справочной литературой, базами данных сайтов заводов-изготовителей.
1. 1. Обоснование целесообразности и необходимости автоматизации технологического процесса
В этом разделе приводится информация об области применения производимой продукции, а также информация об ее назначении (например: керамические изделия). Перечисляются этапы технологического процесса (например: производство керамических изделий состоит из нескольких этапов): - процесс приготовления шихты; - сушка керамического порошка; - формовка и прессование керамических изделий; - обжиг керамических изделий). Описываются методы изготовления продукта и исходные материалы производства (например: пластичное формование керамических изделий, или другой метод, который применяют для формования изделий сложной формы, - метод шликерного литья). Исходные материалы (например: глинистые и тонкомолотые материалы, каолин, глины, отощающие компоненты и плавни). Перечисляются контролируемые параметры и допустимые пределы отклонения значений параметров (например: влажность массы для пластического формования должна быть в пределах 18-25%; влажность литейного шликера - в пределах 31-35%; отклонение влажности пластической массы от заданной средней величины не должна превышать ± 0,5%, шликера - соответственно ± 0,8%). Делается вывод о необходимости применения автоматизированной системы контроля и управления технологическим процессом (по показателям экономичности, точности функционирования, быстродействия, инерционности, безопасности и др). Выбирается этап технологического процесса производства, подлежащий автоматизации, обеспечивающей устойчивую работу технологического оборудования и осуществляющей управляющие воздействия для компенсации изменений в технологическом процессе (например: автоматизация процесса сушки исходного материала). Контроль влажности изделий позволяет корректировать режим сушки и поддерживать влажность керамической массы в заданных пределах). 1. 2. Описание технологического процесса и производственного оборудования
Рассматриваются различные современные устройства, используемые для реализации выбранного процесса производства. Приводится их структура и описание этапов функционирования. Приводится мнемоническая схема автоматического регулирования процесса производства (например: для рассматриваемого примера сушки исходного материала используются распылительные сушилки). Распылительные сушилки применяют для снижения влажности массы до 7- 9% перед ее прессованием. Математическое описание звеньев системы автоматизации следует начинать с ТОУ. В технической литературе тепловые объекты автоматизации (например, распылительная сушилка) с достаточной степенью точности описываются последовательным соединением звена чистого запаздывания и апериодического звена первого порядка. Значения постоянных времени и времени запаздывания определяются по переходным характеристикам. Однако в ряде случаев, когда невозможно получить переходную характеристику при составлении математической модели ТОУ следует использовать статистические данные по их характеристикам, полученные экспериментально в ходе штатной работы установки методом пассивного эксперимента, когда через определенные промежутки времени фиксируются значения входной и выходной величины ТОУ. Такой путь называется идентификацией объектов автоматизации. 1. 3. Требования к системе автоматизации технологического процесса
Анализ технологического процесса позволяет построить структуру системы автоматизации и сформулировать требования, предъявляемые к системе автоматизации технологического процесса. В приведенном выше примере применение автоматического регулирования влажности шликера по температуре отходящих газов позволяет: - уменьшить среднеквадратическое отклонение влажности шликера; - увеличить качество керамических изделий; - уменьшить брак при прессовании. Для обеспечения положительного эффекта использования системы автоматизации, к ней предъявляются следующие требования: - статическая ошибка: не более ± 5 %; - перерегулирование: не более 10 %; - время переходного процесса: от 0, 1 до 0, 2 с; - запас устойчивости по амплитуде: не менее 20 дБ; - запас устойчивости по фазе: от 20 до 80 градусов. ГЛАВА 2. ИДЕНТИФИКАЦИЯ ОБЪЕКТА АВТОМАТИЗАЦИИ
2. 1. Особенности построения моделей технологических объектов управления
Сложность идентификации технологических процессов во многом зависит от наличия априорной информации о технологических объектах управления, их статических и динамических характеристик. Определение характеристик объекта управления выполняется различными способами, например, могут быть рассмотрены методы, связанные с проведением физичес В самом общем случае, связь между входным и «теоретическим» выходным сигналами может быть задана в виде некоторого оператора Ψ. При этом наблюдаемый выходной сигнал объекта может быть описан на основе соотношения: y
(
t
) =
Ψ[u
(
t
)
] + e
(
t
).
Принцип суперпозиции позволяет объединить все действующие помехи в одну общую e
(
t
)
и приложить ее к выходу линейной модели. При рассмотрении задач идентификации все помехи считают статически независимыми, что позволяет моделировать их в виде гауссовского процесса (шума). Перед началом экспериментальных исследований проводят априорный анализ перечня входных переменных с целью отбора и включения в состав модели информативных параметров, т. е. оказывающих наиболее сильное воздействие на выходные переменные y
(
t
)
. В первую очередь в их состав включают управляющие входные переменные, с помощью которых осуществляется регулирующее воздействие на ТОУ. Если в процессе идентификации структура модели не меняется, то выполняется только оценивание параметров модели (идентификация в узком смысле). Однако можно менять и структуру модели, подбирая наиболее адекватную описываемому процессу. При этом вид модели, ее структура выводится из физических представлений о сути процессов в ТОУ. Например, простейший сглаживающий фильтр (RC-цепь) описывается известными законами электротехники, для него можно записать: u(t) = RCdy(t)/dt + y(t),
где Uin
(t) = u(t), Uout
(t) = y(t).
Если такая структура (с точностью до вектора коэффициентов β
) известна, то при известном входном сигнале u
(
t
)
описание объекта можно представить в виде: y
(
t
) =
F
(
β
,
t
) +
e
(
t
),
где F
–
функция известного вида, зависящая от β
и времени t
.
Последнее уравнение позволяет после проведения эксперимента, заключающегося в фиксации входного и выходного сигналов на каком-то интервале времени, провести обработку экспериментальных данных и каким-либо методом (например, методом наименьших квадратов) найти оценку вектора параметров β
. Отметим, что при экспериментальном определении параметров модели необходимо обеспечить: ● подбор адекватной структуры модели; ● выбор такого входного сигнала, чтобы по результатам эксперимента можно было найти оценки всех параметров модели. Наиболее просто задача определения параметров решается для линейных объектов, для которых выполняется принцип суперпозиции. В задачах идентификации под линейными объектами чаще понимаются объекты, линейные по входному воздействию. Как правило, идентификация – многоэтапная процедура, состоящая из этапов: 1. Структурная идентификация, включающая определение структуры математической модели на основании теоретических соображений. 2. Параметрическая идентификация включает в себя процедуру оценивания параметров модели по экспериментальным данным. 3. Проверка адекватности – проверка качества модели в смысле выбранного критерия близости выходов модели и объекта. Следует отметить, что в связи с многообразием объектов и различных подходов к их моделированию существует множество вариантов решения задачи идентификации.
Линейные непрерывные стационарные динамические объекты могут быть представлены (без учета действия шума e
(
t
)
) в виде: Дифференциального уравнения
. Наиболее универсальная модель, имеющая форму
где na
– порядок модели (na
> nb
); ai
и bj
– постоянные коэффициенты (параметры модели); u
(
j
)
(
t
)
и y
(
i
)
(
t
)
– производные, соответственно, входного и выходного сигналов. Передаточной функции.
Модель определяется как отношение преобразования Лапласа выходного сигнала к преобразованию Лапласа входного сигнала где L
{●}
– символ преобразования Лапласа, р –
переменная (оператор Лапласа). Импульсной характеристики
w
(
t
) и
переходной функции
h
(
t
)
. Импульсная характеристика определяется как реакция объекта на входной сигнал в виде δ-функции. Переходная функция h(t) определяется как реакция объекта на входной сигнал в виде единичного скачка. Соотношения между этими характеристиками имеют следующий вид: L
{
w
(
t
)}=
W
(
p
),
w
(
t
)=
h
’
(
t
) ,
L
{
h
(
t
)}=
W
(
p
)/
p
При нулевых начальных условиях связь между выходными и входными сигналами описывается интегралом свертки: или в операторной форме: Частотной характеристики.
Частотные характеристики объекта определяются его комплексным коэффициентом передачи W
(
jω
)
. Модуль комплексного коэффициента передачи │W
(
jω
)
│= A
(
ω
)
представляет собой амплитудно-частотную характеристику (АЧХ) объекта с передаточной функцией W
(
p
)
, а аргумент arg(W
(
jω
))=
φ
(
ω
)
– фазочастотную характеристику (ФЧХ). Графическое представление W
(
jω
)
, на комплексной плоскости при изменении ω
от 0 до ∞, то есть график амплитудно-фазовой характеристики (АФХ) в полярных координатах в отечественной литературе называется годографом, а в англоязычной – диаграммой Найквиста. В теории автоматического управления часто используется логарифмическая амплитудно-частотная характеристика (ЛАЧХ), равная 20 lg │W
(
jω
)
│. В 70-е годы 20 века Розенброком был создан метод «размытых» частотных характеристик, предназначенный для автоматизированного проектирования систем с несколькими входами и выходами, ориентированный на использовании средств вычислительной техники и названный в последствие методом переменных состояния (МПС). В основе этого метода лежит представление дифференциальных уравнений в нормальной форме Коши, которое дополняется алгебраическими уравнениями, связывающими выходные переменные с переменными состояния: x
’ =
Ax
+
Bu
y
=
Cx
, где u
– вектор входных воздействий; y
– вектор выходных воздействий; x
– вектор переменных состояния; A, B, C
– матрицы коэффициентов размерности n
x n
,
n
x m
,
r
x n
соответственно; n
– число переменных состояния или максимальная степень производной исходного дифференциального уравнения; m
– число входов; r
– число выходов. Для систем с одним входом и одним выходом уравнения переменных состояния можно сформулировать следующим образом. При выборе n
координат системы (объекта) в качестве переменных ее состояния (такими координатами, например, могут быть выходной сигнал y
(
t
)
и n
-1
его производных) принимаем xi
,
i
=
1,2,…, n
и данную систему можно описать следующими уравнениями для переменных состояния: х
′
(t
) =
A
х
(
t
)
+ Bu
(
t
),
y
(
t
) =
C
х
(
t
) +
Du
(
t
),
где х
(
t
) =
[x
1
(
t
),
x
2
(
t
),…,
xn
(
t
)
]t
– вектор-столбец переменных состояния; A
,
B
,
C
,
и D
при скалярных u
(
t
)
и y
(
t
)
– соответственно матрица размера n
Для дискретных объектов, функционирование которых представляется дискретным временем tk
=
kT
(T
– интервал дискретизации), наиболее общим видом описания является разностное уравнение (аналог дифференциального): yk
+a1
yk-1
+ ... +ana
yk – na
= b1
uk
+ b2
uk – 1
+ b3
uk - 2
+ ... + bnb
uk – nb + 1
,
где yk – i
= y
[(k – i
)T
] , uk – j
= u
[(k – j
)T
] Связь между входом и выходом может быть отражена следующими соотношениями: • через дискретную свертку: где ω –
ординаты решетчатой весовой функции объекта, или, с использованием аппарата Z
–
преобразования: • через дискретную передаточную функцию: которая определяется на основании разностного уравнения после применения к обеим частям этого уравнения Z
– преобразования: На практике в большинстве случаев измерение непрерывных сигналов производится в дискретные моменты времени, что представляет определенное удобство при последующей обработке данных на ЭВМ. Поэтому представление непрерывных объектов дискретными моделями является актуальной задачей. Хотя такое представление может быть осуществлено с некоторой степенью приближенности. 2. 3.
Виды
моделей
пакета
System Identification Toolbox
Одним из расширений MATLAB является пакет System Identification Toolbox, который содержит средства для создания математических моделей линейных динамических систем, на основе наблюдаемых входных и выходных данных. Он имеет удобный графический интерфейс, позволяющий организовывать данные и создавать модели. Приведем несколько распространенных моделей дискретных объектов, используемых в пакете System Identification Toolbox для временной области, учитывающих действие шума наблюдения: · Модель авторегрессии AR (AutoRegressive) – считается простым описанием: A(z) y(t) = e(t)
, где A(z) =
1 + a
1
z
– 1
+ a
2
z
– 2
+...+ a na
z – na
. · ARX – модель (Autoregressive with eXternal input) – более сложная модель: A
(z
)y(t
) = B
(z
) u
(t
) + e
(t
) , Здесь и ниже e
(t
) – дискретный белый шум. · ARMAX где nk
– величина задержки (запаздывания), · Модель «вход-выход» (в иностранной литературе такая модель называется «Output-Error», то есть «выход-ошибка», сокращенно ОЕ): где · Модель Бокса-Дженкинса (BJ): полиномы B
(z
), F
(z
), C
(z
) определены ранее, а полином D
(z
) определяется по формуле: Данные модели можно рассматривать, как частные случаи обобщенной параметрической линейной структуры: при этом все они допускают расширение для многомерных объектов (имеющих несколько входов и выходов). · Модель для переменных состояния (State-space): где A
,
B
,
C
,
D
– матрицы соответствующих размеров, v
(t
) – коррелированный белый шум наблюдений. Возможна и другая (так называемая о где К – некоторая матрица (вектор столбец), е(t) – дискретный белый шум (скаляр). В своей работе пакет System Identification Toolbox использует три внутренних вида матричного представления моделей, которые с помощью операторов и функций пакета преобразуются во все выше перечисленные виды моделей объектов: ● так называемый тета-формат (для временных моделей); ● частотный формат (для частотных моделей); ● формат нулей и полюсов. 2. 4. Основные операторы и функции пакета
System
Identification
Toolbox
Приведем основные операторы и функции пакета System Identification Toolbox, которые набираются в командной строке MATLAB или могут быть использованы при написании программ для m-файлов. Наиболее полную информацию о содержании, написании и использовании этих функций можно получить в литературе [2] и справочной части help MATLAB. idhelp
– используется для вызова подсказаки о возможностях пакета. iddemo
– используется для вызова демонстрационных примеров. ident
– команда вызова графического интерфейса пользователя. m
idprefs
– команда задает (изменяет) директорию для файла midprefs.mat,хранящего информацию о начальных параметрах графического интерфейса мользователя при его открытии. predict
– команда осуществляет прогноз выхода объекта по его ттета-модели и сучетом информации о его предыдущих фактических значениях выхода ( рекомендуется для расчета прогноза значений временной последовательности). id
sim
– возвращает выход модели тета-формата. iddata
– создает файл объекта данных. detrend –
удаляет тренд из набора данных. idfilt
– команда фильтрует данные с помощью фильтра Баттерворта. idinput
- команда генерирует входные сигналы для идентификации. merge
– объединяет несколько экспериментов. misdata
– оценивает и заменяет потерю входных и выходных данных в файле созданном с помощью команды iddata
. esample
– восстанавливает форму квантованного сигнала данных прореживанием и интерполяцией и изменяет частоту дискретизации. Функции непараметрического оценивания
: covf
– выполняет расчет авто - и взаимных корреляционных функций совокупности экспериментальных данных. cra –
определяет оценку импульсной характеристики методом коррелированного анализа для одномерного (один вход – один выход) объекта. etfe
– возвращает оценку дискретной передаточной функции для обобщенной линейной модели одномерного объекта в частотной форме. impulse
– выводит на дисплей импульсную характеристику модели. spa
– возвращает частотные характеристики объекта и оценки спектральных плотностей его сигналов для обобщенной линейной модели объекта (возвращает модель объекта в частотном формате). step
– выводит на дисплей переходную (временную) характеристику модели объекта (реакция на единичное ступенчатое воздействие). Функции параметрического оценивания
: ar
– оценивает параметры модели авторегрессии (AR), то есть коэффициенты полинома A
(z),
при моделировании скалярных временных последовательностей. arx
– оценивает параметры моделей ARX и AR. bj
– оценивает параметры модели Бокса-Дженкинса. Ivar
– оценивает параметры скалярной AR-модели. iv4
– оценивает параметры для моделей ARX с использованием четырехступенчатого метода инструментальной переменной. n4sid
– используется для оценивания параметров моделей для переменных состояния в канонической форме при произвольном числе входов и выходов. ivx
– оценивает параметры ARX-моделей методом инструментальной переменной. oe
– оценивает параметры ОЕ-модели. pem
– оценивает параметры обобщенной многомерной линейной модели. Функции задания структуры модели
: idpoly
– создавать модель объекта в виде полинома. idss
– создает модель объекта в виде переменных состояния. i
darx
– создает многопараметрическую ARX-модель объекта. idgrey
– создает модель объекта, созданную пользователем. arx
2
th
– создает матрицу модели тета-формата по полиномам ARX-модели многомерного объекта. canform
– создает каноническую форму модели для переменных состояния многомерного объекта. mf
2
th
– преобразует структуру модели для пременных состояния в тета-формат. poly
2
th
– создает модель тета-формата из исходной модели «вход-выход». Функции извлечения данных о моделях
: arxdata
– возвращает матрицы коэффициентов полиномов ARX-моделей, а также их среднеквадратические отклонения. polydata
– возвращает матрица коэффициентов полиномов. ssdata
– функция возвращает матрицы(и величину интервала дискретизации в дискретном случае) ss-моделей (моделей переменных состояния). tfdata
– данная функция возвращает числитель и знаменатель передаточной функции. zpkdata
– функция возвращает нули, полюсы и обобщенные коэффициенты передачи для каждого канала модели тета-формата или LTI-модели (если используется пакет Control System Toolbox c именем sys). idfrd
– данная функция создает частотную модель объекта в frd-формате. idmodred
– функция уменьшает порядок модели объекта. c2d, d2c
– первая функция преобразует непрерывную модель в дискретную. Вторая – наоборот. ss, tf, zpk, frd
– функции создания моделей стационарных систем в виде модели переменных состояния (ss), передаточной функции по ее заданным нулям и полюсам (zpk), передаточной функции, записанной в операторном виде (tf) и в частотном виде (frd). Функции отображения модели
: bode
,
bodeplot
,
ffplot
– функции отображения логарифмических частотных характеристик. plot
– отображение входных - выходных данных для данных объекта. present
– функция отображения вида модели тета-формата с оценкой среднеквадратического отклонения, функции потерь и оценки точности модели.. pzmap
– отображает нули и полюсы модели (с областями неопределенности). nyquist
– отображает диаграмму Найквиста (гадограф АФХ) передаточной функции. view
– отбражение LTI-модели (при использовании пакета Control System Toolbox). compare
– функция позволяет сравнить выходы модели и объекта с выводом на дисплей сравнительных графиков и указанием оценки адекватности модели. resid
– функция вычисляет остаточную ошибку для заданной модели и соответствующие корреляционные функции. Функции выбора структуры модели
: aic, fpe
– функции вычисляют информационный критерий AIC и конечную ошибку модели. arxstruc
– функция вычисляет функции потерь для ряда различных конкурирующих ARX-моделей с одним выходом. selstruc
– функция осуществляет выбор наилучшей структуры модели В состав библиотеки System ID Blocks блоков Simulink входят блоки, позволяющие производить оценивание ряда типовых моделей: ● модели авторегрессии AR (AutoRegressive model estimator); ● ARX-модели (AutoRegressive Moving Average with eXternal input model estimator); ● ARX-модели (AutoRegressive Moving Average with eXternal input model estimator); ● модели Бокса-Дженкинса BJ (Box-Jenkins model estimatjr); ● обобщенной линейной модели (General model estimator using Predictive Error Method); ● модели «вход-выход» OE (Output-error model estimator). Правила работы с данными блоками аналогичны правилам для других блоков Simulink. Полученная модель отображается в основном окне MATLAB. 2. 5. Пример использования пакета
System
Identification
Toolbox
для идентификации технологических объектов управления
· априорный анализ ТОУ с целью выбора структуры модели; · проведение предварительного (небольшого по объему) исследования объекта с целью уточнения оценки структуры модели (этот этап желателен, особенно при отсутствии априорной информации о ТОУ); · разработка методики основного экспериментального исследования ТОУ, составление плана эксперимента; · проведение основного экспериментального исследования для получения массива данных (ui
, yi
); · математическая обработка массива данных (с использованием пакета System Identification Toolbox) с целью определения параметров модели и ее адекватности, доверительных границ параметров и выходной координаты модели. При этом в процессе исследования ТОУ необходимо принять некоторые допущения, позволяющие применить хорошо отработанный аппарат анализа стационарных, линейных объектов: · технологический объект управления является системой с сосредоточенными параметрами; · технологический объект управления стационарен, т. е. статические и динамические свойства ТОУ неизменны во времени; · уравнения моделей ТОУ линеаризуются в малом, т.е. при небольших отклонениях ± ∆
yi
от выбранной "рабочей" точки (рабочего режима) ТОУ. Массив данных [u
i
, y
i
] образуется в результате трудоемкой операции расшифровки регистрограмм по приборам измерительной системы. Однако широкое развитие микропроцессорной и вычислительной техники и внедрение ее в производственные технологические процессы позволили существенно усовершенствовать техническое обеспечение идентификации ТОУ. Обработка массива данных с помощью пакета System Identification Toolbox предполагает следующие этапы: · обработка и преобразование данных с целью создания файла данных; · анализ экспериментальных данных с целью предварительного определения основных характеристик ТОУ; · параметрическое оценивание данных с целью создания различных видов моделей (описанных во втором разделе) в тета-формате; · задание структуры модели; · изменение и уточнение структуры модели (если это необходимо); · проверка адекватности и сравнение различных моделей с целью выбора наилучшей; · преобразование модели из тета-формата в вид удобный для дальнейшего использования при анализе и синтезе системы управления. На каждом этапе идентификации имеется возможность графического отображения результатов моделирования и извлечения необходимой информации об объекте. 2. 6. Обработка данных при построении модели объекта управления
Начнем процедуру построения аналитической модели технического объекта управления (ТОУ) – сушилки. В сушилку подводится шликер, где он распыляется. Инжекционные горелки создают высокую температуру в зоне распыления материала. Распыленные частицы, теряя влагу, уже в виде порошка собираются в днище сушилки, откуда поступают непосредственно с бункера над прессами. Основная задача системы регулирования состоит в поддержании характеристик продукта (влажность шликера) в заданных пределах. При этом, необходимость организации процесса автоматического регулирования обусловлена и наличием также случайных возмущений (т. е. неотвратимых в реальной обстановке факторов), воздействие которых может привести к нарушению технологического процесса и отклонению характеристик продукта от заданных значений. В частности, к таким возмущениям могут относиться случайные изменения интенсивности Как было отмечено выше, аналитическая модель сложных ТОУ может быть построена на основе массива входных и выходных данных, полученных в результате физического эксперимента, проведенного над ТОУ. При этом необходимо учитывать и воздействие случайных факторов, поэтому при проведении физического эксперимента случайное воздействие должно быть смоделировано. Для дальнейших вычислений будем использовать статистические данные, полученные при изучении теплового объекта, и содержащиеся в файле project
4,
который включает массив данных, состоящий из 159 значений входного параметра – расход газа, выраженный в м3
/час и 159 значений выходного параметра – температуры в объекте в градусах Цельсия. Для проведения модельного эксперимента и формирования массива входных и выходных данных об объекте автоматизации следует в Simulink построить S-модель установки для снятия массивов входных и выходных данных, изображенную на рис. 2. 0. Рис.
2.0. S-модель объекта автоматизации
Набрать в блоке Transfer Fcn передаточную функцию составляющей компоненты объекта автоматизации своего варианта и запустить моделирование. Перед запуском моделирования в раскрывающемся меню Simulation следует открыть окно Configuration Parameters и установить Stop time равным 999 (или набрать 999 непосредственно в окне модели). Для загрузки файла в рабочую среду Workspace системы MATLAB необходимо в режиме командной строки выполнить следующую команду: В результате выполнения команды в рабочей области Workspace появится массив входных переменных u2 и массив выходных переменных y2. Интервал дискретизации (промежутки времени, через которые производились измерения входных и выходных величин) в ходе эксперимента был принят равным 0. 08 с, который необходимо указать дополнительно в командной строке: >> ts=0.08 Этот массив данных при использовании в дальнейшем в пакете System
Identification
Toolbox
необходимо объединить в единый файл, содержащий необходимую информацию о входных и выходных параметрах объекта, их значениях и единицах измерения. >> dan=iddata(y2,u2,ts). Результат выполнения команды комментируется следующей фразой MATLAB: Time domain data set with 1097 samples. Sampling interval: 0.08 Outputs Unit (if specified) y1 Inputs Unit (if specified) u1 Сформированный файл dan. m указывает, что он содержит результаты 1097 измерений с интервалом дискретизации 0. 08 с. Входными переменными является массив u
1, а выходным параметром – y
1. Для наглядности сформированного файла необходимо в его структуру ввести обозначения входных и выходных данных c указанием размерностей параметров: >> dan.outputn='температура'; >> dan.inputn='расход газа'; >> dan.inputUnit='м3/ч'; >> dan.outputUnit='гр.С 100' В конечном итоге сформированный файл данных dan.m имеет следующий вид: Time domain data set with 1097 samples. Sampling interval: 0.08 Outputs Unit (if specified) температура гр.С 100 Inputs Unit (if specified) расход газа м3/ч Полную информацию о файле dan.m можно получить воспользовавшись командой: >> get(dan) Результат выполнения команды комментируется следующей фразой MATLAB: ans = Domain: 'Time' Name: [] OutputData: [1097x1 double] y: 'Same as OutputData' OutputName: {'температура'} OutputUnit: {'гр.С 100'} InputData: [1097x1 double] u: 'Same as InputData' InputName: {'расход газа'} InputUnit: {'м3/ч'} Period: Inf InterSample: 'zoh' Ts: 0.0800 Tstart: [] SamplingInstants: [1097x0 double] TimeUnit: '' ExperimentName: 'Exp1' Notes: [] UserData: [] ExperimentName: 'Exp1' Notes: [] UserData: [] Для графического представления данных используется команда plot (dan), либо команда idplot (datta), однако в последнем случае графики не будут содержать информации о названии переменных и их размерностях. Исходные данные с использованием команды plot (dan) приведены на рис. 2. 1. >> plot(dan) Для дальнейшего использования полученных исходных данных необходимо провести предварительную обработку этих данных с цель удаления тренда (постоянной соста Для удаления тренда пользуются функцией dtrend
: >> zdan=dtrend(dan) Исполнение функции приведет к появлению в командной строке записи: Time domain data set with 1097 samples. Sampling interval: 0.08 Outputs Unit (if specified) температура гр.С 100 Inputs Unit (if specified) расход газа м3/ч а в рабочей области Workspace системы MATLAB файла zdan. Рис. 2. 1. Исходные данные для идентификации технического объекта управления.
>> plot(zdan)). Файл в дальнейшем будет использован для построения моделей ТОУ. Рис.2. 2. Исходные данные для идентификации технического объекта управления при отсутствии тренда
Кроме указанной команды удаления тренда в пакете System Identification Toolbox имеются другие функции обработки данных эксперимента, которые приведены в описании пакета System Identification Toolbox. Применение этих функций производится в тех случаях, когда проведен предварительный анализ ТОУ и определены возможные помехи либо некоторые другие динамические характеристики, либо появляется необходимость изменить интервал дискретизации в случае повышенной погрешности представления модели ТОУ в ходе параметрического оценивания ТОУ. Следующим этапом идентификации является определение статистических и частотных характеристик массивов исходных данных. 2
. 7. Оценивание статистических и частотных характеристик исходных данных
Как уже отмечалось выше, при формировании массива исходных данных с использованием физического эксперимента над техническим объектом управления, воздействующий на объект входной сигнал был представлен случайным процессом с нулевым математическим ожиданием (т. е. центрированный после исключения тренда). Процесс будем считать эргодическим, что необходимо для практических приложений теории случайных процессов т. к. дает возможность по одной достаточно продолжительной реализации случайного процесса судить о его статистических характеристиках. В соответствие с свойствами стационарного эргодического процесса любая статистическая характеристика, полученная усреднением по ансамблю возможных реализаций, с вероятностью сколь угодно близкой к единице, может быть получена усреднением за достаточно большой промежуток времени из одной единственной реализации случайного процесса. Поэтому любая реализация исходных данных может быть использована нами для получения статистических характеристик массивов исходных данных т. к. в ходе планирования и проведения эксперимента сказать заранее, по какой реализации пойдет процесс, невозможно. Для характеристики связи между значениями случайного процесса в различные моменты времени, вводятся понятия корреляционной (автокорреляционной)
функции и спектральной плотности случайного процесса. Корреляционной функцией случайного процесса X
(t
) называют неслучайную функцию двух аргументов R
(t
1
; t
2
), которая для каждой пары произвольно выбранных значений моментов времени t
1
и t
2
равна математическому ожиданию произведения двух случайных величин случайного процесса, соответствующих этим моментам времени. Между дисперсией случайного процесса и корреляционной функцией существует прямая связь – дисперсия случайного стационарного процесса равна значению корреляционной функции. Статистические свойства связи двух случайных процессов X
(t
) и G
(t
) можно охарактеризовать взаимной корреляционной функцией R
xg
(t
1
, t
2
). Взаимная корреляционная функция Rxg
(τ
) характеризует взаимную статистическую связь двух случайных процессов X
(t
) и G
(t
) в разные моменты времени, отстоящие друг от друга на промежуток времени τ
. Если случайные процессы X
(t
) и G
(t
) статистически не связаны друг с другом и имеют равные нулю средние значения, то их взаимная Анализируя свойства корреляционной функции можно сделать вывод: чем слабее взаимосвязь между предыдущим X
(t
) и последующим X
(t
+
τ
) значениями случайного процесса, тем быстрее убывает корреляционная функция Rx
(τ
). Случайный процесс, в котором отсутствует связь между предыдущими и последующими значениями, называют чистым случайным процессом или белым шумом. В случае белого шума время корреляции τ
R
= 0 и корреляционная функция представляет собой δ-функцию. При исследовании автоматических систем управления удобно пользоваться еще одной характеристикой случайного процесса, называемой спектральной плотностью. Спектральная плотность S
x
(ω) случайного процесса X
(t
) определяется как преобразование Фурье корреляционной функции Rx
(τ
). Физический смысл спектральной плотности состоит в том, что она характеризует распределения мощности сигнала по частотному спектру. В пакете System Identification Toolbox имеется четыре функции cra
,
etfe
,
covf
,
и spa
непараметрического оценивания совокупности экспериментальных данных. Функция cra
выполняет расчет авто- и взаимных корреляционных функций, оценку импульсной характеристики методом корреляционного анализа для одномерного объекта массива экспериментальных данных. Написание этой функции следующее: cra(z); [ir,R,cl] = cra(z, M, na, plot); cra(R) Аргументы: · z – матрица экспериментальных данных вида z = [y2 u2], где y2 - вектор – столбец, соответствующий выходным данным; · u2 - вектор – столбец, соответствующий входным данным; · М – максимальное значение дискретного аргумента для которого производится расчет оценки импульсной характеристики (по умолчанию М = 20); · na – порядок модели авторегрессии (порядок многочлена), которая используется для расчета параметров отбеливающего фильтра (по умолчанию na = 10). При na = 0 в качестве идентифицирующего используется не преобразованный входной сигнал; · Если plot = 0, то график отсутствует, если plot = 1, то график полученной оценки импульсной характеристики вместе с 99% - м доверительным коридором, если plot = 2, то выводятся графики всех корреляционных функций. Возвращаемые величины: ir – оценка (вектор значений) импульсной характеристики; R – матрица, элементы первого столбца которой – значения дискретного аргумента, элементы второго столбца – значения оценки автокорреляционной функции выходного сигнала, элементы третьего столбца – значения оценки автокорреляционной функции входного сигнала, элементы четвертого столбца – значения оценки взаимной корреляционной функции. Для примера сушилки шликера эти величины имеют следующие значения: М и na приняты по умолчанию [], []. >> [ir,R,cl]=cra(zdan,[],[],2) 0.0134 0.1469 0.2256 0.1864 0.0956 0.0634 0.0457 0.0168 0.0066 0.0053 0.0046 0.0029 0.0068 -0.0068 -0.0099 -0.0099 -0.0017 0.0058 0.0150 0.0053 0.0081 R = -20.0000 0.0011 0.0015 -0.0123 -19.0000 0.0015 -0.0021 -0.0221 -18.0000 0.0017 0.0007 -0.0370 -17.0000 0.0017 0.0069 -0.0287 -16.0000 0.0013 0.0123 0.0080 -15.0000 0.0005 0.0074 0.0289 -14.0000 -0.0003 0.0051 0.0470 -13.0000 -0.0010 0.0092 0.0236 -12.0000 -0.0018 -0.0070 0.0419 -11.0000 -0.0019 0.0064 0.0221 -10.0000 -0.0010 -0.0008 0.0000 -8.0000 0.0001 0.0005 0.0018 -7.0000 0.0011 -0.0003 -0.0124 -6.0000 0.0031 0.0001 -0.0299 -5.0000 0.0065 0.0005 -0.0161 -4.0000 0.0110 0.0001 -0.0167 -3.0000 0.0163 -0.0001 0.0021 -2.0000 0.0261 -0.0007 0.0152 -1.0000 0.0393 0.0001 0.0259 0 0.0479 0.2477 0.0304 1.0000 0.0393 0.0001 0.3341 3.0000 0.0163 -0.0001 0.4239 4.0000 0.0110 0.0001 0.2174 5.0000 0.0065 0.0005 0.1442 6.0000 0.0031 0.0001 0.1040 7.0000 0.0011 -0.0003 0.0382 8.0000 0.0001 0.0005 0.0150 9.0000 -0.0005 0.0004 0.0121 10.0000 -0.0010 -0.0008 0.0105 11.0000 -0.0019 0.0064 0.0066 12.0000 -0.0018 -0.0070 0.0154 13.0000 -0.0010 0.0092 -0.0155 14.0000 -0.0003 0.0051 -0.0225 15.0000 0.0005 0.0074 -0.0224 16.0000 0.0013 0.0123 -0.0038 17.0000 0.0017 0.0069 0.0131 18.0000 0.0017 0.0007 0.0341 19.0000 0.0015 -0.0021 0.0119 20.0000 0.0011 0.0015 0.0185 cl = 0.0343 в г б а Рис. 2. 3. Графики функций: а) автокорреляционная функция выходного сигнала; б) автокорреляционная функция входного сигнала; в) взаимная корреляционная функция; г) импульсная характеристика.
Можно получить более подробный график импульсной характеристики, если выполнить функцию cra
с одним аргументом zdan (Рис. 2. 4): >> cra(zdan) ans = 0.0134 0.1469 0.2256 0.1864 0.0956 0.0634 0.0457 0.0168 0.0066 0.0053 0.0046 0.0029 0.0068 -0.0068 -0.0099 -0.0099 -0.0017 0.0058 0.0150 0.0053 0.0081 Рис. 2. 4. Импульсная характеристика
Необходимо отметить, что на графиках по оси абсцисс откладываются промежутки времени τ
= t
i
– t
i
-1
, а по оси ординат значения корреляционных функций для входного u
2
и выходного у2
сигналов; значения взаимокорреляционой функции и импульсной характеристики. Из полученных характеристик следует, что с увеличением τ
наблюдается резкий спад корреляционной зависимости входного сигнала, что свидетельствует о слабой взаимосвязи между сечениями процесса, соответствующими произвольным моментам времени (процесс более близок к белому шуму, а автокорреляционная функция к дельта-функции). Выходная величина наоборот более плавно изменяет свои состояния от одного момента времени к другому и, следовательно, взаимосвязь между предыдущим и последующим значениями выходного сигнала более тесная, чем у входного. Для получения частотных характеристик экспериментальных данных воспользуемся функциями оценивания частотных характеристик. Функция spa
возвращает частотные характеристики одномерного объекта [g, phiv] = spa(z) [g, phiv, z_spe] = spa(z,M,w,maxsize,T) Аргументы: · z – матрица исходных данных; · М – ширина временного окна (по умолчанию М = min(30, length(z) /10), где length(z) – число строк матрицы z); · w – вектор частот для расчета частотных характеристик (по умолчанию [1: 128]/128*pi/T); · Т – интервал дискретизации; · maxsize – параметр, определяющий максимальный размер матриц, создаваемых в процессе вычислений. Возвращаемые величины: · g – оценка W(e jωT
) в частотном формате; · phiv – оценка спектральной плотности шума v(t); · z_spe – матрица спектральных плотностей входного и выходного сигналов. Построим диаграмму Боде (АЧХ, ФЧХ), используя функции spa и bodeplot и данные, полученные при изучении теплового объекта, и содержащиеся в файле dryer2 >> load project24; >> z=[y2 u2]; >> g=spa(z); >> bodeplot(g) Результаты моделирования (АЧХ построена в логарифмическом масштабе) без доверительного коридора представлены на рис. 2. 5. Рис. 2. 5. Частотные характеристики технического объекта управления
Полученные зависимости подтверждают высокочастотную составляющую значений входного и выходного сигналов. Границы изменения частот на графиках установлены по умолчанию. Для получения частотных характеристик вместе с доверительным коридором шириной в три среднеквадратических отклонения в пакете System Identification Toolbox MATLAB имеется следующие возможности: · установка границ изменения частот с помощью команды >> w=logspace(w1
,w2
,N), где w1
– нижняя граница диапазона частот (10w
1
), w2
– верхняя граница диапазона частот (10w
2
) и N – количество точек графика. · построение АФХ, ФЧХ и S
(ω
) – функции спектральной плотности шума e(
t
)
· вычисление g – оценки АФХ и ФЧХ в частотном формате и phiv – оценки спектральной плотности шума с помощью команды >> [g,phiv]=spa(z,[],w); Графики АФХ, ФЧХ и S
(ω
) строятся с доверительным коридором в три среднеквадратических отклонения с помощью команды >> bodeplot([g p],'sd',3,'fill'), где 'sd' – указывает на сплошную линию доверительного коридора (по умолчанию эта линия штриховая); 3 – величина доверительного коридора в три среднеквадратических отклонения; 'fill' – способ заливки доверительного коридора (желтым цветом). Построим АЧХ, ФЧХ, используя функции spa, bodeplot, logspace и данные, полученные в файле Project24 с соответствующим доверительным коридором: >> w = logspace(-2,pi,128); >> [g,phiv]=spa(z,[],w); >> bodeplot([g,phiv],3,'fill') Результаты моделирования представлены на рис. 2. 6. Рис. 2.6. Оценки АЧХ и ФЧХ вместе с доверительным коридором.
>> bodeplot([phiv],'sd',3,'fill') Результаты моделирования представлены на рис. 2. 7. Рис.2. 7. График оценки
S(
ω) вместе с доверительным коридором
Полученный график оценки спектральной плотности шума с доверительным коридором показывает наличие равномерного распределения мощности сигнала по частотному спектру с последующим спадом мощности на частоте выше 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); >> 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); 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. Окно выбора структуры модели
Взамен строковой переменной '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 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 >> 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 В качестве аргумента функции может выступать любая из рассмотренных ранее моделей, например 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 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 Функция th2zp
рассчитывает нули, полюса и статические коэффициенты передачи (коэффициенты усиления) модели тета - формата zn4s многомерного объекта: >> [zepo,k]=th2zp(zn4s) 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
| ||||||||||||||||||||||||||||