Набор данных

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

Загрузите данные с сайта supportfiles 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-м циклами, $$\Delta Q_{100-10}(V)$$.

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)');

Вычисление сводной статистики с помощью $$\Delta Q_{100-10}(V)$$ кривые каждой камеры в качестве индикатора состояния. Важно иметь в виду, что эта сводная статистика не должна иметь четкого физического смысла, они просто должны показать прогнозирующие возможности. Для примера отклонения $$\Delta Q_{100-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]:

The 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 

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

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

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

где $$P_{\alpha }(w)=\frac{(1-\alpha )}{2}\Vert w{\Vert }_2^2+\alpha \Vert w{\Vert }_1$$. Упругая сеть приближается к регрессии хребта, когда $$\alpha$$ близок к 0, и это то же самое, что и регуляризация лассо, когда $$\alpha =1$$. Обучающие данные используются для выбора гиперпараметров $$\alpha$$ и $$\lambda$$, и определить значения коэффициентов, $$w$$. The lasso функция возвращает установленные коэффициенты регрессии методом наименьших квадратов и информацию о подгонке модели в качестве вывода для каждого $$\alpha$$ и $$\lambda$$ значение. The lasso функция принимает вектор значений для $$\lambda$$ параметр. Поэтому для каждого значения $$\alpha$$, коэффициенты модели $$w$$ и $$\lambda$$ вычисляют, получая наименьшее значение 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

Чтобы сделать подобранную модель устойчивой, используйте данные валидации, чтобы выбрать окончательные значения $$\alpha$$ и $$\lambda$$ гиперпараметры. С этой целью выберите коэффициенты, соответствующие четырем самым низким значениям 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] Severson, K.A., Attia, PM, 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