Создание структурированных и пользовательских моделей Используя System Identification Toolbox™

Открыть скрипт

В этом примере показано, как оценить параметры в пользовательских структурах модели. Такие структуры заданы IDGREY (линейное пространство состояний) или IDNLGREY (нелинейное пространство состояний) модели. Мы займемся расследованиями, как присвоить структуру, зафиксировать параметры и создать зависимости среди них.

Данные об эксперименте

Мы исследуем данные, произведенные (симулированным) двигателем постоянного тока. Мы сначала загружаем данные:

load dcmdata
who

Your variables are:

text  u     y

Матричный y содержит эти два выходных параметров: y1 угловое положение вала двигателя и y2 скорость вращения. Существует 400 выборок данных, и шаг расчета составляет 0,1 секунды. Вход содержится в векторном u. Это - входное напряжение к двигателю.

z = iddata(y,u,0.1); % The IDDATA object
z.InputName = 'Voltage';
z.OutputName = {'Angle';'AngVel'};
plot(z(:,1,:))

Фигура: данные об измерении: напряжение, чтобы повернуть

plot(z(:,2,:))

Фигура: данные об измерении: напряжение к скорости вращения

Выбор структуры модели

     d/dt x = A x + B u + K e
        y   = C x + D u + e

Мы создадим модель двигателя постоянного тока. Динамика двигателя известна. Если мы выбираем x1 в качестве углового положения и x2 как скорость вращения, легко настроить модель в пространстве состояний следующих воздействий пренебрежения символа: (см. Пример 4.1 в Ljung (1999):

         | 0     1  |      | 0   |
d/dt x = |          | x  + |     | u
         | 0   -th1 |      | th2 |

        | 1  0 |
   y =  |      | x
        | 0  1 |

Параметр th1 здесь инверсия, постоянная во времени из двигателя и th2 таково что th2/th1 статическое усиление от входа до скорости вращения. (См. Ljung (1987) для как th1 и th2 относитесь к физическим параметрам двигателя). Мы оценим эти два параметра из наблюдаемых данных. Структура модели (параметрированное пространство состояний) описанный выше может быть представлена в MATLAB® с помощью моделей IDSS и IDGREY. Эти модели позволяют вам выполнить оценку параметров с помощью экспериментальных данных.

Спецификация номинальной (первоначальной) модели

Если мы предполагаем, что th1=1 и th2 = 0.28 получаем номинальную или первоначальную модель

A = [0 1; 0 -1]; % initial guess for A(2,2) is -1
B = [0; 0.28]; % initial guess for B(2) is 0.28
C = eye(2);
D = zeros(2,1);

и мы группируем это в объект модели IDSS:

ms = idss(A,B,C,D);

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

ms.Structure.a

 
ans =
 
       Name: 'A'
      Value: [2x2 double]
    Minimum: [2x2 double]
    Maximum: [2x2 double]
       Free: [2x2 logical]
      Scale: [2x2 double]
       Info: [2x2 struct]

 
1x1 param.Continuous

ms.Structure.a.Value
ms.Structure.a.Free

ans =

     0     1
     0    -1


ans =

  2x2 logical array

   1   1
   1   1

Спецификация свободных (независимых) параметров Используя модели IDSS

Таким образом, мы должны теперь отметить это, это только (2,2) и B (2,1), которые являются свободными параметрами, которые будут оценены.

ms.Structure.a.Free = [0 0; 0 1];
ms.Structure.b.Free = [0; 1];
ms.Structure.c.Free = 0; % scalar expansion used
ms.Structure.d.Free = 0;
ms.Ts = 0;  % This defines the model to be continuous

Первоначальная модель

ms % Initial model

ms =
  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
   x1   0   1
   x2   0  -1
 
  B = 
         u1
   x1     0
   x2  0.28
 
  C = 
       x1  x2
   y1   1   0
   y2   0   1
 
  D = 
       u1
   y1   0
   y2   0
 
  K = 
       y1  y2
   x1   0   0
   x2   0   0
 
Parameterization:
   STRUCTURED form (some fixed coefficients in  A, B, C).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

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

Оценка свободных параметров модели IDSS

Ошибка предсказания (наибольшее правдоподобие) оценка параметров теперь вычисляется:

dcmodel = ssest(z,ms,ssestOptions('Display','on'));
dcmodel

dcmodel =
  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
   x1       0       1
   x2       0  -4.013
 
  B = 
       Voltage
   x1        0
   x2    1.002
 
  C = 
           x1  x2
   Angle    1   0
   AngVel   0   1
 
  D = 
           Voltage
   Angle         0
   AngVel        0
 
  K = 
        Angle  AngVel
   x1       0       0
   x2       0       0
 
Parameterization:
   STRUCTURED form (some fixed coefficients in  A, B, C).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 2
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                       
Estimated using SSEST on time domain data "z".
Fit to estimation data: [98.35;84.42]%        
FPE: 0.001071, MSE: 0.1192

Ориентировочные стоимости параметров вполне близко к используемым, когда данные были симулированы (-4 и 1). Чтобы оценить качество модели, мы можем симулировать модель с фактическим входом и сравнить его с фактическим выходом.

compare(z,dcmodel);

Мы можем теперь, например, построить нули и полюса и их области неопределенности. Мы будем чертить области, соответствующие 3 стандартным отклонениям, поскольку модель довольно точна. Обратите внимание на то, что полюс в начале координат является абсолютно бесспорным, поскольку это - часть структуры модели; интегратор от скорости вращения до положения.

clf
showConfidence(iopzplot(dcmodel),3)

Теперь мы можем сделать различные модификации. С 1,2 элементами из матрицы A (зафиксированный к 1) говорит нам тот x2 производная x1. Предположим, что датчики не калибруются, так, чтобы мог быть неизвестный коэффициент пропорциональности. Чтобы включать оценку такой константы, мы только "выпускаем" A(1,2) и переоценка:

dcmodel2 = dcmodel;
dcmodel2.Structure.a.Free(1,2) = 1;
dcmodel2 = pem(z,dcmodel2,ssestOptions('Display','on'));

State-space Model Identification
                                
Estimation data: Time domain data z          
Data has 2 outputs, 1 inputs and 400 samples.
Number of states: 2                          

Algorithm: Nonlinear least squares with automatically chosen line search method
 <br>
------------------------------------------------------------------------------------------
 <br>
                          Norm of      First-order    Improvement (%) <br> Iteration       Cost       step         optimality     Expected   Achieved    Bisections <br>------------------------------------------------------------------------------------------
     0     0.00106082          -             40      0.0256           -         -                                                                                                                                                                                                                                                                                                                                                                         
     1     0.00106055     0.0038           8.21      0.0256      0.0255         0
     2     0.00106055   5.48e-05       0.000879    1.13e-05    1.13e-05         0
------------------------------------------------------------------------------------------
Estimating parameter covariance...
done.
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 2, Number of function evaluations: 5    
                                                              
Status: Estimated using PEM                                   
Fit to estimation data: [98.35;84.42]%, FPE: 0.00107658

Получившаяся модель

dcmodel2

dcmodel2 =
  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
   x1       0  0.9975
   x2       0  -4.011
 
  B = 
       Voltage
   x1        0
   x2    1.004
 
  C = 
           x1  x2
   Angle    1   0
   AngVel   0   1
 
  D = 
           Voltage
   Angle         0
   AngVel        0
 
  K = 
        Angle  AngVel
   x1       0       0
   x2       0       0
 
Parameterization:
   STRUCTURED form (some fixed coefficients in  A, B, C).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 3
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                     
Estimated using PEM on time domain data "z".
Fit to estimation data: [98.35;84.42]%      
FPE: 0.001077, MSE: 0.1192

Мы находим что предполагаемый A(1,2) близко к 1. Чтобы сравнить две модели, мы используем compare команда:

compare(z,dcmodel,dcmodel2)

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

Предположим, что мы точно знаем статическое усиление двигателя постоянного тока (от входного напряжения до скорости вращения, например, из предыдущего эксперимента переходного процесса. Если статическим усилением является G, и постоянная времени двигателя является t, затем модель в пространстве состояний становится

          |0     1|    |  0  |
d/dt x =  |       |x + |     | u
          |0  -1/t|    | G/t |

          |1   0|
   y   =  |     | x
          |0   1|

С G известный, существует зависимость между записями в различных матрицах. Для того, чтобы описать, что, ранее используемый путь "Свободными" параметрами не будет достаточен. Мы таким образом должны записать файл MATLAB, который производит A, B, C, и D, и опционально также K и матрицы X0 как выходные параметры, для каждого данного вектора параметра как вход. Это также берет вспомогательные аргументы в качестве входных параметров, так, чтобы пользователь мог изменить определенные вещи в структуре модели, не имея необходимость редактировать файл. В этом случае мы позволяем известному статическому усилению G будьте введены как таковые аргумент. Файл, который был записан, имеет имя 'motorDynamics.m'.

type motorDynamics

function [A,B,C,D,K,X0] = motorDynamics(par,ts,aux)
%MOTORDYNAMICS ODE file representing the dynamics of a motor.
%
%   [A,B,C,D,K,X0] = motorDynamics(Tau,Ts,G)
%   returns the State Space matrices of the DC-motor with
%   time-constant Tau (Tau = par) and known static gain G. The sample
%   time is Ts.
%
%   This file returns continuous-time representation if input argument Ts
%   is zero. If Ts>0, a discrete-time representation is returned. To make
%   the IDGREY model that uses this file aware of this flexibility, set the
%   value of Structure.FcnType property to 'cd'. This flexibility is useful
%   for conversion between continuous and discrete domains required for
%   estimation and simulation.
%
%   See also IDGREY, IDDEMO7.

%   L. Ljung
%   Copyright 1986-2015 The MathWorks, Inc.

t = par(1);
G = aux(1);

A = [0 1;0 -1/t];
B = [0;G/t];
C = eye(2);
D = [0;0];
K = zeros(2);
X0 = [0;0];
if ts>0 % Sample the model with sample time Ts
   s = expm([[A B]*ts; zeros(1,3)]);
   A = s(1:2,1:2);
   B = s(1:2,3);
end

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

par_guess = 1;

Мы также даем значение 0.25 вспомогательной переменной G (усиление) и шаг расчета.

aux = 0.25;
dcmm = idgrey('motorDynamics',par_guess,'cd',aux,0);

Постоянная времени теперь оценивается

dcmm = greyest(z,dcmm,greyestOptions('Display','on'));

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

dcmm

dcmm =
  Continuous-time linear grey box model defined by "motorDynamics" function:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
           x1      x2
   x1       0       1
   x2       0  -4.006
 
  B = 
       Voltage
   x1        0
   x2    1.001
 
  C = 
           x1  x2
   Angle    1   0
   AngVel   0   1
 
  D = 
           Voltage
   Angle         0
   AngVel        0
 
  K = 
        Angle  AngVel
   x1       0       0
   x2       0       0
 
  Model parameters:
   Par1 = 0.2496
 
Parameterization:
   ODE Function: motorDynamics
   (parametrizes both continuous- and discrete-time equations)
   Disturbance component: parameterized by the ODE function
   Initial state: parameterized by the ODE function
   Number of free coefficients: 1
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                         
Estimated using GREYEST on time domain data "z".
Fit to estimation data: [98.35;84.42]%          
FPE: 0.00107, MSE: 0.1193

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

compare(z,dcmm,dcmodel)

Они ясно очень близки.

Оценка многомерных моделей ARX

Часть пространства состояний тулбокса также обрабатывает многомерный (несколько выходных параметров) модели ARX. Многомерной моделью ARX мы означаем следующее:

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

Здесь (q) ny | нью-йоркская матрица, записи которой являются полиномами в операторе задержки 1/q. k-l элемент обозначается:

$a_{kl}(q)$

где:

$a_{kl}(q) = 1 + a_1 q^{-1} + .... + a_{nakl} q^{-nakl} q$

Это - таким образом полином в 1/q из степени nakl.

Так же B (q) является ny | матрица ню, kj-элемент которой:

$b_{kj}(q) = b_0 q^{-nkk}+b_1 q^{-nkkj-1}+ ... +b_{nbkj} q^{-nkkj-nbkj}$

Существует таким образом задержка nkkj от входного номера j выводить номер k. Наиболее распространенный способ создать тех состоял бы в том, чтобы использовать команду ARX. Порядки заданы как: nn = [na nb nk] с na быть ny-by-ny матрица, чей kj- записью является nakj; nb и nk заданы так же.

Давайте протестируем некоторые модели ARX на dc-данных. Сначала мы могли просто создать общую модель второго порядка:

dcarx1 = arx(z,'na',[2,2;2,2],'nb',[2;2],'nk',[1;1])

dcarx1 =
Discrete-time ARX model:                                                   
  Model for output "Angle": A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + e_1(t)
    A(z) = 1 - 0.5545 z^-1 - 0.4454 z^-2                                   
                                                                           
    A_2(z) = -0.03548 z^-1 - 0.06405 z^-2                                  
                                                                           
    B(z) = 0.004243 z^-1 + 0.006589 z^-2                                   
                                                                           
  Model for output "AngVel": A(z)y_2(t) = - A_i(z)y_i(t) + B(z)u(t) + e_2(t)
    A(z) = 1 - 0.2005 z^-1 - 0.2924 z^-2                                    
                                                                            
    A_1(z) = 0.01849 z^-1 - 0.01937 z^-2                                    
                                                                            
    B(z) = 0.08642 z^-1 + 0.03877 z^-2                                      
                                                                            
Sample time: 0.1 seconds
  
Parameterization:
   Polynomial orders:   na=[2 2;2 2]   nb=[2;2]   nk=[1;1]
   Number of free coefficients: 12
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using ARX on time domain data "z".             
Fit to estimation data: [97.87;83.44]% (prediction focus)
FPE: 0.002157, MSE: 0.1398

Результат, dcarx1, хранится как модель IDPOLY, и все предыдущие команды применяются. Мы могли, например, явным образом перечислить полиномы ARX:

dcarx1.a

ans =

  2x2 cell array

    {[1 -0.5545 -0.4454]}    {[0 -0.0355 -0.0640]}
    {[ 0 0.0185 -0.0194]}    {[1 -0.2005 -0.2924]}

как массивы ячеек, где e.g. {1,2} элементом dcarx1.a является полином (1,2) описанный ранее, имея отношение y2 к y1.

Мы могли также протестировать структуру, где мы знаем тот y1 получен путем фильтрации y2 через фильтр первого порядка. (Угол является интегралом скорости вращения). Мы могли затем также постулировать динамику первого порядка от входа, чтобы вывести номер 2:

na = [1 1; 0 1];
nb = [0 ; 1];
nk = [1 ; 1];
dcarx2 = arx(z,[na nb nk])

dcarx2 =
Discrete-time ARX model:                                                   
  Model for output "Angle": A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + e_1(t)
    A(z) = 1 - 0.9992 z^-1                                                 
                                                                           
    A_2(z) = -0.09595 z^-1                                                 
                                                                           
    B(z) = 0                                                               
                                                                           
  Model for output "AngVel": A(z)y_2(t) = B(z)u(t) + e_2(t)
    A(z) = 1 - 0.6254 z^-1                                 
                                                           
    B(z) = 0.08973 z^-1                                    
                                                           
Sample time: 0.1 seconds
  
Parameterization:
   Polynomial orders:   na=[1 1;0 1]   nb=[0;1]   nk=[1;1]
   Number of free coefficients: 4
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using ARX on time domain data "z".             
Fit to estimation data: [97.52;81.46]% (prediction focus)
FPE: 0.003452, MSE: 0.177

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

compare(z,dcmodel,dcmm,dcarx2)

Наконец, мы могли сравнить bodeplots, полученный из входа, чтобы вывести один для различных моделей при помощи bode: Сначала выход:

dcmm2 = idss(dcmm); % convert to IDSS for subreferencing
bode(dcmodel(1,1),'r',dcmm2(1,1),'b',dcarx2(1,1),'g')

Second output:

bode(dcmodel(2,1),'r',dcmm2(2,1),'b',dcarx2(2,1),'g')

Две первых модели находятся более или менее в точном соглашении. Модели ARX не так хороши, из-за смещения, вызванного цветным ошибочным шумом уравнения. (У нас был белый шум измерения в симуляциях).

Заключения

Оценка моделей с предварительно выобранными структурами может быть выполнена с помощью тулбокса System Identification. В форме пространства состояний параметры могут быть зафиксированы к их известным значениям или ограничены лечь в предписанной области значений. Если отношение между параметрами или другими ограничениями должно быть задано, объекты IDGREY могут использоваться. Модели IDGREY оценивают заданный пользователями файл MATLAB для оценки параметров системы в пространстве состояний. Многомерные модели ARX предлагают другую опцию для того, чтобы быстро оценить мультивыходные модели с заданной пользователями структурой.

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