二次遅れ系システムに​対するフィードフォワ​ード制御波形の導出

既知の二次遅れ系の伝達関数を有するシステムの逆関数を求め、目標信号（のこぎり波）に対してフィードフォワード制御信号を導出するプログラムを作成しております。

プログラムとしては、目標波形y(t)と出力波形u(t)の関係をy(t)=u(t)とすべく、目標波形y(t)のラプラス変換Y(s)に二次遅れ系の伝達関数の逆関数1/G(s)を乗じて、逆ラプラス変換するものです。しかし、Y(s)に1/G(S)を乗じた場合と何も乗じない場合（すなわち、元の目標波形の戻る）で逆ラプラス変換の結果（数式）こそ変わるのですが、fplotのグラフが変わりません。この原因がわからず困っております。

ご助力いただけますようよろしくお願いいたします。

T=0.005;
Tr=0.001;

syms t;
U(t)=2*16*triangularPulse(0,Tr,T,t)-16;

aFreq=3330;
K=1;
ext=0.015
sys=tf(K*aFreq^2,[1,2*ext*aFreq,aFreq^2]);
[num,den]=tfdata(sys);

% Convert transfer function to symbolic function
syms s;
sys_sym=poly2sym(cell2mat(num),s)/poly2sym(cell2mat(den),s);
% Generate invert transfer symbolic function
inv_sys_sym=poly2sym(cell2mat(den),s)/poly2sym(cell2mat(num),s);

Y(t)=ilaplace(laplace(U));
inv_Y(t)=ilaplace(laplace(U)*inv_sys_sym);

fplot(Y,[0 T]);
hold on;
fplot(inv_Y,[0 T]);

Answer 1

[Hiroumi Mita]

伝達関数の逆問題を解ける条件としては

上記の伝達関数sysは厳密にプロパーでなくてはいけないと思います。

上記スクリプトでは

sys =

1.109e07

-----------------------

s^2 + 99.9 s + 1.109e07

ですが、これは厳密にプロパーではありません。そこで

sys=1.109e7*(1e-6*s^2+1e-6*s+1)/( s^2 + 99.9 s + 1.109e07)

とします。これは厳密にプロパーとし、分子の2次と一次の係数の影響を微小にする意図です。

11.09 s^2 + 11.09 s + 1.109e07

------------------------------

s^2 + 99.9 s + 1.109e07

また、スクリプトの最後を少し変えます。

figure(1)

hold on;

fplot(Y,[0 T])

fplot(inv_Y,[0 T])

legend('Y','Inv Y')

hold off

[x1,y1]=fplot(Y,[0 T]);

[x2,y2]=fplot(inv_Y,[0 T]);

添付図を参照すると、この場合は正しく求められます。

[Daisuke Kusakari]

ご回答いただきありがとうございます。

sys=1.109e7*(1e-6*s^2+1e-6*s+1)/( s^2 + 99.9 s + 1.109e07)

の分子のs^2とsにかかる係数1e-6の定数を小さくするほど、プロパーを維持しつつ、影響を微小にできると解釈しました。しかし、定数を小さく変化（1e-6-->1e-9-->1e-12）させると出力されるグラフが全く変わってしまいます。これはどのように理解すればよろしいでしょうか。

何卒宜しくお願い致します。

Answer 2

[Hiroumi Mita]

sys=1.109e7*(1e-6*s^2+1e-6*s+1)/( s^2 + 99.9 s + 1.109e07)

1e-6を順次1e-9,1e-12に変えたときのsysは同じものではありません。

周波数応答を見ると顕著いが現れます。

逆問題では

1/sys=( s^2 + 99.9 s + 1.109e07)/(1.109e7*(1e-6*s^2+1e-6*s+1))

これを扱うことになりますが、分母側を見ると

1.109e7*(1e-6*s^2+1e-6*s+1))->1.109e7*(1e-9*s^2+1e-9*s+1)->1.109e7*(1e-12*s^2+1e-12*s+1)

と変遷するごとに、分母の伝達関数は定数（ゲイン）に近くなります。

そこで、1.109e7*(1e-9*s^2+1e-9*s+1)->1.109e7*(1e-12*s^2+1e-12*s+1)は、低周波領域の特性は、ゲイン＝定数、位相=0になり似ています。

しかし、1.109e7*(1e-6*s^2+1e-6*s+1))≒1e1s^2+1e1s+1.109e7

は、2次、1次の項が無視できない大きさなので

1e-6を1e-9,1e-12に変えたものとは大きく異なります。

[Daisuke Kusakari]

1e-12-->低周波領域の特性は、ゲイン＝定数、位相=0になり似ています。

指数の12を大きくすればするほど、本来求めたいシステムの逆関数に近づくことが分かりました。実験的に指数を大きくしたところ、300までは伝達関数inv_sys3分母のs^2, sの項として反映され、400を超えると伝達関数inv_sys3分母のs^2, sの項は無視されることが分かりました。ただ、いずれの指数でも逆ラプラス変換（ilaplace(laplace(U)*inv_sys_sym)）は、私の想定した波形とはかけ離れていることが不可解です。

逆ラプラス変換ilaplace(laplace(U)*inv_sys_sym)　指数300

逆ラプラス変換ilaplace(laplace(U)*inv_sys_sym)　指数400

T=0.005;
Tr=0.001;

syms t;
U(t)=2*16*triangularPulse(0,Tr,T,t)-16;

aFreq=3330;
K=1;
ext=0.0147;
sys=tf(K*aFreq^2,[1,2*ext*aFreq,aFreq^2]);
[num,den]=tfdata(sys);

%inv_sys1=tf([1,2*ext*aFreq,aFreq^2],K*aFreq^2*[1E-6,1E-6,1]);
%inv_sys2=tf([1,2*ext*aFreq,aFreq^2],K*aFreq^2*[1E-9,1E-9,1]);
inv_sys3=tf([1,2*ext*aFreq,aFreq^2],K*aFreq^2*[1E-300,1E-300,1]);%1E-300の指数300を変更！
%bode(inv_sys1,inv_sys2,inv_sys3,{1e-1,1e7});
[inv_num,inv_den]=tfdata(inv_sys3);

% Convert transfer function to symbolic function
syms s;
sys_sym=poly2sym(cell2mat(num),s)/poly2sym(cell2mat(den),s);
% Generate invert transfer symbolic function
inv_sys_sym=poly2sym(cell2mat(inv_num),s)/poly2sym(cell2mat(inv_den),s);

Y(t)=ilaplace(laplace(U));
inv_Y(t)=ilaplace(laplace(U)*inv_sys_sym);

figure(1)
hold on;
fplot(Y,[0 T]);
fplot(inv_Y,[0 T])

Answer 3

[Hiroumi Mita]

訂正します。厳密にプロパーではないですね。

フィードフォワードのような逆問題が解けるタイプの伝達関数は分母と分子の次数が同じが正です。

さて、MATLABは数値計算機ですので、MATLAB自身が扱える数値の最小値は

eps

ans =

2.2204e-16

です。これより、小さい表現は無意味です。

指数の12を大きくすればするほど、より真値に近いということですが、

これは、以下の係数が小さくなれば、sysは本来のsys1に近くなるので当然です。

sys=1.109e7*(係数*s^2+係数*s+1)/( s^2 + 99.9 s + 1.109e07)

sys1=1.109e7( s^2 + 99.9 s + 1.109e07)

しかし、sysは、sys1を近似した表現で、sysそのものではありません。

或る程度、誤差が乗ることは前提です。係数を1e-300のような、MATLABの数値表現を範囲を超える値に

しても意味がありません。

この方法は、フィードフォワードのような逆問題を解くための近似的手段とご理解ください。

真値を求められる類の問題ではないと思います。