Жизненное предсказание цикла батареи из начальных данных об операции

Скрипт Open Live Script

В этом примере показано, как предсказать остающуюся жизнь цикла быстрого заряженного литий-ионного аккумулятора с помощью линейной регрессии, алгоритма машинного обучения с учителем. Жизненное предсказание цикла литий-ионного аккумулятора с помощью основанного на физике подхода моделирования является очень комплексным из-за различных условий работы и значительной изменчивости устройства даже с батареями от того же производителя. Для этого сценария основанные на машинном обучении подходы обеспечивают многообещающие результаты, когда достаточные тестовые данные доступны. Точное жизненное предсказание цикла батареи на ранних стадиях времени работы аккумулятора допускало бы быструю валидацию новых производственных процессов. Это также позволяет конечным пользователям идентифицировать ухудшенную эффективность с достаточной задержкой, чтобы заменить неисправные батареи. С этой целью только первые 100 основанных на циклах функций рассматриваются для предсказания остающейся жизни цикла, и ошибка предсказания, как показывают, составляет приблизительно 10% [1].

Набор данных

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

Описательные данные: штрихкод Ячейки, заряжая политику, жизнь цикла
Сводные данные на цикл: номер Цикла, способность выброса, внутреннее сопротивление, время зарядки
Данные собраны в цикле: Время, температура, линейно интерполировало способность выброса, линейно интерполировал напряжение

Загрузите данные MathWorks supportfiles сайт (это - большой набор данных, ~1.7GB).

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

Извлечение признаков

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

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]. В частности, изменение в напряжении выброса изгибается между циклами, $Δ Q (V)$ , является очень эффективным при диагнозе ухудшения. Поэтому выберите различие в способности выброса в зависимости от напряжения между 100-м и 10-м циклом, $Δ 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)');

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

$Δ Q_{100 - 10} (V)$ регистрируйте отклонение [DeltaQ_var]
$Δ Q_{100 - 10} (V)$ регистрируйте минимум [DeltaQ_min]
Наклон линейной подгонки к способности исчезает кривая, циклы 2 - 100 [CapFadeCycle2Slope]
Точка пересечения линейной подгонки к способности исчезает кривая, циклы 2 - 100 [CapFadeCycle2Intercept]
Разрядите способность в цикле 2 [Qd2]
Среднее время зарядки по первым 5 циклам [AvgChargeTime]
Минимальное Внутреннее Сопротивление, циклы 2 - 100 [MinIR]
Внутреннее различие в Сопротивлении между циклами 2 и 100 [IRDiff2And100]

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

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

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

$y_{i} = w^{T} x_{i} + β$

где $y_{i}$ количество циклов для ячейки $i$ , $x_{i}$ isa $p$ - размерный характеристический вектор для ячейки $i$ и $w$ $p$ - размерный вектор коэффициента модели и $β$ скалярная точка пересечения. Линейная модель упорядочена с помощью эластичной сети, чтобы обратиться к высоким корреляциям между функциями. Для $α$ строго между 0 и 1, и неотрицательный $λ$ , коэффициент регуляризации, эластичная сеть решает задачу:

$w_{}^{ˆ} = \min_{w} ‖ y - X w - β ‖_{2}^{2} + λ P_{α} (w)$

где $P_{α} (w) = \frac{(1 - α)}{2} ‖ w ‖_{2}^{2} + α ‖ w ‖_{1}$ . Эластичная сетевая регрессия гребня подходов, когда $α$ близко к 0, и это совпадает с регуляризацией лассо когда $α = 1$ . Обучающие данные используются, чтобы выбрать гиперпараметры $α$ и $λ$ , и определите значения коэффициентов, $w$ . lasso функция возвращает адаптированные коэффициенты регрессии наименьших квадратов и информацию о припадке модели, как выведено для каждого $α$ и $λ$ значение. lasso функция принимает вектор из значений для $λ$ параметр. Поэтому для каждого значения $α$ , коэффициенты модели $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 и 1 200 циклами. Вне 1 200 циклов ухудшается производительность модели. Модель последовательно недооценивает остающуюся жизнь цикла в этой области. Это, в основном, потому что тест и наборы данных валидации содержат значительно больше ячеек с общей жизнью приблизительно 1 000 циклов.

Вычислите RMSE предсказанной остающейся жизни цикла.

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

rmseTestModel = 211.6148

Другая метрика, рассмотренная для оценки результатов деятельности, является средней погрешностью в процентах, заданной [1] как

$% e r r = \frac{1}{n} Σ_{i = 1}^{n} \frac{| y_{i} - {y_{}^{ˆ}}_{i} |}{y_{i}} \times 100$

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, пополудни, Чжин, N. и др. "Управляемое данными предсказание жизни цикла батареи перед полным ухудшением". Туземная энергия 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

Смотрите также

lasso

Документация