System Identification Toolbox のデータおよびモデルオブジェクト

スクリプトを開く

この例は、System Identification Toolbox™ で使用できるデータおよびモデルオブジェクトを管理する方法を示しています。システム同定とは、データからモデルを構築することです。データセットは、複数の情報で特性が決まります。こうした情報には、入出力信号、サンプル時間、変数名、単位などがあります。同様に、推定モデルには、推定パラメーター、その共分散行列、モデル構造など、異なる種類の情報が含まれます。

これはつまり、データおよびモデル周辺の適切な情報をオブジェクトにパッケージ化するときに最適であり、必要であることを意味します。System Identification Toolbox には、このようなオブジェクトが多数含まれています。これらのオブジェクトの基本機能についてもこの例で説明します。

IDDATA オブジェクト

最初にいくつかのデータを作成します。

u = sign(randn(200,2)); % 2 inputs
y = randn(200,1);       % 1 output
ts = 0.1;               % The sample time

入力と出力を 1 つのオブジェクトで収集するには、以下を実行します。

z = iddata(y,u,ts);

名前を入力すると、データ情報が表示されます。

z =

Time domain data set with 200 samples.
Sample time: 0.1 seconds               
                                       
Outputs      Unit (if specified)       
   y1                                  
                                       
Inputs       Unit (if specified)       
   u1                                  
   u2

データは、plot(z) と同様に plot コマンドで iddata としてプロットされます。サブプロットを継続して進めるには、いずれかのキーを押します。ここでは、個別のチャネルをプロットします。

plot(z(:,1,1)) % Data subset with Input 1 and Output 1.

plot(z(:,1,2)) % Data subset with Input 2 and Output 1.

出力と入力を取得するには、以下を使用します。

u = z.u;   % or, equivalently u = get(z,'u');
y = z.y;   % or, equivalently y = get(z,'y');

データの一部を選択するには

zp = z(48:79);

最初の出力と 2 番目の入力を選択するには

zs = z(:,1,2);  % The ':' refers to all the data time points.

サブ選択は統合できます。

plot(z(45:54,1,2)) % samples 45 to 54 of response from second input to the first output.

チャネルには、既定の名前 'y1'、'u2' などが付けられます。この名前は、次のコマンドで任意の値に変更できます。

set(z,'InputName',{'Voltage';'Current'},'OutputName','Speed');

同様に、以下のように指定できます。

z.inputn = {'Voltage';'Current'}; % Autofill is used for properties
z.outputn = 'Speed';    % Upper and lower cases are also ignored

データ記録とプロットの場合、単位も設定できます。

z.InputUnit = {'Volt';'Ampere'};
z.OutputUnit = 'm/s';
z

z =

Time domain data set with 200 samples.
Sample time: 0.1 seconds               
                                       
Outputs       Unit (if specified)      
   Speed         m/s                   
                                       
Inputs        Unit (if specified)      
   Voltage       Volt                  
   Current       Ampere

現在のプロパティはすべて、(あらゆるオブジェクトで) 以下のように取得します。

get(z)

ans = 

  struct with fields:

              Domain: 'Time'
                Name: ''
          OutputData: [200x1 double]
                   y: 'Same as OutputData'
          OutputName: {'Speed'}
          OutputUnit: {'m/s'}
           InputData: [200x2 double]
                   u: 'Same as InputData'
           InputName: {2x1 cell}
           InputUnit: {2x1 cell}
              Period: [2x1 double]
         InterSample: {2x1 cell}
                  Ts: 0.1000
              Tstart: 0.1000
    SamplingInstants: [200x1 double]
            TimeUnit: 'seconds'
      ExperimentName: 'Exp1'
               Notes: {}
            UserData: []

前述のプロパティに加えて、'Period' は入力の周期を表し、Period = inf の場合は非周期的入力であることを意味します。

z.Period

ans =

   Inf
   Inf

入力のサンプル間動作は、'zoh' (ゼロ次ホールド、すなわち、区分的定数) または'foh' (1 次ホールド、すなわち、区分的線形) として得られます。同定ルーチンでは、この情報を使ってアルゴリズムを調整します。

z.InterSample

ans =

  2x1 cell array

    {'zoh'}
    {'zoh'}

チャネルは (入力と出力の両方)、"水平連結" (すなわち、z = [z1 z2]) で追加できます。

z2 = iddata(rand(200,1),ones(200,1),0.1,'OutputName','New Output',...
    'InputName','New Input');
z3 = [z,z2]

z3 =

Time domain data set with 200 samples.
Sample time: 0.1 seconds               
                                       
Outputs          Unit (if specified)   
   Speed            m/s                
   New Output                          
                                       
Inputs           Unit (if specified)   
   Voltage          Volt               
   Current          Ampere             
   New Input

z3 のいくつかのチャネルをプロットしてみましょう。

plot(z3(:,1,1)) % Data subset with Input 2 and Output 1.

plot(z3(:,2,3)) % Data subset with Input 2 and Output 3.

入力の生成

コマンド idinput は、一般的な入力信号を生成します。

u = idinput([30 1 10],'sine'); % 10 periods of 30 samples
u = iddata([],u,1,'Period',30) % Making the input an IDDATA object.

u =

Time domain data set with 300 samples.
Sample time: 1 seconds                 
                                       
Inputs       Unit (if specified)       
   u1

iddata の入力が SIM に適用され、iddata の出力が得られます。sim を使用して、入力 u による推定モデル m の応答を取得します。また、モデルのノイズダイナミクスに従ってノイズをモデルの応答に追加します。この操作は "AddNoise" シミュレーションオプションを使用して行います。

m = idpoly([1 -1.5 0.7],[0 1 0.5]);  % This creates a model; see below.
options = simOptions;
options.AddNoise = true;
y = sim(m,u,options) % simulated response produced as an iddata object

y =

Time domain data set with 300 samples.
Sample time: 1 seconds                 
Name: m                                
                                       
Outputs      Unit (if specified)       
   y1

シミュレーション入力 u および出力 y は、以下のように 1 つの iddata オブジェクトに結合できます。

z5 = [y u] % The output-input iddata.

z5 =

Time domain data set with 300 samples.
Sample time: 1 seconds                 
Name: m                                
                                       
Outputs      Unit (if specified)       
   y1                                  
                                       
Inputs       Unit (if specified)       
   u1

iddata オブジェクトの詳細については、help iddata に記載されています。

線形モデルオブジェクト

すべてのモデルは、MATLAB® オブジェクトとして提供されます。使用するモデルの種類に応じて、いくつか異なるオブジェクトがありますが、ほとんどが透過的です。

load iddata1
m = armax(z1,[2 2 2 1]);  % This creates an ARMAX model, delivered as an IDPOLY object

このモデルに関連するプロパティはすべて 1 つのオブジェクト (ここでは idpoly) にパッケージ化されています。詳細を表示するには、その名前を入力します。

m =
Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t)
  A(z) = 1 - 1.531 z^-1 + 0.7293 z^-2                    
                                                         
  B(z) = 0.943 z^-1 + 0.5224 z^-2                        
                                                         
  C(z) = 1 - 1.059 z^-1 + 0.1968 z^-2                    
                                                         
Sample time: 0.1 seconds
  
Parameterization:
   Polynomial orders:   na=2   nb=2   nc=2   nk=1
   Number of free coefficients: 6
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using ARMAX on time domain data "z1".  
Fit to estimation data: 76.38% (prediction focus)
FPE: 1.127, MSE: 1.082

    Many of the model properties are directly accessible

m.a    % The A-polynomial

ans =

    1.0000   -1.5312    0.7293

プロパティのリストは、以下のように取得します。

get(m)

                 A: [1 -1.5312 0.7293]
                 B: [0 0.9430 0.5224]
                 C: [1 -1.0587 0.1968]
                 D: 1
                 F: 1
    IntegrateNoise: 0
          Variable: 'z^-1'
           IODelay: 0
         Structure: [1x1 pmodel.polynomial]
     NoiseVariance: 1.1045
        InputDelay: 0
       OutputDelay: 0
         InputName: {'u1'}
         InputUnit: {''}
        InputGroup: [1x1 struct]
        OutputName: {'y1'}
        OutputUnit: {''}
       OutputGroup: [1x1 struct]
             Notes: [0x1 string]
          UserData: []
              Name: ''
                Ts: 0.1000
          TimeUnit: 'seconds'
      SamplingGrid: [1x1 struct]
            Report: [1x1 idresults.polyest]

推定されたパラメーターの共分散を各パラメーターについて標準偏差 +/- 1 の不確かさの値で表示するには、present を使用します。

present(m)

                                                                              
m =                                                                           
Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t)                     
  A(z) = 1 - 1.531 (+/- 0.01801) z^-1 + 0.7293 (+/- 0.01473) z^-2             
                                                                              
  B(z) = 0.943 (+/- 0.06074) z^-1 + 0.5224 (+/- 0.07818) z^-2                 
                                                                              
  C(z) = 1 - 1.059 (+/- 0.06067) z^-1 + 0.1968 (+/- 0.05957) z^-2             
                                                                              
Sample time: 0.1 seconds                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   na=2   nb=2   nc=2   nk=1                             
   Number of free coefficients: 6                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Termination condition: Near (local) minimum, (norm(g) < tol)..                
Number of iterations: 3, Number of function evaluations: 7                    
                                                                              
Estimated using ARMAX on time domain data "z1".                               
Fit to estimation data: 76.38% (prediction focus)                             
FPE: 1.127, MSE: 1.082                                                        
More information in model's "Report" property.

すべてのモデルパラメーター (または自由モデルパラメーターのみ) と、各モデルパラメーターの不確かさのフラットリストを取得するには、getpvec を使用します。共分散行列全体を取得するには、getcov を使用します。

[par, dpar] = getpvec(m, 'free')
CovFree = getcov(m,'value')

par =

   -1.5312
    0.7293
    0.9430
    0.5224
   -1.0587
    0.1968


dpar =

    0.0180
    0.0147
    0.0607
    0.0782
    0.0607
    0.0596


CovFree =

    0.0003   -0.0003    0.0000    0.0007    0.0004   -0.0003
   -0.0003    0.0002   -0.0000   -0.0004   -0.0003    0.0002
    0.0000   -0.0000    0.0037   -0.0034   -0.0000    0.0001
    0.0007   -0.0004   -0.0034    0.0061    0.0008   -0.0005
    0.0004   -0.0003   -0.0000    0.0008    0.0037   -0.0032
   -0.0003    0.0002    0.0001   -0.0005   -0.0032    0.0035

nf = 0、nd = 0 は、一般的な線形モデルの次数を示し、ARMAX モデルの場合は特殊例となります。

レポートには、推定プロセスに関する情報が含まれます。

m.Report
m.Report.DataUsed       % record of data used for estimation
m.Report.Fit            % quantitative measures of model quality
m.Report.Termination    % search termination conditions

ans = 

              Status: 'Estimated using ARMAX with prediction focus'
              Method: 'ARMAX'
    InitialCondition: 'zero'
                 Fit: [1x1 struct]
          Parameters: [1x1 struct]
         OptionsUsed: [1x1 idoptions.polyest]
           RandState: [1x1 struct]
            DataUsed: [1x1 struct]
         Termination: [1x1 struct]


ans = 

  struct with fields:

            Name: 'z1'
            Type: 'Time domain data'
          Length: 300
              Ts: 0.1000
     InterSample: 'zoh'
     InputOffset: []
    OutputOffset: []


ans = 

  struct with fields:

    FitPercent: 76.3807
       LossFcn: 1.0824
           MSE: 1.0824
           FPE: 1.1266
           AIC: 887.1256
          AICc: 887.4123
          nAIC: 0.1192
           BIC: 909.3483


ans = 

  struct with fields:

                 WhyStop: 'Near (local) minimum, (norm(g) < tol).'
              Iterations: 3
    FirstOrderOptimality: 7.2436
                FcnCount: 7
              UpdateNorm: 0.0067
         LastImprovement: 0.0067

最小化に関するオンライン情報を入手するには、"Display" 推定オプションを使用します (指定可能な値は "off"、"on" および "full")。これにより、モデル推定の進捗状況を示す進捗状況ビューアーが表示されます。

Opt = armaxOptions('Display','on');
m1 = armax(z1,[2 2 2 1],Opt);

線形モデルのバリアント - IDTF、IDPOLY、IDPROC、IDSS および IDGREY

線形モデルにはいくつかの種類があります。上記に多項式タイプのモデルで idpoly バージョンの例を示します。Box-Jenkins モデル、出力誤差モデル、ARMAX モデルなど、多項式タイプのモデルの各バリアントは、対応する推定器 bj, oe, armax, arx などを使用して取得します。これらはすべて、idpoly オブジェクトとして表されます。

その他のバリアントは、状態空間モデルの idss、ユーザー定義の構造化状態空間モデルの idgrey、伝達関数モデルの idtf、プロセスモデルの idproc です (ゲイン + 遅延 + 静的ゲイン)。

モデルを評価するコマンド bode, step, iopzmap, compare などはすべて、たとえば以下のように、モデルオブジェクトで直接機能します。

compare(z1,m1)

状態/空間への変換、伝達関数と零点/極は、idssdata、tfdata および zpkdata で取得します。

[num,den]  = tfdata(m1,'v')

num =

         0    0.9430    0.5224


den =

    1.0000   -1.5312    0.7293

'v' は、cell 配列ではなくベクトルとして、num と den が返されることを意味します。cell 配列は、多変数システムを処理する場合に便利です。また、num と den の値で標準偏差 1 の不確かさを取得するには、次のように記述します。

[num, den, ~, dnum, dden] = tfdata(m1,'v')

num =

         0    0.9430    0.5224


den =

    1.0000   -1.5312    0.7293


dnum =

         0    0.0607    0.0782


dden =

         0    0.0180    0.0147

同定されたモデルを Control System Toolbox の数値 LTI に変換

オブジェクトは、tf、ss および zpk のように、Control System Toolbox™ モデルオブジェクトに直結しています。また、Control System Toolbox が使用可能であれば、これらの LTI オブジェクトに変換できます。たとえば、tf は、idpoly オブジェクトを tf オブジェクトに変換します。

CSTBInstalled = exist('tf','class')==8;
if CSTBInstalled % check if Control System Toolbox is installed
    tfm = tf(m1) % convert IDPOLY model m1 into a TF object
end

tfm =
 
  From input "u1" to output "y1":
    0.943 z^-1 + 0.5224 z^-2
  ----------------------------
  1 - 1.531 z^-1 + 0.7293 z^-2
 
Sample time: 0.1 seconds
Discrete-time transfer function.

IDLTI モデルを Control Systems Toolbox の LTI モデルに変換する場合、ノイズ成分は保持されません。ノイズチャネルを LTI モデルの標準入力として含めるには、"augmented" フラグを使用します。

if CSTBInstalled
    tfm2 = tf(m1,'augmented')
end

tfm2 =
 
  From input "u1" to output "y1":
    0.943 z^-1 + 0.5224 z^-2
  ----------------------------
  1 - 1.531 z^-1 + 0.7293 z^-2
 
  From input "v@y1" to output "y1":
  1.051 - 1.113 z^-1 + 0.2069 z^-2
  --------------------------------
    1 - 1.531 z^-1 + 0.7293 z^-2
 
Input groups:           
      Name      Channels
    Measured       1    
     Noise         2    
                        
Sample time: 0.1 seconds
Discrete-time transfer function.

モデル tfm2 ではノイズチャネルは v@y1 という名前になっています。

System Identification Toolbox のデータおよびモデル オブジェクト

IDDATA オブジェクト

線形モデル オブジェクト

線形モデルのバリアント - IDTF、IDPOLY、IDPROC、IDSS および IDGREY

同定されたモデルを Control System Toolbox の数値 LTI に変換

System Identification Toolbox のデータおよびモデルオブジェクト

線形モデルオブジェクト