コンテンツへスキップ

短期集中:振動論と制御理論 工学系の数学入門 の感想です。

振動を例にモデル化しグラフを描くまで、制御の一番はじめの体験ができる本。

短期集中:振動論と制御理論 工学系の数学入門 日本評論社 吉田 勝俊

期間: ~2019年2月11日の3週間程度

内容(個人的まとめ)

  1. 力学の運動方程式
  2. 単振動,減衰振動,強制振動
  3. スケール変換
  4. ハーモニックバランス
  5. 複素数,ラプラス変換,伝達関数
  6. 解析力学入門
  7. 棒立て制御(倒立振子)
  8. 線形化

特徴

単振り子や倒立振子を取り上げて、ソフト「MATX」で計算を体験しながら理解できる本でした。全体的に簡単に書いているので、なんとなく理解できます。 キーワードが載っているので、次につなげることができます。

感想

(参考)読む前: 電気科3年後期のテスト期間でした。
力学は入試+学部1年でやったきりで、さっぱりわかりませんでした。 制御工学の初歩的知識は持っており、なんとなく伝達関数を記述できました。 また、私はpythonで計算しましたが、pythonはこれで初めてさわりました(最下部のリンク)。

全体

この本の内容は

うちの卒研生が配属直後の1カ月以内に本書の全課題をこなし、

本書 「はじめに」より

と書かれている通りの難易度です。

力学、振動の運動方程式の解き方など、非常に基本的な部分からスタートして倒立振子の姿勢制御までを紹介しています。 数値計算ソフトMATXを使用しながら演習できます。

私は MATX ではなくPythonでやりました。pythonは方程式を解いたりグラフを描いたりすることが簡単にできるのでおすすめです。MATXのコードが載っているので、コーディングの参考にできました。

力学がさっぱりわからない状態で、倒立振子の制御がやりたかった私に最適な本でした。参考にした本、この後に読むといい本が載っているのもかなり有難かったです。

1. 力学の運動方程式

よくある内容です。

2. 単振動,減衰振動,強制振動

電気科の私にはLRC回路が馴染みがありましたが、ばねとダンパーで復習できました。

3. スケール変換

2次方程式を解くのは面倒なので、パラメータを減らして1次方程式にします。このような操作をスケール変換と呼びます。 1次方程式にすると、プログラムで解くのが簡単になります。

4. ハーモニックバランス法

解き方です。

5. 複素数,ラプラス変換,伝達関数

基本的な内容です。

6. 解析力学入門

この本で、解析力学というものを初めて知りました。それだけでもこの本を読んで良かったと思いました。

解析力学では、運動の座標を変化させたり一般化して表現したりするみたいです。ラグランジュの運動方程式の知識があると、一般化された座標の運動方程式が比較的簡単に立てられます。

ラグランジュの運動方程式をこの本だけで理解するのは厳しかったです。

7. 棒立て制御(倒立振子)

ラグランジュの運動方程式は使わずに倒立振子の運動方程式を立てて、グラフに書きました(下記の関連記事)。途中式は(よくあることですが)書かれていないので大変でした。それっぽい波形が得られたので満足です。

8. 線形化

あまりよくわからなかったのでさらっと読みました。

次に読みたい

解析力学キャンパス・ゼミ 改訂2 マセマ出版社 馬場 敬之 (著)

力学 裳華房 原島 鮮 (著)

Pythonによる制御工学入門 オーム社 南 裕樹 (著)

関連記事

Jupyter Notebookで倒立振子の振動波形を表示する

Scilab/Xcosが必要になったので入れて使ってみます。動作の確認として、MATLABでやったのと同じ解析をしてみました。

Scilab/Xcosが必要になったので入れて使ってみます。windows10,Scilab6.0.2です。動作の確認として、MATLABでやったのと同じ解析をしてみました。

目次

  • 導入
    • scilabとは
    • Xcosとは
  • インストール
  • 使ってみる
    • ブロックを並べる
    • 設定
    • シミュレーション
    • 構文
      • polyとは
      • syslinとは
      • bodeとは
      • dampとは
  • MATLAB/SimlinkとScilab/Xocsの比較
    • MATLAB/Simlink
    • Scilab/Xcos

導入

Scilabとは

Scilabはフリーのオープンソフトウェアで、端的に言ってしまえばMATLABのフリー版です。様々な数値解析ができます。

詳しくは公式にて

MATLABに関してはこちらの記事をご覧ください。=>MATLAB使い方個人的まとめ

Xcos

動的な連続/離散時間システムのモデリングやシミュレーションができます。 MATLABにとってのSimlinkです。グラフィカルにブロックを並べることでモデリングでき、コンパイルやシミュレーションもまとめてできる優れものです。

X2C(https://www.lcm.at/virtuelle-entwicklung/x2c/)

インストール

公式ページ: Download Scilab 6.0.2 -scilab

言語設定でjapaneseを選択し、デフォルトのままfullバージョンでインストールします。

インストールページのスクショ

使ってみる

MATLABと同様に、コマンドウィンドウ、Xcosウィンドウ、SimNoteのウィンドウをそれぞれ使って解析を行います。アプリケーション/Xcosで起動するとXcosウィンドウとライブラリが現れます。

Scilabの画面
Xcos

ブロックを並べてブロック線図を作り、sin波を表示させてみます。

ブロックを並べる

まずは信号源、出力/表示などからsin波発生器、CSCOPE、CLOCKを見つけて右クリックかドラッグでXcos上に配置し、矢印をドラッグして配線します。 スコープは赤い矢印があるものを選びましょう。これはイベントシミュレーションを駆動するイベント信号になります。

ブロック図の作成

設定

それぞれのブロックで右クリックすると項目が現れて、その中からブロックパラメータを開きます。 CFSCOPEを設定します。[参考1]

  • Color or mark vector: スコープで表示するグラフ線の色またはマークの連続によるグラフを指定する。入力信号数にかかわらず8セットの数字が必要。数字が正のとき、連続線で色が指定でき、負のとき、マーク形状が指定できる。(plot2dコマンドに準ずる)
  • Output window number: グラフ・ウィンドウ番号の指定。
  • Output window position: グラフ・ウィンドウの表示位置指定。左上の位置をx,yのドットで示す。
  • Output window size: グラフ・ウィンドウの表示サイズをドット数で指定。
  • Ymin: グラフ・ウィンドウのY軸最小値
  • Ymax: グラフ・ウィンドウのY軸最大値
  • Refresh period: グラフ・ウィンドウのX軸(時間)最大値 シミュレーション設定時間がこの値より大きい場合、グラフをクリアし再度続きの時間から描画する。基本的使い方では、シミュレーション時間と一致させる。
  • Buffer size: 入力データを一時的に蓄えておくメモリ領域のサイズを指定。グラフ表示速度に大きく影響するので、変化を観察したい場合以外は、少なくともRefresh period以上に設定する。最小は2。
  • Accept herited events: (inherited?) 継承の無効(0)・有効(1)を指定。1にするとイベント・ポート赤が削除される。アクティベーションが上流で行われる場合、継承が可能。

SCOPEの積算時間を10に、最小値~最大値を-10~10にしました。 クロックとsin発生器のパラメータも設定します。初期化時間を0.0に、ゲイン(振幅)5にします。

シミュレーションの設定

シミュレーション

シミュレーションを開始します。シミュレーション/パラメータ設定から、[参考1]

  • Final integration time: シミュレーション終了時間[単位は任意]
  • Real time scaling: 正の実数nを入力すると、シミュレーション設定単位時間tがn×t秒としてシミュレーションが実行される。n=1とした場合、シミュレーションの設定単位時間が1秒として解析表示される。
  • Integrator absolute tolerance: ソルバーの絶対誤差公差
  • Integrator relative tolerance: ソルバーの相対誤差公差
  • Tolerance on time: ODE(常微分方程式)ソルバーの最小時間公差
  • Max integration time interval: ソルバー呼び出しの最大時間間隔。"too many calls"というエラーメッセージが出た場合、値を小さくする。(解析精度を落とす。)
  • Solver 0(CVODE) - 100(IDA): 数値ソルバーの選択。
  • maximum step size: ソルバーがとる最大時間ステップ

最終積算時間を1.0E1.0、すなわち10sにしました。 AEBならA*10^Bになります。

document/test_0620/image.scgはこのようになりました。指定した通りに動いていることが確認できます。

document/ test_0620/image.scg の出力

構文

ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 3. 付録: Scilab言語の基本構文 -SimCircuit

Scilab コンソールの使い方の基本

「;」を後ろにつけると結果を表示しないみたいですね。

polyとは

[参考2]によれば、ploy関数は以下のように書きます。

p=poly(実数もしくは行列,"変数")
-->a=[-3 -4];
-->p=poly(a,"x")
 p  = 

            2
   12 +7x +x 

コンソールに上のように入力すると、12+7x+x^2を返します。((x+3)(x+4)=0, のときx=-3,-4) pを多項式変数として扱う事ができます。

-->s=poly(0,'s')
s =

    s
polyをコンソールに入力したときの出力

syslinとは

[参考3]によれば、syslin関数は以下のように書きます。

q=syslin(連続時間システムc or離散時間システムd,多項式)
--> s=poly(0,'s');

--> P=syslin('c',1/(s+1))
 P  = 

           
     1     
   ------  
           
   1 + s   

bodeとは

[参考4]によれば、syslin関数は以下のように書きます。

bode(syslinで定義された線形システムの多項式,x軸の最小値,x軸の最大値,x軸の刻み,['コメント'])

コメントにタイトルなどをかきます。 こだわらないなら関数のみ与えればOKです。

dampとは

[参考5]によれば、syslin関数は以下のように書きます。

[wn,z] = damp(連続時間のsys関数)

MATLAB/SimlinkとScilab/Xcosの比較

某授業のレポート1でやった問題の伝達関数を使って確認していきます。

問題: オープンループの伝達関数L(s)が以下の時の、回ループと閉ループでのボーデ線図をそれぞれ表示します。

$$L(s)=\frac{0.1}{s(s+0.5)(s+1)}$$

MATLAB

  • オープンループ

$$L(s)=\frac{0.1}{s^3+1.5s^2+0.5s}$$

MATLABオープンループでのボーデ線図

位相余裕:59° 利得余裕:17.5dB

  • フィードバック $$L(s)=\frac{0.1}{s^3 + 1.5s^2 + 0.5s + 0.1}$$

A. 周波数応答

MATLABフィードバックでのボーデ線図

B. ステップ応答

設定がうまくいかず、ぎこちないですがこんな感じ。最終値は1です。

MATLABフィードバックでのステップ応答

Scilab

  • オープンループ

複数行のコマンドを何度も使用する場合、アプリケーション/SciNoteに書くとスクリプトのようにまとめて実行できます。 見やすさのために出力なしで実行しました。

SciNoteを使用している様子

A. 周波数応答

まずはコマンドで実行してみます。

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
bode(p)
show_margins(p)
[gm,fp]=g_margin(p) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(p) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 

出力

ゲイン余裕[dB]

   17.501225

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   59.289834

 ゲイン交点角周波数[Hz]

   0.0293667

ゲイン余裕:位相が-180°と交差する周波数(位相交差角周波数)においてゲインが0dBより何dB小さいか

位相余裕:ゲインが0dBと交差する周波数(ゲイン交差角周波数)

オープンループボーデ線図は

Scilab オープンループでのボーデ線図

MATLABと一致してますね。

B. ステップ応答

Scilab オープンループでのステップ応答
  • フィードバック

A. 周波数応答

コマンドからやります。フィードバックした場合の伝達関数G(s)は以下のようになるので

$$G(s)=\frac{L(s)}{1+L(s)}$$

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
g=p/(1+p)
disp(g)
bode(g)
show_margins(g)
[gm,fp]=g_margin(g) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(g) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 
[wn,z] = damp(g); disp("固有周波数[Hz]") disp(wn) disp("減衰係数") disp(z) 

出力

                           
            0.1            
   ----------------------  
                    2   3  
   0.1 + 0.5s + 1.5s + s   

 ゲイン余裕[dB]

   16.258267

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   114.71226

 ゲイン交差角周波数[Hz]

   0.0313482

 固有周波数[Hz]

   0.2964606
   0.2964606
   1.1378001

 減衰係数

   0.6108736
   0.6108736
   1.

同じ曲線を描いていることを確認しました。

Scilab フィードバックでのボーデ線図

B. ステップ応答

最終値1で、100秒まで計算しています。

Scilab フィードバックでのステップ応答

コードは[参考6]を参考にさせて頂きました。

雑記

  • Xcosのライブラリで固まったらEscキーを押しましょう!ドラッグより右クリックの方がいいです。
  • 周波数応答を調べるとき、関数に伝達関数を与えています。なのでXcosでブロック線図を描こうと思ったらブロックを並べて伝達関数を計算させて、それを取得して関数に入れる、という手順を踏むと思うので回りくどいです。ここでやった例はすぐ伝達関数を計算できるので、コマンドでやった方が早いと思いました。
  • 何気に変数ブラウザ便利ですね。

参考

1.ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 2.3 Xcosによるシミュレーション・ツールとしての使い方 (GUI型ブロック線図入力モジュール) -SimCircuit

2. Scilab関数の説明:polyとは

3. Scilab関数の説明:syslinとは

4. Scilabでのボード線図の描き方

5. damp -scilab

6. scilabで制御系設計

MATLABの使い方を個人的にメモします。 ボード線図と時間応答。ごく簡単な部分だけです。

本家様

MATLAB 入門 -MathWorks

目次

  • MATLABとは
  • simlinkとは
  • 構成
  • 使い方

MATLABとは

技術計算のための有料ソフトで、基本的な機能の他はさらに課金する必要があります。その拡張機能によって様々な場面で使用されており、企業で使用されるケースも多々あります。学生は研究室で使えたりします。個人で購入する場合でも学割があり、生協でかなり安く買えるそうです。

行列と配列の数学を直接表現するプログラム言語と反復的に実施する分析、設計プロセスに適したデスクトップ環境がセットになっています。

基本的には3つのウインドウをつかう。Simlinkは拡張機能のようなもの。

  • コマンドウィンドウ
  • m-fileウィンドウ
  • Simlinkウィンドウ

MATLAB -MathWorks

(追記)無料の似たようなソフトScilabというものを見つけました。=>Scilab/Xcosのインストールと使い方 MATLABとの比較

simlinkとは

Simlinkウインドウを使用すると、グラフィカルにブロック図を書けばコードに変換してくれます。制御工学とかと相性が良いですね。

Simlink -MathWorks

構成

コマンドウィンドウ

会話的で、すぐに結果が出ます。

>> a = [1 2 3; 4 5 6; 7 8 9]
a =
        1 2 3
        4 5 6
        6 7 8
コマンドウィンドウ

こんな感じ。関数も使えます。範囲で指定する方法がわからなかったので直接配列を書きました。

>> x = [0 pi*1/6 pi*1/3 pi*2/3 pi*5/6 pi pi*4/3 pi*3/2 pi*5/3 2*pi]
>> plot(x,sin(x))
sin(x)

linespce(始点、終点、点数)を使えば楽ですが、研究室に合ったやつには入っていない関数らしい。

関数を定義する事も出来ますが、とくに使う予定がないので省略します。

m-fikeウィンドウ

shellみたいなものだと思います。上とほとんど同じです。

>> x = [0 pi*1/6 pi*1/3 pi*2/3 pi*5/6 pi pi*4/3 pi*3/2 pi*5/3 2*pi]
>> plot(x,cos(x))

保存するとき、名前がモデルと完全に同じだとエラーを吐くらしいので、微妙に変えておくと無難です(user/tutorial01_p.mにしました)。メニューバーのデバッグ/実行で実行します。この時たぶん、カレントディレクトリにありません、などエラーを吐いたら、コマンドウィンドウ上の方のカレントディレクトリを変更します。

デバック/実行か実行アイコンから実行するとこんな感じ。

cos(x)

Simlinkウィンドウ

ドラッグアンドドロップでブロック線図を作成します。割愛。

ちなみに、m-fileのエディタでmdlファイルを開くとめちゃくちゃいろいろ書いてます。

使い方

簡単な周波数解析

MATLABを起動し、メニューバーにあるSimlinkアイコン(ごちゃついたヤツ)からSimlinkを起動します。

Simlink

Simlinkを起動すると、自動的にlibrary Browserが開くと思います。

Simlinkのファイル/新規作成/モデルをクリックして、出てきたウィンドウを保存(user/tutorial01.mdlにしました)。

Simlink library Browser には関数やinputなどが入っていて、libraryからmodelへドラッグアンドドロップすることでブロック線図を完成させます。

ひとまずこのような図を作ります。

  • Continuousには1次遅れ系の伝達関数っぽいブロック(Transfer Fcn)などが入っています。
  • Sinksにはoutput(Out1)やオシロスコープ的なもの(Scope)などが入っています。
  • SoucesにはInput(In1)や時間測定器(Clock)
continuous
sink
souces

ブロックの端に「>」というのがありますので、それをドラッグして他のブロックに接続します。

開ループ

これでひとまず開ループの完成。矢印上で右クリック・ドラッグすると→を分岐できます。 Math Operatioのsumと合わせて使用して閉ループを作ることもできます。

関数をバブルクリックするとパラメータを変えるウィンドウが出てきますので、適宜変えていきます。

閉ループ

Numeratorを[1]、Denominator[2 3] にすると、

$ \frac{1}{2s+3}$

$$ \frac{1}{2s+3}$$

Numeratorを[1]、Denominator[2 -3 4] にすると、

$$ \frac{1}{2s^2-3s+4}$$

になります。ブロック図ができたらプロットのためのプログラムを組みます。

m-fileを作成します。

コマンドウィンドウから、新規作成/M-fineで、M-fileウィンドウを開きます。

T = 1
[num, den]= linemod('tutorial01') 
w=logspace(-3,3,100); 
figure 
bode(tf(num,den), w) 
grid on 

これで、実行すると、ボード線図が描画されます。このfigureで保存すると、画像を保存できます。

bode線図

簡単な時間解析

時間情報をグラフに書き込むときは、出力情報に紐づけられた時間の配列を変数に格納してplot(t,y)などを実行することになります。

同様にモデルを作成します。

  • Sourceからステップ信号(Step) 、時間(Clock)
  • Continuousから伝達関数(Transfer Fcn)
  • Sinkから「外部へ」(simout)

ブロック線図を作成したら、simoutをダブルクリックして、Clockとつながっている方は変数名をt、保存フォーマットをArray(配列)に、transfer Fcnとつながっている方は変数名をy、保存フォーマットをArray(配列)に変更します。

シミュレーションパラメータ

メニューバーのシミュレーション/シミュレーションパラメータから開始時間を0、終了時間を10*T、ソルバオプションのタイプを固定ステップ、ode1(Euler) ステップサイズをT/10に変更したら、メニューバーの下の△で実行。

コマンドウィンドウで、

>>figure
>>pot(t,y)
>>xlabel('Time[s]'),ylabel('Output')

とすると、こんな感じ。

時間解析

これを、m-fileのコマンドでやると以下のようになります。

T = 1;
figure
plot(t,y)
xlabel('Time[s]'),ylabel('Output')
grid on

以上です。あとはその都度組み合わせて使っていきましょう。

(追記 2019/6/24)Scilabに関して追記しました。