Работа с пакетом в режиме командной строки

Работа с пакетом в режиме командной строки

Возможности работы в режиме командной строки

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

• графического интерфейса;

• обработки и преобразования данных;

• отображения модели;

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

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

 

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

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

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

• выбора структуры модели;

• преобразования модели;

• извлечения информации о модели;

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

• прочие;

• демонстрации возможностей пакета.

Команды графического интерфейса

Данная группа содержит две команды: ident и midprefs.

Команда ident — команда открытия рассмотренного выше графиче­ского интерфейса пользователя. Она записывается в виде:

ident

ident(session, directory)

Команда открывает окно графического интерфейса. При второй фор­ме записи открытие интерфейса сопровождается загрузкой сессии из файла с именем session, находящегося в директории directory. При­меры:

ident('iddata1, sid')

ident(‘mydata,sid', '\matlab\data\cdplayer\')

Команда midprefs задает (изменяег) директорию для файла midprefs.mat, хранящего информацию о начальных параметрах графического ин­терфейса пользователя при его открытии. Пример:

midprefs('/home/l jung/mat1ab/')

Функции обработки и преобразования данных

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

Функция detrend

Функция detrend удаляет тренд из набора данных. Она используется в виде:

zd - detrend(z)

zd - detrend(z,o,brkp)

Ее аргументы: z — матрица исходных данных, в которой столбцы со­ответствуют векторам данных; o (по умолчанию 0) — параметр, оп­ределяющий порядок тренда (нулевой или первый); если аргумент brkp не определен, из вектора данных вычитается линейная функция, соответствующая тренду; если данный аргумент определен, из исход­ныхданных вычитается «ломаная линия» с точками разрыва произ­водной при элементах, номера которых содержатся в brkp (только для линейных трендов). Возвращаемая величина: zd — матрица преобра­зованных данных.

В следующем примере:

 

zd(:,1) = detrend(z(:,1),1,119);

zd(:,2) = detrend(z(:,2));

 

удаляется V-образный тренд из данных первого столбца с пиком на 119-м элементе и тренд нулевого порядка (постоянная составляю­щая) — для элементов второго столбца матрицы исходных данных.

 

Функция idfilt

 

zf = idfilt(z,ord,Wn)

[zf, thf] = idfilt(z,ord,Wn,hs)

 

возвращает данные, фильтрованные с помощью фильтра Баттерворта. Эта функция создает фильтр Баттерворта порядка ord и осуществляет фильтрацию данных (по столбцам) матрицы данных z = [у u]. Если аргумент hs не задан и Wn содержит только один элемент, создается фильтр низких частот (ФНЧ) с частотой среза Wn.

Если hs -'high', создается фильтр верхних частот (ФВЧ) с частотой среза Wn.

Если Wn - [Wnl Wnh], то есть является двухэлементным вектором, и аргумент hs не определен, создается полосовой фильтр (ПФ) с полосой пропускания от Wnl до Wnh. Если же в этом случае задан аргумент hs = 'stop', создается заграждающий полосовой фильтр. Порядок фильтра в этих случаях — 2*ord.

Отметим, что частоты Wnl и Wnh задаются в относительных единицах, в долях от так называемой частоты Найквиста, определяемой как ω= π/Т (или f = 1/2T), где Т — выбранный интервал дискретизации.

Возвращаемые величины: zf — матрица отфильтрованных данных, thf — описание созданного фильтра в так называемом тета-формате (см. ниже).

 

Функция idinput

Функция idinput в формах

u = idinput(N)

u =idinput(N,type,band,levels),

u = idinput(N,'sine',band,levels,auxvar)

генерирует входной идентифицирующий сигнал. Она создает несколько разновидностей входного сигнала, наиболее часто используемых в задачах идентификации при возможности организации активного эксперимента. Аргумент N задает количество значений сигнала. Строковая переменная type задает тип сигнала и может принимать значения (по умолчанию — 'rbs'):

• type = 'rs' — задание случайного гауссова сигнала;

• type - 'rbs' — задание случайного бинарного (двоичного) сигнала;

• type = ‘prbs' — задание псевдослучайного двоичного сигнала (ПСДС);

• type * 'sine' — задание полигармонического сигнала (в виде сум­мы синусоид).

Двухэлементный вектор-строка band * [wlow, whigh] задает диапазон частот сигнала (для сигналов типа 'rs', *rbs' и 'sine'), при этом wlow — нижняя частота, a whigh — верхняя частота (частоты задаются в долях от частоты Найквиста — см. выше). Белый шум, таким обра­зом, характеризуется значением band = [0, 1], которое является и значением по умолчанию.

Для сигнала типа ПСДС ('prbs') используется значение band = [twologp, М], при этом 1/М — длительность интервала постоянства сигнала, а значение 2^twologp-1 определяет период ПСДС. Макси­мальное значение периода соответствует twologp=18 (или twologp=0); в этом случае по умолчанию band = [0, 1].

Двухэлементный вектор levels - [minu. maxu] задает минимальный и максимальный уровни для сигналов типа 'rs¹, 'rbs' и 'sine'. Для сигнала типа 'rs' величина minu имеет значение, равное математическому ожиданию минус среднее квадратическое (или стандартное) отклонение, a maxu — математическому ожиданию плюс среднее квадратическое отклонение.

Для сигнала типа 'sine' возможно задание четвертого аргумента — auxvar = [no_of_sinusoids, no_of_trials], где no_of_sinusoids — число гармоник, образующих сигнал, a no_of_trials определяет, сколько итераций должно быть сделано, чтобы минимизировать амплитуду си­нала путем задания различным синусоидам произвольных начальных фаз. Значение но умолчанию — auxvar = [10, 10].

Возвращаемая величина — N-элементный вектор и.

Функция idresamp

Функция idresamp, записываемая в формах

zr = idresamp(z,R)

[zr,R_act] = idresamp(z,R,filter_order,tol)

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

Аргументы: z — матрица исходных данных (каждый столбец соответствует вектору данных), R — коэффициент изменения интервала дискретизации; при R>1 происходит прореживание исходного сигна­ла, при R<1 — его интерполяция. Допустимо любое положительное значение R, но оно будет заменено достигнутым приближением R_ act.

• filterorder — порядок фильтра, который будет использован для предварительной фильтрации исходного сигнала (по умолча­нию — 8);

• tol — допустимое относительное отклонение R_act от R. По умолчанию — 0.1.

Возвращаемые величины — матрица преобразованных сигналов тс и достигнутое значение R_act.

Функции отображения модели

 

Эта группа образована шестью функциями.

Функция bodeplot

bodeplot — функция отображения логарифмических частотных ха­рактеристик (диаграмма Боде):

hodeplot(g)

bodeplot(g1 g2... gn])

bodeplot(g,sd,C,mode)

Она имеет следующие аргументы:

• g, gl, g2,..., gn — имена оценок частотных характеристик;

• sd — при sd>0 к графикам добавляются соответствующие доверительные коридоры (по умолчанию sd=O);

• аргумент С задает выводимые характеристики:
o С = 'А' — выводится только график АЧХ;
o С = 'Р — выводится только график ФЧХ;
o С = 'В' — выводятся оба графика;

• для получения всех графиков на одной и той же диаграмме зада­ется mode = 'same'.

Возвращаются графики частотных характеристик с логарифмиче­скими масштабами для АЧХ и оси частот.

Функция ffplot

Функция ffplot имеет такое же назначение, как предыдущая. Запи­сывается в виде:

ffpiot(g)

ffplot([gl g2... gn])

ffplot(g,sd,C,mode)

Аргументы — как у предыдущей функции. Возвращаются графики частотных характеристик с логарифмическим масштабом для АЧХ, но линейным для оси частот (единицы частоты — Гц).

Функция idplot

 

Функция idplot графически отображает входные и выходные данные:

idplot(z)

idplot(z,int,T,ny,pc)

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

• z =[у х] — матрица данных; первые столбцы соответствуют выхо­дам, последние — входам объекта;

• ny — количество выходов, по умолчанию nу =1;

• int — вектор-строка, задающий количество точек графика. По умолчанию выводятся все исходные данные;

• Т — интервал дискретизации (для задания правильного масштаба времени);

• рс — строковая переменная, определяющая вид представления входного сигнала (входных сигналов). При рс = 'рс' (значение по умолчанию) данный сигнал полагается кусочно-постоянным между моментами дискретизации, при рс = '11' — линейно изменяю­щимся.

Функция nyqplot

Функция nyqplot осуществляет построение годографа (диаграммы Найквиста):

nyqplot(g)

nyqplot([gl g2... gn])

nyqplot(g,sd,mode)

Аргументы g, gl, g2,..., gn, sd, mode — такие же, как у функции bodeplot.

Возвращается график АФХ (годограф).

Функция present

Функция present возвращает информацию о модели (представленной в так называемом тета-формате):

present(th)

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

Функция zpplot

Функция zpplot возвращает графическое изображение нулей и по­люсов моделей объекта:

zpplot(zepo)

zpplot(zpform(zepo1,zepo2,...,zepon))

zpplot (zepo, sd, mode, axis)

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

• zepo — матрица, содержащая информацию о нулях и полюсах модели объекта и получаемая с помощью функции th2zp (см. ниже);

• zpform(zepol,zepo2,.-, zepon) — функция объединения информа­ции о нулях и полюсах нескольких моделей одного и того же объ­екта с одними и теми же входами;

• при sd>0 к нулям и полюсам добавляются соответствующие дове­рительные области (по умолчанию sd=0);

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

• Когда в zepo содержится информация о нескольких входах, характер отображения определяется строковой переменной mode. При mode = 'sub' открываются одновременно графики, относящиеся ко всем входам; при mode = 'same' информация представляется в од­
ном графическом окне, при этом переход от изображения нулей и полюсов для одного входа к изображению нулей и полюсов для другого входа производится при нажатии клавиши Enter; mode = ’sep' — то же, но предыдущий график убирается до обработки сле­дующего входа. По умолчанию mode = 'sub';

• axis = [xl x2 yl y2] — задает масштабы осей. Задание axis = m — то же самое, что и axis = [-m m -m m].

Функция covf

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

R = covf(z,M)

R = covf(z,M,maxsize)

Аргументы:

• z— матрица данных размером Nxnz, каждый столбец которой co­ответствует входному или выходному сигналу (обычно z= [y u]);

• М — максимальная величина дискретного аргумента, для которой pассчитываются корреляционные функции, минус единица;

• maxsize — параметр, определяющий максимально допустимый размер матриц.

Возвращаемая величина — матрица R размером nz²xM с элементами

 

 

Функция сra

Функция сra определяет оценку ИХ методом корреляционного ана­лиза для одномерного (один вход — один выход) объекта:

cra(z);

[ir,R,cl] - cra(z,M,na,p1ot);

cra(R):

Аргументы:

• z — матрица экспериментальных данных вида z = [у и], где у —вектор-столбец, соответствующий выходным данным, u — вектор-столбец входных данных;

• М — максимальное значение дискретного аргумента, для которого производится расчет оценки ИХ, по умолчанию М=20;

• na — порядок модели авторегрессии (степени многочлена A(z)), которая используется для расчета параметров «обеляющего» фильт­ра Ф(г), по умолчанию nа=10. При nа=0 в качестве идентифици­рующего используется непреобразованный входной сигнал;

• plot — plot=0 означает отсутствие графика, plot =1 (по умолчанию) —график полученной оценки ИХ вместе с 99%-м доверительным коридором, plot =2 - выводятся графики всех корреляционных функций.

Возвращаемые величины: ir — оценка ИХ (вектор значений); cl — 99%-й доверительный коридор для оценки ИХ; R — матрица, эле­менты первого столбца которой — значения дискретного аргумента, элементы второго столбца — значения оценки автокорреляционной функции выходного сигнала (возможно, отфильтрованного), элемен­ты третьего столбца — значения оценки автокорреляционной функ­ции входного сигнала (возможно, «обеленного»), элементы четвер­того столбца — значения оценки взаимной корреляционной функции.

Функция spa

Функция 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 — вектор частот, для которых производится расчет частотных характеристик, по умолчанию [l:128]/128*pi/T;

• Т — интервал дискретизации;

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

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

• g — оценка W{e^jωТ)в частотном формате;

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

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

Рассмотрим пример. Пусть исходные данные содержатся в файле dryer2.mat. Воспользуемся функцией spa для нахождения оценок ам-плитудно- и фазочастотных характеристик объекта с выводом резуль-1лта в форме графиков.

»% Загрузка данных
» load dryer2
» z=[y2 u2];

» g - spa(z); %Оценива ние модели

» bodeplot(g) %Построение диаграммы Боде

Результат представлен на рисунке (оценка АЧХ построена в лога­рифмическом масштабе).

В продолжение примера ниже приведены функции, обеспечиваю­щие вывод графиков АЧХ, ФЧХ и S_v(ω) с доверительными коридорами шириной в три среднеквадратических отклонения

» w = logspace(-2,pi,l28);

» [g,phiv] = spa(z,[ ],w);

» % (пустая матрица означает значение по умолчанию )

» bodeplot ([g phiv], 3)

 

Функция etfe

Функция etfe так же, как предыдущая функция, возвращает оценку дискретной передаточной функции для обобщенной линейной моде­ли (см. выше) одномерного объекта в частотной форме. Рекомендуется к использованию для узкополосных объектов. Записывается в виде:

g = etfe(z)

g = etfe(z,M,N,T)

Аргументы:

• z — как для предыдущей функции;

• М — как для предыдущей функции;

• N — определяет диапазон частот для расчета (эта величина долж­на быть степенью 2) значений частотных характеристик, согласно формуле w = [1:N]/N*pi/T. По умолчанию N=128;

• Т — интервал дискретизации (по умолчанию Т=1).
Возвращаемая величина g = W(e^jωТ).

 

В продолжение предыдущего примера (с его исходными данными) проведем сравнение результатов использования функций spa и etfe. Соответствующая программа и графики приведены ниже.

» ge = etfe(z):

» gs - spa(z);

» bodeplot([ge gs])

 

Функция ar

 

Функция аг оценивает параметры модели авторегрессии (AR), то есть коэффициенты полинома А(z), при моделировании скалярных вре­менных последовательностей:

 

th - ar(y,n)

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

 

Аргументы:

• у — вектор-столбец данных, содержащий N элементов;

• n — порядок модели (число оцениваемых коэффициентов);

• аргумент approach (строковая переменная) определяет метод оце­нивания:

o 'fb'— прямой-обратный метод (разновидность метода наимень­ших квадратов). Используется по умолчанию;

о 'ls' — метод наименьших квадратов (МНК);

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

o 'burg' — метод Бэрга (комбинация МНК с минимизацией гар­монического среднего);

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

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

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

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

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

o window = 'pow* — последующие отсутствующие данные заменя­ются нулями, так что суммирование расширяется до момента времени N+n;

o window = 'ppw' — и начальные, и последующие отсутствующие данные заменяются нулями. Используется в алгоритме Юла— Уокера;

• Аргумент maxsize определяет максимальную размерность задачи;

• Т — интервал дискретизации.

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

• th — информация о модели в так называемом тета-формате (внут­реннем матричном формате представления параметрических моделей пакета System Identification);

• refl — информация о коэффициентах и функции потерь.

Пример использования функции аг для оценивания модели авторег­рессии 4-го порядка приведен ниже.

» load dryer2; %3aгрузка данных

» у = [у2(1:300)];

» th = ar(y,4); %Оценивание AR-модели

» present(th) %Вывод информации о модели

Функция armax

Функция armax оценивает параметры ARMAX-модели:

th = armax(z,nn)

th = armax(z,nn,'trace')

[th, iter_info] = armax(z,nn,maxiter,tol,lim,maxsize,T, 'trace')

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

• z = [у u] — матрица экспериментальных данных; в случае многомерного входа и представляет собой матрицу с числом столбцов, равным числу входов;

• nn = [na nb nc nk] — задаваемые параметры ARMAX-модели (сте­пени полиномов и величина задержки);

• 'trace' — наличие этого аргумента задает выдачу информации о ходе итерационного процесса вычислений (трассировка);

• maxiter — определяет максимальное число итераций при решении оптимизационной задачи, по умолчанию maxiter =10;

• tol — определяет заданную погрешность вычислений при поиске
минимума, по умолчанию tol = 0.01;

• lim — так называемый критерий робастности, по умолчанию lim = 1.6.
При lim = 0 используется обычный МНК;

• maxsize — определяет максимальную размерность задачи, по умол­
чанию

• Т — интервал дискретизации, но умолчанию равен 1. Возвращаемые величины:

• th — информация о модели в тета-формате;

• iterinfo = [номер последней итерации, значение функции потерь на данной итерации] — информация о выполнении итерационно­го процесса оценивания модели.

Использование функции агтах иллюстрируется следующим приме­ром:

» load dryer2 %3aгрузка данных

» y =[у2 u2];

» t = armax(z,[2 2 2 1]); %Оценивание ARMAX-иодели

» present(th) %Вывод информации о модели

 

This matrix was created by the command ARMAX on 1/19 2001 at 9:56

Loss fen: 0.0071844 Akaike's FPE: 0.0072712 Sampling interval 1

The polynomial coefficients and their standard deviations are

В =

0 0.0073 0.0451

0 0.0021 0.0025

A =

1.00 -1.6490 0.7026

0 0.0097 0.0088

C =

1.00 -0.2741 -0.4759

0. 0.0289 0.0285

В данном случае задана структура модели, все полиномы которой имеют 2-й порядок. Для выводимой информации, по-видимому, осо­бых разъяснений не требуется.

Функция агх

Функция агх оценивает параметры моделей AR и ARX:

th = arx(z,nn)

th = arx(z,nn,maxsize,T)

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

nn - [nа nb nк]

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

A{z)y{t) = B_{(z)u_{(t-nk_x) +

при этом

• в представлении z = [y u] экспериментальных данных u = [ul, u2,..., unu] есть матрица, столбцы которой соответствуют входам;

• nb и nk — векторы-столбцы, определяющие порядки полиномов B_j (z) и задержки nk_j для входов объекта.

Функция bj

Функция bj оценивает параметры модели Бокса—Дженкинса:

th = bj(z,nn)

th = bj(z,nn,'trace')

[th, iter_info] = bj(z,nn,maxiter,tol.,lim,maxsize,T,'trace')

Аргументы и возвращаемые величины идентичны приведенным для функции аrmах за тем исключением, что в данном случае аргумент nn содержит пять параметров: nn = [nb nc nd nf nk], то есть степени че­тырех полиномов и величину задержки.

Функция canstart

Функция canstart оценивает коэффициенты канонической формы модели для переменных состояния:

th = canstart(z,orders,nu)

th = canstart(z,orders,nu,dkx)

Аргументы функции следующие:

• z — матрица данных;

• orders — скаляр, задающий порядок п модели;

• nu— число входов;

• dkx = [d, k, x] — трехэлементный вектор, в котором:,

o k = l означает, что матрица К в модели полностью параметризо­вана (то есть подлежит оцениванию); к=0 означает, что К = 0;

o d имеет тот же смысл, но по отношению к матрице модели D;

o х=1 означает, что параметризован вектор начальных значений переменных состояния; при х = 0 данные значения полагаются нулевыми.

По умолчанию dkx = [0, 1, 1].

Возвращаемая величина th — модель в тета-формате.

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

» load dryer2 % 3агрузка данных

» z=[y2 u2];

» th = canstart(z,2,1); %Оценивание mодели

» present(th) %Вывод информации омодели

Функция ivar

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

A(z)y(t) = v(t)

где шум v(t) предполагается коррелированным случайным процес­сом типа процесса скользящего среднего, то есть имеющим представ­ление v(t) = H(z)e(t) (здесь e(t) — белый шум), с использованием квазиоптимального метода инструментальной переменной (так на­зываемого IV-метода). Функция записывается в виде:

th = ivar(y,na)

th = ivar(y,na,nc,maxsize,T)

где

• у — вектор данных;

• na — порядок модели авторегрессии (степень полинома A(z));

• nc - порядок полинома H(z), по умолчанию na=nc;

maxsize и Т имеют тот же смысл, что и для функций аrmах и bj.

Возвращаемая величина th — модель в тета-формате.

Функция ivx

Функция ivx оценивает параметры ARX-модели методом инструмен­тальной переменной:

th = ivx(z.nn.x)

th = ivx(z,nn,x.maxsize,T)

Функция аналогична функции агх и является ее альтернативой для случая, когда шум наблюдений является коррелированным. Аргумент х — матрица задаваемых инструментальных переменных. Остальные аргументы и возвращаемая величина аналогичны рассмотренным для функции arx.АМЕЧАНИЕ: Вместо функции агх рекомендуется использовать функцию iv4.

Функция iv4 оценивает параметры ARX-модели с использованием четырехступенчатого метода инструментальной переменной. Может применяться для идентификации многомерных объектов:

th = iv4(z,nn)

th = iv4(z,nn,maxsize,T)

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

Ниже приведен пример использования данной функции для иденти­фикации объекта с одним выходом и двумя входами при различных задержках по входам (0 и 2), но при одинаковых степенях полино­мов для данных входов (2):

» z = [у u1 u2]:

» nb = [2 2]; % Задание степеней полиномов Bi(z)

» nk = [0 2]; %Задание задержек по входам

» %Оценивание модели

» th = iv4(z,[2 nb nk]);

Функция n4sid

Функция n4sid используется для оценивания параметров моделей для переменных состояния в канонической форме при произвольном чис­ле входов и выходов:

h = n4sid(z)

[th,A0] = n4sid(z,rder,ny,auxord,dkx,maxsize,T,'trace')

где:

• z — матрица экспериментальных данных;

• аргумент order задает порядок модели. Если данный аргумент вводится как вектор-строка, например, order = [1:10], то предва­рительные расчеты выполняются по моделям всех заданных порядков (по умолчанию — от 1-го по 10-й) с выводом графика, позволяющсго выбрать оптимальный порядок. Если order = 'best' (по умолчанию), выбирается модель «наилучшего» порядка (пу­тем выбора наилучшей среди моделей 1-го— 10-го порядков);

• nу — количество выходов, по умолчанию nу = 1;

• auxord - так называемый дополнительный порядок, используемый алгоритмом oценивания. Должен быть больше, чем порядок, зада­ваемый параметром order (по умолчанию auxord = 1.2*order+3;
Если auxord вводится как вектор-строка (так же, как order), выби­рается модель наилучшего порядка;

• dkx, maxsize, T, 'trace' — как для функций bj, canstart, ivar.
Возвращаемые величины:

• th — имя (идентификатор) построенной модели в тета-формате;

• А0 — наилучший порядок модели (определяемый с помощью зада­ния auxord).

Функция ое

Функция ое оценивает параметры ОЕ-модели:

th = oe(z,nn)

th = oe(z,nn,'trace')

[th, iter_info] = oe(z,nn,maxiter,tol,lim,maxsize,T, 'trace')

Здесь аргумент nn задается в форме nn = [nb nf nk] и содержит ин­формацию о степенях полиномов модели В(г), F{z) и о величине за­держки nk.

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

Функция pem

Функция реm оценивает параметры обобщенной многомерной линей­ной модели вида

 

 

Функция записывается в виде:

th = pem(z,nn)

th = pem(z,nn, 'trace')

[th, iter_info] = pem(z,nn,index,maxiter,tol,lim,maxsize,T,'trace')

где:

• аргумент nn имеет представление nn = [na nb nc nd nf nk], где па, nb,
nc, nd, nf, nk в многомерном случае являются векторами-столбцами;

• index — вектор-строка, содержащий индексы коэффициентов, ко­торые должны быть оценены (по умолчанию оцениваются все ко­эффициенты).

Остальные аргументы и возвращаемые величины соответствуют рас­смотренным выше.

Функция rarmax

Функция rarmax реализует процедуру рекуррентного оценивания ко­эффициентов ARMAX-модели:

thm = rarmax(z,nn,adm,adg)

[thm,yhat,P,phi,psi] = rarmax(z,nn,adm,adg,th0,P0,phi0,psi0)

Здесь:

• z — матрица данных;

• nn — аргумент, задающий параметры ARMAX-модели, nn = [na nb nc nk];

• adm и adg — аргументы, задающие вид процедуры оценивания, на­пример, значения adm ='ff' и adg = lam задают рекуррентный ме­тод наименьших квадратов, adm = 'ug' и adg = gam — градиентный метод оценивания и т. п. (дополнительную информацию можно
получить, используя команду help rarmax);

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

• Р0 — начальная матрица ковариаций оцениваемых коэффициен­тов, по умолчанию равна единичной матрице, умноженной на 10⁴;

• phi0 и psi0 — соответственно, начальные значения вектора данных и вектора-градиента (по умолчанию — нулевые). Их размеры дос­таточно сложным образом зависят от параметров модели.

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

• thm — матрица оценок коэффициентов модели; k-я строка матрицы содержит оценки, соответствующие моменту времени k в сле­дующем порядке:

thm(k,:) = [al,a2,...ana,bl,...,bnb];

• yhat — вектор прогнозируемых значений выхода (соответствующий вектору z(k,l));

• Р, phi, psi — конечные значения матрицы ковариаций оценок коэффициентов, вектора данных и вектора-градиента.

Функция гагх

Функция гагх реализует процедуру рекуррентного оценивания коэф­фициентов AR и ARX-моделей:

thm = rarx(z,nn,adm,adg)

[thm,yhat,P,phi] = rarx(z,nn,adm,adg,h0,P0,phi0)

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

Функция rbj

Функция rbj — то же, что предыдущая функция, но для модели Бокса-Дженкинса:

thm = rbj(z,nn,adm,adg)

[thm,yhat,P,phi,psi] =...

rbj(z,nn,adm,adg,th0,P0,phi0,psi0)

Аргументы и возвращаемые величины — как у функции гагтах, с уче­том того, что оценивается модель Бокса—Дженкинса.

Функция rое

Функций rое — то же, что предыдущая функция, но для ОЕ-модели:

thm = roe(z,nn,adm,adg)

[thm,yhat,P,phi,psi] = roe(z,nn,adm,adg,th0,P0,phi0,psi0)

Аргументы и возвращаемые величины те же, что у функции rаrmах.

Функция rреm

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

thm = rpem(z,nn,ddm,adg)

[thm,yhat,P,phi,psi] = rpem(z,nn,adm,adg,th0,P0,phi0,psi0)

Аргументы и возвращаемые величины те же, что у функции rarmax.

Функция rplr

Функция rplr имеет то же назначение, что и предыдущая функция. Отличие состоит в алгоритме оценивания. Функция записывается в виде:

thm = rplr(z,nn,adm,adg)

[thm,yhat,P.phi] = rplr(z,nn,adm,adg,th0,P0,phi0)

Функция segment

Функция segment осуществляет разбиение (сегментацию) данных и оценивание коэффициентов моделей для каждого сегмента в пред­положении, что в пределах каждого сегмента коэффициенты модели являются постоянными, но могут изменяться (скачкообразно) при переходе от сегмента к сегменту. Оцениваемые модели — AR, ARX, ARMAX. Функция записывается в виде:

segm = segment(z,nn)

[segm,V,thm,R2e] = segment(z,nn,R2,q,Rl,M,th0,P0,ll,mu)

Здесь:

• z — матрица данных;

• аргумент пп определяет порядок модели:
o для ARMAX-модели nn = [nа nb nc nk];
o для ARX-модели nn = [nа nb nк];

o для AR-модели nn = nа;

• R2 — оценка дисперсии шума наблюдений;

• q — вероятность скачкообразного изменения модели в момент вре­мени k_y по умолчанию 0.01;

• R1 — предполагаемая матрица ковариаций оценок коэффициентов при их возможном скачкообразном изменении, по умолчанию —единичная матрица;

• М — число параллельных моделей, используемых в алгоритме, по умолчанию 5;

• th0 и Р0 — смысл данных аргументов рассмотрен выше, при описа­нии функции rаrmах (здесь начальное значение Р0 — единичная мат­рица, умноженная на 10);

• l l — аргумент, определяющий существование параллельных моделей (каждая модель не исключается, по крайней мере, в течение 11 итераций), по умолчанию l l = 1;

• mu — параметр алгоритма оценивания матрицы R2, но умолчанию 0.97.
Возвращаемые величины:

• segm — матрица, k-я строка которой содержит оценки коэффици­ентов, соответствующие моменту времени k (как возвращаемая матрица в функциях rаrх и rаrmах);

• thm — матрица, содержащая оценки коэффициентов, соответствую­щие моменту времени k без учета сегментации данных;

• V — сумма квадратов ошибок сегментированной модели;

• R2e — полученная оценка матрицы R2 для момента времени k.

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

y(t) = b₁* 1

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

» у = sin([l:50]/3)'; %Задание гармонического сигнала

» thm = segment([y,ones(size(y))],[0 1 1],0,1); %Сегментация данных

» plot ([thm,у]) %Вывод результатов

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

 

В эту группу входят следующие 6 функций.

Функция arx2th

Функция arx2th конструирует матрицу модели так называемого тета-формата по полиномам ARX-модели многомерного объекта:

th = arx2th(A,B,ny,nu)

th = arx2th(A,B,ny,nu,lam,T)

Исходная ARX-модель многомерного объекта имеет вид:

y(t) + A₁y(t-1) + A₂y(t-2) +... +A_nay(t-na) = В₀ и(t) + B₁ u(t -1) +...+ B_nb u(t - nb) + e(t),

где y(t) — вектор-столбец с nу элементами, u(t) — вектор-столбец с nu элементами (то есть объект имеет nu входов и nу выходов), А_к и В_к — матрицы размеров nyxny и nuxnu соответственно, e(t) — вектор шума.

Аргументы функции следующие:

• А = [I А₁ А₂... А_na], где I — единичная матрица;

• В = [В_о В₁... В_nЬ];

• ny и nu — количество выходов и входов модели (объекта);

• lam — матрица ковариаций шума e(t), no умолчанию — единичная матрица;

• Т — интервал дискретизации.

Функция canform

Функция canform конструирует каноническую форму модели для пе­ременных состояния многомерного объекта:

ms = canform(orders,nu)

ms = canform(orders,nu,dkx)

Используется совместно с функцией mf2th (см. ниже). Аргументы ана­логичны аргументам описанной выше функции canstart. Возвращае­мая величина ms — матрица, отражающая информацию о структуре модели.

Функция mf2th

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

th = mf2th(model,cd,parval)

th = mf2th(model,cd,parval,aux,lambda,T)

где:

• model — имя m-файла, содержащего информацию о том, как матрицы модели зависят от свободных (оцениваемых) параметров;

• cd — строковая переменная, задающая тип модели:
o cd='c' — модель для непрерывного времени,

o cd='d' — модель для дискретного времени;

• parval — вектор-строка, содержащий номинальные значения пара­метров;

• Т — интервал дискретизации;

• aux — матрица вспомогательных переменных, используемых для различных целей и задаваемых в файле model;

• lambda — интенсивность шума.
Возвращаемая величина th — модель в тета-формате.

 

Функция modstruc

Функция modstruc определяет структуру модели для переменных со­стояния с частично незаданными коэффициентами матриц модели:

ms = modstruc(A,B,C,D,K)

ms = modstruc(A,B,C,D,K,x0)

Аргументы А, В, С, D, К, х0 — соответственно матрицы модели и вектор начальных условий (по умолчанию — нулевой). Известные коэффи­циенты матриц задаются числами, неизвестные (свободные), подле­жащие оцениванию, — как NaN. Возвращаемая матрица ms содержит информацию о структуре модели для дальнейшего использования с помощью функции ms2th.

В следующем примере системы 2-го порядка неизвестные элементы содержатся в матрицах А и В:

» A = [NaN,0;0,NaN]; % Задание матрицы А

» B = [NaN;NaN]; % Задание матрицы В

» О = [1,1]: % Задание матрицы С

» D = 0; % Задание матрицы D

» К = [0;0]; % Задание матрицы К

» ms = modstruc(A,B,C,D,K) % Формирование матрицы структуры модели

 

ms =

NaN 0 NaN 1 0 0 0 1

0 NaN NaN 1 0 0 0 2

Функция ms2th

Функция ms2th преобразует матрицу структуры в модель тета-фор-мата:

th = ms2th(ms)

th = ms2th(ms,cd,parval,lambda,T)

Аргументы и возвращаемая величина соответствуют описанным для предыдущих функций данной группы. Для примера определим сле­дующую структуру модели (для непрерывного времени):

. Г⁰» °1 Г^ез1 [о ej [e₄j

у =[1 \]х + е

с начальными значениями коэффициентов θ = [-0.2 -0.3 2 4] и оце­ним соответствующие свободные параметры:

» load dryer2 %Загрузка данных

» z = [у2 и2];

» %Задание структур матриц модели

» A = [NaN,0;0,NaN];

» B = [NaN;NaN];

» С = [1,1];

» ms=modstruc(A,B,C,0,[0;0]);

» th=ms2th(ms,'c',[-0.2,-0.3,2,4]);

» th=pem(z,th):

» present(th) %Отображение информации о модели

This matrix was created by the command РЕМ on 1/20 2001 at 12:22

Loss fen: 0.17279 Akaike's FPE: 0.17418 Continuous time model estimated using

sampling interval 1

The state-space matrices with standard deviations given as imaginary parts are

a =

-0.2585 + 0.2844i 0

0 - 0.2816 + 0.33231 i

 

b =

3.5019 + 93.6630i

-3.5374 + 93.6538i

 

 

c = 1 1

d =

k =

 

x0 =

lambda =

0.1728

Функция poly2th

Функция poly2thконструирует модель тета-формата из исходной мо­дели «вход-выход»:

th = poly2th(A,B)

th = poly2th(A,B,C,D,F,lam,T)

Исходная модель такая же, как и для рассмотренной функции реm (см. выше). Аргументы функции следующие:

• А, В, С, D, F — векторы-строки, содержащие коэффициенты соответ­ствующих полиномов;

• Т — интервал дискретизации;

• 1am — дисперсия шума модели.

Возвращаемая величина th — модель в тета-формате.

 

Функция fixpar

Функция fixpar фиксирует коэффициенты в структурах тета-моделей для переменных состояния и ARX, созданных с помощью функ­ций ms2th и arx2th:

thn = fixpar(tho,matrix)

thn = fixpar(tho,matrix,elements,parval)

где:

• tho — исходная тета-модель;

• аргумент matrix (строковая переменная):

o если исходная модель задана через переменные состояния, при­нимает одно из значений: 'А', 'В', 'С', 'D', 'К' или 'х0';

o если исходная модель — ARX-типа, то параметр matrix прини­мает одно из знач