Документация

Системное описание

Это тематическое исследование касается данных, собранных от лабораторных весов "фен". (Преподаватель Процесса обратной связи PT326; См. также страницу 525 в Ljung, 1999). Процесс работает можно следующим образом: Воздух вентилируется через трубу и нагревается во входе. Температура воздуха измеряется термопарой при выходе. Вход является напряжением по нагревающемуся устройству, которое является только сеткой проводов резистора. Выход является температурой воздуха выхода, представленной измеренным напряжением термопары.

Подготовка Данных для Анализа

Сначала мы загружаем данные ввода - вывода к MATLAB® Workspace.

load dryer2;

Векторный y2, выход, содержит 1 000 измерений напряжения термопары, которое пропорционально температуре в воздушном потоке выхода. Векторный u2 содержит 1 000 точек входных данных, состоящих из напряжения, применился к нагревателю. Вход был сгенерирован как бинарная случайная последовательность, которая переключается от одного уровня до другого с вероятностью 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: Импульсная характеристика модели FIR, оцененной с помощью 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;

Смотрите "Объекты данных и Объекты модели в System Identification Toolbox" пример для получения дополнительной информации относительно объектов модели. Чтобы узнать, какие свойства объекта модели могут быть получены, используйте 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
           Report: [1x1 idresults.ssest]
       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]

Чтобы выбрать значения матриц пространства состояний и их 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+15 *

    0.1470    0.0957    0.0732
    0.2935    0.1998    0.2061
    1.2756    0.6196    0.2394


dB =

   1.0e+13 *

    0.3388
    1.1938
    4.0974


dC =

   1.0e+14 *

    2.2060    1.5977    0.3666


dD =

     0


dK =

   1.0e+13 *

    1.3862
    4.9801
    7.2037

Неопределенность является довольно большой даже при том, что модель соответствует данным об оценке хорошо. Это вызвано тем, что модель сверхпараметризована, то есть, она имеет более свободные параметры, чем, что могло быть однозначно определено. Отклонение параметров в таких случаях не четко определено. Однако это не подразумевает, что модель ненадежна. Мы можем построить время - и частотная характеристика этого графика и просмотреть отклонение как области уверенности, как обсуждено затем.

Анализ предполагаемой модели

Диаграмма Боде сгенерированной модели может быть получена с помощью 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 довольно надежно.

Оценка моделей с предписанной структурой

System Identification Toolbox может также использоваться, чтобы получить модель с предписанной структурой. Например, модель разностного уравнения с 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 = 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: Предвещайте ответы m1, m2 и m3 сравненный с непараметрической спектральной моделью gs оценки.

Частотные характеристики из этих трех моделей/методов очень близки. Это указывает, что этот ответ надежен.

Кроме того, годограф Найквиста может анализироваться с областями неопределенности, отмеченными на определенных частотах:

showConfidence(nyquistplot(m1,'b',m2,'r',m3,'c',gs,'g'),3)

Рисунок 13: годографы Найквиста моделей m1, m2, m3 и gs.

Непараметрическая модель gs показывает большую часть неопределенности в ответ.

Дополнительная информация

Для получения дополнительной информации об идентификации динамических систем с System Identification Toolbox посещают страницу информации о продукте System Identification Toolbox.