最適化関数を利用したシステム同定

Global Optimization Toolbox の産業応用

このワークシートでは、システム同定による自動車のサスペンションシステムのチューニングを行います。

路上の凹凸に応じてある指定された振る舞いを示すサスペンションを設計したいと考えます。
問題の変数は、
バネ定数 K  と  ダンパ変数 B  です。 質量に M  を与え、
標準的な凹凸の期待される振幅に対して、目的の反応と一致する反応のシステムを作成するために
 バネ定数 K  と ダンパ変数 B  の値を見つけます。

>    restart:

目的とするデータの作成

目的とする時系列のデータが無いため、下のような応答を示す関数からデータを作成します。
この関数は指数関数で作成されており、衝撃の振幅が 0.1 メートルに達する場合、その振幅は指数関数的に減衰します。

>    Target := t->(0.14*exp(-4.9*t)*cos(5*t -.7754));

Target := proc (t) options operator, arrow; .14*exp(-4.9*t)*cos(5*t-.7754) end proc

グラフ表示して反応を確かめます。

>    plot( Target(t),t=0..2,title="目的の反応",labels=["秒","メーター"],axes=frame);

[Maple Plot]

上記の反応の  0.5 秒毎 のデータを取り、目的のデータとして利用します。

>    time_sample:=[seq(i/2,i=1..4)];

time_sample := [1/2, 1, 3/2, 2]

それではグラフ表示するために、時間軸データとして作成します。

>    xydata:=ListTools[Transpose]([time_sample,[seq(Target(t),t=time_sample)]]);

xydata := [[1/2, -.1850800734e-2], [1, -.4886026617e-3], [3/2, .8133982026e-4], [2, -.7608201709e-5]]

データをグラフ表示します。

>    plot([xydata,Target(t)],t=0..2,style=[point,line],color=[black,red],title="目的の反応と使用するデータ",labels=["秒","メーター"],axes=frame);

[Maple Plot]

このデータと一致する反応をシステム同定を行って作り出すことが目的です。

システム同定

マス・バネ・ダンパ系のモデル式を作成します。
この式は減衰振動方程式と言われ、時間が経過するにつれて減衰していく微分方程式です。

>    suspension:=M*diff(x(t),t,t)+B*diff(x(t),t)+K*x(t)=0;

suspension := M*diff(x(t),`$`(t,2))+B*diff(x(t),t)+K*x(t) = 0

質量 M  に 450  と与えます。

>    M:=450:

それでは初期値に x(0) = 0.1, x'(0)=0  を与えこの式の解を求めます。

>    ans:=dsolve({suspension,x(0)=0.1,D(x)(0)=0});

ans := x(t) = 1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp((-1/900*B+1/900*(B^2-1800*K)^(1/2))*t)+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp((-1/900*B-1/900*(B^2-1800*K)^(1/2))*...
ans := x(t) = 1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp((-1/900*B+1/900*(B^2-1800*K)^(1/2))*t)+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp((-1/900*B-1/900*(B^2-1800*K)^(1/2))*...

得られた微分方程式の解を B,K,t の関数として定義します。

>    F:=unapply(rhs(ans),B,K,t):

B,K に適当な値を値を代入して、グラフ上で目的の反応と比較してみます。

>    plot([xydata,F(1000,10000,t),Target(t)],t=0..2,style=[point,line,line],color=[black,blue,red],title="目的の反応と B=1000,K=10000 の時の応答比較",labels=["秒","メーター"],axes=frame);

[Maple Plot]

>    plot([xydata,F(2000,20000,t),Target(t)],t=0..2,style=[point,line,line],color=[black,blue,red],title="目的の反応と B=2000,K=20000 の時の応答比較",labels=["秒","メーター"],axes=frame);

[Maple Plot]

今回は B  と K  にそれぞれ 1000  と 10000  等の適当な値を代入してグラフ化しましたが、グラフを確認すると適切なパラメータを代入できていると言えません。
真の目標は、目的の反応と
同じ応答を示す B,K の値を見つけることです

>    Error:=add((F(B,K,time_sample[i])-Target(time_sample[i]))^2,i=1..nops(time_sample));

Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...
Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...
Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...
Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...
Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...
Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...
Error := (1/20*(B*(B^2-1800*K)^(1/2)+B^2-1800*K)/(B^2-1800*K)*exp(-1/1800*B+1/1800*(B^2-1800*K)^(1/2))+1/20*(B^2-B*(B^2-1800*K)^(1/2)-1800*K)/(B^2-1800*K)*exp(-1/1800*B-1/1800*(B^2-1800*K)^(1/2))+.1850...

この2乗差の総和の式は、非常に複雑で非線形な式ですが、 この式を評価関数として GlobalSolve 関数を利用し最適化を行い、総和の値が最小となるような B,K  の値を発見します。

>    sol := GlobalOptimization[GlobalSolve]( Re(Error),B=100..10000,K=100..40000);

sol := [.412502906919710016e-7, [B = 4531.02176962435806, K = 22691.7138449499362]]

2乗差の総和が 0 に非常に近くなっており、その値となる B,K も求められました。

それでは理想のデータを黒い点、目的の反応の式を赤線、GlobalSolve で求められた式を緑線でグラフ表示します。

>    plot([xydata,subs({op(sol[2])},F(B,K,t)),Target(t)], t=0..2,
        labels=["秒","メーター"], title = "目的と実際の反応", thickness=1, color=[black,green,red],style=[point,line,line],axes=frame);

[Maple Plot]

赤線がほとんど目立ちませんので、非常によく再現できていることが視覚的にもわかります。

それでは下に、目的の反応の式と、GlobalSolve で求めた式の差の絶対値を取り、その差も確認してみます。

>    plot(abs(Target(t) - subs({op(sol[2])},F(B,K,t))), t=0..2,
        title="`目的の反応 - 実際の反応`の値の絶対値", labels=["秒","メーター"], thickness=2);

[Maple Plot]

非常に小さな値ですが最も差が大きい所は 0.2 秒あたりであることが確認できます。
これを改善したい場合、もっと細かくデータを取り、もう一度 GlobalSolve で最適化を行うことも有効かもしれませんので、試してみましょう。

より精度よく再現するための再チューニング

グラフの曲率が高い部分、 0 秒から 1 秒までを 0.05 秒づつデータを取り、その他に 1.2 秒と 2 秒の時のデータを利用します

>    time_sample:=[seq(i/20,i=1..20),seq(i/5,i=[6,10])];

time_sample := [1/20, 1/10, 3/20, 1/5, 1/4, 3/10, 7/20, 2/5, 9/20, 1/2, 11/20, 3/5, 13/20, 7/10, 3/4, 4/5, 17/20, 9/10, 19/20, 1, 6/5, 2]

>    xydata:=ListTools[Transpose]([time_sample,[seq(Target(t),t=time_sample)]]):

作成したデータ点をグラフ表示します。

>    plot([xydata,Target(t)],t=0..2,style=[point,line],color=[black,red],title="目的の反応と使用するデータ",labels=["秒","メーター"],axes=frame);

[Maple Plot]

それでは、評価関数を作成します。(評価関数はより複雑なため、表示はしません)

>    Error:=add((F(B,K,time_sample[i])-Target(time_sample[i]))^2,i=1..nops(time_sample)):

GlobalSolve 関数で最適化を行います。オプションで、初期値では一回目のGlobalSolve で得られた結果を利用し、
B に 4531、K に 22691  を使い、計算時間の制限を1000秒と指定します。

>    sol := GlobalOptimization[GlobalSolve](Re(Error),B=100..10000,K=1000..40000,initialpoint=[B=4531,K=22691],timelimit=1000);

sol := [.410973146301959979e-9, [B = 4411.71231561056902, K = 22062.5862395228141]]

誤差の総和が前回実行時よりも小さくなっており、より素晴らしい値が求まっていることが期待できます。
それではグラフ表示して確認してみましょう。

>    plot([xydata,subs({op(sol[2])},F(B,K,t)),Target(t)], t=0..2,
        labels=["秒","メーター"], title = "目的と実際の反応", thickness=1, color=[black,green,red],style=[point,line,line],axes=frame);

[Maple Plot]

下記で誤差の絶対値をグラフ化しその差を確認します。

>    plot(abs(Target(t) - subs({op(sol[2])},F(B,K,t))), t=0..2,
        title="`目的の反応 - 実際の反応`の値の絶対値", labels=["秒","メーター"], thickness=2);

[Maple Plot]

結論:
Maple と GlobalOptimization Toolbox を用いることで、モデル式からデータとのシステム同定を行いました。
この GlobalSolve 関数の力は如何でしょうか?大変精度よく、システム同定を行うことができているといえます。
今回はデータの与え方や初期値を利用することで、より精度が高くなるということを示しましたが、
その他のオプションも利用することで、より最適なパラメータを導くことができるかもしれません。