メインコンテンツ

Simulink における LSTM ネットワークを使用した物理システムのモデリング

この例では、長短期記憶 (LSTM) ニューラル ネットワークを使用して、Simulink® モデルでバーチャル センサーとして機能する低次元化モデル (ROM) を作成する方法を示します。

物理モデルの代わりに ROM を使用すると、元の物理モデルの精度を損なうことなく必要な計算量を減らすことができます。組み込みシステムでモデルを高サンプル レートで実行する必要があるバーチャル センサーなどの用途において、ROM は元の物理システムを適切に近似できます。ROM を構築するにはさまざまな手法がありますが、この例では、LSTM-ROM (LSTM ネットワークを活用する ROM の一種) を構築し、それを Deep Learning Stateful Predict ブロックの一部として Simulink モデルで使用します。

この例では、LSTM ネットワークに学習させるために、元のモデルを使用して学習データを生成します。ROM 内の学習済み LSTM は、負荷シャフトの B 信号と F 信号、ならびに制御圧力を入力として受け取り、負荷シャフトの B 信号と F 信号の次の値を予測します。この例では、モデルの学習後に LSTM-ROM コンポーネントを作成し、Simulink モデル内の物理コンポーネントと置き換えます。次の図は、Simulink モデル内の物理コンポーネントを LSTM-ROM サブコンポーネントに置き換える様子を示しています。

LSTM-ROM 用の LSTM ネットワークの学習は大量の計算を必要とするタスクであり、実行に長い時間がかかることがあります。この例を高速化するため、この例では学習をスキップして事前学習済みのネットワークを読み込みます。代わりにネットワークに学習させるには、doTraining 変数を true に設定します。

doTraining = false;

学習データの生成

この例にサポート ファイルとして添付されている Simulink モデル ex_sdl_flexible_shaft を使用して、学習データを生成します。このモデルにアクセスするには、この例をライブ スクリプトとして開きます。この Simulink モデルはフレキシブル シャフトをシミュレートします。このモデルには、集中定数法を使用してモデル化された 2 つのアルミニウム製フレキシブル シャフト (20 個のセグメントから成るモーター シャフトと 5 個のセグメントから成る負荷シャフト) があります。どちらのシャフトにも、慣性、減衰、および剛性ねじりバネが含まれています。シミュレーションの開始時、クラッチのロックは解除されており、被駆動シャフトは自由に動きます。モーター シャフトの初期速度は 200 rad/s に設定されており、システムは定常状態で起動します。このモデルは、クラッチに与えられる圧力を制御パラメーターとして使用し、モデルのダイナミクスを決定します。

この Simulink モデルは次の 4 つの値を出力します。

  1. 負荷シャフトのベース (B) 信号

  2. 負荷シャフトの従属 (F) 信号

  3. 制御圧力

  4. モーター シャフトの F 信号

最初の 3 つの状態は、LSTM-ROM の学習に使用されます。4 番目の状態は、学習済みのモデルの精度を評価するために使用されます。

次の図は、このモデルの構造を示しています。

この Simulink モデルは、ワークスペース変数 stopTimetimeInterval に依存します。これらの変数は、それぞれシミュレーションの最終タイム ステップと出力タイム ステップ間の間隔を指定します。シミュレーションの初期タイム ステップは 0 です。

停止時間を 0.2 に指定し、時間間隔を 5×10-5 に指定します。

stopTime = 0.2;
timeInterval = 5e-5;

このモデルの Pressure ブロックは、クラッチに与えられる最大圧力の値を定義するワークスペース変数 maxPressure に依存します。105106 の間で等間隔に設定された 20 個の異なる maxPressure の値について、モデルを実行します。cell 配列 data に出力データを収集します。各要素は、指定された圧力プロファイルを使用して計算された時系列観測値に対応します。

numObservations = 20;
maxPressures = linspace(1e5,1e6,numObservations);

data = cell(numObservations,1);

for i = 1:numObservations
    maxPressure = maxPressures(i);
    simout = sim("ex_sdl_flexible_shaft");
    data{i} = simout.simout.Data';
end

シミュレーションの時間ステップを抽出します。

times = simout.simout.Time;
numTimeSteps = length(times);

最初の 5 つのシミュレーションの制御圧力をプロットします。

figure
for i = 1:5
    pressure = data{i}(3,:);
    plot(times,pressure);
    hold on
end

title("Input Pressure")
legend("Observation " + (1:5))
xlabel("Time")
ylabel("Pressure (Pa)")
hold off

シミュレーションの 1 つにおける負荷シャフトの B 信号と F 信号、ならびにモーター シャフトの F 信号をプロットします。

idx = 4;
BLoadShaft = data{idx}(1,:);
FLoadShaft = data{idx}(2,:);
FMotorShaft = data{idx}(4,:);

figure
plot(times,BLoadShaft, ...
    times,FLoadShaft, ...
    times,FMotorShaft)

legend("B - Load Shaft", "F - Load Shaft", "F - Motor Shaft")
title("Model Dynamics (Maximum Pressure =  " + maxPressures(idx) + " Pa)")

学習用データの準備

非常に長いシーケンスで LSTM に学習させる場合、各タイム ステップで計算された勾配の累積によって勾配が消失し、学習プロセスが準最適の結果に収束する可能性があります。勾配の消失を防ぐには、学習データをダウンサンプリングし、情報の大幅な喪失を抑えながらシーケンスを大幅に短くします。

データをダウンサンプリングするには、学習データのサンプル時間よりも大きいサンプル時間を指定します。Stateful Predict ブロックで LSTM ネットワークを使用する場合は、[サンプル時間] パラメーターに同じ値を指定しなければなりません。

サンプル時間を 10-3 に指定します。

sampleTime = 1e-3;

サンプル時間をシミュレーションの時間間隔で割った固定間隔でタイム ステップを抽出し、学習データをダウンサンプリングします。

intervalDownsampled = sampleTime / timeInterval;
timeStepsDownsampled = 1:intervalDownsampled:numTimeSteps;

for i = 1:numObservations
    dataDownsampled{i} = data{i}(:,timeStepsDownsampled);
end

この例にサポート ファイルとして添付されている trainingPartitions 関数を使用して、学習データを学習用とテスト用の区画に均等に分割します。このファイルにアクセスするには、例をライブ スクリプトとして開きます。

[idxTrain,idxTest] = trainingPartitions(numObservations,[0.5 0.5]);

maxPressuresTrain = maxPressures(idxTrain);
maxPressuresTest = maxPressures(idxTest);
dataTrain = dataDownsampled(idxTrain);
dataTest = dataDownsampled(idxTest);

学習データから予測子とターゲットを抽出します。予測子は、負荷シャフトの B 信号と F 信号、ならびに制御圧力です。ターゲットは、1 タイム ステップ分シフトされた負荷シャフトの B 信号と F 信号です。予測子は、dataTrain の最初の 3 つのチャネルに対応します。ターゲットは、1 タイム ステップ分シフトされた dataTrain の各要素における最初の 2 つのチャネルに対応します。

inputStatesTrain  = [1 2 3];
outputStatesTrain = [1 2];

numObservationsTrain = numel(dataTrain);
for i = 1:numObservationsTrain
    XTrain{i} = dataTrain{i}(inputStatesTrain,1:end-1);
    TTrain{i} = dataTrain{i}(outputStatesTrain,2:end);
end

ネットワーク アーキテクチャの定義

次の B 信号と F 信号の値を予測する LSTM ネットワークを以下のように定義します。

  • シーケンス入力の場合は、入力サイズが入力数に一致するシーケンス入力層を指定します。入力を再スケーリングし、値が 0 から 1 の間になるように正規化します。

  • 入力特徴間の交互作用を学習させるには、出力サイズが 200 の全結合層と、それに続く ReLU 層を含めます。

  • シーケンス データの長期的な依存関係を学習させるには、200 個の隠れユニットをもつ 2 つの LSTM 層と、それに続く ReLU 層を含めます。

  • 正しいサイズの予測を出力するには、応答の数に一致するサイズの全結合層を含めます。

  • 回帰用にネットワークに学習させるには、回帰層を含めます。

numHiddenUnits = 200;

numFeatures = numel(inputStatesTrain);
numResponses = numel(outputStatesTrain);

layers = [
    sequenceInputLayer(numFeatures,Normalization="rescale-zero-one")
    fullyConnectedLayer(numHiddenUnits)
    reluLayer
    lstmLayer(numHiddenUnits)
    lstmLayer(numHiddenUnits)
    reluLayer
    fullyConnectedLayer(numResponses)
    regressionLayer];

学習オプションの指定

学習オプションを指定します。

  • Adam ソルバーを使用して、ネットワークの学習を 3×104 エポック行います。

  • 勾配が発散するのを防ぐには、しきい値 1 で勾配をクリップします。

  • 学習を改善するには、学習率係数が段階的に減少するようにスケジュールします。初期学習率 5×10-3 を使用し、104 回の反復ごとに学習率を 40% 減らします。

  • プロットに学習の進行状況を表示し、詳細出力を表示しないようにします。

  • GPU が利用できる場合、GPU で学習を行います。既定では、trainNetwork 関数は利用可能な GPU で学習を行います。GPU を使用するには、Parallel Computing Toolbox とサポートされている GPU デバイスが必要です。サポートされているデバイスの詳細については、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。

options = trainingOptions("adam", ...
    MaxEpochs=3e4, ...
    GradientThreshold=1, ...
    InitialLearnRate=5e-3, ...
    LearnRateSchedule="piecewise", ...
    LearnRateDropPeriod=1e4, ...
    LearnRateDropFactor=0.6, ...
    Verbose=0, ...
    Plots="training-progress");

LSTM ネットワークの学習

trainNetwork 関数を使用して LSTM ネットワークに学習させます。Simulink モデルで LSTM ネットワークを使用するには、ネットワークを MAT ファイルに保存します。

LSTM-ROM 用の LSTM ネットワークの学習は大量の計算を必要とするタスクであり、実行に長い時間がかかることがあります。この例を高速化するため、この例では学習をスキップして事前学習済みのネットワークを読み込みます。代わりにネットワークに学習させるには、doTraining 変数を true に設定します。

filename = "flexibleShaftLoadNet.mat";
if doTraining
    net = trainNetwork(XTrain, TTrain, layers, options);
    save(filename,"net")
else
    load(filename);
end

Simulink モデルにおける LSTM-ROM の使用

学習済みの LSTM ネットワークを使用して、負荷シャフトの予測された B 信号と F 信号を出力する Simulink サブコンポーネントを以下のように作成します。このブロックは、フィードバック ループを通して次のタイム ステップの予測を使用します。

Stateful Predict ブロックで、[ファイル パス] パラメーターを保存された学習済みネットワークの場所に設定し、学習データの生成に使用したのと同じサンプル時間を指定します。Stateful Predict ブロックの入力に接続された Rate Transition ブロックと、Unit Delay ブロックで、同じサンプル時間を使用します。

元の Simulink モデルの負荷シャフトを、LSTM-ROM サブコンポーネントに置き換えます。結果のモデルは ex_sdl_flexible_shaft_lstm に保存されます。このファイルにアクセスするには、例をライブ スクリプトとして開きます。

モデルのテスト

ホールドアウトされたテスト データ セットを使用して、モデルの精度を評価します。

LSTM ネットワークだけでなく、モデル全体の精度をテストするには、元の Simulink モデルによって生成されたモーター シャフトの出力とシミュレートされた F 信号を比較します。

テスト データからターゲットを抽出します。テスト ターゲットは、元のモデルを使用してシミュレートされたモーター シャフトの F 信号を 1 タイム ステップ分シフトした信号です。

numObservationsTest = numel(dataTest);

outputStatesTest = 4;

for i = 1:numObservationsTest
    TTest{i} = dataTest{i}(outputStatesTest,2:end);
end

テスト データ セットの最大圧力値ごとにシミュレーションを実行し、シミュレートされたモーター シャフトの F 信号を cell 配列 YTest に保存します。

errs = [];
for i = 1:numObservationsTest
    maxPressure = maxPressuresTest(i);
    simout = sim("ex_sdl_flexible_shaft_lstm");

    YTest{i} = simout.simout.Data(1:end-1,4)';
end

散布図で各タイム ステップにおける予測を可視化します。

figure
scatter([TTest{:}], [YTest{:}])
xlabel("Target")
ylabel("Prediction")

m = min([TTest{:} YTest{:}]);
M = max([TTest{:} YTest{:}]);
xlim([m M])
ylim([m M])

hold on
plot([m M],[m M],"r--")

予測誤差をヒストグラムで可視化します。

figure
histogram([TTest{:}] - [YTest{:}])
xlabel("Error")
ylabel("Frequency")

ゼロに近い値は、予測が正確であることを示します。

新しいデータを使用したシミュレーション

未確認の最大圧力値 2.7×105 Pa でモデルを実行します。

maxPressure = 2.7e5;
simout = sim("ex_sdl_flexible_shaft_lstm");
Y = simout.simout.Data(:,4)';

予測された F 信号をプロットで可視化します。

figure
plot(simout.simout.Time, Y)
ylim([0 inf])
xlabel("Time")
ylabel("F Signal - Motor Shaft")
title("Model Predictions (Maximum pressure = " + maxPressure + " Pa)")

参考

| | | | |

トピック