自作のVVVFインバーターの惰行時の制御

フィードバックなしに空回しモーターの惰行時の回転速度を推定する方法を紹介します。

前置き

こんにちは。自作のVVVFインバーターでモーターを回すのは楽しいですよね。
f:id:yuppi5:20210220165241j:plain

ところで、電車用のインバーターでは回転中に電圧を0にして惰行（だこう）する制御があります。普通のモーター用のインバーターでは回転中に電圧が0になることはないので、惰行は電車用のインバーター特有の制御です。

私が作っているインバーターでは電車のインバーターの再現をするのが目的なので、当然惰行もします。加減速はフィードバック無しでもある程度平気(V/f一定制御)ですが、惰行はフィードバックが無いと厳しいです。インバーター周波数と実際のモーターの回転数がかけ離れた状態で電圧を上げると、過大な電流が流れます。スイッチング素子が壊れるかもしれません。惰行をする場合、フィードバックは必須です。

しかし、私のインバーターでは負荷をかけてモーターを回すことはなく、空回しをさせています。よって計算することでフィードバックなしにモーターの回転数を推定できると考えました。

モーターを惰性で回して回転数を測定する

f:id:yuppi5:20170901100753j:plain
フォトインタラプターとマイコンで測定しました。マイコンはLPC1114を使用しました。またこのマイコンのGPIOにはヒステリシス機能がついていました。この機能を有効にして外部割り込みでうまく回転数を測定しました。
回転数のグラフは以下のようになりました。横軸が時間、縦軸が回転数です。
f:id:yuppi5:20170901213728p:plain
インバーター周波数を180Hzにして落ち着いたあとに惰行を開始しました。この瞬間のモーターの回転数は約90rpsですが、都合上2倍にして表示しています。
180Hz(90rps)から、指数関数のような形で回転数が下がっていき、28秒あたりで回転が止まりました。
この線は4回測定した中央値ですが、線がガタガタですね。200回くらい試行しないといけなかったでしょうか。

曲線を近似する

上のようなグラフを曲線で近似して、回転数推定に必要なパラメーターを算出します。空回しなので、急に止まったりすることはなく、自然な減速をしていくはずです。結果だけ見たい方は近似した曲線と必要なパラメーターを算出するスクリプトに飛んでください。値を突っ込むと必要なパラメーターが出てくるスクリプトがあります。

微分方程式を立てる

モーターの構造を単純化します。
f:id:yuppi5:20170903003005p:plain
回転数を {\displaystyle f} と置き、

  • 摩擦抵抗 { \displaystyle = -\mu N = -T } (定数なのでまとめて { \displaystyle T} と置いた)
  • 回転周波数に比例する空気抵抗 { \displaystyle = -kf } ( { \displaystyle k} は定数)

があると仮定します。
回転周波数の2乗に比例する抵抗は面倒なので計算していません。
空気抵抗は、モーターの外扇の羽や回転子から生まれると考えます。

角加速度 { \displaystyle \frac{df}{dt} } は上記2つの抵抗の合計になります：
\[
\frac{df}{dt} = - kf - T
\]
特定したいのは、 { \displaystyle T} と { \displaystyle k} の値です。
私のインバーターは角加速度でモーターの回転数を増減させているので、この2つの値がわかれば十分です。
[追記]本来の運動方程式は { \displaystyle m \frac{df}{dt} = - kf - T } ( { \displaystyle m} は質量)と質量が必要ですが、 { \displaystyle k} と { \displaystyle T} は { \displaystyle m} ですでに除算されたことにすれば問題ないです。完全に質量のことを忘れていました。

微分方程式を解く・数値解析

※対数の底は { \displaystyle e}とします。
{ \displaystyle \frac{df}{dt} = - kf - T } を変形します：
\[
\frac{df}{dt} = - kf - T \\
\frac{ 1 }{ f + \frac{T}{k} } \frac{df}{dt} = -k
\]
両辺を { \displaystyle t } で積分します：
\[
\int \frac{ 1 }{ f + \frac{T}{k} } \frac{df}{dt} dt = - \int k dt \\
\int \frac{ 1 }{ f + \frac{T}{k} } df = - \int k dt \\
\log \left| f + \frac{T}{k} \right| = - k t + C'
\]

C'は積分定数

\[
f + \frac{T}{k} = \pm e^{-kt + C'} \\
f + \frac{T}{k} = C e^{-kt} \\
(C = \pm e^{C'}) \\
f = C e^{-kt} - \frac{T}{k} \ ...\ (1)
\]
ここで初期回転数 { \displaystyle f_0 } を代入して { \displaystyle C } を消します。
回転数 { \displaystyle f } は { \displaystyle t=0 } のとき { \displaystyle f_0 } だから：
\[
f_0 = - \frac{T}{k} + C \\
C = f_0 + \frac{T}{k} \ ...\ (2)
\]
式(1)に式(2)を代入して：
\[
f(t) = \left( f_0 + \frac{T}{k} \right) e^{-kt} - \frac{T}{k}
\]
とりあえず式は解けました。
ここから変形をして { \displaystyle T } を消します。
モーターの回転が止まった { \displaystyle (f(t)=0)} ときの時間を { \displaystyle t_{end} } と置くと
\[
0 = \left( f_0 + \frac{T}{k} \right) e^{-k t_{end} } - \frac{T}{k} \\
\frac{T}{k} = \frac{ f_0 e^{-k t_{end}} }{ 1 - e^{-k t_{end}} } \ ...\ (3)
\]
・{ \displaystyle T } は式(3)に { \displaystyle k } を掛けると算出できます：
\[
T = \frac{ f_0 e^{-k t_{end}} }{ 1 - e^{-k t_{end}} } k \ ...\ (4)
\]

式(1)に式(3)を代入：
\[
f(t) = \frac{ f_0 }{ 1 - e^{-k t_{end}} } ( e^{-k t} - e^{-k t_{end}} )
\]
ここから、適当な実測点 { \displaystyle t_1, f_1(=f(t_1))} を取って { \displaystyle k } に関する公式を出しますが…
\[
f_1 = \frac{ f_0 }{ 1 - e^{-k t_{end}} } ( e^{-k t_1} - e^{-k t_{end}} ) \\
0 = \frac{f_0}{f_1} \frac{ e^{-k t_1} - e^{-k t_{end}} }{ 1 - e^{-k t_{end}} } - 1 \ ...\ (5)
\]
式(5)が { \displaystyle 0 } になるような { \displaystyle k } を出せばいいんですが、これを解析的に解く({ \displaystyle k=□} の形にする)ことは不可能です。
{ \displaystyle k } に適当な値を当てはめて、数値的に解きます。

式(5)を頑張って変形すると：
\[
k_{n+1} = - \frac{1}{t_1} \log \left( \frac{ f_1 }{ f_0 } (1 - e^{-k_n t_{end}}) + e^{-k_n t_{end}} \right)
\]
初期値 { \displaystyle k_0 } (適当に 0.1 を代入しています)を入れて、出てきた { \displaystyle k_1 } を右辺に代入し計算して { \displaystyle k_2 } を出します。これを右辺に代入し計算して { \displaystyle k_3 } を出します。
という風に、何回か繰り返して { \displaystyle k } の値を数値的に算出します。15回ほど繰り返すといい具合に収束します。
この数値計算プログラムは次の近似した曲線と必要なパラメーターを算出するスクリプトに置いておきます。

近似した曲線と必要なパラメーターを算出するスクリプト

必要な値：
f:id:yuppi5:20170902215840p:plain

  • { \displaystyle f_0 } : 初期回転数
  • { \displaystyle t_1} : ある点(最初と最後以外の点・真ん中の少し左側あたりが良い)
  • { \displaystyle f_1} : { \displaystyle t_1}のときの回転数
  • { \displaystyle t_{end} } : 回転が止まったときの時間

{ \displaystyle k,\ T } は：

  • { \displaystyle k_{n+1} = - \frac{1}{t_1} \log \left( \frac{ f_1 }{ f_0 } (1 - e^{-k_n t_{end}}) + e^{-k_n t_{end}} \right) }
  • { \displaystyle T = \frac{ f_0 e^{-k t_{end}} }{ 1 - e^{-k t_{end}} } k }

近似する曲線の式は：

  • { \displaystyle f(t) = \frac{ f_0 }{ 1 - e^{-k t_{end}} } ( e^{-k t} - e^{-k t_{end}} ) }

「微分方程式を解く・数値解析」より、{ \displaystyle k } は数値計算で求めます：
・ { \displaystyle f_0 } :
・ { \displaystyle t_1} :
・ { \displaystyle f_1} :
・ { \displaystyle t_{end} }

・ { \displaystyle k} :
・ { \displaystyle T } :
入っている値は、今回測定したものです。

測定した回転数から必要なパラメーターを求める

近似した曲線と必要なパラメーターを算出するスクリプトのスクリプトを使います。
もう値がスクリプトに入力されていますが、改めて示します。これが測定した値です：

  • { \displaystyle f_0 = 180 }
  • { \displaystyle t_1 = 12.4 }
  • { \displaystyle f_1 = 41.25 }
  • { \displaystyle t_{end} = 28 }

計算すると、パラメーターは次のようになりました：

  • { \displaystyle k = 0.105929 }
  • { \displaystyle T = 1.03544 }

グラフを書きます：
f:id:yuppi5:20170902221026p:plain
うまくフィットしました。高回転時ほど回転数の増減が激しいので、 { \displaystyle f_1 } は { \displaystyle f_0 } に近いところ、半分くらいのところを選ぶといいと思います。{ \displaystyle f_1 = 41.25 } は低すぎたかもしれません。