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

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

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

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

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

次の図は、Simulink モデルで置き換える物理コンポーネントを示しています。

次の図は、物理コンポーネントを LSTM-ROM サブコンポーネントで置き換えた Simulink モデルを示しています。

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

doTraining = false;

学習データの生成

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

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

負荷シャフトのベース (B) 信号
負荷シャフトの従属 (F) 信号
制御圧力
モーターシャフトの F 信号

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

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

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

停止時間を 0.2 に指定し、時間間隔を $5 \times 1 0^{- 5}$ に指定します。

stopTime = 0.2;
timeInterval = 5e-5;

このモデルの Pressure ブロックは、クラッチに与えられる最大圧力の値を定義するワークスペース変数 maxPressure に依存します。 $1 0^{5}$ と $1 0^{6}$ の間で等間隔に設定された 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 に学習させる場合、各タイムステップで計算された勾配の累積によって勾配が消失し、学習プロセスが準最適の結果に収束する可能性があります。勾配の消失を防ぐには、学習データをダウンサンプリングし、情報の大幅な喪失を抑えながらシーケンスを大幅に短くします。

データをダウンサンプリングするには、学習データのサンプル時間よりも大きいサンプル時間を指定します。LSTM ネットワークを層ブロックとして Simulink にエクスポートする場合は、exportNetworkToSimulink の名前と値の引数 SampleTime に同じ値を指定しなければなりません。

サンプル時間を $1 0^{- 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}(1:end-1,inputStatesTrain);
    TTrain{i} = dataTrain{i}(2:end,outputStatesTrain);
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)];

学習オプションの指定

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

Adam ソルバーを使用してネットワークに 500 エポック学習させます。
勾配が発散するのを防ぐには、しきい値 1 で勾配をクリップします。
学習を改善するには、学習率係数が段階的に減少するようにスケジュールします。初期学習率 $5 \times 1 0^{- 3}$ を使用し、 $1 0^{4}$ 回の反復ごとに学習率を 40% 減らします。
プロットに学習の進行状況を表示し、詳細出力を表示しないようにします。
GPU が利用できる場合、GPU で学習を行います。既定では、trainnet 関数は利用可能な GPU で学習を行います。GPU を使用するには、Parallel Computing Toolbox とサポートされている GPU デバイスが必要です。サポートされているデバイスの詳細については、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。

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

LSTM ネットワークの学習

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

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

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

exportNetworkToSimulink を使用して、ネットワークを Simulink にエクスポートします。

モデルに "lstm_model" という名前を付ける。
Simulink モデルがネットワークの独自のコピーを確実に保持するように、ネットワークをモデルワークスペースに保存する。
元のモデルと同じサンプル時間を使用する。

exportNetworkToSimulink(net, ...
    ModelName="lstm_model", ...
    SaveNetworkInModelWorkspace=true, ...
    SampleTime=string(sampleTime));

モデルのサブシステムを表示します。サブシステムは、ニューラルネットワークを深層学習層ブロックのモデルとして表現します。

LSTM-ROM コンポーネントを作成します。エクスポートされたモデルをサブシステムとして使用して、負荷シャフトの予測された B 信号と F 信号を出力する Simulink コンポーネントを作成します。モデルはフィードバックループを通じて次のタイムステップの予測を使用します。

元の 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}(2:end, outputStatesTest);
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
x = cat(1,TTest{:});
y = cat(1,YTest{:});
scatter(x,y)
xlabel("Target")
ylabel("Prediction")

m = min([x y],[],"all");
M = max([x y],[],"all");
xlim([m M])
ylim([m M])

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

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

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

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

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

未確認の最大圧力値 $2.7 \times 1 0^{5}$ 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)")

参考