メディアチーム【突入電流測定編】Pythonによるオシロスコープ波形データ解析の秘訣
メディアチームIT技術
Pythonで突入電流波形解析してみた！
今回は、オシロスコープで得られた csv の突入電流波形データに対して、Python で自動評価するテクニックを紹介します。
電子部品ヒューズの選定手順に含まれる、I2t-T カーブへの適用確認をジュール積分計算とグラフ描画によって行います。
電源回路の評価で必要不可欠な突入電流波形解析を Python で行います。
記事の最後にコピペで使えるサンプルコードも付属しています！
突入電流波形について
突入電流とは、電気回路の電源オン時に瞬間的に流れる電流となっており、短い時間で急速に増加したあと急速に減衰し、定常値に収束するような挙動をとります。
主に回路内のコンデンサ部品の合計容量によってその大きさは変わりますが、一般的には机上計算が難しい複雑な振る舞いとなります。
そのため、オシロスコープを使って実測値を測定しなければ評価することはできません。
突入電流を評価する意義
では、突入電流を評価する意義は何なのでしょうか。
最もポピュラーなケースが電子部品ヒューズの選定です。
過電流を検出し、安全のため自動で溶断するヒューズは、瞬間的な突入電流で誤って溶断してしまう可能性があります。
そのリスクを判断するためには突入電流の波形評価が必要不可欠なのです。
突入電流波形の評価方法
評価の具体的な方法としては、各社ヒューズメーカーが詳しい説明資料を用意しています。
ここでは、評価の中でも少々ややこしい I2t-T 特性の評価について簡単に説明していきます。
下図は Python で生成した大雑把な突入電流波形の例です。
まず、得られた突入電流波形の始めから時間単位でジュール積分値を計算します。
積分値は各時間ごとに累積させていき、横軸を時間、縦軸を累積積分値としたグラフを描きます。
このグラフは I2t-T 曲線と呼ばれ、この曲線がメーカー側が提示する範囲に収まっていれば合格となります。
各ヒューズのモデルごとに I2t-T 曲線は用意されていますので、エンジニアはヒューズ本来の溶断電流値を満たしつつ、合格範囲内となるモデルを選定していくことになります。
I2t-T 曲線は例えば、下記のメーカーのヒューズモデルのデータシートを閲覧することで確認できます。
SOC 株式会社：DC35V 3216サイズ 速断 表面実装型プロテクター
https://www.socfuse.com/ja/products/detail/dc35vp11cf/
本記事での実施内容
本記事では、上記で述べた I2t-T 特性の評価に注目し、Python で突入電流波形から I2t-T 曲線を描くことを行います。
突入電流波形は今までのように、オシロスコープから Pyvisa などを使って得られたことを想定し、csv 形式のものを用意します。
今回は、Python を使って突入電流を擬似的に再現した波形データを取り扱います。
コード内では、csv の読み出しに始まり、逐次処理方式で突入電流の始まりを自動検出し、ジュール積分値を累積する形で計算します。
最後に matplotlib を使い、グラフを作成する形です。
準備
今回のコード実行に必要な環境は下記になります。
インストールが必要なもの
- Python の開発環境（Spyderなど）
- Pandas モジュール（波形データ生成に必要）
- math モジュール（波形データ生成に必要）
- random モジュール（波形データ生成に必要
筆者の実行環境
- Python（3.6）
- Spyder（3.2.6）
解析する波形データ
波形データは Python で簡単に生成できますので、下記のコードを実行してみてください。
上の方にある突入電流波形の例通りのものが生成されるはずです。
1import pandas as pd
2import random
3import math
4import matplotlib.pyplot as plt
5
6X_SAMPLE = 3000
7NOISE_SEED_OFF = 123
8
9#rise_wave
10RISE_TIME=300
11PEAK_TIME = 400
12DAMP_TIME = 1200
13PEAK_CUR = 3.2
14STABLE_CUR =1.5
15NOISE_AMP_2=0.05
16RC=300
17PEAK_INT = (PEAK_TIME-RISE_TIME)*4
18DAMP_INT = (DAMP_TIME-PEAK_TIME)*4
19
20xx = [ i for i in range(3000)]
21
22analog_wave_rush=[]
23noise_rush=[]
24for i in range(len(xx)):
25 random.seed(NOISE_SEED_OFF+i)
26 noise_rush.append(NOISE_AMP_2*random.randint(-100,100)/100)
27
28 if i<RISE_TIME:
29 analog_wave_rush.append(noise_rush[i])
30 elif i<PEAK_TIME:
31 rush_dt = PEAK_CUR*math.sin(2/PEAK_INT*xx[i-RISE_TIME]*math.pi)
32 analog_wave_rush.append(rush_dt+noise_rush[i])
33 elif i<DAMP_TIME:
34 damp_dt = PEAK_CUR - (PEAK_CUR-STABLE_CUR)*math.cos(2/DAMP_INT*xx[i-PEAK_TIME]*math.pi + math.pi/2*3)
35 analog_wave_rush.append(damp_dt+noise_rush[i])
36 else:
37 rise_dt = STABLE_CUR*(1 - math.exp(-1/RC*xx[i-RISE_TIME]))
38 analog_wave_rush.append(STABLE_CUR+noise_rush[i])
39
40plt.figure(figsize=(10,8))
41plt.plot(xx,analog_wave_rush)
42plt.savefig("./analog_wave_rush.jpg")
43plt.show()
44
45wave_data= pd.DataFrame({"wave_rush":analog_wave_rush})
46wave_data.to_csv("analog_wave_rush.csv")
コード本文
それでは、解析をスタートしましょう！
今回も波形データに特化した時系列逐次処理で処理を行っていきます。
詳細については前回記事（記事ID26995へのリンク）を参照してみてください。
①モジュールのインポートとcsvデータの読み込み
まずは下記のようにモジュールのインポートと csv の読み込みを行います。
ここは前回記事とほとんど変わりません。
1#変数の定義
2ch1_rise_st=[]
3ch1_rise_st_buf =[ 0 for i in range(10)]
4TimeCnt=[]
5JouleIntg=[]
6prevIntg=0
7
8#立ち上がり検出を1回だけ実施するためのフラグ
9rise_st_f = True
10start_measure_f = False②使用する変数の定義とフラグの初期化
次に変数とフラグを定義して初期化します。
使われ方については後述のコードで説明していきます。
1#変数の定義
2ch1_rise_st=[]
3ch1_rise_st_buf =[ 0 for i in range(10)]
4TimeCnt=[]
5JouleIntg=[]
6prevIntg=0
7
8#立ち上がり検出を1回だけ実施するためのフラグ
9rise_st_f = True
10start_measure_f = False③波形データの時系列逐次処理
前回から使っている時系列逐次処理の部分です。
こちらも従来通りですので、説明は割愛します。
メインは次から説明する while 文の中身です。
1#データの読み込み開始
2##最初の１行目はループ前に読み込む
3row_data = wave_file.readline()
4split_data = row_data.split(",")
5
6dt_cnt=0
7while row_data:
8 index_dt = float(split_data[0])
9 ch1_dt = float(split_data[1])
10 ###########この間に解析処理を入れる##########
11～～～～～
12
13 #######################################
14
15 ##次の行を読み込む。ここで最終行なら空白が返ってくるためループ終了
16 dt_cnt+=1
17 row_data = wave_file.readline()
18 split_data = row_data.split(",")
19wave_file.close()④突入電流開始時間の検出
メインの処理の前に突入電流の検出を行います。
あらかじめ設定した閾値を超過した立ち上がりタイミングから検出処理を開始します。
コードを見ていきましょう！
1 #測定開始の判定
2 #バッファを使って閾値の超過タイミングを検出する
3 ##サンプル間隔を開けてノイズの影響を軽減する
4 if index_dt % DET_SMP_INT == 0 and rise_st_f:
5 #リングバッファへの保存
6 for i in range(len(ch1_rise_st_buf)-1):
7 ch1_rise_st_buf[len(ch1_rise_st_buf)-1-i] = ch1_rise_st_buf[len(ch1_rise_st_buf)-1-i-1]
8 ch1_rise_st_buf[0]=ch1_dt
9 #5回連続で値が閾値を超えた場合に閾値超過と判定する
10 for i in range(len(ch1_rise_st_buf)):
11 if ch1_rise_st_buf[i] > START_TH:
12 if i >4 :
13 start_measure_f = True
14 rise_st_f=False
15 StartPnt = index_dt
16 break
17 else:
18 break今までと同様に、リングバッファを使って、連続して閾値を超えたタイミングを検出しています。
検出時には StartPnt に開始時のサンプル位置を記録し、グラフ描画のスタート地点としています。
⑤ジュール積分値の計算とグラフ描画用の配列作成
メインとなるジュール積分値の計算です。
基本的には突入電流の検出タイミングから指定したサンプル間隔ごとに計算していく形になります。
コードを見ていきましょう！
1 #ジュール積分の計算をする処理
2 if index_dt % DET_SMP_INT == 0 and start_measure_f:
3 #検出タイミングを0secとしてサンプル間隔毎にグラフ用の配列に追加する
4 TimeCnt.append((index_dt - StartPnt)*SMP_TIME)
5 #現在の時間でI2tを計算し、前回の値を加算することでジュール積分値を算出し、グラフ用の配列に追加する
6 integ_calc = ch1_dt*ch1_dt*(SMP_TIME*DET_SMP_INT) + prevIntg
7 JouleIntg.append(integ_calc)
8 #ジュール積分値は前回の値を累積するので今回の値を保持しておく
9 prevIntg=integ_calcグラフを描画するのに必要なリストデータは時間 TimeCnt と積分値 JouleIntg なのでこの２つに要素を追加していく形になります。
まず、
1TimeCnt.append((index_dt - StartPnt)*SMP_TIMEで実時間を検出タイミングを起点に算出します。
ここで、波形データ1点ごとのサンプリング時間 SMP_TIME を正しく把握しておくことが重要です。
次に
1integ_calc = ch1_dt*ch1_dt*(SMP_TIME*DET_SMP_INT) + prevIntg前回時間におけるジュール積分値に現在時間における積分値を加算して累積していきます。
積分値の計算は、離散データでは
時間間隔 × 電流値の2乗
で計算ができるため、データのサンプリング時間（SMP_TIME）と処理のサンプリング間隔（DET_SMP_INT）から正しい値を算出しましょう。
上記を繰り返し処理していくことでグラフ描画用のデータが生成されます。
⑥ I2t-t グラフの描画
最後に波形データ処理の結果得られたジュール積分値を時間軸にプロットします。
ポイントとしては、比較するメーカーのグラフと描画方法と範囲を揃えることです。
メーカーのグラフは図データのみになるので、最終的には目視で比較することとなります。
1#I2t-tグラフを描画する
2plt.xscale("log")
3plt.yscale("log")
4plt.grid(which='major')
5plt.grid(which='minor')
6plt.xlim(GRAPH_XLIM)
7plt.ylim(GRAPH_YLIM)
8plt.scatter(TimeCnt, JouleIntg,s=3)
9plt.show()下図のグラフが生成されました。
このグラフとメーカーのグラフを見比べてみましょう。
こちらのリンクの pdf 形式のデータシートからグラフを参照することができます。
プロットしたデータよりマージンがある曲線は3150のようです。
したがって、型番 DC35VP11CF では定格電流3.15A を選べばよいことが分かりました。
サンプルスクリプト
最後に今回紹介した解析コードの全文を下記に記載します。
今までの説明で省かれていた閾値などの設定値も記載しており、解析したい対象に応じて値を変えれば、同様の解析が手持ちの波形データでも行えるはずです。
1import matplotlib.pyplot as plt
2
3wave_path = "./analog_wave_rush.csv"
4wave_file = open(wave_path,'r')
5skip_row_count = 1
6#設定値関連
7SMP_TIME = 100 * 10**(-6) #取得した波形データのデータごとのサンプリング時間（今回は100us)
8START_TH = 0.1 #突入電流の検出を開始する電流閾値[A]
9DET_SMP_INT = 2 #電圧閾値を検出する際にデータ比較をするサンプル間隔
10GRAPH_XLIM = [100 * 10**(-6),1.0] #I2t-tグラフのX軸の描画範囲（比較するメーカーのグラフに合わせる）
11GRAPH_YLIM = [10**(-2),10**2] #I2t-tグラフのY軸の描画範囲（比較するメーカーのグラフに合わせる）
12
13#変数の定義
14ch1_rise_st=[]
15ch1_rise_st_buf =[ 0 for i in range(10)]
16TimeCnt=[]
17JouleIntg=[]
18prevIntg=0
19
20#立ち上がり検出を1回だけ実施するためのフラグ
21rise_st_f = True
22start_measure_f = False
23
24#ヘッダ行等のスキップ
25while skip_row_count !=0:
26 row_data = wave_file.readline()
27 skip_row_count-=1
28
29#データの読み込み開始
30##最初の１行目はループ前に読み込む
31row_data = wave_file.readline()
32split_data = row_data.split(",")
33
34dt_cnt=0
35while row_data:
36 index_dt = float(split_data[0])
37 ch1_dt = float(split_data[1])
38 ###########この間に解析処理を入れる##########
39
40 #測定開始の判定
41 #バッファを使って閾値の超過タイミングを検出する
42 ##サンプル間隔を開けてノイズの影響を軽減する
43 if index_dt % DET_SMP_INT == 0 and rise_st_f:
44 #リングバッファへの保存
45 for i in range(len(ch1_rise_st_buf)-1):
46 ch1_rise_st_buf[len(ch1_rise_st_buf)-1-i] = ch1_rise_st_buf[len(ch1_rise_st_buf)-1-i-1]
47 ch1_rise_st_buf[0]=ch1_dt
48 #5回連続で値が閾値を超えた場合に閾値超過と判定する
49 for i in range(len(ch1_rise_st_buf)):
50 if ch1_rise_st_buf[i] > START_TH:
51 if i >4 :
52 start_measure_f = True
53 rise_st_f=False
54 StartPnt = index_dt
55 break
56 else:
57 break
58
59 #ジュール積分の計算をする処理
60 if index_dt % DET_SMP_INT == 0 and start_measure_f:
61 #検出タイミングを0secとしてサンプル間隔毎にグラフ用の配列に追加する
62 TimeCnt.append((index_dt - StartPnt)*SMP_TIME)
63 #現在の時間でI2tを計算し、前回の値を加算することでジュール積分値を算出し、グラフ用の配列に追加する
64 integ_calc = ch1_dt*ch1_dt*(SMP_TIME*DET_SMP_INT) + prevIntg
65 JouleIntg.append(integ_calc)
66 #ジュール積分値は前回の値を累積するので今回の値を保持しておく
67 prevIntg=integ_calc
68
69 #######################################
70
71 ##次の行を読み込む。ここで最終行なら空白が返ってくるためループ終了
72 dt_cnt+=1
73 row_data = wave_file.readline()
74 split_data = row_data.split(",")
75wave_file.close()
76
77#I2t-tグラフを描画する
78plt.xscale("log")
79plt.yscale("log")
80plt.grid(which='major')
81plt.grid(which='minor')
82plt.xlim(GRAPH_XLIM)
83plt.ylim(GRAPH_YLIM)
84plt.scatter(TimeCnt, JouleIntg,s=3)
85plt.savefig("./I2t-t_graph.jpg")
86plt.show()こちらの記事もオススメ！
ライトコードでは、エンジニアを積極採用中！
ライトコードでは、エンジニアを積極採用しています！社長と一杯しながらお話しする機会もご用意しております。そのほかカジュアル面談等もございますので、くわしくは採用情報をご確認ください。
採用情報へ
「好きを仕事にするエンジニア集団」の（株）ライトコードです！ ライトコードは、福岡、東京、大阪、名古屋の４拠点で事業展開するIT企業です。 現在は、国内を代表する大手IT企業を取引先にもち、ITシステムの受託事業が中心。 いずれも直取引で、月間PV数1億を超えるWebサービスのシステム開発・運営、インフラの構築・運用に携わっています。 システム開発依頼・お見積もり大歓迎！ また、現在「WEBエンジニア」「モバイルエンジニア」「営業」「WEBデザイナー」を積極採用中です！ インターンや新卒採用も行っております。 以下よりご応募をお待ちしております！ https://rightcode.co.jp/recruit
