最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

摩擦モデリング: 静的 SISO システムの MATLAB ファイル モデリング

この例では、MATLAB 関数を ODE ファイルとして使用した、静的な単入力、単出力のシステムのグレー ボックス モデリングを説明します。

システム同定では通常、動的モデルのパラメーターを同定します。ただし、静的モデルもそれ自身、またはより大規模なモデルのサブモデルとして、対象となります。後者の例を「産業用ロボット アーム」(idnlgreydemo13.m) の事例として検討します。ここでは、静的摩擦モデルがロボット アーム モデルの固定 (推定済み) コンポーネントとして導入されています。

ここでは、静摩擦現象を IDNLGREY モデルとして表し、対応する係数を推定する方法を説明します。

連続微分可能摩擦モデル

不連続および区分的連続摩擦モデルは、高性能の連続コントローラーにとって問題となることがよくあります。この事実に対して、Makkar、Dixon、Sawyer、Hu は、実際に発生する摩擦現象の大部分を表現する連続的で微分可能な摩擦モデルを提唱しました。新しい摩擦モデル構造は、次の文献で説明されています。

C. Makkar, W. E. Dixon, W. G. Sawyer, and G.Hu "A New Continuously

Differentiable Friction Model for Control Systems Design", IEEE(R)/ASME

International Conference on Advanced Intelligent Mechatronics,

Monterey, CA, 2005, pages 600-605.

また、ここで行う静的同定実験のベースにもなります。

Makkar などが提唱する摩擦モデルは、別の本体と接触する本体のすべり角度 v(t) と摩擦力 f(t) とを、静的な関係によって結び付けます。

f(t) = g(1)*(tanh(g(2)*v(t) - tanh(g(3)*v(t))

+ g(4)*tanh(g(5)*v(t)) + g(6)*v(t)

ここで、g(1)、...、g(6) 不明な正のパラメーターです。このモデル構造には、実世界のアプリケーションから発生する優れたプロパティがいくつか示されます。

  1. 摩擦モデルは原点に関して対称的である

  2. 静的な摩擦係数は g(1)+g(4) によって近似される

  3. 方程式の第 1 項 tanh(g(2)*v(t) - tanh(g(3)*v(t) は、いわゆるストライベック効果を捉えており、ここで摩擦項はかなり急速に減衰して原点付近ですべり速度が上昇する

  4. クーロン摩擦効果は項 g(4)*tanh(g(5)*v(t)) によってモデル化される

  5. 粘性摩擦の散逸は最後の項 g(6)*v(t) によって反映される

摩擦の全般について、および提唱されたモデル構造の詳細については、上記で挙げた文献を参照してください。

IDNLGREY 摩擦モデリング

次に、静的摩擦を記述する IDNLGREY モデル オブジェクトを作成します。通常どおり、IDNLGREY モデリング ファイルを作成することから始めます。ここで、MATLAB ファイル friction_m.m を次の内容で作成します。

function [dx, f] = friction_m(t, x, v, g, varargin)

%FRICTION_M Nonlinear friction model with Stribeck, Coulomb and viscous

% dissipation effects.

% Output equation.

f = g(1)*(tanh(g(2)*v(1))-tanh(g(3)*v(1))) ... % Stribeck effect.

+ g(4)*tanh(g(5)*v(1)) ... % Coulomb effect.

+ g(6)*v(1); % Viscous dissipation term.

% Static system; no states.

dx = [];

状態更新 dx は常にモデル ファイルから返され、静的モデリングの場合は空 ([]) でなければならないことに注意してください。

次のステップは、モデル ファイル、モデル順序についての情報、推測されたパラメーター ベクトルなどを、IDNLGREY コンストラクターへの入力引数として渡すことです。また、入力および出力、すべてのモデル パラメーターが正でなければならない状態の名前と単位を指定します。

FileName      = 'friction_m';          % File describing the model structure.
Order         = [1 1 0];               % Model orders [ny nu nx].
Parameters    = {[0.20; 90; 11; ...
                  0.12; 110; 0.015]};  % Initial parameters.
InitialStates = [];                    % Initial initial states.
Ts            = 0;                     % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts,    ...
                'Name', 'Static friction model',                   ...
                'InputName', 'Slip speed', 'InputUnit', 'm/s',     ...
                'OutputName', 'Friction force', 'OutputUnit', 'N', ...
                'TimeUnit', 's');
nlgr = setpar(nlgr, 'Minimum', {zeros(5, 1)});   % All parameters must be >= 0.

これらの処理の後、プロパティが次のような初期摩擦モデルが得られます。

present(nlgr);
                                                                              
nlgr =                                                                        
Continuous-time nonlinear static model defined by 'friction_m' (MATLAB file): 
                                                                              
    y(t) = H(t, u(t), p1) + e(t)                                              
                                                                              
 with 1 input(s), 0 state(s), 1 output(s), and 6 free parameter(s) (out of 6).
                                                                              
 Inputs:                                                                      
    u(1)  Slip speed(t) [m/s]                                                 
 Outputs:                                                                     
    y(1)  Friction force(t) [N]                                               
 Parameters:    Value                                                         
   p1(1)   p1       0.2      (estimated) in [0, Inf]                          
   p1(2)             90      (estimated) in [0, Inf]                          
   p1(3)             11      (estimated) in [0, Inf]                          
   p1(4)           0.12      (estimated) in [0, Inf]                          
   p1(5)            110      (estimated) in [0, Inf]                          
   p1(6)          0.015      (estimated) in [0, Inf]                          
                                                                              
Name: Static friction model                                                   
                                                                              
Status:                                                                       
Created by direct construction or transformation. Not estimated.              
More information in model's "Report" property.                                

同定実験では、摩擦モデル全体だけでなく、縮小した摩擦モデルの動作の調査も関心の対象になっています。縮小モデルとは、完全なモデルの 3 つの項のうちの 2 つの項をもつ摩擦モデルを指します。これを調べるため、完全なモデル構造のコピーを 3 つ作成し、各コピーで 2 つの項のみが使用されるようにパラメーター ベクトルを固定します。

nlgr1 = nlgr;
nlgr1.Name = 'No Striebeck term';
nlgr1 = setpar(nlgr1, 'Value', {[zeros(3, 1); Parameters{1}(4:6)]});
nlgr1 = setpar(nlgr1, 'Fixed', {[true(3, 1); false(3, 1)]});
nlgr2 = nlgr;
nlgr2.Name = 'No Coulombic term';
nlgr2 = setpar(nlgr2, 'Value', {[Parameters{1}(1:3); 0; 0; Parameters{1}(6)]});
nlgr2 = setpar(nlgr2, 'Fixed', {[false(3, 1); true(2, 1); false]});
nlgr3 = nlgr;
nlgr3.Name = 'No dissipation term';
nlgr3 = setpar(nlgr3, 'Value', {[Parameters{1}(1:5); 0]});
nlgr3 = setpar(nlgr3, 'Fixed', {[false(5, 1); true]});

入出力データ

入力すべり速度が -10 m/s から 10 m/s の範囲でランプ型の方法で変化する、2 種類の (シミュレーションを実行した) データセットを自由に使用できます。データを読み込んで同定実験用に 2 つの IDDATA オブジェクトを作成します。ze は推定用、zv は検証用です。

load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'frictiondata'));
ze = iddata(f1, v, 1, 'Name', 'Static friction system');
set(ze, 'InputName', 'Slip speed', 'InputUnit', 'm/s',  ...
        'OutputName', 'Friction force', 'OutputUnit', 'N', ...
        'Tstart', 0, 'TimeUnit', 's');
zv = iddata(f2, v, 1, 'Name', 'Static friction system');
set(zv, 'InputName', 'Slip speed', 'InputUnit', 'm/s',  ...
        'OutputName', 'Friction force', 'OutputUnit', 'N', ...
        'Tstart', 0, 'TimeUnit', 's');

推定に使用される入出力データがプロット ウィンドウに表示されます。

figure('Name', ze.Name);
plot(ze);

図 1: 摩擦を示すシステムからの入出力データ

初期摩擦モデルの性能

入出力データと 4 種類の初期摩擦モデルを使用できるようになり、次の疑問はこれらのモデルの妥当性です。COMPARE で実行されるシミュレーションによって、推定データセットについて調べてみます。

set(gcf,'DefaultLegendLocation','southeast');
compare(ze, nlgr, nlgr1, nlgr2, nlgr3);

図 2: 4 つの初期摩擦モデルの実際の出力とシミュレーションの出力の比較

パラメーター推定

初期モデルはいずれも、実際の出力を適切に記述できません。これを解決するため、4 つのすべてのモデル構造のモデル パラメーターを推定します。最大で 30 回の反復を実行し、許容誤差がほぼゼロ (現実には実データでは発生しない) の場合にのみ反復を停止するよう、すべての推定を設定します。これらの計算には時間がかかります。

opt = nlgreyestOptions('Display', 'on');
opt.SearchOptions.MaxIterations = 30;
opt.SearchOptions.FunctionTolerance = eps;
opt.EstimateCovariance = false; 

nlgr  = nlgreyest(nlgr,  ze, opt);
nlgr1 = nlgreyest(nlgr1, ze, opt);
nlgr2 = nlgreyest(nlgr2, ze, opt);
nlgr3 = nlgreyest(nlgr3, ze, opt);

推定された摩擦モデルの性能

COMPARE を使用して得られた 4 つのモデルの実際の出力をシミュレーションの出力と比較して、モデルの性能をもう一度確認します。ただし、今回は比較は検証データセット zv に基づきます。

compare(zv, nlgr, nlgr1, nlgr2, nlgr3);

図 3: 4 つの推定された摩擦モデルの実際の出力とシミュレーションの出力の比較

このシステムでは、完全なモデルのほうが縮小されたモデルよりも性能がよいことがはっきりとわかります。ただし、縮小されたモデルは、モデル化している効果を取り込むことができ、それぞれの場合、推定の一致率は高くなります。粘性損失項が省略されているモデルでは、適合が最低になっています。完全なモデルの高い一致率は当然のことで、モデル構造は実際の出力データの生成に使用されたものだからです。完全なモデルのパラメーターも実際のモデル出力の生成に使用されたパラメーターと非常に似通っています。

disp('   True       Estimated parameter vector');
   True       Estimated parameter vector
ptrue = [0.25; 100; 10; 0.1; 100; 0.01];
fprintf('   %7.3f    %7.3f\n', [ptrue'; getpvec(nlgr)']);
     0.250      0.249
   100.000    106.637
    10.000      9.978
     0.100      0.100
   100.000     87.699
     0.010      0.010

4 つのすべての摩擦モデルに適用された最終予測誤差 (FPE) 基準 (低いほうが良い) によって、完全な摩擦モデルが優位であることが示されます。

fpe(nlgr, nlgr1, nlgr2, nlgr3);
   1.0e-03 *

    0.0002    0.1665    0.1584    0.7931

動的システムについては、PE を使用して静的モデルの予測誤差を調べることもできます。4 つの完全な摩擦モデルについてこれを実行し、残差がランダムであることが最後にわかります。

pe(zv, nlgr);

図 4: 推定された完全な摩擦モデルで取得された予測誤差

さらに、完全な摩擦モデルの残差 (残り) を調べてランダムであることを確認します。

figure('Name', [nlgr.Name ': residuals of estimated IDNLGREY model']);
resid(zv, nlgr);

図 5: 推定された完全な摩擦モデルで取得された残差

静的モデルのステップ応答も計算してプロットできます。単位ステップを適用して、4 つのすべての推定された摩擦モデルについてこれを実行します。

figure('Name', [nlgr.Name ': step responses of estimated models']);
step(nlgr, nlgr1, nlgr2, nlgr3);
legend(nlgr.Name, nlgr1.Name, nlgr2.Name, nlgr3.Name, 'location', 'SouthEast');

図 6: 4 つの推定された摩擦モデルの単位ステップ応答

最後に、パラメーターの推定された標準偏差、損失関数など、完全な摩擦モデルのいくつかのプロパティを示します。

present(nlgr);
                                                                                                          
nlgr =                                                                                                    
Continuous-time nonlinear static model defined by 'friction_m' (MATLAB file):                             
                                                                                                          
    y(t) = H(t, u(t), p1) + e(t)                                                                          
                                                                                                          
 with 1 input(s), 0 state(s), 1 output(s), and 6 free parameter(s) (out of 6).                            
                                                                                                          
 Inputs:                                                                                                  
    u(1)  Slip speed(t) [m/s]                                                                             
 Outputs:                                                                                                 
    y(1)  Friction force(t) [N]                                                                           
 Parameters:    Value                                                                                     
   p1(1)   p1      0.249402      (estimated) in [0, Inf]                                                  
   p1(2)            106.637      (estimated) in [0, Inf]                                                  
   p1(3)            9.97835      (estimated) in [0, Inf]                                                  
   p1(4)          0.0999916      (estimated) in [0, Inf]                                                  
   p1(5)            87.6992      (estimated) in [0, Inf]                                                  
   p1(6)          0.0100019      (estimated) in [0, Inf]                                                  
                                                                                                          
Name: Static friction model                                                                               
                                                                                                          
Status:                                                                                                   
Termination condition: Change in parameters was less than the specified tolerance..                       
Number of iterations: 9, Number of function evaluations: 10                                               
                                                                                                          
Estimated using Solver: FixedStepDiscrete; Search: lsqnonlin on time domain data "Static friction system".
Fit to estimation data: 99.68%                                                                            
FPE: 2.428e-07, MSE: 2.414e-07                                                                            
More information in model's "Report" property.                                                            

まとめ

この例では、静的システムの IDNLGREY モデリングを実行する方法を説明しました。この手順は基本的に、動的システムのモデル化と同じです。