Этот пример показывает, как использовать нейронные сети с длительной кратковременной памятью (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 дБ от ее значения постоянного тока.
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 не напрямую связано с долгосрочным поведением, а скорее добавляет гибкость для корректировки оценки с первого уровня.
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