Itamika logo Itamika
分解加速度制御

分解加速度制御

ロボットアームの制御について

 ロボットアームの制御をする場合、足回りと同じように、目標の座標から目標の角度を逆運動学で解析し、現在の角度との誤差に対してPID制御などのフィードバック制御を適用する しかし、これを行うには以下のような問題がある

  • 三自由度以下のロボットアームでなければならない 1
  • ロボットアームは基本、非線形なため、静止動作含め制御するのが難しい

分解加速度制御は、この2つ目の非線形なモデルを制御する方法であるといえる 簡潔に説明すると、ロボットアームの運動方程式に基づいてトルクを計算することで、非線形状態フィードバック制御を行うような制御方法である よって、運動方程式の計算にラグランジュの運動方程式を使用するため、これから説明する

※7月19日追記。最近、twitterで聞かれたので、これを機会に書き直してみることにする。 ※2022年12月8日変更。ちゃんと書き直した。


ラグランジュの運動方程式

ラグランジュの運動方程式とか一般人は知らないので説明する まず、以下の式が公式である

ddt(Lx˙)Lx=0\frac{d}{dt}\left(\frac{\partial \boldsymbol{L}}{\partial \dot{x}}\right)-\frac{\partial \boldsymbol{L}}{\partial x}=0
  • tは時間
  • xは位置
  • Lはラグランジアン

を意味している

ラグランジアンは以下を示すものである L = T ー U Tは運動エネルギー Uはポテンシャルエネルギー を意味している

ポテンシャルエネルギーは、位置エネルギーを含むものなので、とりあえず位置エネルギーだと思っておけばよい

導出

ニュートンの運動方程式

ここから、ラグランジュの運動方程式がニュートンの運動方程式と同じになることを説明する 大学物理ではニュートンの運動方程式は以下のように示す

mx¨=Fm\ddot{x}=\boldsymbol{F}
  • mは質量
  • Fは力

を意味している ここでFはポテンシャルエネルギーから受ける力であるため、以下のようにあらわせる

F=Ux\boldsymbol{F}=-\frac{\partial \boldsymbol{U}}{\partial x}

ラグランジュの運動方程式

大学物理では運動エネルギーは以下のように示す

T=12mx˙2\boldsymbol{T}=\frac{1}{2}m\dot{x}^2

よって、運動エネルギーはXの一階微分(速度)の関数といえるため、ラグランジアンは

L=T(x˙)U(x)\boldsymbol{L}=\boldsymbol{T}(\dot{x})-U(x)

と表現できる。 これを偏微分すると以下のようになる

Lx˙=Tx˙=mx˙,  Lx=Ux\frac{\partial \boldsymbol{L}}{\partial \dot{x}}=\frac{\partial \boldsymbol{T}}{\partial \dot{x}}=m\dot{x},\;\frac{\partial \boldsymbol{L}}{\partial x}=-\frac{\partial U}{\partial x}

右の式はニュートンの運動方程式の右辺と同じである。 また、左の式を t で微分すると以下のようになる。

ddt(Lx˙)=mx¨\frac{d}{dt}\left(\frac{\partial \boldsymbol{L}}{\partial \dot{x}}\right)=m\ddot{x}

これは、ニュートンの運動方程式の左辺と同じである ここで、ラグランジュの運動方程式は、以下であったため、

ddt(Lx˙)Lx=0\frac{d}{dt}\left(\frac{\partial \boldsymbol{L}}{\partial \dot{x}}\right)-\frac{\partial \boldsymbol{L}}{\partial x}=0

移項すれば、ラグランジュの運動方程式とニュートンの運動方程式は同じであることがわかる。

分解加速度制御 ロボットアームの運動方程式と制御

運動方程式の導出

ロボットアームの運動方程式は、ラグランジュの運動方程式を使用して計算する

※今回は簡単のために、1関節のロボットアームを考えることにする。 ※上の説明で x (位置) を示していたものは、ロボットアームで扱うのは角度であるため、θで表現する

ラグランジュの運動方程式は、以下の式になる

ddt(Lθ˙)Lθ=τ\frac{d}{dt}\left(\frac{\partial \boldsymbol{L}}{\partial \dot{\theta}}\right)-\frac{\partial \boldsymbol{L}}{\partial \theta}=\tau
  • θはアーム関節の角度
  • τ(タウ)は加えられる力 (目標トルクとなる)

を意味している

また、ラグランジアンは以下のようになる (関節数が増えると複雑になる)

L=12ml2θ˙2+12Iθ˙2mglcosθ\boldsymbol{L}=\frac{1}{2}ml^2\dot{\theta}^2+\frac{1}{2}I\dot{\theta}^2-mglcos\theta
  • Iは慣性モーメント
  • l(エル)はアームの長さ

を意味している

偏微分すると、以下のようになる

Lθ˙=ml2θ˙+Iθ˙,  Lθ=mglsinθ\frac{\partial \boldsymbol{L}}{\partial \dot{\theta}}=ml^2\dot{\theta}+I\dot{\theta},\;\frac{\partial \boldsymbol{L}}{\partial \theta}=mglsin\theta

また、左の式は t で微分しないといけないので、

ddt(Lθ˙)=ml2θ¨+Iθ¨\frac{d}{dt}\left(\frac{\partial \boldsymbol{L}}{\partial \dot{\theta}}\right)=ml^2\ddot{\theta}+I\ddot{\theta}

よって、このロボットアームの運動方程式は、

ml2θ¨+Iθ¨mglsinθ=τml^2\ddot{\theta}+I\ddot{\theta}-mglsin\theta=\tau

となる

よく見ると、左辺は慣性力と重力を表していることに気が付くだろう これは、一般的に以下のような式で表現される

J(θ)θ¨+G(θ)=τJ(\theta)\ddot{\theta} + G(\theta) = \tau
  • Jは慣性モーメントとかの行列
  • Gは重力

を意味している

2関節以上になると遠心力・コリオリ力が出てくる さらに、摩擦力を考慮してラグランジアンを立てれば同じように計算できる

角加速度の計算

θの二階微分以外は、慣性モーメントと重力なので、アーム角度から容易に計算することができる つまり、θの二階微分つまり角加速度を計算すれば目標トルクを求めることが可能となる よって、以下のような非線形状態フィードバック補償を考える

J(θ)τ+G(θ)=τJ(\theta)\tau^{\prime}+ G(\theta) = \tau

(θの二階微分がτプライムに置き換わっている)

これにより線形な系まで落とし込むことができた しかし、運動方程式とのモデル誤差は避けられないため、これを補償するために、サーボ補償器というものを定義する

サーボ補償器は以下のようなものである

τ=θd¨+Kd(θd˙θ˙)+Kp(θdθ)\tau^{\prime} = \ddot{\theta_{d}} + K_{d}(\dot{\theta_{d} }- \dot{\theta}) + K_{p}(\theta_{d} - \theta)
  • θdは目標角度
  • KdはDゲイン
  • KpはPゲイン

を意味している

θdーθ=eと置くと、ゲインを適切に求めることで、eを0に収束させることができることがわかる つまりPD制御ってこと

座標系の変換(上平のオリジナル)

この状態で制御をするには、すべてθベース、つまりアームの関節角度ベースで制御しなければならない 低能な僕では、ゲインをどのくらいにすればよいのか見当もつかないし、 経路生成的問題でも、普通直交座標系の経路を生成するはずなので、 直交座標系でPD制御を行いたい

よって、以下のようにサーボ補償器を定義し直す

τ=x¨d+Kd(x˙dx˙)+Kp(xdx)\tau^{\prime\prime}=\ddot{x}_{d}+K_{d}(\dot{x}_{d}-\dot{x})+K_{p}(x_{d}-x)

しかし、これでは座標系が違うので、これを座標変換することを考える 3ステップあるため、順に説明する まず、角速度を速度に変換する式は、

x˙(t)=J(t)θ˙(t)\dot{x}(t) = J(t)\dot{\theta}(t)
  • J(t)はヤコビ行列

を意味している

つぎに、両辺をtで微分すると

x¨(t)=J˙θ˙(t)+Jθ¨(t)\ddot{x}(t) = \dot{J}\dot{\theta}(t) + J\ddot{\theta}(t)

これを移項とかすると

τ=θ¨(t)=J1(x¨(t)J˙θ˙(t))\tau^{\prime} = \ddot{\theta}(t) = J^{-1}(\ddot{x}(t) - \dot{J}\dot{\theta}(t))

になる よって、上のサーボ補償器で計算したxの二階微分をこれに代入して、目標角加速度を計算する

これで一通り終わりだが、実際にはモデル化誤差が問題になってくる。そのため、式に積分項を追加することで誤差を吸収する方法が考えられる。よって、フィートバックパラメータKiK_{i}を用いて再定義する。

τ=xd¨+Kd(xd˙x˙)+Kp(xdx)+Ki(xdx)dt\tau^{\prime\prime} = \ddot{x_{d}} + K_{d}(\dot{x_{d}} - \dot{x}) + K_{p}(x_{d} - x) +K_{i}\int(x_{d} - x)dt

僕はこれでいけたと思っていますが、本当はあっているかわからないので、もし間違っているところに気付いたら教えてください。お願いします。

この計算方法は、目標トルクを計算するため、計算トルク法と呼ばれる。また、加速度を用いて計算するため分解加速度制御とも呼ばれる。 ちなみに、運動方程式を G(θ)=τG(\theta) = \tau の形にすると重力の分しか出力をしていないので重力のみを補償する形になり、静止することができる。

モータを制御する

これで、アームの目標トルクは計算することができた。 このトルクを制御する方法として、以下の2つが考えられえる

  • 電流制御によるトルク制御を行い、モーターを制御する
  • PWMによって制御するために、モーターをモデル化し、トルクをPWM値へと変換する

電流制御器は簡単に設計できるため、目標トルクを電流に変換して、適当にPIDとかで電流制御やってればいける モーターのモデル化は次に示す

モーターのモデル化

 まず、モーターのデータシートより、モーターのトルク定数、端子間抵抗、誘導起電力定数を読みるか、実測する。そして、計算トルク法により計算した目標トルクをトルク定数で割り、目標の電流を計算する。これと端子間抵抗をかけたものを目標電圧1とする。また、誘導起電力定数とモーターの現在の角速度をかけたものを目標電圧2とする。目標電圧1と目標電圧2を足して、電源の電圧で割ると、目標のPWM値が計算できる。これをモータードライバの出力にすることで、ロボットアームを制御することができる。ちなみに、トルク定数をKτK_{\tau}、端子間抵抗をRR、誘導起電力定数をKEK_{E}、電源電圧をV0V_{0}、角速度をω\omega、与えるトルクをτref\tau^{ref}とすると、

PWM=RτrefKτ+KEωV0\Large PWM = \frac{R\frac{\tau^{ref}}{K_{\tau}} + K_{E}\omega}{V_{0}}

こんな感じになる。

 この制御方法はロボットアームとモーターの両方のモデル化が必要になるため、モデルをきちんと計算することが重要になってくる。しかし、(4)の式はPIDになっているため、このパラメータしだいで意外と綺麗に動かすことができるため、おすすめの制御方法である。


軌道生成

 軌道生成には三次以上のスプライン曲線を用いることで加速度まで連続性が保たれた、滑らかな軌道をつくる。加速度でも滑らかな曲線をつくるために五次スプライン曲線を使うことも考えられるが、今回は三次スプライン曲線を用いた。  そもそも、スプライン曲線とはある多項式によって描かれる、複数の点を通るなめらかな曲線のことで、iiを点の順番とすると、三次のスプライン曲線の公式は、

Si(x)=ai+bi(xxi)+ci(xxi)2+di(xxi)3S_{i}(x) = a_{i} + b_{i}(x - x_{i}) + c_{i}(x - x_{i})^2 + d_{i}(x - x_{i})^3

このa,b,c,da,b,c,dを求めるだけで軌道が生成できるため、三次スプライン曲線はよく使われる。a,b,c,da,b,c,dは、三次スプライン曲線の条件式5つを連立させて解く。しかし、この公式のxxはロボコンではおそらく時間を示すttを用いて使う(自分はそうした)ため、これからはttで置き換えることにし、xxは座標のx,y,zx,y,zとかのxxとして扱うことにすると、公式はこう置き換えることができる…と思う。

Si(t)=ai+bi(tti)+ci(tti)2+di(tti)3S_{i}(t) = a_{i} + b_{i}(t - t_{i}) + c_{i}(t - t_{i})^2 + d_{i}(t - t_{i})^3

また、点の数をnnとすると、条件は以下になる。

条件1

Si(ti)=xiS_{i}(t_{i}) = x_{i}

条件2

Si(ti+1)=Si+1(ti)=xi+1S_{i}(t_{i+1}) = S_{i+1}(t_{i}) = x_{i+1}

条件3

Si(ti+1)=Si+1(ti+1)S^{\prime}_{i}(t_{i+1}) = S^{\prime}_{i+1}(t_{i+1})

条件4

Si(ti+1)=Si+1(ti+1)S^{\prime\prime}_{i}(t_{i+1}) = S^{\prime\prime}_{i+1}(t_{i+1})

条件5

S0(0)=Sn1(xn)=0S^{\prime\prime}_{0}(0) = S^{\prime\prime}_{n-1}(x_{n}) = 0

a,b,c,da,b,c,dを点と点の間の数だけ解くことでスプライン曲線を生成することができる。また、上の公式を使うと目標の座標(位置)が求められるが、あとは、公式を時間ttで微分すれば速度の公式になって、もう一回微分すれば加速度の式になるだけだが、一応実装用に書いておく。ttは点と点の間を一秒としたときの時間である。

x(t)=a+bt+ct2+dt3x(t)=b+2ct+3dt2x(t)=2c+6dtx(t) = a + b*t + c*t^2 + d*t^3 \\ x^{\prime}(t) = b + 2c*t + 3d*t^2 \\ x^{\prime\prime}(t) = 2c + 6d*t

(y,zも同様)(y, zも同様)

ちなみに、ロボコンでバスタオルと取った時の軌道は コメント 2019-10-29 184855.png (19.5 kB) こんな感じでした


注釈

Footnotes

  1. 三自由度までしか逆運動学が求められず、解析的に求めることができないため、ニュートン・ラフソン法などで近似するしかない