Создание структурированных и пользовательских моделей с использованием 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

Поэтому теперь мы должны отметить, что только A (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 известно, что существует зависимость между записями в различных матрицах. В порядок для описания этого, ранее использованный способ с параметрами «Free» будет недостаточным. Таким образом, мы должны записать файл 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

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

compare(z,dcmm,dcmodel)

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

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

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

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

Здесь A (q) является ny | 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 | nu матрицей, чей 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-entry 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]}

как массивы ячеек, где, например, элемент {1,2} dcarx1.a является полиномом 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)

Наконец, мы могли бы сравнить бодеплоты, полученные из входа, чтобы вывести одну для различных моделей при помощи 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-модели не так хороши, из-за смещения, вызванного небелым шумом ошибки уравнения. (У нас был белый шум измерения в симуляциях).

Заключения

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

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