モータ制御の引き出し

モータ制御の仕事をしています。

PI制御電流帰還における伝達関数のエクセルマクロシミュレーション

今回はベクトルモータ制御の電流帰還制御に使われるPI制御の伝達関数周波数応答シミュレーションをエクセルマクロを用いて行う方法を説明しようと思います。

 

①PI制御電流帰還オープンループ伝達関数の周波数応答エクセルマクロシミュレーション

モータのベクトル制御でのPI制御を用いた電流帰還は非干渉制御を行うと、下図のようになり、オープンループゲインは以下の式になる。                           f:id:runge0704:20210816172829p:plain

f:id:runge0704:20210816173205p:plain

X(s):電流指令値(入力)

Y(s):モータ電流値(出力)

Ca(s):PI制御部関数

Cm(s):モータ負荷関数

R:モータの抵抗値

L:モータのインダクタンス

K:比例制御の比例定数

Ta:積分制御の時定数

Tm:モータの時定数(L/R)

sはラプラス演算子 

 

オープンループなので、帰還は考えないようにすると、伝達関数(出力/入力)は

Ca(s)・Cm(s)となる。

本来はd軸、q軸と別々に制御するが、今回はそれを省略する。

モータ負荷部Cm(s)ブロックの入力はモータ印加電圧であり、そこからモータ電流を導くCm(s)はすなわちアドミタンスであり、それは1/(R+sL)の一次遅れになり、モータの時定数Tm(L/R)を用いて表すと、1/(R(1+Tms))となる。

では、ボード線図を作成するために伝達関数の絶対値を求める。

絶対値計算は分母と分子をそれぞれ計算する。

伝達関数を変形して以下の通りにする。

f:id:runge0704:20210817232200p:plain

分子、分母それぞれの絶対値を算出する。

〇分子 s=jωを分子の式に代入すると

 K*(Tas+1)=K*(1+jTa*ω)

 絶対値を計算する。

 |K*(1+jTa*ω)|=K√(1+(Ta*ω)^2)

 

〇分母 s=jωを分子の式に代入すると

 Tas*R*(1+Tms)=Ta*R*(jω-Tm*ω^2)

 絶対値を計算する。

  |Ta*R*(jω-Tm*ω^2)|=Ta*R*ω√((Tm*ω)^2+1)

 

オープンループゲインの絶対値は

Ca(s)・Cm(s)=|分子|/|分母|にて計算することができる。

 

次にモータ電気定数R,Lと制御器の設定応答角周波数ω0から

最適比例定数Kと制御器時定数Taを算出する。

①モータ定数RとLを入力する。

   例えばR=0.05Ω、L=0.002Hとする。

②比例定数は制御器の応答角周波数(ω0)から決まる。

   応答角周波数(ω0)と設定すると、比例定数Kは

   K=ω0Lと導くことができる。

  これは積分制御なしの比例制御のみを考えたときに、

  応答角振動数ω0と時定数Tは次の関係になる。

   T=1/ω0 

<K=ω0Lの証明>

        例えばω0 を500rad/s、L=0.002Hとして、K=ω0Lより

   K=1としてステップ応答のシミュレーションをすると

   下図のような立ち上がりとなり、時間0.02sec(1/ω0)のとき、

   値が6.59 (指令値が10)となり、ほぼ時定数(1-1/e   0.632ですが)

           になる。

f:id:runge0704:20210912125614p:plain

 ③積分制御の時定数Taの求め方

 Ta=Tmとすると、オープンループゲインの式

 f:id:runge0704:20210817232200p:plain

  においてTaS+1 と1+TmS が同じになるので、

 オープンループゲインはK/TaSRとsに対して1次の簡単な式になる。

 そのためオープンループゲイン周波数応答シミュレーションは下図のように

 極、零の無い直線になる。対数軸基準なので、直線=反比例となる。

 このようなTaの求め方を極零相殺法と言う。

 Tm=L/Rなので数値を代入して、Tm=0.04となるので、

 Taも0.04とした時の周波数応答シミュレーション結果を下図に示す。

f:id:runge0704:20210830002326p:plain

 さらにTaを10倍、0.1倍と変化させた時の周波数応答シミュレーション結果を下図に示す。

また、このゲイン曲線の違いがクローズループになるとどうなるかを後程シミュレーションする。

f:id:runge0704:20210901010710p:plain

Ta=0.4の時はTmによる極点が25rad/secに見られ、Taによる零点2.5rad/secと10倍離れていることが分かる。(極・零はω=1/Taから算出できる。)

オープンループゲインとクローズループゲインの周波数応答エクセルマクロシミュレーションはここからダウンロードしてください。

(右クリックで「新しいウインドウで開く」を行い、ダウンロードしてください。) 

エクセルマクロ

 

PI制御電流帰還クローズループ伝達関数周波数応答エクセルマクロシミュレーション 

 

非干渉制御を行うと、PI制御を用いた電流帰還は下図のようになり、

オープンループゲインは以下になる。 

f:id:runge0704:20210829231057p:plain

 

f:id:runge0704:20210829231243p:plain


K:比例制御の比例定数

Ta:積分制御の時定数
R:モータの抵抗値               

Tm:モータの時定数(L/R)

そのクローズループゲインは上のブロック図より

Y(s)/(X(s)-Y(s))=Ca(s)*Cm(s)

Y(s)=Ca(s)*Cm*1

Y(s)(1+Ca(s)*Cm(s))=Ca(s)*Cm(s)*X(s)

Y(s)/X(s)=Ca(s)*Cm(s)/(1+Ca(s)*Cm(s))

となるので、クローズループゲインは

f:id:runge0704:20210829231851p:plain

f:id:runge0704:20210829231914p:plainTas(R(1+Tms)を分母分子に掛ける。

 

 

絶対値計算は分母と分子をそれぞれ計算する。

|分子| K√((1)^2+(Taω)^2) S=jωを代入する。

|分母| √(K-TaTmRω2)2+(KTaω+RTaω)2)

クローズループのゲイン絶対値は|分子|/|分母|となる。

 

オープンループゲインの時と同じ定数にて周波数応答のシミュレーションを行った結果を下図に示す。Taは標準の0.04と0.4,0.004とした。

f:id:runge0704:20210901005914p:plain

Ta=0.004の時、ω=500rad/sec付近でピークを持ってしまい、

発振しやすくなってしまう。

また、初めにω0を500rad/secと設定したが、シミュレーション結果より

カットオフ角振動数は500rad/secと読み取れる。

それに付加してTa=0.04、ω0を250、500、1000つまりKを0.5、1、2とした時のシミュレーションを下図に示す。ω0を倍、半分にすることでカットオフ角周波数も変化する。

 

f:id:runge0704:20210904084000p:plain

但し、ここで注意することはモータ制御のキャリア周波数とカットオフ角周波数
の関係について、デジタルサンプリング誤差を考慮すると、カットオフ周波数(f0=ω0/2π)*10<キャリア周波数 以上の差が必要になる。

 

クローズループの時間応答シミュレーションはここからダウンロード出来ます。

時間応答のエクセルマクロ

(右クリックで「新しいウインドウで開く」を行い、ダウンロードしてください。) 

マクロはTaを3パターンとしたシミュレーションです。

 

下図にクローズループの時間応答シミュレーション結果を示す。

f:id:runge0704:20210912142727p:plain

Ta(積分時間)を標準から1/100(ゲインは100倍)をした状態のTa3では

起動波形が不安定になっており、発振の可能性がある。

 

以上

*1:s)X(s)-Y(s

VBAによるモータベクトル制御シミュレーション

モータベクトル制御のシミュレーションとしてはPSIMやLTspice,simulinkなどGUIを使ったものがありますが、これらは結局C言語ソフトウェアでも実施でき、それならばエクセルに付属のVBA(エクセルマクロ)でもできるのではないかと思い、作成してみました。

VBAにてシミュレーションする利点は、ソフトのインストールなど必要なく、エクセルさえあれば簡単に実行できることです。

例えばソフト作成者はベクトル制御の各ステップの計算値を確認したり、グラフ化することでベクトル制御への理解を深めることができるでしょう。

今回は条件として、

Id=0 Iq:電流指令での電流ベクトル制御としました。

マクロ込みのエクセルファイルはここからダウンロードできます。

エクセルマクロ

(右クリックし、新しいタブで開くで開き、ダウンロードする。)

 

0.使い方

事前準備としてエクセルソフトを開き、ファイル⏎その他/オプション⏎リボンのユーザ設定ウインドウの「開発」にチェックを入れる。

f:id:runge0704:20210307230040p:plain

ダウンロードしたマクロ付きエクセルファイルを開き、「開発」タブを開き、

「マクロ」ボタンでマクロウインドウを開き、ID_0_ベクトル制御を実行します。

f:id:runge0704:20210307230358p:plain

結果がデータシートに表示されますので、範囲を選択する形でグラフを表示します。




1.ダイアグラム

f:id:runge0704:20210307192925p:plain



電流指令をId=0、Iq=iref として入力し、UVW相のモータ電流を検出したiu1,iv1,iw1をdq変換し、フィードバックしています。d軸サーボコントローラ、q軸サーボコントローラではPI制御を行い、その出力を逆dq変換し3相モータ出力電圧u1,v1,w1として出力します。

インバータ実物は効率の面でPWM出力ですが、シミュレーションでは制御を簡単にするため省略しています。モータ電流をシミュレートするに必要なモータ逆起電力とdq変換、dq逆変換に必要な位相(theta)はモータモデルから速度の積分として得られるようにしています。

モータモデルはモータに発生するトルクはId=0なので磁石トルクのみと考えられるので、T=KT*Iとしてシミュレートしています。(KTはトルク定数)

 

 2.コード

〇変数、定数の宣言

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の三相電流からαβ軸変換波形は下図のようになります。

f:id:runge0704:20210306162741p:plain

iu1とa1は重なります。b1はa1から90度遅れになります。

 

αβ軸からdq軸への変換は下図の波形になります。

f:id:runge0704:20210306163602p:plain

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を算出する。コードを下記に示す。

 

dq→αβ 逆dq変換


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は重なっています。

f:id:runge0704:20210307104134p:plain

本来はこの後で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=0制御では磁気トルクのみとなるので、単純に

コイル電流m3とトルク定数KTから求められます。

m4=KT×m3

 

モータモデルからモータ速度m7、モータ位相m8、モータ逆起電圧thetaを算出するのは下記のコードとなります。

m1 = Vq - m8 '出力電圧ー誘起電圧
m2 = m1 / Ra
m3 = rungekut(m2, m3, taue, dt) '4次元ルンゲクッタ法 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) '4次元ルンゲクッタ法 La,Raによるフィルタ

iu1=iu3’結果をiu1に代入する。

〇電源電圧について

実際のハードではインバータの電源は片電圧ですが、シミュレーション上は

両電源としています。

両電源とすることで、モータ中点電圧=0Vとしてシミュレーションすることができ、オフセット電圧を考慮しないで済みます。

〇シミュレーションステップについて

本シミュレーションは時間軸を基にしたシミュレーションですが、ステップ時間dt

を一定にすると位相急変時に位相変化量が1ステップで0.1°以上になることがあり、そのことでシミュレーションが発振する可能性があるため下記コードにより制限しています。

If m7 < 10 And m7 > -10 Then
carF = 600
Else: carF = m7 * 60

End If
dt = 1 / carF

〇モータパラメータ

 本シミュレーションに用いるモータパラメータはY結線の1相分を基準とします。

トルク定数:発生トルク/コイル電流(相電流)の0-peak値

コイルインダクタンス:Y結線の1相分の値

コイル抵抗:Y結線の1相分の値

逆起電圧定数:誘起電圧(相電圧)の0-peak値/角速度

例で示しているパラメータは出力15kWのサーボモータのものです。 

 

以上、今後は速度ループの追加や、d軸電流固定値の追加など行っていきたいと思います。