ArduinoとMATLABを利用したDCモータの制御系設計

モデルベースデザイン

ArduinoとMATLAB/Simulinkを連携させてDCモータを速度制御してみる

manao55.hatenablog.com

元の記事になります。

manao55.hatenablog.com

上記事の続きとなります。

実機を使った制御設計において, PIDゲインをいろいろ変化させながら実験を繰り返したり, ジーグラ・ニコルスの限界感度法といった持続振動条件を発見し, ゲイン設計を行うといった手法は, ”動き”を"推定"し"獲得"してゆくことになりますが, 時間がかかってしまうことが多いかと思います。

そのため, 制御プラントの数式モデルを獲得し, それに基づいて設計を行うモデルベース設計が行われることが望ましいと思われます。（システム制御の分野に限らず, 回路設計等広範な分野でも採用可かと）試行錯誤を伴う実験を回避したり, 仕様をみたす制御系設計が可能になったり, 新しい制御理論の実験検証がしやすくなります。製品開発期間の短縮, 問題の早期発見といったメリットもあります。

前回では, システム同定によって, 1次遅れ系制御対象（DCモータの速度制御系）のパラメータを同定し, 数式モデルを獲得しました。そこで今回は, 制御プラントの数式モデルを用いつつ, 制御特性の向上を狙って (なるべく簡単に) フィードバック制御設計を行います。

よろしくお願いします。

極配置法による制御ゲイン設計

ここでは, 閉ループ伝達関数の極を配置することによる「極配置法」によってPIゲインを決めていきます。プラントの伝達関数を

${ \displaystyle P(s) = \frac{K}{Ts+1} }$

PIコントローラの伝達関数を

${ \displaystyle C(s) = K_{P} + \frac{K_{I}}{s} }$

と置くことにより, 下図のフィードバック制御系の目標値 $r$ から出力 $y$ までの閉ループ伝達関数を求めてみます。

閉ループ伝達関数は以下の公式で計算できます。

${ \displaystyle G(s) = \frac{C(s)P(s)}{1+C(s)P(s)} }$

上式分子は”前向き伝達関数”で, $r$ から $y$ へ直線的にたどったときの伝達関数を掛け合わせたものなので $C(s)P(s)$ となり, また, 一巡伝達関数はフィードバックループを一巡した伝達関数を掛け合わせたものなので, これも $C(s)P(s)$ となります。したがって $r$ から $y$ までの伝達関数は,

${ \displaystyle G(s) = \frac{C(s)P(s)}{1+C(s)P(s)} }$

${ \displaystyle =\frac{\frac{K}{Ts+1}(K_{P}+K_{I}/s)}{1+\frac{K}{Ts+1}(K_{P}+K_{I}/s)} }$

${ \displaystyle = \frac{\frac{KK_{P}}{T}s+\frac{KK_{I}}{T}}{s^2+\frac{1+KK_{P}}{T}s+\frac{KK_{I}}{T}} }$

となります。この伝達関数の, 分母式である,

${ \displaystyle s^2+\frac{KK_{P}+1}{T}s+\frac{KK_{I}}{T} = 0 }$

の根が閉ループの極となります。極と $p1$ , $p2$ が一致するようにゲイン設計をします。 $p1$ と $p2$ を根に持つ方程式は,

${ \displaystyle (s-p_1)(s-p_2) = s^2 -(p_1+p_2)s + p_1p_2 = 0 }$

上2式を比較すると, 次の式が得られます

${ \displaystyle K_P = -\frac{(p_1+p_2)T+1}{K}, K_I = \frac{p_1p_2T}{K} }$

$p_1$ と $p_2$ は負の実数や共役複素数を設定します。極と応答の理論背景についての考察, ,ステップ応答は, 逆ラプラス変換により

${ \displaystyle y(t) = \mathcal{L}^{-1}[G(s)u(s)]\ }$

${ \displaystyle =\mathcal{L}^{-1}[G(s)\frac{1}{s}]\ }$

ここで $G(s)$ が極 $a_i$ を持つ場合,

${ \displaystyle G(s)=\frac{n(s)}{(s-a_1)(s-a_2)...(s-a_n)} }$

となり,

${ \displaystyle y(t) = \mathcal{L}^{-1} [\frac{c_0}{s}+\frac{c_1}{s-a_1}+\frac{c_2}{s-a_2}+...+\frac{c_n}{s-a_n}]\ }$

${ \displaystyle = c_0+c_1\mathrm{e}^{a_1t} + c_2\mathrm{e}^{a_2t} + ... +c_n\mathrm{e}^{a_nt} }$

となり, 指数関数の和が応答になることが分かります。また, 極が虚部を持たなければ振動せず, 実部が絶対値の大きい負であれば応答は早く収束するイメージもつくかと思います。

それでは実験をしていきます。また今回使用するDCモータをミニ四駆で用いられるレブチューン2モータを使用しました。ちなみに筆者はサイクロンマグナム†が好きです。サイクロンマグナムに使用するモータを同定し, 制御するという体裁で進めさせていただければと思います。マグナム.jpg 注† : サイクロンマグナム装置外観

前回の記事の手法でシステム同定を行いましたがパラメータは

$K : 1.02$

$T : 0.74$

となり, 波形は,

上図のような波形となりました。（すこしずれてしまいましたが...）では制御実験を行います。

モータの電圧指令値を 1V 初期値, 0.5V 振幅のステップ電圧として, 制御器をPI制御器, 設定する極を-2の重根にした場合で実験を行いました。mファイルでSimulinkファイルを呼び出します。

% velo_pi_mbd.m

% Initialize & load data
close all
clear all
load model_data

% PI design by pole placement
p1 = input('p1 = ');
p2 = p1; % 重根の場合

Kp = -((p1+p2)*T + 1)/K;
Ki = p1*p2*T/K;

% Open simulink model
open_system('velo_pi_mbd_sl');
open_system('velo_pi_mbd_sl/Output');

ts = 1/50;
sim('velo_pi_mbd_sl')

% EOF of velo_pi_mbd.m

次に波形となります。上図で黄色が指令値,青が実機応答で, オレンジがシミュレーションによる数理応答になります。下図は青が指令値と数理応答の誤差, 黄色が実機応答と指令値の誤差になります。

$K_p = 1.92$

$K_i = 2.90$

極指定を $p1=p2=-3$ の重根で設定します。

$K_p = 3.37$

$K_i = 6.52$

極を左半面に寄せるほど応答は早くなります。ちなみに $p1=p2=-4$ に設定すると,

応答が早くなりましたが, オーバーシュートが生じます。PゲインもIゲインも増加しますので, 直感的にゲインを上げ続けると発振していきそうですね。（安定限界は調べていません）

$s$ 領域下における極配置設計, お手軽でなかなか良いと個人的に感じます。今回 D 制御は入れていませんが, 古典制御におけるPID制御則はやはり偉大です。。。さすがこの世界を支えている制御則です。

z 変換とデジタル制御について

ATMega328Pというマイクロプロセッサを用いたArduinoでの実験を行いました。制御器は $s$ 領域における連続時間の下での設計でしたが, 本来マイコンが行っている計算は, 周期的に測定量をサンプリングし, デジタル回路のなかで計算を行ってから制御量を決定するので, その動作は断続的なものとなります。（デジタル制御と呼ばれる所以です）また, 今回のシステムは連続ダイナミクス（DCモータは連続時間において1次遅れ系）と離散ダイナミクス（デジタル論理のこと。コンピュータ, オートマトン）が混在したシステムのなので, ハイブリッド制御系と言ったりすることもできます。厳密にはデジタル制御器を実現し, 離散時間下で性能を評価したほうがよいですね。離散信号を扱うのなら " $z$ 変換" を用いるのが便利です。

$k$ を整数として,

${ \displaystyle x(k)=0 (k<0) }$

を満たす $k$ の関数について,このような $x$ は数列

${ \displaystyle {{x(0),x(1),x(2),...}} }$

を表し, $x(k)$ に対して $z$ の関数 $X(z)$ を

${ \displaystyle X(z)=x(0)+x(1)z^{-1}+x(2)z^{-2}+... }$

${ \displaystyle =\sum_{k=0}^{∞}x(k)z^{-k} }$

で定義します。この（正則）関数を $x$ の $z$ 変換と呼びます。離散的時刻 $k$ に対して $z^{-k}$ を割り当てると解釈するといいでしょう。また数列 $x(k)$ に対応するインパルス列 $x(t)$ を

${ \displaystyle x(t)=\sum_{k=0}^{∞}x(k)\delta(t-kT_{s}) }$

で, サンプリング信号を模擬します。時刻 $kT_s$ に発生する単位インパルス[tex:\delta(t-kT{s})]のラプラス変換は, $e^{-skT{s}}$ であり, したがって $x(t)$ のラプラス変換 $X(s)$ は,

${ \displaystyle X(s)= \mathcal{L}[x(t)]\ =\sum_{k=0}^{∞}x(k)e^{-skT_s} }$

となり, ここで,

${ \displaystyle e^{sT_s} = z }$

と置き換えると. $X(z)$ の右辺を得ることができます。つまり, 「インパルス列 $x(t)$ の $z$ 変換は $x(t)$ のラプラス変換において $e^{sT_s}$ を $z$ で置き換えたもの」と解釈できます。この式から, $z$ が1ステップ（時間 $T_s$ ）だけ信号を進める作用を持ち, $z^{-1}$ は1ステップだけ遅らせる作用をすることが示唆されます。