ワットチェッカー的なものを自作する (その1): エンジニア徒然草

2014年8月 1日 (金)

ワットチェッカー的なものを自作する (その1)

前回の記事にワットチェッカーを買って冷蔵庫の消費電力を測った話を書いたが、その際「自分で作った方がたぶん安い」と思わなくもなかった。正確さなどをあまり欲張らなければマイクロコントローラーのAD変換機能などを使って割と簡単に実現できるはずだ。「ワットチェッカー自作」でGoogle検索すると実例がいくつも見つかる。

冷蔵庫の消費電力測定はそれなりに急いでいたので結局既製品を買ったが、これには経時変化を記録したりする機能が無いので、この辺の補完を目指して検討を進めてみた。

/*----------- -----------*/

電流・電圧信号の取り出し

ワットチェッカーでは抵抗を使って取り出しているようだ。

完成された製品なので回路全体を絶縁した状態に保てるからこれで良い。しかし自作する場合、途中の信号を調べる時に感電したり漏電遮断器をトリップさせてしまったりしないよう入口で絶縁してしまうのが安全である。なお上の回路の電圧信号は負荷電圧を正しく反映していない。電流信号、つまりRsによる電圧降下分が入っていないのだが、その振幅を十分小さく (100mV位) すれば無視できる。あるいは信号処理で引き算することも可能なはずだ。

交流信号を絶縁して取り出すにはトランスを使うのが簡単である。電圧信号には小型の電源トランス、電流信号には電流トランスを使う。

簡単な商用電源の極性確認回路を追加している。常用するのであればケースに収めたり保護用ヒューズを入れたりする必要があるが、実験用なのでバラックセットである。

分電盤などの既存配線の電流を測る場合はクランプ型の電流トランスを使うが、これは単価1千円程度と比較的高価である。今回のように後から配線を通す場合はモールド型が使えて、これは単価3百円程度で入手できる。なお高周波用の電流トランスがあり、これはさらに安価で売られている場合が多いが、最低周波数が数kHz以上と商用電源周波数では使えないので間違えて買わないよう注意する必要がある。

いわゆるユニバーサル基盤を使って組んでいるが、そのまま使ったのではAC100Vがかかる部分の絶縁に不安がある。接地線を除くACラインに接続された部分の周囲のランド1列をドリルでさらって絶縁距離を確保するようにした。

デジタルオシロを使って信号を確認するついでに電流、電圧、そして有効電力を計算させてみた。

定格電力10WのLED電球を点灯した場合の信号だ。この場合の元の信号の値を得るにはトランスのゲイン、つまり変換係数を使う必要がある。

電流ゲイン: Ki = 0.1V/A
電圧ゲイン: Kv = 0.11V/V

元の信号の値は以下のように計算できる。

電流RMS値: 0.0158V ÷ Ki = 0.158A
電圧RMS値: 11.4V ÷ Kv = 104V
有効電力: 0.0912VV ÷ Ki ÷ Kv = 8.29 AV = 8.29W
皮相電力: 電流RMS値 x 電圧RMS値 = 0.158 x 104 = 16.4VA
力率: 有効電力 ÷ 皮相電力 = 8.29 ÷ 16.4 = 0.505

校正していないので正確な数字ではないが、個人的に使う分には問題ない程度の確度は得られているはずだ。しかし計算された有効電力の数字は定格電力に比べて少し小さ過ぎるかもしれない。これは波形の画面上の振幅が小さ過ぎたことによる誤差が影響している可能性がある。

電流波形を見ると酷く歪んでいて多量の高調波成分を含んでいる事が分かる。このRMS値を正確に測るにはそれなりに速いサンプルレートが必要だ。電圧波形にもひずみが認められるが、その程度は僅かである。

信号の高調波成分を調べる

高調波成分を調べるのはデジタル信号処理を行う場合の必要なサンプルレートを見積もるためだ。あまり周波数の高い成分を捉える必要が無ければサンプルレートは低くて良い。

オシロにFFT機能が無いので、PCのオーディオ機能を利用して信号のスペクトラムを調べた。トランスの出力信号をPCのオーディオ入力につなぐにはレベルに調整する必要がある。

この回路なら電圧信号は約100mV RMSでPCに入力される。また電流トランスの負担抵抗の合成値が25Ωになるので、電流信号は4Aの時に100mVになる。まずWaveSpectraを使って電圧信号を調べてみた。

存在する高調波は奇数倍のものだけである。最もレベルの高い3倍と5倍の高調波が基本波に対して-30dbなので歪率は5%程度だ (-30 + 3 = -27db → 約1/20)。つまり歪を完全に無視しても、そのことによる誤差は5%程度である。また1kHzまでの高調波を捉えられれば、それより周波数の高い高調波を取りこぼすことによる誤差が1%に達することは無いはずだ。

その一方電流は、予想通り高調波が非常に多い。

こうなるとWaveSpectraでは歯が立たない。今欲しいのは、ある周波数よりも上にある高調波成分のパワーが占める割合である。WaveSpectraは手軽で使いやすいのだが、このような場合を想定して作られてはいない。こうなると自分でプログラムを組むのが結局は早道だったりする。

WaveSpectraでファイルに落としたサンプルをPythonで作ったプログラムで加工したものである。フーリエ変換を行ってパワースペクトラムを求め、全体の合計を1.0に正規化した。周波数の高い方からの累積曲線を見れば欲しかった答えが求められる。

例えば、無視する高調波電力の合計を1/10,000以下 (1/100の二乗) にすればそのことによる誤差は1%以下になる。上のグラフから、この条件だと8kHzまでの高調波は無視できないことが分かる。また2kHz以上の高調波を無視すると10%程度の誤差になる。

今話題にしている「誤差」が直接影響するのはRMS値である。有効電力については、基本的に電流波形の基本波成分が正しく捉えられていれば良いのだが、この「正しく」が曲者だ。捉える必要が無くても高い周波数の成分が存在する場合、マイクロコントローラーで良く使われているSAR (Successive Approximation Register、逐次比較) 方式のADCでサンプルレートを下げてしまうとエイリアス信号の影響が大きくなるからである。ただしエイリアス信号の大きさはパワースペクトラムの累積値ではなくサンプル周波数付近の高調波信号の大きさが影響する。

以上から今回試した電流波形をSAR方式のADCで測る場合の最低サンプルレートについて以下の事が言えそうである。

電流のRMS値の精度を10%以内にするには4ksps以上必要
電流のRMS値の精度を1%以内にするには16ksps以上必要
有効電力の精度を10%以内にするには900sps以上必要
有効電力の精度を1%以内にするには4ksps以上必要

手っ取り早く測ってみる

最終的にはマイクロコントローラーを使った仕組みを考えているが、とりあえずPCを使って測ってみることにした。サンプルレートは44.1kspsにしたので不足するおそれはない。後はプログラムを組めば済む話である。

# -*- coding: cp932 -*-
""" WattSimple.py: Realtime Watt Meter. """

import pyaudio
import time
import numpy as np
from WattParam import Cv, Ci

# Audio設定
FMT = pyaudio.paInt16
NCH = 2
SPS = 44100

# __main__で定義していないグローバル変数
tot_sec = 0
tot_kwh = 0.0

def callback(in_data, frame_count, time_info, status):
    global tot_sec, tot_kwh
    data = np.fromstring(in_data, np.int16)
    i_data = Ci * data[::2]     # 偶数番目のデータ (L ch) は電流
    v_data = Cv * data[1::2]    # 奇数番目のデータ (R ch) は電圧
    i = i_data.std()            # RMS電流の計算
    v = v_data.std()            # RMS電圧の計算
    va = i * v                  # 皮相電力の計算
    p = (i_data * v_data).sum() / i_data.size   # 有効電力の計算
    i = 0. if i < 0.01 else i   # 以下4行はノイズ対策
    va = 0. if i < 0.01 else va
    p = 0. if i < 0.01 else p
    pf = 0. if i < 0.01 or va < 0.1 else p/va
    pf = 1.0 if pf > 1. else pf # 誤差対策
    tot_sec += 1
    tot_kwh += p / 3600
    strtime = time.strftime("%Y/%m/%d %H:%M:%S")
    print ("[%s] %.3f Wh in %d:%02d:%02d  [Ctrl+C to stop]"
           % (strtime, tot_kwh,
              tot_sec//3600, (tot_sec%3600)//60, tot_sec%60) )
    print ("\t\t%5.1fV %6.3fA %6.1fVA %6.1fW (PF=%.3f)" % (v, i, va, p, pf) )
    return (in_data, pyaudio.paContinue)

if __name__ == "__main__":
    pa = pyaudio.PyAudio()
    st = pa.open(format = FMT, channels = NCH, rate = SPS,
                 input = True, output = False,
                 frames_per_buffer = SPS, stream_callback = callback)
    try:
        while True:
            time.sleep(0.1)
    except KeyboardInterrupt:           # Ctrl+Cで終了
        pass
    except:
        raise
    finally:
        st.stop_stream()
        st.close()
        pa.terminate()

手っ取り早く組むにはPythonが一番だ。標準装備ではないモジュールpyaudioとnumpyを使っているので必要なら別途インストールする。WattParamモジュールは校正で変更する可能性のある値を集めた一種の設定ファイルで、WattSimple.pyと同じフォルダーに入れておく。

# -*- coding: cp932 -*-
""" WattParam.py: Analog Parameter of Watt Meter. """

# 校正方法:
#   - Vvtin =とVvtout =の右辺を実測値に修正する
#   - ADJ = の右辺をVvtin =の右辺と同じ値にする
#   - WattSimple.pyを実行して表示される電圧RMS値が100近くになるよう録音レベル調整
#   - ADJ = の右辺を最終的に表示された電圧RMS値に修正する

ADJ = 110.1 

# アナログインターフェイス・パラメーター
Vvtin = 100                 # 電圧トランスの入力電圧実測値
Vvtout = 11                 # 電圧トランスの出力電圧実測値
Vdiv = 1. / (100 + 1)       # 電圧信号の減衰比
Nct = 1000                  # 電流トランスの二次巻き数
Rct = 25                    # 電流トランスの負担抵抗値
Kv = Vdiv * Vvtout / Vvtin  # 電圧信号経路のゲイン [V/V]
Ki = 1.0 * Rct / Nct        # 電流信号経路のゲイン [V/A]
Cv = 1. / ADJ               # 電圧換算係数: 入力換算で10mV/LSB相当
Ci = Cv * Kv / Ki           # 電流換算係数: 入力換算で436μA/LSB相当

DMM (Digital Multi Meter、デジタルテスターのこと) があれば、DMMと部品の誤差程度の確度に校正できる手順をコメントに記入した。録音レベルがプログラムから設定できればもう少しスマートにできたのだが、簡単な方法が見つからなかった。

WattSimple.pyを実行すると、Ctrl+Cが押されるまでコマンドプロンプトに消費電力などを1秒毎に延々と表示し続ける。

情報をprintする部分で、例えば1分毎とか5分毎にデータベースを更新するような処理を行うようにすると実用的なプログラムにつながるのではないかと思う。しかしデスクトップコンピューターは数十Wの電力を消費するし、ノートでも20W程度は消費しているから、定格電力10WのLED電球の消費電力を測るためだけにこのプログラムを動かすのはナンセンスだ。多分唯一有意義な使い道はこのプログラムを動かすコンピューターの消費電力を測ることだろう。デスクトップで動かしたときのCPU使用率は1%未満だったので、このプログラムを動かすことによる消費電力の増加は僅かである。

機能確認できたのでこれで済ませるはずだった。しかし実際に動かしてみるとプロットが欲しくなる。

ウインドウサイズを変更したときなどに若干不審な挙動を示す場合があるが、一応それらしくプロットできることが確認できた。

# -*- coding: cp932 -*-
""" WattScope.pyw: GUI version of Realtime Watt Meter. """

import pyaudio
import time, threading
import numpy as np
import matplotlib.pyplot as plt
import wx
from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg
from WattParam import Cv, Ci

POINTS = 300        # プロットのポイント数
PLTSTYLE = 'r.-'    # プロットのスタイル

# Audio設定
FMT = pyaudio.paInt16
NCH = 2
SPS = 44100

# __main__で定義していないグローバル変数
v = 0.0
i = 0.0
va = 0.0
pf = 0.0
p = 0.0
tot_sec = 0
tot_kwh = 0.0
wav = []
waa = []

class MyFrame(wx.Frame):
    def __init__(self):
        super(MyFrame,self).__init__(None, wx.ID_ANY, size=(800, 600) )
#        self.SetIcon(wx.Icon('WattScope.ico', wx.BITMAP_TYPE_ICO) )
        self.SetTitle("Watt Scope")
        self.values = []
        self.times = []
        self.bg = None
        self.lock = threading.Lock()
        self.fig = plt.Figure()
        self.ax = self.fig.add_subplot(111, xlim=(0, POINTS) )
        self.canv = FigureCanvasWxAgg(self, wx.ID_ANY, self.fig)
        self.line, = self.ax.plot([], [], PLTSTYLE, animated=True,
                                  label="Power [W]")
        self.wava, = self.ax.plot([], [], 'g', animated=True,
                                  label="Current Waveform", alpha=0.5)
        self.wavv, = self.ax.plot([], [], 'y', animated=True,
                                  label="Voltage Waveform", alpha=0.5)
        self.ax.get_yaxis().set_animated(True)
        self.ax.legend(loc="upper left")
        self.ax.set_xlabel("Time [sec]")
        self.ax.set_ylabel("Power Consumption [W]")
        self.time_template = '0 sec = %s'
        self.time_text = self.ax.text(0.01, 0.02, '', animated=True,
                                      transform=self.ax.transAxes)
        self.v_template = '%5.1f V'
        self.v_text = self.ax.text(0.85, 0.9, 'v', animated=True,
                                   transform=self.ax.transAxes)
        self.a_template = '%6.3f A'
        self.a_text = self.ax.text(0.85, 0.85, 'a', animated=True,
                                   transform=self.ax.transAxes)
        self.va_template = '%6.1f VA'
        self.va_text = self.ax.text(0.85, 0.8, 'va', animated=True,
                                    transform=self.ax.transAxes)
        self.p_template = '%6.1f W'
        self.p_text = self.ax.text(0.85, 0.75, 'p', animated=True,
                                   transform=self.ax.transAxes)
        self.pf_template = 'PF=%.3f'
        self.pf_text = self.ax.text(0.85, 0.7, 'pf', animated=True,
                                    transform=self.ax.transAxes)
        self.kwh_template = '%.3f Wh in %d:%02d:%02d'
        self.kwh_text = self.ax.text(0.7, 0.95, 'kwh', animated=True,
                                     transform=self.ax.transAxes)
        self.ax.grid(True)
        self.Bind(wx.EVT_SIZE, self.on_size)
        self.Bind(wx.EVT_ICONIZE, self.on_size)
        self.Bind(wx.EVT_MAXIMIZE, self.on_size)

    def on_size(self, event):
        event.Skip()
        self.bg = None

    def animation(self):
        if self.bg == None:
            self.lock.acquire()
            self.fig.canvas.draw()
            self.bg = self.canv.copy_from_bbox(self.ax.get_figure().bbox)
            self.lock.release()
        else:
            self.canv.restore_region(self.bg)
        self.values.append(p)
        t = time.time()
        self.times.append(t)
        if len(self.values) > POINTS:
            self.values = self.values[-POINTS:]
            self.times = self.times[-POINTS:]
        atimes = np.array(self.times) + POINTS - t - 1
        time0 = time.localtime(t - POINTS + 1)
        self.time_text.set_text(self.time_template %
                                time.strftime("%Y/%m/%d %H:%M:%S", time0) )
        self.v_text.set_text(self.v_template % v)
        self.a_text.set_text(self.a_template % a)
        self.va_text.set_text(self.va_template % va)
        self.p_text.set_text(self.p_template % p)
        self.pf_text.set_text(self.pf_template % pf)
        self.kwh_text.set_text(self.kwh_template % (tot_kwh,
                                                    tot_sec//3600,
                                                    (tot_sec%3600)//60,
                                                    tot_sec%60) )
        maxvalue = 1.2 * max(self.values)
        ytop = 5000
        for x in [100, 200, 500, 1000, 2000]:
            if x > maxvalue:
                ytop = x
                break
        self.line.set_data(atimes, self.values)
        wa = (0.5 + 0.4 * waa
              / (10 if a < .01 else max(1, waa.max() - waa.min() ) ) ) * ytop
        wv = (0.5 + wav / 500 ) * ytop
        self.wava.set_data(np.arange(waa.size), wa)
        self.wavv.set_data(np.arange(waa.size), wv)
        self.ax.set_xlim([0, POINTS] )
        self.ax.set_ylim([0, ytop] )

        self.ax.draw_artist(self.ax.get_yaxis() )
        self.ax.draw_artist(self.time_text)
        self.ax.draw_artist(self.v_text)
        self.ax.draw_artist(self.a_text)
        self.ax.draw_artist(self.va_text)
        self.ax.draw_artist(self.p_text)
        self.ax.draw_artist(self.pf_text)
        self.ax.draw_artist(self.kwh_text)
        self.ax.draw_artist(self.line)
        self.ax.draw_artist(self.wava)
        self.ax.draw_artist(self.wavv)
        
        self.canv.blit(self.ax.clipbox)

def callback(in_data, frame_count, time_info, status):
    global v, a, va, pf, p, tot_sec, tot_kwh, waa, wav
    data = np.fromstring(in_data, np.int16)
    a_data = Ci * data[::2]     # 偶数番目のデータ (L ch) は電流
    v_data = Cv * data[1::2]    # 奇数番目のデータ (R ch) は電圧
    t100ms = SPS / 10
    if POINTS < t100ms:
        waa = np.array(a_data[:t100ms:t100ms//POINTS] )
        wav = np.array(v_data[:t100ms:t100ms//POINTS] )
    else:
        waa = np.array(a_data[:POINTS] )
        wav = np.array(v_data[:POINTS] )
    a = a_data.std()            # RMS電流の計算
    v = v_data.std()            # RMS電圧の計算
    va = a * v                  # 皮相電力の計算
    p = (a_data * v_data).sum() / a_data.size   # 有効電力の計算
    a = 0. if a < 0.01 else a   # 以下4行はノイズ対策
    va = 0. if a < 0.01 else va
    p = 0. if a < 0.01 else p
    pf = 0. if a < 0.01 or va < 0.1 else p/va
    pf = 1.0 if pf > 1. else pf # 誤差対策
    tot_sec += 1
    tot_kwh += p / 3600
# Debug用 -- 拡張子を.pyに変更する必要あり
##    print ("[%s] %5.1fV %6.3fA %6.1fVA %6.1fW (PF=%.3f)"
##           % (time.strftime("%Y/%m/%d %H:%M:%S"),
##              v, a, va, p, pf) )
    frm.animation()             # アニメーションを1フレーム描画
    return (in_data, pyaudio.paContinue)

if __name__ == "__main__":
    # sources: http://stackoverflow.com/a/1552105/674475
    import ctypes
    myappid = 'mitaka1954.wattscope.wattscope.000' # arbitrary string
    ctypes.windll.shell32.SetCurrentProcessExplicitAppUserModelID(myappid)

    wxa = wx.App(redirect=False)
    frm = MyFrame()
    pa = pyaudio.PyAudio()
    st = pa.open(format = FMT, channels = NCH, rate = SPS,
                 input = True, output = False,
                 frames_per_buffer = SPS, stream_callback = callback)
    try:
        frm.Show(True)
        wxa.MainLoop()          # 終了するにはプロットウインドウを[X]で閉じる
    except KeyboardInterrupt:   # コマンドプロンプトでCtrl+Cが押されても正常終了する
        pass
    except:
        raise
    finally:
        st.stop_stream()
        st.close()
        pa.terminate()

少し余計な処理も行っているが、私のデスクトップ環境でこのプログラムを動かした場合のCPU使用率は3%程度である。

今回のまとめ

ワットチェッカーは電圧や電流を測るのだが、その内容はオーディオ周波数帯域の信号処理だと考えた方が実情に合っている。

電流信号と電圧信号を適切に前処理すれば、ワットチェッカー的なものはPCのオーディオ機能を使って比較的簡単に実現できる。しかしPCの計算処理能力のごく一部だけしか使わないので、ワットチェッカー的なもの専用にPCを使うのはエネルギー効率がとても悪い。やはりマイクロコントローラーなどもっと規模の小さい計算機を使うべきだ。これについて引き続き検討する。

2014年8月 1日 (金) パソコン・インターネット, 趣味 | 固定リンク
Tweet