024_DynamicSystem.mw

ダイナミックシステムの解析と設計

イントロダクション

このワークシートでは、Maple で行うダイナミックシステムの設計と解析例を紹介します。

システムの定義

[Inserted Image]

> restart;

> with(plots):
with(inttrans):

Warning, the name changecoords has been redefined

線形システム解析

フォワードパス伝達関数を定義します。

> G := K/(s*(M*s+C));

G := K/(s*(M*s+C))

基準信号(ステップ関数)

> R := 3/s;

R := 3/s

閉ループ伝達関数で計算し、式を整理します。

> X := R*(G/(1+G));
X := simplify(X);

X := 3*K/(s^2*(M*s+C)*(1+K/(s*(M*s+C))))

X := 3*K/(s*(M*s^2+s*C+K))

逆ラプラス変換を計算し、過渡応答を求めます。

> xt := simplify(invlaplace(X,s,t));

xt := -3*(-(C^2-4*M*K)^(1/2)+cosh(1/2*t*(C^2-4*M*K)^(1/2)/M)*exp(-1/2*t*C/M)*(C^2-4*M*K)^(1/2)+C*sinh(1/2*t*(C^2-4*M*K)^(1/2)/M)*exp(-1/2*t*C/M))/(C^2-4*M*K)^(1/2)

値を代入します。 M = 1, C = 1.2, K = 3.

> x1 := subs(M=1,C=1.2,K=3,xt);

x1 := .9231861822*I*(-3.249615362*I+3.249615362*I*cosh(1.624807681*I*t)*exp(-.6000000000*t)+1.2*sinh(1.624807681*I*t)*exp(-.6000000000*t))

応答をプロットします。

> linresp := plot(x1, t=0..10, axes=boxed, color=black): display(linresp);

[Plot]

値を代入します (M=1, K=3)。今度は x(t,C) です。サーフェイスをプロットします。

> x2 := subs(M=1, K=3, xt);

x2 := -3*(-(C^2-12)^(1/2)+cosh(1/2*t*(C^2-12)^(1/2))*exp(-1/2*t*C)*(C^2-12)^(1/2)+C*sinh(1/2*t*(C^2-12)^(1/2))*exp(-1/2*t*C))/(C^2-12)^(1/2)

> plot3d(x2, C=.1..3, t=0..10, axes=boxed);

[Plot]

値を代入します (M=1, C=1.2)。今度は x(t,K) です。サーフェイスをプロットします。

> x3 := subs(M=1, C=1.2, xt);

x3 := -3*(-(1.44-4*K)^(1/2)+cosh(1/2*t*(1.44-4*K)^(1/2))*exp(-.6000000000*t)*(1.44-4*K)^(1/2)+1.2*sinh(1/2*t*(1.44-4*K)^(1/2))*exp(-.6000000000*t))/(1.44-4*K)^(1/2)

> plot3d(x3, K=.5..3, t=0..10, axes=boxed);

[Plot]

非線形ばねモデル

非線形ばねの解析モデル K(w) = K1w + K2w3 を考えます。

支配微分方程式

> gov := M*diff(x(t),t,t) + C*diff(x(t),t) + stiff = ref;

gov := M*(diff(x(t), `$`(t, 2)))+C*(diff(x(t), t))+stiff = ref

非線形のばねを定義します。

> nlstiff := K1*w - K2*w^3;

nlstiff := K1*w-K2*w^3

初期条件を定義します。

> inits:= x(0)=0,D(x)(0)=0 ;

inits := x(0) = 0, D(x)(0) = 0

非線形部分、基準信号 (r(t)=3H(t))、パラメータ値を代入します。

> de1:=subs(stiff=subs(w=x(t),nlstiff), ref=subs(w=xref(t),nlstiff),gov):
de2:=subs(xref(t)=mag*Heaviside(t),de1):
de3:=subs(M=1,C=1.2,K1=3,K2=0.05,mag=3,de2);

de3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+3*x(t)-0.5e-1*x(t)^3 = 9*Heaviside(t)-1.35*Heaviside(t)^3

物理的に t はつねに 0 または 正なので、Heavisede ステップ関数は 1 で置き換えることができます。

> de4 := subs(Heaviside(t)=1, de3):

数値的に非線形微分方程式を解きます。

> sol2:=dsolve({de4,inits}, x(t), type=numeric, method=classical[rk4], stepsize=.05);

sol2 := proc (x_classical) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14...

> nlresp:=plots[odeplot](sol2, [t, x(t)], 0..10, color=blue, axes=boxed):
plots[display](nlresp);

[Plot]

実験データに基づく非線形モデル

実験データを定義します。

> datax:=[0,.5,1.1,1.88,2.08,2.65,3.9,5.1];
datay:=[0,2.1,4.3,6.8,7.1,8.29,9.7,10.3];

datax := [0, .5, 1.1, 1.88, 2.08, 2.65, 3.9, 5.1]

datay := [0, 2.1, 4.3, 6.8, 7.1, 8.29, 9.7, 10.3]

> pldata := zip((x,y)->[x,y], datax, datay):
dataplot := pointplot(pldata, symbol=BOX):

キュービックスプラインを計算し、プロットします。

> stiff3:=spline(datax, datay, w, cubic);

stiff3 := PIECEWISE([4.34588552399999984*w-.583542095700000041*w^3, w < .5], [.145885524+3.90822895300000006*w-.875313143513378856*(w-.5)^2+.787848890699999970*(w-.5)^3, w < 1.1], [.220397020+3.708729...

> stplot3:=plot(stiff3, w=0..5, color=red):

> display(dataplot, stplot3, axes=boxed);

[Plot]

スプライン、基準信号、パラメータ値を 支配方程式に代入します。

> dex1:=subs(stiff=subs(w=x(t), stiff3), ref=subs(w=xref(t),stiff3),gov):
dex2:=subs(xref(t)=mag*Heaviside(t), dex1):
dex3:=subs(M=1, C=1.2, K1=3, K2=0.05, mag=3, dex2);

dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...dex3 := (diff(x(t), `$`(t, 2)))+1.2*(diff(x(t), t))+(PIECEWISE([4.34588552399999984*x(t)-.583542095700000041*x(t)^3, x(t) < .5], [.145885524+3.90822895300000006*x(t)-.875313143513378856*(x(t)-.5)^2+.7...

数値的に解を求め、線形モデルと非線形解析モデルとともにプロットします。上と同様に、t は 0 または正の値なので、Heaviside 関数を 1 に置き換えます。

> dex4:=subs(Heaviside(t)=1, %):

> sol3:=dsolve({dex4,inits}, x(t), type=numeric, method=classical[rk4], stepsize=.05);

sol3 := proc (x_classical) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14...

> nldatresp :=plots[odeplot](sol3, [t,x(t)], 0..10, color=red):

> display({linresp, nlresp, nldatresp});

[Plot]

剛性関数を比較します。

> stplot1:=plot(subs(K1=3,K2=.05,nlstiff), w=0..5, color=blue):

> stplot2:=plot(3*w, w=0..5, color=black):

> display(stplot1, stplot2, stplot3, axes=boxed);

[Plot]