031_ControlBridge.mw

ローディングブリッジの制御装置

イントロダクション

ローディング・ブリッジは負荷をある地点から他の地点に運びます。台車の加速と減速は好ましくない振動を発生してしまいます。そのため負荷はさらに安定し難くなり、時間もかかってしまいます。負荷がある地点から他の地点へ素早く移動し、すぐに安定するように、制御装置を設計します。

[Inserted Image]

> restart;

問題の定義

ローディング・ブリッジのパラメータは以下の通りです。

x[K] - 台車位置 [ m ]m[K]
- 台車の質量 [ kg ]
F[K]
- 台車の駆動力 [ N ]
x[G]
,
z[G] - 負荷の位置 [ m ]m[G]
- 負荷の質量 [ m ]

r - ロープの長さ [ m ]
theta
- ロープの角度,

g - 重力による加速度 [
m/(s^2) ]

台車の駆動は一次遅れの装置のように振る舞い、その特性は増幅 K[M] [ N/V ] と時定数 T[M] [ s ]  により決まります。パラメータの値は以下のとおりです。

> val:= m[K]=1000, m[G]=4000, r=10, g=9.81, K[M]=100, T[M]=1:

運動の方程式の決定

台車の加速度は力を考慮することで得られます。 S  はロープを通して台車に伝播する力です。

> eq1 := m[K]*diff(x[K](t),t,t) = F[K](t) + S*sin(theta(t));

eq1 := m[K]*(diff(x[K](t), `$`(t, 2))) = F[K](t)+S*sin(theta(t))

負荷についての垂直方向と水平方向の力は以下の通りです。

> eq2 := m[G]*diff(x[G](t),t,t) = -S*sin(theta(t));

eq2 := m[G]*(diff(x[G](t), `$`(t, 2))) = -S*sin(theta(t))

> eq3 := m[G]*diff(z[G](t),t,t) = m[G]*g - S*cos(theta(t));

eq3 := m[G]*(diff(z[G](t), `$`(t, 2))) = m[G]*g-S*cos(theta(t))

負荷の座標は

> eq4 := x[G](t) = x[K](t) + r*sin(theta(t));

eq4 := x[G](t) = x[K](t)+r*sin(theta(t))

> eq5 := z[G](t) = r*cos(theta(t));

eq5 := z[G](t) = r*cos(theta(t))

これらの方程式を2回微分 (ロープの長さは一定と仮定) します。そして eq2 および eq3 の中へ入れます。これで以下が解ります。

> s1 := diff(eq4,t,t);

s1 := diff(x[G](t), `$`(t, 2)) = (diff(x[K](t), `$`(t, 2)))-r*sin(theta(t))*(diff(theta(t), t))^2+r*cos(theta(t))*(diff(theta(t), `$`(t, 2)))

> s2 := diff(eq5,t,t);

s2 := diff(z[G](t), `$`(t, 2)) = -r*cos(theta(t))*(diff(theta(t), t))^2-r*sin(theta(t))*(diff(theta(t), `$`(t, 2)))

> s3 := subs(s1, eq2);

s3 := m[G]*((diff(x[K](t), `$`(t, 2)))-r*sin(theta(t))*(diff(theta(t), t))^2+r*cos(theta(t))*(diff(theta(t), `$`(t, 2)))) = -S*sin(theta(t))

> s4 := subs(s2, eq3);

s4 := m[G]*(-r*cos(theta(t))*(diff(theta(t), t))^2-r*sin(theta(t))*(diff(theta(t), `$`(t, 2)))) = m[G]*g-S*cos(theta(t))

最後の方程式を解き S を求めます。そして s3 に代入します。

> s5:=S=solve(s4,S);

s5 := S = -(-m[G]*r*cos(theta(t))*(diff(theta(t), t))^2-m[G]*r*sin(theta(t))*(diff(theta(t), `$`(t, 2)))-m[G]*g)/cos(theta(t))

> tmp:=subs(s5,s3);

tmp := m[G]*((diff(x[K](t), `$`(t, 2)))-r*sin(theta(t))*(diff(theta(t), t))^2+r*cos(theta(t))*(diff(theta(t), `$`(t, 2)))) = (-m[G]*r*cos(theta(t))*(diff(theta(t), t))^2-m[G]*r*sin(theta(t))*(diff(the...tmp := m[G]*((diff(x[K](t), `$`(t, 2)))-r*sin(theta(t))*(diff(theta(t), t))^2+r*cos(theta(t))*(diff(theta(t), `$`(t, 2)))) = (-m[G]*r*cos(theta(t))*(diff(theta(t), t))^2-m[G]*r*sin(theta(t))*(diff(the...

この方程式に cos(theta(t))/m[G] をかけ、方程式の右辺を左辺に移項します。これで、位置 x[K] と角度 theta にのみ依存する微分方程式が得られます。

> tmp2 := expand((tmp*cos(theta(t)))/m[G]);

tmp2 := cos(theta(t))*(diff(x[K](t), `$`(t, 2)))-cos(theta(t))*r*sin(theta(t))*(diff(theta(t), t))^2+r*cos(theta(t))^2*(diff(theta(t), `$`(t, 2))) = -cos(theta(t))*r*sin(theta(t))*(diff(theta(t), t))^...tmp2 := cos(theta(t))*(diff(x[K](t), `$`(t, 2)))-cos(theta(t))*r*sin(theta(t))*(diff(theta(t), t))^2+r*cos(theta(t))^2*(diff(theta(t), `$`(t, 2))) = -cos(theta(t))*r*sin(theta(t))*(diff(theta(t), t))^...

> deq1:=simplify(lhs(tmp2)-rhs(tmp2))=0;

deq1 := cos(theta(t))*(diff(x[K](t), `$`(t, 2)))+r*(diff(theta(t), `$`(t, 2)))+sin(theta(t))*g = 0

方程式 eq2 を方程式 eq1 に代入します。その結果を方程式 s1 に代入します。 位置 x[K] と角度 theta にのみ依存する2番目の微分方程式を得ます。

> tmp4:=algsubs(rhs(eq2)=lhs(eq2),eq1);

tmp4 := m[K]*(diff(x[K](t), `$`(t, 2))) = F[K](t)-m[G]*(diff(x[G](t), `$`(t, 2)))

> deq2:=expand(subs(s1,tmp4));

deq2 := m[K]*(diff(x[K](t), `$`(t, 2))) = F[K](t)-m[G]*(diff(x[K](t), `$`(t, 2)))+m[G]*r*sin(theta(t))*(diff(theta(t), t))^2-m[G]*r*cos(theta(t))*(diff(theta(t), `$`(t, 2)))

sin(theta(t)) theta(t) および cos(theta(t)) を 1 に設定することで微分方程式を線形化します。

角速度は小さいので diff(theta(t), t) の2次の項は無視できます。

> lindeq1 := subs(sin(theta(t))=theta(t), cos(theta(t))=1, (diff(theta(t),t))^2=0, deq1);

lindeq1 := (diff(x[K](t), `$`(t, 2)))+r*(diff(theta(t), `$`(t, 2)))+theta(t)*g = 0

> lindeq2 := subs(sin(theta(t))=theta(t), cos(theta(t))=1, (diff(theta(t),t))^2=0, deq2);

lindeq2 := m[K]*(diff(x[K](t), `$`(t, 2))) = F[K](t)-m[G]*(diff(x[K](t), `$`(t, 2)))-m[G]*r*(diff(theta(t), `$`(t, 2)))

台車を駆動すると一次遅れの装置のように振る舞うので、次の微分方程式によって記述することができます。(ここで K[M] は増幅で T[M] は時定数です。)

> deq3 := T[M]*diff(F[K](t),t) + F[K](t) = K[M]*u;

deq3 := T[M]*(diff(F[K](t), t))+F[K](t) = K[M]*u

次の連立線形微分方程式を得ます。

> sys:={lindeq1, lindeq2, deq3};

sys := {(diff(x[K](t), `$`(t, 2)))+r*(diff(theta(t), `$`(t, 2)))+theta(t)*g = 0, m[K]*(diff(x[K](t), `$`(t, 2))) = F[K](t)-m[G]*(diff(x[K](t), `$`(t, 2)))-m[G]*r*(diff(theta(t), `$`(t, 2))), T[M]*(dif...sys := {(diff(x[K](t), `$`(t, 2)))+r*(diff(theta(t), `$`(t, 2)))+theta(t)*g = 0, m[K]*(diff(x[K](t), `$`(t, 2))) = F[K](t)-m[G]*(diff(x[K](t), `$`(t, 2)))-m[G]*r*(diff(theta(t), `$`(t, 2))), T[M]*(dif...

システムの過渡応答のグラフィカルな表示

パラメータに値を代入し、以下の初期条件の下で連立微分方程式を解きます。

> sysval := subs(val, u=10, sys);

sysval := {(diff(F[K](t), t))+F[K](t) = 1000, 1000*(diff(x[K](t), `$`(t, 2))) = F[K](t)-4000*(diff(x[K](t), `$`(t, 2)))-40000*(diff(theta(t), `$`(t, 2))), (diff(x[K](t), `$`(t, 2)))+10*(diff(theta(t),...sysval := {(diff(F[K](t), t))+F[K](t) = 1000, 1000*(diff(x[K](t), `$`(t, 2))) = F[K](t)-4000*(diff(x[K](t), `$`(t, 2)))-40000*(diff(theta(t), `$`(t, 2))), (diff(x[K](t), `$`(t, 2)))+10*(diff(theta(t),...

> init := {x[K](0) = 2, D(x[K])(0)=0, theta(0)=0, D(theta)(0)=0, F[K](0)=0};

init := {x[K](0) = 2, D(x[K])(0) = 0, theta(0) = 0, D(theta)(0) = 0, F[K](0) = 0}

> sol := dsolve(sysval union init, [x[K](t), theta(t), F[K](t)], type=numeric);

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

システムの位置、台車の速度とロープの角度の過渡応答 を時間に関してプロットします。

> plot(['op(2,sol(t)[2])'], t=0..30, axes=boxed,
title=" Position as a function of Time",
labels=["Time [s]", "x [K]"]);

[Plot]

> plot(['op(2,sol(t)[3])'], t=0..30, axes=boxed,
title=" Speed as a function of Time",
labels=["Time [s]", "Speed [K]"]);

[Plot]

> plot(['op(2,sol(t)[4])'], t=0..30, axes=boxed,
title=" Angle as a function of Time",
labels=["Time [s]", "Theta "]);

[Plot]

システム行列の設定

DEtools および linalgパッケージを使って、これまでの連立微分方程式を、一次の連立微分方程式に変更します。

> with(DEtools):
with(linalg):

Warning, the assigned name adjoint now has a global binding

Warning, the protected names norm and trace have been redefined and unprotected

convertsys コマンドを使って一時の連立方程式に変換します。

> syst := convertsys(sys, init, [x[K](t), theta(t), F[K](t)], t, X, X_p);

syst := [[X_p[1] = X[2], X_p[2] = (X[3]*g*m[G]+X[5])/m[K], X_p[3] = X[4], X_p[4] = -(X[3]*g*m[K]+X[3]*g*m[G]+X[5])/(r*m[K]), X_p[5] = (K[M]*u-X[5])/T[M]], [X[1] = x[K](t), X[2] = diff(x[K](t), t), X[3...syst := [[X_p[1] = X[2], X_p[2] = (X[3]*g*m[G]+X[5])/m[K], X_p[3] = X[4], X_p[4] = -(X[3]*g*m[K]+X[3]*g*m[G]+X[5])/(r*m[K]), X_p[5] = (K[M]*u-X[5])/T[M]], [X[1] = x[K](t), X[2] = diff(x[K](t), t), X[3...

genmatrix コマンドで連立方程式を行列の形式に変換します。 ここで diff(X(t), t) = A*x+b*u です。

> A:=genmatrix(map(rhs=0,syst[1]), [ X[1],X[2],X[3],X[4],X[5]], 'inhom');

A := matrix([[0, 1, 0, 0, 0], [0, 0, m[G]*g/m[K], 0, 1/m[K]], [0, 0, 0, 1, 0], [0, 0, -(g*m[K]+m[G]*g)/(r*m[K]), 0, -1/(r*m[K])], [0, 0, 0, 0, -1/T[M]]])

> b := map(x->-(1*x)/u, inhom);

b := vector([0, 0, 0, 0, K[M]/T[M]])

極配置による制御装置の設計

システムが可制御であるかの判定

最初に駆動装置の電圧により、システムが可制御であるか調べます。 制御行列を計算し、その行列式がゼロでないことを確かめます。

> Q:=concat(b, multiply(A,b), multiply(A^2,b), multiply(A^3,b), multiply(A^4,b));

Q := matrix([[0, 0, K[M]/(m[K]*T[M]), -K[M]/(m[K]*T[M]^2), (-m[G]*g/(m[K]^2*r)+1/(m[K]*T[M]^2))*K[M]/T[M]], [0, K[M]/(m[K]*T[M]), -K[M]/(m[K]*T[M]^2), (-m[G]*g/(m[K]^2*r)+1/(m[K]*T[M]^2))*K[M]/T[M], (...

> det(Q);

K[M]^5*g^2/(T[M]^5*m[K]^4*r^4)

もし K[M] がゼロでなく、T[M] , m[K] そして r が有限値であるならシステムは可制御です。

極を決定

制御回路を開発するために 制御システムの極を計算し表示します。

> pole:=eigenvalues(A);

pole := 0, 0, -1/T[M], (-r*m[K]*g*(m[K]+m[G]))^(1/2)/(r*m[K]), -(-r*m[K]*g*(m[K]+m[G]))^(1/2)/(r*m[K])

> pole:=subs(val,[pole]);

pole := [0, 0, -1, 2.214723459*I, -2.214723459*I]

> plot(map([Re, Im],pole), -2..0, -2.5..2.5, style=point, symbol=circle);

[Plot]

制御回路の目標は、制御エラーがゼロで、素早く、良く減衰する 安定な過渡効果を生成することです。 虚軸の極はできるだけ左にあったほうが良いので、結果として極は負の実部を持ちます。 干渉の感度をできるだけ低く押さえ、駆動がオーバーステアリング(過度に操作する)ことを防ぐためには システムの速度を十分に速くするべきです。すなわち、極を必要な限り左になるようにすれば良いのです。 制御回路の応答伝達関数が、素早く、良く減衰する2次遅れの装置となり、他の極の影響が無視できるように、 極を配置します。

2次遅れ装置の応答伝達関数 G  およびステップ関数応答 h  は次のようになります。

> G := s -> K/(1+2*d*T*s+T^2*s^2);

G := proc (s) options operator, arrow; K/(1+2*d*T*s+T^2*s^2) end proc

> h := t -> K-(K*exp(-(d*t)/T)*sin((sqrt(1-d^2)*t)/T))/sqrt(1-d^2);

h := proc (t) options operator, arrow; K-K*exp(-d*t/T)*sin(sqrt(1-d^2)*t/T)/sqrt(1-d^2) end proc

減衰振動の包絡は

> he := t -> K-(K*exp(-(d*t)/T))/sqrt(1-d^2);

he := proc (t) options operator, arrow; K-K*exp(-d*t/T)/sqrt(1-d^2) end proc

限界 he(infinity ) = K との差は

> abs(he-K) = (K*exp(-(d*t)/T))/sqrt(1-d^2);

abs(-he+K) = K*exp(-d*t/T)/(1-d^2)^(1/2)

減衰は d =0.7 が適当です。t =25s 後に、距離が 2%を過ぎないよう、その値を固定します。これにより T の値が、次のように計算されます。

> r1 := (K*2)/100 = (K*exp(-(d*t)/T))/sqrt(1-d^2);

r1 := 1/50*K = K*exp(-d*t/T)/(1-d^2)^(1/2)

> r2:=solve(subs(d=0.7,t=25,r1),{T});

r2 := {T = 4.118911534}

これは応答伝達関数を導き出します。

> G_ := subs(d=0.7,t=25,r2,G(s));

G_ := K/(1+5.766476148*s+16.96543222*s^2)

もう一つのの極を-1に設定します。これで特性多項式を得ます。

> p := (s+1)^3*simplify(denom(G_)/lcoeff(denom(G_)));

p := (s+1)^3*(0.5894338482e-1+.3398956226*s+s^2)

これで極の分布は次の様になります。

> plot(map([Re, Im],[fsolve(p=0, s, complex)]), -2..0, -1..1, style=point, symbol=circle);

[Plot]

制御回路のパラメータの計算

制御回路のパラメータを計算するために、 可制御行列の逆行列の最後の行が必要です。

> Q_:=inverse(Q);

Q_ := matrix([[0, (m[K]+m[G])/K[M], 0, m[G]*r/K[M], T[M]/K[M]], [(m[K]+m[G])/K[M], T[M]*(m[K]+m[G])/K[M], m[G]*r/K[M], T[M]*m[G]*r/K[M], 0], [T[M]*(m[K]+m[G])/K[M], r*m[K]/(K[M]*g), T[M]*m[G]*r/K[M], ...

> Q_5:=row(Q_,5);

Q_5 := vector([T[M]*m[K]*r/(K[M]*g), 0, T[M]*m[K]*r^2/(K[M]*g), 0, 0])

システム行列を特性多項式に入れ Q_5 をかけて、制御行列 R が得られます。

> R:=scalarmul(Q_5,coeff(p,s,0)):
for i from 1 to 5 do
 R:=matadd(R,multiply(Q_5,A^i),1,coeff(p,s,i)):
end do:
i:='i':
R:=subs(val,evalm(R));

R := vector([.6008499983, 5.267337178, 1424.575602, 135.3103468, 0.2339895620e-1])

フィルターの増幅

> TMP := array(1..5,1..5,[[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[b[5]*R[i]$i=1..5]]);

TMP := matrix([[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [.6008499983*K[M]/T[M], 5.267337178*K[M]/T[M], 1424.575602*K[M]/T[M], 135.3103468*K[M]/T[M], 0.2339895620e-1*K[M]/T[M...

> s[1] := 1/(multiply(multiply([1,0,0,0,0], inverse(matadd(TMP, A, 1, -1))), b));

s[1] := .6008499985

制御回路に対する方程式の計算

シミュレーションを行うために、制御回路に対する方程式を計算します。

> X := [x[i](t)$i=1..5];

X := [x[1](t), x[2](t), x[3](t), x[4](t), x[5](t)]

> Xp := [diff(x[i](t),t)$i=1..5];

Xp := [diff(x[1](t), t), diff(x[2](t), t), diff(x[3](t), t), diff(x[4](t), t), diff(x[5](t), t)]

> R := convert(evalm(R), vector);

R := vector([.6008499983, 5.267337178, 1424.575602, 135.3103468, 0.2339895620e-1])

> r2 := scalarmul(b, innerprod(R, X));

r2 := vector([0, 0, 0, 0, (.6008499983*x[1](t)+5.267337178*x[2](t)+1424.575602*x[3](t)+135.3103468*x[4](t)+0.2339895620e-1*x[5](t))*K[M]/T[M]])

> SYS := geneqns(A, X, matadd(matadd(Xp, scalarmul(b,s[1]*w), 1, -1
##                                Xp, scalarmul(b,s[1]*w), 1, -1
), r2, 1, 1));

SYS := {x[4](t) = diff(x[3](t), t), -(g*x[3](t)*m[K]+m[G]*g*x[3](t)+x[5](t))/(r*m[K]) = diff(x[4](t), t), -x[5](t)/T[M] = (diff(x[5](t), t))-.6008499985*w*K[M]/T[M]+(.6008499983*x[1](t)+5.267337178*x[...SYS := {x[4](t) = diff(x[3](t), t), -(g*x[3](t)*m[K]+m[G]*g*x[3](t)+x[5](t))/(r*m[K]) = diff(x[4](t), t), -x[5](t)/T[M] = (diff(x[5](t), t))-.6008499985*w*K[M]/T[M]+(.6008499983*x[1](t)+5.267337178*x[...SYS := {x[4](t) = diff(x[3](t), t), -(g*x[3](t)*m[K]+m[G]*g*x[3](t)+x[5](t))/(r*m[K]) = diff(x[4](t), t), -x[5](t)/T[M] = (diff(x[5](t), t))-.6008499985*w*K[M]/T[M]+(.6008499983*x[1](t)+5.267337178*x[...

> SYS:=simplify(subs(val, w=10, SYS));

SYS := {x[4](t) = diff(x[3](t), t), -4.905000000*x[3](t)-1/10000*x[5](t) = diff(x[4](t), t), -x[5](t) = (diff(x[5](t), t))-600.8499985+60.08499983*x[1](t)+526.7337178*x[2](t)+142457.5602*x[3](t)+13531...SYS := {x[4](t) = diff(x[3](t), t), -4.905000000*x[3](t)-1/10000*x[5](t) = diff(x[4](t), t), -x[5](t) = (diff(x[5](t), t))-600.8499985+60.08499983*x[1](t)+526.7337178*x[2](t)+142457.5602*x[3](t)+13531...SYS := {x[4](t) = diff(x[3](t), t), -4.905000000*x[3](t)-1/10000*x[5](t) = diff(x[4](t), t), -x[5](t) = (diff(x[5](t), t))-600.8499985+60.08499983*x[1](t)+526.7337178*x[2](t)+142457.5602*x[3](t)+13531...

> INIT := {x[1](0)=2, x[2](0)=0, x[3](0)=0, x[4](0)=0, x[5](0)=0};

INIT := {x[5](0) = 0, x[1](0) = 2, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0}

> sol := dsolve(SYS union INIT, [x[1](t),x[2](t),x[3](t),x[4](t),x[5](t)], type=numeric);

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

> plot(['op(2,sol(t)[2])'], t=0..30, axes=boxed,
title=" Position as a function of Time",
labels=["Time [s]", "Position [K]"]);

[Plot]

> plot(['op(2,sol(t)[3])'], t=0..30, axes=boxed,
title=" Speed as a function of Time",
labels=["Time [s]", "Speed [K]"]);

[Plot]

> plot(['op(2,sol(t)[4])'], t=0..30, axes=BOXED,
title=" Angle as a function of Time",
labels=["Time [s]", "Theta"]);

[Plot]

結論

Maple を使用することで、ローディング・ブリッジの微分方程式を設定しそれを線形化しました。そして、この微分方程式から、極配置によって制御装置を設計しました。この制御装置により、負荷は、事前に指定した時間後に、指定した誤差内で安定するようになります。

参考文献

Otto F?llinger, Regelungstechnik, H?thig.