ThingSpeak を使用したサーボのギア列のライブ RUL 推定

この例では次を使用します。

この例では、Arduino ベースのデータ収集システムから ThingSpeak™ へ、また ThingSpeak から MATLAB® で実行中の RUL 推定エンジンへのサーボモーターデータのリアルタイムストリーミングを通して、サーボモーターのギア列の残存耐用期間 (RUL) を推定する方法を説明します。

ホビーサーボモーターを駆動する電流信号のモーター電流シグネチャ解析 (MCSA) を使用して、モーターとギア列の故障を示すいくつかの周波数の関心領域から、周波数領域 (スペクトル) の特徴を抽出します。特徴の組み合わせを使用して、以降の RUL 推定に用いる健康インジケーター (HI) を構成します。

MCSA は、サーボのギア列でトルクまたは回転数の変動を誘発し、相関しているモーター電流の変化につながる故障を診断するのに役立つ手法です。MCSA はモーターの電流信号のみが解析に必要とされ、高価な検出ハードウェアが追加で必要にならないため、モーターの故障解析に最適であると実証されています。特に、加速度計や他の振動センサーをもつ計器でのギア列へのアクセスが容易ではない場合、従来の振動センサーを使用したギア故障の検出は困難です。

この例では、教育的な実験室実習や産業用途のプロトタイピングに適した単純な市販のコンポーネントを使用して、リアルタイムのデータストリーミング、特徴抽出、および RUL 推定システムを作成する方法を説明します。ハードウェアなしでこの例を実行するには、関数 sendSyntheticFeaturesToThingSpeak を使用して合成データを生成します。

データストリーミングシステムと RUL 推定エンジンを構成するための簡略化されたワークフローには、以下の手順が含まれます。

市販のコンポーネントを使用して、ハードウェアとデータ収集システムを開発します。
リアルタイムデータを、Arduino® UNO のようなマイクロコントローラーから MATLAB へとストリーミングします。
MATLAB でリアルタイムデータを処理して特徴を抽出し、ThingSpeak を使用してその特徴をクラウドにストリーミングします。
クラウドデータを MATLAB の RUL 推定および可視化エンジンへとストリーミングします。

ハードウェアセットアップとデータ収集

ハードウェアの概要

この例では、Futaba の標準的なホビー用サーボ S3003 を連続回転用に改造したものから電流データを収集しました。サーボは、内部の DC モーターの高い回転数を、出力シャフトの高いトルクに変換します。これを実現するため、サーボは DC モーター、ナイロン製または金属製ギアのセット、および制御回路から構成されています。その制御回路を取り外して、5 V の定電圧を DC モーターのリード線に直接印加できるようにし、DC モーターを通る電流を電流検出抵抗器により直接監視できるようにしています。

反対側にある 2 つ目のサーボは、その DC モーターターミナルが低抵抗の抵抗器を介して互いに短絡されており、これを使用して駆動サーボに対する高トルクの逆負荷を生成することで、ギア列の劣化を速め、MCSA 解析のためにモーター電流信号に大きな変動を引き起こします。

サーボモーターとギア列

下の図に示すように、Futaba の S3003 サーボは、4 組の噛み合うナイロン製ギアで構成されています。DC モーターシャフト上のピニオン P1 は段付きギア G1 と噛み合います。ピニオン P2 は段付きギア G1 の成形部分であり、段付きギア G2 と噛み合います。ピニオン P3 はギア G2 の成形部分であり、段付きギア G3 と噛み合います。ピニオン P4 は G3 とともに成形されており、出力シャフトに取り付けられた最終ギア G4 と噛み合います。G1 と P2、G2 と P3、G3 と P4 の段付きギアの組は自由回転ギアです $-$ すなわち、これらは各シャフトに固定されていません。ギアのセットは 278:1 の減速を行い、DC モーターを 5 V で駆動する場合、約 6117 rpm のモーター回転数を出力シャフトで約 22 rpm に減速します。

次の表に、歯数と、各ギア噛み合いにおける出力速度、ギアの噛み合い周波数、累積ギア減速比の理論値を示します。

サーボの出力シャフトの回転速度は、赤外線フォトインタラプタを、直径 40 mm の 3-D プリントによる 16 スロットタコメーターホイールとともに使用して測定されました。ホイールの 16 個のスロットは等間隔であり、IR フォトインタラプタはホイールの 1 回転につき正確に 16 パルスとなるよう配置されました。サーボとフォトインタラプタは、3-D プリントしたマウントによって保持されました。

駆動 DC モーターは 5 V の定電圧で駆動され、4 組のギアで 278:1 の減速を行うため、出力シャフトの回転数は平均値で 22 rpm になりました。ギアの噛み合いとギアの歯の経時的劣化による出力回転数の変動、およびそれに対応するモーター電流の変化は、MATLAB での将来の解析のためセンサーによって取得されました。

フィルター処理と増幅の回路

DC モーター巻線の消費電流は、0.5 Ω 0.5 W の抵抗器両端の電圧低下を測定することにより、オームの法則を使用して計算されました。電流測定値の変化は小さすぎて高精度では検出できないため、単電源センサーインターフェイスアンプ (Analog Devices AD22050N) を使用して電流信号が 20 倍に増幅されました。増幅された電流信号は、アンチエイリアシング 5 次楕円ローパスフィルター (Maxim MAX7408) を使用して平滑化とノイズ除去を行ってから、Arduino Uno のアナログデジタルコンバーター (ADC) ピンに送られました。アンチエイリアシングフィルターのコーナー周波数は、アナログデジタルコンバーターのナイキスト周波数 (750 Hz) で十分に減衰するよう、約 470 Hz に定められました。

電子回路は、下の写真に示される追加の受動コンポーネントをいくつか使用して、プロトタイププリント配線基板 (PCB) 上に作成されました。

以下のフローチャートが示すように、Arduino Uno は ADC 入力を通して送られた電流信号を 1500 Hz のサンプリング周波数でサンプリングし、250000 bps のボーレートで動作するシリアル通信チャネルで、タコメーターのパルスとともにコンピューターにストリーミングしました。

overview flowchart-01.png

この例のデータの場合、コンピューターとのシリアル通信の効率を高めるため、シリアル信号は Arduino UNO 上で、データサンプル (電流 + タコメーター) あたり 3 つの 8 ビットバイナリ値として符号化されました。

リアルタイムのデータストリーミング

Arduino UNO から MATLAB へのデータのストリーミング

Arduino UNO は、効率化のためにモーター電流とタコメーターパルス値をバイナリバッファーに符号化し、5 分ごとに 16384 (2^14) サンプルからなるバッチデータをコンピューターに送信しました。データは、MATLAB コードでのバッチデータの抽出が簡単になるよう、データバイト前後の開始マーカーと終了マーカーで構成される単純なデータプロトコルを使用して、1500 Hz のレートで Arduino UNO から送信されました。

詳細と、ハードウェアセットアップから MATLAB へのデータのストリーミングについては、添付の Arduino スケッチファイル servo_data_producer.ino を参照してください。

MATLAB でのライブデータストリームの受信

MATLAB スクリプトを使用して、適切なボーレートとシリアル COM ポート番号で Arduino Uno からシリアルデータのバッチを取得することができます。

まず、Arduino または同等のハードウェアとの通信用にシリアルポートを開きます。serialportlist コマンドを使用して、シリアル COM ポートを見つけることができます。この例では、Arduino UNO は COM4 に接続されています。250000 bps のボーレートを使用します。

port = "COM4";
baudRate = 250000;

適切なタイムアウトを設定し、ファイル servo_data_producer.ino にある Arduino コードからの正しいデータ終端子を使用します。この例では、6 分のタイムアウトを使用します。次に、flush を使用してシリアルポートのバッファーをクリアします。

s = serialport(port,baudRate,'Timeout',360); 
configureTerminator(s,"CR/LF");
flush(s);

ここで、適切な数の読み取りサイクルを使って、符号化されたデータストリームを Arduino から読み取り、処理します。この例では、10 回の読み取りサイクルを使用します。詳細については、「サポート関数」 節にあるサポート関数 readServoDataFromArduino を参照してください。

count = 0;
countMax = 10;
while (count < countMax)
  try
    count = count + 1;
    fprintf("\nWaiting for dataset #%d...\n", count);
    readServoDataFromArduino(s);
  catch e
    % Report any errors.
    clear s;
    throw(e);
  end
end

ここで、clear コマンドを使用してシリアルポートを閉じます。

clear s;

MATLAB でのデータ処理および特徴抽出

ライブストリームデータの処理

データが MATLAB にストリーミングされると、readServoDataFromArduino スクリプトを使用してデータ行列を構成し、サーボモーターデータの最新のバッチを格納することができます。詳細については、「サポート関数」 節にあるサポート関数 readServoDataFromArduino を参照してください。

readServoDataFromArduino(s)

結果の timetable には、Arduino UNO からストリーミングされた最新データセットのタイムスタンプ、モーター電流、およびタコメーターのパルスの値が含まれています。

関数 processServoData を使用して、適切な物理単位が伴った timetable を出力します。

T = processServoData(X)

結果の timetable は、モーターのデータから予測的特徴を抽出する際の便利なコンテナーとなります。

必要なハードウェアがない場合は、関数 sendSyntheticFeaturesToThingSpeak を使用してこの例を実行し、合成データを生成することができます。

特徴抽出

新しいデータセットそれぞれから、定格回転数を計算してギア列での対象とする周波数を検出し、パワースペクトル上の周波数と正しく照合します。

F = extractFeaturesFromData(T);

1500 Hz のサンプリング周波数値をtachorpmコマンドで使用して、出力シャフトの定格回転数を計算します。このインスタンスでは、RPM 値は 22 rpm 付近で一定です。

tP = T.TachoPulse;
rpm = tachorpm(tP,Fs,'PulsesPerRev',16,'FitType','linear');
rpm = mean(rpm);

モーター回転数をギアセットの物理パラメーターとともに用いると、故障周波数帯域の構成が可能になりますが、これはスペクトルメトリクスを計算するための重要な前提条件です。まず、ギア列の駆動ギアの歯数と定格 rpm を使用して、対象周波数を計算します。対象周波数は実際の出力回転数の値 (Hz 単位) で、次の表に示す理論値と近い値になりました。

次にfaultBandsコマンドを使用して、次の対象周波数を含む、すべての出力回転数に対する周波数帯域を構成します。選択された基本周波数、高調波、および側波帯は、ライブエディタータスクを使用したギア列データの解析およびスペクトル特徴の抽出の例でより詳しく説明されています。特に、サーボシステムで監視すべき重要な故障周波数は、以下のように特定されました。

FS1 の最初の 6 つの高調波と、FS2 の 0:1 側波帯
FS2 の 1 番目、2 番目、4 番目の高調波と、FS3 の 0:1 側波帯
FS3 の 3 番目の高調波

% Number of gear and pinion teeth.
G4 = 41; G3 = 35; G2 = 50; G1 = 62;
P4 = 16; P3 = 10; P2 = 10; P1 = 10;

% Shaft speeds in Hz.
FS5 = rpm / 60;
FS4 = G4 / P4 * FS5;
FS3 = G3 / P3 * FS4;
FS2 = G2 / P2 * FS3;
FS1 = G1 / P1 * FS2;

% Generate the fault bands of interest.
FB_S1 = faultBands(FS1, 1:6, FS2, 0:1);
FB_S2 = faultBands(FS2, [1 2 4], FS3, 0:1);
FB_S3 = faultBands(FS3, 3);
FB = [FB_S1; FB_S2; FB_S3];

faultBandMetricsコマンドをモーター電流信号のパワースペクトル密度 (pwelchコマンド) とともに使用して、合計で 85 のスペクトルメトリクスを計算します。

モーターの電流データのパワースペクトルを計算します。

mC = T.MotorCurrent;
[ps,f] = pwelch(mC,[],[],length(mC),Fs);

すべての故障帯域についてスペクトルメトリクスを計算します。

metrics = faultBandMetrics(ps,f,FB);
metrics = metrics{:, [4 13 85]}; % Select significant metrics.

計算したメトリクスから、追跡する特徴を選択します。シャフト 1 の回転数 (FS1) の 1 番目の高調波と 2 番目の高調波におけるスペクトルのピーク振幅のような、重要な特徴をいくつかを選択します。

features = [metrics,mean(mC),rpm];

モーターの平均電流と平均回転数 (rpm 単位) も、追加の特徴として記録されることに注意してください。

ライブデータストリーミングおよび RUL 推定

クラウドでの特徴データの保存

5 分ごとのバッチでリアルタイムのデータからさまざまな特徴が計算されたら、同じ MATLAB スクリプトの一部として、クラウドストレージと残存耐用期間の解析用に ThingSpeak へと送信します。

関数 sendFeaturesToThingSpeak を使用して、特徴値を別々のフィールドとして ThingSpeak チャネルに送信します。

sendFeaturesToThingSpeak(features)

以下に示すのは、特徴値を保存するいくつかのフィールドで構成されている ThingSpeak チャネルの例です。チャネルは、便利なデータリポジトリおよび可視化プラットフォームとして機能します。

ThingSpeak チャネルは、特徴値の変化を追跡するための、リアルタイムの特徴監視ダッシュボードを提供します。また、残存耐用期間の推定などの特徴処理に用いるデータソースとしても機能します。

サーボモーターから特徴を生成するためのハードウェアセットアップがない場合は、次のコード例を使って、以降の RUL 推定の計算を記述するための合成の特徴値を生成できます。異なったマシンで実行されている別々の MATLAB インスタンスを設けて、ThingSpeak からのデータ読み取りと、そこへの書き込みを行うことができます。

hasArduino = false; % Set to true when using Arduino UNO.
if ~hasArduino
  % Set timer to send new feature data every minute.
  tmr1 = timer('ExecutionMode','fixedSpacing','Period',60);
  tmr1.TimerFcn = @(~,~) sendSyntheticFeaturesToThingSpeak();
  tmr1.start();
  return;
end

クラウドからの特徴データの読み取り

コンパイル済みの MATLAB アプリとして、または展開済みの Web アプリとしても展開できる個別の MATLAB スクリプトを使用して、ThingSpeak からサーボモーターの特徴を読み取り、RUL メトリクスのリアルタイム推定を計算することができます。このインスタンスでは、初期の指数劣化モデルを構成し、新たな特徴値が利用可能になるにつれ定期的にモデルを更新します。

% Set timer to check for new data every minute and update the RUL model, mdl.
tmr2 = timer('ExecutionMode','fixedSpacing','Period',60);
tmr2.TimerFcn = @(~,~) readFeaturesFromThingSpeak(mdl,channelID,readKey);
tmr2.start();

新しい特徴セットはそれぞれ、指数劣化モデルの更新と、RUL のリアルタイム推定の計算に使用されます。

残存耐用期間のライブ推定

この例では、学習データが履歴データではなく、前の手順で計算された特徴で取得したような、コンポーネントの状態のリアルタイムの観測値であると仮定します。指数劣化モデル (exponentialDegradationModel) は、こうしたタイプのリアルタイム RUL 推定と可視化に適しています。RUL モデルでは、健全性インデックスをしきい値として使用し、システムの RUL を推定します。この例の場合、しきい値は、前の手順で計算されたサーボモーターの特徴 Total Band Power に関するものです。合計帯域パワーは、サーボモーターギア列の重要な周波数領域に関連付けられた故障帯域周波数内における、スペクトルエネルギーの変化を捉えます。

任意の事前分布データと指定されたノイズ分散を使って、指数劣化モデルを作成します。

ThingSpeak から読み取ったデータ table のサブセットの観測データに対し、ライフタイムとデータの変数名を指定します。次に、それぞれの特徴値を使用し、モデルに保存されている現在のライフタイム値によってコンポーネントの RUL を予測します。

channelID = 1313749; % Replace with the channel ID to write data to.
readKey  = 'KYIDUZ1ENDT3TGG0'; % Replace with the read API key of your channel.

% Read most recent degradation data from ThingSpeak.
T = thingSpeakRead(channelID,'ReadKey',readKey,'OutputFormat','timetable');

% Construct the exponential degradation model using Total Band Power as the
% health index.
mdl = exponentialDegradationModel( ...
  'Theta', 1, 'ThetaVariance', 100, ...
  'Beta', 1, 'BetaVariance', 100, ...
  'NoiseVariance', 0.01, ...
  'LifeTimeVariable', T.Properties.DimensionNames{1}, ...
  'DataVariables', T.Properties.VariableNames{3}, ...
  'LifeTimeUnit', "hours");

% Set timer to check for new data every minute.
tmr2 = timer('ExecutionMode','fixedSpacing','Period',60);
tmr2.TimerFcn = @(~,~) readFeaturesFromThingSpeak(mdl,channelID,readKey);
tmr2.start();

MATLAB スクリプトは ThingSpeak チャネルで使用できる特徴データを追跡し、新たな特徴値が利用可能になると RUL モデルを更新します。これにより、リアルタイムの RUL 推定と可視化の機能が提供されます。

サポート関数

サポート関数 readServoDataFromArduino は、Arduino UNO から符号化されたデータストリームを読み取り、処理します。

function readServoDataFromArduino(s)
% Reads and processes encoded data streams from Arduino.

% Find start of encoded data stream.
str = readline(s);
N = sscanf(str, "BEGIN:%d"); % Expected number of data points (should be 16384).

% Acquire new data until the end of data stream.
if ~isempty(N) && N > 0
  X = zeros(N,2);

  % Read and store encoded binary data.
  disp("Reading servo motor data...");
  for k = 1:N
    % Read three bytes each time containing motor current and tacho pulse values.
    x = [read(s,1,"uint16"), read(s,1,"uint8")];
    X(k,:) = double(x);
  end

  % Find end of encoded data stream.
  readline(s); % Remove CR/LF from data stream.
  str = readline(s);
  if str == "END"
    disp("Processing servo data...");
    T = processServoData(X);

    disp('Extracting features...');
    F = extractFeaturesFromData(T);

    disp('Sending features to ThingSpeak...');
    sendFeaturesToThingSpeak(F);
  else
    error('Error reading end-of-file marker from serial port.');
  end
end
end

サポート関数 processServoData は、生のデータ行列を、適切な変数名と物理単位を伴う timetable に変換します。

function T = processServoData(X)
% Converts raw data matrix to a timetable with appropriate variable names and
% physical units.
Fs = 1500; % Sampling rate.
Rsens = 0.5; % Current sense resistor value in ohm.
Kamp = 20; % Sensor amplifier gain.
count2volt = 5/1024; % Scaling factor between digital reading to analog voltage.

T = timetable('SampleRate', Fs);
T.MotorCurrent = X(:,1) * count2volt / Kamp / Rsens * 1000; % Convert to mA.
T.TachoPulse   = X(:,2); % On/off tacho pulse values.
T.Properties.VariableUnits = {'mA', 'level'};
end

サポート関数 extractFeaturesFromData は、時間領域のサーボモーターデータから特徴を table 形式で抽出します。

function features = extractFeaturesFromData(T)
% Extracts up to 8 features from time-domain servo data provided in table
% format.
Fs = 1500; % Sampling rate.

% Average motor speed.
tP = T.TachoPulse;
rpm = tachorpm(tP, Fs, 'PulsesPerRev', 16, 'FitType', 'linear');
rpm = mean(rpm);

% Number of gear and pinion teeth.
G4 = 41; G3 = 35; G2 = 50; G1 = 62;
P4 = 16; P3 = 10; P2 = 10; P1 = 10;

% Shaft speeds in Hz.
FS5 = rpm / 60;
FS4 = G4 / P4 * FS5;
FS3 = G3 / P3 * FS4;
FS2 = G2 / P2 * FS3;
FS1 = G1 / P1 * FS2;

% Generate the fault bands of interest.
FB_S1 = faultBands(FS1, 1:6, FS2, 0:1);
FB_S2 = faultBands(FS2, [1 2 4], FS3, 0:1);
FB_S3 = faultBands(FS3, 3);
FB = [FB_S1; FB_S2; FB_S3];

% Compute the power spectrum of the motor current data.
mC = T.MotorCurrent;
[ps,f] = pwelch(mC, [], [], length(mC), Fs);

% Compute spectral metrics for all fault bands.
metrics = faultBandMetrics(ps, f, FB);
metrics = metrics{:, [4 13 85]}; % Select significant metrics.

% Pick features to track.
features = [metrics, mean(mC), rpm];
end

サポート関数 sendFeaturesToThingSpeak は、特徴値を別々のフィールドとして ThingSpeak チャネルに送信します。

function sendFeaturesToThingSpeak(features)
% Sends feature values to a ThingSpeak channel as separate fields.

% Configure the Channel ID and the Write API Key values.
channelID = 1313749; % Replace with your channel ID to write data to.
writeKey = 'X96MRW1TTZC1XNV3'; % Replace with the write API key of your channel.

% Write the features to the channel specified by the 'channelID' variable.
fields = 1:numel(features); % Up to 8 features
thingSpeakWrite(channelID, features, 'Fields', fields, 'WriteKey', writeKey);
end

サポート関数 readFeaturesFromThingSpeak は、最新のデータの読み取りのタイムスタンプを追跡することで、新たなデータが利用可能になると ThingSpeak チャネルから特徴値を読み取ります。この関数はまた、新しい特徴値を使用して RUL モデルを更新し、予測された RUL 値をプロットします。

function readFeaturesFromThingSpeak(mdl, channelID, readKey)
% Read feature values from ThingSpeak to update the RUL estimation.
persistent timestamp
persistent hplot

if isempty(timestamp)
  % Start with last ten days of data.
  timestamp = dateshift(datetime, "start", "day", -10);
end

% Read recent data from all fields.
disp('Reading features from ThingSpeak...');
T = thingSpeakRead(channelID, 'ReadKey', readKey, ...
  'OutputFormat', 'timetable', 'DateRange', [timestamp, datetime]);
if ~isempty(T)
  T = T(T.Timestamps > timestamp, :); % Only keep most recent data.
end

if ~isempty(T)
  T = T(T.Timestamps > timestamp, :); % Only keep most recent data.
  timestamp = T.Timestamps(end); % Update time stamp for most recent data.

  % Set RUL degradation threshold.
  threshold = 10;

  % The training data is not historical data, but rather real-time observations
  % of the servo motor features.
  N = height(T);
  estRUL = hours(zeros(1,N));
  for i = 1:N
    update(mdl, T(i,:))
    estRUL(i) = predictRUL(mdl, threshold);
  end

  % Visualize RUL over time.
  if isempty(hplot) || ~isvalid(hplot)
    hplot = plot(T.Timestamps(:)', estRUL, 'r-o');
    title('Estimated RUL at Time t')
    xlabel('Time, t')
    ylabel('Estimated RUL')
  else
    hplot.XData = [hplot.XData, T.Timestamps(:)'];
    hplot.YData = [hplot.YData, estRUL];
  end
end
end

特徴値の生成にサーボモーターと Arduino UNO が利用できない場合は、サポート関数 sendSyntheticFeaturesToThingSpeak を使用して合成の特徴データを用意し、ThingSpeak への特徴値のストリーミングの記述に用いてください。

function sendSyntheticFeaturesToThingSpeak()
% Generate synthetic features and send them to ThinkSpeak.
persistent bandPower
if isempty(bandPower)
  bandPower = 5;
end

% Simulate motor features and fault band metrics.
motorCurrent = 308+12*rand(1); % between [308 320] mA
rpm = 20+3*rand(1); % between [20 23] rpm
bandPower = bandPower + rand(1)/2; % bandPower degradation
metrics = [rand(1), rand(1), bandPower];
features = [metrics, motorCurrent, rpm];

disp('Sending features to ThingSpeak...');
sendFeaturesToThingSpeak(features)
end

参考

faultBands | faultBandMetrics | exponentialDegradationModel | pwelch | tachorpm

ThingSpeak を使用したサーボのギア列のライブ RUL 推定

ハードウェア セットアップとデータ収集