制御工学における無駄時間要素をパデ近似(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.14

noname#101087 · ベストアンサー

>多項式の正負をどのように判定されているのか.....

単純な目視に過ぎません。
有理式のままではゴチャゴチャして焦点が定まらないので、分子多項式だけ書き出してみましょう。

L=1のとき、原式を因数分解。
　　-(s^3/120) + (s^2/10) - (s/2) + 1 = -(s-a)(s^2-bs+c)/120
この両辺に 120 を乗算すれば、
　　-s^3 + 12s^2 - 60s + 120 = -(s-a)(s^2-bs+c)　　(ac=120)

一方、mathstudy さんの因数分解の結果は、
　　s^3 - 12s^2 + 60s - 120 = (s-a)(s-σ-jω)(s-σ+jω) = (s-a)(s^2-bs+c)
これは原式を正負逆転したものですね。

お礼

ご回答頂きありがとうございます。
返信遅くなり申し訳有りません。
ご指摘頂いた分子式で眼が覚めました。
もっとも単純なミスをしていたわけですね。
求めようとしたランプ応答が描けました。
ありがとうございました。
結局小生質問の解は以下のようになりました。
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}
分母分子の３次式をそれぞれカルダノの方法でとくと１つの実根と
２つの虚根（共役複素数）が得られる。これらを、複素平面上で
分母は左半面、分子は右半面で分母分子は左右対称であることから
一般解を実根a、複素根をσ±jωとおく、分子多項式の係数の正負に
注意すると伝達関数は次のようにおける。
G(s)=-(s-a)(s-σ-jω)(s-σ+jω)/(s+a)(s+σ-jω)(s+σ+jω)
ランプ応答をr(t)とするとr(t)=L^-1[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の各係数を部分分数展開して求める。
　上式を満たす連立方程式は次式のとおり。
　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σ((a^2)*σ+2aX-2aσ^2-σX)/{ω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)

回答No.13

noname#101087

G(s)/s^2 = {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}s^2]　　…(1)
　　　　　　　　={(s-a)((s-σ)^2+ω^2)}/{s^2(s+a)((s+σ)^2+ω^2)}　　　　　　　　　…(2)
　　　　　　　　=A/s+B/s^2+C/(s+a)+{D(s+σ)+Eω}/{(s+σ)^2+ω^2)}　　　　　　 …(3)
から、
　　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σ)}
が得られた。

… (1) と (2) では、分子多項式の正負が逆転してます。
(2) と (3) では、正負のつじつまがあってます。(B=-1)
(3) の逆 Laplace 変換をグラフにすれば、(1) のグラフを正負逆転したものになるのは当然だと思います。

お礼

ご回答頂きありがとうございます。
せっかくご指摘いただいているのに、
多項式の正負をどのように判定されているのか小生、無知で申し訳
ありませんが分かりません。
判定方法をご教示いただければ幸いです。
また、
G(s)={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次方程式を分子、分母毎に解いた結果、一例としてL=1のとき
分母は
x1=-4.6443707092516、x2,x3=-3.67781464537391±3.50876191956745i
分子は
x1=4.644371011781658、x2,x3=3.677814734101803±3.508761733998084i
となることから、一般解の場合
分子x=a、σ±jω、分母x=-a、-σ±jωとおけるので
G(s)={(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)}
={(s-a)((s-σ)^2+ω^2)}/{(s+a)((s+σ)^2+ω^2)}
と表わせると考えたのですが、元の
G(s)={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}
と正負が逆転してしまったのは、なぜなのでしょうか。
ご教示いただければ幸いです。

回答No.12

noname#101087

>２つ目のF(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)
　　F(s) = G(s)/s^2 = {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}s^2]　　……(2)

分子の最低次係数は、(1) では - 、(2) では + 。分母のほうはどちらも + 。
つまり、(1) では C2 < 0 、(1) では C2 > 0 。

…というだけのハナシでした。

お礼

ご回答頂きありがとうございます。
さらに誤記訂正ありがとうございます。

F(s) = G(s)/s^2 = {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}s^2]　　……(2)
について、小生がANo.7のお礼欄で記載させていただいたように、
G(s)={(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)}
とおいて
L^-1[G(s)/s^2]=A+B*t+C*EXP(-at)+D*EXP(-σt)cos(ωt)＋E*EXP(-σt)sin(ωt)
となるよう各係数A,B,C,D,Eを連立方程式を立てて求めたのですが
どこか、間違っているところなどお気づきの点がありますでしょうか。
ご教示いただければ幸いです。

回答No.11

noname#101087

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)
これだと Ramp の符号が反転。

　　F(s) = G(s)/s^2 = {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}
これなら Ramp は符号反転しない。

…ということでした。

補足

ここで、２つ目のF(s)の右辺分母にs^2が記載されていないのは誤記でしょうか。
正しくは
F(s) = G(s)/s^2 = {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}/s^2
ですよね。
ご教示いただきたくお願いいたします。

回答No.10

noname#101087

>..... パデ近似したものも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/numeanal2/node19.html
また、3次/3次のパデ近似は載っていませんが一般式が以下のURLに
載っています。ご参考まで
http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/numeanal2/node19.html

以上より以下の有理関数の近似式が正しいとすると
ランプ応答が正負反転してしまうのは、小生の計算方法がどこか
おかしいということと考えられますが、いかがでしょうか。
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)を上式のようにおくことは間違いないと考えたのですが。
ご指摘ご指導をよろしくお願いいたします。

回答No.9

noname#101087

[さらなる蛇足]

>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.8

noname#101087

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

[部分分数展開]
確かに「連立方程式」を立てて解くのも一策です。
当方のやり方は、展開結果を想定して各項ごとに定数 {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 にコメント)
>(奇数次なので、極性は反転してます)

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

お礼

ご回答頂きありがとうございます。
なるほど代数的な手法と留数定理による２つの方法を利用するわけですね。
大変参考になりました。
小生も最初、連立方程式を解くのが面倒になり、途中で留数を
使用して解こうとしましたが、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までは多少振動しても右上がりになる
と考えられますがいかがでしょうか。

補足

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

回答No.7

noname#101087

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

なるほど、スケーリングの問題ですね。いろんな処理策があります。
一例だけ。
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 多項式)

お礼

小生の方法で実施したところ問題が発生してしまい、こちらのほうに
ついてもアドバイスをいただければ幸いです。
小生が最初考えたやり方で以下のように答えを求めてみました。
　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)
　各係数をマイナスにすると求めようとするグラフになりました。
　各係数をマイナスするには、上記の連立方程式の右辺にマイナスを
　付加しなければ連立方程式が成立ちません。
　原因が分からず四苦八苦しております。お気づきになられた点など
　アドバイスいただければ幸いです。

補足

ご回答頂きありがとうございます。
申し訳有りませんが、ご教示頂いた内容が小生乏しい知識のため、
理解できません。申し訳ありませんが初学者にも分かりやすく
ご教示いただければ幸いです。分からないところを以下のように整理しました。
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次方程式の解を求めるとき
　ご使用されていると理解しましたが、その理解でよろしいでしょうか。
　カルダノの方法で解いても問題ないものと理解しましたがよろしいでしょうか。

回答No.6

noname#101087

>ランプ応答なので
>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に対しても部分分数展開の各係数が求まる。
よって、無駄時間係数がどんな場合でもランプ応答を
グラフに描けるようになると考えてます。
お分かりいただけたでしょうか。

回答No.5

noname#101087

>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.4

noname#101087

ネットの神「Beesel って何じゃ。Bessel だぞよ！」
小生「ヘヘエーッ！」