Используйте сеть LSTM для идентификации линейной системы

В этом примере показано, как использовать нейронные сети долгой краткосрочной памяти (LSTM), чтобы оценить линейную систему и сравнивает этот подход к оценке передаточной функции.

В этом примере вы исследуете способность сети LTSM получить базовую динамику смоделированной системы. Для этого вы обучаете сеть 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. The axes with title Fourth-order mixed step response contains an object of type line.

Диаграмма Боде показывает пропускную способность системы, которая измеряется как первая частота, где усиление опускается ниже 70,8% или приблизительно 3 дБ, его значения DC.

bodeplot(fourthOrderMdl)

Figure contains 2 axes. Axes 1 contains an object of type line. This object represents fourthOrderMdl. Axes 2 contains an object of type line. This object represents fourthOrderMdl.

fb = bandwidth(fourthOrderMdl)
fb = 62.8858

Сгенерируйте обучающие данные

Создайте набор данных сигналов ввода и вывода, что можно использовать, чтобы обучить сеть LSTM. Для входа сгенерируйте случайный Гауссов сигнал. Симулируйте ответ передаточной функции fourthOrderMdl к этому входу, чтобы получить выходной сигнал.

Гауссовы шумовые данные

Задайте свойства для случайного Гауссова шумового учебного сигнала.

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.mlx для деталей). Несмотря на то, что несколько архитектур могут работать, эта оптимизация обеспечивает наиболее в вычислительном отношении эффективный. Процесс оптимизации также показал это сложностью увеличений передаточной функции при применении LSTM к другим линейным передаточным функциям, архитектура сети значительно не изменяется. Скорее номер эпох должен был обучить сетевые увеличения. Количество скрытых модулей, требуемых для моделирования системы, связано с тем, сколько времени движущие силы берут к ослаблению. В этом случае существует две отличных части к ответу: высокочастотный ответ и низкочастотный ответ. Более высокое количество скрытых модулей требуется, чтобы получать низкочастотный ответ. Если более низкое количество модулей выбрано, высокочастотный ответ все еще моделируется. Однако оценка низкочастотного ответа ухудшается.

Создайте сетевую архитектуру.

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

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

Начальная скорость обучения влияет на успех сети. Используя начальную скорость обучения, которая является слишком высокими результатами в высоких градиентах, которые приводят к более длительным учебным временам. Более длительные учебные времена могут привести к насыщению полносвязного слоя сети. Когда сеть насыщает, выходные параметры отличаются и сетевые выходные параметры NaN значение. Следовательно, используйте значение по умолчанию 0,01, который является относительно низкой начальной скоростью обучения. Это приводит к монотонно уменьшающейся невязке и кривым потерь. Используйте кусочное расписание уровня, чтобы помешать алгоритму оптимизации то, чтобы быть захваченным в локальных минимумах в начале стандартной программы оптимизации.

options = trainingOptions('adam', ...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'GradientThreshold',10, ...
    'Shuffle','once', ...
    'Plots','training-progress',...
    'ExecutionEnvironment','gpu',...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropPeriod',100,...
    'Verbose',0,...
    'ValidationData',[{xval'} {yval'}]);

loadNetwork = true; % Set to true to train the network using a parallel pool.
if loadNetwork
    load('fourthOrderMdlnet','fourthOrderNet')
else
    poolobj = parpool;
    fourthOrderNet = trainNetwork(urgs,yrgs,Networklayers,options);
    delete(poolobj)
    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. The axes 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. Axes 1 with title Hidden State contains 200 objects of type line. Axes 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. The axes 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. The axes 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. The axes 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. The axes contains 2 objects of type line.

Амплитудный сдвиг

Затем увеличьте амплитуду ступенчатой функции, чтобы исследовать сетевое поведение как системное входное перемещение за пределами области значений обучающих данных. Чтобы измерить дрейф за пределами области значений обучающих данных, можно измерить функцию плотности вероятности амплитуд в Гауссовом шумовом сигнале. Визуализируйте амплитуды в гистограмме.

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

Figure contains an axes. The axes 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. The axes with title Prediction Error as a Function of Deviation from Training Rrange contains an object of type line.

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. Axes 1 with title Best Performance contains 2 objects of type line. Axes 2 with title Worst Performance 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. The axes 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. The axes 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. The axes 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. The axes 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
     LSTMMSE        TFMSE   
    __________    __________

    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

% Constuct 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
Для просмотра документации необходимо авторизоваться на сайте