基礎から学ぶDCモータの1Dモデル化② ～MATLAB/Simulink～

本記事は「基礎から学ぶDCモータの1Dモデル化～MATLAB/Simulink～」の２回目です。

前回の記事ではフレミング左手の法則について復習し、DCモータの回転原理について解説してそれを数理式（数理モデル）に置き換えてMATLAB/Simlinkでモデル化する方法を紹介しました。

DCモータの数理式は線形の微分方程式ですので、行列を用いるとシンプルに記述することができます。また、行列を用いて各状態量のつり合いを取りながら解くことで数値計算の精度向上も期待できます。（力のつり合い計算への応用等で有用）

本記事では、DCモータの1Dモデルを現代制御理論の分野で用いる状態方程式に直すことで、MATLAB/Simlinkのモデルをシンプルに記述する方法について解説します。

前回の記事

本記事はこちらの記事の続きとなります。

2023年1月8日基礎から学ぶDCモータの1Dモデル化① ～MATLAB/Simulink～

制御理論の概要

DCモータを1Dモデル化する目的としては、以下のものが考えられるかと思います。

ロボットなどでアクチュエータとして使われているDCモータを制御する制御モデルを構築したい。
モデルベース開発（MBD）による開発効率向上のため、プラントモデルを構築したい。

どちらの場合においても、MATLABでモデル化をするときは制御理論（Control theory）が基礎になりますので押さえておく必要があります。

ここで制御とは、対象を希望通りに「制（せい）し、御（ぎょ）する」ことを言い、電流・電圧・変位・速度・加速度などの物理量を目標値（希望する値）に一致させる操作のことを言います。

制　勝手なふるまいを押さえる
御　思うように操る

なお、制御モデルとプラントモデルの違いについては、以下の記事をご覧ください。

2022年7月13日 MBD(ﾓﾃﾞﾙﾍﾞｰｽ開発)におけるﾌﾟﾗﾝﾄﾓﾃﾞﾙ I/F ｶﾞｲﾄﾞﾗｲﾝ

制御理論の歴史

制御理論には大きく分けて古典制御理論（Classical control theory）と現代制御理論（Modern control theory）の２種類があります。

古典制御理論は１８世紀後半に、蒸気機関の回転速度を一定にするためにジェームズ・ワットが実用化した遠心調速機が始まりとされています。

それから理論的研究が進んでいき、古典制御理論は1950年代ごろに確立していきました。

1950年代に入ると第二次世界大戦で戦勝国となったアメリカとソ連は宇宙開発競争でしのぎを削ります。

ロケットを制御するときに位置や姿勢の測定にはセンサーが用いられますが、センサーの値はノイズだらけです。その値をそのまま用いると、宇宙飛行のような長距離移動では誤差が蓄積していき、目的地にたどり着くことは到底できません。

しかし、1960年になるとそうしたノイズをキャンセルするカルマンフィルタ（Kalman filter）が発明されます。これにより、1969年7月20日にアポロ11号は人類初の月面着陸に成功します。

Apollo Guidance Computer — Fig. 2 カルマンフィルタが搭載された
アポロ・ガイダンス・コンピュータ（AGC）

カルマンはたくさんの現代制御理論の重要な論文を発表し、現代制御理論の基本的な枠組みをほとんど１人で作り上げました。そのため、カルマンは「現代制御理論の父」と呼ばれています。

古典制御理論と現代制御理論

数理モデルの表現には伝達関数表現と状態空間表現があります。

伝達関数とは、入出力を表す関数におけるすべての変数の初期値がゼロであるとして、ラプラス変換したときの入力と出力の比のことを言います。

出力を$y(t)$のラプラス変換を$Y(s)$、入力$u(t)$のラプラス変換を$U(s)$とすると伝達関数$G(s)$は下式で表されます。

$$ G(s) = \frac{Y(s)}{U(s)} = \frac{\mathcal{L}[y(t)]}{\mathcal{L}[u(t)]} 　\tag{1} $$

よって、伝達関数表現は下図のように入力変数$U(s)$から直接、出力変数$Y(s)$へ影響を与える式で表されます。

Transfer function representation — Fig. 3 伝達関数表現

これに対し、状態空間表現は入力変数$\boldsymbol{u}(t)$から状態変数$\boldsymbol{x}(t)$までの影響を状態方程式で表し、その後、状態変数$\boldsymbol{x}(t)$と出力変数$\boldsymbol{y}(t)$とのつながりを出力方程式で表す方法となります。

State-space representation — Fig. 4 状態空間表現

古典制御理論では伝達関数表現を用い、現代制御理論では状態空間表現を用います。これにより古典制御理論では１入力１出力システム（SISOシステム）が基本となります。

現代制御理論では行列やベクトルを用いた式を使うことで、多入力多出力システム（MIMOシステム）を扱うことができます。

状態空間表現を用いたDCモータの1Dモデル化

では、前回の記事で導出した数理式をもとに、実際に状態方程式と出力方程式を求めたいと思います。

今回は例題として、電圧$E_{mt}$と負荷トルク$T_{mt}$を入力とし、電流$i$を出力とするモデルを作成したいと思います。

カヲル

本記事では線形代数の慣例に従って、小文字で太字の量記号をベクトル、大文字の量記号を行列とするよ。

状態方程式の導出

モータの端子間電圧を$E_{mt}$、電流を$i$、コイルの抵抗を$R$、インダクタンスを$L$、逆起電力定数を$K_E$、トルク定数を$K_T$、慣性モーメントを$J_{rt}$、回転速度を$\omega$、粘性抵抗係数を$D$、負荷トルクを$T_{mt}$、時間を$t$と置くと、DCモータの数理式は以下の通りでした。

$$\begin{eqnarray} &&E_{mt} = Ri + L \frac{di}{dt} + K_E \omega \tag{2}\\ &&J_{rt} \frac{d\omega}{dt} = K_T i \, – \, D \omega \, – \, T_{mt} \tag{3} \end{eqnarray}$$

ここでは、Fig. 4で示したように下式のような状態方程式へ変換します。

$$ \dot{\boldsymbol{x}} = A \boldsymbol{x} + B \boldsymbol{u} \tag{4} $$

入力変数$\boldsymbol{u}$（ベクトル）の成分を$E_{mt}$と$T_{mt}$とすると

$$ \boldsymbol{u} = \begin{bmatrix} E_{mt}\\ T_{mt}\\ \end{bmatrix} \tag{5} $$

また、現代制御理論にて状態変数と呼ばれる内部変数$\boldsymbol{x}$（ベクトル）の成分を$i$と$\omega$とすると

$$ \boldsymbol{x} = \begin{bmatrix} i\\ \omega\\ \end{bmatrix} \tag{6} $$

以上の変数の定義を頭に入れながら、式（２）、（３）を式変形すると

$$\begin{eqnarray} \frac{d}{dt}i &=& \frac{E_{mt} \, – \, Ri \, – \, K_E \omega}{L} \\ &=&\left( -\,\frac{R}{L}i \,-\,\frac{K_E}{L}\omega \right)+\left(\frac{1}{L}E_{mt} \right) \tag{7}\\ \\ \frac{d}{dt}\omega &=& \frac{K_T i \, – \, D \omega \,-\, T_{mt}}{J_{rt}} \\ &=& \left( \frac{K_T}{J_{rt}}i \,-\, \frac{D}{J_{rt}} \omega \right) + \left(-\frac{1}{J_{rt}}T_{mt} \right) \tag{8} \end{eqnarray}$$

式（７）、（８）を行列の形で書くと

$$ \frac{d}{dt} \begin{bmatrix} i\\ \omega\\ \end{bmatrix} = \begin{bmatrix} -\frac{R}{L} & -\frac{K_E}{L}\\ \frac{K_T}{J_{rt}} & -\frac{D}{J_{rt}}\\ \end{bmatrix} \begin{bmatrix} i\\ \omega\\ \end{bmatrix} + \begin{bmatrix} \frac{1}{L} & 0\\ 0 & -\frac{1}{J_{rt}}\\ \end{bmatrix} \begin{bmatrix} E_{mt}\\ T_{mt}\\ \end{bmatrix} \tag{9} $$

よって、状態方程式の行列$A$と$B$は以下のようになります。

$$\begin{eqnarray} A &=& \begin{bmatrix} -\frac{R}{L} & -\frac{K_E}{L}\\ \frac{K_T}{J_{rt}} & -\frac{D}{J_{rt}}\\ \end{bmatrix} \tag{10}\\ \\ B &=& \begin{bmatrix} \frac{1}{L} & 0\\ 0 & -\frac{1}{J_{rt}}\\ \end{bmatrix} \tag{11} \end{eqnarray}$$

出力方程式の導出

下図のような状態変数$\boldsymbol{x}$の成分である電流$i$を出力とする出力方程式を求めます。

$$\begin{eqnarray} \boldsymbol{y} &=& \boldsymbol{c} \, \boldsymbol{x} \tag{12} \\ &=& \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} i\\ \omega\\ \end{bmatrix} \tag{13} \end{eqnarray}$$

よって、出力方程式のベクトル$\boldsymbol{c}$は以下のようになります。

$$ \boldsymbol{c} = \begin{bmatrix} 1 & 0 \end{bmatrix} \tag{14} $$

mファイルの作成

前回の記事で作成したmファイルをベースに行列$A$、$B$とベクトル$\boldsymbol{c}$を追加します。

追加分は13～19行目になります。今回はカタログ（RE-260RA-2295）に記載のある$3\,V$印加時の最大効率点（$1.31 \, mN \cdot m$）の電流を計算するために、10～11行目も変更しました。

clear; clc;

Kt  = 2.54 * 10^-3;     % トルク定数 [N-m/A]
Ke  = 2.88 * 10^-3;     % 逆起電力定数 [V/(rad/s)]
R   = 1.11;             % 抵抗 [ohm]
L   = 1.4 * 10^-4;      % インダクタンス [H]
Jrt = 1.4 * 10^-5;      % 慣性モーメント [kg-m^2]
D   = 4 * 10^-7;        % 粘性抵抗係数 [N-m/(rad/s)]

Emt = 3;                % モータの端子間電圧 [V]
Tmt = 1.31 * 10^-3;     % 負荷トルク [N-m]

% 状態方程式
A   = [-R/L   -Ke/L  ;
       Kt/Jrt -D/Jrt];
B   = [1/L    0      ;
       0      -1/Jrt];
% 出力方程式
c   = [1 0];

sim('MotorModel.slx');