レギュレータ

2020年4月29日

はじめに

バネ-ダンパー-質点系について、レギュレータによる制御を行う。

バネ-ダンパー-質点系のシミュレーション

2 自由度のバネ-ダンパー-質点系を線形システム状態方程式として解く。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)

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()
As, Bs, u = LM.linearize(q_ind=q, qd_ind=q.diff(t), A_and_B=True)

具体的な数値を設定する。

import numpy as np

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

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)

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

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

線形システム状態方程式をシミュレートする。ここでは SciPy を用いる。

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

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

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()

レギュレータ

線形システム状態方程式を次式で表す。

これに以下のようなフィードバック形の入力を与える。

状態方程式は次式のような形になる。

このようなシステムをレギュレータという。フィードバックゲイン \(F\) はシステムが安定になるように選ぶ。すなわち、\(A_0 + B_0 F\) の固有値の実数部分が負になるようにする。

\(b_0\)、\(k\) をベクトルとし、\(F = b_0 k^T\)、\(b = -B_0 b_0\) として、\(A_0 + B_0 F = A_0 - b k^T\) と表す。SciPy モジュールに \(A - B K\) の形の行列について、指定した固有値になる \(K\) を求める関数 signal.place_poles() があるので、これを用いる。これは、\(K\) によってシステムの極 (固有値) を配置するということで極配置と呼ばれる。ここでは 2 つめの質点にのみフィードバックを効かせるために \(b_0 = (0\ 1)^T\) とする。

b0 = np.array([[0.], [1.]])
b = -B0.dot(b0)
poles = [-2.+0.5j, -2.-0.5j, -2., -1.]

bunch = signal.place_poles(A0, b, poles)
k = bunch.gain_matrix
print(k)

A0f = A0 - b.dot(k)

from scipy import linalg
print(linalg.eigvals(A0f))

出力は以下のようになる。

[[ -4.5 -18.9  14.   -5.5]]
[-1.+0.j  -2.+0.5j -2.-0.5j -2.+0.j ]

指定した固有値になっているのがわかる。

SciPy の代わりに Python Control Systems Library を用いる場合は、place() 関数を用いる。

from control import matlab
k = matlab.place(A0, b, poles)

レギュレーションのシミュレーションを行うと、次のようになる。

sys = signal.lti(A0f, 0.*B0, C, D)
t, y, x = signal.lsim(sys, 0., t, x0)

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()

各変数が 0 に向かって制御されていることがわかる。

すべての質点にフィードバックを効かせられるとして、\(F\) を直接求めることもできる。

bunch = signal.place_poles(A0, -B0, poles)
F = bunch.gain_matrix
print(F)

A0f = A0 + B0.dot(F)

この場合のシミュレーション結果は以下のようになる。

参考

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