Pythonで考えるDCモータの制御（３）ランプ入力に対する１型ＰＩ制御の応答

はじめに

こんにちは、こうへいです。

前回まで、DCモータの速度PI制御について考えてみました。

rikei-tawamure.com

これからは、さらに進んでマイクロマウスやロボットでは速度を台形状に制御する事が多いので、その検証をできるように考えてみたいと思うのですが。まずは、基本のステップ入力とランプ入力について見ていきます。

制御におけるステップ入力とランプ入力

f:id:kouhei_ito:20200130193651j:plain — ステップ入力とランプ入力

階段状の信号を入力として捉えたものがステップ入力、直線の坂道状の信号を入力と捉えたものがランプ入力と呼ばれています。

ステップ入力について

制御において、この値を目標にしてそこまで速やかに達して欲しい等と考える場合はステップ入力と考えてよく、階段の高さが目標値になり基本的にはその値は変わりません。

単位ステップ入力のラプラス変換は次式で表されます。

$\displaystyle{ R(s)=\frac{1}{s}\\ }$

単位なので、１以外の大きさの場合は任意の数値をかけて表現します。

これまでの、記事では特に説明してなかったかもしれませんがステップ入力を入力として考えていたわけです。

ランプ入力について

一方でだんだん大きくなっていく（または、小さくなっていく）目標を設定したい場合は、線形に増大していく場合はランプ入力を使います。

$\displaystyle{ R(s)=\frac{1}{s^2}\\ }$

こちらも傾きが１の場合を表し任意の傾きのランプを表現したければ、任意の数値をかけて表現します。

前回のPI制御システム（１型）のランプ応答

モータのモデルは以下の１次遅れモデルで、機械的時定数を $\tau_m$ 、モータゲイン $K_m$ をとします。

$\displaystyle{ G_m(s)=\frac{K_m}{\tau_m s + 1}\\ }$

それぞれのパラメータについては、モータの各種値を用いて

${\displaystyle K_m = \frac{Ke}{DR+K^2} \\ }$

${\displaystyle \tau_m = \frac{JR}{DR+K^2}\\ }$

以上のようになります。

これらついては以前の記事を確認してみてください。

kouhei-no-homepage.hatenablog.com

比例ゲインを $K_p$ 、積分ゲインを $K_i$ とすると制御系の伝達関数は以下のようになります。

$\displaystyle{ W(s)=\frac{K_m K_p s + K_m Ki}{\tau_m s^2 +(1+K_m K_p )s + K_m Ki}\\ }$

このシステムの極を $\alpha$ 、 $\beta$ とすると傾き $a$ のランプ入力を入れた際の応答を表す式は、以下のようになります。

$\displaystyle{ \omega(t)=K_1e^{\alpha t}+K_2 e^{\beta t} + at - \frac{a}{K_m K_i}\\ }$

ここで $K_1$ 、 $K_2$ については、以下のようになります。

$\displaystyle{ K_1=\frac{a(K_m K_p \alpha + K_i)}{\tau_m (\alpha - \beta) \alpha^2}\\ }$

$\displaystyle{ K_2=\frac{a(K_m K_p \beta + K_i)}{\tau_m (\beta - \alpha) \beta^2}\\ }$

1型サーボと定常偏差

前回のPI制御器は積分器が一つしか入っていないので、これを１型サーボと呼びます。二つ入ると２型サーボとなります。

通常、ステップ入力は積分器一つと同じ形をしているので、１型にすると定常偏差がなくなりますが、ランプ入力では２型サーボにしなければ定常偏差が残ります。

したがって、今回のように１型サーボにランプ入力を加えると、上の応答の式からわかるのは

$\displaystyle{ \frac{a}{K_m K_i}\\ }$

の定常偏差が残るということです。

分母の積分ゲイン $K_i$ を大きくすると、定常偏差を小さくできる事がわかりますが。 $K_i$ を大きくすると、応答が振動的になります。

　ランプ入力に対する応答

それでは、目標の傾きを１として、いくつかゲインを決めて、実際に計算した結果を比較してみましょう。

PI制御システムの特性方程式（伝達関数の分母＝０）は次のようになります。

$\displaystyle{ \tau_m s^2 + (1+K_m K_p)s + K_m K_i=0\\ }$

この式の判別式から、極が異なる実数解、共益複素解、重解になるかならないかの境界値は重解になる場合の値です。判別式は次のようになります。

$\displaystyle{ D=(1+K_m K_p)^2 -4 \tau_m K_m K_i\\ }$

上式から $K_p$ の値から極が重解になる $K_i$ を求めるには次式を計算すると良いです。

$\displaystyle{ K_i=\frac{(1+K_m K_p)^2 }{4 \tau_m K_m}\\ }$

この時に $K_p=\frac{6}{500}$ のときには境界値の値が $K_i=1.6827$ 程度となり。この値より小さいと異なる二つの実数解、大きいと共役複素解になります。

異なる実数解の場合

kp=6/500, ki=1として計算するとランプ応答は次のようになります。

f:id:kouhei_ito:20200131155706p:plain — Ki=1の場合のランプ応答

誤差を計算してみると次のようになります。

f:id:kouhei_ito:20200131160034p:plain — Ki=1の場合のランプ応答の誤差

十分に時間が経つと、定常偏差が残る事がわかります。

共役複素解の場合

kp=6/500, ki=100として計算すると、極が共役複素解になり、ランプ応答は次のようになります。

f:id:kouhei_ito:20200131160328p:plain — Ki=100の場合のランプ応答

ki=1の場合とは時間軸の範囲が違うので注意です。

こちらも、誤差のグラフを示します。

f:id:kouhei_ito:20200131161446p:plain — Ki=100の場合のランプ応答の誤差

理屈どうりに、定常偏差は小さくなりますが、応答は振動的になります。

ランプ応答計算プログラム

import numpy as np
import matplotlib.pyplot as plt

#1724DCモータの諸元
R=3.41
K=6.59e-3
L=75e-6
D=1.4e-7
J=1e-7

#1次遅れモデルのパラメータ
K_m=K/(D*R+K**2)
tau_m=J*R/(D*R+K**2)


# 比例ゲインと積分ゲイン
kp=6/500
ki=100

# 傾きの目標値
ref=1

# 特性方程式の係数
A=tau_m
B=1+K_m*kp
C=K_m*ki

# 極の計算
p=np.roots([A, B, C])
print('Pole={}'.format(p))

# 応答式の係数
k1=ref*K_m*(kp*p[0]+ki)/tau_m/(p[0]-p[1])/p[0]/p[0]
k2=ref*K_m*(kp*p[1]+ki)/tau_m/(p[1]-p[0])/p[1]/p[1]
k32=ref
k31=-ref/K_m/ki
print(k1,k2,k32,k31)

# 応答の計算
t=np.linspace(0, 0.1, 15000)
y=k1*np.exp(p[0]*t) + k2*np.exp(p[1]*t) + k32*t + k31
#y=k32*t - (k1*np.exp(p[0]*t) + k2*np.exp(p[1]*t) + k32*t + k31) #誤差を計算する場合はこちら

# 可視化
plt.figure(figsize=(8,4))
plt.plot(t, y, linewidth=1, label='1stOrderModel')
plt.plot(t, k32*t, label='Reference') #誤差を表示するときはこちらをコメントにする
plt.legend(loc='lower right')
plt.xlabel('Time(s)')
plt.ylabel('Anglurr Velocity(rad/s)')
plt.grid()
plt.show()