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

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

compare(z,dcmm,dcmodel)

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

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

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

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

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

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

Заключения

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

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