[Inserted Image]

『Global Optimization Toolbox による最適化』

数式処理システム Maple で動く実践的最適化アプリケーション

Global Optimization による産業適用例

〜自動車サスペンションシステムのチューニング〜

> restart;

ここでは、自動車のサスペンション・システムの設計における、バネ定数、ダンパー定数などの最適値を探し出します。

理想的応答例

ある衝撃に対する応答が、以下のように指数関数を用いて減衰する形で表現できる場合を考えます。

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

> plot( Target(t), t = 0..1, title="Target Response", labels=["seconds","metres"], thickness=2 );

[Plot]

>

実応答の導出

いま考えているサスペンションシステムの数学モデルは、運動方程式により以下で表現できます。

> System_Equation := m*diff(x(t),t,t) + b*diff(x(t),t) + k*x(t) = 0;

System_Equation := m*(diff(x(t), `$`(t, 2)))+b*(diff(x(t), t))+k*x(t) = 0

ここで、m は各ホイール上での車の質量とします。係数 b 及び k によりシステムの振幅の大きさが左右されます。

なお、初期条件として x(0) = Float(1, -1)とします。

いま、m = 450kg の場合での k, b, t の関数としてシステムの応答を計算します。

> m := 450:

> sol := dsolve( {System_Equation, x(0)=.1, D(x)(0)=0} );

sol := x(t) = 1/20*(b*(b^2-1800*k)^(1/2)+b^2-1800*k)*exp((-1/900*b+1/900*(b^2-1800*k)^(1/2))*t)/(b^2-1800*k)+1/20*(b^2-b*(b^2-1800*k)^(1/2)-1800*k)*exp((-1/900*b-1/900*(b^2-1800*k)^(1/2))*t)/(b^2-1800...

これを k, b, t の関数として定義します。

> Actual := unapply(rhs(sol), k, b, t);

Actual := proc (k, b, t) options operator, arrow; 1/20*(b*(b^2-1800*k)^(1/2)+b^2-1800*k)*exp((-1/900*b+1/900*(b^2-1800*k)^(1/2))*t)/(b^2-1800*k)+1/20*(b^2-b*(b^2-1800*k)^(1/2)-1800*k)*exp((-1/900*b-1/...

もしも k = 10^4, b = 10^3 のとき、理想応答と実応答を比べてみると、以下のグラフのようになります。理想応答にはほど遠く、大きく振幅してしまっていることがわかります。

> plot(
 [ Target(t), Actual(10^4, 10^3, t) ],
 t=0..1,
 thickness=3, labels=["seconds","metres"],
 title="Target Response and Actual Response for k = 10^4 and b = 10^3",
 legend=["Target Response", "Actual Response"]
);

[Plot]

>

実応答と理想応答の差

原理的には、理想応答と実応答の差は平方和の積分値として測定できますが、実応答関数が複雑であるために、この計算を厳密に行うことは現実的ではありません。従って、積分値を計算するために適当なサンプリング時間値に対しての誤差値を定義し、その誤差値が最小となるような係数を考えることにします。ここでは、 t = 1/2, 1, 3/2, 2 でのサンプリングを用います。

> time_sample:= 1/2, 1, 3/2, 2;

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

理想応答と実応答の誤差は、二乗誤差として定義します。

> Error := add( (Actual(k, b, time_sample[i]) - Target(time_sample[i]))^2, i = 1..4);

Error := (1/20*(b*(b^2-1800*k)^(1/2)+b^2-1800*k)*exp(-1/1800*b+1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+1/20*(b^2-b*(b^2-1800*k)^(1/2)-1800*k)*exp(-1/1800*b-1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+0.18...Error := (1/20*(b*(b^2-1800*k)^(1/2)+b^2-1800*k)*exp(-1/1800*b+1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+1/20*(b^2-b*(b^2-1800*k)^(1/2)-1800*k)*exp(-1/1800*b-1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+0.18...Error := (1/20*(b*(b^2-1800*k)^(1/2)+b^2-1800*k)*exp(-1/1800*b+1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+1/20*(b^2-b*(b^2-1800*k)^(1/2)-1800*k)*exp(-1/1800*b-1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+0.18...Error := (1/20*(b*(b^2-1800*k)^(1/2)+b^2-1800*k)*exp(-1/1800*b+1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+1/20*(b^2-b*(b^2-1800*k)^(1/2)-1800*k)*exp(-1/1800*b-1/1800*(b^2-1800*k)^(1/2))/(b^2-1800*k)+0.18...

誤差関数を k, b を変数として曲面プロットしてみます。

> plot3d(
 Re(Error),
 b=100..5000, k=1000..25000,
 axes=boxed, shading=zhue,
 transparency=.35, style=patchnogrid,
 orientation=[75,75]
);

[Plot]

>

Global Optimizationによる誤差の最小化

Global Optimization Toolboxの機能を用いて、誤差関数 Error を最小にする係数を定めます。

> with( GlobalOptimization );

[GlobalSolve, Interactive]

> infolevel[GlobalOptimization] := 2:
st := time():
sol := GlobalSolve( Re(Error), b=100..10000, k=100..40000 );
time() - st;

GlobalSolve:, `calling NLP solver`
SolveGeneral:, `calling global optimization solver`

SolveGeneral:, `number of problem variables`, 2

SolveGeneral:, `number of nonlinear inequality constraints`, 0

SolveGeneral:, `number of nonlinear equality constraints`, 0

SolveGeneral:, `trying evalhf mode`

SolveGeneral:, `trying evalf mode`

LGORUN:, `total number of function evaluations`, 2044
LGORUN:, `runtime in external solver`, 31.

LGORUN:, `maximum constraint infeasibility`, 0.

sol := [0.412502906919710016e-7, [b = 4531.02176962435806, k = 22691.7138449499362]]

31.085

最小値付近での誤差関数の挙動を確認してみます。

> plot3d(
 Re(Error),
 b=4300..4600, k=21000..23000,
 axes=boxed, shading=zhue, transparency=.35,
 style=patchnogrid, orientation=[75,75]
);

[Plot]

>

最適値の適用ー理想応答にどれだけ近いか?

前節で得られた最適係数 b, k をサスペンションモデルに適用してみます。

> assign( sol[2] );

> plot(
 [Target(t), Actual(k, b, t)+.05, .05  ],
 t=0..1,
 labels=["seconds","metres"],
 title = "Target and Actual Responses Overlaid (with vertical offset)",
 thickness=2, color=[red,green,black]
);

[Plot]

ほぼ理想的な振舞いとなっていることがわかります。

>