皆さんこんにちは。AmedTech代表の天野です。
伝達関数という言葉を聞いたことがございますでしょうか？これは制御工学の分野で入力に対する出力の応答を数式化したものです。今回はこの伝達関数についてご紹介するとともに、excelを使ってこのシミュレーションをしてしまう手法をご紹介いたします。
伝達関数とは
世の中にはいろいろな機器・装置の類があふれています。昔はほとんどのからくりは機械仕掛けだったので、どれだけ複雑であっても分解することによってどうやって動くのかは分析が可能でた。ところがマイクロコンピュータの発達により多くの動作がソフトウェアによってなされるようになり分解しただけではどうやって動作しているのかがわからなくなってきました。その代わり、以前は目的によってハードウェアの設計を個々に行っていた開発作業が、ソフトウェアを変更するだけでハードウェアは大きく変えずに実施できるようになり、機器・装置の値段は劇的に低下しました。伝達関数とは装置をブラックボックスと考えた時に入力に対する出力を関数として数学的に表現したものです。以前のハードウェアに依存した装置開発の時代には数学的な興味から発達した分野ですが、マイクロコンピュータの発達に伴い急速にソフトウェアに組み込まれて、現在では多くの機器制御に使われています。
科学の発展は自然現象をいかに記述するかの歴史です。ここの現象をより一般的に記述しようとした学問が物理学であり、より個々の減少に注目した学問が科学です。そしてそのすべてに共通する言語として数学が生まれました。さらに産業革命とともにそれまでに学んできた事実をもとに自分たちの目的に合ったシステムを作り上げようとしたのが工学、すなわちエンジニアです。伝達関数も自然現象を記述ための数学から発展してきたものです。自然現象はその多くが微分方程式で記述できます。例を上げれば放射能の減少、人口の拡大、刺激に対する応答、惑星の運動などすべて微分方程式であらわされます。そうなってきますと、ではこの微分方程式を利用して自分たちの欲しいシステムを作ってしまえ、となるのはまあ自然の成り行きかもしれません。あるシステム入力に対する入力と出力をそれぞれy(t)、x(t)(注意：入力x、出力yではない！)とすると、一般に次のような定係数線形常微分方程式で記述されます
$$a_{n}\frac{d^nx}{dt^n}+a_{n-1}\frac{d^{n-1}x}{dt^{n-1}}+…+a_{1}\frac{dx}{dt}+a_{0}x=b_{m}\frac{d^my}{dt^m}+b_{m-1}\frac{d^{m-1}y}{dt^{m-1}}+…+b_{1}\frac{dy}{dt}+b_{0}y$$
$$n\geq m$$
このように微分方程式で記述されたシステムにある入力を入れた時にどのような出力が出てくるかは、この方程式を一定の初期条件で説いてやればよいわけですが、エンジニアは数学者ではないので、このままではとても大変、というより不可能です。
ラプラス変換
そこで、ラプラス変換という、とても便利な手法を利用します。ラプラス変換については、これもよく解説されたウェッブページがたくさんありますので、割愛しますが、一番重要な事項は上の定係数線形常微分方程式で一般化されたシステムがラプラス変換によって以下の式になることです。
$$G(s)=\frac{X(s)}{Y(s)}=\frac{b_{m}s^m+b_{m-1}s^{m-1}+…+b_{1}s+b_{0}}{a_{n}s^n+a_{n-1}s^{n-1}+…+a_{1}s+a_{0}}$$
ここでY(s)とX(s)は入力と出力のラプラス変換値でありG(s)がラプラス変換で記述されたシステムの特性、すなわち伝達関数となることです。微分方程式は時間領域と呼ばれ時間tの関数として記述されますが、ラプラス変換式ではs領域と呼ばれる以下の変換式に伴うパラメータで記述されます。
$$F(s)=\int_0^\infty f(t)e^{-st}dt$$
なんでこんな面倒くさい方法を取るかと言えば、s領域ではt領域での微分・積分が、除算・積算ですむからです。こうして計算しておいてからラプラス逆変換をすれば、面倒くさい微分方程式の解の計算が、ほぼ機械的な変換によって可能になるからです。
一時遅れ要素
この伝達関数のもっとも単純な形である
$$G(s)=\frac{X(s)}{Y(s)}=\frac{1}{Ts+1}$$
を1次遅れ要素と呼びます。ちなみに自然界ではほとんどの場合分母の次数が分子の次数より高くなることがほとんどで、分子のほうが高くなるものはわざわざ作ってやらなくてはならないそうです。したがって、分母にsが一つ、分子にsが0個のものがもっとも単純な伝達関数となります。実は、もっと単純なF(s)=1であるとかF(s)=1/sであるとかもあるのですが、前者はインパルス関数といって、時間が０で無限大の大きさを持ち、それ以外では出力0の関数をさし、後者は時間０で大きさが０から１へ変化するステップ関数を指します。これらは挙動のわからないシステムの応答を調べるのによく使われます。
さて上の1次遅れ要素ですが、これを逆ラプラス変換して微分方程式で書くと次のようになります。
$$T\frac{dx}{dt}+x=y$$
この式の意味するところは出力の時間変化分に定数をかけたものと出力を加えると入力になる、です。なぜか出力をx、入力をyと記述するので、こんがらがってしまいますよね。ではこの回路のインパルス応答とステップ応答はどうなるでしょうか？インパルス応答とは、時間0で無限大の入力を入れた時の出力の変化、ステップ応答とは時間0で1の大きさの入力を入れた時の出力の変化を指します。現実世界の電気回路では無限大の入力や瞬時に立ち上がる電源などはありませんが、ある程度コンデンサの容量があれば瞬間的にスイッチをONにすることでインパルス入力の代わりにしたり、スイッチをONにしてからの出力の変化を観察することでステップ応答を測定したりできます。数学的な解は以下の通りです。
インパルス応答
$$g(t)=\mathcal{L}^{-1}[X(s)]=\mathcal{L}^{-1}[Y(s)G(s)]=\mathcal{L}^{-1}[1\frac{1}{1+Ts}]=\mathcal{L}^{-1}[\frac{1/T}{s+(1/T)}]=\frac{1}{T}e^{-\frac{t}{T}}$$
ステップ応答
$$h(t)=\mathcal{L}^{-1}[X(s)]=\mathcal{L}^{-1}[Y(s)G(s)]=\mathcal{L}^{-1}[\frac{1}{s}\frac{1}{1+Ts}]=\mathcal{L}^{-1}[\frac{1}{s}+\frac{1}{s+\frac{1}{T}}]=1-e^{-\frac{t}{T}}$$
この結果を見るとなぜこの伝達関数が一時遅れ要素と呼ばれるかわかりやすいですね。青い線がインパルス応答なので、瞬間的に入力があって出力が上がるものの、入力はすぐになくなってしまうので、少しずつ出力が減少していき、最後にはゼロになります。一方赤い線はステップ応答なので、入力に対し出力がある一定の遅れをもって追随することがわかります。ちなみにこの、どれくらい遅れるか、を規定するTを時定数と呼びます。
エクセルによるシミュレーション
さて、このままでは単に数式のお話になってしまうので、最近のPCではほぼすべてにインストールされているマイクロソフトのエクセルでこの伝達関数をシミュレーションしてみましょう。基本的な考え方は、微分は差分で表せる性質を利用します。微分とはある時点におけるグラフの傾きを表します。傾きが大きければ微分値も大きく、傾きが0であれば微分値も0です。この性質を利用して極大値や極小値を計算したりしますね。正確な微分値はある時点における接線の傾きですが、この”ある時点”より微小に前の時間や後の時間における値を結んだ直線の傾きも接線にかなり近い値になります。
この図では赤い線が接線ですが、赤い点と青い点（微小に後ろの点+δ）を結んだ線も赤い点と緑の点（微小に前の点－δ）を結んだ点も、接線に近いといえます青い点と緑の点を結んだ線が接線に最も近いですが、2つの微小に時間のずれた値が必要になります。測定した値を後から解析するのであればどの方法を使ってもよいのですが、リアルタイムで解析、あるいは制御をしようとすると、過去の値はわかっても未来の値はわかりません。したがってこの図でいう緑の点と赤い点を利用して微分をすることになります。このように1点における微分を2点を通る直線で近似することを差分化といいます。
それでは1時遅れ要素を差分化してみましょう。
$$X(x)=Y(s)G(s)=Y(s)\frac{1}{1+Ts}$$
$$X(s)(1+Ts)=Y(s)$$
$$X(s)+TsX(s)=Y(s)$$
$$\mathcal{L}^{-1}[X(s)+TsX(s)]=\mathcal{L}^{-1}[Y(s)]$$
sX(s)は時間領域の微分に相当します
$$x(t)+T\frac{dx(t)}{dt}=y(t)$$
$$x(t)+T\frac{x(t)-x(t-1)}{\delta t}=y(t)$$
$$x(t)(1+\frac{T}{\delta t})=y(t)+\frac{T}{\delta t}x(t-1)$$
$$x(t)=\frac{1}{(\frac{T}{\delta t}+1)}(y(t)+\frac{T}{\delta t}x(t-1))$$
$$x(t)=\frac{\delta t}{T+\delta t}y(t)+\frac{T}{T+\delta t}x(t-1)$$
エクセル化においてはデータを離散化しδt秒ごとのデータを観察します。ここでx(t-1)あるいはy(t-1)とは現在の値x(t)、y(t)よりも微小時間（δt）前の値です。それでは具体的な入力方法を見てみましょう。まずはA列、B列の表題をそれぞれ”Time”と”y(t)と入力し、初期値として時間を0（秒）、y(t)とx(t)を０と入力します。また、E1にdt（刻み時間）をE2にT（時定数）を入力します。（今回はそれぞれ0.1（秒）と10にしました）
また、式で使いやすいように刻み時間dtを名前定義します。数式→名前の定義をクリックすると以下のようなウィンドウが出てきますので、このままOKを押します。時定数のTも同様に定義してください。
A3には次の数式を入れます =A2+dt そして、C3には上で作成した伝達関数の出力式を入れます=dt/(T+dt)*B3+T/(T+dt)*C2 そしてA3からC3までマウスを左クリックしたままポインタを移動させると太枠が現れるので、右下の小さな四角をクリックし、マウス左ボタンを押したままず～～～っと下（500行もあれば十分だと思います）までスクロールします。B列は入力なので、最初は0、A列が1の時に1を入れ（B12)、以降はすべて1にします。この結果のグラフがこれです。
ちゃんと１秒でステップ入力が入り、それに相当する出力が出ていますね。では、これは正しい出力か検証してみましょう。上でやった1次遅れ要素へのステップ応答と重ねてみればわかります。そのためには、立ち上がりは１秒目ではなく０秒目からにしなくてはならないので、Ｂ列はB2からすべて1にします。そしてこのグラフと前の1次遅れ要素のステップ応答のグラフを重ねてみたのがこれです。
茶色い線は差分法の値から微分方程式の解を引いたものですが、最も大きいところで2%以下です。決して小さくはありませんが、エクセルという表計算ソフトで伝達関数をシミュレートするお手軽さを考慮すると悪くないと思います。