MuPAD入門書 - JMAAB（Japan MATLAB Automotive Advisory Board

プラントモデリングのための MuPAD® 入門
【 本書作成の狙い ・・・本文 2.4 章（p.7）から抜粋・・・ 】
自動車制御系開発の現場では、数式処理ツールはまだほとんど利用されていない。そのため、いくら数式処理がプラン
トモデリングのあるべき姿の実現に有効と言っても、実際に数式処理ツールを使おうとすると初歩的な使い方から学習が
必要となり、効果的に使えるようになるまでには時間を要してしまう。
また、既存の数式処理ツールの解説記事には、制御系開発のためのプラントモデリングという立場で書かれたものは存
在していないと思われる。そこで今回、誮でも自由に入手・再配布が可能な入門書を作成し、公開することとした。
JMAAB Plant Model Workshop 3
1
目次
1
入門書の取り扱い ................................................................................................................................... 5
1.1 著作権 ................................................................................................................................................................. 5
1.2 本書の取り扱い .................................................................................................................................................... 5
1.3 内容に関する問い合わせ先 ................................................................................................................................ 5
1.4 注意事項.............................................................................................................................................................. 5
1.5 本書において必要となる MathWorks 製品 ....................................................................................................... 5
2
はじめに ................................................................................................................................................... 6
2.1
JMAAB および Plant Model Workshop 3 （PMWS3） について......................................................................... 6
2.2
プラントモデリングのあるべき姿 ............................................................................................................................ 6
2.3
プラントモデリングに係る MathWorks 製品について .......................................................................................... 7
2.3.1 Simulink® と SimscapeTM の関係 ................................................................................................................ 7
2.3.2 Symbolic Math ToolboxTM と MuPAD® の関係 .......................................................................................... 7
2.4 本書作成の狙い .................................................................................................................................................. 7
2.5
3
プラントモデリングと本書のカバー範囲 ................................................................................................................ 8
MuPAD® の使い方 ............................................................................................................................... 10
3.1 MuPAD® の起動................................................................................................................................................ 10
3.2 数式処理............................................................................................................................................................ 13
3.2.1 例題:
3.2.1.1
回転摩擦トルクモデル........................................................................................................................ 13
回転摩擦トルクモデルの数式処理と解説 ............................................................................................... 13
3.2.2 例題: DC モーターモデル............................................................................................................................ 22
3.2.2.1 DC モーターモデルの数式処理と解説 ................................................................................................... 22
Appendix ..................................................................................................................................................... 32
A.1
MATLAB® と MuPAD® の関係と違い .............................................................................................................. 33
A.2
MuPAD® ヘルプの使い方 ................................................................................................................................. 35
A.2.1 ヘルプウィンドウの起動 ............................................................................................................................... 35
A.2.2 ヘルプウィンドウの使い方 ............................................................................................................................ 36
A.3
Notebook の実行方法...................................................................................................................................... 40
A.4
代入 （:=） と 等号 （=） の考え方 .................................................................................................................. 41
JMAAB Plant Model Workshop 3
2
A.5
piecewise と if … then … の違い...................................................................................................................... 43
A.6
数式処理による状態変数の取得 ...................................................................................................................... 47
A.7
MuPAD® mu ファイルについて .......................................................................................................................... 49
参考情報..................................................................................................................................................................... 50
JMAAB Plant Model Workshop 3
3
Revision
Rev
1
JMAAB
内容
Plant Model Workshop 3
日時
2013 – 06
JMAAB Plant Model Workshop 3
備考
4
1
入門書の取り扱い
1.1 著作権

本書の著作権は JMAAB （Japan MBD Automotive Advisory Board） に帰属します。

JMAAB は、本書の内容に関し、いかなる保証もするものではありません。万一本書を利用して丌具合等があった場合
でも、JMAAB は一切責任を負いかねます。また、本書に記載されている事項は予告なしに変更または廃止されること
がありますので、あらかじめご了承ください。
1.2 本書の取り扱い

本書は、非営利目的、または利用者内部で使用する場合に限り、複製が可能です。

また、本書を引用する場合は、本書からの引用であることを明示し、引用された著作物の題号や著作者名を明示する
等の引用の要件を満たす必要があります。
1.3 内容に関する問い合わせ先
本書に関するお問い合わせは、JMAAB 事務局（[email protected]）までご連絡下さい。
1.4 注意事項
本書におけるツールのバージョンは、R2012b を使用しています。他のバージョンでは、操作環境、操作結果、表示などが
異なることがあります。
1.5 本書において必要となる MathWorks 製品

MATLAB®

Symbolic Math Toolbox

SimscapeTM
MathWorks 社における基幹アプリケーション。
TM
MATLAB のオプション製品。本製品に MuPAD® が含まれる。
プラントモデリングに特化した Simulink のオプション製品。
JMAAB Plant Model Workshop 3
5
2
はじめに
2.1 JMAAB および Plant Model Workshop 3 （PMWS3） について
JMAAB （Japan MBD Automotive Advisory Board） は、「開発環境では協調し、競争は製品で」を基本理念として 2001 年
に日本のいくつかの自動車メーカと ECU サプライヤによって MathWorks 製品のユーザー会として発足した。日本の自動
車業界におけるモデルベース開発技術の普及・推進を加速すべく、2011 年に JMAAB の M を MATLAB の M から
Model-Based Development の M へ変更した。2013 年 3 月現在、自動車メーカ 12 社、 ECU サプライヤ 12 社の合計 24
社がコアメンバーとして活動している。
プラントモデル（制御対象モデル）に関わる JMAAB での取組みとしては、モデルの効率的な再利用方法や、組織間でモ
デルを容易に交換するために必要な要件の明確化、モデリングルールの標準化を目指した活動などを行っている。
2012 年 6 月から 2013 年 3 月までの間、Plant Model Workshop 3 （PMWS3）として数式処理ツール MuPAD® （Symbolic
Math ToolboxTM に含まれる） とプラントモデリングツール SimscapeTM を対象にしたワークフローの検討などに取組み、本書
もその一環として作成された。
PMWS3 メンバー会社を以下に示す。
【PMWS3 メンバー会社】
カルソニックカンセイ株式会社, 株式会社デンソー, トヨタ自動車株式会社, トヨタテクニカルディベロップメ
ント株式会社, 日立オートモティブシステムズ株式会社, 株式会社本田技術研究所, マツダ株式会社, 三菱ふそ
うトラック・バス株式会社
2.2 プラントモデリングのあるべき姿
大規模化と複雑化が進んでいる自動車制御系の開発において、制御対象が未完成の状態であっても制御系の開発を進
めるフロントローディング開発が重要性を増している。ハード未完の状態で制御系の開発を進めるには、ハードの代わりとなる
モデル（プラントモデル）を必要に応じてタイムリーに作成または入手して利用することが必要である。
しかしながら、制御系に対する要求・制約の高度化に伴ってプラントモデルに対する要求・制約も高度化しており、制御開
発プロジェクトでタイムリーにプラントモデルを利用していくには、既存のモデルを繰り返し使うだけでは必ずしも十分ではない
場合が増えている。
そのような場面では、新しい要求に応えられる新しいプラントモデルを納期に遅れることなく作らなければならない。また、
一種類の自動車の開発にも多数の会社が関わっており、会社の壁を越えたプラントモデルの流通によって制御系の開発を
促進することも重要である。
従って、プラントモデリングのあるべき姿は、単にモデルライブラリのようなものを充実させれば良いといった単純な話ではな
く、モデルを短期間で作成したり修正したりする技術、モデルを流通させる仕組み、システム全体のシミュレーションをコントロ
ーラーとの閉ループで実施する技術など、様々な側面を持っている。
JMAAB Plant Model Workshop 3
6
2.3 プラントモデリングに係る MathWorks 製品について
2.3.1
Simulink® と SimscapeTM の関係
MATLAB 製品群のうち、プラントモデルの作成には主に、 Simulink または Simscape が用いられる。Simscape はプラント
モデリング専用として 2007 年にリリースされた比較的新しい Simulink のオプション製品であり、 Simscape がリリースされる
までは、一般的なモデリング環境である Simulink が制御モデルの作成やプラントモデリングに用いられてきた。しかしながら、
Simulink のソルバーは基本的に常微分方程式を対象とするため、一般的に、微分代数方程式で表現される物理システム
をモデル化するには、モデル作成者があらかじめシステムを大幅に単純化しなければならない場合がある。モデル化の対象
が比較的小規模なシステムや、大規模であっても物理的には単純なシステムの場合には問題は無いかもしれないが、大規
模かつ複雑な物理システムを Simulink によりモデル化するのは困難である。
一方、 Simscape は微分代数方程式ソルバーを備え、Simscape のコンポーネント内部には微分代数方程式を直接的に
記述することができる。また、コンポーネント間の接続には保存則や釣り合いの式が適用される仕組みが搭載されており、大
規模で複雑な物理システムのモデル化にも対応可能である。更に、Simscape のコンポーネントと Simulink のブロックをフレ
キシブルに接続する仕組みも提供されており、コントローラーとの閉ループ構築も容易に可能である。従って、先述したプラン
トモデリングのあるべき姿の観点からは、プラントモデリングには Simulink よりも Simscape を用いることが望ましい。
ただし、Simscape はあくまでも微分代数方程式系の数値シミュレーション専用ツールであり、Simscape のコンポーネント
内部に記述する式の物理的な正しさを保障したり、既に記述されている式の物理的な意味を理解したりすることは、ユーザ
ーの責任である。これらの用途には、微分積分や代入消去、更にグラフ化など、数式を対話的に処理することが必要なため、
数式処理専用ツール MuPAD® を利用するのが便利である。
2.3.2
Symbolic Math ToolboxTM と MuPAD® の関係
MATLAB® では数式処理を行う製品として Symbolic Math Toolbox (SMT) を提供している。
SMT には 2 つの利用方法がある。1 つは MATLAB のプロンプトや M スクリプトで数式処理コマンドとして利用する方法で
ある。もう 1 つは MuPAD と呼ばれる独立した数式処理専用 GUI アプリケーションとして利用する方法である。いずれの場
合も数式処理を実際に実行するエンジンそのものは同じである。
物理システムをプラントモデルとして数式で表現する場合、時間のみを独立変数とする集中定数系として表現することが
多い。MuPAD では、任意の物理量を時間の関数として定義し、数式処理や数値計算を適用することが可能である。
MATLAB と MuPAD の関係と違いについては Appendix A.1 に記載した。
2.4 本書作成の狙い
自動車制御系開発の現場では、数式処理ツールはまだほとんど利用されていない。そのため、いくら数式処理がプラント
モデリングのあるべき姿の実現に有効と言っても、実際に数式処理ツールを使おうとすると初歩的な使い方から学習が必要
となり、効果的に使えるようになるまでには時間を要してしまう。
また、既存の数式処理ツールの解説記事には、制御系開発のためのプラントモデリングという立場で書かれたものは存在
していないと思われる。そこで今回、誮でも自由に入手・再配布が可能な入門書を作成し、公開することとした。
JMAAB Plant Model Workshop 3
7
2.5 プラントモデリングと本書のカバー範囲
プラントモデルには、大きく分けて 2 種類の形態がある。1 つは部品単位でのモデル（コンポーネントモデル）であり、もう 1 つ
はコンポーネントモデルを組み合わせて構成されるモデルである。
ただし、モデル化する対象を、1 つのコンポーネントとしてモデル化するべきなのか、或いは、複数のコンポーネントを用いて
モデル化するべきなのかは実際の制御系開発の要求に応じて決められるべきものであり、モデル化対象だけで決めることは
できない。また、モデル化の際に必要となる技術という観点からは、1 つのコンポーネントモデルの作成と、複数のコンポーネ
ントから成るモデルの作成において、いずれの目的にも共通のツールが有効に使えるとは限らない。例えば、1 つのコンポー
ネントモデルを作成する場面では、基本的な物理法則の式から定式化を始めて最終的にシミュレーション可能な式を導出
することになるが、その時には数式処理技術を対話的に使うことが想定される。これには先述の通り MuPAD が適している。
一方、複数のコンポーネントを組み合わせてシステムを構築するときには、コンポーネント間の接続に対して物理量の釣り合
いの式や保存の式を自動的に生成する機能が備わっている SimscapeTM を用いるのがふさわしい。
本入門書のカバー範囲としては、1 つのコンポーネントモデルを作成する場面を想定しており、そこで用いる数式処理
（MuPAD®） の基礎について例題を用いて解説を行っている（下図 A 参照）。
物理システムのダイナミクスをプラントモデルとして数式で表現し MuPAD で取り扱うにあたって、本入門書では以下の 3
つのケースを考えることにする。
Case 1.
解析解
Case 2.
定常状態
Case 3.
数値積分（シミュレーション）
1 つ目の解析解は、モデルによって導出可能な場合とそうでない場合がある。解析解を得ることができる場合には、シミュ
レーション結果を確かめるための比較対象として利用する。
モデルが複雑で解析解が得られない場合でも、何らかの定常入力に対してシステムが到達する定常状態ならば、求める
ことができるはずである。そこで、 MuPAD で定常状態を求めることを 2 つ目のケースとする。下図 A・B において、これら
Case 1 と Case 2 は 「① 数式的アプローチ」 に相当する。
実開発で用いるプラントモデルは通常、解析解の導出や定常状態の分析だけで十分に検討しきれるほど単純ではなく、
最終的にはシミュレーションをしなければならないことが多い。そこで 3 つ目のケースとして、 MuPAD でシミュレーションを実
施する場合についても取り上げる。これは図中の 「② 数値的アプローチ」 に相当する。
解析解や定常状態をあらかじめ調べた上でシミュレーションを行えば、シミュレーションの結果が妥当であるかを確実に判
断することができるので、本入門書ではそのような流れを想定した MuPAD の使い方を紹介する。
尚、非線形なモデルに対して動作点を仮定して線形化し制御設計を実施するなどといった内容は本入門書では取り扱わ
ない。
JMAAB Plant Model Workshop 3
8
図 A： プラントモデリングと本書のカバー範囲
図 B： プラントモデリングのワークフロー
JMAAB Plant Model Workshop 3
9
3
MuPAD® の使い方
ここでは、MuPAD の起動方法、及び、 MuPAD の基本的な使い方について、回転摩擦トルクと DC モーターを具体例と
して取り上げながら説明する。
3.1 MuPAD® の起動
MuPAD は次の方法で起動できる。
【コマンドによる起動】
MATLAB コマンドウィンドウで mupad とタイプイン
【マウス操作による起動】
□ R2012a 以前 □
1) MATLAB 起動画面のスタートボタンから Symbolic Math  MuPAD を選択
（スタートボタンは R2012b 以降廃止）
JMAAB Plant Model Workshop 3
10
2) Welcome to MuPAD ウィンドウ から New Notebook を選択
（Welcome to MuPAD ウィンドウ は R2012b 以降のバージョンにおいても mupadwelcome コマンドで起動する）
□ R2012b 以降 □
1) MATLAB 起動画面のアプリケーションタブを選択しプルダウンをクリック
2) MuPAD Notebook を選択
JMAAB Plant Model Workshop 3
11
上記の操作より、Notebook と呼ばれる次の図のような専用エディターが起動する。
Notebook は以下の 3 つの入力領域から構成される。
テキスト入力領域
コマンド入力領域
実行結果出力領域
グラフ表示領域
JMAAB Plant Model Workshop 3
12
3.2 数式処理
MuPAD による数式処理について、実システムにおける物理現象を用いて解説する。
ここでは MuPAD で検討した数式を Simscape 言語として出力し、 Simulink による制御設計において活用することを念
頭においている。したがって、以下 2 項では、回転摩擦トルクモデルと DC モーターモデルを例題とし、 Simscape での活用
を前提とした数式処理、及び、操作手順について説明している。
3.2.1
例題: 回転摩擦トルクモデル
回転摩擦モデルは 2 つの並置された回転体が接触することにより発生する摩擦トルクを算出するモデルである。
回転摩擦モデルの概念図を以下（左図）に示す。尚、 Simscape の Foundation Library には、回転摩擦を模擬する
Rotational Friction ブロック（コンポーネント）が用意されている（右図）。以下では、この Rotational Friction ブロックのソース
コードに記述されているモデル式と同じものを MuPAD で分析する。
回転摩擦トルクモデルの概念図 及び Rotational Friction ブロック
3.2.1.1
回転摩擦トルクモデルの数式処理と解説
回転摩擦トルクモデルの数式は、並置された 2 つの回転体の間に発生する摩擦トルクを、回転体間の角速度差の代数式
として表現している。回転摩擦トルクモデル（代数式）の数式処理について以下に示す。また、本例題の流れは、2.5 節で述
べたプラントモデリングのアプローチ「① 数式的アプローチ」の Case 1 に該当する。
尚、ここでは、各コマンドの詳細な説明は省略し、ヘルプの記載箇所を示した。
となる。ヘルプウィンドウの詳細な使い方は Appendix
は、ヘルプウィンドウ起動アイコンの意
A.2 に記載している。また、以下の段落において、左の段落は
MuPAD Notebook の各入力領域をイメージしている。記述されるコマンドを、実際の MuPAD Notebook にコピー＆ペースト
して実行できるよう考慮した。右の段落は MuPAD コマンドに対する解説である。
Rotational Friction
1. Initialize MuPAD session
1. MuPAD の初期化
reset():
Notebook に対する初期化は reset コマンドを用
いる。このコマンドは開いている Notebook 内の数式
や変数を初期化する。
Notebook の初めに reset コマンドを実行するこ
とは必頇ではないが、この方法を推奨する。
JMAAB Plant Model Workshop 3
13
：（コロン）はステートメントの行末を示し、実行結果
が非表示となる。実行結果を表示させる場合は、;
（セミコロン）を使用する。
Notebook では、コマンド入力領域に記述した各
ステートメントは Enter キーを押すことで実行できる。
Enter
キー以外の実行方法もあるが、詳細は
Appendix
A.3 に記載している。
尚、コマンド入 力領域に おけ る 改行は Shift +
Enter を用いる。
補足として、Ctrl + i でコマンド入力領域を、また、
Ctrl + t でテキスト入力領域を挿入することができ
る。
sysname(Arch);
Notebook の流通や再利用を考慮し、ここでは、
Notebook の作成環境情報を取得する。
sysname コマンドは、OS 名を出力する。
version();
version コマンドは、MuPAD カーネルバージョンを
出力する。
【R2012a 以前】
> Commands Tab > “reset” とタイプイン
> Commands Tab > “ sysname ” とタイプイン
> Commands Tab > “version” とタイプイン
> Commands Tab > “ : ” とタイプイン
> Commands Tab > “ ; ” とタイプイン
【R2012b 以降】
> MuPAD Functions > reset を選択
> MuPAD Functions > sysname を選択
> MuPAD Functions > version を選択
> MuPAD Functions > ; を選択
(表示内容に : についても記載されている)
JMAAB Plant Model Workshop 3
14
2. Define parameter symbols
2. ギリシャ文字の定義
w := Symbol::omega; // Rotational speed [rad/s]
回転体間の角速度差を Symbol::omega コマンド
により、ギリシャ文字 ω として定義する。
MuPAD は ギ リ シ ャ 文 字 を
Symbol::alpha,
Symbol::beta などのように Symbol コマンドで記述
できる。
尚、 Symbol コマンドはギリシャ文字だけではなく
特殊文字や添え字なども設定できる。
コマンド入力時、Symbol:: や Symbol::om のよう
に途中までタイプしたところで Ctrl + Space を押すと、
それ以降に続く MuPAD コマンドの補間リストが表示
される。
:= （コロン・イコール）は MuPAD における代入処
理を意味し、この記号の右にあるものが、この記号の
左にある MuPAD 変数に代入される。
ここでは可読性の向上と、後の数式処理で ω を
容易に使えるよう、変数として w に代入している。
Notebook コマンド入力領域において、// （スラッ
シュ・スラッシュ） を使用することで、// 以降から改
行までをコメント化できる。
【R2012a 以前】
> Commands Tab > “Symbol” とタイプイン
【R2012b 以降】
>
> MuPAD > MuPAD 言語の基礎
> 記号と特殊文字 > 例および使用方法
> Non-Greek Characters in Text Regions
JMAAB Plant Model Workshop 3
15
3. Define parameter values
3. パラメータの定義
params :=
四角括弧 [ ] とカンマで、list を定義できる。list は
[// Breakaway friction torque [N*m]
BrkwyTrq = 25
要素の並び順も定義の一部となる。また、
params[1] で
// Coulomb friction torque [N*m]
ColTrq = 20
, ColTrq = 20
BrkwyTrq = 25,
params[2] で
というように list に定義した各要素を
取得できる。
// Viscous friction coefficient [N*m*s/rad]
, ViscCoef = 0.001
尚、コロンを伴わない等号記号 = は MuPAD で
// Transition approximation coefficient [s/rad]
, TransCoef = 10
は数学的な等号を意味し、MuPAD 変数への代入記
号 := （コロン・イコール）とは異なる。
// Linear region velocity threshold [rad/s]
数式処理において代入と等号を使い分けることは
, VelThr = 1e-4
重要になる。代入（:=）と等号（=）についての基本的
];
な考え方を Appendix
A.4 に記載した。
【R2012a 以前】
> Commands Tab > “ [ ] ” とタイプイン
【R2012b 以降】
> MuPAD Functions > _index を選択
4. Define equations
4. モデル方程式の定義
eqs := table(): // Define table for equations
モデルの方程式として最終的に Simscape で使
うものは list ではなく table にストアしておくことを推
奨する。table 変数は、table コマンドで定義できる。
尚、本例題では式が 1 つしかなく、また、式変形せ
ずに最終的な式を直接入力しているため特に table
を使う利点は分かりにくいが、説明の意味を込め、敢
えて table を使用している。
JMAAB Plant Model Workshop 3
16
// Equation of rotational friction torque
ここでの数式は、piecewise コマンドを用いた条件
eqs[`1.1`] :=
付きの数式として定義している。if 文も利用可能だ
T(t) = piecewise(
が、数式処理を念頭においた場合、piecewise コマ
[ (-VelThr<=w(t)) and (w(t)<=VelThr)
, BrkwyTrq*w(t)/VelThr ]
ンドの使用を推奨する。
尚、piecewise コマンドと if 文の違いについて
, [ w(t)>VelThr
Appendix
A.5 に記載している。
, ViscCoef*w(t)+ColTrq+(BrkwyTrq-ColTrq)*
exp(-TransCoef*w(t)) ]
piecewise コマンド中の w(t), T(t) の記述は時間
, [ Otherwise
関数を意味する。このように 「 記号（記号） 」 と記
, ViscCoef*w(t)-ColTrq-(BrkwyTrq-ColTrq)*
exp(-TransCoef*abs(w(t))) ]
述することによって、MuPAD はこれらを関数として
取り扱う。
):
ここでは、定義した数式を table に eqs[`1.1`] と
して保存している。list とは異なり table の添え字は
ユーザーがフレキシブルに決めることができる。
記号「`」バックティックは日本語キーボードの
SHIFT + @ で入力できる。このようにバックティックで
添 え 字 を 括 る こ と で 、 `1.1` の 部 分 を 例 え ば
`01.01.01` や `1.a` のように比較的自由な文字列
として使用できる。
eqs;
table にストアした式は、整った形で表示でき
る。
【R2012a 以前】
> Commands Tab > “table” とタイプイン
> Commands Tab > “piecewise” とタイプイン
【R2012b 以降】
> MuPAD Functions > table を選択
> MuPAD Functions > piecewise を選択
JMAAB Plant Model Workshop 3
17
5. Symbolic manipulation and evaluation
// Substitute params into eqs[1.1]
5. 数式処理と数式の評価
subs コマンドを用いて、定義した数式に対しパラメ
tmpeq := subs(eqs[`1.1`], params);
ータ値を代入する。
table へストアした方程式にはパラメータが値では
なく記号として含まれており、式が大きくなるほど、ま
たは記号が増えるほど、どのような関数形なのか把
握することが難しくなる。しかし、パラメータ値を代入
することで、関数形をより把握しやすくすることができ
る。
尚、この式は Simscape へ出力しないため table
へはストアせず、適当なテンポラリ変数 tmpeq に代
入している。
// Evaluation at an operating point
evalAt(tmpeq, w(t)=1.0);
続 い て 、 数 式 の 評 価 と し て evalAt コ マ ン ド 、
assuming コマンド、assume コマンドを使用する。
evalAt コマンドでは w(t) に対し一定値を代入し、
その値における式の評価結果を確認できる。
// assume w(t) > 10 with the assuming command
tmpeq
assuming(w(t) > 10);
assuming コマンド、及び、assume コマンドで、分
岐などを含む式に対し仮定を付不し、その条件下に
おける式を確認できる。
assuming コマンドは１つの処理に対して条件を設
定したい場合に使用する。
// assume w(t) > 10 with the assume command
複数の処理に対し継続的に条件を適用したい場
assume(w(t) > 10):
合は、assume コマンドと unassume コマンドを使用
tmpeq;
する。assume コマンドを用いた場合は、その後、
unassume コマンドにより条件を無効化することがで
きる。
// unassume for w(t)
unassume(w(t)):
【R2012a 以前】
tmpeq;
> Commands Tab > “subs” とタイプイン
> Commands Tab > “evalAt” とタイプイン
> Commands Tab > “assuming” とタイプイン
> Commands Tab > “assume” とタイプイン
> Commands Tab > “unassume” とタイプイン
【R2012b 以降】
> MuPAD Functions > subs を選択
> MuPAD Functions > evalAt を選択
> MuPAD Functions > assuming を選択
> MuPAD Functions > assume を選択
> MuPAD Functions > unassume を選択
JMAAB Plant Model Workshop 3
18
6. Visualize
6. グラフ表示
// Advance preparation for tmpeq before visualize
tgt := [lhs(tmpeq), rhs(tmpeq)]:
ここではトルクを角速度差の関数としてプロットす
る。
vars := select(select(indets(tgt, PolyExpr)
, testtype, "function")
, _not@testtype, "exp"):
repl := map(vars, v -> v=op(v,0));
先ず、時間の関数として定義されている記号から、
(t) を取り除く。この処理は、プロットしたい式に含ま
れている関数に対し MuPAD が割り当てている各種
の内部情報を利用し、関数 f(t) をただの記号 f で
置き換えるための式 f(t)=f を構築する方法である。
ｌhs コマンドで、等式の左辺を取得することができ
る。また、rhs コマンドで、等式の右辺を取得すること
ができる。
// Alternative procedure for above statement
// repl = { w(t)=w, T(t)=t };
//*
上記の手法は入門書で取り扱うレベルを超えてい
るため、一般性は失われるものの、代わりの簡易な
手法として、左記のコメント化しているステートメント（*
viseq := subs(tmpeq, repl);
記述箇所）のように、時間関数を変数に置き換える
式を直接入力しても良い。こちらの方法を使う場合
は、* の行の先頭にあるコメント記号 // を削除し、
実行する。
g1 := plot::Function2d( rhs(viseq)
グラフの表示には、plot::Function2d コマンド、及
, w=0..PI, LineStyle=Solid, LineWidth=0.5*unit::mm
, Color=RGB::Blue, LegendText="Rotational velocity"
):
び、plot コマンドを使用する。
plot::Function2d コマンドはグラフとして表示する
式に対し、軸の範囲・線種・線幅・凡例などを設定で
きる。
尚、plot::Function2d コマンドの中で、w の範囲
を w=0..PI として設定しているが、この PI は円周
率πを示す。
MATLAB では円周率を小文字で pi と入力する
が、MuPAD では大文字 PI と入力する。MATLAB
と MuPAD のコマンドの違いなどについて Appendix
A.1 に記述している。
JMAAB Plant Model Workshop 3
19
plot( g1
plot コマンドは、plot::Function2d コマンドで設定
, AxesTitles=["Rotational velocity[rad/s]"
,"Torque[N*m]"]
したグラフの表示や、タイトル・凡例の表示位置など
座標軸に関する設定を行う。
, YAxisTitleOrientation=Vertical
, XTicksNumber = Low, YTicksNumber = Normal
【R2012a 以前】
, LegendVisible=TRUE, LegendPlacement=Top
> Commands Tab > “lhs” とタイプイン
, LegendAlignment=Right
> Commands Tab > “rhs” とタイプイン
, GridVisible=TRUE, SubgridVisible = TRUE
> Commands Tab > “select” とタイプイン
, Header="Rotational Friction Torque"
> Commands Tab > “indets” とタイプイン
, BorderWidth=0.15*unit::mm
> Commands Tab > “map” とタイプイン
):
> Commands Tab > “op” とタイプイン
> Commands Tab > “->” とタイプイン
> Commands Tab
> “plot::Function2d” とタイプイン
> Commands Tab > “plot” とタイプイン
【R2012b 以降】
> MuPAD Functions > lhs を選択
> MuPAD Functions > rhs を選択
> MuPAD Functions > select を選択
> MuPAD Functions > indets を選択
> MuPAD Functions > map を選択
> MuPAD Functions > op を選択
> MuPAD Functions > proc を選択
(表示内容に
-> についても記載される)
> MuPAD Functions >plot::Function2d を選択
> MuPAD Functions > plot を選択
7. Generate Simscape code
7. Simscape 言語の出力
// Advance preparation before generating Simscape code
検討した数式を Simscape 言語として出力する。
ssceq := subs(eqs[`1.1`], repl);
先ず、Simscape 言語に変換できるよう、方程式
eqs[`1.1`] から、(t) を取り除いている。repl は、19
ページのグラフ表示で作成した置き換え用の変数で
ある。
尚、将来のリリースでは、このような処理は必要な
くなる予定である。
JMAAB Plant Model Workshop 3
20
print(Unquoted
MuPAD か ら Simscape 言 語 へ の 変 換 に は 、
,generate::Simscape(ssceq, NoWarning));
T == if (VelThr < omega),
generate::Simscape コマンドを用いる。
ただし、現在変換できるのは、Simscape 言語の
ColTrq+omega*ViscCoef+exp(-omega*TransCoef)
equations セクションに記述できる式のみである。
*(BrkwyTrq-ColTrq);
else
左記の出力結果を、Simscape Foundation Library
if ((0.0 <= omega+VelThr) && (omega <= VelThr)),
(omega*BrkwyTrq)/VelThr;
に含まれる Rotational Friction ブロック のソース
（friction.ssc）の中の equations セクションと比べる
else
と、完全に同一の表現ではないが、計算内容として
-ColTrq-exp(-TransCoef*abs(omega))*(BrkwyTrq-ColTrq)
+omega*ViscCoef;
は同じものが出力されていることが分かる。
尚、friction.ssc は、Rotational Friction ブロックを
end;
Simulink エディター上に配置し、このブロック上で右
end;
クリックして表示されるコンテキストメニューの中から
View Simscape source を選択することで開くことがで
きる。
また、Simscape 言語に関する詳細は、下記 URL
から参照できる。
http://www.mathworks.com/help/pdf_doc/physm
od/simscape/simscape_lang.pdf
fprint( Unquoted, Text, "FromMuPAD_01.txt"
出力 は、上 記 のように print コマンド を用いて
, generate::Simscape(ssceq, NoWarning)
Notebook に出力する方法と、左記のように fprint
):
コマンドを用いてテキストファイルに出力する方法が
ある。
【R2012a 以前】
> Commands Tab
> “generate::Simscape” とタイプイン
> Commands Tab > “print” とタイプイン
> Commands Tab > “fprint” とタイプイン
【R2012b 以降】
> MuPAD Functions > generate::Simscape を選択
> MuPAD Functions > print を選択
> MuPAD Functions > fprint を選択
JMAAB Plant Model Workshop 3
21
3.2.2 例題: DC モーターモデル
DC モーターは電気エネルギーと機械エネルギーを相互に変換する装置である。DC モーターの概略図を以下（左図）に
示す。尚、Simscape のオプション製品となる SimElectronics® では、DC モーターを模擬したブロックが用意されている（右
図）。また、Simscape においては、DC モーターに関する例題が用意されており、MATLAB コマンドウィンドウで ssc_dcmotor
とタイプインすることで該当モデルを開くことができる。更に、MathWorks 社ホームページのユーザーコミュニティサイト
（MATLAB Central）では、 Simscape 言語を用いて DC Motor をモデリングした例（下記 URL）が公開されている。
http://www.mathworks.com/matlabcentral/fileexchange/25014
DC モーター概略図 及び DC Motor ブロック
3.2.2.1
DC モーターモデルの数式処理と解説
DC モーターモデルの式は、キルヒホッフの法則から導出される電圧と電流に関する電気系の式と、ニュートンの法則から
導出されるトルクに関する機械系の式から成り、常微分方程式として表現される。
DC モーターモデル（微分方程式）の数式処理について以下に示す。また、本例題の流れは、プラントモデリングのアプロ
ーチ「① 数式的アプローチ」の Case 2、及び、「② 数値的アプローチ」の Case 3 に該当する。
尚、前節の回転摩擦トルクモデルで解説した内容はここでは省略している。また、前節同様に、左の段落は MuPAD
Notebook の各入力領域を意味し、右の段落は MuPAD コマンドに対する解説である。
DC Motor
Initialize MuPAD session
reset():
sysname(Arch);
version();
Define parameter symbols
w := Symbol::omega; // Motor angular velocity [rad/s]
tau := Symbol::tau; // Motor torque [N*m]
JMAAB Plant Model Workshop 3
22
1. Define parameter values
1. パラメータの定義
params :=
MuPAD では、D は微分のコマンドとして D(x) な
[ R
= 3.9
// Armature resistance [Ohm]
どのように使われるため、ユーザー定義の記号として
, L
= 1.2e-5
// Armature inductance [H]
D の利用は推奨できない。ここでは、モーター軸の
, J
= 1e-6
// Inertia [kg*m^2]
粘性減衰係数を、 Dr という記号を用いて定義し
, Dr = 3e-6
// Rotor damping [N*m/(rad/s)]
た。
, Ki = 7.2e-5
// Torque constant [N*m/A]
, Kb = 7.2e-5
// Back-emf constant [V/(rad/s)]
];
【R2012a 以前】
> Commands Tab > “D” とタイプイン
【R2012b 以降】
> MuPAD Functions > D を選択
2. Define equations
2. モデル方程式の定義
eqs := table(): // Define table for equations
複数の方程式で物理システムを表現する場合、
物理法則をそのまま記述したような基本的な式に加
え、それらの式を変形した後の式もシステムの方程
式として用いたい場合がある。
式変形を伴う場合には、table オブジェクトだけで
なく作業用に別の MuPAD 変数も作成し、変形中
の式は作業用の MuPAD 変数へ代入することを推
奨する。
ただし、本例題では導出過程を省いて最終的な
モデル式を定義しているため、table へ直接代入を
行っている。
// Dynamics of rotational motion
DC モーターの式は、トルクに関する機械系の式
eqs[`1.1`] := J*diff(w(t),t) = - Dr*w(t) + tau(t);
eqs[`1.1`], eqs[`1.2`] と、電圧と電流に関する電気
系の式 eqs[`2.1`], eqs[`2.2`] として定義できる。
ここでは、 eqs[`1.2`] で定義しているように、モー
eqs[`1.2`] := tau(t) = Ki * i(t);
ター発生トルクは電流に比例するとした。
JMAAB Plant Model Workshop 3
23
// Dynamics of electrical circuit
eqs[`2.1`], eqs[`2.2`] において、V(t) は印加電圧、
eqs[`2.1`] := L*diff(i(t),t) + R*i(t) = V(t) - Vb(t);
Vb(t) はモーター逆起電力を表す。
これらの式は回転角速度 ω と電流 i に関する
eqs[`2.2`] := Vb(t) = Kb*w(t);
微分方程式として表されるが、微分は diff コマンド
を用いて表現できる。
【R2012a 以前】
> Commands Tab > “diff” とタイプイン
【R2012b 以降】
> MuPAD Functions > diff を選択
3. Symbolic manipulation
3. 数式処理
前項で定義したモデル式に対して数式処理を適
用する場合、table ではなく list を使う方が便利で
ある。この場合 table から list への変換作業が必
要になるものの、添え字の使い方に関して自由度の
高い table と、数式処理の適用に関して自由度の
高い list を使い分ける方が、取り扱うモデルの規模
が大きくなっても効率的に作業ができる。
// Get the index from equation table
idx := map([op(eqs)], lhs);
table から list へ変換するため、先ず map コマ
ンドと op コマンドを用いて、 table オブジェクト eqs
の添え字を取得し idx に代入している。
map コマンドは、第 1 引数のオブジェクトに対して、
// Convert table to list with idx
第 2 引数の処理を順次適用する。
eqs2 := []:
op コマンドは、数式や記号の構成要素から任意
for k in idx do
の要素を取り出すなど、MuPAD が構成する様々な
eqs2 := [op(eqs2), eqs[k]]:
オブジェクトから必要な要素を取り出す。
end_for:
table か ら list へ は 、 for…end_for ル ー プ と 、
// Check the number of equations
idx を用いて変換し、 eqs2 に代入している。また、
nops(eqs2) * ` equations`;
ここでは list に変換された数式の数を nops コマン
ドを用いて確認している。
JMAAB Plant Model Workshop 3
24
eqs2
list は table と異なり、要素ごとに改行されること
なく表示される。
eqs3 :=
subs コマンドを用いて、中間変数（トルク τ(t) と
[ subs(eqs2[1], eqs2[2])
内部起電力 Vb(t)）を消去し、状態変数（角速度
, subs(eqs2[3], eqs2[4])
ω(t) と電流 i(t)）および外部入力（電圧 V(t)）だけ
];
の式に変換する。
【R2012a 以前】
> Commands Tab > “map” タイプイン
> Commands Tab > “op” とタイプイン
> Commands Tab > “for” とタイプイン
> Commands Tab > “nops” とタイプイン
【R2012b 以降】
> MuPAD Functions > map を選択
> MuPAD Functions > op を選択
> MuPAD Functions > for を選択
> MuPAD Functions > nops を選択
4. Analyze steady-state
4. 定常状態の分析
// Advance preparation
定常状態では物理量が時間とともに変化せず、
tgt := [map(op(eqs3), lhs), map(op(eqs3), rhs)]:
vars := select(select(indets(tgt, PolyExpr)
, testtype, "function")
, _not@testtype, "diff"):
一定の値をとる。
そこで、定常状態を分析するための準備として、
先ず、時間の関数を定数とみなし、時間関数から定
数に置き換えるための変数 SS を用意する。
SS := map(vars, v -> v=op(v,0));
ここでは、前節のグラフ表示において用いた方法
と同様に、時間関数として定義している変数から (t)
を取り除く方法を用いている。
// Alternative procedure for above statement
// SS := [ w(t)=w , i(t)=i , V(t)=V ];
//*
上記の手法は入門書で取り扱うレベルを超えてい
るため、簡易な代替え手法として、左記のコメント化
しているステートメント（* 記述箇所）のように、時間
関数を変数に置き換える式を直接入力してもよい。
JMAAB Plant Model Workshop 3
25
// Get steady-state equations as tmp1
tmp1 := evalAt( eqs3, SS );
SS を eqs3 に代入し、evalAt コマンドにより評価
すると、微分項はゼロとなり定常状態の方程式
tmp1 が得られる。
// Substitute params into tmp1
得られた定常状態の方程式に、パラメータ値を代
tmp2 := evalAt( tmp1, params );
入すると、定常状態の状態変数と外部入力の関係
式 tmp2 が得られる。
// Simultaneous equation with tmp2
この関係式から状態変数について陽に解くと、状
SSans := op(solve(tmp2, [w,i]));
態変数の定常値は外部入力の関数 SSans のよう
に求められる。この結果を、数値計算（シミュレーショ
ン）の正しさを確認する手段として使うことができる。
【R2012a 以前】
> Commands Tab > “solve” とタイプイン
【R2012b 以降】
> MuPAD Functions > solve を選択
5. Numerically integrate differential equations
5. 微分方程式の数値解法
数式処理ソフトを用いて微分方程式を解く場合、
解析解を得ることができる場合もあるが、一般的に
は実際の開発で用いるプラントモデルは非線形であ
ることが多く、外部入力に対するシステムの挙動を
解析解という形で得られるのは稀である。そのため、
数値計算（シミュレーション）により解を求める必要が
あり、MuPAD では numeric::ode2vectorfield コマン
ドと numeric::odesolve2 コマンドを使用する。
// Define input
ここでは先ず、外部入力となるモーターへの印加
Input := V(t) = 1.0;
電圧 V(t) が一定であると仮定し、Input を定義す
る。
// Define initial conditions
状態変数に対し、初期条件 (initial conditions)
ICs := w(0)=0 , i(0)=0 ;
を ICs として設定する。
JMAAB Plant Model Workshop 3
26
// Substitute params and initial conditions into eqs3
syseqs := op(evalAt(eqs3, params, Input));
ここで、システム方程式 eqs3 にパラメータ値
params、及び、外部入力 Input を代入し、 syseqs
を得る。
// Define differential variables
numeric::ode2vectorfield コマンドには、状態変数
DVs := [w(t),i(t)];
(differential variables) を入力引数として渡す必要
があるため、DVs を設定する。尚、ここでは状態変
数を直接入力しているが、数式処理を用いてより一
般的に設定する方法を、Appendix
A.6 に記載し
た。
// Convert the initial value problem and the fields
numeric::ode2vectorfield コマンドに、システム方
// into a procedure acceptable by numeric::odesolve2.
程式 syseqs 、初期値 ICs 、状態変数 DVs を渡し、
dsys1 := numeric::ode2vectorfield([syseqs,ICs], DVs):
戻り値を dsys1 にストアする。
// Get a function of numeric approximate solution
dsys2 := numeric::odesolve2(dsys1, GAUSS(3)):
次 に 、 numeric::odesolve2
コ マ ン ド に 、
numeric::ode2vectorfield コ マ ン ド の 戻 り 値 dsys1
と数値計算のアルゴリズム（ソルバー）を設定する。
ここではスティッフソルバーである GAUSS を指定
している。ソルバーを指定しない場合 Runge-Kutta
に基づく DOPRI78 が自動的に選択されるが、物理
システムの数値計算に使用するには必ずしも適切と
は限らない。
// Evaluate numerically
numeric::odesolve2 コマンドの戻り値 dsys2 を用
DVs;
いて、指定する時刻（例として 0.1 秒、2 秒、8 秒）に
tt := time():
おける数値計算結果を得ることができる。
dsys2(0.1);
dsys2 か ら 得 ら れ る 各 要 素 の 並 び 順 は 、
dsys2(2);
numeric::ode2vectorfield に対して指定した状態変
dsys2(8);
数の list （DVs） に対応している。そこで、最初に
print( (time() - tt)/1000, " seconds");
DVs を表示するようにしている。また、計算時間を計
測するため、time() コマンドを用いている。
JMAAB Plant Model Workshop 3
27
// Show again SSans (analytic solution for steady state)
SSans;
定常状態の解析解と、数値解を比較するため、
定常状態の解析解 SSans を再度表示する。この結
果より、入力電圧 V=1 に対する数値計算は正しい
定常状態に到達していることが分かる。
【R2012a 以前】
> Commands Tab
> “numeric::odesolve2” とタイプイン
> Commands Tab
> “numeric::ode2vectorfield” とタイプイン
【R2012b 以降】
> MuPAD Functions > numeric::odesolve2 を選択
> MuPAD Functions
> numeric::ode2vectorfield を選択
6. Visualize
6. グラフ表示
g1 := plot::Function2d( dsys2(t)[1]
ここでは、ここまでの分析で得られた電機子電流
, t=0..2, LineStyle=Solid, LineWidth=0.5*unit::mm
i(t) と回転角速度 ω(t) の時間発展と、それぞれ
, Color=RGB::Blue, LegendText="w(t)"
の定常状態をグラフとして表示する。
, Mesh=121, AdaptiveMesh=0
):
plot::Function2d コマンドの戻り値を plot コマンド
に複数設定することで、1 つのプロット上に複数のグ
g2 := plot::Function2d( subs(rhs(SSans[1]), V=1)
ラフを表示できる。
, t=1..2, LineStyle=Dotted, LineWidth=0.7*unit::mm
, Color=RGB::Red, LegendText="w(t): steady state"
, Mesh=121, AdaptiveMesh=0
こ こ で は 、 plot::Function2d コ マ ン ド に 対 し て 、
AdaptiveMesh オプションを設定している。
):
この AdaptiveMesh オプションは設定するメッシ
ュサイズ（Mesh オプションで設定）に対し、補間アル
ゴリズムを適用するか否か、及び、補間レベルを設
定するオプションとなる。
数式の振る舞いにもよるが、ここでは補間アルゴリ
ズムの副作用を防ぐため、AdaptiveMesh=0 と設定
し、補間アルゴリズムを適応しないようにしている。
JMAAB Plant Model Workshop 3
28
g3 := plot::Function2d( dsys2(t)[2]
, t=0..2, LineStyle=Solid, LineWidth=0.5*unit::mm
, Color=RGB::Green, LegendText="i(t)"
, Mesh=121, AdaptiveMesh=0
):
g4 := plot::Function2d( subs(rhs(SSans[2]), V=1)
, t=1..2, LineStyle=Dotted, LineWidth=0.7*unit::mm
, Color=RGB::Orange, LegendText="i(t): steady state"
, Mesh=121, AdaptiveMesh=0
):
plot(g1, g2
, AxesTitles=["time [s]", "Angular velocity [rad/s]"]
, YAxisTitleOrientation=Vertical
, XTicksNumber = Low, YTicksNumber = Normal
, LegendVisible=TRUE,LegendPlacement=Top
, LegendAlignment=Right
, GridVisible=TRUE, SubgridVisible = TRUE
, Header="DC Motor - angular velocity -"
, BorderWidth=0.15*unit::mm
):
plot( g3,g4
, AxesTitles=["time [s]", "Current [A]"]
, YAxisTitleOrientation=Vertical
, XTicksNumber = Low, YTicksNumber = Normal
, LegendVisible=TRUE,LegendPlacement=Top
, LegendAlignment=Right
, GridVisible=TRUE, SubgridVisible = TRUE
, Header="DC Motor - current -"
, BorderWidth=0.15*unit::mm
):
JMAAB Plant Model Workshop 3
29
7. Generate Simscape code
7. Simscape 言語の出力
// Define procedure in order to generate Simscape code
PREP := proc(eqs)
ここでは
eqs[`1.1`] ～ eqs[`2.2`] の 4 式 を
Simscape 言語に変換する。
begin
前節同様に、Simscape 言語に変換できるよう、
tmpeq :=
数式を書き換えている。詳細な説明はここでは省略
eval(text2expr("eqs"."[".expr2text(idx[NN])."]"))
するが、書き換え用のプロシージャ（サブ関数のよう
$NN=1..nops(eqs):
map(tmpeq
なもの）を proc…end_proc コマンドで定義している。
尚、将来的には、このような処理についての記述
, n -> text2expr( stringlib::subs( expr2text(n)
は必要なくなる予定である。
, "`&"=""
, ";`"="", "(t)"=""
)
)
)
end_proc:
JMAAB Plant Model Workshop 3
30
print(Unquoted
Simscape 言語へ変換した結果を出力するには、
, generate::Simscape(PREP(eqs), NoWarning)):
前節同様に print コマンド、もしくは fprint コマンド
J*omega.der == tau-Dr*omega;
を使用する。また、これらのコマンドの中で PREP プ
tau == Ki*i;
ロシージャを呼ぶ。
R*i+L*i.der == V-Vb;
Vb == Kb*omega;
FileName := "FromMuPAD_02.txt":
fprint コマンドで FromMuPAD_02.txt ファイルへ
fid := fopen(FileName, Append, Text):
変 換 結 果 を 保 存 す る 際 、 こ こ で は equations と
fprint(Unquoted, fid,"equations"):
end キーワードも合わせて出力している。
fprint( Unquoted, fid
, generate::Simscape(PREP(eqs), NoWarning)):
fprint(Unquoted, fid,"end"):
fclose(fid):
3.2.2 節で紹介した MATLAB Central に掲載され
ている
DC Motor モ デ ル （ DCMotor.ssc ） の
equations セクションは、本例題で出力した当該箇
所と比較して、回転数 ω(t) の微分項が含まれて
いない（別途イナーシャコンポーネントを接続する設
計になっている）ことを除き、同じ計算内容になって
いる。
【R2012a 以前】
> Commands Tab > “proc” とタイプイン
> Commands Tab > “system” とタイプイン
> Commands Tab > “fopen” とタイプイン
> Commands Tab > “fclose” とタイプイン
【R2012b 以降】
> MuPAD Functions > proc を選択
> MuPAD Functions > system を選択
> MuPAD Functions > fopen を選択
> MuPAD Functions > fclose を選択
JMAAB Plant Model Workshop 3
31
Appendix
JMAAB Plant Model Workshop 3
32
A.1 MATLAB® と MuPAD® の関係と違い
MuPAD は、MATLAB のオプション製品となる Symbolic Math Toolbox に含まれる数式処理専用 GUI である。しかしな
がら現状は MATLAB から独立したアプリケーションのため、類似するコマンドの定義方法や使い方が一部異なっている。例
えば、円周率πは、MATLAB では小文字 pi で表現し、 MuPAD では大文字 PI で表現する。
MATLAB と MuPAD の違いに関する詳細は以下から確認できる。尚、ヘルプの使い方については、Appendix A.2 に記
述している。
【R2012a 以前】
MATLAB GUI
> コンテンツタブ > Symbolic Math Toolbox > User’s Guide
> MuPAD in Symbolic Math Toolbox Software > Integration of MuPAD and MATLAB
【R2012b 以降】
>
> MATLAB と MuPAD の統合 > MuPAD 関数の呼び出し> MATLAB 構文と MuPAD 構文の違い
また、MATLAB と MuPAD ではプロンプトについても違いがある。
MATLAB では、図 A1 に示すように、パラメータや関数、或いは、式を定義するには = （イコール）を用い、C 言語などの
一般的な手続き型プログラミング言語と同様に、右辺を左辺に代入するという考えに基づいてプログラミングを行う。
図 A1 :
MATLAB プログラミング例
MATLAB の基本的な操作方法などに関しては、下記 URL から参照できる。
http://www.mathworks.com/help/matlab/getting-started-with-matlab.html
JMAAB Plant Model Workshop 3
33
一方、MuPAD において変数や数式などをある変数に代入する（定義する）には、 := （コロン・イコール）を使用する。尚、
MuPAD において = （イコール）は数学的な等号を意味する。
図 A1 で示した MATLAB のコードに対し、数式処理ソフトである MuPAD を用いて同様の処理を行った場合、図 A2 に
示すように、数式をシンボリックに扱うことができる。
図 A2 :
MuPAD プログラミング例
図 A1 と図 A2 を比べると分かるように、数式処理を使うと、定数（A, B, f）と関数（x(t), x1(t), x2(t)）の区別が記号を見る
だけで明らかになる。また、記号としてパラメータが含まれる式（Eq2, Eq3）や、パラメータを値に置き換えた式（Eq4）のように、
直接的に式を確認することができる。
JMAAB Plant Model Workshop 3
34
A.2 MuPAD® ヘルプの使い方
ここでは、MuPAD のヘルプウィンドウの起動方法、及び、使い方について説明する。
尚、R2012a 迄のバージョンでは、Help Browser と呼ばれる MuPAD 専用のヘルプウィンドウが用意されている。R2012b
以降は仕様変更に伴い Documentation Center と呼ばれるヘルプウィンドウが起動する。
A.2.1
ヘルプウィンドウの起動
次の 3 つの方法で MuPAD のヘルプウィンドウを起動できる。
□
MuPAD Notebook の Help メニューバーから、Open Help を選択
□
MuPAD Notebook の Standard toolbar から、アイコン
□
MuPAD Notebook が起動している状態で F1 キーを押す
を選択
JMAAB Plant Model Workshop 3
35
A.2.2
ヘルプウィンドウの使い方
R2012a 迄のバージョンと R2012b 以降のバージョンではヘルプウィンドウが異なる。ここでは、R2012a 以前のバージョン
とそれより後のバージョンに別けて使い方を説明する。
【R2012a 以前】
R2012a 以前のバージョンでは、次のようなヘルプウィンドウ（Help Browser）が起動する。
Help Browser の具体的な使用方法として、Commands タブの使い方について説明する。
Commands タブは、あらかじめコマンド名を把握している場合に使用する。例えば、solve コマンドの詳細な利用方法を確
認するには、以下のようになる。
1) Help Browser の Commands Tab を選択
2) Enter command to search for: エディタブルテキストボックスに “solve” を入力
3) Select command to display: に表示する関連項目から solve を選択
JMAAB Plant Model Workshop 3
36
その他、MuPAD の Help および Help Browser に関する情報は以下から確認できる。
> Contents Tab > Getting Started > Accessing MuPAD Help
【R2012b 以降】
R2012b 以降のバージョンでは、次のようなヘルプウィンドウ（Documentation Center）が起動する。
Documentation Center の具体的な使用方法について説明する。
前述のように予め solve コマンドを把握している前提で、Documentation Center の利用方法について説明する。
JMAAB Plant Model Workshop 3
37
1) Documentation Center の エディタブルテキストボックスに solve と入力し Enter を押す
2) カテゴリで絞り込む から MuPAD を選択
3) 表示される該当項目から solver を選択
JMAAB Plant Model Workshop 3
38
尚、Documentation Center を 3 の状態に保つことで、MuPAD に関連する項目に絞り検索できる。
その他、MuPAD の Help に関する情報は以下から確認できる。
>
> MuPAD > 例および使用方法 > Access Help for Particular Command
また、MuPAD のコマンド一覧が、ヘルプウィンドウの下方
JMAAB Plant Model Workshop 3
から確認できる。
39
A.3 Notebook の実行方法
MuPAD Notebook 内で、ある変数の値を変更しその変更を Notebook 全体に伝播させるには、 Notebook メニューか
らアクセス可能な次に示す方法を利用すると便利である。

Evaluate
Enter と等価

Evaluate From Beginning
Notebook の最初からカーソル位置まで実行

Evaluate To End
カーソル位置から Notebook の最後まで実行

Evaluate All
Notebook 全てを実行
JMAAB Plant Model Workshop 3
40
A.4 代入 （:=） と 等号 （=） の考え方
MuPAD を用いたモデルの定義や式の処理では、代入（:=）と等号（=）を使い分けることが重要である。特にモデルパラメ
ータの扱いでは、等号を使用することが強く推奨される。以下に、例を示しながら代入と等号の考え方について説明する。
Notebook を初期化する。
reset(): sysname(Arch), version();
パラメータを定義する。パラメータのシンボルと値は等号（=）を使って表現することを推奨する。次の例では、 J_1 と Dr_1 がシンボル、 1e-6
と 3e-6 が値である。
params1 :=
[ J_1 = 1e-6
, Dr_1 = 3e-6
];
パラメータを定義するもう１つの方法として、以下のように MuPAD 変数をパラメータシンボルとして使い、値を代入（:=）する方法もある。
しかしながら、後述する一連の理由により、この方法は推奨しない。
J_2 := 1e-6;
Dr_2 := 3e-6;
等式（=）によるアプローチでは、次のようにパラメータシンボルを使って式を定義したとき、 MuPAD からの出力表示は、パラメータがシンボル
のまま保持される。
eq1 := J_1*diff(w(t),t) = - Dr_1*w(t) + trq(t);
代入（:=）によるアプローチでは、 MuPAD が返す式の中のパラメータは、自動的に値に置き換えられる。
eq2 := J_2*diff(w(t),t) = - Dr_2*w(t) + trq(t);
等式（=）によるアプローチでは、式の中のパラメータシンボルを値に置き換えるには、代入コマンド （evalAt） を明示的に実行することになる。
evalAt(eq1, params1);
JMAAB Plant Model Workshop 3
41
次に示すように、外部ファイルにパラメータ情報を保存する場合、代入（:=）によるアプローチでは値だけが保存され、対応するシンボルは保存
されないため情報が欠損してしまう （params2.txt の内容を参照）。しかし、等式（=）によるアプローチでは、シンボルと値が等号で結ばれた形で
保存され、情報の欠損は発生しない （params1.txt の内容を参照）。
fid := fopen("params1.txt", Write, Text):
fprint(fid, params1, eq1);
fclose(fid);
等式によるアプローチの出力結果 （param1.txt）
fid := fopen("params2.txt", Write, Text):
fprint(fid, J_2, Dr_2, eq2);
fclose(fid);
代入によるアプローチの出力結果 （param2.txt）
上記の例を見ると、代入（:=）によるアプローチの方が手数が尐なく一見効率的に見えるかもしれない。しかし、代入アプロ
ーチでは、ユーザーの入力と MuPAD の出力を両方見なければ一連の処理の流れを理解することができない。
一方、等式（=）によるアプローチでは MuPAD の出力だけを見て処理を理解することができる。これは、より複雑な式を扱
ったり、より多くの式を扱ったりするときほど重要となる。また、数式処理を適用する上では、式の中でパラメータが記号として
保持されていることが重要になる場合がある。
以上の理由から、等式（=）によるアプローチを一貫して採用することを強く推奨する。
JMAAB Plant Model Workshop 3
42
A.5 piecewise と if … then … の違い
MuPAD における条件分岐の記述方法について補足する。
MuPAD における条件分岐の表現方法としては、if…then… コマンドを用いた方法と piecewise コマンドを用いる方法があ
る。if 文は条件分岐プログラム用途に向いており、数式処理を前提とする場合 piecewise コマンドの使用を推奨する。これ
らのコマンドの違いについて以下に示す。
Non-linear spring model
Initialize notebook
reset(): sysname(Arch), version();
Define parameters
params := [
K1=10
// spring stiffness (Low rate)
,K2=1000
// spring stiffness (High rate)
,SWP = 0.05
// Switching point between K1 and K2
];
Define nonlinear spring equations
1) Define with if…then…elif…end_if
spring1 :=
x -> if abs(x)<SWP then
K1*x
elif x>SWP then
sign(x)*(K1*SWP+K2*(x-SWP)^2)
elif x<-SWP then
sign(x)*(K1*SWP+K2*(x+SWP)^2)
else
0
end_if;
x -> (if abs(x) < SWP the...
JMAAB Plant Model Workshop 3
43
The default output just shows a part of the entire procedure when defining MuPAD procedure. Therefore,
if you want to see all of the procedure, you can see it with expose command.
expose(spring1);
x -> (if abs(x) < SWP then
K1*x
else
if SWP < x then
sign(x)*(K1*SWP + K2*(x - SWP)^2)
else
if x < -SWP then
sign(x)*(K1*SWP + K2*(x + SWP)^2)
else
0
end_if
end_if
end_if)
2) Define with piecewise function
spring2 := piecewise(
[ abs(x)<=SWP, K1*x]
, [ x>SWP, sign(x)*(K1*SWP+K2*(x-SWP)^2)]
, [ x<-SWP, sign(x)*(K1*SWP+K2*(x+SWP)^2)]
, [ Otherwise, 0]
);
Define list of equations
ListSpringEqs := [spring1, spring2]:
JMAAB Plant Model Workshop 3
44
Evaluate the equations
Integrate equations
int(ListSpringEqs[1],x);
Error: The integrand is invalid. [int]
int(ListSpringEqs[2],x);
Substitute params to equations
tmpSpringEqs := subs(ListSpringEqs, params):
Evaluate equation when assuming x>0.1
assume(x>0.1):
expose(tmpSpringEqs[1]);
x -> (if abs(x) < 0.05 then
10*x
else
if 0.05 < x then
sign(x)*(10*0.05 + 1000*(x - 0.05)^2)
else
if x < -0.05 then
sign(x)*(10*0.05 + 1000*(x + 0.05)^2)
else
0
end_if
end_if
end_if)
tmpSpringEqs[2];
unassume(x):
JMAAB Plant Model Workshop 3
45
以上のように、piecewise を利用した場合のみ、int コマンドによる積分などの数式処理や、assume コマンドで仮定条件を
設定したときの数式の評価が可能となる。
尚、ここではグラフ化までは行っていないが、どちらもグラフ化は可能である。
【R2012a 以前】
> Commands Tab > “if” とタイプイン > if を選択
> Commands Tab > “->” とタイプイン > ->を選択
> Commands Tab > “expose” とタイプイン > ->を選択
【R2012b 以降】
>
> if を選択
>
> proc を選択
>
> expose を選択
JMAAB Plant Model Workshop 3
46
A.6 数式処理による状態変数の取得
DC モーターモデルの例では、微分方程式の解法において、状態変数 DVs を直接入力し設定を行った（p.27）。しかし、
状態変数の数が多い大規模なモデルなどでは直接入力する方法は現実的ではない。そこで、数式処理コマンドを用いて、
普遍的に状態変数を抽出する方法を以下に示す。尚、以下の方法は DC モーターモデル例における 「5.微分方程式の解
法」 の箇所で実行することを前提に記述している。
準備としてシステム方程式を、方程式の右辺と左辺を個別に要素として持つ list に変換する。
tmp1 := map(eqs3, eq -> (lhs(eq),rhs(eq)) );
続いて、因子を要素として持つ list に変換する。
tmp2 := map(tmp1, eq -> op(indets(eq,RatExpr)));
更に、重複する要素が存在する場合を考慮し、list を set （集合）に変換して重複要素を無くす。この時点で、システム方程式に含まれてい
る全てのパラメータと関数が抽出される。
tmp3 := coerce(tmp2,DOM_SET);
関数だけを取り出し、更に微分演算子を含む要素を取り出す。
vars := select(tmp3, testtype, "function");
diffs := select(vars, has, diff);
その結果から微分演算子を除去することで、状態変数 （differential variables） を得ることができる。
DVs := map(diffs, op) minus {t};
JMAAB Plant Model Workshop 3
47
参考までに、代数変数 （algebraic variables） は以下の集合演算によって得ることができる。
AVs := vars minus diffs minus DVs;
上述の通り、数式処理を用いることで、モデル規模によらず確実に状態変数を抽出することができる。
【R2012a 以前】
> Contents Tab > Getting Started > Programing Basics > Using Procedures
> Contents Tab > Getting Started > Using the MuPAD Editor to Write New Functions > Calling Functions Created in the MuPAD Editor
【R2012b 以降】
>
> MATLAB と MuPAD の統合 > MuPAD 関数の呼び出し> 独自の MuPAD プロシージャの使用
>
> MATLAB と MuPAD の統合 > MuPAD のファイルとインターフェイス
JMAAB Plant Model Workshop 3
48
A.7 MuPAD® mu ファイルについて
MuPAD の Notebook ファイルの拡張子は mn である。これに対し、MuPAD プロシージャ（proc…end_proc） を保存す
るときには、拡張子を mu とする。この mu ファイルはテキストファイルであり、MATLAB エディターやメモ帳などで編集でき
る。また、このファイルは、Notebook （mn ファイル） から呼び出して使う。
mu ファイルに関する詳細は以下から確認できる。
【R2012a 以前】
> Contents Tab > Getting Started > Programing Basics > Using Procedures
> Contents Tab > Getting Started > Using the MuPAD Editor to Write New Functions > Calling Functions Created in the MuPAD Editor
【R2012b 以降】
>
> MATLAB と MuPAD の統合 > MuPAD 関数の呼び出し> 独自の MuPAD プロシージャの使用
>
> MATLAB と MuPAD の統合 > MuPAD のファイルとインターフェイス
JMAAB Plant Model Workshop 3
49
参考情報
【関連図書】
生越 茂樹 : 基礎からの MuPAD, 工学社 (2004)
Christopher Creutzig, Walter Oevel : MuPAD Tutorial, Springer (2004)
William J. Palm III : System Dynamics -2nd Edition-, McGraw-Hill (2009)
【Web セミナ】
https://www.mathworks.co.jp/company/events/webinars/wbnr59783.html
https://www.mathworks.co.jp/company/events/webinars/wbnr56081.html
JMAAB Plant Model Workshop 3
50