利用LSTM网络进行线性系统辨识
这个例子展示了如何使用长短期记忆(LSTM)神经网络来估计线性系统，并将这种方法与传递函数估计进行了比较。
在本例中，您将研究LSTM网络捕获建模系统的底层动态的能力。为了做到这一点，你训练LSTM网络的输入和输出信号从一个线性传递函数，并测量网络响应的精度阶跃变化。
传递函数
本例使用四阶传递函数，具有混合快慢动态和中等阻尼。适度的阻尼使系统动力学在较长的时间范围内衰减，并显示了LSTM网络捕获混合动力学的能力，而不需要一些重要的响应动力学衰减。通过指定系统的极点和零点来构造传递函数。
fourthOrderMdl = zpk(4(9 + 5我;9-5i; 2 + 50我;2-50i), 5 e5);[stepResponse, stepTime] = (fourthOrderMdl)步;
画出传递函数的阶跃响应。
情节(stepTime stepResponse)网格在轴紧标题(“四阶混合步响应”）
波德图显示了系统的带宽，这是测量增益低于其直流值的70.8%(约3 dB)的第一个频率。
bodeplot (fourthOrderMdl)
fb =带宽(fourthOrderMdl)
fb = 62.8858
生成的训练数据
构造一个输入和输出信号的数据集，可以用来训练LSTM网络。对于输入，生成一个随机的高斯信号。模拟传递函数的响应fourthOrderMdl向此输入获得输出信号。
高斯噪声数据
指定随机高斯噪声训练信号的属性。
signalType =“该公司”；%高斯signalLength = 5000;信号中的点数fs = 100;%采样频率signalAmplitude = 1;最大信号振幅
生成高斯噪声信号使用idinput然后缩放结果。
开始= idinput (signalLength signalType);开始= (signalAmplitude / max(开始))*开始的;
根据采样率生成时间信号。
采用= 0:1 / fs:长度/ fs-1 / fs(开始);
使用lsim函数生成系统的响应，并将结果存储在yrgs．对模拟输出进行转置，使其对应于LSTM数据结构，该结构需要行向量而不是列向量。
yrgs = lsim (fourthOrderMdl开始,采用);yrgs = yrgs ';
类似地，创建一个较短的验证信号在网络训练期间使用。
signalType xval = idinput (100);yval = lsim (fourthOrderMdl xval,采用(1:10 0));
创建和培训网络
下面的网络架构是通过使用贝叶斯优化例程确定的，其中贝叶斯优化代价函数使用独立的验证数据(参见所附的bayesianOptimizationForLSTM.mlx的详细信息)。尽管多种体系结构可以工作，但这种优化提供了计算效率最高的一种。优化过程还表明，当将LSTM应用于其他线性传递函数时，随着传递函数复杂度的增加，网络的体系结构并没有发生显著变化。相反，训练网络所需的时数增加了。建模系统所需的隐藏单元的数量与动态衰减所需的时间有关。在这种情况下，响应有两个不同的部分:高频响应和低频响应。需要更多的隐藏单元来捕获低频响应。如果选择较低数量的单元，仍然模拟高频响应。然而，低频响应的估计变差。
创建网络架构。
numResponses = 1;featureDimension = 1;numHiddenUnits = 100;maxEpochs = 1000;miniBatchSize = 200;Networklayers = [sequenceInputLayer (featureDimension).．.lstmLayer (numHiddenUnits).．.lstmLayer (numHiddenUnits).．.fullyConnectedLayer (numResponses).．.regressionLayer];
初始学习率影响网络的成功。使用过高的初始学习率会导致高梯度，从而导致较长的训练时间。较长的训练时间会导致网络全连接层的饱和。当网络饱和时，输出发散，网络输出a南价值。因此，使用默认值0.01，这是一个相对较低的初始学习率。这导致了一个单调递减的残差和损失曲线。使用分段速率调度来避免优化算法在优化程序开始时陷入局部极小值。
选择= trainingOptions (“亚当”，.．.“MaxEpochs”maxEpochs,.．.“MiniBatchSize”miniBatchSize,.．.“GradientThreshold”10.．.“洗牌”，“一次”，.．.“阴谋”，“训练进步”，.．.“ExecutionEnvironment”，“图形”，.．.“LearnRateSchedule”，“分段”，.．.“LearnRateDropPeriod”, 100,.．.“详细”0,.．.“ValidationData”, {xval '} {yval}]);loadNetwork = true;%设置为true，使用并行池训练网络。如果loadNetwork负载(“fourthOrderMdlnet”，“fourthOrderNet”）其他的poolobj = parpool;fourthOrderNet = trainNetwork(开始、yrgs Networklayers,选项);删除(poolobj)保存(“fourthOrderMdlnet”，“fourthOrderNet”，“开始”，“yrgs”）;结束
评估模型的性能
如果网络能够成功捕获系统的动态行为，那么它的性能就会很好。通过测量网络准确预测系统对步进输入响应的能力来评估网络性能。
构造一个步进输入。
stepTime = 2;%在几秒钟内stepAmplitude = 0.1;stepDuration = 4;%在几秒钟内构造步进信号和系统输出。时间= (0:1 / fs: stepDuration) ';stepSignal = [0 (sum(时间< = stepTime) 1); stepAmplitude *的(和(时间> stepTime), 1)];systemResponse = lsim (fourthOrderMdl、stepSignal、时间);对网络输入的输入输出信号进行转置。stepSignal = stepSignal ';systemResponse = systemResponse ';
利用训练过的网络对系统响应进行评估。在一个图中比较系统和估计的响应。
fourthOrderMixedNetStep =预测(fourthOrderNet stepSignal);图的标题(“阶跃响应估计”)图(时间、systemResponse“k”四thordermixednetstep)网格在传奇(“系统”，“估计”)标题(“四阶一步”）
这幅图显示了合身的两个问题。首先，网络的初始状态不是平稳的，这导致了信号开始时的瞬态行为。第二，网络的预测有轻微的偏移。
初始化网络和调整适合度
为了将网络状态初始化到正确的初始条件，您必须更新网络状态以对应于测试信号开始时的系统状态。
通过比较系统在初始状态下的估计响应和系统的实际响应，可以调整网络的初始状态。利用网络对初始状态的估计与初始状态的实际响应之间的差值来修正系统估计中的偏移量。
设置网络初始状态
当网络使用从0到1的步进输入执行估计时，LSTM网络的状态(LSTM层的单元状态和隐藏状态)会向正确的初始条件漂移。方法可以在每个时间步骤提取单元格和网络的隐藏状态predictAndUpdateState函数。
只使用步骤之前的单元格和隐藏状态值，该步骤在2秒内发生。定义一个2秒的时间标记，并提取该标记之前的值。
stepMarker = time <= 2;yhat = 0 (sum (stepMarker), 1);hiddenState = 0(金额(stepMarker), 200);% 200 LSTM单位cellState = 0(金额(stepMarker), 200);为ntime = 1:sum(stepMarker) [fourthOrderNet,yhat(ntime)] = predictAndUpdateState(fourthOrderNet,stepSignal(ntime)');hiddenState (ntime:) = fourthOrderNet.Layers .HiddenState (2, 1);cellState (ntime:) = fourthOrderNet.Layers .CellState (2, 1);结束
接下来，绘制步骤前一段时间内的隐藏状态和单元状态，并确认它们收敛到固定值。
图subplot(2,1,1) plot(time(1:20 00)，hiddenState(1:20 00，:))网格在轴紧标题(“隐蔽的状态”) subplot(2,1,2) plot(time(1:200)，cellState(1:200，:))网格在轴紧标题(“细胞状态”）
为了初始化零输入信号的网络状态，选择一个零输入信号，并选择持续时间，使信号足够长，使网络达到稳态。
initializationSignalDuration = 10;%在几秒钟内initializationValue = 0;initializationSignal = initializationValue * 1 (1, initializationSignalDuration * fs);fourthOrderNet = predictAndUpdateState (fourthOrderNet initializationSignal);
验证初始条件为零或接近零。
图zeroMapping = predict(fourthOrderNet,initializationSignal);情节(zeroMapping)轴紧
现在网络已经正确初始化，再次使用网络预测步长响应并绘制结果。最初的干扰消失了。
fourthOrderMixedNetStep =预测(fourthOrderNet stepSignal);图的标题(“阶跃响应估计”)图(时间、systemResponse“k”，.．.时间、fourthOrderMixedNetStep“b”网格)在传奇(“系统”，“估计”)标题(“四阶阶跃-调整状态”）
调整网络抵消
即使在您设置网络初始状态以补偿测试信号的初始条件之后，在预测响应中仍然可以看到一个小的偏移。这是因为LSTM网络在训练过程中学习到错误的偏见术语。您可以通过使用用于更新网络状态的相同初始化信号来修复偏移量。初始化信号被期望将网络映射为零。零和网络估计之间的偏移量是网络学习到的偏差项中的误差。在每一层计算的偏差项的总和接近于在响应中检测到的偏差。然而，在网络输出处调整网络偏置项比在网络的每一层调整单个偏置项更容易。
偏见=意味着(预测(fourthOrderNet initializationSignal));fourthOrderMixedNetStep = fourthOrderMixedNetStep-bias;图的标题(“阶跃响应估计”)图(时间、systemResponse“k”、时间、fourthOrderMixedNetStep“b -”)传说(“系统”，“估计”)标题(“四阶步长-调整偏移量”）
离开训练场地
所有用于训练网络的信号的最大振幅为1，阶跃函数的振幅为0.1。现在，研究一下网络在这些范围之外的行为。
时移
通过调整步进的时间引入时移。将该步骤的时间设置为3秒，比训练集中的时间长1秒。绘制结果网络输出，注意输出被正确地延迟了1秒。
stepTime = 3;%在几秒钟内stepAmplitude = 0.1;stepDuration = 5;%在几秒钟内(stepSignal、systemResponse、时间)= generateStepResponse (fourthOrderMdl、stepTime stepAmplitude, stepDuration);fourthOrderMixedNetStep =预测(fourthOrderNet stepSignal);bias = fourthOrderMixedNetStep(1) - initializationValue;fourthOrderMixedNetStep = fourthOrderMixedNetStep-bias;图绘制(时间、systemResponse“k”、时间、fourthOrderMixedNetStep“b”网格)在轴紧
振幅键
接下来，增加阶跃函数的振幅，以研究当系统输入移动到训练数据范围之外时的网络行为。要测量训练数据范围之外的漂移，可以测量高斯噪声信号中振幅的概率密度函数。在直方图中可视化振幅。
图直方图(开始,“归一化”，“pdf”网格)在
根据分布的百分位数设置阶跃函数的振幅。将错误率绘制成百分比的函数。
pValues = [60:2:98, 90:1:98, 99:0.1:99.9 99.99];stepAmps = prctile(开始,pValues);%振幅stepTime = 3;%在几秒钟内stepDuration = 5;%在几秒钟内stepMSE = 0(长度(stepAmps), 1);fourthOrderMixedNetStep =细胞(长度(stepAmps), 1);步骤=细胞(长度(stepAmps), 1);为nAmps = 1:长度(stepAmps)%四阶混合(stepSignal、systemResponse、时间)= generateStepResponse (fourthOrderMdl、stepTime stepAmps (nAmps) stepDuration);fourthOrderMixedNetStep {nAmps} =预测(fourthOrderNet stepSignal);bias = fourthOrderMixedNetStep{nAmps}(1) - initializationValue;fourthOrderMixedNetStep {nAmps} = fourthOrderMixedNetStep {nAmps}偏见;stepMSE (nAmps) =√总和(systemResponse-fourthOrderMixedNetStep {nAmps}) ^ 2));步骤{nAmps 1} = systemResponse;结束图绘制(pValues stepMSE,“波”)标题(“预测误差是偏离训练范围的函数”网格)在轴紧
次要情节(2,1,1)情节(时间、步骤{1}“k”、时间、fourthOrderMixedNetStep {1},“b”网格)在轴紧标题(“最佳表现”)包含(“时间”) ylabel (系统响应的次要情节(2,1,2)情节({}结束时间、步骤,“k”、时间、fourthOrderMixedNetStep{}结束,“b”网格)在轴紧标题(表现最差的)包含(“时间”) ylabel (系统响应的）
当阶跃响应的幅度超出训练集的范围时，LSTM尝试估计响应的平均值。
这些结果表明，使用训练数据的重要性，在相同的范围内的数据将用于预测。否则，预测结果不可靠。
改变系统带宽
通过对四种不同网络的四阶混合动态传递函数建模，研究系统带宽对LSTM网络选择的隐藏单元数量的影响:
具有5个隐藏单元和单个LSTM层的小型网络
具有10个隐藏单元和单个LSTM层的中型网络
全网络包含100个隐藏单元和一个LSTM层
具有2个LSTM层的深度网络(每层有100个隐藏单元)
加载经过训练的网络。
负载(“variousHiddenUnitNets.mat”）
产生一个步进信号。
stepTime = 2;%在几秒钟内stepAmplitude = 0.1;stepDuration = 4;%在几秒钟内构造步进信号。时间= (0:1 / fs: stepDuration) ';stepSignal = [0 (sum(时间< = stepTime) 1); stepAmplitude *的(和(时间> stepTime), 1)];systemResponse = lsim (fourthOrderMdl、stepSignal、时间);对网络输入的输入输出信号进行转置。stepSignal = stepSignal ';systemResponse = systemResponse ';
使用各种训练过的网络估计系统响应。
smallNetStep =预测(smallNet stepSignal) -smallNetZeroMapping(结束);medNetStep =预测(medNet stepSignal) -medNetZeroMapping(结束);fullnetStep = predict(fullNet,stepSignal) - fullNetZeroMapping(end);doubleNetStep = predict(doubleNet,stepSignal) - doubleNetZeroMapping(end);
画出估计的响应。
图的标题(“阶跃响应估计”)图(时间、systemResponse“k”，.．.时间、doubleNetStep‘g’，.．.时间、fullnetStep“r”，.．.时间、medNetStep“c”，.．.时间、smallNetStep“b”网格)在传奇({“系统”，“双网”，“充分净”，“地中海净”，“小净”}，“位置”，“西北”)标题(“四阶一步”）
注意所有的网络都很好地捕捉了响应中的高频动态。但是，为了比较系统缓慢变化的动态，可以绘制响应的移动平均值。LSTM捕捉线性系统长期动态(低频动态)的能力与系统的动态和LSTM中隐藏单元的数量直接相关。LSTM中的层数与长期行为没有直接关系，而是增加了从第一层调整估计的灵活性。
图的标题(“慢动力学组件”)图(时间、movmean (systemResponse, 50),“k”)举行在情节(时间、movmean (doubleNetStep, 50),‘g’)图(时间、movmean (fullnetStep, 50),“r”)图(时间、movmean (medNetStep, 50),“c”)图(时间、movmean (smallNetStep, 50),“b”网格)在传奇(“系统”，“双网”，“充分净”，“地中海净”，“小净”，“位置”，“西北”)标题(“四阶一步”）
在测量的系统响应中添加噪声
在系统输出中加入随机噪声，探讨噪声对LSTM性能的影响。为此，在测量的系统响应中添加水平为1%、5%和10%的白噪声。利用噪声数据对LSTM网络进行训练。对于相同的噪声数据集，利用特遣部队．模拟这些模型，并使用模拟的响应作为性能比较的基线。
使用与前面相同的步长函数:
stepTime = 2;%在几秒钟内stepAmplitude = 0.1;stepDuration = 4;%在几秒钟内(stepSignal、systemResponse、时间)= generateStepResponse (fourthOrderMdl、stepTime stepAmplitude, stepDuration);
加载训练过的网络并估计系统响应。
负载(“noisyDataNetworks.mat”) netNoise1Step = predictAndAdjust(netNoise1,stepSignal,initializationSignal,initializationValue);netNoise5Step = predictAndAdjust (netNoise5 stepSignal、initializationSignal initializationValue);netNoise10Step = predictAndAdjust (netNoise10 stepSignal、initializationSignal initializationValue);
传递函数估计器(特遣部队)用于估计在上述噪声水平下的函数，以比较网络对噪声的弹性(见附带的noiseLevelModels.m更多的细节)。
负载(“noisyDataTFs.mat”) tfStepNoise1 = lsim(tfNoise1,stepSignal,time);tfStepNoise5 = lsim (tfNoise5、stepSignal、时间);tfStepNoise10 = lsim (tfNoise10、stepSignal、时间);
绘制生成的响应。
图绘制(时间、systemResponse“k”，.．.时间、netNoise1Step.．.时间、netNoise5Step.．.时间,netNoise10Step)网格在传奇(系统响应的，“1%噪声”，“5%噪声”，“10%噪声”)标题(“带噪声数据的深度LSTM”）
现在，画出估计的传递函数。
图绘制(时间、systemResponse“k”，.．.时间、tfStepNoise1.．.时间、tfStepNoise5.．.时间,tfStepNoise10)网格在传奇(系统响应的，“1%噪声”，“5%噪声”，“10%噪声”)标题(“适合噪声数据的传递函数”）
计算均方误差，以更好地评估不同模型在不同噪声水平下的性能。
msefun = @ (y, yhat)意味着(倍根号((y-yhat) ^ 2) /长度(y));% LSTM错误lstmMSE (1) = msefun (systemResponse netNoise1Step);lstmMSE (2) = msefun (systemResponse netNoise5Step);lstmMSE (3) = msefun (systemResponse netNoise10Step);%传递函数错误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表LSTMMSE TFMSE __________ __________ 5.1791 9.9064 2.5577 8.8621 1.0115 e-05 e-07 e-05 e-06 e-05 e-05 3.6831
噪声对LSTM和传递函数估计结果都有相似的影响。
辅助函数
函数(stepSignal、systemResponse、时间)= generateStepResponse(模型、stepTime stepAmp signalDuration)为给定的模型生成步长响应。％%检查模型类型modelType =类(模型);如果nargin < 2 stepTime = 1;结束如果nargin < 3 stepAmp = 1;结束如果nargin < 4 signalDuration = 10;结束构造步进信号如果模型。Ts = 0 Ts = 1e-2;时间= (0:Ts: signalDuration) ';其他的时间= (0:model.Ts: signalDuration) ';结束stepSignal = [0 (sum(时间< = stepTime) 1); stepAmp *的(和(时间> stepTime), 1)];开关modelType情况下｛“助教”，“zpk”} systemResponse = lsim(model,stepSignal,time);情况下“idpoly”systemResponse = sim(模型、stepSignal时间);否则错误(“不支持传递的模型”）结束stepSignal = stepSignal ';systemResponse = systemResponse ';结束
您也可以从以下列表中选择网站: