В этом примере показано, как использовать нейронные сети долгой краткосрочной памяти (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')
Диаграмма Боде показывает пропускную способность системы, которая измеряется как первая частота, где усиление опускается ниже 70,8% или приблизительно 3 дБ, его значения DC.
bodeplot(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')
График показывает две проблемы с подгонкой. Во-первых, начальное состояние сети не является стационарным, который приводит к переходному поведению в начале сигнала. Во-вторых, предсказание сети имеет небольшое смещение.
Чтобы инициализировать сетевое состояние к правильному начальному условию, необходимо обновить сетевое состояние, чтобы соответствовать состоянию системы в начале тестового сигнала.
Можно настроить начальное состояние сети путем сравнения предполагаемого ответа системы при начальном условии к фактическому ответу системы. Используйте различие между оценкой начального состояния сетью и фактическим ответом начального состояния, чтобы откорректировать для смещения по системной оценке.
Когда сеть выполняет оценку с помощью входа шага от 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')
Чтобы инициализировать сетевое состояние для нулевого входного сигнала, выберите входной сигнал нуля и выберите длительность так, чтобы сигнал был достаточно длинен для сетей, чтобы достигнуть устойчивого состояния.
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
Теперь, когда сеть правильно инициализируется, используйте сеть, чтобы предсказать переходной процесс снова и построить результаты. Начального воздействия не стало.
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')
Даже после того, как вы устанавливаете сетевое начальное состояние компенсировать начальное условие тестового сигнала, маленькое смещение все еще отображается в предсказанном ответе. Это вызвано тем, что неправильного срока смещения, что сеть 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')
Все сигналы, используемые, чтобы обучить сеть, имели максимальную амплитуду 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 histogram(urgs,'Normalization','pdf') grid on
Установите амплитуду ступенчатой функции согласно процентили распределения. Постройте коэффициент ошибок как функцию процентили.
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
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')
Как амплитуда перемещений переходного процесса за пределами области значений набора обучающих данных, 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')
Обратите внимание, что все сети получают высокочастотную динамику в ответе хорошо. Однако постройте скользящее среднее значение ответов для того, чтобы сравнить медленную различную динамику системы. Способность 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')
Добавьте случайный шум в систему выход, чтобы исследовать эффект шума на эффективности 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 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')
Вычислите среднеквадратическую ошибку, чтобы лучше оценить эффективность различных моделей на различном уровне шума.
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