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

Цели

Оцените и проверьте линейные модели из данных multi-input/single-output (MISO), чтобы найти ту, которая лучше всего описывает динамику системы.

После выполнения этого руководства вы сможете выполнить следующие задачи с помощью командной строки:

Описание данных

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

В первом эксперименте оператор ввел импульсную волну на оба входа. Во втором эксперименте оператор ввел импульсную волну на первый вход и шаговый сигнал на второй вход.

Загрузка данных в рабочее пространство MATLAB

Загрузите данные.

load co2data

Эта команда загружает следующие пять переменных в Рабочее пространство MATLAB:

Для обоих экспериментов входные данные состоят из двух столбцов значений. Первый столбец значений - темп потребления химикатов (в килограммах в минуту), а второй столбец значений - ток (в амперах). Выходными данными является один столбец скорости производства углекислого газа (в миллиграммах в минуту).

Графическое изображение входных/выходных данных

Постройте график входных и выходных данных из обоих экспериментов.

plot(Time,Input_exp1,Time,Output_exp1)
legend('Input 1','Input 2','Output 1')
figure
plot(Time,Input_exp2,Time,Output_exp2)
legend('Input 1','Input 2','Output 1')

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

Второй график показывает второй эксперимент, где оператор применяет импульсную волну к первому входу и шаговый сигнал ко второму входу.

Удаление равновесных значений из данных

Построение графика данных, как описано в Графическом изображении Входных/Выходных Данных, показывает, что входы и выходы имеют ненулевые равновесные значения. В этом фрагменте руководства вы вычитаете равновесные значения из данных.

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

Увеличьте изображение графиков, чтобы увидеть, что самый ранний момент, когда оператор применяет нарушение порядка к входам, происходит после 25 минут установившихся условий (или после первых 50 выборок). Таким образом, среднее значение первых 50 выборок представляет условия равновесия.

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

Input_exp1 = Input_exp1-...
   ones(size(Input_exp1,1),1)*mean(Input_exp1(1:50,:));
Output_exp1 = Output_exp1-...
   mean(Output_exp1(1:50,:));
Input_exp2 = Input_exp2-...
   ones(size(Input_exp2,1),1)*mean(Input_exp2(1:50,:));
Output_exp2 = Output_exp2-...
   mean(Output_exp2(1:50,:));

Использование объектов для представления данных для системы идентификации

Объекты System Identification Toolbox™ данных, iddata и idfrd, инкапсулируйте значения данных и свойства данных в одну сущность. Можно использовать команды System Identification Toolbox, чтобы удобно манипулировать этими объектами данных как отдельными сущностями.

В этом фрагменте руководства вы создаете два iddata объекты, по одному для каждого из двух экспериментов. Вы используете данные Эксперимента 1 для оценки модели и данные Эксперимента 2 для валидации модели. Вы работаете с двумя независимыми наборами данных, потому что вы используете один набор данных для оценки модели, а другой для валидации модели.

Создание объектов iddata

Вы, должно быть, уже загрузили выборочные данные в MATLAB^® рабочая область, как описано в разделе Загрузка данных в рабочее пространство MATLAB.

Используйте эти команды для создания двух объектов данных, ze и zv :

Ts = 0.5; % Sample time is 0.5 min
ze = iddata(Output_exp1,Input_exp1,Ts);
zv = iddata(Output_exp2,Input_exp2,Ts);

ze содержит данные эксперимента 1 и zv содержит данные эксперимента 2. Ts - шаг расчета.

iddata конструктор требует три аргумента для данных временной области и имеет следующий синтаксис:

data_obj = iddata(output,input,sampling_interval);

Просмотр свойств iddata объект, использовать get команда. Для примера введите эту команду, чтобы получить свойства данных оценки:

get(ze)

ans = 

  struct with fields:

              Domain: 'Time'
                Name: ''
          OutputData: [2001x1 double]
                   y: 'Same as OutputData'
          OutputName: {'y1'}
          OutputUnit: {''}
           InputData: [2001x2 double]
                   u: 'Same as InputData'
           InputName: {2x1 cell}
           InputUnit: {2x1 cell}
              Period: [2x1 double]
         InterSample: {2x1 cell}
                  Ts: 0.5000
              Tstart: []
    SamplingInstants: [2001x0 double]
            TimeUnit: 'seconds'
      ExperimentName: 'Exp1'
               Notes: {}
            UserData: []

Чтобы узнать больше о свойствах этого объекта данных, смотрите iddata страница с описанием.

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

% Set time units to minutes
ze.TimeUnit = 'min';
% Set names of input channels
ze.InputName = {'ConsumptionRate','Current'};
% Set units for input variables
ze.InputUnit = {'kg/min','A'};
% Set name of output channel
ze.OutputName = 'Production';
% Set unit of output channel
ze.OutputUnit = 'mg/min';

% Set validation data properties
zv.TimeUnit = 'min';
zv.InputName = {'ConsumptionRate','Current'};
zv.InputUnit = {'kg/min','A'};
zv.OutputName = 'Production';
zv.OutputUnit = 'mg/min';

Можно проверить, что InputName свойство ze изменяется или «index» в это свойство:

ze.inputname;

Графическое изображение данных в объекте данных

Можно построить график iddata объекты, использующие plot команда.

Постройте график данных оценки.

plot(ze)

Нижние оси показывают входы ConsumptionRate и Current, и верхние оси показывают выходные ProductionRate .

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

figure   % Open a new MATLAB Figure window
plot(zv) % Plot the validation data

Увеличьте изображение графиков, чтобы увидеть, что процесс эксперимента усиливает первый вход (ConsumptionRate ) в 2 раза и усиливает второй вход (Current ) в 10 раз.

Выбор подмножества данных

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

Ze1 = ze(1:1000);
Zv1 = zv(1:1000);

Для получения дополнительной информации об индексации в iddata объекты соответствующую страницу с описанием.

Зачем оценивать модели степенной и частотной характеристики?

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

В этом фрагменте руководства вы оцениваете эти модели с помощью набора данных ze. Вы, должно быть, уже создали ze, как описано в разделе Создание объектов iddata.

Графики отклика из этих моделей показывают следующие характеристики системы:

Оценка частотной характеристики

Продукт System Identification Toolbox обеспечивает три функции для оценки частотной характеристики:

Использование spa для оценки частотной характеристики.

Ge = spa(ze);

Постройте график частотной характеристики как диаграмма Боде.

bode(Ge)

Амплитудные достигать максимума на частоте 0,54 рад/мин, что предполагает возможное резонансное поведение (комплексные полюсы) для первой комбинации вход-выход - ConsumptionRate на Production .

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

Оценка эмпирической переходной характеристики

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

% model estimation
Mimp = impulseest(Ze1,60);

% step response
step(Mimp)

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

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

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

Зачем оценивать задержки?

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

Если вы не знаете задержек ввода/вывода для вашей системы из эксперимента, можно использовать программное обеспечение System Identification Toolbox, чтобы оценить задержку.

Оценка задержек с использованием структуры модели ARX

В случае систем с одним входом можно считать задержку на графике импульсной характеристики. Однако в случае систем с несколькими входами, таких как в этом руководстве, вы можете не сказать, какой вход вызвал начальное изменение выхода, и вы можете использовать delayest вместо этого команда.

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

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

y (t) представляет выход в момент t, u (t) представляет вход в момент t, na - количество полюсов, nb - количество b параметров (равное количеству нулей плюс 1), nk - количество выборок перед тем, как вход влияет на выход системы (называется задержкой или потерей времени модели), и e (t) является нарушением порядка белого шума.

delayest принимает, что _na = nb = 2 и что шум e белый или незначительный, и оценивает _nk.

Чтобы оценить задержку в этой системе, введите:

delayest(ze)

ans =

     5    10

Этот результат включает два числа, потому что существует два входа: предполагаемая задержка для первого входа 5 выборки данных и предполагаемая задержка для второго входа 10 выборки данных. Потому что шаг расчета для экспериментов 0.5 min, это соответствует 2.5 -min задержка перед тем, как первый вход влияет на выход, и 5.0 -min задержка перед тем, как второй вход влияет на выход.

Оценка задержек с помощью альтернативных методов

Существует два альтернативных метода оценки временной задержки в системе:

Зачем оценивать порядок модели?

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

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

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

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

Команды для оценки порядка модели

В этом фрагменте руководства вы используете struc, arxstruc, и selstruc чтобы оценить и сравнить простые полиномиальные (ARX) модели для области значений порядков моделей и задержек, и выбрать лучшие заказы на основе качества модели.

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

Порядок модели для первой комбинации вход-выход

В этом руководстве существует два входа в систему и один выход, и вы оцениваете порядки модели для каждой комбинации вход/выход независимо. Можно либо оценить задержки от двух входов одновременно, либо по одному входу за раз.

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

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

Чтобы оценить порядок модели для первой комбинации вход/выход:

Порядок модели для второй комбинации вход-выход

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

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

Чтобы оценить порядок модели для второй комбинации вход/выход:

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

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

Можно оценить передаточную функцию из этих порядков, используя tfest команда. Для оценки можно также выбрать просмотр отчета о прогрессе путем установки Display опция для on в наборе опций, созданном tfestOptions команда.

Opt = tfestOptions('Display','on');

Соберите порядки моделей и задержки в переменные, которые будут переданы tfest.

np = [2 1];
ioDelay = [2.5 5];

Оцените передаточную функцию.

mtf = tfest(Ze1,np,[],ioDelay,Opt);

Просмотрите коэффициенты модели.

mtf

mtf =
 
  From input "ConsumptionRate" to output "Production":
                  -1.382 s + 0.0008305
  exp(-2.5*s) * -------------------------
                s^2 + 1.014 s + 5.447e-12
 
  From input "Current" to output "Production":
                 5.93
  exp(-5*s) * ----------
              s + 0.5858
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: [2 1]   Number of zeros: [1 0]
   Number of free coefficients: 6
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                         
Estimated using TFEST on time domain data "Ze1".
Fit to estimation data: 78.92%                  
FPE: 14.22, MSE: 13.97                          

Отображение модели показывает более 85% подгонки к данным оценки.

Валидация модели

В этом фрагменте руководства вы создаете график, который сравнивает фактический выход и выход модели с помощью compare команда.

compare(Zv1,mtf)

Сравнение показывает, что подгонка составляет около 60%.

Остаточный анализ

Используйте resid команда для оценки характеристик невязок.

resid(Zv1,mtf)

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

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

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

Вы, должно быть, уже подготовили свои данные, как описано в разделе Подготовка данных.

Для определения порядков модели можно использовать следующие результаты оценочных заказов модели:

Используйте idproc команда для создания двух структур модели, по одной для каждой комбинации вход/выход:

midproc0 = idproc({'P2ZUD','P1D'}, 'TimeUnit', 'minutes');

Массив ячеек {'P2ZUD','P1D'} задает структуру модели для каждой комбинации вход/выход:

В примере также используется TimeUnit параметр для определения модуля времени, используемой.

Просмотр структуры модели и значений параметров

Просмотрите две получившиеся модели.

midproc0

midproc0 =
Process model with 2 inputs: y = G11(s)u1 + G12(s)u2
  From input 1 to output 1:                         
                       1+Tz*s                       
  G11(s) = Kp * ---------------------- * exp(-Td*s) 
                1+2*Zeta*Tw*s+(Tw*s)^2              
                                                    
         Kp = NaN                                   
         Tw = NaN                                   
       Zeta = NaN                                   
         Td = NaN                                   
         Tz = NaN                                   
                                                    
  From input 2 to output 1:                         
             Kp                                     
  G12(s) = ---------- * exp(-Td*s)                  
            1+Tp1*s                                 
                                                    
        Kp = NaN                                    
       Tp1 = NaN                                    
        Td = NaN                                    
                                                    
Parameterization:
    {'P2DUZ'}    {'P1D'}
   Number of free coefficients: 8
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

Значения параметров установлены в NaN потому что они еще не оценены.

Определение начальных догадок для временных задержек

Установите свойство задержки по времени для объекта модели в 2.5 мин и 5 min для каждой входной/выходной пары в качестве начальных догадок. Кроме того, установите верхний предел задержки, потому что доступны хорошие начальные предположения.

midproc0.Structure(1,1).Td.Value = 2.5;
midproc0.Structure(1,2).Td.Value = 5;
midproc0.Structure(1,1).Td.Maximum = 3;
midproc0.Structure(1,2).Td.Maximum = 7;

Оценка параметров модели с помощью procest

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

В зависимости от его аргументов, procest оценивает различные модели черного полинома. Можно использовать procestдля примера оценить параметры для структур линейного непрерывного времени передаточной функции, пространства состояний или полиномиальной модели. Как использовать procestнеобходимо предоставить структуру модели с неизвестными параметрами и данные оценки в качестве входных параметров.

Для этого фрагмента руководства необходимо уже определить структуру модели, как описано в разделе «Определение структуры модели процесса». Использование midproc0 как структура модели и Ze1 как данные оценки:

midproc = procest(Ze1,midproc0);
present(midproc)

                                                                  
midproc =                                                         
Process model with 2 inputs: y = G11(s)u1 + G12(s)u2              
  From input "ConsumptionRate" to output "Production":            
                       1+Tz*s                                     
  G11(s) = Kp * ---------------------- * exp(-Td*s)               
                1+2*Zeta*Tw*s+(Tw*s)^2                            
                                                                  
         Kp = -1.1807 +/- 0.29987                                 
         Tw = 1.6437 +/- 714.6                                    
       Zeta = 16.036 +/- 6958.9                                   
         Td = 2.426 +/- 64.276                                    
         Tz = -109.19 +/- 63.738                                  
                                                                  
  From input "Current" to output "Production":                    
             Kp                                                   
  G12(s) = ---------- * exp(-Td*s)                                
            1+Tp1*s                                               
                                                                  
        Kp = 10.264 +/- 0.048404                                  
       Tp1 = 2.049 +/- 0.054901                                   
        Td = 4.9175 +/- 0.034374                                  
                                                                  
Parameterization:                                                 
    {'P2DUZ'}    {'P1D'}                                          
   Number of free coefficients: 8                                 
   Use "getpvec", "getcov" for parameters and their uncertainties.
                                                                  
Status:                                                           
Termination condition: Maximum number of iterations reached..     
Number of iterations: 20, Number of function evaluations: 279     
                                                                  
Estimated using PROCEST on time domain data "Ze1".                
Fit to estimation data: 86.2%                                     
FPE: 6.081, MSE: 5.984                                            
More information in model's "Report" property.                    

В отличие от моделей полинома в дискретном времени, модели в непрерывном времени позволяют вам оценить задержки. В этом случае предполагаемые значения задержки отличаются от начальных догадок, которые вы задали 2.5 и 5 , соответственно. Большие неопределенности в оценочных значениях параметров G_1(s) указать, что динамика от входных 1 к выходу не захватываются хорошо.

Чтобы узнать больше об оценке моделей, см. «Модели процесса».

Валидация модели

В этом фрагменте руководства вы создаете график, который сравнивает фактический выход и выход модели.

compare(Zv1,midproc)

График показывает, что результат модели разумно согласен с измеренным выходом: существует согласие 60% между моделью и данными валидации.

Использование resid для выполнения остаточного анализа.

resid(Zv1,midproc)

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

Измените алгоритм итерационной оценки параметра на Левенберг-Марквардт.

Opt = procestOptions;
Opt.SearchMethod = 'lm';
midproc1 = procest(Ze1,midproc,Opt);

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

Оценка модели процесса с помощью модели шума

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

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

Opt = procestOptions;
Opt.DisturbanceModel = 'ARMA1';
midproc2 = procest(Ze1,midproc0,Opt);

Можно вводить 'dist' вместо 'DisturbanceModel'. Имена свойства не чувствительны к регистру, и вы должны включать только фрагмент имени, которая уникально идентифицирует свойство.

Сравните полученную модель с измеренными данными.

compare(Zv1,midproc2)

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

Выполните остаточный анализ.

resid(Zv1,midproc2)

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

Модели порядков для оценки полиномиальных моделей

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

Вы, должно быть, уже подготовили свои данные, как описано в разделе Подготовка данных.

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

Оценка линейной модели ARX

О моделях ARX.  Для системы с одним входом/одним выходом (SISO), структура модели ARX:

y (t) представляет вывод в момент t, u (t) представляет вход в момент t, na - количество полюсов, nb - количество нулей плюс 1, nk - количество выборок перед тем, как вход влияет на выход системы (называемый задержкой или потерей времени модели), и e (t) - нарушение порядка белого шума.

Структура модели ARX не различает полюсы для отдельных входных/выходных путей: деление уравнения ARX на A, которое содержит полюсы, показывает, что A появляется в знаменателе для обоих входов. Поэтому можно задать _na в сумме полюсов для каждой пары вход/выход, которая равна 3 в этом случае.

Продукт System Identification Toolbox оценивает параметры $$a_1\dots a_n$$ и $$b_1\dots b_n$$ используя заданные данные и порядок модели.

Оценка моделей ARX с использованием arx.  Использовать arx чтобы вычислить полиномиальные коэффициенты с помощью быстрого, нежесткого метода arx:

marx = arx(Ze1,'na',3,'nb',[2 1],'nk',[5 10]);
present(marx) % Displays model parameters
              % with uncertainty information

                                                                              
marx =                                                                        
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)                           
                                                                              
  A(z) = 1 - 1.027 (+/- 0.02907) z^-1 + 0.1678 (+/- 0.042) z^-2 + 0.01289 (   
                                                         +/- 0.02583) z^-3    
                                                                              
  B1(z) = 1.86 (+/- 0.189) z^-5 - 1.608 (+/- 0.1888) z^-6                     
                                                                              
  B2(z) = 1.612 (+/- 0.07392) z^-10                                           
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   na=3   nb=[2 1]   nk=[5 10]                           
   Number of free coefficients: 6                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Estimated using ARX on time domain data "Ze1".                                
Fit to estimation data: 90.7% (prediction focus)                              
FPE: 2.768, MSE: 2.719                                                        
More information in model's "Report" property.                                

MATLAB оценивает полиномы A , B1 , и B2. Неопределенность для каждого из параметров модели вычисляется как 1 стандартное отклонение и появляется в круглых скобках рядом с каждым значением параметров.

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

marx = arx(Ze1,[3 2 1 5 10]);

Доступ к данным моделям.  Модель, которую вы оценили, marx, является дискретным временем idpoly объект. Чтобы получить свойства этого объекта модели, можно использовать get функция:

get(marx)

                 A: [1 -1.0267 0.1678 0.0129]
                 B: {[0 0 0 0 0 1.8599 -1.6084]  [0 0 0 0 0 0 0 0 0 0 1.6118]}
                 C: 1
                 D: 1
                 F: {[1]  [1]}
    IntegrateNoise: 0
          Variable: 'z^-1'
           IODelay: [0 0]
         Structure: [1x1 pmodel.polynomial]
     NoiseVariance: 2.7436
        InputDelay: [2x1 double]
       OutputDelay: 0
                Ts: 0.5000
          TimeUnit: 'minutes'
         InputName: {2x1 cell}
         InputUnit: {2x1 cell}
        InputGroup: [1x1 struct]
        OutputName: {'Production'}
        OutputUnit: {'mg/min'}
       OutputGroup: [1x1 struct]
             Notes: [0x1 string]
          UserData: []
              Name: ''
      SamplingGrid: [1x1 struct]
            Report: [1x1 idresults.arx]

Вы можете получить доступ к информации, хранимой этими свойствами, используя запись через точку. Для примера можно вычислить дискретные полюсы модели, найдя корни A полином.

marx_poles = roots(marx.a)

marx_poles =

    0.7953
    0.2877
   -0.0564

В этом случае вы получаете доступ к A полином с использованием marx.a.

Модель marx описывает динамику системы с помощью трех дискретных полюсов.

Подробнее.  Чтобы узнать больше об оценке полиномиальных моделей, см. «Входно-выходные полиномиальные модели».

Для получения дополнительной информации о доступе к данным данных моделей Раздел «Извлечение данных».

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

О моделях пространства состояний.  Общая структура модели пространства состояний:

y (t) представляет выход в момент t, u (t) представляет вход в момент t, x (t) является вектором состояния в момент t, а e (t) является нарушением порядка белого шума.

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

Продукт System Identification Toolbox оценивает матрицы А, B, C, D и K пространства состояний, используя порядок модели и заданные вами данные.

Структура модели пространства состояний является хорошим выбором для быстрой оценки, потому что она содержит только два параметра: n количество полюсов (размер матрицы A) и nk - задержка.

Оценка моделей пространства состояний с использованием n4sid.  Используйте n4sid команда для определения области значений порядков моделей и оценки эффективности нескольких моделей пространства состояний (порядки 2-8):

mn4sid = n4sid(Ze1,2:8,'InputDelay',[4 9]);

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

Вертикальная ось является относительной мерой того, насколько каждое состояние способствует поведению ввода/вывода модели (журнал сингулярных значений ковариационной матрицы). Горизонтальная ось соответствует порядку модели n. Этот график рекомендует n=3, обозначенный красным прямоугольником.

В поле порядка Chosen Order отображается рекомендуемый порядок модели, 3 в этом случае по умолчанию. Выбор порядка можно изменить с помощью раскрывающегося списка Chosen Order. Примените значение в поле Chosen Order и закройте окно выбора заказа нажав Apply.

По умолчанию n4sid использует свободную параметризацию формы пространство состояний. Чтобы оценить каноническую форму, установите значение SSParameterization свойство к 'Canonical' . Можно также задать задержку ввода-вывода (в выборках), используя 'InputDelay' свойство.

mCanonical = n4sid(Ze1,3,'SSParameterization','canonical','InputDelay',[4 9]);
present(mCanonical);  %  Display model properties

                                                                              
mCanonical =                                                                  
  Discrete-time identified state-space model:                                 
    x(t+Ts) = A x(t) + B u(t) + K e(t)                                        
       y(t) = C x(t) + D u(t) + e(t)                                          
                                                                              
  A =                                                                         
                       x1                  x2                  x3             
   x1                   0                   1                   0             
   x2                   0                   0                   1             
   x3  0.0737 +/- 0.05919  -0.6093 +/- 0.1626    1.446 +/- 0.1287             
                                                                              
  B =                                                                         
             ConsumptionR             Current                                 
   x1   1.844 +/-   0.175   0.5633 +/-  0.122                                 
   x2   1.063 +/-  0.1673    2.308 +/- 0.1222                                 
   x3  0.2779 +/- 0.09615    1.878 +/- 0.1058                                 
                                                                              
  C =                                                                         
               x1  x2  x3                                                     
   Production   1   0   0                                                     
                                                                              
  D =                                                                         
               ConsumptionR       Current                                     
   Production             0             0                                     
                                                                              
  K =                                                                         
               Production                                                     
   x1  0.8674 +/- 0.03169                                                     
   x2  0.6849 +/- 0.04145                                                     
   x3  0.5105 +/- 0.04352                                                     
                                                                              
  Input delays (sampling periods): 4  9                                       
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   CANONICAL form with indices: 3.                                            
   Feedthrough: none                                                          
   Disturbance component: estimate                                            
   Number of free coefficients: 12                                            
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Estimated using N4SID on time domain data "Ze1".                              
Fit to estimation data: 91.39% (prediction focus)                             
FPE: 2.402, MSE: 2.331                                                        
More information in model's "Report" property.                                

mCT1 = n4sid(Ze1, 3, 'Ts', 0, 'InputDelay', [2.5 5])
mCT2 = ssest(Ze1, 3,'InputDelay', [2.5 5])

Подробнее.  Чтобы узнать больше об оценке моделей пространства состояний, смотрите Модели пространства состояний.

Оценка модели Бокса-Дженкинса

О моделях Бокса-Дженкинса.  Общая структура Box-Jenkins (BJ):

Чтобы оценить модель BJ, необходимо задать параметры _nb, _nf, _nc, _nd и _nk.

В то время как структура модели ARX не различает полюсы для отдельных входных/выходных путей, модель BJ обеспечивает большую гибкость в моделировании полюсов и нулей нарушения порядка отдельно от полюсов и нулей динамики системы.

Оценка модели BJ с использованием полиеста.  Вы можете использовать polyest для оценки модели BJ. polyest является итерационным методом и имеет следующий общий синтаксис:

polyest(data,[na nb nc nd nf nk]);

Чтобы оценить модель BJ, введите:

na = 0;
nb = [ 2 1 ];
nc = 1;
nd = 1;
nf = [ 2 1 ];
nk = [ 5 10];
mbj = polyest(Ze1,[na nb nc nd nf nk]);

Эта команда задает nf=2 , nb=2 , nk=5 для первого входа и nf=nb=1 и nk=10 для второго входа.

Отобразите информацию о модели.

present(mbj)

                                                                              
mbj =                                                                         
Discrete-time BJ model: y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t)              
  B1(z) = 1.823 (+/- 0.1792) z^-5 - 1.315 (+/- 0.2367) z^-6                   
                                                                              
  B2(z) = 1.791 (+/- 0.06431) z^-10                                           
                                                                              
  C(z) = 1 + 0.1068 (+/- 0.04009) z^-1                                        
                                                                              
  D(z) = 1 - 0.7452 (+/- 0.02694) z^-1                                        
                                                                              
  F1(z) = 1 - 1.321 (+/- 0.06936) z^-1 + 0.5911 (+/- 0.05514) z^-2            
                                                                              
  F2(z) = 1 - 0.8314 (+/- 0.006441) z^-1                                      
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   nb=[2 1]   nc=1   nd=1   nf=[2 1]                     
   nk=[5 10]                                                                  
   Number of free coefficients: 8                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Termination condition: Near (local) minimum, (norm(g) < tol)..                
Number of iterations: 6, Number of function evaluations: 13                   
                                                                              
Estimated using POLYEST on time domain data "Ze1".                            
Fit to estimation data: 90.75% (prediction focus)                             
FPE: 2.733, MSE: 2.689                                                        
More information in model's "Report" property.                                

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

Полиномы C и D дать числитель и знаменатель шумовой модели, соответственно.

mbj = bj(Ze1,[2 1 1 1 2 1 5 10]);

Подробнее.  Чтобы узнать больше об идентификации моделей полинома ввода-вывода, таких как BJ, смотрите Модели Полинома ввода-вывода.

Сравнение выхода модели с измеренным выходом

Сравните выходные данные моделей ARX, пространства состояний и Box-Jenkins с измеренным выходом.

compare(Zv1,marx,mn4sid,mbj)

compare строит графики измеренного выхода в наборе данных валидации относительно моделируемого выхода моделей. Входные данные из набора данных валидации служат входными данными для моделей.

Выполните остаточный анализ на моделях ARX, пространства состояний и Box-Jenkins.

resid(Zv1,marx,mn4sid,mbj)

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

Симуляция выхода модели

В этом фрагменте руководства вы симулируете выход модели. Вы, должно быть, уже создали модель непрерывного времени midproc2, как описано в Оценке Моделей Процесса.

Симуляция выхода модели требует следующей информации:

Например, следующие команды используют iddata и idinput команды для создания набора входных данных и использования sim для симуляции выхода модели:

% Create input for simulation
U = iddata([],idinput([200 2]),'Ts',0.5,'TimeUnit','min');
% Simulate the response setting initial conditions equal to zero
ysim_1 = sim(midproc2,U);

Чтобы максимизировать подгонку между моделируемой характеристикой модели и измеренным выходом для того же входа, можно вычислить начальные условия, соответствующие измеренным данным. Лучшие подходящие начальные условия могут быть получены при помощи findstates на версии предполагаемой модели в пространстве состояний. Следующие команды оценивают начальные состояния X0est из набора данных Zv1:

% State-space version of the model needed for computing initial states
midproc2_ss = idss(midproc2);
X0est = findstates(midproc2_ss,Zv1);

Затем симулируйте модель с помощью начальных состояний, оцененных из данных.

% Simulation input
Usim = Zv1(:,[],:);
Opt = simOptions('InitialCondition',X0est);
ysim_2 = sim(midproc2_ss,Usim,Opt);

Сравните моделируемый и измеренный выходы на графике.

figure
plot([ysim_2.y, Zv1.y])
legend({'model output','measured'})
xlabel('time'), ylabel('Output')

Предсказание выхода

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

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

predict(midproc2,Ze1,5)

Сравните предсказанные выходные значения с измеренными выходными значениями. Третий аргумент compare задает пятиэтапное предсказание. Когда вы не задаете третий аргумент, compare принимает бесконечный горизонт предсказания и симулирует выход модели вместо этого.

compare(Ze1,midproc2,5)

Использовать pe чтобы вычислить ошибку предсказания Err между предсказанным выходом midproc2 и измеренный выход. Затем постройте график спектра ошибок с помощью spectrum команда.

[Err] = pe(midproc2,Zv1);
spectrum(spa(Err,[],logspace(-2,2,200)))

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