2.3.Программные средства идентификации и моделирования систем При исследовании идентификационных моделей

Привет, Вы узнаете о том , что такое 2.3.Программные средства идентификации и моделирования систем При исследовании идентификационных моделей, Разберем основные их виды и особенности использования. Еще будет много подробных примеров и описаний. Для того чтобы лучше понимать что такое 2.3.Программные средства идентификации и моделирования систем При исследовании идентификационных моделей , настоятельно рекомендую прочитать все из категории Идентификация систем управления.

Пакет System Identification Toolbox содержит средства для создания математических моделей динамических систем на основе наблюдаемых входных и выходных данных. Он имеет гибкий графический интерфейс, помогающий организовать данные и создавать модели [2, 12, 29, 30]. Методы идентификации, входящие в пакет, применимы для решения широкого класса задач: от проектирования систем управления и обработки сигналов до анализа временных рядов и вибрации. Основные свойства пакета:

•простой и гибкий интерфейс;

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

•методы авторегрессии;

•анализ отклика во временной и частотной области;

•отображение нулей и полюсов передаточной функции системы;

•анализ невязок при тестировании модели;

•построение сложных диаграмм, таких как диаграмма Найквиста,

и др.

Графический интерфейс упрощает предварительную обработку данных, а также диалоговый процесс идентификации модели. Возможна также работа с пакетом в командном режиме и с применением расширения Simulink. Операции загрузки и сохранения данных, выбора диапазона, удаления смещений и трендов выполняются с минимальными усилиями и находятся в главном меню.

Представление данных и идентифицированных моделей организовано графически таким образом, что в процессе интерактивной идентификации пользователь легко может вернуться к предыдущему шагу работы. Для новичков существует возможность просматривать

следующие возможные шаги: графические средства позволяют отыскать любую из ранее полученных моделей и оценить ее качество в сравнении с другими моделями.

Начав с измерения выхода и входа, можно создать параметрическую модель системы, описывающую ее поведение в динамике. Пакет поддерживает все традиционные структуры моделей, включая авторегрессию, структуру Бокса–Дженкинса и др. Он поддерживает линейные модели пространства состояний, которые могут быть определены как в дискретном, так и в непрерывном пространстве. Эти модели могут включать произвольное число входов и выходов. В пакете есть функции, которые можно использовать как тестовые данные для идентифицированных моделей. Идентификация линейных моделей широко используется при проектировании систем управления, когда требуется создать модель объекта. В задачах обработки сигналов модели могут быть использованы для адаптивной обработки сигналов. Методы идентификации успешно применяются и для финансовых приложений.

2.3.1. Виды моделей пакета System Identification Toolbox

Дискретные модели, используемые в пакете System Identification Toolbox для временной области, учитывающие действие шума наблюдения, могут быть представлены следующим образом:

• модель авторегрессии AR (AutoRegressive):

A(z) y(t) = e(t),

где e(t) – дискретный белый шум; A(z) =1+ a1z−1 + a2 z−2 +K+ ana z−na ;

• ARX-модель (AutoRegressive with external input):

A(z) y(t) = B(z)u(t) +e(t),

где B(z) =1+b1z−1 +b2 z−2 +K+bnb z−nb ;

• ARMAX-модель (AutoRregressive-Moving Average with external input) – модель авторегрессии скользящего среднего:

y(t) +a1 y(t −1) +K+ ana y(t −n) =

= b1u(t) +b2u(t −1) +K+bnb y(t −m) +e(t), A(z) y(t) = B(z)u(t −nk) +C(z)e(t),

81

где nk – величина задержки(запаздывания);

C(z) =1+c1z−1 +c2 z−2 +K+cnc z−nc ;

• модель «вход–выход» («Output–Error»):

y(t) =

B(z)

u(t −nk) +e(t) ,

F(z)

где F(z) =1+ f z−1

+ f

2

z−2 +K+ f

nf

z−nf ;

1

• модель Бокса–Дженкинса (BJ):

y(t) =

B(z)

u(t

−nk) +

C(z)e(t) ,

F(z)

D(z)

где D(z) =1+ d1z−1 + d2 z−2 +K+dnd z−nd .

Данные модели можно рассматривать как частные случаи обобщенной параметрической линейной структуры:

A(z) y(t) = FB((zz))u(t −nk) + CD((zz))e(t).

При этом все они допускают расширение для многомерных объектов (имеющих несколько входов и выходов);

• модель для переменных состояния (State-space): x(t +1) = Ax(t) + Bu(t),

y(t) = C(x(t) + Du(t) +v(t),

где A, B, C, D – матрицы коэффициентов, входа, выхода и обхода системы; v(t) – коррелированный белый шум наблюдений. Возможна и другая (каноническая) форма представления данной модели:

x(t +1) = Ax(t) + Bu(t) + Ke(t),

y(t) = C(x(t) + Du(t) +e(t),

где K – некоторая весовая матрица; e(t) – дискретный белый шум. В своей работе пакет System Identification Toolbox использует три

внутренних вида матричного представления моделей, которые с помощью операторов и функций пакета преобразуются во все вышеперечисленные виды моделей объектов:

82

•тета-формат (для временных моделей);

•частотный формат (для частотных моделей);

•формат нулей и полюсов.

2.3.2.Использование пакета System Identification Toolbox

для построения модели системы

Основные операторы и функции пакета System Identification Toolbox набираются в командной строке MatLab или могут быть использованы при написании программ для т-файлов.

Список основных операторов и функций и их описание приведены в приложении 1.

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

восновном окне MatLab.

Вкачестве примера использования пакета System Identification Toolbox для идентификации технических систем рассмотрена распылительная сушилка [12].

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

ром входным воздействием на систему является расход газа, м3/ч, а выходным регулируемым параметром – температура, °С. Процесс идентификации включает следующие этапы:

•априорный анализ объекта с целью выбора структуры модели;

•проведение предварительного исследования объекта с целью уточнения оценки структуры модели;

•разработка методики основного экспериментального исследования объекта, составление плана эксперимента;

•проведение основного экспериментального исследования для получения массива данных ( ui , yi );

•математическая обработка массива данных (с использованием пакета System Identification Toolbox) с целью определения параметров модели и ее адекватности, доверительных границ параметров и выходной координаты модели.

При этом в процессе исследования системы необходимо принять

некоторые допущения, позволяющие применить хорошо отработанный аппарат анализа стационарных, линейных объектов:

• технический объект является системой с сосредоточенными параметрами;

83

•технический объект стационарен, т.е. его статические и динамические свойства неизменны во времени;

•уравнения моделей изучаемой системы линеаризуются в малом, т.е. при небольших отклонениях ±∆yi от выбранной «рабочей» точки

(рабочего режима) объекта.

Массив данных ui , yi образуется в результате трудоемкой операции

расшифровки регистрограмм по приборам измерительной системы. Однако широкое развитие микропроцессорной и вычислительной техники и внедрение ее в производственные технологические процессы позволили существенно усовершенствовать техническое обеспечение идентификации системы. Обработка массива данных с помощью пакета System Identification Toolbox предполагаетследующие этапы:

•обработка и преобразование данных с целью создания файла;

•анализ экспериментальных данных с целью предварительного определения основных характеристик объекта;

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

•задание структуры модели;

•изменение и уточнение структуры модели (если это необходимо);

•проверка адекватности и сравнение различных моделей с целью выбора наилучшей;

•преобразование модели из тета-формата в вид, удобный для дальнейшего использования при анализе и синтезе системы управления.

На каждом этапе идентификации имеется возможность графического отображения результатов моделирования и извлечения необходимой информации об объекте.

2.3.2.1. Обработка данных при построении модели объекта

Построение аналитической модели технического объекта начинается с описания принципа функционирования объекта – сушилки. Всушилку подводится шликер, где он распыляется. Инжекционные горелки создают высокую температуру в зоне распыления материала. Распыленные частицы, теряя влагу, уже в виде порошка собираются в днище сушилки, откуда поступаютнепосредственно с бункера над прессами.

84

Основная задача состоит в построении модели в виде передаточной функции. При этом необходимо учитывать наличие также случайных возмущений (т.е. неотвратимых в реальной обстановке факторов), воздействие которых должно учитываться в моделировании. В частности, к таким возмущениям могут относиться случайные изменения интенсивности подачи шликера в сушилку или случайные изменения качества топлива в горелках и т.д. Эти изменения необходимо учесть при построении модели сушилки. Теоретически задача будет сводиться к изучению поведения системы с учетом воздействия на нее случайных факторов.

Аналитическая модель сложных систем может быть построена на основе массива входных и выходных данных, полученных в результате физического эксперимента, проведенного над системой. При этом воздействие случайных факторов должно быть смоделировано.

Статистические данные, полученные при изучении теплового объекта, помещаются в отдельный файл (под именем Project), который включает массив данных, состоящий из 1000 значений входного параметра – расхода газа, м3/ч, и 1000 значений выходного параметра – температуры в объекте, °С.

Для загрузки файла в рабочую среду Workspace системы MatLab необходимо в режиме командной строки выполнить следующую команду:

>>load Project.

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

менных y2. Интервал дискретизации (T0 = 0,08 ) необходимо указать дополнительно в командной строке:

>> ts = 0,08.

Этот массив данных при использовании в дальнейшем в пакете System Identification Toolbox необходимо объединить в единый файл, содержащий необходимую информацию о входных и выходных параметрах объекта, их значениях и единицах измерения.

Для объединения исходных данных в единый файл dan.m пользуются командой

>> dan=iddata(y2,u2,ts).

85

Результат выполнения команды комментируется следующей фра-

зой MatLab:

Тimе domain data set »ith 1000 samples. Sampling interval: 0,08 Outputs Unit (if specified) у1 Inputs Unit (if specified) u1.

Сформированный файл dan.m указывает, что он содержит результаты 1000 измерений с интервалом дискретизации 0,08 с. Входной переменной является массив u1, а выходным параметром – y1.

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

>>dan.outputn='температура'; >>dan.input='расход газа'; >>dan.inputUnit='м3/ч', >>dan.outputUnit='гр,С 100'.

В итоге сформированный файл данных dan.m имеет следующий

вид:

Time domain data set with 1000 samples. Sampling interval: 0,08 Оutputs Unit (if specified) температура гр.С 100 Inputs Unit (if specified) расход газа м3/ч.

Полную информацию о файле dan.m можно получить, воспользовавшись командой

>>get(dan).

Для графического представления данных используется команда plot (dan) либо команда idplot (datta), однако в последнем случае графики не будут содержать информации о названии переменных и их размерностях. Исходные данные с использованием команды >>plot(dan) приведены на рис. 2.2.

Для дальнейшего использования полученных исходных данных необходимо провести их предварительную обработку с целью удаления тренда (постоянной составляющей) из набора данных и, если необходимо, отфильтровать данные с помощью имеющихся средств в пакете System Identification Toolbox.

Для удаления тренда пользуются функцией dtrend:

>> zdan=dtrend(dan).

86

Рис. 2.2. Исходные данные для идентификации объекта управления

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

Тimе domain data set with 1000 samples. Sampling interval: 0.08 Outputs Unit (if specified) температура гр.С 100 Inputs Unit (if specified) расход газа м3/ч.

В рабочей области Workspace системы появится файл zdan. Получен новый файл zdan.m, в котором отсутствует постоянная

составляющая сигналов (рис. 2.3). Рисунок получен после выполнения команды

>> plot(zdan)).

Файл в дальнейшем будет использован для построения моделей объекта управления.

Кроме указанной команды удаления тренда в пакете System Identification Toolbox имеются другие функции обработки данных эксперимента, которые приведены в описании пакета System Identification Toolbox. Применение этих функций производится в тех случаях, когда проведен предварительный анализ объекта идентификации и определены возможные помехи либо некоторые другие динамические харак-

87

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

Рис. 2.3. Исходные данные для идентификации объекта при отсутствии тренда

Следующим этапом идентификации является определение статистических и частотных характеристик массивов исходных данных.

2.3.2.2.Оценивание статистических

ичастотных характеристик исходных данных

Как уже отмечалось выше, при формировании массива исходных данных с использованием физического эксперимента над техническим объектом управления воздействующий на объект входной сигнал был представлен случайным процессом с нулевым математическим ожиданием (т.е. центрированный после исключения тренда). Процесс будем считать эргодическим, что необходимо для практических приложений теории случайных процессов, так как дает возможность по одной достаточно продолжительной реализации случайного процесса судить о его статистических характеристиках. В соответствии со свойствами стационарного эргодического процесса любая статистическая характеристика, полученная усреднением по ансамблю возможных

88

реализаций, с вероятностью, сколь угодно близкой к единице, может быть получена усреднением за достаточно большой промежуток времени из одной единственной реализации случайного процесса. Поэтому любая реализация исходных данных может быть использована для получения статистических характеристик массивов исходных данных, поскольку в ходе планирования и проведения эксперимента сказать заранее, по какой реализации пойдет процесс, невозможно.

Для характеристики связи между значениями случайного процесса в различные моменты времени применяются корреляционная функция Rx (τ) и спектральная плотность случайного процесса Sx (ω) .

Статистические свойства связи двух случайных процессов X (t) и G(t) можно охарактеризовать взаимной корреляционной функцией

Rxg (τ) .

В пакете System Identification Toolbox имеется четыре функции cra, etfe, conf и spa непараметрического оценивания совокупности экспериментальных данных. Функция сrа выполняет расчет авто- и взаимных корреляционных функций, оценку импульсной характеристики методом корреляционного анализа для одномерного объекта массива экспериментальных данных. Написание этой функции следующее:

cra(z);

[ir,R,cl] = cra(z, М, na, plot), cra(R).

Здесь z – матрица экспериментальных данных вида z = [y2 u2]; у2 – вектор-столбец, соответствующий выходным данным; u2 – векторстолбец, соответствующий входным данным; М – максимальное значение дискретного аргумента, для которого производится расчет оценки импульсной характеристики (по умолчанию М = 20); na – порядок модели авторегрессии (порядок многочлена), которая используется для расчета параметров отбеливающего фильтра (по умолчанию na = 10). При na = 0 в качестве идентифицирующего используется непреобразованный входной сигнал; plot – наличие графика:

если plot = 0, то график отсутствует,

eсли plot = 1, то график полученной оценки импульсной характеристики вместе с 99%-м доверительным интервалом,

89

если plot = 2, то выводятся графики всехкорреляционных функций. Возвращаемые величины:

ir – оценка (вектор значений) импульсной характеристики;

R – матрица, элементы первого столбца которой – значения дискретного аргумента, элементы второго столбца – значения оценки автокорреляционной функции выходного сигнала, элементы третьего столбца – значения оценки автокорреляционной функции входного сигнала, элементы четвертого столбца – значения оценки взаимной корреляционной функции.

Например, для сушилки шликера эти величины имеют следующие значения.

Ми na приняты по умолчанию [], [].

>>[ir,R,cI]=cra(zdan,[], [],2).

На рис. 2.4 приведены результаты расчета автокорреляционной функции выходного сигнала (Covf for filtered y); автокорреляционной функции входного сигнала (Covf for prewhitened u); взаимной корреляционной функции (Correlation from u to у); импульсной характери-

стики (Impulse response estimate).

Можно получить более подробный график импульсной характеристики, если выполнить функцию cra с одним аргументом zdan

(рис. 2.5):

>> cra(zdan).

Необходимо отметить, что на графиках по оси абсцисс откладываются промежутки времени τ = ti −ti−1 , а по оси ординат – значения кор-

реляционных функций для входного u2 и выходного у2 сигналов; значения взаимокорреляционой функциии импульсной характеристики.

Из полученных характеристик следует, что с увеличением τ наблюдается резкий спад корреляционной зависимости входного сигнала, что свидетельствует о слабой взаимосвязи между сечениями процесса, соответствующими произвольным моментам времени (процесс более близок к белому шуму, а автокорреляционная функция к дельтафункции). Выходная величина, наоборот, более плавно изменяет состояния от одного момента времени к другому, и, следовательно, взаимосвязь между предыдущим и последующим значениями выходного сигнала более тесная, чем у входного.

аб

вг

Рис. 2.4. Графики функций: а – автокорреляционная функция выходного сигнала; б – автокорреляционная функция входного сигнала; в – взаимная корреляционная функция; г – импульсная характеристика

Рис. 2.5. Импульсная характеристика

Для получения частотных характеристик экспериментальных данных воспользуемся функциями оценивания частотных характеристик. Об этом говорит сайт https://intellect.icu . Функция spa возвращает частотные характеристики одномерного объекта и оценки спектральной плотности его сигналов для обобщенной линейной модели объекта:

91

[g, phiv]= spa(z);

[g, phiv, z_spe]=spa(z, М, maxsize, Т);

здесь z – матрица исходных данных; М – ширина временного окна (по умолчанию М = min(30, length(z)/10), где length(z) – число строк мат-

рицы z); Т – интервал дискретизации; maxsize – параметр, определяющий максимальный размер матриц, создаваемых в процессе вычислений.

Возвращаемые величины:

g – оценка W (e jωT ) в частотном формате;

phiv – оценка спектральной плотности шума v(t) ;

z_spe – матрица спектральных плотностей входного и выходного сигналов.

Диаграмма Боде (АЧХ, ФЧХ) строится с использованием функций spa и bodeplot и данных, полученных при изучении объекта и содержащихся в файле Project:

>>load Project;

>>х=[у2 u2];

>>S=spa(z);

>>bodeplot(g).

Результаты моделирования без доверительного интервала представлены на рис. 2.6.

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

Для получения частотных характеристик вместе с доверительным интервалом шириной в три среднеквадратических отклонения в пакете System Identification Toolbox MatLab имеются следующие возможности:

•установка границ изменения частот с помощью команды

>>w=logspace(w1, w2, N),

где w1 – нижняя граница диапазона частот (10w1 ); w2 – верхняя граница диапазона частот (10w2 ); N – количество точек графика;

92

Рис. 2.6. Частотные характеристики объекта управления

•построение АФХ, ФЧХ и S(ω) – функции спектральной плотности шума e(t) ;

•вычисление g-оценки АФХ и ФЧХ в частотном формате и phivоценки спектральной плотности шума с помощью команды

>>[g, phiv]=spa(z, [], w).

Графики АФХ, ФЧХ и S(ω) строятся с доверительным интервалом в три среднеквадратических отклонения с помощью команды

>> bodeplot([g р], 'sd',3,'fill'),

где 'sd' указывает на сплошную линию доверительного интервала (по умолчанию эта линия штриховая); 3 – величина доверительного коридора в три среднеквадратических отклонения; 'fill' – способ заливки доверительного коридора (желтым цветом).

Построение АЧХ, ФЧК с использованием функций spa, bodeplot, logspace и данных, полученных в файле Project с соответствующим доверительным коридором, осуществляется следующим образом:

>> w =logspace(-2,pi,128); >>[jg, phiv]=spa(z, [], w);

93

>> bodeplot([g, phiv], 3, 'fill').

Результаты моделирования представлены на рис. 2.7.

Рис. 2.7. Оценки АЧХ и ФЧХ вместе с доверительным интервалом

Для построения графика оценки спектральной плотности шума

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

>>bodeplot([phiv], 'sd', 3, 'fill').

Результаты моделирования показаны на рис. 2.8.

Рис. 2.8. Оценки S(ω) вместе с доверительным интервалом

94

Полученный график оценки спектральной плотности шума с доверительным интервалом свидетельствует о наличии равномерного распределения мощности сигнала по частотному спектру с последующим спадом мощности на частоте выше 1,1 рад/с.

Далее необходимо выполнить параметрическое оценивание.

2.3.3.3. Параметрическое оценивание данных

Параметрическое оценивание экспериментальных данных проводится с целью определения параметров модели заданной структуры путем минимизации выбранного критерия качества модели (чаще всего – среднего квадрата рассогласования выходов объекта и его постулируемой модели).

Для проведения параметрического оценивания массив экспериментальных данных необходимо разделить условно на две части (необязательно равные):

>>zdanv=zdan(I:500);

>>zdane=zdan(501:1000).

Первая часть массива данных будет использоваться для параметрического оценивания и построения модели системы. Вторая часть будет необходима для верификации (проверки качества) модели, оценки адекватности полученной модели и определения погрешностей идентификации. Необходимо отметить, что параметрическая идентификация в пакете System Identification Toolbox выполняется в дискретном виде и полученные модели являются дискретными.

В пакете System Identification Toolbox рассмотрены различные виды моделей, которые с разной степенью достоверности описывают объект идентификации. Для выбора наиболее приемлемой структуры и вида моделей при параметрическом оценивании экспериментальных данных в пакете System Identification Toolbox имеются специальные функции:

•параметрического оценивания;

•задания структуры модели;

•изменения и уточнения структуры модели и выбора структуры модели.

Функция оценивания ar определяет параметры авторегрессии:

95

A(z) y(t) = e(t),

где A(z) =1+ a z−1

+ a z−2

+Ka

na

z−na , т.е. коэффициенты полинома

1

2

A(z) при моделировании скалярных временных последовательностей. Функция имеет синтаксис

th = аr(у, n)

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

[th, refl]=ar(y, n, approach, win, maxsize, Т),

где у – вектор-столбец данных, содержащий N элементов; n – порядок модели (число оцениваемых коэффициентов); approach – аргумент (строковая переменная), определяющая метод оценивания:

•'1s' – метод наименьших квадратов;

•'yw' – метод Юла-Уокера;

•'burg' – метод Бэрга (комбинация метода наименьших квадратов

сминимизацией гармонического среднего);

•'gl' – метод с использованием гармонического среднего.

Если любое из данных значений заканчивается нулем (например, burg0), то вычисление сопровождается оцениванием корреляционных функций.

win – (строковая переменная) используется в случае отсутствия части данных:

•win='now' – используются только имеющиеся данные (применяется по умолчанию, за исключением случая approach = 'yw');

•win='prw' – отсутствующие начальные данные заменяются нулями, так что суммирование начинается с нулевого момента времени;

•win='prw' – последующие отсутствующие данные заменяются нулями, так что суммирование расширяется до момента времени N+n;

•win='ppw' – и начальные, и последующие отсутствующие данные заменяются нулями (используется в алгоритме Юла-Уокера);

maxsize определяет максимальную размерность задачи; Т – интервал дискретизации.

Возвращаемые величины:

th – полученная модель авторегрессии в тета-формате (внутреннем матричном формате представления параметрических моделей па-

кета System Identification);

96

relf – информация о коэффициентах и функции потерь.

Для использования функции параметрического оценивания ar необходимо из массива экспериментальных данных, записанных в файле dan, выделить выходную переменную у с помощью команды

>>y=dan.у;

что равносильно команде

>>y=get(dan,'у'); >>th=ar(y,4);

Discrete timе IDPOLY model: A(q)y(t)=e(t)

A(q)=1 -2,195q^-1+1.656q^-2-0.4167q^-3-0.04415q^- 4 Estimated using AR ('fb'/'now') from data set y

Loss function 0.0104048 and FРЕ 0.0104884 Sampling interval: 1.

Полная информация о модели авторегрессии th может быть получена с помощью команды

>> prcsent(th).

Discrete-time IDPOLY model: A(q)y(t)=e(t) A(q)=1-2.195(±0.03165 q^-1+1.656(±0.07533 q^-2- -0.4167(±0.07542)q^-3-0.04415(±0.03174)q^- 4 Estimated using АК ('fb'/'now') from data set y Loss function 0.0104048 and FРЕ 0.0104884 Sampling interval: 1

Created: 07-Jan-2012 08:47:27

Last modified: 07-Jan-2012 08:47:27.

В информации приведены сведения о том, что модель является дискретной и для оценивания ее параметров используется прямой обратный метод (разновидность метода наименьших квадратов), на это указывает строковая переменная 'fb' (используется по умолчанию); для построения модели необходимы только имеющиеся данные у, на что указывает строковая переменная 'now' (используется по умолчанию); определены: функция потерь Loss function как остаточная сумма квадратов ошибки и так называемый теоретический информационный кри-

терий Акейке (Akaike's Information Theoretic Criterion – AIC) FPE; ин-

тервал дискретизации Sampling interval.

97

Следующая функция arx оценивает параметры модели AR и ARX, представленной зависимостью

A(z) y(t) = B(z)u(t) +e(t)

или в развернутом виде:

y(t) + a1 y(t −1) +K+ ana y(t −n) =

= b1u(t) +b2u(t −1) +K+bnb y(t −m) +e(t).

Здесь и ниже e(t) – дискретный белый шум.

B(z) =1+b1z−1 +b2 z−2 +K+bnb z−nb .

Функция имеет следующий синтаксис: dar=arx(z,nn)

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

dar=arx(z,nn,maxsize,Т),

где z – экспериментальные данные; nn – задаваемые параметры модели (аргумент nn содержит три параметра: na – порядок (число коэффициентов) полинома A(z); nb – порядок полинома В(z); nk – величина задержки); maxsize – максимальная размерность задачи; Т – интервал дискретизации.

При выборе полинома возникает вопрос о его степени. Известно, что с увеличением порядка полиномов улучшается степень адекватности модели реальному объекту. Однако при этом получаются громоздкие выражения и увеличивается время моделирования. Поэтому для нахождения оптимального порядка полиномов можно воспользоваться функциями выбора структуры модели.

Функция arxstruc вычисляет функции потерь для ряда различных конкурирующих ARX-моделей с одним выходом:

v=arxstruc(ze,zv,NN)

или

v=arxstruc(ze,zv,NN,maxsize);

где ze, zv – соответственно матрицы экспериментальных данных для оценивания и верификации моделей; NN – матрица задания конкурирующих структур со строками вида nn=[na nb nk]; maxsize – максимальная размерность задачи.

98

Возвращаемая величина v – матрица, верхние элементы каждого ее столбца (кроме последнего) являются значениями функции потерь для ARX-моделей, структура которых отображается последующими элементами столбцов (т.е. каждый столбец соответствует одной модели). Первый элемент последнего столбца – это число значений экспериментальных данных для верификации моделей.

Функция selstruc осуществляет выбор наилучшей структуры модели из ряда возможных вариантов:

[nn,vmod]=selstruc(v), [nn,vmod]=selstruc(v,с),

где v – матрица, возвращаемая функцией arxstruc; с – строковая переменная, определяющая вывод графика или критерий отбора наилучшей структуры:

•с='plot' – выводится график зависимости функции потерь от числа оцениваемых коэффициентов модели;

•с='log' – выводится график логарифма функции потерь;

•с='aic' – график не выводится, но возвращается структура, минимизирующая теоретический информационный критерий Акейке;

•с='mdl' – возвращается структура, обеспечивающая минимум критерия Риссанена минимальной длины описания;

•с=a – выбирается структура, которая минимизирует значение функции потерь:

vmod = v(1+a(d / N)),

где v – значение функции потерь; d – число оцениваемых коэффициентов модели; N – объем выборки экспериментальных данных, используемых для оценивания.

Возвращаемые величины: nn – выбранная структура;

vmod – значение соответствующего критерия. Например, для данных Project можно:

• задать пределы изменения порядка модели: >>NN=struc(1:10,1:10,1);

• вычислить функции потерь: >>v=arxstruc(zdane,zdаnv,NN);

99

• выбрать наилучшую структуру порядков полиномов: >>[nn,vmod]=selstruc(v,'plot'),

где 'plot' – строковая переменная, определяющая вывод графика зависимости функции потерь от числа оцениваемых коэффициентов моде-

ли (рис. 2.9).

Рис. 2.9. Окно выбора структуры модели

В появившемся окне столбики указывают на величину функции потерь. При подведении курсора к соответствующему столбику в правом поле окна отразятся значения порядков полиномов na, nb, nk. В поле графика появятся рекомендации по выбору цвета столбика.

Взамен строковой переменной 'plot' возможны варианты:

•'log' – выводится график логарифма функции потерь;

•'aic' – график не выводится, но возвращается структура, минимизирующая теоретический информационный критерий Акейке

(Akaike's Information Theoretic Criterion – AIC) FPE vmod:

vmod =

v

,

+2(d

))

(1

N

где v – значение функции потерь; d – число оцениваемых коэффициентов модели; N – объем экспериментальных данных, используемых для оценивания;

• 'mdl' – возвращается структура, обеспечивающая минимум критерия Риссанена минимальной длины описания (Rissanen's Minimum Description Lеngth):

vmod =

v

;

+log(N)(d

))

(1

N

• при строковой переменной, равной некоторому численному значению а, выбирается структура, которая минимизирует:

vmod =

v

.

+a(d

))

(1

N

Выбор наилучшей структуры порядков полиномов можно осуществить и с помощью более простой команды:

>>nn=selstruc(v,0)

MatLab возвращает

nn=7 9 1.

С учетом выбранной структуры модели определим вид модели ARX, выполнив функцию arx:

>>darx=arx(zdanv,nn).

Возвращается матрица из 100 столбцов и 4 строк со значениями различных критериев:

vmod=…………. .

Возвращается дискретная модель, представленная в тета-формате (внутренний вид матричных моделей).

Discrete-time IDPOLY model: А(q)у(t)=B(q)u(t)+e(t)

A(q)=1- 0.9601q^-1+0.03384q^-2-0.1035 q^-3+0.089q^-4- 0.02827^q-5+

+0.1383q^-6-0.06615q^-7 B(q)=0.0009708q^-1+0.003863q^-2+0.06428q^-3+0.06245q^-4+ +0.01997q^-5-0.007896q^-6-0.01869q^-7-0.01931 q^-8-0.007944q^-9 Estimated using ARX from data set zdanv

Loss function 0.00139284 and FPE 0.00148493 Sampling interval: 0.08.

101

Функция armax оценивает параметры ARMAX-модели: >>darmax=armax(zdanv,[2 2 2 1]).

Аргументы функции:

zdanv – вектор экспериментальных данных;

[nа nb nc nk] – степени полиномов и величина задержки. Возвращается дискретная модель в тета-формате:

Discrete-time IDPOLY model: A(q)y(t=B(q)u(t)+C(q)e(t) A(q=1-1.648q^-1+0.7054q^'-2 B(q)=0.005017q^-1+0.04556q^-2 C(q)=1-0.278q^-1-0.4333q^-2

Estimated using ARMAX from data set zdanv Loss function 0.00723841 and FPE 0,00741677 Sampling interval: 0.08.

Функция ое оценивает параметры ОЕ-модели: >>zoe=oe(zdanv,[2 2 I]).

Возвращается дискретная модель, представленная в тета-формате:

Discrete time IDPOLY

продолжение следует...

Продолжение:

Часть 1 2.3.Программные средства идентификации и моделирования систем При исследовании идентификационных моделей
Часть 2 - 2.3.Программные средства идентификации и моделирования систем При исследовании идентификационных