Scilab、XCOS

Scilab、Xcosの使い方について。

Bode線図やNyquist線図の作成や、単位ステップ応答の挙動を見るのに手軽で使い勝手が良いです。

(勿論コスト面で余裕があればMatlabの方が良いですが、、、)

Scilabの基本機能

基本コマンド

・コンソール画面への出力

disp()コマンドを使用

disp("Hello")

出力 : Hello

・コメントアウト

入力 : //

・変数の代入

 例えば、a、bという2変数にそれぞれ1、3を代入する場合、

Variables

 

・基本の論理演算子

 

・複素数の表現  

複素数、共役複素数、複素数の絶対値の入力例(上記のa,bの入力に続けて以下入力下さい)

 

↓虚数部の表現 : %i、 共役複素数 : conj()、 絶対値 : abs()

Complexcalc  

↓絶対値の計算確認ついでに平方根の入力確認 : sqrt()

Complexcalc2  

↓指数計算もできます。 : exp()

Complexcalc3  

ここまで分かれば周波数特性G(jω)の計算も楽々です。例えば、1次遅れ+無駄時間を確認してみます。

 

1次遅れ時間τ=10sec、無駄時間:DT=5secとして確認してみます。

 

1次遅れ+無駄時間の伝達関数は、

$$G(s)=\frac{K \cdot exp\left(-DT \cdot s\right)}{\tau \cdot s +1}$$ $$G(s)=\frac{exp\left(-5 \cdot s\right)}{10 \cdot s +1}$$

sにjwを入れると、

$$G(jw)=\frac{exp\left(-5 \cdot j \omega \right)}{10 \cdot j \omega +1}$$

ω=0.1, 1, 10, 100の時のG(jw)を見てみます。

ω=0.1, 1, 10, 100を行列に格納します。配列は以下のように簡単に宣言できます。

Complexcalc4

各ωに対するG(jω)は、for文で求めます。for文は、for a:b ~ endとします。

Complexcalc5

 

・伝達関数のプロット  

例として、以下の伝達関数について、①根軌跡、②Bode線図、③Nyquist線図、④Nichols線図をプロットする方法を示します。

$$f(s)=\frac{3\left(s-1\right)}{s(s+2)(s+5)}$$  

まず上記の伝達関数を定義しおきます。

s=poly(0,'s');
P = syslin('c',3*(s-1)/(s*(s+2)*(s+5)));
 

①根軌跡

evans(P);
Root  

②Bode線図

bode(P);
Bode  

③Nyquist線図

nyquist(P);
Nyquist  

④Nichols線図

black(P);
Nichols  

・伝達関数の時間応答

例として、以下の1次遅れプロセスについて、impulse応答、step応答をプロットする方法を示します。

$$f(s)=\frac{1}{s+2}$$

Impulse応答

s=poly(0,'s');
P=syslin('c',1/(s+2));
y=csim('impulse',t,P);
plot2d(t,y);
Time response imp

Step応答

s=poly(0,'s');
P=syslin('c',1/(s+2));
y=csim('step',t,P);
plot2d(t,y);
Time response step

時間応答を見るには、やはりMatalabの方が洗練されていて手軽です。(Matalabはstep(伝達関数)、impulse(伝達関数)でグラフ出力までできるので)

・外部ファイルの読み込み  

例えばcsvファイルに基本物性データを入力しておき、read関数で読みこむ場合、以下のように設定します。(パスは絶対パスでもOKです)

readfunction  

以下のようにcsvファイルが読み込まれます。

readresult  

例えばωのデータをOMEGAに代入したい場合、OMEGA=[Sheets(1,1) Sheets(1,2) Sheets(1,3)]のように、該当する行列番号を指定してやります。

readresulte

 

・sceファイルをサブルーチンとして呼び出す方法  

excec('ファイルパス')の形で呼び出しが可能です。

 

例)excec('vpcalc.sce')

 

例えば、蒸気圧推算は、成分によって式形が異なるので、メインプログラム側では記載せず、サブルーチン側で処理した方がコードが長くならずスッキリします。

 

状態方程式をScilabで解き、結果をグラフ描画する

状態方程式法の解法の詳細についてはこちらをご参照ください。

ここでは既に状態方程式法を3次方程式形まで持っていけることを前提に話を進めます。

Scilabでは3次方程式の解はroots()関数を使うことで簡単に求められます。

以下はSRKのzを求める場合の例になります。

y1=x.^3-x.^2+(Aa(1)-Bb(1)-Bb(1).^2)*x-Aa(1)*Bb(1)
roots(y1)

以下のように複素数解も求まります。

readresulte

グラフ表示させたい場合は、以下のようにplot2dで表示できます。

以下の例では、3通りの圧力、温度でz-f(z)のプロットを表示しています。

plot2d(x,[y1' y2' y3'],[1 2 3]);
xgrid()
xtitle('z diagram','z','f(z)')
hl=legend(['case1';'case2';'case3'],,3)
readresulte

Bode線図、Nyquist線図の作成

ScilabでのBode線図、Nyquist線図の作成方法についてはこちらをご参照ください

最適化関数の活用

Scilabのオンラインヘルプを参照しました。

Xcosの基本機能

1次遅れのステップ応答

まずは最小構成の1次遅れプロセスへのステップ応答です。

メニューバーのアプリケーション→XCOSを起動します。

パレットブラウザから以下のユニットを配置し、フローシートを作成してきいます。

なお、以下の構成は10秒間の応答を見るだけの構成になっておりますので、必要に応じてCSCOPEのRefresh periodや積分時間は見直してください。

グループユニット名設定1設定2
信号源Step_FunctionLステップ時間:0初期値:0、最終値:1
連続時間システムCLRNumerator:1Denominator:1+T*s
出力/表示CSCOPERefresh period:10-
信号源CLOCK_c間隔:0.1初期化時間:0.1

シミュレーション→設定→積分終了時間を1.0E01に変更します。

※Defaultでは1.0E5となっているので、このままでは単位ステップの応答は確認できません

例えば遅れ時間10secの場合、Scilabのコンソール画面でT=10とし、開始ボタン(再生マークです)を押して計算を開始します。

遅れ時間として設定した10secで、出力が1-1/eとなっていることが確認できます。

2次遅れのステップ応答

続いて2次遅れプロセスへのステップ応答です。

折角なので、1次遅れとの違いが見やすいように、入力信号と、1次遅れ、2次遅れの出力をまとめて一つのプロットで出してみます。(失敗過程も含めて以下記載します)

・1回目のトライアル(CMSCOPEを使ってみる)

xcosMul1

・2回目のトライアル、、、(CSCOPEの出力ウィンドウを揃えてみる)

xcosMul2

・3度目の正直!(MUXを使う)

xcosMul3

XcosからのBode線図出力

XcosのマニュアルにXcosのバッチランについてのコマンドがありましたのでご紹介します。

XcosをScilabからバッチ実行することで、Bode線図をXcosファイルから作成することができるようです。

詳細は、ブログにまとめたのでこちらをご参照ください。

以下のように、XcosからBode線図を作成することができます。

位相進み補償

微分方程式の解法

・ode関数

ode関数のタイプとして選択できるのは、以下の通り。

"adams" "stiff" "rk" "rkf" "fix" "discrete" "roots"

例題として、Pythonでも確認した硬い問題を試してみます。(Van Der Pol方程式)

比較前提

Time step:0.1[sec]、積分時間:2000[sec]、μ:1000

van der pol方程式は以下の通り。

$$\frac{d^2x}{dt^2}-\mu \left(1-x^2\right)\frac{dx}{dt}+x=0$$

初期値定義

X0 = [2.0; 1.0];

ode関数定義.以下のvanderpolが関数定義部になります。

[X] = ode("stiff",X0,0,T,vanderpol);

当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

Page Top