Matlabでシステム同定

www.mathworks.com

data objectをインポート

iddata形式をインポート

iddata形式の表現
www.mathworks.com

例：の場合

s=tf('s');
G=1/(s^2+s+1);
ts=0.1;
t=0:ts:20;
u=ones(length(t),1)';
[y,t]=lsim(G,u,t);
z=iddata(y,u',0.1); %uを転置してカラムベクトルに直す
plot(z);

Options = tfestOptions;      
Options.Display = 'on';      
Options.WeightingFilter = [];
tf1 = tfest(z, 2, 1, Options)                              

実行結果（極2, 零1を仮定）

tf1 =
 
  From input "u1" to output "y1":
  -2.694e-10 s + 1
  ----------------
    s^2 + s + 1
 
Continuous-time identified transfer function.

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

Status:                                       
Estimated using TFEST on time domain data "z".
Fit to estimation data: 100%                  
FPE: 3.567e-32, MSE: 3.36e-32   

実際には零点はないので、零を0でやってみる。

Options = tfestOptions;      
Options.Display = 'on';      
Options.WeightingFilter = [];
tf1 = tfest(z, 2, 0, Options)    

実行結果

tf1 =
 
  From input "u1" to output "y1":
       1
  -----------
  s^2 + s + 1
 
Continuous-time identified transfer function.

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

Status:                                       
Estimated using TFEST on time domain data "z".
Fit to estimation data: 100%                  
FPE: 4.913e-32, MSE: 4.674e-32  

完全に復元（同定）できた。
ただし、事前知識（極が２個で零点が０個）を仮定したのがいただけない。ハイパーパラメータも含めて自動で推定したい。
Estimation reportを使えばできそう。フィット率やMSEが取得できる。AICを使えばモデルの良さを検討できそう。
www.mathworks.com
ちょっとやってみる。

２つの同定結果を格納

tf1_zero1 = tf1_zero1 = tfest(z, 2, 1, Options);
tf1_zero0 =  tf1_zero0 = tfest(z, 2, 0, Options);

フィット結果を確認　

>> tf1_zero1.report.Fit;
ans = 
  struct with fields:

    FitPercent: 100.0000
       LossFcn: 3.1981e-32
           MSE: 3.3602e-32
           FPE: 3.5670e-32
           AIC: -1.3984e+04
          AICc: -1.3984e+04
          nAIC: -72.4110
           BIC: -1.3964e+04

>> tf1_zero0.report.Fit
ans = 
  struct with fields:

    FitPercent: 100.0000
       LossFcn: 3.7043e-32
           MSE: 4.6743e-32
           FPE: 4.9127e-32
           AIC: -1.3920e+04
          AICc: -1.3920e+04
          nAIC: -72.0909
           BIC: -1.3903e+04

どっちもフィット率は100%。MSEは零点が１個（複雑なほうのモデル）のほうが高くなっている。（過学習的な）
AICについても、零点が１のほうが若干ではあるが低くなっているため、AICを使えばベストのモデルが求まるわけではなさあそうだが、アタリをつけることはできそう。（極6, 零6の複雑なモデルで試したらAICがかなり高くなったので）

以下は、極が最大m個、零点が最大n個の選択肢の中から、最も良いモデルを調べるプログラム
m=6, n=6くらいでやってみれば実用上はいけるかな？

clear; close all;

s=tf('s');
G=1/(s^2+s+1);
ts=0.1;
t=0:ts:20;
nx=10;
u=[zeros(nx,1)',ones(length(t)-nx,1)'];
[y,t]=lsim(G,u,t);
z=iddata(y,u',0.1); %uを転置してカラムベクトルに直す
plot(z);

Options = tfestOptions;      
Options.Display = 'off';      
Options.WeightingFilter = [];

bestAIC = 99999;
m = 6;
n = 6;
bestm = [];
bestn= [];
for idx=1:m
    for jdx=0:n        
        if jdx>idx 
            continue;
        end
        tf = tfest(z, idx, jdx, Options);
        currentAIC = tf.report.Fit.AIC;
        fprintf('pole:%d, zero:%d, current AIC:%d, bestAIC: %d\n',idx,jdx,currentAIC,bestAIC);
        if currentAIC  < bestAIC
            bestAIC = currentAIC;
            bestm = idx;
            bestn = jdx;
        end        
    end   
end


besttf = tfest(z, bestm, bestn, Options)
figure();title('System ID result');
[y,t]=lsim(G,u,t);
plot(t,y);
hold on;
[y2,t2]=lsim(besttf,u,t);
plot(t2,y2);
legend('original','systemID');

gist.github.com

www.youtube.com

飛行機の例でやってみる
ctms.engin.umich.edu

エレベータ角δからピッチ角θへの伝達関数は以下

上記プログラム（findbest.m）の実行結果

besttf =
 
  From input "u1" to output "y1":
  8.169e-19 s^2 - 1.059e-17 s + 1
  -------------------------------
            s^2 + s + 1
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 2
   Number of free coefficients: 5
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                       
Estimated using TFEST on time domain data "z".
Fit to estimation data: 100%                  
FPE: 6.975e-38, MSE: 6.505e-38                

分子のs^2の項とsの項は微小だから無視すると、もとの伝達関数に完全に一致している。すごすぎ。
ステップ応答がわかれば何でも伝達関数求まるんじゃないかという気がしてきた。

オーダの小さい係数(のオーダ以下)を無視するように変更した
https://www.mathworks.com/matlabcentral/answers/516746-remove-small-numbers-from-transfer-function

clear; close all;

s=tf('s');
%G=1/(s^2+s+1);
G = (1.151*s+0.1774)/(s^3+0.739*s^2+0.921*s); %http://ctms.engin.umich.edu/CTMS/index.php?example=AircraftPitch&section=SystemModeling
ts=0.1;
t=0:ts:20;
nx=10;
u=[zeros(nx,1)',ones(length(t)-nx,1)'];
[y,t]=lsim(G,u,t);
z=iddata(y,u',0.1); %uを転置してカラムベクトルに直す
plot(z);

Options = tfestOptions;      
Options.Display = 'off';      
Options.WeightingFilter = [];

bestAIC = 99999;
m = 6;
n = 6;
bestm = [];
bestn= [];
for idx=1:m
    for jdx=0:n        
        if jdx>idx 
            continue;
        end
        current_tf = tfest(z, idx, jdx, Options);
        currentAIC = current_tf.report.Fit.AIC;
        fprintf('pole:%d, zero:%d, current AIC:%d, bestAIC: %d\n',idx,jdx,currentAIC,bestAIC);
        if currentAIC  < bestAIC
            bestAIC = currentAIC;
            bestm = idx;
            bestn = jdx;
        end        
    end   
end

besttf = tfest(z, bestm, bestn, Options);
[n,d] = tfdata(besttf);
n = cellfun(@(x) {x.*(abs(x)>1e-7)}, n);
d = cellfun(@(x) {x.*(abs(x)>1e-7)}, d);
besttf = tf(n, d)

figure();title('System ID result');
[y,t]=lsim(G,u,t);
plot(t,y);
hold on;
[y2,t2]=lsim(besttf,u,t);
plot(t2,y2);
legend('original','systemID');

実行結果

pole:1, zero:0, current AIC:-2.369800e+01, bestAIC: 99999
pole:1, zero:1, current AIC:-4.843524e+01, bestAIC: -2.369800e+01
pole:2, zero:0, current AIC:9.121234e+02, bestAIC: -4.843524e+01
pole:2, zero:1, current AIC:6.022577e+02, bestAIC: -4.843524e+01
pole:2, zero:2, current AIC:6.261909e+02, bestAIC: -4.843524e+01
pole:3, zero:0, current AIC:-8.554004e+02, bestAIC: -4.843524e+01
pole:3, zero:1, current AIC:-1.365531e+04, bestAIC: -8.554004e+02
pole:3, zero:2, current AIC:-1.367792e+04, bestAIC: -1.365531e+04
pole:3, zero:3, current AIC:-1.377390e+04, bestAIC: -1.367792e+04
pole:4, zero:0, current AIC:7.872736e+02, bestAIC: -1.377390e+04
pole:4, zero:1, current AIC:-5.806630e+03, bestAIC: -1.377390e+04
pole:4, zero:2, current AIC:-1.333585e+04, bestAIC: -1.377390e+04
pole:4, zero:3, current AIC:4.716128e+02, bestAIC: -1.377390e+04
pole:4, zero:4, current AIC:-1.086533e+04, bestAIC: -1.377390e+04
pole:5, zero:0, current AIC:-2.063570e+02, bestAIC: -1.377390e+04
pole:5, zero:1, current AIC:6.484087e+02, bestAIC: -1.377390e+04
pole:5, zero:2, current AIC:7.562263e+02, bestAIC: -1.377390e+04
pole:5, zero:3, current AIC:-1.236406e+04, bestAIC: -1.377390e+04
pole:5, zero:4, current AIC:-1.062038e+02, bestAIC: -1.377390e+04
pole:5, zero:5, current AIC:-1.820478e+03, bestAIC: -1.377390e+04
pole:6, zero:0, current AIC:8.013079e+02, bestAIC: -1.377390e+04
pole:6, zero:1, current AIC:8.634470e+02, bestAIC: -1.377390e+04
pole:6, zero:2, current AIC:5.983620e+02, bestAIC: -1.377390e+04
pole:6, zero:3, current AIC:-4.926007e+02, bestAIC: -1.377390e+04
pole:6, zero:4, current AIC:-1.035856e+04, bestAIC: -1.377390e+04
pole:6, zero:5, current AIC:-2.579125e+03, bestAIC: -1.377390e+04
pole:6, zero:6, current AIC:-3.631200e+03, bestAIC: -1.377390e+04

besttf =
 
      1.151 s + 0.1774
  -------------------------
  s^3 + 0.739 s^2 + 0.921 s
 
Continuous-time transfer function.

バッチリだ。
実用的すぎる