Pythonで考えるDCモータの制御（５）制御の基礎の基礎とPython

はじめに
ラプラス変換と逆変換
道具としてのラプラス変換：微分方程式を解く道具
覚えておきたいラプラス変換
- 覚えておきたいラプラス変換ベスト２
- 覚えておきたいラプラス変換　よく入力信号として登場するもの
とりあえずf(t)のラプラス変換はF(s)とする
- 小文字は大文字へ
覚えておきたいラプラス変換の性質
有利関数
有利関数の極と零点
有利関数の逆ラプラス変換と部分分数展開
- 重複する極がない場合
- 重複する極がある場合
部分分数展開における係数の求め方
- 極に重複がない場合
- 極に重複がある場合
伝達関数
DCモータの定常角速度の計算
過渡応答（時間応答）の計算
ブロック線図と伝達関数
フィードバック制御
Pythonでの多項式関数の取り扱い
おわりに
参考文献

はじめに

これまでの記事を読み進めるための制御工学に関する基礎的なことをまとめておこうと思い、執筆しました。

制御工学に関する計算とPython

「Pythonで考えるDCモータの制御」シリーズなのにあまりPythonに触れていませんでした。 PythonにはMatlabライクに使えるパッケージが存在しています。

僕はあまり知りませんが、他にも制御に関するモジュールはあると思います。上記のPython Control Systems Libraryはかなり充実していると思います。

ということで本シリーズもPython-control等を使えば簡単に結果が出せると思うのですが仕事ではなく興味関心に基づいてやっていることなので、効率を追求する必要もなくちょっと計算がしんどいところをコンピュータ様に最低限のお手伝いをしてもらう考え方で当面は記事を書いていこうと思います。 MatlabやScilabなどがやっていることを類推するようなもので、それも楽しいです。

更に進んで、これまでやっていたことの土台の上に何か積み重ねていこうとしたときは大いに、素晴らしいツールを活用していきたいです。その側面でも、いつかは記事を書いてみたいです。

余談

はじめてMatlabを触った２５年ぐらい昔、すごく感動したものです。

ラプラス変換もろくにできなかったですが Simulinkで伝達関数のブロックを挿入してブロック線図を作ったら、簡単にステップ応答が現れる！興奮しました。数式追いかけて現象理解していなくても、見た目で出てくるというのは良いものです。

制御工学にものすごく入っていきやすくなったような気がしています。まあ、同じ環境があっても入ってこない人も大勢いますが・・・

その後、研究室に転がっていたAD/DA変換ボードがSimulink上で動かせて倒立振子をそれで立てたときには鳥肌立ったのを覚えています。（RealTimeWorkShopの初期のころのバージョンだったと記憶しています。）出来たときに先輩や先生や同期に自慢しました。

初めて勉強するときにMatlabのようなツールがあってラッキーでした。（先生や先輩に感謝！）

Matlabはマウサーは自由に使えます。良い時代です。

jp.mathworks.com

今回は制御についてです

今までの記事で説明不足感が有ったところを補うために制御について、書いてみようかなと思います。

ラプラス変換と逆変換

ラプラス変換を深く語るほどの力量は無いので、制御の道具としての使い方の一端に触れさせてもらおうと思います。

ラプラス変換

関数 $f(t)$ に対して次の定積分操作をすることで $F(s)$ にする操作をラプラス変換と呼びます。

$\displaystyle{ F(s)=\int_0^\infty f(t) e^{-st} dt\\ }$

上の定積分は時間が０から始まっています。ラプラス変換の世界にはマイナスの時間はありません。

以下のようにも記述します。

$\displaystyle{ \mathcal{L}\left[f(t)\right]=F(s) }$

$f(t)$ は原関数とも呼ばれます。

ラプラス逆変換

逆に $F(s)$ から元の $f(t)$ を求める操作をラプラス逆変換と言います。

ラプラス変換はしんどい？！

ラプラス変換の定義にしたがって積分計算をするには部分積分というテクニックをつかって「しこしこ」することになるのですが、高校や大学以来、積分なんてした事が無い場合、「あんまりしたく無い気分」が盛り上がりますよね。

そうして、制御では道具として、ハンダゴテやスパナやドライバのように使うこととなるラプラス変換が嫌になります。

ハンダ付もハンダゴテが嫌いだから、嫌いというのは稀で、使いこなせないか、下手くそだから嫌になるんですよね。

道具が使いこなせないから、制御が嫌いになるのは残念でならないんですが、その先を見据えて行くしか無いのかなあ。

本当は、「ここはいっちょ、やっていやるか！」と１０００回ぐらい積分計算の練習をすれば、何のことは無いと実感し「ストン」と腑に落ちます。（自分は教える事があるので、学生さんたちに出す課題を作って、それを解く作業をするので１０００回に近い回数は練習していると思います。ようやく慣れてきました（笑））手を動かすのが唯一の正解だと実感しています。（それが、中々できないからねえ・・・）

いや！今はMatlabやPythonがありますから、まずはそっちから入るのも手ですよ。（たぶん・・・）「教え方もっと工夫しなきゃなあ。」と反省することしきりです。

道具としてのラプラス変換：微分方程式を解く道具

微分方程式を代数方程式へ

結局南極、ラプラス変換は何に使うのかという事ですが、線形の微分方程式をうまいこと解く事ができます。ラプラス変換を微分方程式に施すと、微分方程式は代数方程式の形になります。代数演算の操作で以後制御対象のモデルが扱う事ができます。その世界では微分や積分の演算をすることはなくシンプルになります。

シンプルな操作の代償としてラプラス変換の手間をかける

シンプルに作業を進めるために、ラプラス変換（積分の計算）で一手間かけると言ったイメージです。しかし、あとで触れますがその積分計算自体をほとんどする必要がなく、更に手軽になります。

入力も掛け算するだけで表現可能

制御対象を微分方程式で表してラプラス変換すると制御対象の入出力関係をあらわす 伝達関数がえられます。制御対象に何か入力を加えることを考えた時、ラプラス変換された世界ならば、伝達関数に入力関数のラプラス変換表現を 掛けるだけで出力を表す式が得られます。

他方、制御対象を微分方程式で表したまま、これをやろうとすると、入力の項に、入力を表す関数を代入し式変形し、微分方程式を解かなければなりません。制御屋はここで、結局のところ慣れ親しんだラプラス変換を使いますから、初めからラプラス変換を使っていた方が良いということになります。

覚えておきたいラプラス変換

制御工学で活用する際はラプラス変換の積分計算を実際にできる必要性はあまりありません。そのかわり、いくつかの原関数のラプラス変換の結果や、ラプラス変換の性質を覚えておいた方が良いです。

覚えておきたいラプラス変換ベスト２

覚えとおきたいランキングNo .1

$\displaystyle{ \mathcal{L}\left[ e^{-a t} \right] = \frac{1}{s+a}\\ }$

覚えとおきたいランキングNo .2

$\displaystyle{ \mathcal{L}\left[ \frac{1}{(n-1)!}t^{n-1}e^{-a t} \right] = \frac{1}{(s+a)^{n}}\\ }$

覚えておきたいラプラス変換　よく入力信号として登場するもの

入力の関数として頻繁に登場する、インパルス関数、ステップ関数、ランプ関数、正弦波関数（余弦波関数）は抑えて欲しいものです。

インパルス関数（ディラックのデルタ関数）

$\displaystyle{ \mathcal{L}\left[ \delta(t) \right] = 1\\ }$

ステップ関数

$\displaystyle{ \mathcal{L}\left[ u(t)\ \rm{or}\ 1\right] = \frac{1}{s}\\ }$

ランプ関数

$\displaystyle{ \mathcal{L}\left[ t \right] = \frac{1}{s^2}\\ }$

正弦波関数（角周波数 $\omega$ ）

$\displaystyle{ \mathcal{L}\left[ \sin(\omega t) \right] = \frac{\omega}{s^2+\omega^2}\\ }$

余弦波関数（角周波数 $\omega$ ）

$\displaystyle{ \mathcal{L}\left[ \cos(\omega t) \right] = \frac{s}{s^2+\omega^2}\\ }$

とりあえずf(t)のラプラス変換はF(s)とする

原関数 $f(t)$ のラプラス変換を $F(s)$ として以下の話は進みます。

僕は

「原関数の正体は判らないけれど、とりあえず $f(t)$ としてそのラプラス変換は $F(s)$ としておきましょう。」

これが、ラプラス変換を制御の道具として使っていく際に、大変重要だと思います。

関数の記号は $f$ ですが、これが $g$ でも $h$ でも良いのです。

$\displaystyle{ \mathcal{L} [ g(t) ]=G(s)\\ }$

$\displaystyle{ \mathcal{L} [ h(t) ]=H(s)\\ }$

小文字は大文字へ

「小文字の時間関数は大文字の $s$ の関数に変換される！」

と機械的に思っておくと便利です。慣れたら、大文字、小文字変換もなしで $(t)$ 、 $(s)$ で区別したり、それも省略して文脈で判断するなど、していけます。

覚えておきたいラプラス変換の性質

線形性

ラプラス変換しても足し算は足し算、定数倍は定数倍の関係が保持されます。

$\displaystyle{ \mathcal{L}\left[ \alpha f(t) \pm \beta g(t) \right] = \alpha F(s) \pm \beta G(s)\\ }$

これは当たり前に使えないとならない、大事な大事な大事な性質です。

微分のラプラス変換

「原関数を微分したもののラプラス変換はどうなるの？」ということで、ラプラス変換で微分方程式を解くアルゴリズにとって大事な性質です。

$f(t)$ の $n$ 階微分を以下のように表現するとすると

$\displaystyle{ \frac{d^{(n)} f(t)}{d t^{(n)}}\\ }$

微分のラプラス変換は以下のようになります。

$\displaystyle{ \mathcal{L}\left[ \frac{d^{(n)} f(t)}{d t^{(n)}} \right] = s^{n} F(s) - s^{n-1} f(0) - s^{n-2} \dot{f} (0) - \cdots - \left.\frac{d f^{(n-1)}(t)}{d t^{(n-1)}} \right|_{t=0} \\ }$

積分のラプラス変換

$\displaystyle{ \mathcal{L}\left[ \int_0^t f(u) du \right] = \frac{1}{s} F(s)\\ }$

ラプラス変換の世界で積分は $\frac{1}{s}$ であると言うところが、大変重要です。

制御では積分器を

$\displaystyle{ \frac{1}{s}\\ }$

で表します。

最終値の定理

たとえば、モータに５V電圧をかけたら最終的に角速度（回転数）はいくらになるのか知りたい場合に、モータのステップ応答を表すラプラス変換された式があれば以下で示す最終値の定理を使えば簡単に求まります。

最終値の定理はラプラス変換された入出力関係がわかれば、その時間関数が判らなくても、時間領域で十分に時間が経った後の出力の定常値が以下の計算でわかると言うものです。

ただし、知りたいシステムが時間が立つと発散する場合は使えません。

$\displaystyle{ \lim_{t \to \infty} f(t)=\lim_{s \to 0} s F(s) }$

$F(s)$ の前に $s$ がかけられているので、ご注意を！

初期値の定理

最終値があれば初期値もあると言うわけで、以下のようになります。

$\displaystyle{ \lim_{t \to 0} f(t)=\lim_{s \to \infty} s F(s) }$

こちらの方は、覚えるための動機付けを説明するのが難しいですね。

推移定理

前回の台形の入力の話題の際に触れました

時間関数をT時間遅れにする場合は

$\displaystyle{ \mathcal{L} [f(t-T)]=F(s) e^{-Ts}\\ }$

ちなみに、 $T\gt0$ です。

合成積（コンボリュージョン）のラプラス変換

「おお！コンボリュージョン！」（おおげさに）

「あなたは、どうして、コンボリュージョンなの！」（おおげさに）

（※注　この記事をずーっと書き続けて、つかれてナチュラルハイになって来ています。）

ラプラス変換の説明をしていて、このコンボリュージョンを、説明するのはかなり難しく困難な仕事です。説明しないと言う選択肢もありかもしれません。

関数 $g(t)$ と $u(t)$ の下の式に示すような掛け算をして時刻０から所望する時間tまでの積分を考えたときに、これを合成積または畳み込み積分（コンボリュージョン）と言います。積分を書くのが大変なので $g*f$ で表します。

式そのまま、日本語で説明しただけです。すみません。

$\displaystyle{ g*u:=\int_0^t g(t-\tau) u(\tau) d \tau\\ }$

その、ラプラス変換はどうなるのか？という問に対する答えが以下です。

$\displaystyle{ \mathcal{L}\left[\int_0^t g(t-\tau) u(\tau) d \tau \right]=G(s) U(s)\\ }$

ラプラス変換の性質として合成積のラプラス変換はこうなると言うことです。

卵と鶏論争になりそうですが・・・・（今の僕の限界）

関数 $g(t)$ は $G(s)$ の逆ラプラス変換と考えますと時刻０より前はありません。従って、 $g(t-\tau)$ の値は時刻 $\tau$ より前は無ということになります。そして、時刻 $\tau$ の時に $g(0)$ の値が出力されます。つまり $g(t-\tau)$ は時間軸上で、 $\tau$ 時間だけ $g(t)$ が遅れて出てくることになります。

そこで、以下の式について考えると

$\displaystyle{ g(t-\tau) u(\tau) \\ }$

これは、関数 $g(t-\tau)$ を時刻 $\tau$ の時の $u(t)$ の大きさである、 $u(\tau)$ 倍するということです。

これを $\tau$ を小さく刻みながら、０からtまで積分するのが合成積ということになります。

後ほど明らかにしたいと思いますが、制御工学では $g(t)$ はシステムの単位インパルス応答を表します。また、 $u(t)$ は任意の入力で時々の入力の大きさを表します。つまり、合成積は単位インパルス応答に時々の入力の大きさをかけて全時間について足し合わせ（積分）たものとなりシステムの任意入力に対する応答を表します。

そして、重要なのはスステムの任意入力に対する応答のラプラス変換が

$\displaystyle{ G(s) U(s) \\ }$

ということです。

$G(s)$ は伝達関数とよばれ、システムの性質を表し、システムの数学モデルである微分方程式から導かれます。また、 $U(s)$ は任意入力の時間関数のラプラス変換です。

合成積はシステムの応答を表します。それは、難しそうな積分です。しかし、ラプラス変換された世界では、システムと入力の単なる積になります。

これが伝達関数を使った制御工学の大事なところとかなり密接に関係していると思うのですが、まだ、伝達関数やその時間応答についての話をしていないので、中途半端な説明になってしまいました。

合成積については、今後もう少しきっちり説明できたらなと思います。

以上で、ラプラス変換の性質についてはおわります。

嗚呼！結局

結局、教科書的に網羅した感じになってしまいました。当、「Pythonで考えるDCモータの制御」的に大事なところに焦点を当ててすっきりしたかったかなと思いますが今回はこんな感じでお茶を濁します。

有利関数

以下のように、多項式関数の分数式を有利関数と呼びます。

$\displaystyle{ F(s) = \frac{a_1 s^m + a_2 s^{m-1} + \cdots + a_m }{b_1 s^n + b_2 s^{n-1} + \cdots + b_n } \\ }$

有利関数の極と零点

上記のように一般に有利関数は分母分子ともにsの多項式になり、分母＝０の式を特性方程式と呼び、その根を極と呼びます。分子＝０の式には名称はありませんが、その根は零点と呼ばれます。分子と分母と次数は通常

分母の次数≧分子の次数

となります。

極や零点は多項式方程式の解なので複数あります。零点を $z_x$ 、極を $p_x$ 分子をm次、分母をn次とすると、有利関数 $F(s)$ を次のように表す事ができます。

${\displaystyle F(s)=\frac{a_1(s-z_1)(s-z_2) \cdots (s-z_m)}{b_1(s-p_1)(s-p_2)\cdots(s-p_n)}\\ }$

この際、重複する極や重複する零点がある場合があります。また共通の極と零点がある場合は約分されているとします。

有利関数の逆ラプラス変換と部分分数展開

制御工学では有利関数 $F(s)$ の逆ラプラス変換が求められる事が重要になります。

そして、有利関数は必ず部分分数展開ができます。

重複する極がない場合

もし、重複する極がないとしますと、上記の式は、以下のような形にする事ができます。

${\displaystyle F(s)=\frac{K_1}{s-p_1}+\frac{K_2}{s-p_2} + \cdots + \frac{K_n}{s-p_n}\\ }$

$K_x$ は求めなければなりませんが、ヘビサイドの展開定理により容易に求める事が可能ですので、求められたとします。

部分分数展開したのは、上式の右辺の各項が覚えておきたいいラプラス変換のランキングNo.1が適用できるからです。

この式を逆ラプラス変換する際に、No.1を適用すると、簡単に以下の時間関数が得られます。

${\displaystyle f(t)=K_1 e^{p_1 t} + K_2 e^{p_2 t} + \cdots + K_n e^{p_n t}\\ }$

この様に、極が重複しない場合の有利関数の逆ラプラス変換は必ず、指数関数の和になります。

重複する極がある場合

有利関数の極 $p_i$ が $j$ 重解だったとします。すると有利関数は次のように表現できます。

${\displaystyle F(s)=\frac{a_1(s-z_1)(s-z_2) \cdots (s-z_m)}{b_1(s-p_1)(s-p_2)\cdots(s-p_{n-i})(s-p_i)^j}\\ }$

上式の部分分数展開は以下のようになります。

${\displaystyle F(s)=\frac{K_1}{s-p_1}+ \cdots + \frac{K_{i-1}}{s-p_{i-1}} + \frac{K_{ij}}{(s-p_i)^j}+ \frac{K_{i(j-1)}}{(s-p_i)^{(j-1)}}+\cdots+\frac{K_{i1}}{s-p_i}\\ }$

係数 $K_{xy}$ についても容易に求める事ができます。故に、上式に逆ラプラス変換を施すと時間応答の式が求められ、その際には覚えておきたいいラプラス変換のランキングNo.2を適用することで簡単に求められます。

逆ラプラス変換した結果を以下に示します。

${\displaystyle f(t)=K_1 e^{p_1 t} + \cdots + K_{i-1} e^{p_{i-1} t} \\+ K_{ij} \frac{1}{(j-1)!}t^{j-1}e^{p_j t}+ K_{i(j-1)} \frac{1}{(j-2)!}t^{j-2}e^{p_j t}+\cdots+ K_{i1} e^{p_j t}\\ }$

部分分数展開における係数の求め方

極に重複がない場合

極に重複がない場合は極 $p_i$ に対応する係数 $K_i$ は以下のようにして求めます。

$\displaystyle{ K_i=\lim_ {s \to p_i}(s-p_i)F(s)\\ }$

極限内は $\frac{0}{0}$ になりますが、分母分子に $s-p_i$ が存在して約分できるので、 $\frac{0}{0}$ は回避できます。

しかし、これは手計算している場合や数式処理では容易にできる事ですが、数値計算上はあまりシンプルには行かないので、 ロピタルの定理を適用すると数値演算でもわりとシンプルに計算をする事ができます。

どういうことかと言いいますと、有利関数を次のように簡単に表すとします。

$\displaystyle{ F(s)=\frac{N(s)}{D(s)}\\ }$

この場合は $D(s)$ は特性方程式です。係数を求める式に適用すると

$\displaystyle{ K_i=\lim_ {s \to p_i}\frac{(s-p_i)N(s)}{D(s)}\\ }$

この場合も極限内は $\frac{0}{0}$ なります。この、 $\frac{0}{0}$ になる場合に、分母と分子を両方微分して極限をとったときに、元の極限計算と結果が一致するというのがロピタルの定理です。

$\displaystyle{ K_i=\lim_ {s \to p_i}\frac{(s-p_i)N(s)}{D(s)}=\lim_ {s \to p_i}\frac{\{(s-p_i)N(s)\}^\prime}{\{D(s)\}^\prime}\\ }$

１回の微分では $\frac{0}{0}$ が解消しない場合は、微分をさらにします。

極に重複がある場合

極 $p_i$ が $k$ 重解だったとします。その時に係数 $K_{ij}$ を求める方法は以下のようになります。

$\displaystyle{ K_{ij}=\frac{1}{(k-j)!} \lim_{s \to p_i} \frac{d^{(k-j)} \{(s-p_i)F(s)\}}{d s^{(k-j)}}\\ }$

$j$ は $j=k\sim1$ の範囲の整数です。

ここで $0!=1$ とし $\frac{d^{(n)} F(s)}{d s^{(n)}}$ は $F(s)$ の $n$ 階微分とします。特殊ケースの $\frac{d^{(0)} F(s)}{d s^{(0)}}$ は微分しないということです。これは、 $j=k$ の時に起こります。

これら、係数を求める手順も、このように求めると言われても、やってみないとチンプンカンプンなものです。慣れるとパズルを解くような感覚にもなりますが、計算機にやって欲しいところかもしれません。本シリーズではPython先生にやってもらおうと考えています。

伝達関数

ここから制御の基礎の話に進みます。

伝達関数とは

ラプラス変換された世界でシステムの入出力関係を表す関数を伝達関数と呼びます。もう少し具体的にいうと、出力のラプラス変換を入力のラプラス変換で割った物です。

伝達関数を $G(s)$ 、入力を $U(s)$ 、出力を $Y(s)$ としたときそれらの間には以下のような関係があります。

$\displaystyle{ Y(s)=G(s) U(s) \\ }$

あるいは

$\displaystyle{ \frac{Y(s)}{U(s)}=G(s)\\ }$

伝達関数の求め方

伝達関数を求める手順は以下の通りです

与えられた制御対象の微分方程式を初期値０としてラプラス変換する。
ラプラス変換された式はsに関する代数方程式なので、そこから出力／入力を求める。

例題　モータの伝達関数を求める

では、実際に伝達関数を求めてみます。

もうおなじみだと思いますがDCモータの微分方程式は以下の二本です。（負荷省略）

$\displaystyle{ L\frac{d i(t)}{d t}+ Ri(t) + K\omega(t) = e(t)\\ }$

$\displaystyle{ J\frac{d \omega(t)}{d t}+ D \omega(t) = K i(t)\\ }$

以上の式をこれまでの議論に基づいてラプラス変換してみます。

$\displaystyle{ Ls I(s)+ RI(s) + K\Omega(s) = E(s)\\ }$

$\displaystyle{ Js\Omega(s)+ D \Omega(s) = K I(s)\\ }$

小文字の関数は大文字へ直せば良いんでした。また、微分は微分のラプラス変換の性質を適用しましたが、初期値０なので、大変簡単です。

さらに、整理を進めましょう

$\displaystyle{ (Ls +R )I(s) = E(s)- K\Omega(s)\\ }$

$\displaystyle{ (Js+ D) \Omega(s) = K I(s)\\ }$

電流と角速度について解くと

$\displaystyle{ I(s) = \frac{E(s)- K\Omega(s)}{Ls +R} \\ }$

$\displaystyle{ \Omega(s) = \frac{K}{Js+ D} I(s)\\ }$

１式目を２式目に代入して電流を消去します。

$\displaystyle{ \Omega(s) = \frac{K}{(Js+ D)(Ls +R)} E(s)　- \frac{K^2}{(Js+ D)(Ls +R)} \Omega(s)\\ }$

$\Omega$ について解くと

${\displaystyle \Omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} E(s)\\ }$

これがDCモータの入出力関係を表す式となります。

出力／入力に書き直すと

${\displaystyle \frac{\Omega(s)}{E(s)} =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \\ }$

上式の右辺が伝達関数です。

DCモータの定常角速度の計算

伝達関数がわかりましたので、ラプラス変換の性質にある最終値の定理を用いてDCモータの定常角速度を求めてみましょう。

DCモータに電圧e(V)引加した時、十分に時間がたったら角速度はどのくらいになるのか解いてみます。

DCモータに電圧e(V)引加は制御工学的にはどのように考えるかですが、一度e(v)を印加したら後はずっと印加し続けるので、大きさe(v)のステップ入力を与えることになります。

ステップ入力のラプラス変換表現を考慮すると入力のラプラス変換表現 $E(s)$ は以下のようになります。

$\displaystyle{ E(s)=\frac{e}{s} \\ }$

これを入出力の式に代入してみると以下のようになります。

${\displaystyle \Omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \frac{e}{s}\\ }$

最終値の定理を適用してみます。

${\displaystyle \omega(\infty) =\lim_{s \to 0 }s\Omega(s)=\lim_{s \to 0} s \frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \frac{e}{s}\\ }$

最終値の定理で掛けられるsとステップ入力の分母のsがうまく約分されて、極限が不定形にならず計算できそうです。

${\displaystyle \omega(\infty) =\lim_{s \to 0} \frac{Ke}{JL s^2 + (DL+JR)s + DR+K^2} \\ }$

計算を進めると

${\displaystyle \omega(\infty) = \frac{Ke}{DR+K^2} \\ }$

過渡応答（時間応答）の計算

定常値の計算は最終値の定理を用いると良い事がわかりました。途中の状態はどうなるのか、すなわち過渡応答、各時刻にたいする角速度はどうなるのかを調べてみることにします。

モータについてのみならず、他の対象についても制御工学では代表的な入力信号に、インパルス信号、ステップ信号、ランプ信号があります。それらを入力とした場合について検証を進めていきます。

単位インパルス応答

インパルス信号のラプラス変換は１ですのでシステムの応答式のラプラス変換は

$\displaystyle{ Y(s)=G(s) \\ }$

となって、とどのつまり伝達関数そのものになります。

つまり、伝達関数の逆ラプラス変換がインパルス応答の時間関数です。

モータの伝達関数に当てはめてみます。

${\displaystyle \Omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2} \\ }$

ここから、極を求めたときに分母は２次式ですから極は二つ求められます。２次方程式の解の公式でもPythonを用いても、容易に解は求められます。その時、極が実数でも複素数であっても今後の議論は成り立ちます。

さて、極を $p_1$ 、 $p_2$ としましょう。すると応答式は因数定理を用いて以下のように書き直せます。

${\displaystyle \Omega(s) =\frac{K}{JL(s -p_1)(s-p_2)} \\ }$

次に、この式を逆ラプラス変換して時間応答を表す式を導けば完了です。

有利関数の説明のところを参考にして、部分分数展開します。

${\displaystyle \Omega(s) =\frac{K_1}{s -p_1}+ \frac{K_2}{s-p_2} \\ }$

有利関数の逆ラプラス変換の話を振り返って、係数 $K_1$ と $K_2$ について考えてみましょう。

DCモータの場合は重解になることはまず考えられませんので重解にならない場合を適用すればOKです。

${\displaystyle K_1 =\lim_{s \to p_1}\frac{(s -p_1)K}{JL(s -p_1)(s-p_2)}=\lim_{s \to p_1}\frac{K}{JL(s-p_2)}=\frac{K}{JL(p_1-p_2)} \\ }$

${\displaystyle K_2 =\lim_{s \to p_2}\frac{(s -p_2)K}{JL(s -p_1)(s-p_2)}=\lim_{s \to p_2}\frac{K}{JL(s-p_1)}=\frac{K}{JL(p_2-p_1)} \\ }$

ここまでできてしまえば、終わったようなもので、逆ラプラス変換した時間応答の式を示します。

${\displaystyle \omega(t) =K_1 e^{p_1} + K_2 e^{p_2} \\ }$

以上からインパルス応答を求める計算で少し大変なのは、極を求める作業になります。

２次方程式

${\displaystyle a s^2 + b s +c =0 \\ }$

の解を求めるのをPython先生にやってもらうには以下のようにnumpyモジュールを使います。

import numpy as np

#a,b,cには適当に数値を代入する

np.roots([a , b, c ])

高次の多項式の解を求めるのも同様に各項の係数をリストとして与えるだけです。

ステップ応答

次はステップ応答です。

ステップ信号は $\frac{1}{s}$ ですが、定常値の時も考えたように大きさe(v)のステップ入力とすると

$\displaystyle{ E(s)=\frac{e}{s} \\ }$

となり

応答を計算する式は

$\displaystyle{ Y(s)=G(s) \frac{e}{s} \\ }$

DCモータの伝達関数を当てはめてみます。

${\displaystyle \Omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2}\frac{e}{s} \\ }$

ここから、極を求めたときに分母は３次式ですから極は３つ求められますが、一つは０なことは明らかです。３次方程式の解の公式でもPythonを用いても、容易に解は求められます。しかし、この場合は極の一つが０であとはインパルス応答で求めた極と同じであることも明らかです。

すると、０以外の極を $p_1$ 、 $p_2$ としましょう。すると応答式は因数定理を用いて以下のように書き直せます。

${\displaystyle \Omega(s) =\frac{Ke}{JL(s -p_1)(s-p_2)s} \\ }$

次に、この式を逆ラプラス変換して時間応答を表す式を導けば完了です。

有利関数の説明のところを参考にして、部分分数展開します。

${\displaystyle \Omega(s) =\frac{K_1}{s -p_1}+ \frac{K_2}{s-p_2} + \frac{K_3}{s}\\ }$

有利関数の逆ラプラス変換の話を振り返って、係数 $K_1$ 、 $K_2$ と $K_3$ について考えてみましょう。

インパルス応答の時と同様に３つの極は重解にならないので、重解にならない場合を適用すればOKです。

${\displaystyle K_1 =\lim_{s \to p_1}\frac{(s -p_1)Ke}{JLs(s -p_1)(s-p_2)s}=\lim_{s \to p_1}\frac{Ke}{JL(s-p_2)s}=\frac{Ke}{JL(p_1-p_2)p_1} \\ }$

${\displaystyle K_2 =\lim_{s \to p_2}\frac{(s -p_2)Ke}{JL(s -p_1)(s-p_2)s}=\lim_{s \to p_2}\frac{Ke}{JL(s-p_1)s}=\frac{Ke}{JL(p_2-p_1)p_2} \\ }$

${\displaystyle K_3 =\lim_{s \to 0}\frac{s Ke}{JL(s -p_1)(s-p_2)s}=\lim_{s \to 0}\frac{Ke}{JL(s-p_1)(s-p_2)}=\frac{Ke}{JLp_1 p_2} \\ }$

続いて、逆ラプラス変換した時間応答の式を示します。 $e^{0}$ は１ですので第三項は少し趣が違いますが、同じことをしています。

${\displaystyle \omega(t) =K_1 e^{p_1} + K_2 e^{p_2} + K_3\\ }$

結局南極、ステップ応答の場合も極がわかれば応答がわかるということになります。

ランプ応答

次はランプ応答です。

ランプ信号は $\frac{1}{s^{2}}$ ですが、傾き $\alpha$ のランプ入力とすると

$\displaystyle{ E(s)=\frac{\alpha}{s^2} \\ }$

となり

応答を計算する式は

$\displaystyle{ Y(s)=G(s) \frac{\alpha}{s^2} \\ }$

DCモータの伝達関数を当てはめてみます。

${\displaystyle \Omega(s) =\frac{K}{JL s^2 + (DL+JR)s + DR+K^2}\frac{\alpha}{s^2} \\ }$

ここから、極を求めたときに分母は4次式ですから極は4つ求められますが、2つは０の重解であることは明らかです。４次方程式の解はPythonを用いて容易に求められます。しかし、この場合は極の２つが０であとはインパルス応答で求めた極と同じであることも明らかです。

すると、０以外の極を $p_1$ 、 $p_2$ としましょう。すると応答式は因数定理を用いて以下のように書き直せます。

${\displaystyle \Omega(s) =\frac{Ke}{JL(s -p_1)(s-p_2)s^2} \\ }$

次に、この式を逆ラプラス変換して時間応答を表す式を導けば完了です。

有利関数の説明のところを参考にして、部分分数展開します。ここでは極に重解が含まれるので、今までとは趣が違いますよ！

${\displaystyle \Omega(s) =\frac{K_1}{s -p_1}+ \frac{K_2}{s-p_2} + \frac{K_{32}}{s^2} ++ \frac{K_{31}}{s} \\ }$

有利関数の逆ラプラス変換の話を振り返って、係数 $K_1$ 、 $K_2$ 、 $K_{32}$ と $K_{31}$ について考えてみましょう。

こんかいは極に２重解が含まれるので、重解にならない場合と重解になる場合の両方考えます。

$K_1$ 、 $K_2$ は重解にならない場合で、今までと同様です

${\displaystyle K_1 =\lim_{s \to p_1}\frac{(s -p_1)K\alpha}{JLs(s -p_1)(s-p_2)s^2}=\lim_{s \to p_1}\frac{K\alpha}{JL(s-p_2)s^2}=\frac{K\alpha}{JL(p_1-p_2)p_1^2} \\ }$

${\displaystyle K_2 =\lim_{s \to p_2}\frac{(s -p_2)K\alpha}{JL(s -p_1)(s-p_2)s^2}=\lim_{s \to p_2}\frac{K\alpha}{JL(s-p_1)s^2}=\frac{K\alpha}{JL(p_2-p_1)p_2^2} \\ }$

$K_{32}$ 、 $K_{31}$ は重解の極に対応するので、重解の場合を適用します。

${\displaystyle K_{32} =\frac{1}{0!}\lim_{s \to 0}\frac{s^2 K\alpha}{JL(s -p_1)(s-p_2)s^2}=\lim_{s \to 0}\frac{K\alpha}{JL(s-p_1)(s-p_2)}=\frac{K\alpha}{JLp_1 p_2} \\ }$

${\displaystyle K_{31} =\frac{1}{1!}\lim_{s \to 0} \left\{\frac{s^2 K\alpha}{JL(s -p_1)(s-p_2)s^2}\right\}^\prime=\lim_{s \to 0} \left\{\frac{K\alpha}{JL(s-p_1)(s-p_2)}\right\}^\prime \\ }$

${\displaystyle =\lim_{s \to 0} \frac{-K\alpha \{JL(s-p_1)(s-p_2)\}^\prime}{\{JL(s-p_1)(s-p_2)\}^2}=\lim_{s \to 0} \frac{-KJL\alpha \{(s-p_2)+(s-p_1)\}}{\{JL(s-p_1)(s-p_2)\}^2}\\ }$

${\displaystyle = \frac{KJL\alpha (p_1+p_2)}{(JLp_1p_2)^2}=\frac{K\alpha (p_1+p_2)}{JL(p_1p_2)^2}\\ }$

続いて、逆ラプラス変換した時間応答の式を示します。

${\displaystyle \omega(t) =K_1 e^{p_1} + K_2 e^{p_2} + K_{32} t + K_{31}\\ }$

ランプ応答の場合も極がわかれば応答がわかるということになります。

解と係数の関係

補足です。

以上のお話の中では、 $p_1$ と $p_2$ は以下の２次方程式の答えです。

${\displaystyle s^2 + \frac{DL+JR}{JL}s + \frac{DR+K^2}{JL} =0 }$

２次方程式の解と係数の関係を適用すると

${\displaystyle p_1+p_2=\frac{DL+JR}{JL} }$

${\displaystyle p_1p_2=\frac{DR+K^2}{JL} }$

これらを、今までの結果に代入しても良いです。

ブロック線図と伝達関数

制御システムを考える時にラプラス変換の世界とベストマッチする図示表現がブロック線図です。ブロック線図を読み解くルールはわずかに３つです。以下にそれを示します。

伝達要素

ブロック線図は四角（ブロック）と矢印または線で構成されます。

ブロックは伝達特性を表して伝達関数や定数を意味していて、上の図のようにブロックが伝達関数 $G(s)$ を表しているとすると上の図の意味するところを数式で表現すると以下のようになります。

${\displaystyle Y(s)=G(s)U(s)\\ }$

ブロックに信号 $U(s)$ が入力され、ブロックが表している伝達特性 $G(s)$ がそれに掛けられたものが出力 $Y(s)$ であることを表しています。

加え合わせ点

白丸に入って行った信号を足したり引いたりしたものが出ていく信号であることを表しています。足すか引くかは、入っていく信号の矢印の白丸に近いところに書き添えます。

例の図では２本の信号が入って１本の信号が出ていますが、２本以上入る場合を考える場合は複数の矢印を白丸につなげて良いです。

例の図が表していることを、数式で表すと以下のようになります。

${\displaystyle e(s)=R(s)-Y(s)\\ }$

分岐

同じ信号を複数に分けたい時に、図のようにします。

これも、分けたいだけ、黒丸から矢印をはやしても良いです。

ブロック線図例題

このブロック線図が表す入力 $R(s)$ から出力 $Y(s)$ までの伝達関数を求めてみましょう。

考える際は、出力に注目して、そこから始めるのが簡単です。

出力 $Y(s)$ を辿ると、加え合わせ点に入っていきます。ここで $R(s)$ と加え合わせが起きます、その結果は

${\displaystyle U(s)=R(s)-Y(s)\\ }$

である事がブロック線図から読み解けます。 $U(s)$ はブロックに入っていきますから、図から以下のことが、わかります。

${\displaystyle Y(s)=G(s)U(s)\\ }$

$U(s)$ に加え合わせ点での結果を代入すると

${\displaystyle Y(s)=G(s)(R(s)-Y(s))\\ }$

この式を $Y(s)$ について解く事で、目的が達成されます。つまり以下のようになります。

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

どんな複雑なブロック線図でも、この方法論で伝達特性を読み解く事ができます。

制御の教科書にたくさん載っているブロック線図の読み替えは知っていると便利ですが、一生懸命に覚える必要はないと思います。やっていれば自然と覚えます。

フィードバック制御

出力を目標と比較して誤差があれば、誤差に応じて何かをして、目標を達成しようとすることをフィードバック制御と呼びます。

制御工学の根本です。

DCモータの制御を極シンプルにして考えると、DCモータの伝達関数を $G$ （ $(s)$ は省略しています。）設計したコントローラーの伝達関数を $C$ として、最も単純なフィードバック制御を下図に表しました。

出力は $\omega(s)$ 、目標値は $R(s)$ 、モータへの制御入力は $U(s)$ 、誤差は $e(s)$ となっています。

伝達関数を考える時に目標値から出力だけではなく、制御入力や誤差についても求めることができます。以下に、これらについて、先に述べた方法を適用して、それぞれの伝達関数を求めた場合を示します。

目標から出力までの伝達関数

${\displaystyle \frac{CG}{1+CG}\\ }$

目標から制御入力までの伝達関数

${\displaystyle \frac{C}{1+CG}\\ }$

目標から誤差までの伝達関数

${\displaystyle \frac{1}{1+CG}\\ }$

これらは、分母は同じですから、全ての場合で極は同じであることが分かります。

これらの伝達関数を展開して、有利関数の逆ラプラス変換をすれば、フィードバック制御の時間応答が求められます。

すべて、同じ工程で検討することができます。

Pythonでの多項式関数の取り扱い

Pythonで有利関数を扱う際に知っておくと便利な多項式関数の扱い方について、説明します。

たとえばsの多項式

${\displaystyle f(s)=s^2 + 2s +3 }$

${\displaystyle g(s)=4s +5 }$

をpythonで表現するときは

import numpy as np

f=np.poly1d([1, 2, 3])
g=np.poly1d([4, 5])

とします。

これでfやgについていろいろな事ができます

f=0の根

fの根を求めるには

f.r

多項式の和や積

期待どうりに和や積が求められます。

f+r
g*r

値を求める

例えばs=2の時の多項式の値を評価したければ

f(2)
g(2)

とできます。

多項式の微分

微分もできます

f.deriv()

多項式の積分

積分もできます

f.integ()

以上を駆使することにより、有利関数の逆ラプラス変換において、係数の算出などで威力を発揮できます。後日、使用例などを記事にしたいと思っています。

おわりに

自問自答しました

すごい時間をかけて書きましたが、制御工学をブログで語る必要があるかと自問自答しはじめました。（教科書見ればいいのでは？という疑問が湧きました。）

しかし、本ブログにおいて、基本的な部分の説明が足りていないので記事が完結していないと感じていました。リファレンス的なものが必要と感じ、頑張って書いてみました。

非常に疲れましたが、自分の頭の整理にはすごくよかったのです。

PI制御は書き加えようかな

あと、PI制御について書き下したものを書き加えるかどうか悩んでいます。 PI制御については、この記事を下敷きにすれば、他の記事で取り扱っているので、それで良いかなと思いもします。

周波数応答も書かなきゃ

野望が湧いてくるのですが、どうしようかな。

Pythonについて

Pythonの部分が少なくて申し訳ないのですが、解析的に作業を進めていくと結局、「極さえ求められばいいよね」ということになってしまい、Pythonの登場の余地が少なくなります。これについてはPythonに慣れていない人のためにグラフを描く方法などを加えたらいいのかなと思っています。

最後の方に書いた、多項式をpythonで取り扱う方法は、解析をすっ飛ばして数値計算でなんとかするための解を与えていると思いますので、もう少し使いこなしていきたいです。

だれが読むのかと、またまた自問自答

しかし、制御を判っている人には、読まなくてもいい記事だし、まだ初心者の人には舌足らずで、ただ長いだけで読みたくない記事になっていると思います。もう少し、よくしたいのですが、今後の課題です。

ようやく台形制御へ！

寄り道しましたが、ようやく台形入力に対するPI制御の応答を書く準備が整いました。

これからも細々と続けたいと思います。

参考文献

あまり関係ないのもありますが、僕の制御の先生や先輩の本なので上げておきます。

最後の本はこんなに高いのかとびっくり

一連のモータ制御シリーズは以下となります。

rikei-tawamure.com

はじめに

制御工学に関する計算とPython

余談

今回は制御についてです

ラプラス変換と逆変換

ラプラス変換

ラプラス逆変換

ラプラス変換はしんどい？！

道具としてのラプラス変換：微分方程式を解く道具

微分方程式を代数方程式へ

シンプルな操作の代償としてラプラス変換の手間をかける

入力も掛け算するだけで表現可能

覚えておきたいラプラス変換

覚えておきたいラプラス変換ベスト２

覚えておきたいラプラス変換 よく入力信号として登場するもの

とりあえずf(t)のラプラス変換はF(s)とする

小文字は大文字へ

覚えておきたいラプラス変換の性質

線形性

微分のラプラス変換

積分のラプラス変換

最終値の定理

初期値の定理

推移定理

合成積（コンボリュージョン）のラプラス変換

有利関数

有利関数の極と零点

有利関数の逆ラプラス変換と部分分数展開

重複する極がない場合

重複する極がある場合

部分分数展開における係数の求め方

極に重複がない場合

極に重複がある場合

伝達関数

伝達関数とは

伝達関数の求め方

例題 モータの伝達関数を求める

DCモータの定常角速度の計算

過渡応答（時間応答）の計算

単位インパルス応答

ステップ応答

ランプ応答

解と係数の関係

ブロック線図と伝達関数

伝達要素

加え合わせ点

分岐

ブロック線図例題

フィードバック制御

目標から出力までの伝達関数

目標から制御入力までの伝達関数

目標から誤差までの伝達関数

Pythonでの多項式関数の取り扱い

f=0の根

多項式の和や積

値を求める

多項式の微分

多項式の積分

おわりに

自問自答しました

PI制御は書き加えようかな

周波数応答も書かなきゃ

Pythonについて

だれが読むのかと、またまた自問自答

ようやく台形制御へ！

参考文献

覚えておきたいラプラス変換　よく入力信号として登場するもの

例題　モータの伝達関数を求める