このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
2 つの剛体間の乾燥摩擦: 複数の実験データを使用したパラメーター推定
この例では、複数の実験データを使用して非線形グレー ボックス モデルのパラメーターを推定する方法を説明します。2 つの剛体間の乾燥摩擦を示すシステムを議論の基礎として使用します。このシステムでは、1 つの剛体は固定され、もう 1 つは固定された剛体上を、図 1 に示すように外力に応じて前後に移動します。
図 1: 2 つの剛体システムの概略図
2 つの剛体間の乾燥摩擦のモデル化
ニュートンの第 3 法則を使用して、移動する剛体の動きを次のように記述します。
F_tot(t) = m*a(t) = m*d/dt v(t) = m*d^2/dt^2 s(t)
ここで F_tot(t) は外力 F(t) から 2 つの剛体の接触による摩擦力を差し引いたものと等しくなります。摩擦力は、すべり摩擦力 F_s(t) と乾燥摩擦力 F_d(t) との合計と仮定されます。前者は通常は速度の線形関数 F_s(t) = k*v(t) としてモデル化されます。ここで、k は不明なすべり摩擦パラメーターです。もう一方の乾燥摩擦はより複雑な現象です。次の文献
A. Clavel、M. Sorine、Q. Zhang 著『Modeling and identification of a leaf spring system』、In third IFAC Workshop on Advances in Automotive Control、2001年
常微分方程式でモデル化されます。
d/dt F_d(t) = -1/e*|v(t)|*F_d(t) + f/e*v(t)
ここで e と f はそれぞれ距離と力の 2 つの不明パラメーターです。
入力信号を u(t) = F(t) [N] として、状態を次のように示します。
x1(t) = s(t) Position of the moving body [m]. x2(t) = v(t) Velocity of the moving body [m/s]. x3(t) = F_d(t) Dry friction force between the bodies [N].
モデル パラメーターは次のとおりです。
m Mass of the moving body [m]. k Sliding friction force coefficient [kg/s]. e Distance-related dry friction parameter [m]. f Force-related dry friction parameter [N].
次の状態空間モデル構造を得られます。
d/dt x1(t) = x2(t) d/dt x2(t) = 1/m*(u(t) - k*x2(t) - x3(t)) d/dt x3(t) = 1/e*(-|x2(t)|*x3(t) + f*x2(t))
y(t) = x1(t)
これらの方程式は C MEX モデル ファイル twobodies_c.c に入力されます。この状態と出力更新方程式の compute_dx と compute_y は次のようになります。
/* State equations. */
void compute_dx(double *dx, double t, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* Retrieve model parameters. */
double *m, *k, *e, *f;
m = p[0]; /* Mass of the moving body. */
k = p[1]; /* Sliding friction force coefficient. */
e = p[2]; /* Distance-related dry friction parameter. */
f = p[3]; /* Force-related dry friction parameter. */
/* x[0]: Position. */
/* x[1]: Velocity. */
/* x[2]: Dry friction force. */
dx[0] = x[1];
dx[1] = (u[0]-k[0]*x[1]-x[2])/m[0];
dx[2] = (-fabs(x[1])*x[2]+f[0]*x[1])/e[0];
}
/* Output equation. */
void compute_y(double *y, double t, double *x, double *u, double **p,
const mxArray *auxvar)
{
/* y[0]: Position. */
y[0] = x[0];
}
モデル構造を記述するファイルを作成した次のステップは、モデル化の状況を反映する IDNLGREY オブジェクトを作成することです。また、入力、出力、状態の名前と単位、モデル構造のモデル パラメーターについての情報を追加します。ここで、パラメーターと InitialState はベクトルとして指定されます。これは、既定では、NLGREYEST が呼び出されたときにすべてのモデル パラメーターが推定され、初期状態ベクトルは推定されないということを意味します。
FileName = 'twobodies_c'; % File describing the model structure. Order = [1 1 3]; % Model orders [ny nu nx]. Parameters = [380; 2200; 0.00012; 1900]; % Initial parameter vector. InitialStates = [0; 0; 0]; % Initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Two body system', ... 'InputName', 'Exogenous force', ... 'InputUnit', 'N', ... 'OutputName', 'Position of moving body', ... 'OutputUnit', 'm', ... 'TimeUnit', 's'); nlgr = setinit(nlgr, 'Name', {'Position of moving body' ... 'Velocity of moving body' ... 'Dry friction force between the bodies'}); nlgr = setinit(nlgr, 'Unit', {'m' 'm/s' 'N'}); nlgr = setpar(nlgr, 'Name', {'Mass of the moving body' ... 'Sliding friction force coefficient', ... 'Distance-related dry friction parameter' ... 'Force-related dry friction parameter'}); nlgr = setpar(nlgr, 'Unit', {'m' 'kg/s' 'm' 'N'});
入出力データ
この時点で、使用できる入出力データ (シミュレーションを実行したもの) を読み込みます。ファイルには 3 つの異なる (シミュレーションを実行した) テスト実行から得られたデータが含まれ、それぞれサンプリング レート (Ts) 0.005 秒を使用して生成されたノイズで破損された入出力サンプルが 1000 個含まれています。入力 u(t) は移動する剛体に作用する外力 [N] です。実験では、入力は対称的なのこぎり波状信号で、波形繰り返し周波数はすべての実験について同一ですが、最大信号振幅はテスト実行によって異なります。出力 y(t) は移動する剛体の位置 [m] です (固定された剛体に対して相対的)。このモデル化の目的のため、3 つの異なる実験データをもつ 1 つの IDDATA オブジェクトを作成します。
load('twobodiesdata'); z = merge(iddata(y1, u1, 0.005), iddata(y2, u2, 0.005), iddata(y3, u3, 0.005)); z.Name = 'Two body system'; z.InputName = nlgr.InputName; z.InputUnit = nlgr.InputUnit; z.OutputName = nlgr.OutputName; z.OutputUnit = nlgr.OutputUnit; z.Tstart = 0; z.TimeUnit = nlgr.TimeUnit;
前方同定実験に使用される入出力データがプロット ウィンドウに表示されます。
figure('Name', [z.Name ': input-output data']); for i = 1:z.Ne zi = getexp(z, i); subplot(z.Ne, 2, 2*i-1); % Input. plot(zi.SamplingInstants, zi.InputData); title([z.ExperimentName{i} ': ' zi.InputName{1}]); if (i < z.Ne) xlabel(''); else xlabel([z.Domain ' (' zi.TimeUnit ')']); end axis('tight'); subplot(z.Ne, 2, 2*i); % Output. plot(zi.SamplingInstants, zi.OutputData); title([z.ExperimentName{i} ': ' zi.OutputName{1}]); if (i < z.Ne) xlabel(''); else xlabel([z.Domain ' (' zi.TimeUnit ')']); end axis('tight'); end
図 2: 2 剛体システムからの入出力データ
初期 2 剛体モデルの性能
4 つの不明パラメーターを推定する前に、初期パラメーター値を使用してシステムを推定します。既定の微分方程式ソルバー (ode45) を既定のソルバー オプションで使用します。2 つの入力引数のみで呼び出すと、COMPARE は初期状態ベクトル全体を推定します。この場合、実験ごとに 1 つずつです。それぞれが 3 行 1 列の状態ベクトルをもつ 3 つの実験で、合計 9 つの推定初期状態が得られます。シミュレーションの出力と実際の出力がプロット ウィンドウに表示され、一致は悪くありませんが、予期したほどはよくありません
clf compare(z, nlgr);
図 3: 初期の 2 つの剛体モデルの実際の出力とシミュレーションの出力の比較
パラメーター推定
一致を高めるため、次に 4 つのパラメーターを推定します。最初の実験と最後の実験のデータを推定フェーズで使用し、検証を実行するために 2 番目の実験は除外します。
opt = nlgreyestOptions('Display', 'on'); nlgr = nlgreyest(getexp(z, [1 3]), nlgr, opt);
推定された 2 剛体モデルの性能
推定されたモデルの性能を調査するため、最後にシミュレーションを実行します。初期状態構造アレイ nlgr を調整して、実験ごとに推定する状態を COMPARE で完全に指定できます。初期状態 x1(0) と x2(0) を実験 1 に、x2(0) を実験 2 に、x3(0) を実験 3 について推定する構造を定義して使用します。この修正により、測定出力とモデル出力との比較がプロット ウィンドウに表示されます。
nlgr.InitialStates = struct('Name', getinit(nlgr, 'Name'), ... 'Unit', getinit(nlgr, 'Unit'), ... 'Value' , zeros(1, 3), 'Minimum', -Inf(1, 3), ... 'Maximum', Inf(1, 3), 'Fixed', ... {[false false true]; [true false true]; [true true false]}); compare(z, nlgr, compareOptions('InitialCondition', 'model'));
図 4: 推定された 2 つの剛体モデルの実際の出力とシミュレーションの出力の比較
特に注目すべきは、2 番目の実験からのデータによる結果で、これはパラメーター推定には使用されていませんでした。真のシステムのダイナミクスはすべての実験について十分に明確なモデル化がされています。推定されたパラメーターは、実験データの生成に使用されたパラメーターと近いものです。
disp(' True Estimated parameter vector');True Estimated parameter vector
ptrue = [400; 2e3; 0.0001; 1700];
fprintf(' %10.5f %10.5f\n', [ptrue'; getpvec(nlgr)']); 400.00000 402.11398
2000.00000 1986.78107
0.00010 0.00010
1700.00000 1705.01863
最後に、PRESENT コマンドを使用して、推定モデルの追加情報を取得します。
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'twobodies_c' (MEX-file):
dx/dt = F(t, u(t), x(t), p1, ..., p4)
y(t) = H(t, u(t), x(t), p1, ..., p4) + e(t)
with 1 input(s), 3 state(s), 1 output(s), and 4 free parameter(s) (out of 4).
Inputs:
u(1) Exogenous force(t) [N]
States: Initial value
x(1) Position of moving body(t) [m] xinit@exp1 0 (estimated) in [-Inf, Inf]
xinit@exp2 0 (estimated) in [-Inf, Inf]
xinit@exp3 0 (fixed) in [-Inf, Inf]
x(2) Velocity of moving body(t) [m/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
xinit@exp2 0 (estimated) in [-Inf, Inf]
xinit@exp3 0 (fixed) in [-Inf, Inf]
x(3) Dry friction force between the bodies(t) [N] xinit@exp1 0 (fixed) in [-Inf, Inf]
xinit@exp2 0 (fixed) in [-Inf, Inf]
xinit@exp3 0 (estimated) in [-Inf, Inf]
Outputs:
y(1) Position of moving body(t) [m]
Parameters: Value
p1 Mass of the moving body [m] 402.114 (estimated) in [-Inf, Inf]
p2 Sliding friction force coefficient [kg/s] 1986.78 (estimated) in [-Inf, Inf]
p3 Distance-related dry friction parameter [m] 0.000104699 (estimated) in [-Inf, Inf]
p4 Force-related dry friction parameter [N] 1705.02 (estimated) in [-Inf, Inf]
Name: Two body system
Status:
Model modified after estimation.
More information in model's "Report" property.
Model Properties
まとめ
この例では、IDNLGREY モデリングを実行する際に複数の実験データを使用する方法を説明しました。任意の数の実験を利用でき、各実験について、NLGREYEST、COMPARE、PREDICT、SIM、RESID などのコマンドで推定する初期状態を完全に指定できます。
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)