数値表現のインピーダンスをα-β座標系からd-q座標系に移す
はじめに
前記事にて三相LCL回路のd-q座標系におけるインピーダンスをHarnefors氏の論文に基づいて(SymPyの力添えを拝借しながら)導出し,LTspiceによる周波数スキャン(frequency response analysis, FRA)の結果と一致することを確認しました。
その際,今後の議論として以下のような宿題が残りました。
もし,(静止座標系でのインピーダンスである)Zs(s)のトポロジーや回路定数が分からず,したがって数式表現が未知であり,例えばインピーダンスアナライザでの測定結果のようにZs(jω)の数値データだけが手元にあった場合,これを使ってd-q座標系でのインピーダンスZd(jω),Zq(jω)を(近似でもよいので)得ることはできるでしょうか?
前記事を書き終えた時点では,「Zs(jω)の数値データからシステム同定の手法にてZs(s)の数式表現を推定するしかないのではないか」と考えておりました。とは言え,筆者自身でもZs(s + jω)をZs(s)の周りでテイラー展開してみたりと試行錯誤を繰り返していました。
そんな中,ノルウェー科学技術大学(Norges teknisk-naturvitenskapelige universitet, NTNU)のAtle Rygg氏の学位論文(下記リンク参照。以下,「Rygg氏の学位論文」と呼びます)において,d-q座標系,修正対称座標系(modified sequence domain),従来の対称座標系で表現したインピーダンスの間の関係が整理されており,それらの間の座標変換について,簡潔にまとまられていることを見つけました。
これを用いればZs(jω)の数値データからd-q座標系でのインピーダンスZd(jω),Zq(jω)を簡便に計算できることが分かりました。前記事にて取り扱った三相LCL回路に適用してみましたので,本記事にて紹介致します。
また,筆者なりにHarnefors氏の論文とRygg氏の学位論文の比較考察をしてみましたので,馬鹿げた記述かもしれませんが一応書いておきます。
結論
まず結論をざっと示します。静止(α-β)座標系でのインピーダンスがZs(s)である場合(三相平衡を仮定します。また上付きsはstationary = 静止座標系を意味します),Harnefors氏の論文に基づいてd-q座標系に移す場合,
と計算することができます。d軸とq軸にインピーダンスを分けるには,(1)式右辺のZs(s + jω1)の実部と虚部を(ラプラス演算子sはひとまず実数と見なして)計算する必要があります。しかし,もともとのZs(s)の数式表現が手元になく,周波数伝達関数としてのZs(jω)の数値データしか入手できない場合,(1)式に相当する計算は難しそうに見えます。
しかし,Rygg氏の学位論文によれば,何とZs(jω)から
のように簡単に計算することができます。そう,Zs(jω)の周波数をω1だけ両方向にシフトしたものを作り,その和を2で除したもの(つまり平均)がd軸成分,その差を2jで除したものがq軸成分となります。以下,実際にやってみましょう。
目次
各インピーダンスの関係
ここではRygg氏の学位論文に記載されているインピーダンスの変換について概観します。
対称座標系でのインピーダンス
三相システムにおける対称座標系については広く知られていることですが,例えば各相の電圧Va(jω),Vb(jω),Vc(jω)から零相電圧V0(jω),正相電圧Vp(jω),逆相電圧Vn(jω)へは次式のように変換します。
同様の変換は電流についても成立します。インピーダンスが三相平衡であり,正相成分と逆相成分の干渉がない場合,正相回路と逆相回路は独立となり,それぞれのインピーダンスをZp(jω),Zn(jω)とすれば,
と書けます。
修正対称座標系でのインピーダンス
Rygg氏の学位論文は修正対称座標系と呼ばれる次式のような表現方法を導入しています。
(5)式は周波数伝達関数として
と書くこともできます。(5)式のZpn(s)の非対角成分Zpn(s)とZnp(s)は異なる周波数成分の間のカップリングを表現しています。Rygg氏の学位論文ではこれをmirror frequency coupling*1と称しています。正相成分と逆相成分のカップリングがない場合,Zpn(s) = Znp(s) = 0となり,Zpn(s)は対角成分のみとなります。この場合,
と書くことができます。Rygg氏の学位論文ではこのような系をmirror frequency decoupled (MFD)と称しています。三相平衡である場合,MFDになります。
修正対称座標系とd-q座標系の関係
d-q座標系での電圧,電流,インピーダンスは次式のように表現できるとします*2。
ここで,MFDな系である場合,Zdq(s)の対角成分は等しく,非対角成分は大きさが等しく逆符号となります。すなわち,
となります*3。
さて,Rygg氏の学位論文には,修正対称座標系とd-q座標系の間で,電圧・電流は次式のように変換できると記載があります。
(10)式を(5)式に代入すれば,
となり,さらに両辺に左から
を掛ければ,
となります。以上より,変換行列AZがユニタリ行列になるようにうまく係数をかけると,Zdq(s)とZpn(s)の変換は次式のように書けます。
ただし,
です。
修正対称座標系とα-β座標系の関係
三相平衡な場合,静止(α-β)座標系では,
です。修正対称座標系では,(5)式の通り,周波数がs + jω1,s − jω1である場合のインピーダンスZpp(s),Znn(s)を考えるため,α-β座標系でのインピーダンスZs(s)の周波数をシフトすればよいと言えます。すなわち,
となり,周波数伝達関数としては,
と書けます。また,三相平衡の場合,Zpn(jω) = Znp(jω) = 0です。(18)式を(14)式に代入して整理すると,
を得ます(これは(2)式の再掲です)。
計算例
(19)式はZs(jω)の数値データがあれば,横軸(周波数)をシフトさせるだけで計算できそうです。そこで,前記事で取り扱った三相LCL回路を例に計算してみました。図1,2に計算例を示します。
前記事にてHarnefors氏の論文に基づいてSymPy-NumPyの合わせ技でanalyticalに計算した結果と,本記事にてRygg氏の学位論文に基づいて(19)式からnumericalに計算した結果が,見事に一致することが確認できます。
2つの方法の比較考察
以下,筆者としていろいろ考えてみた結果を書き散らしています。
Harnefors氏の表現とRygg氏の表現
さて,2つの異なる表現が同じ結果を与えました。Harnefors氏の表現とRygg氏の表現は等しいはずなので,
が成り立っているはずです。(20)式を見ると,どことなくコサインとサインの指数関数表記
にも似ています。むむッ! これって「テイラー展開してみなはれ」という天の思し召しではありませんか?
テイラー展開してみる
という訳で,(20)式右辺に登場するZs(s + jω)とZs(s − jω)をsの周りでテイラー展開してみましょう*4。
(22),(23)式の右辺を見ると,jω1の偶数次の項は実数に,奇数次の項は虚数になることが分かります。 (22)式と(23)式を辺々加えて2で割ると,それぞれの右辺の奇数次の項がキャンセルするため,
となります。また,(22)式から(23)式を辺々引いて2jで割ると,それぞれの右辺の偶数次の項がキャンセルするため,
と書けます。
ところで,(22)式を実部と虚部に分けることができれば,すなわちHarnefors氏の論文におけるZd(s)とZq(s)に分けることができると言えます。(22)式の右辺を見ると分かる通り,偶数次の項が実数に,奇数次の項が虚数となります。実際に分けてみると,
を得ます。項毎に比較してみると,(26),(27)式右辺はそれぞれ(24),(25)式右辺と等しいことが分かります。すなわち,(20)式が成り立っています。 当初筆者は,無限に続いていく(26),(27)式の項を途中で打ち切れば近似できるのではないかとして,Zs(jω)をjωで何度も微分したりと試行錯誤をしていましたが*5,そのような無限の項を扱うことなく,実部と虚部に簡単に分けることができるというのは大変に不思議です。
数学的に厳密かどうかまったく自信がありませんが,筆者としては,Harnefors氏の論文における表現と,Rygg氏の学位論文における表現は等しいと言えるのではないかと考えます。
LTspiceでの波形の観察
さて,(19)式についてはLTspiceでの波形を見ると,自ずと導き出された可能性があります。図3に三相LCL回路のLTspiceスキマティックを示します。
図3は前記事にてFRAを実施した際に使用した回路ですが,今回はd-q座標系において,
フェーザ表現で,
とした場合の挙動の波形を詳しく見てみます。
図4にシミュレーション波形を示します。また,図5に,三相座標系での電圧v(a), v(b), v(c)および電流ia, ib, ic(電圧信号となるので図5ではv(ia), v(ib), v(ic)と表記されています)のスペクトルを示します。また,図6にd-q座標系での時間波形,d-q座標系でのベクトル図,α-β座標系でのベクトル図の関係を筆者なりに考察した図を示します。以下の説明で図6を適宜参照します。
(28)式および図6(a)の通り,d軸電圧を振幅1 V,周波数70 Hzの正弦波*6とし,q軸電圧を零としています。これをd-q座標系のベクトル図としてプロットすると,図6(b)のように正相成分vp,逆相成分vnの重ね合わせとみることができます。vpとvnのq軸成分が打ち消しあうので,電圧ベクトルvはd軸成分のみを持つことになります。
d-q逆変換,α-β逆変換の後に得られる三相座標系での電圧v(a), v(b), v(c)を図4で確認すると,正弦波にはなっていないことが分かります。また,各相で電圧波形が異なることも分かります(ピーク付近を見比べると,山の位置が相毎に異なります)。図5に示すように,LTspiceのFFT機能を使ってスペクトルを分析してみると,v(a), v(b), v(c)には入力した70 Hzとd-q座標系の回転周波数50 Hzの和120 Hzと差20 Hzの成分が含まれていることが判ります。2つの成分の振幅は同じです。これは,図6(c)に示すように,α-β座標系ではvpとvnの周波数がω + ω1 (120 Hz)とω − ω1 (20 Hz)に割れてしまうためです。また,図5を見ると,v(a), v(b), v(c)の120 Hz成分の相順はa, b, cとなっており正相成分,20 Hz成分の相順はa, c, bとなっており逆相であることが判ります。
当然ながら三相LCL回路のアドミタンスはω + ω1 (120 Hz)とω − ω1 (20 Hz)で異なるので,図6(d)に示すように,流れる電流ia, ib, icの120 Hz成分と20 Hz成分の振幅や位相は異なります。なお,三相LCL回路が三相平衡であれば,正相120 Hzの電圧を印加されれば電流も正相120 Hz,逆相20 Hzの電圧を印加されれば電流も逆相20 Hzとなります。本記事では取り扱いませんでしたが,不平衡であれば,正相と逆相のカップリングが起きるはずです。また,線形回路であるため,120 Hzでも20 Hzでもない第3の周波数成分が生まれたりはしません。結果として,今回の回路パラメータや70 Hzという入力周波数の下では,図5の3,4段目のように20 Hz成分が卓越する形となります。
これをα-β変換,d-q変換して再びd-q座標系に戻すと,図6(e)に示すように,正相成分の120 Hz成分は50 Hzが減じられて70 Hzに,逆相成分の20 Hzは50 Hzが加えられて同じく70 Hzに戻ってきます。正相電流ベクトルipと逆相電流ベクトルinを合成した電流ベクトルiの先端は楕円を描きます。これを時間波形として描けば図6(f)となり,シミュレーション波形としては図4の最下段になります。
このように,静止座標系でのインピーダンスを回転(d-q)座標系に移そうとすると,必然的に正相・逆相成分に分けて考える必要があり,しかもそれぞれが静止(α-β)座標系でω + ω1 (120 Hz)とω − ω1 (20 Hz)という異なる周波数に変換されることから,(19)式のようにそれぞれのインピーダンスの和や差として導出されると定性的には理解することができます。
まとめ
本記事では,静止(α-β)座標系でのインピーダンスZs(jω)の数値データから回転(d-q)座標系でのインピーダンスZd(jω),Zq(jω)を簡便に計算する方法を,Rygg氏の学位論文を紹介しながら説明しました。また,Harnefors氏の論文に記載されている数式との関係をテイラー展開を用いて(誤っているかもしれませんが)考察してみました。
さらに,LTspiceを用いたシミュレーション波形から,回転(d-q)座標系から見たインピーダンスは,静止(α-β)座標系における2つの周波数におけるインピーダンスから導かれることを考察しました。
基本波正相成分を直流に見せようと願った分,他の周波数成分を2つに割らずにはいられない…。 d-q座標系ってそういう仕組みだったんだね…。 あたしって,ほんとバカ…。
何と言うか,d-q変換を使った系統連系変換器の制御にこれまで長く親しんできましたが,つい最近までこのような考えに至ることができず,不勉強を恥じてばかりです。
さて,前記事の最後にも書きましたが,今後はパッシブなインピーダンスの方ではなく,三相系統連系変換器の制御系を含めたインピーダンスについて,同じくRygg氏の学位論文やHarnefors氏の論文,また,F. Blaabjerg先生の論文や教科書を参考にしながら勉強していきたいと思っています*7。
付録: Pythonでの計算方法について
本記事のグラフ(図1,2)を計算した際に書いたJupyter NotebookをGitHubに置きました*8。
Rygg氏の学位論文に基づいて(19)式の計算を行うにあたって,ひとつ工夫したところがあります。往々にして,広い周波数範囲のインピーダンスをデータ化しようとした場合,周波数軸は対数にする(片対数グラフとする)のが一般的かと思われます。その状態で,±ω1 = ±2 π × 50 [rad/s]だけ周波数をシフトしようとした場合,つまりZs(jω)のデータからZs(j(ω ± ω1))を得ようとした場合,ちょうどその周波数でのデータが存在しないことが予期されます(特に50 Hzや60 Hzが相対的に小さく見える高周波領域)。そこで,
from scipy import interpolate ## ±jω1を自由に足し引きできるように線形補完を定義する Zs_interp = interpolate.interp1d(_omega, Zs_bode, fill_value = 'extrapolate') ## 周波数をシフトしてインピーダンスの対角・非対角成分を計算する Zd_rygg = (Zs_interp(_omega + _omega1) + Zs_interp(_omega - _omega1)) / 2 Zq_rygg = (Zs_interp(_omega + _omega1) - Zs_interp(_omega - _omega1)) / 2j
のように,scipy.intepolate.interp1d
関数を用いて,元のZs(jω)すなわちZs_bodeを線形補間するようにしました。図1,2は上のZd_ryggとZq_ryggを_omegaに対してプロットしたものです。
*1:Mirror frequency coupling = 鏡像周波数カップリングとでも訳しましょうか。
*2:Harnefors氏の論文とRygg氏の学位論文では非対角成分の添え字の付け方が逆になっています。(8)式ではRygg氏の書き方に従います。
*3:Harnefors氏の論文では対角成分をZd(s),非対角成分を±Zq(s)と表記していますが,Rygg氏の学位論文ではそれぞれZdiag(s),±Zoffdiag(s)と表記しています。(9)式では前記事に合わせてHarnefors氏の表記方法に従います。
*4:sはラプラス演算子ですが,そのsの周りでのテイラー展開などと言うものができるのかどうか,筆者は残念ながら知りません…💧
*5: 三相LCL回路の静止座標系におけるインピーダンスZs(jω)の数値データから回転座標系におけるZd(jω)とZq(jω)を得る試み…。
どうもダメっぽい💧
いったん頭を冷やそう😅 pic.twitter.com/WcOOMoqyvb
*6:コサインなので厳密には「余弦波」と呼ぶべき? 広義にはコサインも「正弦波」に入ると思っていますが…。
*7:このところNucleoやRaspberry Piに触っていないので,再びプログラミングや電子工作でも遊んでみたいという気もしています…💧