使用 LSTM 网络进行线性系统辨识

此示例使用:

本示例演示了如何使用长短期记忆 (LSTM) 神经网络对线性系统进行估计，并将该方法与传递函数估计方法进行了比较。

例如，在此示例中，您将研究 LSTM 网络捕捉建模系统内在动态特性的能力。为此，您需要利用线性传递函数的输入和输出信号对 LSTM 网络进行训练，并测量该网络对阶跃变化的响应精度。

传递函数

本示例使用了一个包含快慢混合动态特性且阻尼适中的四阶传递函数。适度的阻尼会导致系统动态特性在较长的时间尺度上逐渐衰减，同时也展示了 LSTM 网络在捕获混合动态特性的同时，仍能避免某些重要响应动态特性衰减的能力。通过指定系统的零极点来构建传递函数。

fourthOrderMdl = zpk(-4,[-9+5i;-9-5i;-2+50i;-2-50i],5e5);
[stepResponse,stepTime] = step(fourthOrderMdl);

绘制传递函数的阶跃响应曲线。

plot(stepTime,stepResponse)
grid on
axis tight
title('Fourth-order mixed step response')

Figure contains an axes object. The axes object with title Fourth-order mixed step response contains an object of type line.

波德图显示了系统的带宽，其测量值是增益首次降至 DC 值 70.8% 以下（即约 3 dB）时的频率。

bodeplot(fourthOrderMdl)

MATLAB figure

fb = bandwidth(fourthOrderMdl)

fb = 
62.8858

生成训练数据

构建一个包含输入和输出信号的数据集，用于训练 LSTM 网络。作为输入，生成一个随机高斯信号。仿真传递函数 fourthOrderMdl 对该输入的响应，以获得输出信号。

高斯噪声数据

指定随机高斯噪声训练信号的属性。

rng default
signalType = 'rgs'; % Gaussian 
signalLength = 5000; % Number of points in the signal
fs = 100; % Sampling frequency
signalAmplitude = 1; % Maximum signal amplitude

使用 idinput 生成高斯噪声信号，并对结果进行缩放。

urgs = idinput(signalLength,signalType);
urgs = (signalAmplitude/max(urgs))*urgs';

根据采样率生成时序信号。

trgs = 0:1/fs:length(urgs)/fs-1/fs;

使用 lsim 函数生成系统的响应，并将结果存储在 yrgs 中。对仿真输出进行转置，使其符合 LSTM 数据结构的要求 - 该结构需要行向量而非列向量。

yrgs = lsim(fourthOrderMdl,urgs,trgs);
yrgs = yrgs';

同样，创建一个较短的验证信号，供网络训练时使用。

xval = idinput(100,signalType);
yval = lsim(fourthOrderMdl,xval,trgs(1:100));

创建和训练网络

通过贝叶斯优化算法确定了以下网络架构，其中贝叶斯优化成本函数使用了独立的验证数据（详情请参阅随附的 bayesianOptimizationForLSTM.m）。尽管有多种架构可能有效，但这种优化方案在计算效率方面最为出色。优化过程还表明，当将 LSTM 应用于其他线性函数时，尽管函数的复杂度会增加，但网络架构并不会发生显著变化。相反，训练该网络所需的迭代次数会增加。建模一个系统所需的隐藏单元数量与动态特性振荡衰减所需的时间长短有关。在这种情况下，响应分为两个截然不同的部分：高频响应和低频响应。为了捕捉低频响应，需要更多隐藏单元。即使选择较少的单元数，高频响应仍会被建模。然而，低频响应的估计精度会下降。

构建网络架构。

numResponses = 1;
featureDimension = 1;
numHiddenUnits = 200;
maxEpochs = 1000;
miniBatchSize = 300;

Networklayers = [sequenceInputLayer(featureDimension) ...
    lstmLayer(numHiddenUnits) ...
    lstmLayer(numHiddenUnits) ...
    dropoutLayer(0.02),...
    fullyConnectedLayer(numResponses) ...
    regressionLayer];

初始学习率会影响网络的训练效果。如果初始学习率设定得过高，会导致梯度过大，从而延长训练时间。更长的训练时间可能会导致网络的全连接层出现饱和。当网络饱和时，输出值会发散，网络将输出一个 NaN 值。因此，请使用默认值 0.001，这是一个相对较低的初始学习率。这导致残差曲线和损失曲线大多呈单调递减趋势。使用分段速率表，以防止优化算法在优化过程开始时陷入局部极小值。

options = trainingOptions('adam', ...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'GradientThreshold',20, ...
    'Shuffle','once', ...
    'Plots','training-progress',...
    'ExecutionEnvironment','parallel',...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropPeriod',200,...
    'L2Regularization',1e-3,...
    'LearnRateDropFactor',0.5,...
    'Verbose',0,...
    'ValidationData',[{xval'} {yval'}]);

loadNetwork = true; % Set to false to train the network using a parallel pool.
if loadNetwork
    load('fourthOrderMdlnet','fourthOrderNet')
else
    rng('default')
    fourthOrderNet = trainNetwork(urgs,yrgs,Networklayers,options);
    save('fourthOrderMdlnet','fourthOrderNet','urgs','yrgs');
end

评估模型性能

如果一个网络能够成功捕捉系统的动态行为，那么该网络的表现就是良好的。通过测量网络准确预测系统对阶跃输入响应的能力，来评估网络性能。

构建一个阶跃输入。

stepTime = 2; % In seconds
stepAmplitude = 0.1;
stepDuration = 4; % In seconds

% Construct step signal and system output.
time = (0:1/fs:stepDuration)';
stepSignal = [zeros(sum(time<=stepTime),1);stepAmplitude*ones(sum(time>stepTime),1)];
systemResponse = lsim(fourthOrderMdl,stepSignal,time);

% Transpose input and output signal for network inputs.
stepSignal = stepSignal';
systemResponse = systemResponse';

使用已训练好的神经网络来评估系统响应。在图中比较该系统与估计响应。

fourthOrderMixedNetStep = predict(fourthOrderNet,stepSignal);

figure
title('Step response estimation')
plot(time,systemResponse,'k', time, fourthOrderMixedNetStep)
grid on
legend('System','Estimated')
title('Fourth-Order Step')

Figure contains an axes object. The axes object with title Fourth-Order Step contains 2 objects of type line. These objects represent System, Estimated.

该图显示了拟合中存在的两个问题。首先，网络的初始状态并非稳态，这导致信号起始处出现瞬态行为。其次，该网络的预测结果存在轻微偏差。

初始化网络并调整拟合

为了将网络状态初始化为正确的初始条件，必须更新网络状态，使其与测试信号开始时的系统状态相一致。

可以通过将系统在初始条件下的估计响应与系统的实际响应进行比较，来调整网络的初始状态。利用网络对初始状态的估计值与初始状态的实际响应之间的差异，来校正系统估计中的偏移量。

设置网络初始状态

当网络使用从 0 到 1 的阶跃输入进行估计时，LSTM 网络的状态（LSTM 各层的单元格状态和隐藏状态）会向正确的初始条件收敛。为了直观地展示这一点，请使用 predictAndUpdateState 函数提取网络在每个时间步的单元格状态和隐藏状态。

仅使用该步骤（发生在 2 秒时）之前的单元格和隐藏状态值。定义一个 2 秒的时间标记，并提取截至该标记的值。

stepMarker = time <= 2;
yhat = zeros(sum(stepMarker),1);
hiddenState = zeros(sum(stepMarker),200); % 200 LSTM units
cellState = zeros(sum(stepMarker),200);
for ntime = 1:sum(stepMarker)
    [fourthOrderNet,yhat(ntime)] = predictAndUpdateState(fourthOrderNet,stepSignal(ntime)');
    hiddenState(ntime,:) = fourthOrderNet.Layers(2,1).HiddenState;
    cellState(ntime,:) = fourthOrderNet.Layers(2,1).CellState;
end

接下来，绘制该步之前时段内的隐藏状态和单元状态，并确认它们已收敛到固定值。

figure
subplot(2,1,1)
plot(time(1:200),hiddenState(1:200,:))
grid on
axis tight
title('Hidden State')
subplot(2,1,2)
plot(time(1:200),cellState(1:200,:))
grid on
axis tight
title('Cell State')

Figure contains 2 axes objects. Axes object 1 with title Hidden State contains 200 objects of type line. Axes object 2 with title Cell State contains 200 objects of type line.

要初始化零输入信号的网络状态，请选择零输入信号，并设定适当的持续时间，以确保信号足够长，使网络能够达到稳态。

initializationSignalDuration = 10; % In seconds
initializationValue = 0;
initializationSignal = initializationValue*ones(1,initializationSignalDuration*fs);

fourthOrderNet = predictAndUpdateState(fourthOrderNet,initializationSignal);

验证初始条件是否为零或接近零。

figure
zeroMapping = predict(fourthOrderNet,initializationSignal);
plot(zeroMapping)
axis tight

Figure contains an axes object. The axes object contains an object of type line.

既然网络已正确初始化，请使用该网络再次预测阶跃响应，并绘制结果。最初的骚动已经平息。

fourthOrderMixedNetStep = predict(fourthOrderNet,stepSignal);

figure
title('Step response estimation')
plot(time,systemResponse,'k', ...
    time,fourthOrderMixedNetStep,'b')
grid on
legend('System','Estimated')
title('Fourth-Order Step - Adjusted State')

Figure contains an axes object. The axes object with title Fourth-Order Step - Adjusted State contains 2 objects of type line. These objects represent System, Estimated.

调整网络偏移量

即使将网络初始状态设置为补偿测试信号的初始条件，预测响应中仍可见微小的偏移。这是因为 LSTM 网络在训练过程中学习到了不正确的偏置项。您可以使用与更新网络状态时相同的初始化信号来修正偏移量。初始化信号应将网络映射为零。零点与网络估计值之间的偏移量，即为网络所学习的偏置项中的误差。将各层计算出的偏置项相加，其结果与响应中检测到的偏置值非常接近。不过，调整网络输出端的网络偏置项，比调整网络每一层中的单独偏置项要容易一些。

bias = mean(predict(fourthOrderNet,initializationSignal));
fourthOrderMixedNetStep = fourthOrderMixedNetStep-bias;

figure
title('Step response estimation')
plot(time,systemResponse,'k',time,fourthOrderMixedNetStep,'b-')
legend('System','Estimated')
title('Fourth-Order Step - Adjusted Offset')

Figure contains an axes object. The axes object with title Fourth-Order Step - Adjusted Offset contains 2 objects of type line. These objects represent System, Estimated.

移出训练区

用于训练该网络的所有信号的最大振幅均为 1，而阶跃函数的振幅为 0.1。现在，请研究网络在这些范围之外的行为。

时间偏移

通过调整步骤的时间来引入时间偏移。将该步骤的时间设置为 3 秒，比训练集中的时间长 1 秒。将生成的网络输出绘制成图，并注意输出已正确延迟了 1 秒。

stepTime = 3; % In seconds
stepAmplitude = 0.1;
stepDuration = 5; % In seconds
[stepSignal,systemResponse,time] = generateStepResponse(fourthOrderMdl,stepTime,stepAmplitude,stepDuration);

fourthOrderMixedNetStep = predict(fourthOrderNet,stepSignal);
bias = fourthOrderMixedNetStep(1) - initializationValue;
fourthOrderMixedNetStep = fourthOrderMixedNetStep-bias;

figure
plot(time,systemResponse,'k', time,fourthOrderMixedNetStep,'b')
grid on
axis tight

Figure contains an axes object. The axes object contains 2 objects of type line.

振幅偏移

接下来，增加阶跃函数的振幅，以研究当系统输入超出训练数据范围时网络的行为。要测量训练数据范围之外的漂移，可以测量高斯噪声信号中振幅的概率密度函数。将振幅以直方图的形式可视化。

figure
histogram(urgs,'Normalization','pdf')
grid on

Figure contains an axes object. The axes object contains an object of type histogram.

根据分布的百分位数设置阶跃函数的振幅。将误差率作为百分位数的函数绘制出来。

pValues = [60:2:98, 90:1:98, 99:0.1:99.9 99.99];
stepAmps = prctile(urgs,pValues); % Amplitudes
stepTime = 3; % In seconds
stepDuration = 5; % In seconds

stepMSE = zeros(length(stepAmps),1);
fourthOrderMixedNetStep = cell(length(stepAmps),1);
steps = cell(length(stepAmps),1);

for nAmps = 1:length(stepAmps)
    % Fourth-order mixed
    [stepSignal,systemResponse,time] = generateStepResponse(fourthOrderMdl,stepTime,stepAmps(nAmps),stepDuration);
    
    fourthOrderMixedNetStep{nAmps} = predict(fourthOrderNet,stepSignal);
    bias = fourthOrderMixedNetStep{nAmps}(1) - initializationValue;
    fourthOrderMixedNetStep{nAmps} = fourthOrderMixedNetStep{nAmps}-bias;
    
    stepMSE(nAmps) = sqrt(sum((systemResponse-fourthOrderMixedNetStep{nAmps}).^2));
    steps{nAmps,1} = systemResponse;
end

figure
plot(pValues,stepMSE,'bo')
title('Prediction Error as a Function of Deviation from Training Rrange')
grid on
axis tight

Figure contains an axes object. The axes object with title Prediction Error as a Function of Deviation from Training Rrange contains a line object which displays its values using only markers.

subplot(2,1,1)
plot(time,steps{1},'k', time,fourthOrderMixedNetStep{1},'b')
grid on
axis tight
title('Best Performance')
xlabel('time')
ylabel('System Response')
subplot(2,1,2)
plot(time,steps{end},'k', time,fourthOrderMixedNetStep{end},'b')
grid on
axis tight
title('Worst Performance')
xlabel('time')
ylabel('System Response')

Figure contains 2 axes objects. Axes object 1 with title Best Performance, xlabel time, ylabel System Response contains 2 objects of type line. Axes object 2 with title Worst Performance, xlabel time, ylabel System Response contains 2 objects of type line.

当阶跃响应的振幅超出训练集的范围时，LSTM 会尝试估计该响应的平均值。

这些结果表明，使用与用于预测的数据处于同一量级的训练数据非常重要。否则，预测结果将不可靠。

更改系统带宽

通过使用四种不同的网络对四阶混合动态传递函数进行建模，研究系统带宽对 LSTM 网络所选隐藏单元数量的影响：

一个包含 5 个隐藏单元和单层 LSTM 的小型神经网络
一个包含 10 个隐藏单元和单个 LSTM 层的中等规模神经网络
包含 100 个隐藏单元和单个 LSTM 层的完整神经网络
包含 2 个 LSTM 层的深度神经网络（每个层有 100 个隐藏单元）

加载已训练好的神经网络。

load('variousHiddenUnitNets.mat')

生成一个阶跃信号。

stepTime = 2; % In seconds
stepAmplitude = 0.1;
stepDuration = 4; % In seconds

% Construct step signal.
time = (0:1/fs:stepDuration)';

stepSignal = [zeros(sum(time<=stepTime),1);stepAmplitude*ones(sum(time>stepTime),1)];
systemResponse = lsim(fourthOrderMdl,stepSignal,time);

% Transpose input and output signal for network inputs.
stepSignal = stepSignal';
systemResponse = systemResponse';

利用各种已训练好的神经网络对系统响应进行估计。

smallNetStep = predict(smallNet,stepSignal)-smallNetZeroMapping(end);
medNetStep = predict(medNet,stepSignal)-medNetZeroMapping(end);
fullnetStep = predict(fullNet,stepSignal) - fullNetZeroMapping(end);
doubleNetStep = predict(doubleNet,stepSignal) - doubleNetZeroMapping(end);

绘制估计响应值。

figure
title('Step response estimation')
plot(time,systemResponse,'k', ...
    time,doubleNetStep,'g', ...
    time,fullnetStep,'r', ...
    time,medNetStep,'c', ...
    time,smallNetStep,'b')
grid on
legend({'System','Double Net','Full Net','Med Net','Small Net'},'Location','northwest')
title('Fourth-Order Step')

Figure contains an axes object. The axes object with title Fourth-Order Step contains 5 objects of type line. These objects represent System, Double Net, Full Net, Med Net, Small Net.

请注意，所有网络都很好地捕捉到了响应中的高频动态特性。不过，为了比较该系统缓慢变化的动态特性，请绘制响应值的移动平均曲线。LSTM 捕捉线性系统长期动态特性（低频动态特性）的能力，与系统的动态特性以及 LSTM 中的隐藏单元数量直接相关。LSTM 中的层数与长期行为没有直接关系，而是增加了调整第一层估计结果的灵活性。

figure
title('Slow dynamics component')
plot(time,movmean(systemResponse,50),'k')
hold on
plot(time,movmean(doubleNetStep,50),'g')
plot(time,movmean(fullnetStep,50),'r')
plot(time,movmean(medNetStep,50),'c')
plot(time,movmean(smallNetStep,50),'b')
grid on
legend('System','Double Net','Full net','Med Net','Small Net','Location','northwest')
title('Fourth Order Step')

Figure contains an axes object. The axes object with title Fourth Order Step contains 5 objects of type line. These objects represent System, Double Net, Full net, Med Net, Small Net.

向测得的系统响应中添加噪声

在系统输出中添加随机噪声，以探究噪声对 LSTM 性能的影响。为此，请向测得的系统响应中添加强度分别为 1%、5% 和 10% 的白噪声。使用带噪声的数据来训练 LSTM 网络。使用相同的含噪声数据集，通过 tfest 估计线性模型。对这些模型进行仿真，并将仿真结果作为性能比较的基准。

使用与之前相同的阶跃函数：

stepTime = 2; % In seconds
stepAmplitude = 0.1;
stepDuration = 4; % In seconds

[stepSignal,systemResponse,time] = generateStepResponse(fourthOrderMdl,stepTime,stepAmplitude,stepDuration);

加载已训练好的神经网络，并估计系统响应。

load('noisyDataNetworks.mat')
netNoise1Step = predictAndAdjust(netNoise1,stepSignal,initializationSignal,initializationValue);
netNoise5Step = predictAndAdjust(netNoise5,stepSignal,initializationSignal,initializationValue);
netNoise10Step = predictAndAdjust(netNoise10,stepSignal,initializationSignal,initializationValue);

使用传递函数估计器 (tfest) 来估计上述噪声水平下的函数，以比较网络对噪声的鲁棒性（更多细节请参阅随附的 noiseLevelModels.m）。

load('noisyDataTFs.mat')
tfStepNoise1 = lsim(tfNoise1,stepSignal,time);
tfStepNoise5 = lsim(tfNoise5,stepSignal,time);
tfStepNoise10 = lsim(tfNoise10,stepSignal,time);

将生成的响应结果绘制成图。

figure
plot(time,systemResponse,'k', ...
    time,netNoise1Step, ...
    time,netNoise5Step, ...
    time,netNoise10Step)
grid on
legend('System Response','1% Noise','5% Noise','10% Noise')
title('Deep LSTM with noisy data')

Figure contains an axes object. The axes object with title Deep LSTM with noisy data contains 4 objects of type line. These objects represent System Response, 1% Noise, 5% Noise, 10% Noise.

现在，绘制估计的传递函数。

figure
plot(time,systemResponse,'k', ...
    time,tfStepNoise1, ...
    time,tfStepNoise5, ...
    time,tfStepNoise10)
grid on
legend('System Response','1% Noise','5% Noise','10% Noise')
title('Transfer functions fitted to noisy data')

Figure contains an axes object. The axes object with title Transfer functions fitted to noisy data contains 4 objects of type line. These objects represent System Response, 1% Noise, 5% Noise, 10% Noise.

计算均方误差，以便更好地评估不同模型在不同噪声水平下的性能。

msefun = @(y,yhat) mean(sqrt((y-yhat).^2)/length(y));

% LSTM errors
lstmMSE(1,:) = msefun(systemResponse,netNoise1Step);
lstmMSE(2,:) = msefun(systemResponse,netNoise5Step);
lstmMSE(3,:) = msefun(systemResponse,netNoise10Step);

% Transfer function errors
tfMSE(1,:) = msefun(systemResponse,tfStepNoise1');
tfMSE(2,:) = msefun(systemResponse,tfStepNoise5');
tfMSE(3,:) = msefun(systemResponse,tfStepNoise10');

mseTbl = array2table([lstmMSE tfMSE],'VariableNames',{'LSTMMSE','TFMSE'})

mseTbl=3×2 table
    1.0115e-05    8.8621e-07
    2.5577e-05    9.9064e-06
    5.1791e-05    3.6831e-05

噪声对 LSTM 和传递函数估计结果均产生了类似的影响。

辅助函数

function [stepSignal,systemResponse,time] = generateStepResponse(model,stepTime,stepAmp,signalDuration)
%Generates a step response for the given model.
%

%Check model type
modelType = class(model);
if nargin < 2
    stepTime = 1;
end

if nargin < 3
    stepAmp = 1;
end

if nargin < 4
    signalDuration = 10;
end

% Construct step signal
if model.Ts == 0
    Ts = 1e-2;
    time = (0:Ts:signalDuration)';
else
    time = (0:model.Ts:signalDuration)';
end
stepSignal = [zeros(sum(time<=stepTime),1);stepAmp*ones(sum(time>stepTime),1)];
switch modelType
    case {'tf', 'zpk'}
        systemResponse = lsim(model,stepSignal,time);
    case 'idpoly'
        systemResponse = sim(model,stepSignal,time);
    otherwise
        error('Model passed is not supported')
end

stepSignal = stepSignal';
systemResponse = systemResponse';
end