2. Appendix¶
2.1. 倒立振り子の運動方程式¶
倒立振り子装置は水平方向に一次元的に動く台車および一端がその台車に固定された倒立振り子から構成されている。台車と倒立振り子の質量をそれぞれ、 \(M\)および \(m\)で表す。また、倒立振り子の重心 \((x, y)\)の周りの慣性モーメントを \(I\)とする。台車の水平方向の位置および倒立振り子の鉛直方向からはかった角度をそれぞれ \(z\)および \(\theta\)とする。このとき、系のラグランジアンは次のように書き下せる。
此処で \(g\)は重力加速度である。倒立振り子の一端が台車に固定されていることは拘束条件、
で表せる。この拘束条件を解いて、ラグランジアンに代入することで、拘束条件のないラグランジアン、
を得る。
2.1.1. 線型化¶
\(\theta\)の変化が微小な範囲では、運動方程式を線型化近似で取り扱う。式 (51)のラグランジアンを微小な \(\theta\)について展開し、\(\theta\)の二次の範囲までをとると、
となる。これより運動方程式は、
となる。
2.2. 振り上げについての若干の考察。（未完)¶
フィードバックの失敗、大きな擾乱などにより、ある程度以上の振り子の傾きや台車位置がそのリミットに到達するような時には、 フィードバックの異常動作を避けるため、フィードバックの働きを一時的に停止し、台車位置をレール中心に移動させる。 プログラムはさらに台車を振り子の動きに合わせて動かすことで、倒れてしまった振り子を倒立の位置に持ってくる。 効率的にこの操作を行うための考察をここに述べる。 (未完)
2.3. Octave¶
Octaveは GPL(GNU General Public License)にしたがって配付されているフリーソフトウェアである。ほぼ MATLAB互換の処理系であって、制御に関するパッケージも MATLABの制御用パッケージのサブセット (＋独自拡張 )が実装されている。この報告では、 Octave2.1.34を使用した。 Octaveはもともと MATLAB互換なフリーソフトウェアを目指して開発が進められてきたが、一部では独自の発展を遂げている。
Octaveのソースコードなどは、 http://www.octave.org/download.html から入手可能である。
2.1.34は（2001年 5月 30日現在）開発版の状態であり、安定版は 2.0.16であるが、 2.1.34でも特に問題は起きていない。この報告で紹介した octaveプログラムは、すべて 2.1.34でのみ動作を確認している。 Octaveの controlパッケージは 2.0.16にくらべ、 2.1.34は大きく進歩している。 2.0.16で実行するためには、かなりの変更が必要と思われる。
改訂版の執筆時(2014年8月)の最新版は3.8.2となっている。再改訂版執筆時(2019年10月31日)の最新版は、 5.1.0(2019.2.25)となっている。
Linux/Windows/MacOSXにbinary installが提供されている。ソースコードも提供されているので、自ら buildすることも可能である。ビルドには Cコンパイラとfortran コンパイラ（たとえばgfortran)が必要 である。此処で利用しているcontrolパッケージはoctave本体のインストール後に、別途インストールする手続き が必要である。octaveにcontrolパッケージをインストールするには、octaveプロンプトで:
octave:1> pkg install -forge control
を実行すればよい。この他のオプションパッケージについては、http://octave.sourceforge.net/control/function/initial.html <Octave Forge> に詳しい。
octaveのControlシステムライブラリは2.xまでとなっている。新しいライブラリはSLICOTとよばれるライブラリを 利用する形で利用可能のようだが、調査中。SLICOTはPythonラッパーも用意されているらしい(slycot)。ただ、SLICOTのbuildにはFortran compilerが必要。
からgfortranやg77はダウンロード出来るみたいである。残念ながら、Python/slycont/SLICOTのc2d関数はMIMO系 (ここで考えて いるシステムは 入力２、出力１のMIMO系）はまだサポートされていない。(MIMOがサポートされないのは、現代制御理論としては 片手落ちと云う気がするけれど）ということで、現在Freesoftでこれらの係数を計算し直す方法がない。 (Octave2を動くようにする、Octave3でcontrolシステムが使えるようになるのをまつ等の手段を考えるのか？)
2.3.1. 極配置法¶
Octaveの標準制御パッケージには、極配置法によるフィードバックゲインを支援するための関数、 placeが用意されている。 MATLABの placeとは異なり、一入力系にしか適用出来ない (残念ながら）。
.. _CodeOctPoles:
sys=ss2sys(A,B,C,D);
poles=[-240, -5+2 j, -5 -2 j , -5];
K=place(sys,poles);
2.4. Python-control モジュール¶
プログラム言語PythonでもMATLAB/Octaveと同等の制御システム向けライブラリが 開発されている。python-control モジュール[A-1]がそれである。
python -m pip install control scipy slycot
2.4.1. Python control モジュールによる設計例¶
本文中で示された matlab/octaveのプログラムをpython.controlを用いて書き直してみた。 一部を取り出して紹介する。
A=np.matrix(
[[0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
[0.00000000e+00, 0.00000000e+00, -160.22529505e+02, 0.0],
[0.00000000e+00, 43.66133683, 573.55755483, 5.00908711]])
B=np.transpose(np.matrix([[0., 0. , 94.57663738, -337.84229409]]))
# 角度と位置がシステムの出力(観測値)
C=np.matrix([[1, 0, 0, 0],[ 0, 1, 0, 0 ]]);
D=np.matrix([[0], [0]]);
matrix_rank(A)
matrix_rank(ctrb(A,B))
sigw=diag([40,1,1,1])
sigv=eye(2)
# solve Ricatti equation.
#[H,p,e]=lqe(A,G,C,sigw,sigv);
[H,p,e]=lqr(np.transpose(A), np.transpose(C), sigw, sigv)
H=np.transpose(H)
Q=diag([10.0, 4.0, 1.0, 1.0])
R=eye(1)
[Kr,p,e]=lqr(A, B, Q, R)
eig(A-B*Kr)
python-controlのbode()関数はSISO系の場合だけに有効であるMIMOの場合には、
freqresp()メソッドあるいはfreqresp()関数を使って、
周波数応答を求めたあと、Bode図を描く。グラフの作成には、python.matplotlib.pyplotを利用している。
freqs=logspace(-2,1,1024)
mag,phase,omega=sys.freqresp(freqs)
pyplot.subplot(2,1,1)
pyplot.plot(omega, mag[0,0])
pyplot.subplot(2,1,2)
pyplot.plot(omega, phase[0,0])
振幅magおよび 位相phaseは、主力信号、入力信号、
および周波数に対応した指標を持つ三次元の配列として求められる。
上記の状態フィードバックゲインKrと状態推定きのゲインHを使って、コントローラーを定義し、そのコントローラによるフィードバックシステムを構成する。
図. 2.1 Feedback systemのダイアグラム:コントローラは対象システム(plant)の出力から、状態量を推定、推定状態量に基づき、制御量を決定する。
図. 2.2 Feedback systemのダイアグラム：コントローラにシステム入力を接続した場合。状態推定には、plantの出力(y)およびplantの入力(u)を用いる。plantへの参照値(r)が変化した場合にも、plantの状態と推定値の差が大きくならないと期待される。
from control import StateSpace,initial_responcse,feedback
plant=StateSpace(A,B,C,D)
controller=StateSpace(A-B*Kr-H*C, H, Kr, matrix(zeros((1,2))))
fbsys=series(plant, controller).feedback(1)
# or
fbsys=target.feedback(controller)
図. 2.3はこのフィードバックシステムによるインパルス応答を計算した例です。
図. 2.3 状態フィードバックによるインパルス応答。control.impulse_response()を用いた。一段目はコントローラーからの出力(モーターの速度指令値). 二段目と三段目はそれぞれ、台車の位置と振り子の角度を示す。青が対象系の状態値、橙色は状態推定器の推定値、緑は状態値と推定値の差を示している。
プログラムに実装する際には、連続系のパラメータを離散系に変換する。matlab/octaveでは\(c2d\)が使われたが、python.controlモジュールでは, control.sample_system()関数あるいは システムオブジェクトのsample()メソッドを使う。離散化の方法(method)として、"zoh"(default),"bilinear(tustin)", "forward_diff(euler)","backword_diff"が用意されている。応答関数で定義されるSISOシステムでは"matched"もサポートされる。
discrete_sys=sample_system(sys, scanperiod, "bilinear"
2.5. MATLAB から Python-control モジュールへの移行¶
python.controlモジュールには control.matlabモジュールが含まれており、 matlabの制御関係の関数などが含まれている。しかしながら、完全な互換性が あるわけでは無い。 東京電機大学出版会から出版されている「MATLABによる 制御理論の基礎」および「MATLABによる制御系設計」に掲載されているMATLAB プログラムは、「MATLABによる制御系設計」に付属のCD-ROMに収録されている。 これらのプログラムのいくつかをPython.control モジュールを使った Pythonモジュールに書き直してみた。この際に、注意すべき点をいくつか 上げておく。
2.5.1. importするモジュール¶
制御関係の関数などは、control.matlabモジュールをimportすることで、 概ね使えるが、グラフの描画などのためにmatplotlibから必要な関数をimportする必要がある、
from control.matlab import *
from numpy.linalg import matrix_rank as rank, eig, inv, pinv,svd
from numpy import matrix,zeros,diag,eye,sqrt,log10
from matplotlib.pyplot import clf,plot,show,draw,figure,subplot
from matplotlib.pyplot import xlabel,ylabel,axis,loglog,semilogx,semilogy,xlim,ylim,legend,title,grid,suptitle
これらの関数をインポートしておけば、大体の用は足りるでしょう。
numpy, matplotlibに代えて pylab() をimport *してしまうのも一つの考え方です[1]
| [1] | pylabのディクショナリには、900以上の名前が登録されており、これらを全てglobalスコープにimportしてしまうと、思わぬところで名前の衝突が起きるかもしれません。 |
2.5.2. 配列の定義¶
matlabでは行列の表記が文法に組み込まれており、簡単に行列を入力することが出来ます。例えば、
A=[ 1 2 3; 4 5 6]
は\(\begin{pmatrix} 1 & 2& 3 \\ 4 & 5 &6\end{pmatrix}\)を変数Aに割り当てます。
pythonではnumpyの matrix を使って
A=numpy.matrix(
(
(1,2,3),
(4,5,6)
))
と定義することで、変数Aに行列を割り当てます。実は numpy の行列は、
A=numpy.matrix(" 1 2 3; 4 5 6")
といった定義も可能です。しかし、
a=1;b=2
A=numpy.matrix(" a b 3; 4 5 6")
のように変数名をその値に含むような表式は受け付けてもらえません。結局のところ、
a=1;b=2
A=numpy.matrix([[a, b, 3],[4, 5, 6]])
と書き換えることが必要になります。ということで、実用的にはnumpy.matrix()を使うことになるでしょう。
2.5.3. 複素数の取り扱い¶
matlabでの複素数は '1+sqrt(2) i' のように記号 'i' を使って表現されます。pythonでは
複素数はcomplex(1,sqrt(2))()のようにcomplex()関数を使って複素数オブジェクトを
作成することになります。
% in matlab/octave
c= 1+ sqrt(2) i
# in python
c=complex(1, sqrt(2))
2.5.4. 対象システムの指定¶
matlabでは応答関数を表現するのに、s多項式の分子(num)と分母(den)のペア( [num, den] )を使うことが出来ます。
このペアを、bode,stepなどの関数に直接引数として渡すことが出来ます。これらの関数では、StateSpaceを表現する四つの行列のくみ、'A, B, C, D'を引数に指定すること
一方、python.controlではこれらの関数(bode, step ... )では、引数として渡せるシステムはStateStatかTransfer Functionのオブジェクトになり、このペア [num, den] (あるいは行列の組）を渡すことは出来ません。[num, den]のペアから応答関数のオブジェクトをtf()で作成して（あるいは ss関数を使って、行列の組からState Spaceオブジェクトを作って）、これを引数として渡します。
% in matlab/octave
bode([num, den])
step(A,B,C,D)
# in python
bode(tf(num,den))
step(ss(A,B,C,D))
プロパーな応答関数と状態空間表現との間は、tf2ss()あるいはss2tf()でお互いに
変換可能です。
tf2ss()とss2tf()は概ね交換可能で、tf2ss(ss2tf())() あるいは
ss2tf(tf2ss())()元のオブジェクトに戻ることと思いますが、例外的にこれが成り立たない場合があります。
具体的な例を示します。
>>> ss2=ss(-1,2,1,0)
>>> ss2tf(ss2)
ss2tf(ss2)
2
-----
s + 1
ですが、tf2ss(ss2tf(ss1))はValueErrorになってしまいます。これはpython controlでは応答関数tf([2],[1,1])を定義することは出来ますが、これをtf2ss()で変換することが出来ず、ValueErrorになることが原因です。 原理的には、StateSpaceで定義されたシステムを応答関数に変換しても、再び状態空間表示に
戻すことは、いつも可能なのですから、これはtf2ssの実装上の問題と思われます。(python controlではシステムは、
状態空間表示を基本とする方が良さそうです。)
応答関数を零点と極を指定することで指定することが出来ます。MATLABではこれにはzp2tf関数が使われます。
python.contorl.matlab では、scipy.signalモジュールからインポートされたzpk2tfあるいはzpk2ss関数を使用します。これらの関数はscipy.signalの関数であるため, [num,den]ペアあるいは 四つの行列の組みを値として返します。これらを
python controlでのシステムオブジェクト(LTIオブジェクト）とするにはtf()あるいはss()を使う必要があることに注意しましょう。
zpk表示との変換では、以下の例のように、tf2ssでおきた問題を生じないことがわかります。
ss(*zpk2ss(*ss2zpk(ss2.A,ss2.B,ss2.C,ss2.D)))
2.5.5. システムの離散化¶
連続系のシステムを離散化するプログラムc2dはcontrol.matlabにも用意されていますが、MATLABのそれ以外の離散化
関数は用意されていません。また、c2d()はMATLABのそれとは異なり、[num,den] ペアを受け付けることもありません。というわけで、c2dを使うより、StateSpaseやTransfer Functionのオブジェクトとsampleメソッドを使う方が簡単でしょう。また、matlabのc2dなど関数では、"zoh"(zeroth order hold)の他に、"foh"(first order hold), "bilinear"など様々な方法が用意されています。sampleではzohがデフォルトですが、その他には、"bilinear","euler","backward_diff"そして"gbt"が用意されています。これらの変換はzohを除き全てが、gbt(generalized bilinear transformation)のグループに属しており、gbtが要求するパラメータ \(\alpha\)が bilinear(0.5), euler(0), backward_dif(1.0)の場合になっています。前記の教科書を参考にした離散化アルゴリズムによるシステム応答の違いをpython.controlで計算した結果を以下に示します。
図. 2.4 システムの離散化法によるシステム応答(Bode図）の違い
2.5.6. システム応答を求める関数¶
control.matlabにはシステムの応答を求める関数として、MATLABと同じく, bode,impulse,step,initialの関数が用意されています。
また、これらの機能のうち、よく使われるbode図, 根軌跡図、閉ループのステップ応答をまとめて作成するsisotoolが用意されています。 control.sisotoolsはMATLABの同名のパッケージにちなんで導入されたとのことです。
# python
sys = control.tf([1000], [1,25,100,0])
control.sisotool(sys)
を実行することで、 次のような図が出力されます。
図. 2.5 control.sistoolによる出力, Bode図, root locus図、閉ループのステップ応答がまとめて出力されます。root locus図をマウスでクリックすることで、閉ループのゲインが変更されます。
bode, rlocus, sisotoolは内部でmatplotlibのplotコマンドを発行してグラフを描画します(つまりプログラム全体のどこかで、showコマンドが実行されることが必要だということです）。 その他のシステム応答関数 step/impulse/initial/freqrespは描画のためのデータを返しますが、描画自体は別途 明示的にプログラムする必要があります。 これらのシステム応答関数はMATLABのそれと同様のデータを返しますが、戻り値の順番、数が異なっている場合があります。その場合、python control/control.matlab側の関数に 呼び出し時のオプションを追加する必要があります。例えばstep関数では、状態変数 x の変化を出力する場合には、 'return_x'オプションの指定が必要です。
% MATLAB
y,t=step(sys)
y,t,x = step(sys)
# python
t,y = control.step_response(sys)
y,t = control.step(sys)
# pythonでは状態変数のシミュレーション結果は明示的に要求。
t,y,x=control.step_response(sys,return_x=True)
y,t,x=control.step(sys,return_x=True)