運動方程式の導出とシミュレーション

2020年4月21日

はじめに

運動方程式を導出し、シミュレーションを行う。

運動方程式の導出

バネ-ダンパー-質点が 2 つ直列につながった問題を考える。運動方程式は連立方程式の形になる。直接求めてもよいが、系の自由度が増えてくると大変である。ここではラグランジュの方法を用いる。ラグランジュの運動方程式は次式で表される。

$$\frac{d}{d t} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = Q_i - \frac{\partial F}{\partial \dot{q}_i}$$

ここで \(L = T - V\) はラグランジアンで、\(T\) は運動エネルギー、\(V\) はポテンシャルエネルギーである。\(q_i\) は一般化座標で、\(\dot{q}_i\) は一般化速度、\(Q_i\) は一般化力である。\(F\) は次式で表される粘性減衰力を与える散逸関数である。

$$f_{di} = - \frac{\partial F}{\partial \dot{q}_i}$$

質点の位置を \(z_1\)、\(z_2\)、質点にかかる力を \(f_1\)、\(f_2\)、質量を \(m_1\)、\(m_2\)、バネ定数を \(k_1\)、\(k_2\)、ダンパーの減衰係数を \(c_1\)、\(c_2\) とする。1 側のバネ-ダンパーの端が固定されているものとする。運動エネルギーは次式で表される。

$$T = \frac{1}{2} m_1 \dot{z}_1^2 + \frac{1}{2} m_2 \dot{z}_2^2$$

ポテンシャルエネルギーは次式で表される。

$$V = \frac{1}{2} k_1 z_1^2 + \frac{1}{2} k_2 (z_2 - z_1)^2$$

散逸関数は次式で表される。

$$F = \frac{1}{2} c_1 \dot{z}_1 + \frac{1}{2} c_2 (\dot{z}_2 - \dot{z}_1)^2$$

これらをラグランジュの運動方程式に入れればよいのだが、ここでは Python モジュールの SymPy を用いる。

SymPy による運動方程式の導出

SymPy は数式処理を行うためのモジュールであるが、運動方程式を導出する機能があるので、それを用いる。

モジュールを準備。

from sympy import symbols, diff, Matrix, lambdify
from sympy.physics.mechanics import dynamicsymbols
from sympy.physics.mechanics import ReferenceFrame, Point, Particle
from sympy.physics.mechanics import LagrangesMethod, Lagrangian

変数や定数を準備。

t = symbols("t")
z1, z2, f1, f2 = dynamicsymbols("z1 z2 f1 f2")
m1, m2, c1, c2, k1, k2 = symbols("m1 m2 c1 c2 k1 k2")
q = Matrix([z1, z2])

質点の設定。

N = ReferenceFrame("N")

p1 = Point("p1")
v1 = z1.diff(t)
p1.set_vel(N, v1*N.x)
pa1 = Particle("pa1", p1, m1)
pa1.potential_energy = k1*z1**2/2 + k2*(z2 - z1)**2/2

p2 = Point("p2")
v2 = z2.diff(t)
p2.set_vel(N, v2*N.x)
pa2 = Particle("pa2", p2, m2)

ここでは、質点 1 のポテンシャルエネルギーに系全体のポテンシャルエネルギーを設定している。運動エネルギーは勝手に計算される。力の設定は以下のようになる。

F = c1*v1**2/2 + c2*(v2 - v1)**2/2
fc1 = -F.diff(v1)
fc2 = -F.diff(v2)
fl = [(p1, (f1 + fc1)*N.x), (p2, (f2 + fc2)*N.x)]

ラグランジアンを作り、ラグランジュの運動方程式を導く。

L = Lagrangian(N, pa1, pa2)
LM = LagrangesMethod(L, q, forcelist=fl, frame=N)
LM.form_lagranges_equations()

運動方程式を次式の形で表す。

$$\dot{x} = A_0 x(t) + B_0 u(t)$$ $$x = (x_1\ x_2\ x_3\ x_4) = (z_1\ z_2\ \dot{z}_1\ \dot{z}_2)$$

これは線形時不変システム状態方程式である。つまり運動方程式を状態空間表示する。これを線形化といっている。

As, Bs, u = LM.linearize(q_ind=q, qd_ind=q.diff(t), A_and_B=True)

\(A_s\)、\(B_s\)、\(u\) はそれぞれ以下のようになる。

$$A_s = \displaystyle \left[\begin{matrix}0 & 0 & 1 & 0\\0 & 0 & 0 & 1\\\frac{- k_{1} - k_{2}}{m_{1}} & \frac{k_{2}}{m_{1}} & \frac{- c_{1} - c_{2}}{m_{1}} & \frac{c_{2}}{m_{1}}\\\frac{k_{2}}{m_{2}} & - \frac{k_{2}}{m_{2}} & \frac{c_{2}}{m_{2}} & - \frac{c_{2}}{m_{2}}\end{matrix}\right]$$ $$B_s = \displaystyle \left[\begin{matrix}0 & 0\\0 & 0\\\frac{1}{m_{1}} & 0\\0 & \frac{1}{m_{2}}\end{matrix}\right]$$ $$u = \displaystyle \left[\begin{matrix}\operatorname{f_{1}}{\left(t \right)}\\\operatorname{f_{2}}{\left(t \right)}\end{matrix}\right]$$

シミュレーション

導出した運動方程式を用いてシミュレーションを行う。導出した式に具体的な数値を適用する。

m1n = 1.
m2n = 1.
c1n = 0.5
c2n = 0.5
k1n = 0.2
k2n = 0.2

これを次のように \(A_s\)、\(B_s\) に適用する。

As_func = lambdify((m1, m2, c1, c2, k1, k2), As, modules="numpy")
A0 = As_func(m1n, m2n, c1n, c2n, k1n, k2n)

Bs_func = lambdify((m1, m2), Bs, modules="numpy")
B0 = Bs_func(m1n, m2n)

SciPy モジュールを用いて、運動方程式を線形時不変システム方程式としてシミュレーションする。

import numpy as np
from scipy import signal

C = [0., 0., 0., 0.]
D = [0., 0.]
x0 = [0.5, 1., 0., 0.]

dt = 0.1
end_time = 60.
t = np.linspace(0., end_time, int(end_time/dt) + 1)

sys = signal.StateSpace(A0, B0, C, D)
t, y, x = signal.lsim(sys, 0., t, x0)

グラフを次のように描く。

from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

plt.plot(t, x[:, 0], label="$x_1$")
plt.plot(t, x[:, 1], label="$x_2$")
plt.plot(t, x[:, 2], label="$x_3$")
plt.plot(t, x[:, 3], label="$x_4$")
plt.xlim(0., end_time)
plt.xlabel("$t$")
plt.ylabel("$x$")
plt.legend()
plt.show()

SciPy の代わりに Python Control Systems Library を用いるなら、次のようになる。

from control import matlab

sys = matlab.StateSpace(A0, B0, C, D)
y, t, x = matlab.lsim(sys, 0., t, x0)

lsim() の返り値の順番が SciPy の場合と異なるのに注意。

自前で計算するには、次のようにすればよい (ただし \(B\) や \(C\) の項は省略)。

from scipy.linalg import expm

n = len(x0)
A = expm(A0*dt)

steps = int(end_time/dt) + 1
x = np.zeros((steps, n))

x[0, :] = x0

x1 = np.array(x0)
for step in range(1, steps):
    x1[:] = A.dot(x1)
    x[step, :] = x1

参考

  • 大住晃: 構造物のシステム制御, 森北出版 (2013).