Набор данных

Набор данных содержит измерения от 124 литий-ионных элементов номинальной ёмкостью 1,1 А· ч и номинальным напряжением 3,3 В при различных профилях заряда и разряда. Полный набор данных можно найти здесь [2] с подробным описанием здесь [1]. В этом примере данные обрабатываются, чтобы содержать только подмножество измерений, относящихся к извлекаемым элементам. Это уменьшает размер загрузки без изменения производительности рассматриваемой модели машинного обучения. Обучающие данные содержат измерения из 41 ячейки, валидационные данные содержат измерения из 43 ячеек и тестовые данные содержат измерения из 40 ячеек. Данные для каждой ячейки хранятся в структуре, включающей следующую информацию:

Загрузите данные с сайта файлов поддержки MathWorks (это большой набор данных, ~ 1,7 ГБ).

url = 'https://ssd.mathworks.com/supportfiles/predmaint/batterycyclelifeprediction/v1/batteryDischargeData.zip';
websave('batteryDischargeData.zip',url);
unzip('batteryDischargeData.zip')
load('batteryDischargeData');

Извлечение элементов

Разрядная емкость является ключевой особенностью, которая может отражать работоспособность аккумулятора, и ее значение указывает дальность движения в электромобиле. Постройте график разрядной способности для первых 1000 циклов клеток в обучающих данных для визуализации ее изменения в течение срока службы клетки.

figure, hold on;
for i = 1:size(trainData,2)
    if numel(trainData(i).summary.cycle) == numel(unique(trainData(i).summary.cycle))
      plot(trainData(i).summary.cycle, trainData(i).summary.QDischarge);      
    end
end
ylim([0.85,1.1]),xlim([0,1000]);
ylabel('Discharge capacity (Ah)');
xlabel('cycle');

Как видно на сюжете, угасание мощностей ускоряется ближе к концу жизни. Однако затухание емкости ничтожно мало в первых 100 циклах и само по себе не является хорошей особенностью для прогнозирования срока службы аккумулятора. Поэтому подход, основанный на данных, который учитывает кривые напряжения каждого цикла, наряду с дополнительными измерениями, включая внутреннее сопротивление ячейки и температуру, рассматривается для прогнозирования оставшегося срока службы цикла. Важным характерным пространством является межцикловая эволюция Q (V) - кривой напряжения разряда как функции напряжения для данного цикла [1]. В частности, изменение кривых напряжения разряда между циклами, $\Delta Q($V), очень эффективно в диагностике деградации. Поэтому выберите разность разрядной емкости как функцию напряжения между 100-м и 10-м циклом$$,{}_{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}}\mathrm{\Delta Q100-10}$$( В).

figure, hold on;
for i = 1:size(testData,2)
    plot((testData(i).cycles(100).Qdlin - testData(i).cycles(10).Qdlin), testData(i).Vdlin)   
end
ylabel('Voltage (V)'); xlabel('Q_{100} - Q_{10} (Ah)');

Вычислите сводную статистику, используя $${}_{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{\Delta Q100-10}}(V$$) кривые каждой ячейки в качестве индикатора условия. Важно иметь в виду, что эти сводные статистические данные не должны иметь четкого физического смысла, они должны просто демонстрировать возможности прогнозирования. Например, дисперсия $${}_{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}}\Delta Q100-10($$V) и срок службы цикла сильно коррелированы, как видно на этом логарифмическом графике, тем самым подтверждая этот подход.

figure;
trainDataSize = size(trainData,2);
cycleLife = zeros(trainDataSize,1);
DeltaQ_var = zeros(trainDataSize,1);
for i = 1:trainDataSize
    cycleLife(i) = size(trainData(i).cycles,2) + 1;
    DeltaQ_var(i) = var(trainData(i).cycles(100).Qdlin - trainData(i).cycles(10).Qdlin);
end
loglog(DeltaQ_var,cycleLife, 'o')
ylabel('Cycle life'); xlabel('Var(Q_{100} - Q_{10}) (Ah^2)');

Затем извлеките следующие характеристики из необработанных данных измерений (имена переменных приведены в скобках) [1]:

helperGetFeatures функция принимает исходные данные измерений и вычисляет перечисленные выше характеристики. Возвращается XTrain содержащие элементы в таблице и yTrain содержащий ожидаемый оставшийся срок службы цикла.

[XTrain,yTrain] = helperGetFeatures(trainData);
head(XTrain)

ans=8×8 table
    DeltaQ_var    DeltaQ_min    CapFadeCycle2Slope    CapFadeCycle2Intercept     Qd2      AvgChargeTime     MinIR      IRDiff2And100
    __________    __________    __________________    ______________________    ______    _____________    ________    _____________

     -5.0839       -1.9638          6.4708e-06                1.0809            1.0753       13.409        0.016764     -3.3898e-05 
     -4.3754       -1.6928          1.6313e-05                1.0841            1.0797       12.025        0.016098      4.4186e-05 
     -4.1464       -1.5889          8.1708e-06                  1.08            1.0761       10.968        0.015923     -0.00012443 
     -3.8068       -1.4216          -8.491e-06                1.0974            1.0939       10.025        0.016083     -3.7309e-05 
     -4.1181       -1.6089          2.2859e-05                1.0589            1.0538       11.669        0.015963     -0.00030445 
     -4.0225       -1.5407          2.5969e-05                1.0664            1.0611       10.798        0.016543     -0.00024655 
     -3.9697       -1.5077          1.7886e-05                1.0762            1.0721       10.147        0.016238      2.2163e-05 
     -3.6195       -1.3383         -1.0356e-05                1.0889            1.0851       9.9247        0.016205     -6.6087e-05 

Разработка модели: Линейная регрессия с эластичной регуляризацией сетки

Используйте линейную модель с регуляризацией для прогнозирования оставшейся продолжительности цикла ячейки. Линейные модели имеют низкую вычислительную стоимость, но обеспечивают высокую интерпретируемость. Линейная модель имеет вид:

где $${}_{\mathrm{yi}}$$ - число циклов для ячейки $$i$$, $${}_{\mathrm{xi}}$$ - $$$$p-мерный вектор признаков для ячейки $$i$$ и $$w$$ - $$$$p-мерный вектор коэффициентов модели, а $$\beta$$ - скалярный перехват. Линейная модель регуляризуется с использованием упругой сетки для устранения высоких корреляций между признаками. Для $$\alpha$$ строго между 0 и 1 и неотрицательным $$\lambda$$коэффициент регуляризации, упругая сетка решает задачу:

где $${}_{\mathrm{P\alpha }}(w)\frac{\mathrm{=}(}{\mathrm{}}1-\alpha ){}_{\mathrm{}}^{\mathrm{}}{}_{\mathrm{}}$$2‖w‖22+α‖w‖1. Упругая сетка приближается к регрессии гребня, $$$$когда α близок к 0, и это то же самое, что регуляризация лассо, $$\mathrm{когда}$$α = 1. Обучающие данные используются для выбора $$$$гиперпараметров $$\alpha$$и λ и определения значений $$$$коэффициентов, w.lasso функция возвращает аппроксимированные коэффициенты регрессии методом наименьших квадратов и информацию о соответствии модели в качестве выходных данных для каждого $$$$$$$$значения α и λ. lasso функция принимает вектор значений для $$$$параметра λ. Поэтому для каждого значения $$\alpha$$$$$$$$$$вычисляются модельные коэффициенты w и λ, дающие наименьшее значение RMSE. Как предложено в [1], используйте четырехкратную перекрестную проверку с повторением 1 Монте-Карло для перекрестной проверки.

rng("default")

alphaVec = 0.01:0.1:1;
lambdaVec = 0:0.01:1;
MCReps = 1;
cvFold = 4;

rmseList = zeros(length(alphaVec),1);
minLambdaMSE = zeros(length(alphaVec),1);
wModelList = cell(1,length(alphaVec));
betaVec = cell(1,length(alphaVec));

for i=1:length(alphaVec)
    [coefficients,fitInfo] = ...
        lasso(XTrain.Variables,yTrain,'Alpha',alphaVec(i),'CV',cvFold,'MCReps',MCReps,'Lambda',lambdaVec);
    wModelList{i} = coefficients(:,fitInfo.IndexMinMSE);
    betaVec{i} = fitInfo.Intercept(fitInfo.IndexMinMSE);       
    indexMinMSE = fitInfo.IndexMinMSE;
    rmseList(i) = sqrt(fitInfo.MSE(indexMinMSE));    
    minLambdaMSE(i) = fitInfo.LambdaMinMSE;    
end

Чтобы сделать подогнанную модель надежной, используйте данные проверки, чтобы выбрать конечные значения $$$$$$$$гиперпараметров α и λ. Для этого выбирают коэффициенты, соответствующие четырем самым низким значениям RMSE, вычисленным с использованием обучающих данных.

numVal = 4;
[out,idx] = sort(rmseList);
val = out(1:numVal);
index = idx(1:numVal);

alpha = alphaVec(index);
lambda = minLambdaMSE(index);
wModel = wModelList(index);
beta = betaVec(index);

Для каждого набора коэффициентов вычисляют предсказанные значения и RMSE на данных проверки.

[XVal,yVal] = helperGetFeatures(valData);
rmseValModel = zeros(numVal);
for valList = 1:numVal   
   yPredVal = (XVal.Variables*wModel{valList} + beta{valList});
   rmseValModel(valList) = sqrt(mean((yPredVal-yVal).^2));
end

Выберите модель с наименьшим значением RMSE при использовании данных проверки. Модель, связанная с наименьшим значением RMSE для обучающих данных, не является моделью, соответствующей наименьшему значению RMSE для данных проверки. Таким образом, использование данных проверки повышает надежность модели.

[rmseMinVal,idx] = min(rmseValModel);
wModelFinal = wModel{idx};
betaFinal = beta{idx};

Оценка эффективности обучаемой модели

Данные теста содержат измерения, соответствующие 40 ячейкам. Извлеките элементы и соответствующие ответы из testData. Используйте обученную модель для прогнозирования оставшегося срока службы каждой ячейки в testData.

[XTest,yTest] = helperGetFeatures(testData);
yPredTest = (XTest.Variables*wModelFinal + betaFinal);

Визуализация прогнозируемого и фактического графика жизненного цикла для данных теста.

figure;
scatter(yTest,yPredTest)
hold on;
refline(1, 0);
title('Predicted vs Actual Cycle Life')
ylabel('Predicted cycle life');
xlabel('Actual cycle life');

В идеале все точки на приведенном выше графике должны быть близки к диагональной линии. Из вышеприведенного графика видно, что обученная модель хорошо работает, когда оставшийся срок службы цикла составляет от 500 до 1200 циклов. После 1200 циклов производительность модели ухудшается. Модель последовательно недооценивает оставшийся цикл жизни в этом регионе. Это объясняется главным образом тем, что наборы данных тестирования и проверки содержат значительно больше ячеек с общим сроком службы около 1000 циклов.

Вычислите RMSE прогнозируемого оставшегося срока службы цикла.

errTest = (yPredTest-yTest);
rmseTestModel = sqrt(mean(errTest.^2))

rmseTestModel = 211.6148

Другой метрикой, рассматриваемой для оценки эффективности, является средняя процентная ошибка, определенная [1] как

n = numel(yTest);
nr = abs(yTest - yPredTest);
errVal = (1/n)*sum(nr./yTest)*100

errVal = 9.9817

Заключение

Этот пример показал, как использовать модель линейной регрессии с эластичной регуляризацией сетки для прогнозирования срока службы цикла батареи на основе измерений только из первых 100 циклов. Пользовательские особенности были извлечены из необработанных данных измерения, и модель линейной регрессии с эластичной регуляризацией сетки была установлена с использованием обучающих данных. Гиперпараметры были затем выбраны с использованием набора данных проверки. Эта модель использовалась на тестовых данных для оценки эффективности. Используя измерения только для первых 100 циклов, RMSE оставшегося цикла прогнозирования срока службы клеток в наборе тестовых данных составляет 211,6 и получается средняя процентная погрешность 9,98%.

Ссылки

[1] Северсон, K.A., Attia, P.M., Jin, N. et al. «Прогнозирование срока службы аккумулятора на основе данных до снижения емкости». Nat Energy 4, 383-391 (2019). https://doi.org/10.1038/s41560-019-0356-8

[2] https://data.matr.io/1/

Вспомогательные функции

function [xTable, y] = helperGetFeatures(batch)
% HELPERGETFEATURES function accepts the raw measurement data and computes
% the following feature:
%
% Q_{100-10}(V) variance [DeltaQ_var]
% Q_{100-10}(V) minimum [DeltaQ_min]
% Slope of linear fit to the capacity fade curve cycles 2 to 100 [CapFadeCycle2Slope]
% Intercept of linear fit to capacity fade curve, cycles 2 to 100 [CapFadeCycle2Intercept]
% Discharge capacity at cycle 2 [Qd2]
% Average charge time over first 5 cycles [AvgChargeTime]
% Minimum Internal Resistance, cycles 2 to 100 [MinIR]
% Internal Resistance difference between cycles 2 and 100 [IRDiff2And100]

N = size(batch,2); % Number of batteries

% Preallocation of variables
y = zeros(N, 1);
DeltaQ_var = zeros(N,1);
DeltaQ_min = zeros(N,1);
CapFadeCycle2Slope = zeros(N,1);
CapFadeCycle2Intercept = zeros(N,1);
Qd2 = zeros(N,1);
AvgChargeTime = zeros(N,1);
IntegralTemp = zeros(N,1);
MinIR = zeros(N,1);
IRDiff2And100 = zeros(N,1);

for i = 1:N
    % cycle life
    y(i,1) = size(batch(i).cycles,2) + 1;
    
    % Identify cycles with time gaps
    numCycles = size(batch(i).cycles, 2);
    timeGapCycleIdx = [];
    jBadCycle = 0;
    for jCycle = 2: numCycles
        dt = diff(batch(i).cycles(jCycle).t);
        if max(dt) > 5*mean(dt)
            jBadCycle = jBadCycle + 1;
            timeGapCycleIdx(jBadCycle) = jCycle;
        end
    end
    
    % Remove cycles with time gaps
    batch(i).cycles(timeGapCycleIdx) = [];
    batch(i).summary.QDischarge(timeGapCycleIdx) = [];
    batch(i).summary.IR(timeGapCycleIdx) = [];
    batch(i).summary.chargetime(timeGapCycleIdx) = [];
    
    % compute Q_100_10 stats
    DeltaQ = batch(i).cycles(100).Qdlin - batch(i).cycles(10).Qdlin;
    DeltaQ_var(i) = log10(abs(var(DeltaQ)));
    DeltaQ_min(i) = log10(abs(min(DeltaQ)));
    
    % Slope and intercept of linear fit for capacity fade curve from cycle
    % 2 to cycle 100
    coeff2 = polyfit(batch(i).summary.cycle(2:100), batch(i).summary.QDischarge(2:100),1);
    CapFadeCycle2Slope(i) = coeff2(1);
    CapFadeCycle2Intercept(i) = coeff2(2);
    
    %  Discharge capacity at cycle 2
    Qd2(i) = batch(i).summary.QDischarge(2);
    
    % Avg charge time, first 5 cycles (2 to 6)
    AvgChargeTime(i) = mean(batch(i).summary.chargetime(2:6));
    
    % Integral of temperature from cycles 2 to 100
    tempIntT = 0;
    for jCycle = 2:100
        tempIntT = tempIntT + trapz(batch(i).cycles(jCycle).t, batch(i).cycles(jCycle).T);
    end
    IntegralTemp(i) = tempIntT;
    
    % Minimum internal resistance, cycles 2 to 100
    temp = batch(i).summary.IR(2:100);
    MinIR(i) = min(temp(temp~=0));
    IRDiff2And100(i) = batch(i).summary.IR(100) - batch(i).summary.IR(2);
end

xTable = table(DeltaQ_var, DeltaQ_min, ...
    CapFadeCycle2Slope, CapFadeCycle2Intercept,...
    Qd2, AvgChargeTime, MinIR, IRDiff2And100);

end