В этом примере показано, как разрабатывать и анализировать простые модели на основе реальных данных лабораторных процессов. Начнем с небольшого описания процесса, научимся импортировать данные на панель инструментов и предварительно обрабатывать/кондиционировать их, а затем систематически приступаем к оценке параметрических и непараметрических моделей. После определения моделей мы сравниваем оцененные модели, а также проверяем модель с фактическими выходными данными эксперимента.
Это тематическое исследование касается данных, собранных из лабораторной шкалы «фен». (Инструктор по процессам обратной связи PT326; См. также страницу 525 в Ljung, 1999). Процесс работает следующим образом: Воздух раздувается через трубку и нагревается на входе. Температура воздуха измеряется термопарой на выходе. Вход - это напряжение на нагревательном устройстве, которое представляет собой только сетку резисторных проводов. Выходной сигнал представляет собой температуру воздуха на выходе, представляемую измеренным напряжением термопары.
Сначала загружаем данные ввода-вывода в рабочую область MATLAB ®.
load dryer2;
Вектор y2выходной сигнал содержит 1000 измерений напряжения термопары, которое пропорционально температуре в выходном воздушном потоке. Вектор u2 содержит 1000 входных точек данных, состоящих из напряжения, подаваемого на нагреватель. Вход был сформирован как двоичная случайная последовательность, которая переключается с одного уровня на другой с вероятностью 0,2. Время выборки составляет 0,08 секунды.
Следующим шагом является настройка данных как объекта iddata.
dry = iddata(y2,u2,0.08);
Чтобы получить информацию о данных, просто введите имя iddata объект в окне команды MATLAB:
dry
dry =
Time domain data set with 1000 samples.
Sample time: 0.08 seconds
Outputs Unit (if specified)
y1
Inputs Unit (if specified)
u1
Для проверки свойств указанного выше объекта iddata используйте get команда:
get(dry)
ans =
struct with fields:
Domain: 'Time'
Name: ''
OutputData: [1000x1 double]
y: 'Same as OutputData'
OutputName: {'y1'}
OutputUnit: {''}
InputData: [1000x1 double]
u: 'Same as InputData'
InputName: {'u1'}
InputUnit: {''}
Period: Inf
InterSample: 'zoh'
Ts: 0.0800
Tstart: []
SamplingInstants: [1000x0 double]
TimeUnit: 'seconds'
ExperimentName: 'Exp1'
Notes: {}
UserData: []
Для лучшего ведения книги рекомендуется называть каналы ввода и вывода и единицы времени. Эти имена будут распространяться в ходе анализа этого объекта iddata:
dry.InputName = 'Heater Voltage'; dry.OutputName = 'Thermocouple Voltage'; dry.TimeUnit = 'seconds'; dry.InputUnit = 'V'; dry.OutputUnit = 'V';
Теперь, когда набор данных готов, мы выбираем первые 300 точек данных для оценки модели.
ze = dry(1:300)
ze =
Time domain data set with 300 samples.
Sample time: 0.08 seconds
Outputs Unit (if specified)
Thermocouple Voltage V
Inputs Unit (if specified)
Heater Voltage V
Постройте график интервала от 200 до 300:
plot(ze(200:300));
Рис. 1: Снимок измеренных данных фена.
Из вышеприведенного графика можно видеть, что данные не являются нулевыми средними. Давайте уберем постоянные уровни и сделаем данные нулевыми.
ze = detrend(ze);
Тот же набор данных после его уничтожения:
plot(ze(200:300)) %show samples from 200 to 300 of detrended data
Рис. 2: Детерминированные оценочные данные.
Теперь, когда набор данных был уменьшен и нет явных отклонений, давайте сначала оценим импульсную характеристику системы с помощью корреляционного анализа, чтобы получить некоторое представление о временных константах и т.п.:
clf mi = impulseest(ze); % non-parametric (FIR) model showConfidence(impulseplot(mi),3); %impulse response with 3 standard %deviations confidence region
Рис. 3: Импульсная характеристика модели КИХ, оцененная с использованием ze.
Затененная область отмечает 99,7% доверительный интервал. Имеется временная задержка (мертвое время) из 3 выборок, прежде чем выходной сигнал реагирует на входной сигнал (значительный выходной сигнал вне доверительного интервала).
Простейшим способом начала работы над программой параметрической оценки является построение модели состояния-пространства, в которой порядок модели определяется автоматически, с использованием метода ошибки прогнозирования. Давайте оценим модель, используя ssest методика оценки:
m1 = ssest(ze);
m1 - модель состояния-пространства, идентифицированная непрерывно и представленная символом idss объект. Алгоритм оценки выбирает 3 в качестве оптимального порядка модели. Чтобы проверить свойства расчетной модели, просто введите имя модели в окне команды:
m1
m1 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 -0.4839 -2.011 2.092
x2 3.321 -1.913 5.998
x3 1.623 -17.01 -15.61
B =
Heater Volta
x1 -0.05753
x2 0.02004
x3 1.377
C =
x1 x2 x3
Thermocouple -14.07 0.07729 0.04252
D =
Heater Volta
Thermocouple 0
K =
Thermocouple
x1 -0.9457
x2 -0.02097
x3 2.102
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 18
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "ze".
Fit to estimation data: 95.32% (prediction focus)
FPE: 0.001621, MSE: 0.001526
Отображение предполагает, что модель является произвольной (все записи матриц A, B и C рассматривались как свободные параметры) и что оценочная модель вполне подходит для данных (более 90% подходят). Чтобы извлечь свойства этой модели, например, чтобы получить A матрицу генерируемого выше дискретного состояния-космического объекта можно использовать оператор точки:
A = m1.a;
Дополнительные сведения об объектах модели см. в примере «Данные и объекты модели в панели инструментов идентификации системы». Чтобы узнать, какие свойства объекта модели можно извлечь, используйте get команда:
get(m1)
A: [3x3 double]
B: [3x1 double]
C: [-14.0706 0.0773 0.0425]
D: 0
K: [3x1 double]
StateName: {3x1 cell}
StateUnit: {3x1 cell}
Structure: [1x1 pmodel.ss]
NoiseVariance: 1.2587e-04
InputDelay: 0
OutputDelay: 0
Ts: 0
TimeUnit: 'seconds'
InputName: {'Heater Voltage'}
InputUnit: {'V'}
InputGroup: [1x1 struct]
OutputName: {'Thermocouple Voltage'}
OutputUnit: {'V'}
OutputGroup: [1x1 struct]
Notes: [0x1 string]
UserData: []
Name: ''
SamplingGrid: [1x1 struct]
Report: [1x1 idresults.ssest]
Чтобы получить значения матриц состояния-пространства и их неопределенностей 1 стандартного отклонения, используйте idssdata команда:
[A,B,C,D,K,~,dA,dB,dC,dD,dK] = idssdata(m1)
A =
-0.4839 -2.0112 2.0917
3.3205 -1.9135 5.9981
1.6235 -17.0096 -15.6070
B =
-0.0575
0.0200
1.3770
C =
-14.0706 0.0773 0.0425
D =
0
K =
-0.9457
-0.0210
2.1019
dA =
1.0e+14 *
1.1413 1.1945 0.8367
2.0899 1.4613 1.2117
9.8253 7.1186 1.8997
dB =
1.0e+13 *
0.4013
1.5604
4.5776
dC =
1.0e+14 *
2.1049 1.3852 0.4309
dD =
0
dK =
1.0e+13 *
1.3260
4.5015
7.0925
Неопределенности довольно велики, даже несмотря на то, что модель хорошо вписывается в оценочные данные. Это происходит потому, что модель чрезмерно параметризована, то есть имеет больше свободных параметров, чем то, что может быть однозначно идентифицировано. Дисперсия параметров в таких случаях не вполне определена. Однако это не означает, что модель ненадежна. Мы можем построить график временной и частотной характеристики этого графика и рассматривать дисперсию как доверительные области, как описано ниже.
График Бода генерируемой модели может быть получен с помощью bode как показано ниже:
h = bodeplot(m1);
Рисунок 4: График построения расчетной модели.
Щелкните правой кнопкой мыши на графике и выберите «Характеристика» - > «Область доверия». Или используйте showConfidence для просмотра расхождения ответа.
showConfidence(h,3) % 3 standard deviation (99.7%) confidence region
Рисунок 5: График Бода с 3 областью достоверности стандартного отклонения.
Аналогично, мы можем генерировать график шага и связанную с ним 3 доверительную область стандартного отклонения. Мы можем сравнить ответы и связанные отклонения параметрической модели m1 с непараметрической моделью mi:
showConfidence(stepplot(m1,'b',mi,'r',3),3)
Рисунок 6: Пошаговый график моделей m1 и mi с доверительными областями.
Мы также можем рассмотреть график Найквиста и отметить области неопределенности на определенных частотах эллипсами, соответствующими 3 стандартным отклонениям:
Opt = nyquistoptions;
Opt.ShowFullContour = 'off';
showConfidence(nyquistplot(m1,Opt),3)
Рисунок 7: График Найквиста оценочной модели, показывающий области неопределенности на определенных частотах.
Графики ответов показывают, что расчетная модель m1 является достаточно надежным.
Панель инструментов идентификации системы может также использоваться для получения модели с заданной структурой. Например, модель уравнения разности с 2 полюсами, 1 нулем и 3 задержками выборки может быть получена с помощью arx как показано ниже:
m2 = arx(ze,[2 2 3]);
Для просмотра модели введите имя модели в окне команд.
m2
m2 =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)
A(z) = 1 - 1.274 z^-1 + 0.3942 z^-2
B(z) = 0.06679 z^-3 + 0.04429 z^-4
Sample time: 0.08 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nk=3
Number of free coefficients: 4
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARX on time domain data "ze".
Fit to estimation data: 95.08% (prediction focus)
FPE: 0.001756, MSE: 0.001687
Функция непрерывной временной передачи с 2 полюсами, одной нулевой и 0,2 секундной задержкой передачи может быть оценена с использованием tfest команда:
m3 = tfest(ze, 2, 1, 0.2)
m3 =
From input "Heater Voltage" to output "Thermocouple Voltage":
1.183 s + 26.55
exp(-0.2*s) * ---------------------
s^2 + 11.61 s + 28.63
Continuous-time identified transfer function.
Parameterization:
Number of poles: 2 Number of zeros: 1
Number of free coefficients: 4
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using TFEST on time domain data "ze".
Fit to estimation data: 88.79%
FPE: 0.009126, MSE: 0.008768
Насколько хорошая модель? Один из способов выяснить это - смоделировать его и сравнить выходные данные модели с измеренными выходными данными. Выберите часть исходных данных, которая не использовалась при построении модели, например, от образцов 800 до 900. После предварительной обработки данных проверки мы используем compare как показано ниже, для просмотра качества прогноза:
zv = dry(800:900); % select an independent data set for validation zv = detrend(zv); % preprocess the validation data set(gcf,'DefaultLegendLocation','best') compare(zv,m1); % perform comparison of simulated output
Рис. 8: Смоделированный отклик модели в сравнении с выводом данных проверки.
Здесь можно заметить, что соглашение очень хорошее. Показанное значение «Fit» рассчитывается следующим образом:
Fit = 100*(1 - norm(yh - y)/norm(y-mean(y)))
где y - измеренный выходной сигнал (= | zv.y |), иyh - вывод модели m1.
Чтобы сравнить производительность моделей, которые мы оценили, например m1, m2 и m3 с данными проверки zv, мы можем снова использовать compare команда:
compare(zv,m1,'b',m2,'r',m3,'c');
Рис. 9: Сравнение ответов моделей m1, m2, m3 о наборе данных проверки ze.
Графики полюсов и нулей для моделей можно получить с помощью iopzplot:
h = iopzplot(m1,'b',m2,'r',m3,'c');
Рис. 10: Полюса и нули моделей m1, m2 и m3.
Неопределенности в полюсах и нулях также могут быть получены. В следующем операторе '3' относится к числу стандартных отклонений.
showConfidence(h,3);
Рисунок 11: Карта полюсов и нулей с областями неопределенности.
Приведенные выше частотные функции, полученные из моделей, можно сравнить с той, которая получена непараметрическим методом спектрального анализа (spa):
gs = spa(ze);
spa создает модель IDFRD. bode снова может использоваться для сравнения с передаточными функциями полученных моделей.
w = linspace(0.4,pi/m2.Ts,200); opt = bodeoptions; opt.PhaseMatching = 'on'; bodeplot(m1,'b',m2,'r',m3,'c',gs,'g',w,opt); legend('m1','m2','m3','gs')
Рис. 12: Ответы Bode m1, m2 и m3 сравнение с непараметрической моделью спектральной оценки gs.
Частотные характеристики трех моделей/методов очень близки. Это означает, что этот ответ надежен.
Также график Найквиста может быть проанализирован с областями неопределенности, отмеченными на определенных частотах:
showConfidence(nyquistplot(m1,'b',m2,'r',m3,'c',gs,'g'),3)
Рисунок 13: Найквистские графики моделей m1, m2, m3 и gs.
Непараметрическая модель gs проявляет наибольшую неопределенность в ответ.