Pythonで考えるDCモータの制御（10）PWM制御におけるLocked Anti-Phase方式とSign-Magnitude方式について

アイキャッチ画像

はじめに

今回の話題はDCモータをHブリッジにより駆動する際の

Locked Anti-Phase PWM方式と

Sign-Magnitude PWM方式の

それぞれがどの様なもので

どの様な違いがあるのかを説明させて頂きます。

モータ制御に関連して僕が実際に作ったロボットの話題を一つご紹介します。

既出のマイクロマウス「３式改」ですが、モータドライバは自作でした。

CPLDとゲートドライバとFETを組み合わせて

Locked Anti-Phase方式で制御しておりました。

SH2マイコンとCPLDを８ビットバスでつなげて

CW/CCWに２５６段階のPWM制御が可能な仕様でした。

CPLDの役割はマイコンから指令されたDuty幅でPWM信号をつくり、生成した信号はゲートドライバを経てFETをスイッチングします。

その際、ハイサイドFETとローサイドFETの短絡を防ぐための

デッドタイムを作り出します。

マイクロマウス３式改のモータドライバの写真 — マイクロマウス３式改のモータドライバ

３式改については以下でも触れています。

グラフ理論と迷路探索 - 理系的な戯れ

オドメトリによる移動ロボットの自己位置推定 - 理系的な戯れ

当時のモータ制御の参考は以下のページでした。ランサーロボットのモータドライバも自作したのでそのときからですので２０年ほど前になります。何回も読んでモータドライバを作りました。

NOTE - 回路設計編 - PWMドライブの設計

NOTE - 回路設計編 - モータドライバと象限

DCモータ駆動方式に上記のページで紹介されている以下の二つの方式があります

Locked Anti-Phase PWM方式　（LAP方式）
Sign-Magnitude(BREAK) PWM方式　（SMB方式）

これらのうちLAP方式がリニアな制御理論に一致する制御結果が現れると言うことでそちらを採用してモータドライバを動かすことにしました。

今回は、LAP方式とSMB方式の挙動の違いに関してシミュレーションを通じて概観してみたいと思っています。

はじめに
LAP方式とSMB方式
- Locked Anti-Phase PWM方式
- Sign-Magnitude(Break) PWM方式
Hブリッジの動作
シミュレーションの準備
計算結果と考察
おわりに
- 長文で申し訳ありません
- ブログの間隔が開きがちですが・・・
参考文献
付録（プログラム）

LAP方式とSMB方式

Locked Anti-Phase PWM方式

LAP方式はモータに印加する電圧の極性が常に入れ替わる方式です。常に正転と逆転の指令を入れ替えています。正転指令の時間が逆転指令の時間より長ければ正転しますし、短ければ逆転します。 PWM制御ではパルスの幅をDutyと呼びパーセント（％）で表すことが多いのですが、Duty１００％が正転の最大値、Duty０％が逆転の最大値になりDuty５０％が停止となります。

Sign-Magnitude(Break) PWM方式

SMB方式はPWM１周期の中はオンとブレーキの切り替えだけが起こります。オンの時の極性はその都度変更されますが、１周期の間に極性が変わると言うことはありません。

もしかしたら、普通にロボットを動かすとこちらを採用していることが多いかもしれません。

モータに付加する電圧の与え方で両者を比較してみたのが次の図です。電圧の正負はモータに印加される電圧の極性を示し、実際にマイナス電源がいるわけではありません。

以下の図はそれぞれの方式のモータにかかる電圧を表したものです。

Hブリッジの動作

モータ駆動の正転・逆転駆動回路の基本はHブリッジ回路というもので４つのスイッッチ（下図ではNchFET）を下図の様に配置します。 FETとモータでHの字を形成している様に見えるのでHブリッジ回路と呼ばれています。

４つのスイッチのON・OFFのしかたで以下に示す４通りのモードがあります。

Hブリッジのスイッチ状態による４つのモード図 — Hブリッジのスイッチ状態による４つのモード

①はストップモードでFETに内臓のダイオードでダイオードがONならば電流がGNDから電源方向に流れます
②は正転です正転逆転は相対的な基準です。今回は図のようにモータに電流が左から右に流れている状態です
③は逆転です。正転の反対ですね。
④はハイサイドのスイッチが全てOFF、ローサイドのスイッチが全てONでローサイドで電流が還流できるようにしています。

デッドタイムについて

PWM駆動でモータを制御する場合に理屈通りにON・OFFをすると、ON・OFFのディレイのためにハイサイドとローサイドのスイッチが同時にON状態になりFETが発熱し危険な状態になります。

そこで、PWM駆動回路においてはON・OFFの遷移の際にハイサイドとローサイドが両方ともOFFの時間を設けて両方がONになることを防ぎます。この時間をデッドタイムと呼びます。

デッドタイムについてはFETが確実にOFFしきるまでの時間以上の時間を設ける必要があります。 FETのONとOFFについては様々なところで解説されていますが、下図にあります様にゲートにあるとみなせるコンデンサの充電と放電をする必要があります。

さらに放電時間に加えてそれぞれのFET特性のターンオフ時間が加わり確実にOFFしきるまでの時間が算出されます。これに安全率を考慮してデッドタイムを決めるのが流れだと考えます。

ゲートドライブに関して以下の資料は大変勉強になります。

MOSFET ゲート駆動回路

LAP方式のシーケンス

LAP方式は一周期のうちに正転と逆転を切り替えます。その間にデッドタイム（①ストップモード）を挟みます。よって１周期駆動シーケンスは

②正転モード→①ストップモード→③逆転モード→①ストップモード

になります

SMB方式のシーケンス

SMB方式は１周期の間は正転ならばプラス電圧、逆転ならばマイナス電圧（実際は極性を入れ替えるだけ）が印加されます。よってシーケンスは正転と逆転で場合わけされて

正転の場合

②正転モード→①ストップモード→④ブレーキモード→①ストップモード

逆転の場合

③逆転モード→①ストップモード→④ブレーキモード→①ストップモード

ストップモードは電源をそれまでとは逆に接続した状態になったと見做せますので、電流は鯉の滝登り状態となり急激に減衰します。この状態をモータドライバのデータシートのよってはFast Decayと解説しているものもあります。これに対して④ブレーキモードは電源から切り離されますので電流は緩やかに減少しSlow Decayなどと説明されているものもあります。

シミュレーションの準備

モータの回転の運動方程式は各モードで共通なので、各モードの回路方程式を立てて前述したシーケンス順にモードに応じた回路方程式に切り替えてルンゲ・クッタ法で計算していこうと思います。

①ストップモードの方程式（デッドタイム）

このモードが一番難しいです。ただしく表されているか自信はありません。

このモードに入ってくるダイオードのモデルをどの様に置こうかと悩んでいます。

~~動的なモデルはまだ理解できていないので、静的モデルでとりあえず作ってみる予定です。~~ ~~静的モデルについては、折れ線近似にしてみます。~~

調べた結果、ダイオードの静的モデルは以下のような式で表すことができるようです。

$\begin{eqnarray} i = I_0 \left(e^{\frac{v_d}{n v_t}}-1 \right)\\ \end{eqnarray}$

これは電圧 $v_d$ に対して電流 $i$ が求められる式ですが、回路方程式はこの逆関数のほうが都合が良いので、両辺の対数をとって逆関数を求めます。s

$\begin{eqnarray} v_d =n v_t log \left(\frac{i+I_0}{I_0}\right)\\ \end{eqnarray}$

更には、Hブリッジの場合両ブリッジにダイオードが存在するので全てのスイッチがOFFになってもモータの電流は両方向に流れることが可能です。ただし、電圧に対して逆向きに流れようとするので、長時間流れるといずれ０になります。

そこで、今回のシミュレーション用に作ったダイオードモデル式は以下になります。参考文献の「モータ制御」においても同じようなアプローチが書かれているのでまんざらでも無いのでは？と思っています。

$\begin{eqnarray} v_d(i) =sign (i) n v_t log \left(\frac{|i|+I_0}{I_0}\right)\\ \end{eqnarray}$

$sign(x)$ と言うのは、 $x$ の符号と言う意味です。

数式で言ってる事をグラフ化すると以下になります。

この場合、ダイオードに逆電流が流れ始めると、しばらく逆電流がながれてしまう現象があるのですが、それを再現できないので、後日の課題とします。

ようやく、ダイオードのモデルを適用して、ストップモードの回路方程式を立ててみます。他のモードの式もそうですが、モータが正転するトルクが発生する電流の向きが正としています。

ダイオードの工夫については前述しましたが、左右のブリッジのダイオードのお陰でモータには両方向に電流が流れる事が可能ですが、どちらの向きに流れても常に電圧の向きに対しては逆らって流れるはずなので、電源電圧が流れの向きに対して反対になるようなトリックを組み込んでいます。

以上から、①ストップモードの回路方程式は以下となります。

$\begin{eqnarray} L \frac{d i}{d t} + R i + K \omega + v_d(i) = sign(-i) u\\ \end{eqnarray}$

②正転モードの回路方程式

正の極性で電源電圧が印加されると考えて、以下の式になります。

$\begin{eqnarray} L \frac{d i}{d t} + R i + K \omega = u\\ \end{eqnarray}$

③逆転モードの回路方程式

負の極性で電源電圧が印加されると考えて、以下の式になります。

$\begin{eqnarray} L \frac{d i}{d t} + R i + K \omega = -u\\ \end{eqnarray}$

④ブレーキモードの回路方程式

閉回路に電源はないので閉回路が一周した場合の電位差は０と考えられるので。

$\begin{eqnarray} L \frac{d i}{d t} + R i + K \omega = 0\\ \end{eqnarray}$

モータの回転の運動方程式

これまでこのシーリーズでは何回も出てきましたが、モータの回転の運動方程式は以下です。

$\begin{eqnarray} J \frac{d \omega}{d t} + D \omega = K i\\ \end{eqnarray}$

各種記号

記号	意味	単位
i	電流	A
$\omega$	角速度	rad/s
L	インダクタンス	H
R	コイル抵抗	$\Omega$
u	電源電圧	V
K	トルク定数・逆紀電圧定数	Nm/A・Vs/rad
J	慣性モーメント	kgm²
D	動粘性係数	Nms/rad
$v_d$	ダイオード電圧	V
$v_t$	熱電圧	V
n	補正値	1

計算結果と考察

1周期分の波形からデッドタイムの影響を見る

まず最初は、PWM周波数１MHｚの場合の１周期分の計算結果を示します。 PWMのオン時間６００ｎｓ、デッドタイムが１００ｎｓ、オフ時間が２００ｎｓ合計で周期１μｓになります。

上のグラフは角速度、下のグラフが電流です。以下同様になります。

続いて、逆にＯＦＦ時間が短い場合です。オン時間が２００ｎｓ、オフ時間が６００ｎｓ、デッドタイムは同じ１００ｎｓの波形を示します。

デッドタイムの影響がＯＦＦ時間が長いほうで判りやすく見られますが、デッドタイムは流れている電流を減衰させ０に近づけていきますが０を超えて符号が反転することはありません。理想的にはデッドタイムが入らず１周期の電流波形はきれいな「へ」の字であることが期待されますが、ＯＦＦ時間の長いほうの計算結果では明らかに「へ」の字ではなく、波形の最後が反転しており「Ｚ」字になっています。角速度にもデッドタイムの影響があるのが確認できます。

これは、波形全体の時間、すなわち１周期の長さが１μｓに対してデッドタイムが１００ｎｓと１割の時間を占めているため顕著にみられるものです。デッドタイムは上下のFETの短絡を防ぐためには必要ですが、 PWM制御を理想的に行うためには悪影響を及ぼします。

試みに、デッドタイムを全体の１％にあたる１０ｎｓに設定して、ＯＦＦ時間が長い方の計算をやり直してみます。この際、ＯＦＦ時間は同じ６００ｎｓでデッドタイムを変えた分はＯＮ時間を３８０ｎｓにして全体周期は同じにしています。

こちらは相対的にはデッドタイムの影響は軽減できていると考えます。

したがって、デッドタイムは全体の周期の中でできるだけ小さくしたほうがよく、デッドタイムが無い場合の理想のＰＷＭに近づくと言う当然の結論を得ることができました。

しかし、FETの性能を考えた場合、デッドタイム１００ｎｓは現実的ではなく２μｓや３μｓぐらいにしなければならず PWM周期をもっと大きくしなければなりません。

ところが、仮にデッドタイムを2μsとしてそれをPWM周期の１％に相当するとすると、周期は２００μsです。１７１７モータの電気的時定数は16μsなのでそのおおよそ２０倍ぐらいの長さになり長すぎます。

以下は、１７１７モータをLAP方式でこの値で動かした場合のシミュレーションの結果です。激しく電流が上下するのが確認できます。

チョークコイルを加えて電気的時定数を大きくする

問題があるかもしれないのですが、以後、暫定的にデッドタイムを２μsとしPWM周期の１％にすることとします。つまりPWM周期は２００μsということでPWM周期は５KHzと言うことになり、自分的にはPWM制御をするためには５０KHzや１００KHzをイメージしていたので、それに比べれば遥かに遅く違和感がありますがこの件の検討は後日行うものとし、これで話を進めます。

この場合チョークコイルを接続してインダクタンスを大きくすることを想定しています。ここでは仮にインダクタンスが５００μHになったと仮定します。

この場合で電気的時定数が４６７μsになります。

コイルの例：村田製作所のコイル

では５KHzの計算結果です。

リップルが大きく感じられます。インダクタンスを大きくしても周期２００μsはまだ長いのだと思います。角速度の波形はステップ応答の波形に近く線形解析の結果とほぼ同じです。ただし、電流波形に顕著ですが０．０３秒から０．０５秒ぐらいに不自然な段差があります。これはデッドタイムの影響で応答が非線形になっているものと思われます。のちほど考察します。

こうなると、リップルも小さくなるはずですし１０KHzにあげてデッドタイムは全体の２％に相当する様にしてみてもうまくいくかもしれませんので、計算してみました。

０．０４秒から０．１秒にかけての波形に違和感を感じます。やはりデッドタイム２％の影響は大きい様です。今回はPWM周期２００μs（周波数５KHz）で検討を進めたいと思います。

電流が０に近づいた時のデッドタイムの影響

デッドタイムの長短の影響について考えてみたいと思います。

おそらくDC/DCコンバータの同期整流と非同期整流の違いの様なことがモータのPWM制御にも起こるのだと考えています。

同期整流型と非同期整流型の違い

④ブレーキモードはDC/DCコンバータの同期整流しているイメージですが、①ストップモードはDC/DCコンバータがダイオードだけで非同期整流しているイメージです。

FETがオフになってもダイオードを介して電流は流れるのですが徐々に減衰して０Aに達すると０Aに張り付いてしまいます。そうすると、PWM制御はリップル電流（上下に上下する電流）の平均の電流が流れているとみなせるので０Aに張り付いてしまうと例えば電流がプラスの領域では電流が減衰して０Aに近づくにつれて少しだけ大きく電流が流れていることになります。

計算をしてみるとデッドタイムが大きい場合は以下の様な応答を示します。特にリップル電流波形の下端が奇妙にクリップされているのが確認できます。

０．０３６秒付近の電流は０Aにクリップしていますが、拡大してみると。

０．０７５秒から０．１秒の間も違う感じになっていますが、０．０８秒付近の電流波形を拡大してみました。

０．１１秒から最後までの様子もそれまでと違います。拡大してみると

デッドタイムはON時間の最後とOFF時間の最後に入りますが、どのタイミングでデッドタイムに入るかで上の３パターンのような波形を描きます。この様子が、平均電流の波形に歪みを生じさせることがわかりました。

シミュレーションでしかできないことですが、全くデッドタイムを入れないで計算すると次の様になります。

LAP方式とSMB方式の比較

長くなりましたがLAP方式とSMB方式の比較をしてみたいと思います。

単にそれぞれの波形を見ただけでは何がどう違うのかが判定できなかったので、リニア制御と比較してみることにしました。全ての方式に速度目標値が５００rad/sのPI制御をかけて比較しました。ゲインの設定も３方式とも同じものを使いました。

LAP方式がリニア制御に近い特性を持つと言う話を聞いたことがあるのですが、この結果をみるとかなり近いと感じました。ほぼ傾向は一緒です。前半のすべての制御方式で制御入力が飽和してしまっている部分がLAP方式とSMB方式がほぼ同じ波形でリニア制御が飽和部分は少し早く終わるところが少々疑問だったのですが、どうやらここでもデッドタイムが影響しているようです．試しに現実を無視してLAP方式のデッドタイムを小さくしてみた結果が以下です。

ほとんどリニア制御に重なっていて、LAP方式についてはリニア制御のゲイン調整で大体の当たりがつけられそうな感じがしてきました。SMB方式についてはデッドタイムを小さくしても、リニア制御と同じ波形とみなせる様にはなりませんでした。

おわりに

長文で申し訳ありません

長い記事になりました。

ダイオードのモデルで悩んでましたが、FETのモデルは完全に無視していたりして所々方向性が定まっていない感がありますが、ダイオードに関しては計算して大体の傾向が見えて来たのではと思っています。

慣性モーメントに押し込めて負荷モデルが考慮されていないので今後負荷が加わった場合の検討もできたらと思っています。

シミュレーションをやりながら、何をやっているのか判らなくなる場面もありましたが、一晩寝て朝になると執筆が進む場面もありました。計算しては執筆し計算しては執筆し、また少し戻って書き直すといった感じで書きあがりました。

昔からLAP方式が好きで、自分でもモータドライバを自作して実装しているのですが、満足のいくまでできていません。今回の検討でデッドタイムをおろそかにできないと何となく解ってきたので次に自分で作るときに設計に生かせるかと思っています。

LAP方式が好きな理由として制御法則がそのまま適用できると言う噂がありましたが、それについては真実と考えて良いようです。

LAP方式とSMB方式の比較と優劣ですが、それが判るほどにはいったってないかなと思いますが、SMB方式は波形があまり美しくない感じがします。実機では何か悪いことが起こりそうな気がします。確か、実機でやると、速度０付近で思ったのと違う挙動をする記憶があります。

１７１７モータのようなインダクタンスの小さなモータはLAP方式はかなり厳しいかもしれません。FETの性能をよく調べて、チョークコイルの付加等を考慮して設計パラメータを決めなければなりませんね。

やはり昔、モータドライバをCPLD＋ゲートドライバ＋FETで作ったように、また自作したいものです。まだよくわかってませんがGANを使ったら性能上がるのでしょうか？

ブログの間隔が開きがちですが・・・

最近Twitterが面白くなってきて、そちらで遊んでいたのでブログが疎かになっていますが、ブログも頑張っていこうと思います。

僕のペースでは３週間に２本ぐらい更新できれば今のところ精一杯なのかなと思います。頑張れる時は頑張ろうと思いますが、まずは無理しないで継続を目指します。ネタについては、何か考えているとわからないことや試してみたいこと、調べたいことがたくさん出てくるので、尽きないとは思うのですが、車輪の再発明とか世の中で解決済のことに延々と悩むと言うこともあるかもしれません。

では、今回はこれまでにしたいと思います。長い記事をここまで読んで頂きましてありがとうございます。

では、また！

参考文献

DCモータ制御に関する２０年頼の座右の書

実験とシミュレーションで学ぶモータ制御―Visual Basicコード付き

作者:尚志, 見城,久, 高橋,美奈, 寺内,敦, 藤田,茂, 佐渡友,雅俊, 久保,雅行, 高田
日刊工業新聞社

Amazon

以下のページは大変参考になります。もうすこし読み込んで本記事にも反映させたいところです。

Lock Anti-Phase Drive | Modular Circuits Sign-Magnitude Drive | Modular Circuits

日本語に翻訳された方のブログもありました

Sign-Magnitude Drive (Modular Circuits和訳記事) - まんじゅう日記

記事の中で紹介したゲートドライブに関しての東芝さんの資料

MOSFET ゲート駆動回路

記事の中で紹介したデコデコの同期生流と非同期清流の解説

同期整流型と非同期整流型の違い

付録（プログラム）

ルンゲクッタソルバー＆導関数等

まずは，ルンゲクッタ関係の関数の定義です。他の二つの前にこれをコピーして使ってください。

#%matplotlib inline
#%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import tqdm as tqdm

'''
def rk4(func, t, h, y, *x)
ルンゲ・クッタ法を一回分計算する関数
    引数リスト
    func:導関数
    t：現在時刻を表す変数
    h：刻み幅
    y：出力変数（求めたい値）
    *x:引数の数が可変する事に対応する、その他の必要変数
※この関数では時刻は更新されないため、これとは別に時間更新をする必要があります。
'''
def rk4(func, t, h, y, *x):
    #print(t,h,y, *x, y)
    k1=h*func(t, y, *x)
    k2=h*func(t+0.5*h, y+0.5*k1, *x)
    k3=h*func(t+0.5*h, y+0.5*k2, *x) 
    k4=h*func(t+h, y+k3, *x)
    y=y+(k1 + 2*k2 + 2*k3 + k4)/6

    return y

'''
導関数の書き方
def func(t, y, *state):
    func:自分で好きな関数名をつけられます
    t:時刻変数(変数の文字はtで無くても良い) 
    y:出力変数(変数の文字はyで無くても良い)
    *state:その他の必要変数(引数の数は可変可能))
#関数サンプル
def vdot(t, y, *state):
    s1=state[0]
    s2=state[1]
    return t+y+s1+s2
    
'''
#以下ロボットの位置と速度を計算するためルンゲクッタソルバに渡す導関数

def motor_spec():
    #1717モータ
    J=0.59e-7
    L=500e-6 #チョークコイルでインダクタンスを大きくした想定の値
    #L=17e-6 #本体の1717モータのインダクタンス
    R=1.07
    K=1.98e-3
    D=2.36e-8
    tm=J*R/(D*R+K**2)
    te=L/R
    return J, L, R, K, D, tm, te

def v_d(i):
    n=1
    vt=26e-3
    Is=1e-14
    return np.log((i+Is)/Is)*n*vt

def idot_dead(t, i, omega, u):
    J,L,R,K,D,tm,te=motor_spec()
    return (np.sign(-i)*u -R*i -K*omega -np.sign(i)*v_d(np.abs(i)))/L

def idot_cw(t, i, omega, u):
    J,L,R,K,D,tm,te=motor_spec()
    return (u -R*i -K*omega)/L

def idot_ccw(t, i, omega, u):
    J,L,R,K,D,tm,te=motor_spec()
    return (-u -R*i -K*omega)/L

def idot_break(t, i, omega, u):
    J,L,R,K,D,tm,te=motor_spec()
    return (-R*i -K*omega)/L

def omegadot(t, omega, i):
    J,L,R,K,D,tm,te=motor_spec()
    return (K*i - D*omega)/J

LAP方式　PI制御

#LAP1-with-deadtime control
OMEGA=[]
I=[]
T=[]
omega=0
i=0
Power =3.0
t=0.0

#電気的時定数　1.5887850467289717e-05
#刻み幅
h=1e-7
N=100

Ref=500
err=0
kp=1
ki=1
Integ=0

Sampling=1
Duty=1
DeadTime=20
Period =2000

for n in  tqdm(range(Period*N+1)):

    err=Ref - omega
    Integ+=h*err
    control=kp*err + ki*Integ 
    if control>=Power:
        control=Power
    if control<-Power:
        control=-Power

    Duty = Period/2 + Period*control/Power/2

    if Duty<DeadTime:
        Duty=DeadTime
    if Duty>Period-DeadTime:
        Duty=Period-DeadTime
        
    Ton=int(Duty-DeadTime)
    Toff=int((2000-Duty)-DeadTime)
    
    if n%Sampling==0:
        OMEGA.append(omega)
        I.append(i)
        T.append(t)
    oldomega=omega
    oldi=i

    if n%Period<Ton:
        i= rk4(idot_cw, t, h, oldi, oldomega, Power)
    elif n%Period<Ton+DeadTime:
        i= rk4(idot_dead, t, h, oldi, oldomega, Power)
    elif n%Period<Ton+DeadTime+Toff:
        i= rk4(idot_ccw, t, h, oldi, oldomega, Power)
    else:
        i= rk4(idot_dead, t, h, oldi, oldomega, Power)
        
    omega = rk4(omegadot, t, h, oldomega, oldi)
    t=t+h
    #print('i omega',i,omega,e,TL)

plt.figure(figsize=(10,7))    
plt.subplot(211)
plt.plot(T, OMEGA, label='LAP')
#plt.plot(T, OMEGA, label='SMB')
plt.xlabel('Time(s)')
plt.ylabel('Omega(rad/s)')
plt.legend( loc='lower right')
plt.grid()

plt.subplot(212)
#plt.figure(figsize=(15,5)) 
plt.plot(T, I, label='LAP')
#plt.plot(T, I, label='SMB')
plt.xlabel('Time(s)')
plt.ylabel('Current(A)')
#plt.ylim(-1,1)
plt.legend( loc='upper right')
plt.grid()

T1=T
OMEGA1=OMEGA
I1=I

SMB方式　PI制御

#SMB-with-deadtime control
OMEGA=[]
I=[]
T=[]
omega=0
i=0
Power =3.0
t=0.0

#電気的時定数　1.5887850467289717e-05
#刻み幅
h=1e-7
N=100

Ref=500
err=0
kp=1
ki=1
Integ=0

Sampling=1
Duty=1
DeadTime=20
Period =2000

for n in  tqdm(range(Period*N+1)):
    err=Ref - omega
    Integ+=h*err
    control=kp*err + ki*Integ 
    if control>=Power:
        control=Power
    if control<-Power:
        control=-Power

    Duty = Period/2 + Period*control/Power/2

    if Duty<DeadTime:
        Duty=DeadTime
    if Duty>Period-DeadTime:
        Duty=Period-DeadTime
        
    Ton=int(Duty-DeadTime)
    Toff=int((2000-Duty)-DeadTime)
    
    if n%Sampling==0:
        OMEGA.append(omega)
        I.append(i)
        T.append(t)
    oldomega=omega
    oldi=i

    if n%Period<Ton:
        if control>=0:
            i= rk4(idot_cw, t, h, oldi, oldomega, Power)
        else:
            i= rk4(idot_ccw, t, h, oldi, oldomega, Power)
            #print('hoge')
    elif n%Period<Ton+DeadTime:
        i= rk4(idot_dead, t, h, oldi, oldomega, Power)
    elif n%Period<Ton+DeadTime+Toff:
        i= rk4(idot_break, t, h, oldi, oldomega, Power)
    else:
        i= rk4(idot_dead, t, h, oldi, oldomega, Power)
        
    omega = rk4(omegadot, t, h, oldomega, oldi)
    t=t+h
    #print('i omega',i,omega,e,TL)

plt.figure(figsize=(10,7))    
plt.subplot(211)
plt.plot(T, OMEGA, label='LAP')
#plt.plot(T, OMEGA, label='SMB')
plt.xlabel('Time(s)')
plt.ylabel('Omega(rad/s)')
plt.legend( loc='upper left')
plt.grid()

plt.subplot(212)
#plt.figure(figsize=(15,5)) 
plt.plot(T, I, label='LAP')
#plt.plot(T, I, label='SMB')
plt.xlabel('Time(s)')
plt.ylabel('Current(A)')
#plt.ylim(-1,1)
plt.legend( loc='upper right')
plt.grid()

T2=T
OMEGA2=OMEGA
I2=I

リニア制御　PI制御

PWM制御と比較するリニア制御版です。

#liner control
OMEGA=[]
I=[]
T=[]
omega=0
i=0
Power =3.0
t=0.0

#電気的時定数　1.5887850467289717e-05
#刻み幅
h=1e-7
N=100

Ref=500
err=0
kp=1
ki=1
Integ=0

Sampling=1
Duty=1
DeadTime=20
Period =2000

for n in  tqdm(range(Period*N+1)):
    
    err=Ref - omega
    Integ+=h*err
    control=kp*err + ki*Integ 
    if control>=3:
        control=3
    if control<-3:
        control=-3
        
    if n%Sampling==0:
        OMEGA.append(omega)
        I.append(i)
        T.append(t)
    oldomega=omega
    oldi=i

    i= rk4(idot_cw, t, h, oldi, oldomega, control)
    omega = rk4(omegadot, t, h, oldomega, oldi)
    t=t+h
    #print('i omega',i,omega,e,TL)

plt.figure(figsize=(10,7))    
plt.subplot(211)
plt.plot(T1, OMEGA1, label='LAP')
plt.plot(T2, OMEGA2, label='SMB')
plt.plot(T, OMEGA, label='Liner')
plt.xlabel('Time(s)')
plt.ylabel('Omega(rad/s)')
plt.legend( loc='upper left')
plt.grid()

plt.subplot(212)
#plt.figure(figsize=(15,5)) 
plt.plot(T1, I1, label='LAP')
plt.plot(T2, I2, label='SMB')
plt.plot(T, I, label='Liner')
plt.xlabel('Time(s)')
plt.ylabel('Current(A)')
#plt.ylim(-1,1)
plt.legend( loc='upper right')
plt.grid()