固定ステップサイズとソルバーについて

ごく簡単な2次遅れ系（0次/2次）伝達関数を用いて表される安定系システムのシミュレーション・実験を行っていますが可変ステップは実機で対応していないため、いずれもソルバータイプは固定ステップで設定しています（ソルバーは自動選択です）。

clear all;

% 2次遅れ系(1/2)
% a=0.001637;
% b=5.523e-05;
% c=1;
% d=0.01401;
% e=3.68e-05;
% plant = tf([a b],[c d e])

% 2次遅れ系(0/2)
a=0;
b=6.572e-02;
c=1;
d=16.53;
e=3.973-02;
plant = tf([a b],[c d e])

%　1次遅れ系(0/1)
% a=0;
% b=2.258;
% c=0;
% d=359.5;
% e=1;
% plant=tf([a b],[c d e])

% 目標とサンプリングタイム
r=30;
Ts = 1;

if ~mpcchecktoolboxinstalled('simulink')
    disp('Simulink(R) is required to run this example.')
    return
end

slx = 'simple_model';
open_system(slx);
sim(slx);

システム同定で同定した伝達関数を用いてシミュレーション・実験を行っていますが、固定ステップサイズを1にしたときに出力が無限大に発散してしまいシミュレーションが停止してしまいます。2次遅れ系(1次/2次)や1次遅れ系(0次/1次)の伝達関数ではシミュレーションは正常に終了します。

2次遅れ系（0次/2次）モデルにおいて固定ステップサイズのみを変更すると(0.1 , 0.25 , 0.5 , 1)、その値によって結果が収束したり、振動的になったり、無限大に発散し、シミュレーションが停止するなど、固定ステップサイズによって結果が変動します。

エラー: order_0_2_constant_simple_model (line 38)
ブロック 'simple_model/Transfer Fcn' の時間 110.0 における状態 '1' の導関数が有限ではありません。シミュレーシ
ョンを停止します。解に特異点が存在する可能性があります。そうでない場合は、(固定ステップ サイズを小さくする、または許容
誤差を厳しくすることによって) ステップ サイズを減らしてみてください

診断ビューアーでは「解に特異点が存在する場合があります。そうでない場合は、（固定ステップサイズを小さくする、または許容誤差を厳しくすることによって）ステップサイズを減らしてみてください」と表示されます。実験の都合上、固定ステップサイズをできるだけ大きな値を取りたい場合でもシミュレーションが正常に終了するためには、やはり固定ステップサイズを小さくするという方法しかないのでしょうか。ソルバーを変更してみても、固定ステップサイズが1だと、シミュレーションが停止するか、明らかに違う結果が出てしまいます（目標値30に対し、定常偏差が25残ってしまうなど）。もし根本的な解決策があれば教えて頂ければ幸いです。

また診断ビューアーに表示されているようなステップサイズが大きいと状態の導関数が有限でなくなる理由も教えて頂きたいです。

Answer 1

[Toshinobu Shintai]

おっしゃっている2次遅れ系（0次/2次）モデルのステップ応答は、添付ファイル「応答.png」のようになります。固定ステップサイズというのは、信号処理の世界で言う「サンプリングタイムステップ」にあたります。サンプリングタイムステップが1秒のとき、このモデルのダイナミクスは正確に再現することはできません。発散したりするのもそれが理由です。

一般的に、離散時間の計算でそのモデルのダイナミクスを表現するには、立ち上がり時間の10分の1のサンプリングタイムステップが必要とされています。今回の場合は0.02秒ぐらいに設定するのがよいと考えられます。

Simulinkのブロックはサブシステムごとにサンプリングタイムステップを設定することができます。フィルター部分だけを連続状態にして、それ以外のサンプリングタイムステップを例えば1秒に設定する、ということが可能です。添付のモデルを確認してください。

この2次遅れ系のモデルの計算は簡単なので、この部分だけを0.02秒にしても、計算時間にはそれほど影響しないのではないでしょうか。

[Soichiro Kato]

ご丁寧な説明ありがとうございます。

固定ステップサイズがどのようなものか、またその値の設定の目安は理解することができました。

追加の質問になってしまい恐縮ですが、質問した2次遅れ系（0次/2次）モデルのステップ応答は、添付ファイル「応答.png」のようになるとのことですが、こちらはどのように求めることができるのでしょうか。

matlabのコマンドラインからステップ応答を求めるために、

clear all;
% 2次遅れ系(0/2)
a=0;
b=6.572e-02;
c=1;
d=16.53;
e=3.973-02;
plant = tf([a b],[c d e])
step(plant)

と入力すると、添付したファイル 「step_responce.png」のようになりました。立ち上がり時間や時定数、ステップ応答波形が「応答.png」と全く異なった結果になってしまいました。

「step_responce.png」の波形であれば固定ステップサイズは1でもダイナミクスは再現されると思ったのですが、どうしてこのような違いが生じるのでしょうか。

[Toshinobu Shintai]

「simple_model.slx」の中にあるTransfer Fcnブロックの応答をscopeで表示しただけです。

スクリプトで再現したい場合は、以下のコマンドを実行してください。

a = 0;
b = 100;
c = 1;
d = 5;
e = 100;
plant = tf([a b],[c d e]);
step(plant)

[Soichiro Kato]

返信が大変遅くなり、申し訳ありません。

ステップ応答についての問題は解決できました。

丁寧に教えて頂き、ありがとうございました。

Answer 2

[Hiroumi Mita]

連続系を離散系ソルバーで動かすには、連続系モデルを離散化することをお勧めします。

c2d関数で、ココのようにモデルを連続系から離散系にすれば、遅延演算子zの式に変換できるので

離散実装は容易です。

但し、シャノンのサンプリング定理を抑える必要があります。

さて、システムを1秒で離散化する場合

サンプリング時間Ts=1[s]

サンプリング周波数Fs=1/Ts=1[Hz]

ナイキスト周波数Ns=1/2*Fs=0.5[Hz]

なので、システムに含まれる周波数は理想的でも最大　

ナイキスト周波数までしか再現されません。

あまり荒いサンプリングだと、元の波形と乖離するのはこれが理由です。

a = 0;

b = 100;

c = 1;

d = 5;

e = 100;

plant = tf([a b],[c d e]);

plant_disc=c2d(plant,1) %<=ココ

plant_disc =

1.085 z + 0.08071

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

z^2 + 0.1587 z + 0.006738

サンプル時間: 1 seconds

これを、実装するには、分子・分母をz^2で割り、離散系伝達関数の定義を調べてください。

[Soichiro Kato]

返信が大変遅くなり、申し訳ありません。

z変換について、c2dコマンドは拝見したことがありますが、シャノンのサンプリング定理までは考えていませんでした。

離散系伝達関数の定義について、調べます。

ご丁寧に教えて頂き、ありがとうございました！