Основанная на подобии оценка остающегося полезного срока службы
В этом примере показано, как создать полный рабочий процесс оценки оставшейся полезной жизни (RUL), включая шаги для предварительной обработки, выбора тенденциозных функций, построения индикатора работоспособности путем слияния датчиков, обучения оценок RUL подобия и проверки эффективности прогностики. В примере используются обучающие данные из набора данных PHM2008 challenge из https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/ [1].
Подготовка данных
Поскольку набор данных является маленьким, возможно загрузить в память все данные о деградации. Загрузите и разархивируйте набор данных из https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/ в текущую директорию. Используйте helperLoadData вспомогательная функция для загрузки и преобразования обучающего текстового файла в массив ячеек с расписаниями. Обучающие данные содержат 218 симуляций запуска отказа. Эту группу измерений называют «ансамблем».
degradationData = helperLoadData('train.txt');
degradationData(1:5)ans = 5×1 cell array
{223×26 table}
{164×26 table}
{150×26 table}
{159×26 table}
{357×26 table}
Каждый представитель ансамбля представляет собой таблицу с 26 столбцами. Столбцы содержат данные для идентификатора машины, метки времени, 3 условий работы и 21 измерения датчика.
head(degradationData{1})ans=8×26 table
id time op_setting_1 op_setting_2 op_setting_3 sensor_1 sensor_2 sensor_3 sensor_4 sensor_5 sensor_6 sensor_7 sensor_8 sensor_9 sensor_10 sensor_11 sensor_12 sensor_13 sensor_14 sensor_15 sensor_16 sensor_17 sensor_18 sensor_19 sensor_20 sensor_21
__ ____ ____________ ____________ ____________ ________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________
1 1 10.005 0.2501 20 489.05 604.13 1499.5 1310 10.52 15.49 394.88 2318.9 8770.2 1.26 45.4 372.15 2388.1 8120.8 8.6216 0.03 368 2319 100 28.58 17.174
1 2 0.0015 0.0003 100 518.67 642.13 1584.5 1404 14.62 21.61 553.67 2388 9045.8 1.3 47.29 521.81 2388.2 8132.9 8.3907 0.03 391 2388 100 38.99 23.362
1 3 34.999 0.8401 60 449.44 555.42 1368.2 1122.5 5.48 8 194.93 2222.9 8343.9 1.02 41.92 183.26 2387.9 8063.8 9.3557 0.02 334 2223 100 14.83 8.8555
1 4 20.003 0.7005 0 491.19 607.03 1488.4 1249.2 9.35 13.65 334.82 2323.8 8721.5 1.08 44.26 314.84 2388.1 8052.3 9.2231 0.02 364 2324 100 24.42 14.783
1 5 42.004 0.8405 40 445 549.52 1354.5 1124.3 3.91 5.71 138.24 2211.8 8314.6 1.02 41.79 130.44 2387.9 8083.7 9.2986 0.02 330 2212 100 10.99 6.4025
1 6 20.003 0.7017 0 491.19 607.37 1480.5 1258.9 9.35 13.65 334.51 2323.9 8711.4 1.08 44.4 315.36 2388.1 8053.2 9.2276 0.02 364 2324 100 24.44 14.702
1 7 42 0.84 40 445 549.57 1354.4 1131.4 3.91 5.71 139.11 2211.8 8316.9 1.02 42.09 130.16 2387.9 8082 9.3753 0.02 331 2212 100 10.53 6.4254
1 8 0.0011 0 100 518.67 642.08 1589.5 1407.6 14.62 21.61 553.48 2388.1 9050.4 1.3 47.5 521.74 2388 8133.3 8.4339 0.03 391 2388 100 38.98 23.234
Разделите данные о деградации на набор обучающих данных и набор данных валидации для последующей эффективности оценки.
rng('default') % To make sure the results are repeatable numEnsemble = length(degradationData); numFold = 5; cv = cvpartition(numEnsemble, 'KFold', numFold); trainData = degradationData(training(cv, 1)); validationData = degradationData(test(cv, 1));
Задайте интересующие вас группы переменных.
varNames = string(degradationData{1}.Properties.VariableNames);
timeVariable = varNames{2};
conditionVariables = varNames(3:5);
dataVariables = varNames(6:26);Визуализация выборки данных ансамбля.
nsample = 10;
figure
helperPlotEnsemble(trainData, timeVariable, ...
[conditionVariables(1:2) dataVariables(1:2)], nsample)
Кластеризация рабочего режима
Как показано в предыдущем разделе, нет четкого тренда, показывающей процесс деградации при каждом измерении на отказ. В этом и следующем разделах рабочие условия будут использоваться для извлечения более четких трендов деградации из сигналов датчика.
Заметьте, что каждый представитель ансамбля содержит 3 условия работы: «op_setting_1,» «op_setting_2,» и «op_setting_3.» Во-первых, давайте извлекем таблицу из каждой камеры и соединим их в одну таблицу.
trainDataUnwrap = vertcat(trainData{:});
opConditionUnwrap = trainDataUnwrap(:, cellstr(conditionVariables));Визуализируйте все рабочие точки на 3D графике поля точек. В нем четко показано 6 режимов, и точки в каждом режиме находятся в очень близкой близости.
figure helperPlotClusters(opConditionUnwrap)
Давайте будем использовать методы кластеризации, чтобы автоматически найти 6 кластеров. Здесь используется алгоритм K-средних значений. K-средних значений являются одним из самых популярных алгоритмов кластеризации, но это может привести к локальному оптимуму. Это хорошая практика, чтобы повторить алгоритм кластеризации K-означает несколько раз с различными начальными условиями и выбрать результаты с самой низкой стоимостью. В этом случае алгоритм запускается 5 раз, и результаты идентичны.
opts = statset('Display', 'final'); [clusterIndex, centers] = kmeans(table2array(opConditionUnwrap), 6, ... 'Distance', 'sqeuclidean', 'Replicates', 5, 'Options', opts);
Replicate 1, 1 iterations, total sum of distances = 0.279547. Replicate 2, 1 iterations, total sum of distances = 0.279547. Replicate 3, 1 iterations, total sum of distances = 0.279547. Replicate 4, 1 iterations, total sum of distances = 0.279547. Replicate 5, 1 iterations, total sum of distances = 0.279547. Best total sum of distances = 0.279547
Визуализируйте результаты кластеризации и идентифицированные центроиды кластера.
figure helperPlotClusters(opConditionUnwrap, clusterIndex, centers)
Как иллюстрирует график, алгоритм кластеризации успешно находит 6 рабочих режимов.
Нормализация рабочего режима
Давайте выполним нормализацию по измерениям, сгруппированным по различным рабочим режимам. Во-первых, вычислите среднее и стандартное отклонение каждого измерения датчика, сгруппированного по рабочим режимам, идентифицированным в последнем разделе.
centerstats = struct('Mean', table(), 'SD', table()); for v = dataVariables centerstats.Mean.(char(v)) = splitapply(@mean, trainDataUnwrap.(char(v)), clusterIndex); centerstats.SD.(char(v)) = splitapply(@std, trainDataUnwrap.(char(v)), clusterIndex); end centerstats.Mean
ans=6×21 table
sensor_1 sensor_2 sensor_3 sensor_4 sensor_5 sensor_6 sensor_7 sensor_8 sensor_9 sensor_10 sensor_11 sensor_12 sensor_13 sensor_14 sensor_15 sensor_16 sensor_17 sensor_18 sensor_19 sensor_20 sensor_21
________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________
489.05 604.92 1502.1 1311.4 10.52 15.493 394.32 2319 8784.1 1.26 45.496 371.44 2388.2 8133.9 8.6653 0.03 369.74 2319 100 28.525 17.115
518.67 642.71 1590.7 1409.4 14.62 21.61 553.3 2388.1 9062.3 1.3 47.557 521.36 2388.1 8140.9 8.4442 0.03 393.27 2388 100 38.809 23.285
462.54 536.87 1262.8 1050.6 7.05 9.0275 175.4 1915.4 8014.9 0.93989 36.808 164.56 2028.3 7878 10.916 0.02 307.39 1915 84.93 14.262 8.5552
445 549.72 1354.7 1128.2 3.91 5.7158 138.62 2212 8327 1.0202 42.163 130.53 2388 8088.4 9.3775 0.02 331.12 2212 100 10.584 6.3515
491.19 607.59 1486 1253.6 9.35 13.657 334.46 2324 8729.1 1.0777 44.466 314.85 2388.2 8065 9.2347 0.022299 365.45 2324 100 24.447 14.67
449.44 555.82 1366.9 1131.9 5.48 8.0003 194.43 2223 8355.2 1.0203 41.995 183.01 2388.1 8071.1 9.3344 0.02 334.29 2223 100 14.827 8.8966
centerstats.SD
ans=6×21 table
sensor_1 sensor_2 sensor_3 sensor_4 sensor_5 sensor_6 sensor_7 sensor_8 sensor_9 sensor_10 sensor_11 sensor_12 sensor_13 sensor_14 sensor_15 sensor_16 sensor_17 sensor_18 sensor_19 sensor_20 sensor_21
__________ ________ ________ ________ __________ _________ ________ ________ ________ __________ _________ _________ _________ _________ _________ __________ _________ _________ __________ _________ _________
1.4553e-11 0.47617 5.8555 8.3464 1.1618e-12 0.0047272 0.65536 0.094487 18.057 1.51e-13 0.25089 0.52838 0.096183 16.012 0.037598 9.9929e-16 1.4723 0 0 0.14522 0.086753
4.059e-11 0.48566 5.9258 8.8223 3.7307e-14 0.0011406 0.86236 0.068654 20.061 1.2103e-13 0.25937 0.71239 0.069276 17.212 0.036504 9.9929e-16 1.5046 0 0 0.17565 0.10515
2.7685e-11 0.35468 5.2678 6.9664 2.043e-14 0.0043301 0.45074 0.2743 14.741 0.0010469 0.21539 0.33985 0.28932 13.624 0.044142 8.6744e-16 1.3083 0 5.9833e-12 0.11212 0.06722
0 0.44169 5.6853 7.5741 1.9763e-13 0.0049401 0.44198 0.30713 18.389 0.0014951 0.23584 0.3432 0.33144 16.917 0.037135 2.231e-15 1.4174 0 0 0.10778 0.063991
4.4456e-11 0.46992 5.7664 7.8679 8.9892e-13 0.0046633 0.59984 0.13032 17.983 0.0042388 0.2362 0.49055 0.13285 15.792 0.038156 0.0042084 1.4428 0 0 0.13352 0.079731
4.3489e-11 0.44341 5.7224 7.4842 3.7485e-13 0.0017642 0.4734 0.28889 17.608 0.0017978 0.23242 0.38243 0.31 15.91 0.038624 8.5009e-16 1.3974 0 0 0.11311 0.069348
Статистика в каждом режиме может использоваться для нормализации обучающих данных. Для каждого представителя ансамбля извлеките рабочие точки каждой строки, вычислите ее расстояние до каждого центра кластера и найдите ближайший центр кластера. Затем для каждого измерения датчика вычитайте среднее значение и разделите его на стандартное отклонение этого кластера. Если стандартное отклонение близко к 0, установите нормализованное измерение датчика равным 0, потому что почти постоянное измерение датчика не полезно для оставшейся оценки полезного срока службы. Для получения дополнительной информации о функции modeNormalization см. последний раздел «Вспомогательные функции».
trainDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ... trainData, 'UniformOutput', false);
Визуализируйте данные, нормированные рабочим режимом. Тренды деградации для некоторых измерений датчика теперь выявляются после нормализации.
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(1:4), nsample)
Анализ тенденций
Теперь выберите самые модные измерения датчика, чтобы создать индикатор работоспособности для предсказания. Для каждого измерения датчика оценивают линейную модель деградации и ранжируют наклоны сигналов.
numSensors = length(dataVariables); signalSlope = zeros(numSensors, 1); warn = warning('off'); for ct = 1:numSensors tmp = cellfun(@(tbl) tbl(:, cellstr(dataVariables(ct))), trainDataNormalized, 'UniformOutput', false); mdl = linearDegradationModel(); % create model fit(mdl, tmp); % train mode signalSlope(ct) = mdl.Theta; end warning(warn);
Отсортируйте откосы сигнала и выберите 8 датчиков с наибольшими уклонами.
[~, idx] = sort(abs(signalSlope), 'descend');
sensorTrended = sort(idx(1:8))sensorTrended = 8×1
2
3
4
7
11
12
15
17
Визуализируйте выбранные модные измерения датчика.
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(sensorTrended(3:6)), nsample)
Заметьте, что некоторые из наиболее модных сигналов показывают положительные тренды, в то время как другие показывают отрицательные тренды.
Создайте индикатор работоспособности
Этот раздел посвящен сплавлению измерений датчика в один индикатор работоспособности, с которым обучена основанная на подобии модель.
Все данные пробега до отказа приняты, чтобы начать со здорового условия. Условию здоровья в начале присваивается значение 1, и условию здоровья при отказе присваивается значение 0. Условие здоровья принято линейно деградирующим от 1 до 0 с течением времени. Это линейное ухудшение используется, чтобы помочь сплавить значения датчика. Более сложные методы слияния датчиков описаны в литературе [2-5].
for j=1:numel(trainDataNormalized) data = trainDataNormalized{j}; rul = max(data.time)-data.time; data.health_condition = rul / max(rul); trainDataNormalized{j} = data; end
Визуализируйте условие здоровья.
figure
helperPlotEnsemble(trainDataNormalized, timeVariable, "health_condition", nsample)
Условие здоровья всех представителей ансамбля изменяется с 1 до 0 с различными ухудшающимися скоростями.
Теперь подгоните линейную регрессионую модель условия здоровья с наиболее трендовыми измерениями датчика в качестве регрессоров:
Состояние здоровья ~ 1 + Sensor2 + Sensor3 + Sensor4 + Sensor7 + Sensor11 + Sensor12 + Sensor15 + Sensor17
trainDataNormalizedUnwrap = vertcat(trainDataNormalized{:});
sensorToFuse = dataVariables(sensorTrended);
X = trainDataNormalizedUnwrap{:, cellstr(sensorToFuse)};
y = trainDataNormalizedUnwrap.health_condition;
regModel = fitlm(X,y);
bias = regModel.Coefficients.Estimate(1)bias = 0.5000
weights = regModel.Coefficients.Estimate(2:end)
weights = 8×1
-0.0308
-0.0308
-0.0536
0.0033
-0.0639
0.0051
-0.0408
-0.0382
Создайте один индикатор работоспособности, умножив измерения датчика на связанные с ними веса.
trainDataFused = cellfun(@(data) degradationSensorFusion(data, sensorToFuse, weights), trainDataNormalized, ... 'UniformOutput', false);
Визуализируйте слитый индикатор работоспособности для обучающих данных.
figure helperPlotEnsemble(trainDataFused, [], 1, nsample) xlabel('Time') ylabel('Health Indicator') title('Training Data')
Данные с нескольких датчиков объединяются в один индикатор работоспособности. Индикатор работоспособности сглаживается фильтром скользящего среднего. Для получения дополнительной информации см. вспомогательную функцию «degradationSensorFusion» в последнем разделе.
Применить ту же операцию к данным валидации
Повторите процесс нормализации режима и слияния датчиков с набором данных валидации.
validationDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ... validationData, 'UniformOutput', false); validationDataFused = cellfun(@(data) degradationSensorFusion(data, sensorToFuse, weights), ... validationDataNormalized, 'UniformOutput', false);
Визуализируйте индикатор работоспособности для данных валидации.
figure helperPlotEnsemble(validationDataFused, [], 1, nsample) xlabel('Time') ylabel('Health Indicator') title('Validation Data')
Построение модели RUL подобия
Теперь создайте модель RUL подобия на основе остатков, используя обучающие данные. В этой настройке модель пытается подогнать все слитые данные с помощью полинома 2-го порядка.
Расстояние между данными и данные вычисляется 1-нормой невязки
где является индикатором работоспособности машины , - предполагаемый показатель работоспособности машины использование полиномиальной модели 2-го порядка, идентифицированной в машине .
Счет подобия вычисляется следующей формулой
Учитывая один представителя ансамбля в наборе данных валидации, модель найдет ближайший 50 представитель ансамбля в наборе обучающих данных, подбирает распределение вероятностей на основе 50 представителей ансамбля и использует медиану распределения в качестве оценки RUL.
mdl = residualSimilarityModel(... 'Method', 'poly2',... 'Distance', 'absolute',... 'NumNearestNeighbors', 50,... 'Standardize', 1); fit(mdl, trainDataFused);
Оценка эффективности
Чтобы оценить модель RUL подобия, используйте 50%, 70% и 90% выборочных данных валидации, чтобы предсказать ее RUL.
breakpoint = [0.5, 0.7, 0.9];
validationDataTmp = validationDataFused{3}; % use one validation data for illustrationИспользуйте данные валидации перед первой точкой останова, которая составляет 50% от срока службы.
bpidx = 1; validationDataTmp50 = validationDataTmp(1:ceil(end*breakpoint(bpidx)),:); trueRUL = length(validationDataTmp) - length(validationDataTmp50); [estRUL, ciRUL, pdfRUL] = predictRUL(mdl, validationDataTmp50);
Визуализируйте данные валидации, усеченные на 50%, и его ближайшие соседи.
figure compare(mdl, validationDataTmp50);
Визуализируйте предполагаемый RUL по сравнению с истинным RUL и распределением вероятностей предполагаемого RUL.
figure helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)
Существует относительно большая ошибка между предполагаемым RUL и истинным RUL, когда машина находится на промежуточной стадии работоспособности. В этом примере самые похожие 10 кривых близки в начале, но раздваиваются, когда они приближаются к состоянию отказа, получая примерно два режима в распределении RUL.
Используйте данные валидации перед второй точкой останова, которая составляет 70% от срока службы.
bpidx = 2; validationDataTmp70 = validationDataTmp(1:ceil(end*breakpoint(bpidx)), :); trueRUL = length(validationDataTmp) - length(validationDataTmp70); [estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp70); figure compare(mdl, validationDataTmp70);
figure helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)
Когда наблюдается больше данных, улучшается оценка RUL.
Используйте данные валидации перед третьей точкой останова, которая составляет 90% от срока службы.
bpidx = 3; validationDataTmp90 = validationDataTmp(1:ceil(end*breakpoint(bpidx)), :); trueRUL = length(validationDataTmp) - length(validationDataTmp90); [estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp90); figure compare(mdl, validationDataTmp90);
figure helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)
Когда машина близка к отказу, оценка RUL еще больше улучшается в этом примере.
Теперь повторите ту же процедуру оценки для всего набора данных валидации и вычислите ошибку между предполагаемым RUL и истинным RUL для каждой точки останова.
numValidation = length(validationDataFused); numBreakpoint = length(breakpoint); error = zeros(numValidation, numBreakpoint); for dataIdx = 1:numValidation tmpData = validationDataFused{dataIdx}; for bpidx = 1:numBreakpoint tmpDataTest = tmpData(1:ceil(end*breakpoint(bpidx)), :); trueRUL = length(tmpData) - length(tmpDataTest); [estRUL, ~, ~] = predictRUL(mdl, tmpDataTest); error(dataIdx, bpidx) = estRUL - trueRUL; end end
Визуализируйте гистограмму ошибки для каждой точки останова вместе с ее распределением вероятностей.
[pdf50, x50] = ksdensity(error(:, 1)); [pdf70, x70] = ksdensity(error(:, 2)); [pdf90, x90] = ksdensity(error(:, 3)); figure ax(1) = subplot(3,1,1); hold on histogram(error(:, 1), 'BinWidth', 5, 'Normalization', 'pdf') plot(x50, pdf50) hold off xlabel('Prediction Error') title('RUL Prediction Error using first 50% of each validation ensemble member') ax(2) = subplot(3,1,2); hold on histogram(error(:, 2), 'BinWidth', 5, 'Normalization', 'pdf') plot(x70, pdf70) hold off xlabel('Prediction Error') title('RUL Prediction Error using first 70% of each validation ensemble member') ax(3) = subplot(3,1,3); hold on histogram(error(:, 3), 'BinWidth', 5, 'Normalization', 'pdf') plot(x90, pdf90) hold off xlabel('Prediction Error') title('RUL Prediction Error using first 90% of each validation ensemble member') linkaxes(ax)
Постройте график ошибки предсказания на прямоугольном графике, чтобы визуализировать медиану, 25-75 квантиль и выбросы.
figure boxplot(error, 'Labels', {'50%', '70%', '90%'}) ylabel('Prediction Error') title('Prediction error using different percentages of each validation ensemble member')
Вычислите и визуализируйте среднее и стандартное отклонение ошибки предсказания.
errorMean = mean(error)
errorMean = 1×3
-5.8944 3.1359 3.3555
errorMedian = median(error)
errorMedian = 1×3
-4.8538 5.3763 3.6580
errorSD = std(error)
errorSD = 1×3
26.4916 20.0720 18.0313
figure errorbar([50 70 90], errorMean, errorSD, '-o', 'MarkerEdgeColor','r') xlim([40, 100]) xlabel('Percentage of validation data used for RUL prediction') ylabel('Prediction Error') legend('Mean Prediction Error with 1 Standard Deviation Eror bar', 'Location', 'south')
Показано, что ошибка становится более концентрированной около 0 (меньше выбросов), так как наблюдается больше данных.
Ссылки
[1] А. Саксена и К. Гебель (2008). «PHM08 Challenge Data Set», NASA Ames Prognostics Data Repository (http://ti.arc.nasa.gov/project/prognostic-data-repository), Исследовательский центр НАСА Эймса, Моффетт-Филд, Калифорния
[2] Roemer, Michael J., Gregory J. Kacprzynski, and Michael H. Schoeller. «Улучшенные диагностические и прогностические оценки с использованием слияния информации управления здравоохранением». AUTOTESTCON Proceedings, 2001. IEEE Systems Readiness Technology Conference. IEEE, 2001.
[3] Гебель, Кай и Пьеро Бониссоне. «Прогностическое информационное слияние для систем постоянной нагрузки». Information Fusion, 2005 8-я Международная конференция. Том 2. IEEE, 2005.
[4] Ван, Пэн и Дэвид У. Койт. «Надежность предсказания основанная на моделировании деградации для систем с несколькими показателями деградации». Надежность и поддерживаемость, 2004 ежегодный симпозиум-РАМН. IEEE, 2004.
[5] Jardine, Andrew KS, Daming Lin, and Dragan Banjevic. «Обзор диагностики машин и прогнозирования, реализующий техническое обслуживание на основе условий». Механические системы и обработка сигналов 20.7 (2006): 1483-1510.
Вспомогательные функции
function data = regimeNormalization(data, centers, centerstats) % Perform normalization for each observation (row) of the data % according to the cluster the observation belongs to. conditionIdx = 3:5; dataIdx = 6:26; % Perform row-wise operation data{:, dataIdx} = table2array(... rowfun(@(row) localNormalize(row, conditionIdx, dataIdx, centers, centerstats), ... data, 'SeparateInputs', false)); end function rowNormalized = localNormalize(row, conditionIdx, dataIdx, centers, centerstats) % Normalization operation for each row. % Get the operating points and sensor measurements ops = row(1, conditionIdx); sensor = row(1, dataIdx); % Find which cluster center is the closest dist = sum((centers - ops).^2, 2); [~, idx] = min(dist); % Normalize the sensor measurements by the mean and standard deviation of the cluster. % Reassign NaN and Inf to 0. rowNormalized = (sensor - centerstats.Mean{idx, :}) ./ centerstats.SD{idx, :}; rowNormalized(isnan(rowNormalized) | isinf(rowNormalized)) = 0; end function dataFused = degradationSensorFusion(data, sensorToFuse, weights) % Combine measurements from different sensors according % to the weights, smooth the fused data and offset the data % so that all the data start from 1 % Fuse the data according to weights dataToFuse = data{:, cellstr(sensorToFuse)}; dataFusedRaw = dataToFuse*weights; % Smooth the fused data with moving mean stepBackward = 10; stepForward = 10; dataFused = movmean(dataFusedRaw, [stepBackward stepForward]); % Offset the data to 1 dataFused = dataFused + 1 - dataFused(1); end