制御工学における無駄時間要素をパデ近似(3次/3次)したときのランプ応答について

皆さんよろしくお願いいたします。
標題の時間応答を求めようとしています。
無駄時間要素を無駄時間Lとすると伝達関数はG(s)=e^(-Ls)と表わせます。
この伝達関数をパデ近似(3次/3次)すると次式で表わせます。
e^(-Ls)={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}/{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120}
この分母と分子の3次方程式を解くと1つの実根aと2つの虚根σ±jωが得られます。
すると上式は次式で表わせます。
G(s)={(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)}
この場合のランプ応答r(t)={G(s)/s^2}を求めようとしております。
このr(t)の計算結果を計算で求めたご経験のある方がいらっしゃいましたら、
計算方法と結果をご教示いただきたくお願いいたします。
小生の不足している経験でできることして、部分分数に展開したり
留数定理を用いて計算しようとしておりますがうまくいきません。

No.1 回答者： noname#101087 回答日時： 2008/04/17 09:04

>G(s)={(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)} この場合のランプ応答r(t)={G(s)/s^2}を求めようとしております。

....

部分分数展開が山場のようです。その一例。

　　F(s) = G(s)/s^2 = (s-1)(s-1-j)(s-1+j)/{(s+1)(s+1-j)(s+1+j)s^2}　　……(1)
　　　　　= A/(s+1) + B/(s+1-j) + B~/(s+1+j) + C1/s + C2/s^2　　……(2)

まず C2 を求める。
　　C2 = [(s^2)*F(s)]_s=0 = -1
　　差し引き勘定。
　　E(s) = F(s) - C2/s^2 = (2s^2+8)/{(s+1)(s+1-j)(s+1+j)s}

あとは、E(s) の部分分数展開。
　　E(s) = A/(s+1) + B/(s+1-j) + B~/(s+1+j) + C1/s　　……(3)
　　(B~ は B の共役複素数)
たとえば (3) に (s+1) を乗じたあと、s=-1 を代入して A を求める。
　　A = [(s+1)*E(s)]_s=-1 = -10
同様に、
　　B = [(s+1-j)*E(s)]_s=<-1+j> = -4+2j
　　C1 = [s*E(s)]_s=0 = 4

一般論については「Laurant 展開」を検索・参照のほどを..... 。
　

この回答へのお礼

ご回答頂きありがとうございます。
結果のランプ応答についてエクセルでのご確認まで頂きありがとうございます。
しかし小生の質問が悪く、申し訳有りません。
小生が求めたいのは、むだ時間要素の3次/3次のパデ近似した
e^(-Ls)={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}/{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120}
の一般解のランプ応答を求めたいのです。
今分かっているところは、以下のとおりです。

上式を分母=0、分子=0とし3次方程式をカルダノ法でエクセルで
解くと、1つの実根aと2つの虚根σ±jωが得られることが分かりました。
しかも分母の根（極）と分子の根（零点）を比較すると
複素平面上で極が左半面、零点が右半面で左右対称になります。
そこで一般化した解は
G(s)={(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)}
とおけるところまでは分かりました。
このとき3次の各係数を代入すると上式の解を求めるエクセル式
は次の計算方法で行いました。
http://oshiete1.goo.ne.jp/qa3943927.html
一例としてL=1のとき
分母の根は
x1=-4.6443707092516、x2,x3=-3.67781464537391±3.50876191956745i
分子の根は
x1=4.644371011781658、x2,x3=3.677814734101803±3.508761733998084i
となります。

ランプ応答なので
G(s)/s^2={(s-a)(s-σ-jω)(s-σ+jω)}/{s^2(s+a)(s+σ-jω)(s+σ+jω)}
={(s-a)((s-σ)^2+ω^2)}/{s^2(s+a)((s+σ)^2+ω^2)}
=A/s+B/s^2+C/(s+a)+{D(s+σ)+Eω}/{(s+σ)^2+ω^2)}
として計算を行っておりますが、膨大な計算に悪戦苦闘しており
もっと簡単な方法があればと悩んでおります。
なにか良い智慧がございましたら、ご教示いただきたく
よろしくお願いいたします。

お礼日時：2008/04/19 22:58

No.2 回答者： noname#101087 回答日時： 2008/04/17 13:57

ネットの神「Laurant ではないぞよ。

Laurent じゃ！」
小生「ハハアーッ！」

---------------------------
>ローラン級数 - Wikipedia
http://ja.wikipedia.org/wiki/%E3%83%AD%E3%83%BC% …

No.3 回答者： noname#101087 回答日時： 2008/04/19 08:20

乗りかかった船です。

誤算訂正を兼ねて、もよりの船着場まで。

まずは誤算訂正。
　　A = -10
　　B = 3-j
　　C1 = 4
　　C2 = -1

これを使い、
　　r(t) = -t + 4 -10*exp(-t) +2*exp(-t)*{3*cos(t) - sin(t)}
を EXCEL に描かせてみた。
遅延式(オールパス関数)がいい加減なので、立ち上がり前で若干波打ちますが、ラムプを時間シフト(およそ4)したような応答。
(オールパス関数に Beesel 多項式でも使えば、もっと滑らかな立ち上がりになるでしょう。奇数次なので、極性は反転してます)
　

No.4 回答者： noname#101087 回答日時： 2008/04/19 08:55

ネットの神「Beesel って何じゃ。

Bessel だぞよ！」
小生「ヘヘエーッ！」
　

No.5 回答者： noname#101087 回答日時： 2008/04/19 22:21

>e^(-Ls)の伝達関数をパデ近似(3次/3次)すると次式で表わせます。

>e^(-Ls)={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}/{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120} ....

EXCEL で NEWTON 求根して、一次項で規準化(normalize)してみました。Bessel 多項式と一致。
　　(s+1.00)(s+0.792+j*0.756)(s+0.792-j*0.756)
　

No.6 回答者： noname#101087 回答日時： 2008/04/20 07:43

>ランプ応答なので

>G(s)/s^2={(s-a)(s-σ-jω)(s-σ+jω)}/{s^2(s+a)(s+σ-jω)(s+σ+jω)}={(s-a)((s-σ)^2+ω^2)}/{s^2(s+a)((s+σ)^2+ω^2)}
>=A/s+B/s^2+C/(s+a)+{D(s+σ)+Eω}/{(s+σ)^2+ω^2)}
>として計算を行っております

G(s)/s^2 の「何」を「計算」なさっているのですか？

・G(s)/s^2 の部分分数展開なら ANo.1 & 3 なので、これではなさそう。
・原題のランプ応答 r(t)={G(s)/s^2} を、{G(s)/s^2}のラプラス逆変換のことだとすれば、結果は ANo.3 なので、これでもなさそう。

はて..... ?
　

この回答へのお礼

ご回答頂きありがとうございます。
質問が下手で申し訳有りません。
>G(s)/s^2 の「何」を「計算」なさっているのですか？
については、
G(s)/s^2={(s-a)(s-σ-jω)(s-σ+jω)}/{s^2(s+a)(s+σ-jω)(s+σ+jω)}
={(s-a)((s-σ)^2+ω^2)}/{s^2(s+a)((s+σ)^2+ω^2)}
=A/s+B/s^2+C/(s+a)+{D(s+σ)+Eω}/{(s+σ)^2+ω^2)}
として部分分数展開し各係数A,B,C,D,Eを
a,σ，ωで表わす一般解を求めようとしています。
A,B,C,D,Eが求まれば、逆ラプラス変換し
L^-1[G(s)/s^2]=A+B*t+C*EXP(-at)+D*EXP(-σt)cos(ωt)＋E*EXP(-σt)sin(ωt)
として求めようと考えてます。
これにより、パデ近似したe^(-Ls)について、
1)いかる無駄時間Lに対しても3次方程式の各係数が求まる。
2)1)より、いかなるLに対しても3次方程式の極が求まる。
3)2)より、いかなるLに対しても部分分数展開の各係数が求まる。
よって、無駄時間係数がどんな場合でもランプ応答を
グラフに描けるようになると考えてます。
お分かりいただけたでしょうか。

お礼日時：2008/04/22 16:06

No.7 回答者： noname#101087 回答日時： 2008/04/22 18:33

>無駄時間係数がどんな場合でも.....

なるほど、スケーリングの問題ですね。いろんな処理策があります。
一例だけ。
e^(-s) について r(t) の「グラフに描けるように」しておき、e^(-Ls) の場合には、時間軸を r(t) の L 倍にする方法。
これなら、途中をいちいちやり直す手間は不要です。

------------------------------------------------
ANo.5 のコメント。

>>e^(-Ls)={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}/{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120} ....
>EXCEL で NEWTON 求根して、一次項で規準化(normalize)してみました。Bessel 多項式と一致。
>　　(s+1.00)(s+0.792+j*0.756)(s+0.792-j*0.756)

ここで、分母式のスケーリングいじりをしてました。
　　e^(-s)={1-(s/2)+(s^2/10)-(s^3/120)}/{1+(s/2)+(s^2/10)+(s^3/120)}
の分母式について、
　　s^3+12*s^2+60*s+120 = (s+4.644)(s+3.678+j*3.509)(+3.678-j*3.509)
と分解して、因数ごとに s/4.644 を s と置きなおせば、
>　　(s+1.00)(s+0.792+j*0.756)(s+0.792-j*0.756)
になるわけです。(Bessel 多項式)

この回答への補足

ご回答頂きありがとうございます。
申し訳有りませんが、ご教示頂いた内容が小生乏しい知識のため、
理解できません。申し訳ありませんが初学者にも分かりやすく
ご教示いただければ幸いです。分からないところを以下のように整理しました。
1)『e^(-s) について r(t) の「グラフに描けるように」しておき、
　e^(-Ls) の場合には、時間軸を r(t) の L 倍にする方法。
　これなら、途中をいちいちやり直す手間は不要です。』
　L-1[e^(-s)/s^2]=r(t)|_(L=1)とL-1[e^(-Ls)/s^2]=r(t)|＿(L=任意)
　の関係はどのような式で関係付けられるのでしょうか。
2)Bessel 多項式についてネットで検索してみましたが、
　理解できるものが有りませんでした。3次方程式の解を求めるとき
　ご使用されていると理解しましたが、その理解でよろしいでしょうか。
　カルダノの方法で解いても問題ないものと理解しましたがよろしいでしょうか。

補足日時：2008/04/23 12:12

この回答へのお礼

小生の方法で実施したところ問題が発生してしまい、こちらのほうに
ついてもアドバイスをいただければ幸いです。
小生が最初考えたやり方で以下のように答えを求めてみました。
　G(s)/s^2={(s-a)(s-σ-jω)(s-σ+jω)}/{s^2(s+a)(s+σ-jω)(s+σ+jω)}
　={(s-a)((s-σ)^2+ω^2)}/{s^2(s+a)((s+σ)^2+ω^2)}
　=A/s+B/s^2+C/(s+a)+{D(s+σ)+Eω}/{(s+σ)^2+ω^2)}
　についてA,B,C,Dの各係数を部分分数展開して求めました。
　上式を満たすためには以下の連立方程式を解く必要があります。
　A＋C+D=0
　(2σ+a)A+B+2σC+(a+σ)D+ωE=1
　(σ^2+ω^2+2aσ)A+(2σ+a)B+(σ^2+ω^2)C+aσD+aωE=-(2σ+a)
　a(σ^2+ω^2)A+(σ^2+ω^2+2aσ)B=σ^2+ω^2+2aσ
　a(σ^2+ω^2)B=-a(σ^2+ω^2)
　ここでσ^2+ω^2=Xとおくと上記連立方程式の答えは次のようになりました。
　A=2(X+2aσ)/(aX)
　B=-1
　C=-2(X+a^2+2a)/{a(X+a^2-2aσ)}
　D=4σ(X-a^2+2aσ)/{X(X+a^2-2aσ)}
　E=4σ(2aσ^2+σX-(a^2)*σ-2aX)/{ωX(X+a^2-2aσ)}
　これを次式に代入しました。
　L^-1[G(s)/s^2]=A+B*t+C*EXP(-at)+D*EXP(-σt)cos(ωt)＋E*EXP(-σt)sin(ωt)
　ところがエクセルでグラフを描いたら上下反転したグラフになってしまいました。
　(a=4.6443707092516、σ=3.67781464537391、ω=3.50876191956745)
　各係数をマイナスにすると求めようとするグラフになりました。
　各係数をマイナスするには、上記の連立方程式の右辺にマイナスを
　付加しなければ連立方程式が成立ちません。
　原因が分からず四苦八苦しております。お気づきになられた点など
　アドバイスいただければ幸いです。

　

お礼日時：2008/04/23 12:24

No.8 回答者： noname#101087 回答日時： 2008/04/23 13:21

ステップにわけて整理するのが良さそうですね。

[部分分数展開]
確かに「連立方程式」を立てて解くのも一策です。
当方のやり方は、展開結果を想定して各項ごとに定数 {A, B, C1, C2} を別々に求めてます。
コメント (ANo.1) では伝達関数の零点を勝手に決めてスタートしてます。

>　　F(s) = G(s)/s^2 = (s-1)(s-1-j)(s-1+j)/{(s+1)(s+1-j)(s+1+j)s^2}　　……(1)
>　　　　　= A/(s+1) + B/(s+1-j) + B~/(s+1+j) + C1/s + C2/s^2　　……(2)
>まず C2 を求める。
>　　C2 = [(s^2)*F(s)]_s=0 = -1
>　　差し引き勘定。
>　　E(s) = F(s) - C2/s^2 = (2s^2+8)/{(s+1)(s+1-j)(s+1+j)s}

このくだりでは、式(2) の左辺に s^2 を掛けて、
　　s^2*{A/(s+1) + B/(s+1-j) + B~/(s+1+j)}+ s*C1 + C2
としてから s=0 を代入すると、結果が C2 であることを利用するわけです。
これを左辺について勘定して、
　　C2 = [(s^2)*F(s)]_s=0 = [(s-1)(s-1-j)(s-1+j)/{(s+1)(s+1-j)(s+1+j)}]_s=0 = (-1)(-1-j)(-1+j)/{(+1)(+1-j)(+1+j)} = -1
を得ます。

このあと「差し引き勘定」しているのは、s が2位の極なので、1位の極の留数を別途求められるようにするためです。
ほかの {A, B, C1} も同じ考え方で求まりますのでお試しを。
この「項別代入」と、「連立方程式」とどちらがプログラムし易いかもお考えください。

[蛇足]
>エクセルでグラフを描いたら上下反転したグラフになってしまいました。

それが正解。(ANo.3 にコメント)
>(奇数次なので、極性は反転してます)

さしあたり、こんなとこでしょうか。

この回答への補足

推移定理のところ間違えました。ネットの神様ごめんなさい。
『[f(t-L)]=F(s)e^-(Ls)を利用して』ラプラス変換演算子Lが抜けてました。
正しくは『L[f(t-L)]=F(s)e^(-Ls)を利用して』です。

補足日時：2008/04/23 17:44

この回答へのお礼

ご回答頂きありがとうございます。
なるほど代数的な手法と留数定理による２つの方法を利用するわけですね。
大変参考になりました。
小生も最初、連立方程式を解くのが面倒になり、途中で留数を
使用して解こうとしましたが、EXP{(-σ±jω)t}がでてくるので、
これを更に分解して更にオイラーの定理
{EXP(jωt)+EXP(-jωt)}/2=COS(ωt)と
{EXP(jωt)-EXP(-jωt)}/(2j)=SIN(ωt)
にすることを考えると複雑な係数を計算するのが大変で挫折しました。

エクセルで描いたグラフが上下反転しているのが正解とのことですが、
知識の少ない小生の考えですが、違うと思います。
もともとむだ時間要素のランプ応答r(t)は推移定理
[f(t-L)]=F(s)e^-(Ls)を利用して
r(t)=L^-1[e^(-Ls)/s^2]=t-Lとなります。
これは無駄時間t=L後、傾き1で右上がりの直線になります。
よって、パデ近似したものもt=Lまでは多少振動しても右上がりになる
と考えられますがいかがでしょうか。

お礼日時：2008/04/23 17:33

No.9 回答者： noname#101087 回答日時： 2008/04/23 13:45

[さらなる蛇足]

>2)Bessel 多項式についてネットで検索してみましたが、理解できるものが有りませんでした。3次方程式の解を求めるときご使用されていると理解しましたが、その理解でよろしいでしょうか。
>カルダノの方法で解いても問題ないものと理解しましたがよろしいでしょうか。

3次方程式の求根と「Bessel 多項式」とは別問題。
　　1+(s/2)+(s^2/10)+(s^3/120)
に 120 を乗じ、
　　s^3 + 12*s^2 + 60*s + 120
としてから求根すると、
　　(s+1.00)(s+0.792+j*0.756)(s+0.792-j*0.756)
になり、スケーリングしてみると 3次の「Bessel 多項式」と同じなことに気づいただけのハナシです。
(直流近傍で位相直線の周波数特性をもつ)

スケーリングによる s と t の相互関係については、Laplace 変換のテキストを復習してみて...... 。

No.10 回答者： noname#101087 回答日時： 2008/04/23 18:30

>..... パデ近似したものもt=Lまでは多少振動しても右上がりになると考えられますがいかがでしょうか。

・「振動」は、勝手に作ったオールパス関数のせいです。
Bessel 多項式でオールパス関数を作れば、滑らかに立ち上がるはず。

・正負反転は、e^(-Ls) の近似有理式(オールパス関数)が奇数次だからです。
3次の例でいえば、(s-1)/(s+1) の項が正負を反転させます。
(C2 =-1 が Ramp の符号反転を明示。さらに、前便の C2 の計算過程をご覧ください。複素対の項は正負を反転させてません。実数対が正負反転の張本人)
　
>

この回答へのお礼

ご回答とご指摘をありがとうございます。
振動してしまうのは、有理関数近似を利用している都合上仕方ない
と考えてます。教科書にもうそう載っています。
問題は正負反転ですね。

e^(-Ls)の近似有理式も制御工学の教科書から抜粋したもので
以下の式で間違いないです。自分も近似計算結果を確認しました。
計算に用いたパデ近似の方法については以下のURLをご参照ください。
http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/nu …
また、3次/3次のパデ近似は載っていませんが一般式が以下のURLに
載っています。ご参考まで
http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/nu …

以上より以下の有理関数の近似式が正しいとすると
ランプ応答が正負反転してしまうのは、小生の計算方法がどこか
おかしいということと考えられますが、いかがでしょうか。
e^(-Ls)={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}/{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120}
それとも最初に
G(s)={(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)}
とおいたところに間違いがあるのでしょうか。
ANo.1のお礼に書き込みさせていただいたようにL=1のときの極と零点
は、複素平面上で左右対称であり、複素根は共役複素数なので
G(s)を上式のようにおくことは間違いないと考えたのですが。
ご指摘ご指導をよろしくお願いいたします。

お礼日時：2008/04/24 16:26