The Negligible Lab

Astray in the forest of electrical and electronic circuits. Adrift in the gap between time and frequency domains. 独立独歩を是とする。

FFTを訪ねて［前編］電気屋としての理解を探る

ディジタル信号処理 Python Matplotlib JupyterLab

はじめに

FFTと私

高速フーリエ変換(Fast Fourier Transform, FFT)は，通常の離散フーリエ変換(Discrete Fourier Transform, DFT)に比べて非常に高速なアルゴリズムとして知られており，ランダウの記号を用いると，DFTの計算量がO(N2)であるところ，FFTではO(N log N)になるようです。このため，理学・工学の幅広い分野で応用されております。言わば，ディジタル信号処理の中の1つの大きな“柱”となっています。ただし，入力される時間領域のディジタル信号のデータ点数が2のべき乗（例えば64とか1,024とか4,096など）となっている必要があります*1。

図1に，FFTのうち，本記事で取り扱う周波数間引き型Cooley-Tukeyアルゴリズムの計算途中経過を可視化したGIFアニメーションを示します。Python (JupyterLab)とMatplotlibを用いて作成しました*2。青が実部，橙が虚部を表しています*3。

図1: FFT（周波数間引き型Cooley-Tukeyアルゴリズム）の計算途中経過を可視化

ここでは，時間領域波形f(θ)として，

f(\theta) = 0.8 e^{j \theta} + 0.6 e^{j(7 \theta - \pi / 2)} + 0.4 e^{j(12 \theta - \pi / 2)} + 0.3 e^{j(17 \theta + \pi / 6)} \tag{1}

を与えています*4。jは虚数単位です(j2 = −1)。また，θは図1の横軸全体で0から2πまで変化すると考えて下さい*5。FFTの結果，(1)式の通り，1次（基本波），7次，12次，17次にスペクトルが立っていることが分かります。データ点数はN = 64 (= 26)です。

ところで，筆者はかなり以前から電力・電気系の（自称）技術者（俗にいう「電気屋」）として電力系統の高調波(harmonics)と向き合って参りました。当然ながらフーリエ解析*6のお世話になるのですが，MATLABにしてもNumPyにしても，fftやnumpy.fft.fftとコマンドを打つだけでフーリエ解析ができてしまうので*7，その中身については理解が及んでおりませんでした。しかし，どうしてデータ点数が2のべき乗でなければならないのか，どうしてDFTよりも高速なのか，どのように時間波形からスペクトルを得ているか，本来は理解しないといけないのだろうけれども，よく分からないまま使っているという後ろめたさを抱え続けておりました*8。

また，以前に加速度センサを用いて地震計を作ろう思い立った際，加速度から震度を計算するアルゴリズムにFFTが必須であり*9，挫折した記憶があります。FFTをマイコン（例えばM5Stack Grayなど加速度センサのある系）に実装できれば，地震計の実験もできます。しかし，これも何年間もの休眠プロジェクトとなっています。

インターフェース誌の記事をきっかけとして

そんな中，CQ出版社の「インターフェース」誌の2021年9月号にて「数学とプログラミング」特集が組まれ，マイコン（記事ではソニーのSpresense）にFFTを実装する際の考え方とC言語によるプログラミングの例が記載されておりました*10。これはFFTへの理解を深めるための機運かもしれない。そこで，これまでに遊んできたPython (JupyterLab & Matplotlib)とM5Stack (Arduino IDE)にFFTを実装しながら，FFTの本質を探訪してみようと思い立ちました。

*1:後述しますが，実は3のべき乗でもそれ以上の自然数のべき乗でも，アルゴリズムは存在します。ただし若干効率は下がると思います。

*2:図1では，FFTの結果を見やすくするため，それぞれのバタフライ演算の後でデータを2で除算しています。

*3:最近，1つ目の要素を青，2つ目の要素を橙にするのが流行っていませんか？ MatplotlibもMicrosoft Excelもそうですよね。

*4:このデータは正相成分のみをもっています。

*5:ちょっと乱暴な書き方ですが…ご容赦下さい…。

*6:フーリエ級数展開なのかフーリエ変換なのか，はたまた連続なのか離散なのか，ぼやかして言うときに「フーリエ解析」という言い方は便利だと思っています。

*7:Microsoft Excelでも「分析ツール」でFFTが可能です。

*8:本ブログの筆者の数学力は驚くほど低いことが知られています。

*9:気象庁｜強震観測について｜計測震度の算出方法

*10:2021年9月号 ｜ Interface – CQ出版

続きを読む

M5Stackで遊ぶ ～浮遊静電容量で系統周波数を測る～

電子工作 M5Stack Googleスプレッドシート

はじめに

概論

電力系統の周波数は，言わずもがな，東日本では50 Hz，西日本では60 Hzです。しかし，いつもピッタリ50 Hzや60 Hzとなっているわけではなく，それぞれ50 Hz，60 Hzを中心としながらも，常に揺らいでいます。発電が消費よりも大きければ周波数は上昇し，逆に消費が発電よりも大きければ周波数は低下します。現状，水車や蒸気タービンで回転している同期発電機が大勢を占めていることから，負荷が重ければ減速し，軽ければ加速する──とも言い換えられますね*1。

もちろん，50 Hz，60 Hzから離れ過ぎると問題が生じるので，地域によって異なりますが，±0.2 ~ ±0.3 Hzに収まるよう（±0.1 Hzの滞在率が概ね99%以上*2），水車に流す水量や，ボイラに投入する燃料を調整したり，周波数を観測して特定の発電所を制御したり，時刻，曜日，季節に応じて運転・停止する発電所を選択したりといった措置が取られています。

実際どうなの？

諸外国には，電力系統の周波数をリアルタイムにインターネット上に公開しているTSO (transmission system operator)もあります。例えば，スウェーデンのSvenska Kraftnätなどがそうです*3。

しかし，日本ではこれは実施されていません。電力広域的運営推進機関(OCCTO)ホームページでも，地域間連系線（周波数変換所や直流送電を含む）の融通量を見ることはできますが，周波数は見ることができません。実際，どういう風に周波数が変動しているのか？

*1:太陽光や風力など，系統連系インバータを介した発電設備が増えてくるとこの挙動が変わってくるとして，その対策に関する研究が最近では盛んです。

*2:https://www.occto.or.jp/iinkai/chouseiryoku/jukyuchousei/2020/files/jukyu_shijyo_19_02_02.pdf

*3:The control room | Svenska kraftnät

続きを読む

LTspiceによる三相PWM整流器の上位制御とノーズカーブ

はじめに

サイリスタの転流を交流側の電圧に依存している他励インバータでは，電力系統（以下，単に「系統」とも呼びます）の短絡容量が小さくなると（＝系統が弱いと）安定運転が難しくなると言われています。これに対して，IGBT (insulated-gate bipolar transistor)などのオンオフ制御バルブデバイス（自己消弧素子）を用いる自励インバータでは，系統の短絡容量が小さくても運転できるとされています。とは言うものの，系統の短絡容量と自励インバータの定格電力*1の比であるSCR (short-circuit ratio)が1.5 ~ 2を下回ると，やはり安定性に問題が生じる恐れがあることが指摘されています。

本記事では，本ブログでこれまで取り扱ってきたLTspiceによる三相PWM整流器に，有効電力制御，無効電力制御，交流電圧制御の制御ブロックを追加します。これら3つの制御と直流電圧制御とを合わせて，上位制御系と呼ぶことにします。

また，有効・無効電力制御とSCRの関係について，ノーズカーブの観点から若干の考察を加えます。ノーズカーブとは，系統の電圧安定性を検討する際に描かれる図です（P-V曲線，P-Vカーブとも呼ばれます）。SCRが小さい（＝系統が弱い）場合，三相PWM整流器（自励インバータ）の制御がいかに理想的であったとしても，ある有効電力を超えて受電することはできません。大きな有効電力を得ようと電流を増加させていくと，ある点で連系点の電圧が急激に低下してしまいます。これを電圧崩壊(voltage collapse)と呼びます。これは私の単なる憶測に過ぎませんが，上位制御系の不安定と，単なる電圧崩壊が，これまであまり区別されてこなかったのではないかという思いがあります*2。

本記事はそこで終わりですが，今後，上位制御を含めた三相PWM整流器の小信号モデルを立て，安定性の近似理論解析とLTspiceでの挙動の比較考察に移っていきたいと考えています。

*1:本記事では「定格有効電力」としました。定格皮相電力とするなど，考え方は他にもあるのではと考えます。

*2:いや，私自身の中でもこれをきちんと区別できているのか判然としません?? 俺がそう思うんならそうなんだろ。俺ん中ではな…。

続きを読む

数値表現のインピーダンスをα-β座標系からd-q座標系に移す

回路理論 パワエレ LTspice Python

はじめに

前記事にて三相LCL回路のd-q座標系におけるインピーダンスをHarnefors氏の論文に基づいて（SymPyの力添えを拝借しながら）導出し，LTspiceによる周波数スキャン(frequency response analysis, FRA)の結果と一致することを確認しました。

negligible.hatenablog.com

その際，今後の議論として以下のような宿題が残りました。

もし，（静止座標系でのインピーダンスである）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)，従来の対称座標系で表現したインピーダンスの間の関係が整理されており，それらの間の座標変換について，簡潔にまとまられていることを見つけました。

ntnuopen.ntnu.no

https://www.researchgate.net/publication/328231490_Impedance-based_methods_for_small-signal_analysis_of_systems_dominated_by_power_electronics

これを用いればZs(jω)の数値データからd-q座標系でのインピーダンスZd(jω)，Zq(jω)を簡便に計算できることが分かりました。前記事にて取り扱った三相LCL回路に適用してみましたので，本記事にて紹介致します。

また，筆者なりにHarnefors氏の論文とRygg氏の学位論文の比較考察をしてみましたので，馬鹿げた記述かもしれませんが一応書いておきます。

結論

まず結論をざっと示します。静止(α-β)座標系でのインピーダンスがZs(s)である場合（三相平衡を仮定します。また上付きsはstationary = 静止座標系を意味します），Harnefors氏の論文に基づいてd-q座標系に移す場合，

\begin{align} Z_{d}(s) &= \mathrm{Re} \left [ Z^{s}(s + j \omega_{1}) \right ] \\[10pt] Z_{q}(s) &= \mathrm{Im} \left [ Z^{s}(s + j \omega_{1}) \right ] \end{align} \tag{1}

と計算することができます。d軸とq軸にインピーダンスを分けるには，(1)式右辺のZs(s + jω1)の実部と虚部を（ラプラス演算子sはひとまず実数と見なして）計算する必要があります。しかし，もともとのZs(s)の数式表現が手元になく，周波数伝達関数としてのZs(jω)の数値データしか入手できない場合，(1)式に相当する計算は難しそうに見えます。

しかし，Rygg氏の学位論文によれば，何とZs(jω)から

\begin{align} Z_d(j \omega) &= \frac{Z^{s}(j (\omega + \omega_{1})) + Z^{s}(j (\omega - \omega_{1}))}{2} \\[10pt] Z_q(j \omega) &= \frac{Z^{s}(j (\omega + \omega_{1})) - Z^{s}(j (\omega - \omega_{1}))}{2j} \end{align} \tag{2}

のように簡単に計算することができます。そう，Zs(jω)の周波数をω1だけ両方向にシフトしたものを作り，その和を2で除したもの（つまり平均）がd軸成分，その差を2jで除したものがq軸成分となります。以下，実際にやってみましょう。

続きを読む

三相LCL回路のd-q座標系におけるインピーダンスを求める

回路理論 SymPy LTspice Matplotlib パワエレ

はじめに

系統連系インバータと電力系統の相互作用によって，システムが不安定になることがあります。本ブログでは，過去の記事にて電流制御系と系統インピーダンスの相互作用について検討したことがあります。系統連系インバータの外乱応答（インピーダンス）と系統インピーダンスの比のナイキスト線図を描くことによって，電流制御系の安定性を論じることができました。しかし，この記事では単相の場合しか取り扱っておりませんでした。三相の場合はどのように考えるべきでしょうか。

negligible.hatenablog.com

三相系統連系インバータは大体においてd-q（回転）座標系で制御することが一般的であり，かつ，d軸とq軸で別々の制御を行っています。したがって，三相系統連系インバータはd-q座標系で不平衡なインピーダンスとなることでしょう。そのため，一般的に三相平衡に近いと考えられる系統インピーダンスの方を三相（あるいはα-β（静止）座標系）からd-q座標系に換算して表現するのが好適と考えられます。

しかし，不幸にして筆者は，RL直列回路以外の一般のインピーダンスをd-q座標系で表現する方法を，よくよく考えてみると知りませんでした。その必要性があること自体を想像できておらず，自らの不明を恥じ入っております。

そこで，本記事では，任意のインピーダンスをd-q座標系で表現する手法について，文献を紹介しながら説明します。また，一例として三相LCL回路のインピーダンスのd-q座標系における表現を導出し，これをLTspiceでの周波数スキャン(frequency response analysis, FRA)と比較検証することとしました。

任意のインピーダンスの場合と文献紹介

任意のインピーダンスのd-q座標系での表現

まず，ずばり結論を述べます。あくまでも三相平衡の場合を前提としますが，s領域において，三相座標系での自己インピーダンスがZs(s)，相互インピーダンスがZm(s)である場合，次式が成り立ちます。

\begin{bmatrix} V_{a}(s) \\ V_{b}(s) \\ V_{c}(s) \end{bmatrix} = \begin{bmatrix} Z_{s}(s) & Z_{m}(s) & Z_{m}(s) \\ Z_{m}(s) & Z_{s}(s) & Z_{m}(s) \\ Z_{m}(s) & Z_{m}(s) & Z_{s}(s) \end{bmatrix} \begin{bmatrix} I_{a}(s) \\ I_{b}(s) \\ I_{c}(s) \end{bmatrix} \tag{1}

上式をα-β座標系で表現すると，

\begin{bmatrix} V_{\alpha}(s) \\ V_{\beta}(s) \end{bmatrix} = \begin{bmatrix} Z^{s}(s) & 0 \\ 0 & Z^{s}(s) \end{bmatrix} \begin{bmatrix} I_{\alpha}(s) \\ I_{\beta}(s) \end{bmatrix} \tag{2}

となって，α軸とβ軸の干渉はなくなります。ただし，

Z^{s}(s) = Z_{s}(s) - Z_{m}(s) \tag{3}

であり，上付きsはstationary = 静止座標系を意味します。これを角速度ω1で回転するd-q座標系で表現した場合，次式となります。

\begin{bmatrix} V_{d}(s) \\ V_{q}(s) \end{bmatrix} = \begin{bmatrix} Z_{d}(s) & -Z_{q}(s) \\ Z_{q}(s) & Z_{d}(s) \end{bmatrix} \begin{bmatrix} I_{d}(s) \\ I_{q}(s) \end{bmatrix} \tag{4}

ここで，若干数学的に厳密な書き方ではありませんが，

\begin{align} Z_{d}(s) &= \mathrm{Re} \left [ Z^{s}(s + j \omega_{1}) \right ] \\ Z_{q}(s) &= \mathrm{Im} \left [ Z^{s}(s + j \omega_{1}) \right ] \end{align} \tag{5}

と表現できます。ここで，上式で実部と虚部を取り出すReとImを計算する場合，ひとまずラプラス演算子sを実数と見なして計算します。

文献紹介

d-q座標系での任意のインピーダンスの表現方法については，以下の論文にて詳しく述べられております。以下，この論文を「Harnefors氏の論文」と呼ぶことにします。

ieeexplore.ieee.org

Harnefors氏の論文では，上記の(5)式よりも数学的に厳密な書き方をしている他，三相不平衡な場合についての取り扱いについても記載しております。 ただし，三相不平衡な場合については筆者の理解が追い付いていないため煩雑になるため，まず本記事では三相平衡な場合に限定します。

以降，RL直列回路とLCL回路の場合を例に挙げながら，実際に計算していきます。

続きを読む

Raspberry PiとFlaskでスマートリモコンを作る

Raspberry Pi 電子工作 Matplotlib Python Flask

はじめに

今回のテーマは「家族の役に立つ工作」です。

泊まりがけの旅行から帰宅すると，戸締まりした部屋の中は真夏であれば蒸し風呂のように，また，真冬であれば冷蔵庫のように感じられます。帰宅する数時間前に外出先からエアコンを運転できれば，家に着く頃には快適な温度になっていることでしょう。このような機能を持ったエアコンも製品化されているようですし，そうでなくても，赤外線リモコンの信号を利用してこのような機能を後付けするための「スマートリモコン」と呼ばれるデバイスが今日では低価格で手に入ります。

また，昨今では天井照明がかつてのような「ヒモ」（正式にはプルスイッチ）ではなく，リモコンで制御されることが多くなって参りましたが，これを留守中にも点灯・消灯することにより，あたかも人が在宅しているかの如く見せ掛けて空き巣防止を狙うこともできます。

今回は，勉強とリハビリを兼ねてRaspberry Pi Zero WHでエアコンと天井照明を制御するスマートリモコンを自作した経緯や結果について書きます。なお，有名なirrp.pyに頼るのではなく，赤外線リモコン信号を送信するためのPythonモジュールを自作しました（と言ってもpigpioモジュールには同じように依存します）。これによって，例えばエアコンの温度設定や運転モードの組み合わせに応じたフレームを生成して送信することができます（学習リモコンの場合，すべての組み合わせを予めひとつひとつ学習しておく必要がありますよね）。

また，インターネットからエアコンや天井照明を制御するため，Flaskと呼ばれるフレームワークを用いたWebアプリを作りました。

なお，Raspberry Piを使って学習リモコンやスマートリモコン，スマートハウスを作るという先行研究は非常に多いので*1，本記事のオリジナリティを問われると「ぐぬぬ…」となってしまいますが，敢えて言えば自作モジュールirxmit.pyの部分でしょうか。また，pigpioライブラリ（モジュール）のwaveform機能について基礎的な部分のみですが，和文で書いた記事は少ないのではないかと推測しています。

本記事の最後に，赤外線LEDを6個に増量してパワーアップを図ったVersion 2 (V2)の製作過程について追記しました。十字配線ユニバーサル基板を使った例としてももしかすると参考になるかもしれません。

完成イメージ

まず，結論として完成したスマートリモコンのハードウェアとそのWebアプリのイメージを記載します。趣味なので仕様書のようなものは作っていませんが，作り始める前に思い描いていたものに近い形で実装することができました。

ハードウェア

図1に自作スマートリモコンHATを取り付けたRaspberry Pi Zero WHを示します。秋月電子通商のRaspberry Pi Zero用ユニバーサル基板に手持ちの部品を実装して組み上げました。温湿度センサ，赤外線LEDとそれを駆動するパワーMOSFET，赤外線受光モジュールを搭載しております。思ったよりコンパクトに収まっておりますが，望ましくはOLEDディスプレイやキャラクタLCDなど，表示デバイスを搭載できていたらと反省しております。温湿度センサの周りに抵抗器が密集しているのもあまりよろしくないですね。今年中にこれをプリント基板として作り直してみたいと密かに思っています。

図1: 自作スマートリモコンHATを取り付けたRaspberry Pi Zero WH

Webアプリ

図2に完成したWebアプリのキャプチャ画面を示します。ブラウザ上での見た目はBootstrap 4によってそれっぽく仕上がっております。この裏ではFlaskを用いたPythonのプログラムや，自作の赤外線リモコン信号送信モジュールなどが動いております。また，温湿度のトレンドグラフは馴れ親しんだ（？）Matplotlibで描いております。なお，自宅の家電機器を赤の他人に勝手に制御されると困るので，一応はパスワードで保護しています。これはFlaskのsessionとcookieを使って実現しています。

図2: Webアプリ画面キャプチャ

実際に動いている様子は動画の方が分かりやすいと思いますので，図3にYouTubeの動画を貼りつけます。

www.youtube.com
図3: FlaskによるWebアプリの動作の様子

おぉ…自分で作ったWebアプリとは思えないような美しさ…? 何もかもFlaskとBootstapのお陰ですね💦 たまに温湿度が欠測しているのが玉に瑕です。もう少し工夫が必要かもしれませんね。以下，開発の経緯や具体的な回路，ソフトの作り方について書いていきます。

*1:例えば赤外線LEDを積んだRaspberry Piでスマートルームをつくる(chatbot実装編) | うたかたサバイバーなど。

続きを読む

Raspberry PiでBLDCモータを120°通電制御する

パワエレ モータ制御 Raspberry Pi

はじめに

本ブログの前々記事にてブラシレスDC (BLDC)モータの120°通電制御について書きました。STMicroelectronicsのモータ制御開発キットP-NUCLEO-IHM001の制御装置はもちろんSTM32マイコンの載ったNucleo F302R8です。

このNucleoをRaspberry Piに繋ぎ替えてみたら面白いのではないか──。技術的にはほとんど意味のないことかもしれませんが，Raspberry Pi OS（旧称Raspbian）というLinuxが走っている「コンピュータ」であるRaspberry Piでも交流モータを制御できるのかどうか*1，試してみたいじゃありませんか？

インターネット上にはRaspberry PiでBLDCモータを廻している記事も散見されますが，ESC (electronic speed controller)を途中に挿入している場合が多いようです。モータドライブインバータへのゲートパルスをRaspberry Piで直接生成している例は見つけられませんでしたので，ひょっとすると新しい試みなのかもしれません*2。図1にRaspberry Piを接続したBLDCモータ実験システムを示します。

f:id:s-inoue2010:20210321031911j:plain:w400
図1: Raspberry Piを接続したBLDCモータ実験システム

結論として，C言語でpigpioライブラリ（pigpiod_if2.hではなくpigpio.h）を使い，かつ，マルチコアのRaspberry Piでは，Nucleoに肉薄する制御を（120°通電制御という範囲内で）実現できるらしいことが判りました。

当初はPythonで交流モータを廻すという異端を試してみたかったのですが，今のRaspberry Piではまだ難しいようです*3。 以下では，Raspberry Piでモータを廻すに至る経緯や，プログラム実装への試み，実装したC言語プログラム，シングルコアと4コアの比較について述べます。

*1:Raspberry PiにリアルタイムOSを載せたり，ベアメタルで使うという記事もありますね。いろんな技術を試してみたいものです。

*2:まぁ，Nucleoみたいな普通のマイコンボードでやろうよという制御ではあるのですが…。

*3:Raspberry Pi PicoのMicroPythonだったらできるかも──などと変なことを考えています。

続きを読む