• 配布するファイル等の動作保証およびテクニカルサポート,発生した損害に対する責任は負いません.

準 備

 Windows 7/8/8.1/10 に MATLAB R2013a/2013b/R2014a/2014b/R2015a/R2015aSP1/2015b/R2016a がインストールされている環境を想定しています.
 配布する MATLAB/Simulink ファイルを実行するための準備として,以下の操作を行います.

  •  三つの zip ファイル をダウンロードします.
  •  これら zip ファイルを C ドライブのフォルダ
    • hoge ・・・C ドライブに各自で「新規作成」によりフォルダ hoge を生成する
    にコピーしてから解凍し,フォルダ
    • ip_toolbox_1.0.2
    • sedumi-master
    • yalmip
    を生成します.
  •  MATLAB のコマンドウィンドウで
    >> addpath(genpath('C:\hoge\ip_toolbox_1.0.2\iptools'))
    >> addpath(genpath('C:\hoge\ip_toolbox_1.0.2\odqlab_2.1.3'))
    >> addpath(genpath('C:\hoge\sedumi-master'))
    >> addpath(genpath('C:\hoge\yalmip'))
    と入力するか,もしくは,M ファイル を実行し,パスの設定を行います.パスが通っているかどうかはコマンドウィンドウで
    >> path

        MATLABPATH

      C:\hoge\yalmip
      C:\hoge\yalmip\demos
      C:\hoge\yalmip\extras
      C:\hoge\yalmip\modules
      C:\hoge\yalmip\modules\bilevel
      C:\hoge\yalmip\modules\global
      C:\hoge\yalmip\modules\moment
      C:\hoge\yalmip\modules\parametric
      C:\hoge\yalmip\modules\robust
      C:\hoge\yalmip\modules\sos
      C:\hoge\yalmip\operators
      C:\hoge\yalmip\solvers
      C:\hoge\sedumi-master
      C:\hoge\sedumi-master\conversion
      C:\hoge\sedumi-master\doc
      C:\hoge\sedumi-master\examples
      C:\hoge\sedumi-master\o_win
      C:\hoge\ip_toolbox_1.0.2\odqlab_2.1.3
      C:\hoge\ip_toolbox_1.0.2\iptools
      .............................. 略 ..............................
    のように実行することで確認できます.
     なお,管理者権限を持つユーザの場合は,
    • xxxxx/toolbox/local/matlabrc.m
         ・・・ xxxxx は MATLAB がインストールされているフォルダ
    の最下部に M ファイル の内容を追加することによって,MATLAB の起動時に自動的にパスを設定することもできます.


IP Toolbox

ファイルの説明

  舞鶴高専の川田研究室で製作した台車型倒立振子

およびアーム型倒立振子

をもとにして開発された MATLAB/Simulink 用シミュレータ

を提供します.Windows 7 (64bit) にインストールされた MATLAB R2013a/2013b/R2014a/2014b/R2015a/R2015aSP1/2015b/R2016a で動作を確認しています.
  ip_toolbox_1.0.2.zip を解凍すると,以下のファイル群が生成されます. 詳細については

をお読みください.

フォルダ名 ファイル名 説明
iptools ip_model.slx 倒立振子の非線形シミュレーション用の Simulink ライブラリ (画像)
cdip_para.m 台車型倒立振子の物理パラメータの定義 (p. 44 の表 3.5 (a), (c))
cdip_anime.m 台車型倒立振子のシミュレーション結果のアニメーション表示
cdip.jpg 台車型倒立振子の非線形シミュレータ用画像
cdip_photo.jpg 台車型倒立振子の非線形シミュレータ用画像
adip_para.m アーム型倒立振子の物理パラメータの定義 (p. 55 の表 3.7 (a), (c))
adip_anime.m アーム型倒立振子のシミュレーション結果のアニメーション表示
adip.jpg アーム型倒立振子の非線形シミュレータ用画像
adip_photo.jpg アーム型倒立振子の非線形シミュレータ用画像
odqlab_2.1.3 省略 動的量子化器を設計するために必要な MATLAB ツールボックス ODQ Toolbox/Lab (下記参照)
sample_cdip cdip_pdcont_sim.m
cdip_pdcont.slx
台車の P-D 制御の非線形シミュレーション
( Figure 1, Figure 2, Figure 3, アニメーション )
cdip_lqr_sim.m
cdip_lqr.slx
最適レギュレータによる台車型倒立振子の状態フィードバック制御の非線形シミュレーション
( Figure 1, Figure 2, Figure 3 )
cdip_lqr_servo_sim.m
cdip_lqr_servo.slx
最適レギュレータによる台車型倒立振子のサーボ制御の非線形シミュレーション
( Figure 1, Figure 2, Figure 3, アニメーション )
sample_adip adip_pdcont_sim.m
adip_pdcont.slx
アームの P-D 制御の非線形シミュレーション
( Figure 1, Figure 2, Figure 3, アニメーション )
adip_lqr_sim.m
adip_lqr.slx
最適レギュレータによるアーム型倒立振子の状態フィードバック制御の非線形シミュレーション
( Figure 1, Figure 2, Figure 3 )
adip_lqr_servo_sim.m
adip_lqr_servo.slx
最適レギュレータによるアーム型倒立振子のサーボ制御の非線形シミュレーション
( Figure 1, Figure 2, Figure 3, アニメーション )

使用例

 準備が終わったら,MATLAB のコマンドウィンドウで

>> ip_model Simulink モデル "ip_model.slx" の実行

と入力してください. このとき,非線形シミュレーションを行うための Simulink ブロック群 が表示されます.これら Simulink ブロックは

ブロック名 説明
Cart-driven Inverted Pendulum 台車型倒立振子の非線形モデル (3.22) 式 (p.43) \[ \left\{\begin{array}{rcl} \ddot{z}(t) \hspace{-0.3cm}&=&\hspace{-0.3cm} {}- {a}_{\rm c}\dot{z}(t) + {b}_{\rm c}v(t) \\ {m}_{\rm p}{l}_{\rm p}\cos\theta(t)\cdot\ddot{z}(t) + ({J}_{\rm p} + {m}_{\rm p}{l}_{\rm p}^{2})\ddot{\theta}(t) \hspace{-4cm} && \\ \hspace{-0.3cm}&=&\hspace{-0.3cm}{}- {\mu}_{\rm p}\dot{\theta}(t) + {m}_{\rm p}g{l}_{\rm p}\sin\theta(t) \end{array}\right. \] を表現
Cart-driven Inverted Pendulum (with Quantizer) 台車型倒立振子の非線形モデル (3.22) 式 (p.43) を表現 (D/A 変換の出力レンジ,分解能やロータリエンコーダの分解能を考慮)
Arm-driven Inverted Pendulum アーム型倒立振子の非線形モデル (3.58) 式 (p.54) \[ \left\{\begin{array}{rcl} \ddot{\theta}_{1}(t) \hspace{-0.3cm}&=&\hspace{-0.3cm} {}- {a}_{1}\dot{\theta}_{1}(t) + {b}_{1}v(t) \\ {\alpha}_{3}\cos{\theta}_{12}(t)\cdot\ddot{\theta}_{1}(t) + {\alpha}_{2}\ddot{\theta}_{2}(t) \hspace{-3.5cm}&& \\ \hspace{-0.3cm}&=&\hspace{-0.3cm} {\alpha}_{3}\dot{\theta}_{1}(t)^{2}\sin{\theta}_{12}(t) + {\alpha}_{5}\sin{\theta}_{2}(t) \\ \hspace{-0.3cm}&&\hspace{-0.3cm}{} + {\mu}_{2}\dot{\theta}_{1}(t) - {\mu}_{2}\dot{\theta}_{2}(t) \qquad \end{array}\right. \] を表現 (\({\theta}_{12}(t) = {\theta}_{1}(t) - {\theta}_{2}(t)\))
Arm-driven Inverted Pendulum (with Quantizer) アーム型倒立振子の非線形モデル (3.58) 式 (p.54) を表現 (D/A 変換の出力レンジ,分解能やロータリエンコーダの分解能を考慮)

を表現したものとなっています.
 Simulink ブロック

を利用して,Simulink モデル "cdip_pcont_sim.slx"

を作成します. コンフィギュレーションパラメータや Simulink ブロックのパラメータの設定については図 3.17 (p.58) を参照してください. Simulink モデル "cdip_pcont_sim.slx" は, 振子が真下にぶら下がって静止している状態を台車型倒立振子の初期状態とし, P コントローラにより台車を基準位置から目標位置に移動するための非線形シミュレーションを行います.

 以下で説明する操作は M ファイル "p1c34_cdip_plot_pcont.m" により一度に実行することができます (M ファイル "p1c34_cdip_plot_pcont.m" は Simulink モデル "cdip_pcont_sim.slx" と同じフォルダ内に入れてください).

 Simulink モデル "cdip_pcont_sim.slx" を実行するには, まず,台車型倒立振子の物理パラメータを表 3.5 (a), (c) (p. 44) の値に設定します.そのために,MATLAB のコマンドウィンドウで

>> cdip_para M ファイル "cdip_para.m" の実行

と入力します. つぎに,台車位置 \(z(t)\ {\rm [m]}\),台車速度 \(\dot{z}(t)\ {\rm [m/s]}\),(真上基準の) 振子角度 \(\theta(t)\ {\rm [rad]}\),振子角速度 \(\dot{\theta}(t)\ {\rm [rad/s]}\) の初期値を設定するため,

>> z_0 = 0;
>> theta_0 = pi;
>> dz_0 = 0;
>> dtheta_0 = 0;
\(z(0) = 0\ {\rm [m]}\) 
\(\theta(0) = \pi\ {\rm [rad]}\) (振子は真下で静止)
\(\dot{z}(0) = 0\ {\rm [m/s]}\) 
\(\dot{\theta}(0) = 0\ {\rm [rad/s]}\) 

と入力します.
 以上で非線形シミュレーションを行う準備が完了したので,

  • Simulink モデル "cdip_pcont_sim.slx" の実行ボタンをクリックする
  • Simulink モデル "cdip_pcont_sim.slx" のメニューから「シミュレーション/実行」を選択する
  • Simulink モデル "cdip_pcont_sim.slx" が入っているフォルダにカレントディレクトリを変更してから,MATLAB のコマンドウィンドウで
    >> sim('cdip_pcont_sim') Simulink モデル "cdip_pcont_sim.slx" の実行
    と入力する

のいずれかにより,Simulink モデル "cdip_pcont_sim.slx" を実行します.
 最後に,シミュレーション結果を描画するために

>> figure(1); plot(t,z)
>> figure(2); plot(t,phi)
Figure 1 に台車位置 \(z(t)\) を描画
Figure 2 に真下基準の振子角度 \(\phi(t) = \theta(t) - \pi\) を描画

と入力します.このとき,

のようにグラフが描画されます (このグラフは,M ファイル "p1c34_cdip_plot_pcont.m" の 18-26 行目のようにしてフォントなどをカスタマイズしています). さらに,

>> theta = phi + pi;
>> cdip_anime
振子角度を真上基準に変換 (\(\theta(t) = \phi(t) + \pi\))
M ファイル "cdip_anime.m" の実行

と入力することによって,以下のようにアニメーション表示されます.



SeDuMi/YALMIP

SeDuMi/YALMIP について


YALMIP の使用方法

 YALMIP を利用した LMI 求解の詳細は

およびそのサポートページ が参考になります.
 大雑把に説明すると,YALMIP を利用した LMI 求解のステップは以下のようになります.

ステップ 1  システム行列(状態方程式 \[ \dot{x}(t) = Ax(t) + Bu(t) \] における行列 \(A\), \(B\) など)や 重み行列(評価関数 \[ J = \int_{0}^{\infty}({x}(t)^{\top}Qx(t) + Ru(t)^2)dt \] における重み \(Q = {Q}^{\top} \succ 0\), \(R > 0\) など)を定義します.
ステップ 2
(決定変数の定義)
 関数 sdpvar を利用して,決定変数を定義します.
  • スカラーの決定変数 \(\gamma \in {\mathbb R}\) を定義するには,
    • gamma = sdpvar(1);
    • gamma = sdpvar;
    のいずれかを記述します.
  • 対称行列 \(X = {X}^{\top} \in {\mathbb R}^{n \times n}\) の決定変数を定義するには,
    • X = sdpvar(n,n,'symmetric');
    • X = sdpvar(n,n,'sy');
    • X = sdpvar(n,n);
    • X = sdpvar(n);
    のいずれかを記述します.
  • 対称とは限らない正方行列 \(Y \in {\mathbb R}^{n \times n}\) の決定変数を定義するには,
    • Y = sdpvar(n,n,'full');
    • Y = sdpvar(n,n,'f');
    のいずれかを記述します.
  • 長方行列 \(Z \in {\mathbb R}^{n \times m}\) の決定変数を定義するには,
    • Z = sdpvar(n,m,'full');
    • Z = sdpvar(n,m,'f');
    • Z = sdpvar(n,m);  ・・・ \(n \neq m\)
    のいずれかを記述します.
  • \({S}_{11} = {S}_{11}^{\top} \in {\mathbb R}^{n \times n}\), \({S}_{21} \in {\mathbb R}^{m \times n}\), \({S}_{22} \in {\mathbb R}^{m \times m}\) により構成される行列 \[ {S} = \left[\begin{array}{cc} {S}_{11} & 0 \\ {S}_{21} & {S}_{22} \end{array}\right] \] を定義するには,以下のように記述します.
    • S11 = sdpvar(n,n,'sy');
      S21 = sdpvar(m,n,'f');
      S22 = sdpvar(m,m,'f');
      S = [ S11\(\hspace{10pt}\)zeros(n,m)
      \(\hspace{22.5pt}\)S21\(\hspace{10pt}\)S22 ];
ステップ 3
(LMI の定義)
 [] を利用して,LMI を定義します.まず,LMI の定義を初期化するために
  • LMI = [];
と記述します.そのうえで,たとえば, 連立 LMI \[ \left\{\begin{array}{l} P \succ I\ \ \Longrightarrow\ \ P - I \succ 0 \\ PA + {A}^{\top}P \prec 0 \\ \end{array}\right. \tag{a} \] を定義するには,
  • ep = 1e-6;
    M1 = P - eye(n);
    M2 = P*A + A'*P;
    LMI = [LMI, M1 >= ep*eye(n)];
    LMI = [LMI, M2 <= - ep*eye(n)];
のように記述します.つまり,連立 LMI \({\rm (a)}\) の代わりに \[ \left\{\begin{array}{l} {M}_{1} := P - I \succeq {\varepsilon}I \succ 0 \\ {M}_{2} := PA + {A}^{\top}P \preceq {}-{\varepsilon}I \prec 0 \\ \end{array}\right. \tag{b} \] を定義します (\(\varepsilon\) は与えられた微小な正数). これは,最近のバージョンの YALMIP では
  • \({M}_{1} \succ 0\):\({M}_{1}\) が正定
  • \({M}_{2} \prec 0\):\({M}_{2}\) が負定
という記述をサポートしていないからであり,これらとほぼ等価な表現 (\(\varepsilon\) は与えられた微小な正数)
  • \({M}_{1} \succeq {\varepsilon}I \succ 0 \ \ \Longrightarrow\ \ {M}_{1} - {\varepsilon}I \succeq 0\):\({M}_{1} - {\varepsilon}I\) が半正定
  • \({M}_{2} \preceq {}-{\varepsilon}I \prec 0 \ \ \Longrightarrow\ \ {M}_{2} + {\varepsilon}I \preceq 0\):\({M}_{2} + {\varepsilon}I\) が半負定
で記述しています.
ステップ 4
(LMI の求解)
 関数 solvesdp を利用して,LMI の決定変数を求めます.
  • 定義された LMI に対して,デフォルトのソルバにより LMI 可解問題を解くには,
    • solvesdp(LMI)
    のように記述します.
  • 線形目的関数を \(E = \gamma\) (gamma) とした LMI 最適化問題を解くには,
    • solvesdp(LMI,gamma)
    のように記述します.
 SeDuMi がインストールされている場合, SeDuMi がデフォルトのソルバとなります. 問題が可解であった場合,
  • \(\hspace{19pt}\)info: 'Successfully solved (SeDuMi-1.3)'
    problem: 0
のように表示され,可解でなかった場合,
  • \(\hspace{19pt}\)info: 'Infeasible problem (SeDuMi-1.3)'
    problem: 1
と表示されます.
 問題によっては (たとえば,次数の高い問題を取り扱っている場合には),
  • \(\hspace{19pt}\)info: 'Numerical problems (SeDuMi-1.3)'
    problem: 4
と表示されることがあります. このときには,固有値を調べるなどして問題が解けているかどうか (\(M = {M}^{\top} \succ 0\) であることと \(M = {M}^{\top}\) の固有値がすべて正であることは等価です) を別途,吟味する必要があります.
 ソルバの種類,反復計算の最大回数などの設定を行いたい場合,関数 sdpsettings を使用します.たとえば,ソルバの種類を SeDuMi や LMILAB (MATLAB の Robust Control Toolbox に包含されているソルバ),反復計算の最大回数を 1,000 とするには,
  • ops = sdpsettings;
    ops.solver = 'sedumi';
    ops.sedumi.maxiter = 1000;
あるいは
  • ops = sdpsettings;
    ops.solver = 'lmilab';
    ops.lmilab.maxiter = 1000;
のように設定します.そして,LMI 可解問題の場合,
  • solvesdp(LMI,[],ops)
と記述し,LMI 最適化問題の場合,
  • solvesdp(LMI,gamma,ops)
と記述します.
ステップ 5  関数 double を利用して,得られた決定変数を倍精度に変換します. たとえば,関数 solvesdp により得られた決定変数 \(P\) は
  • P_feas = double(P)
のように入力することで,値を確認できます.

 \(\dot{x}(t) = Ax(t)\) の安定性を上述の \({\rm (b)}\) 式により判別する例は

を参照してください.



ODQ Toolbox/Lab

 ODQ Toolbox/Lab は 動的量子化器を設計するために必要な MATLAB ツールボックスであり,本書では

  • 第 II 部 (発展編) / 第 2 章 ディジタル制御 / 2.3 節 量子化入力制御

で利用します. ODQ Toolbox/Lab は 配布する IP Toolbox Ver.1.0.2 に同封されています.
 ODQ Toolbox/Lab に 関する情報は,開発者 (森田亮介,東 俊一,南 裕樹,杉江俊治) の Web ページ

もしくは文献

を参照してください.
 なお,本書ではソルバとして上述の SeDuMi を利用しています.