VBAによるモータベクトル制御シミュレーション
モータベクトル制御のシミュレーションとしてはPSIMやLTspice,simulinkなどGUIを使ったものがありますが、これらは結局C言語ソフトウェアでも実施でき、それならばエクセルに付属のVBA（エクセルマクロ）でもできるのではないかと思い、作成してみました。
VBAにてシミュレーションする利点は、ソフトのインストールなど必要なく、エクセルさえあれば簡単に実行できることです。
例えばソフト作成者はベクトル制御の各ステップの計算値を確認したり、グラフ化することでベクトル制御への理解を深めることができるでしょう。
今回は条件として、
Id=0 Iq:電流指令での電流ベクトル制御としました。
マクロ込みのエクセルファイルはここからダウンロードできます。
(右クリックし、新しいタブで開くで開き、ダウンロードする。）
0．使い方
事前準備としてエクセルソフトを開き、ファイル⏎その他/オプション⏎リボンのユーザ設定ウインドウの「開発」にチェックを入れる。
ダウンロードしたマクロ付きエクセルファイルを開き、「開発」タブを開き、
「マクロ」ボタンでマクロウインドウを開き、ID_0_ベクトル制御を実行します。
結果がデータシートに表示されますので、範囲を選択する形でグラフを表示します。
1.ダイアグラム
電流指令をId=0、Iq=iref として入力し、UVW相のモータ電流を検出したiu1,iv1,iｗ1をdq変換し、フィードバックしています。d軸サーボコントローラ、q軸サーボコントローラではPI制御を行い、その出力を逆dq変換し3相モータ出力電圧u1,v1,w1として出力します。
インバータ実物は効率の面でPWM出力ですが、シミュレーションでは制御を簡単にするため省略しています。モータ電流をシミュレートするに必要なモータ逆起電力とdq変換、dq逆変換に必要な位相（theta)はモータモデルから速度の積分として得られるようにしています。
モータモデルはモータに発生するトルクはId=0なので磁石トルクのみと考えられるので、T=KT*Iとしてシミュレートしています。（KTはトルク定数）
２．コード
〇変数、定数の宣言
VBAでは変数、定数を宣言しなくても使用できますが、
安全のため宣言します。変数の型はsingle(1.4×10(-45乗）～1.8×10(38乗）とします。
〇dq変換
後に説明するモータの三相電流（iu1,iv1,iw1）をdq変換（相対変換）します。
まず、3相電流をα、β軸の直交2軸に変換します。
α軸電流をa1とすると、a1=iu1、β軸電流をb1とすると、b1=1/√3（iv1-iw1）
その後、αβ軸からdq軸（ロータの磁界方向をd軸、それに直交する方向をq軸とする。）に変換します。
d軸電流をIdとすると、Id＝a1 ×Cos(theta) + b1 × Sin(theta)
q軸電流をIqとすると、Iq = -1×a1 × Sin(theta) + b1 × Cos(theta)
但し、thetaはモータ位相。
Id=0、Iqは可変の時、uvwの三相電流からαβ軸変換波形は下図のようになります。
iu1とa1は重なります。b1はa1から90度遅れになります。
αβ軸からdq軸への変換は下図の波形になります。
Iqはa1,b1の正側ピークとなります。これはこのdq変換を相対変換で行っているためです。dq変換には係数が絶対変換のdq変換もあるので注意が必要です。
相対変換と絶対変換では後で説明するモータパラメータの値を方式に応じて変換する必要があります。
〇piサーボコントローラ
次にd軸、q軸それぞれに対してPI制御の電流フィードバックをかけます。
＿q軸の電流フィードバック
まず、q軸の指令電流cq1をシミュレーション初期から入力すると急峻な変化によりシミュレーションが破綻することがあるので、指令値irefを1次フィルタをかけて入力します。
cq1 = rungekut(iref, cq1, tau2, dt)
ここで、1次フィルタに4次ルンゲクッタフィルタのサブルーチンを使用します。
次に指令電流とフィードバック値Iqから差を算出します。
cq2 = cq1 - Iq
これをpi制御サブルーチンに代入する。
cq2pid = pid1(cq2, kp, ki, dt) 'c2pid
この出力にゲインAを掛ける。
cq3 = cq2pid * A
これをリミット関数に代入し、駆動電圧Vqを算出する。
インバータはハードの都合で出力電圧は電源電圧より大きくなることがないので
リミット関数を挿入する。
Vq = limit(cq3, Vmax, Vmin) 'q軸出力電圧
d軸側も指令電流Id=0として同様にコーディングする。
〇逆dq変換
サーボコントローラ出力のVq、Vdに位相信号thetaを合わせて、逆dq変換を行い
電圧出力u1,v1,w1を算出する。コードを下記に示す。
Va = Vd * Cos(theta) - Vq * Sin(theta)
Vb = Vd * Sin(theta) + Vq * Cos(theta)
'α、β→UVW変換
u1 = Va
v1 = -1 * Va / 2 + Sqr(3) / 2 * Vb
w1 = -1 * Va / 2 - Sqr(3) / 2 * Vb
波形は下図のようになります。u1とVaは重なっています。
本来はこの後でPWM変換があるのですが、このシミュレーションでは省略します。
〇モータモデル
インバータ三相出力電圧u1,v1,w1が算出できましたが、一方Vqとモータパラメータを
モータモデルに入力し、モータの出力トルクを求め、そこからモータ速度、モータ位相、モータ逆起電圧を求めて行きます。
モータモデルのダイアグラムは上図に示したようになり、入力電圧Vqから逆起電圧m8を引いたものを電気的時定数R＋Lを通して電流に変換し、それにトルク定数KTをかけてトルクを算出し、負荷トルクTloadを引いた後で、機械的時定数の慣性モーメントJと粘性摩擦係数からなる時定数を通して速度m7が求められる。
速度m7にモータ逆起定数KEをかけると逆起電圧m8を算出でき、これをフィードバックします。
位相thetaは速度m7×ステップ時間dtの積分で算出できます。
theta=∫m7dt
コードではtheta = theta + m7 *dtとなります。
尚、負荷トルクの印加は急峻な変化を避けるため、下記コードにてスロープを付けている。
If t < 1000 Then
Tload = TloadF * 0.001 * t
Else: Tload = TloadF
End If
モータに発生するトルクm4は今回のId＝０制御では磁気トルクのみとなるので、単純に
コイル電流m3とトルク定数KTから求められます。
m4=KT×m3
モータモデルからモータ速度m7、モータ位相m8、モータ逆起電圧thetaを算出するのは下記のコードとなります。
m1 = Vq - m8 '出力電圧ー誘起電圧
m2 = m1 / Ra
m3 = rungekut(m2, m3, taue, dt) '４次元ルンゲクッタ法 m3モータ電流
m4 = m3 * KT 'm4モータトルク
m5 = m4 - Tload 'モータトルクー負荷トルク
m6 = m5 / de 'モータトルク/粘性制動係数
m7 = rungekut(m6, m7, tauj, dt) 'm74次元ルンゲクッタ法フィルタ 角速度
m8 = m7 / pole * KE 'm8誘起電圧
theta = theta + m7 *dt 'theta位相
〇モータ電流検出
最初に示したモータの三相電流(iu1,iv2,iw1)はU,V,W相それぞれのモータ出力電圧からその相の逆起電圧を引き、それにモータパラメータのR＋Lを時定数をとおして、各相のモータ電流が算出されます。
U相のモータ電流iu1を算出するコードは下記になります。
du1 = u1 - m8 * Sin(theta + pi) 'コイルにかかる電圧
iu2 = du1 / Ra
iu3 = rungekut(iu2, iu3, taue, dt) '４次元ルンゲクッタ法 La,Raによるフィルタ
iu1=iu3’結果をiu1に代入する。
〇電源電圧について
実際のハードではインバータの電源は片電圧ですが、シミュレーション上は
両電源としています。
両電源とすることで、モータ中点電圧＝０Vとしてシミュレーションすることができ、オフセット電圧を考慮しないで済みます。
〇シミュレーションステップについて
本シミュレーションは時間軸を基にしたシミュレーションですが、ステップ時間dt
を一定にすると位相急変時に位相変化量が1ステップで0.1°以上になることがあり、そのことでシミュレーションが発振する可能性があるため下記コードにより制限しています。
If m7 < 10 And m7 > -10 Then
carF = 600
Else: carF = m7 * 60
End If
dt = 1 / carF
〇モータパラメータ
本シミュレーションに用いるモータパラメータはＹ結線の1相分を基準とします。
トルク定数：発生トルク/コイル電流（相電流）の0-peak値
コイルインダクタンス：Y結線の1相分の値
コイル抵抗：Y結線の1相分の値
逆起電圧定数：誘起電圧（相電圧）の0-peak値/角速度
例で示しているパラメータは出力15kWのサーボモータのものです。
以上、今後は速度ループの追加や、d軸電流固定値の追加など行っていきたいと思います。
