Feat: implemented battery profile simulation via BatProfile class

.gitignore

@@ -0,0 +1,2 @@
+python/__pycache__/
+python/functions/__pycache__/

README.md

@@ -1,8 +1,45 @@
+### Profile Simulation with BLAST
+
+This project is forked from ![BLAST](https://github.com/NREL/BLAST-Lite) with
+an implementation of a custom profile simulator. Refer to
+`simulation_with_custom_profile.ipynb` for the usage of the custom profile. The
+custom profile is implemented in `battery_profile.py`.
+
 ### BLAST-Lite
 
-Battery Lifetime Analysis and Simulation Toolsuite in the python programming language.
-Provides a library of battery lifetime and degradation models for various commercial lithium-ion batteries.
+Battery Lifetime Analysis and Simulation Toolsuite (BLAST) provides a library of battery lifetime and degradation models for various commercial lithium-ion batteries from recent years. Degradation models are indentified from publically available lab-based aging data using NREL's battery life model identification toolkit. The battery life models predicted the expected lifetime of batteries used in mobile or stationary applications as functions of their temperature and use (state-of-charge, depth-of-discharge, and charge/discharge rates). Model implementation is in both Python and MATLAB programming languages. The MATLAB code also provides example applications (stationary storage and EV), climate data, and simple thermal management options. For more information on battery health diagnostics, prediction, and optimization, see [NREL's Battery Lifespan](https://www.nrel.gov/transportation/battery-lifespan.html) webpage.
 
 ![Example battery life predictions](example_battery_life.png)
 
-NREL SWR-22-69
+## Caveats
+These battery models predict 'expected life', that is, battery life under nominal conditions. Many types of battery failure will not be predicted by these models:
+- Overcharge or overdischarge
+- Impact of physical damage, vibration, or humidity
+- Operating outside of manufacturer performance and environmental limits, such as voltage, temperature, and charge/discharge rate limits
+- Pack performance loss due to cell-to-cell inbalance
+
+Aging models are generally trained on a limited amount of data, that is, there is not enough information to estimate cell-to-cell variability in degradation rates.
+Battery 'warranty life' is generally much more conservative than 'expected life'.
+
+## Citations:
+ - Sony Murata LFP-Gr battery aging data and model
+     - [Calendar aging data source](https://doi.org/10.1016/j.est.2018.01.019)
+     - [Cycle aging data source](https://doi.org/10.1016/j.jpowsour.2019.227666)
+     - [Model identification source](https://doi.org/10.1149/1945-7111/ac86a8)
+ - Nissan Leaf LMO-Gr (2nd life cells) battery aging data
+     - [Calendar aging data source] (https://doi.org/10.1109/EEEIC/ICPSEUROPE54979.2022.9854784)
+     - [Cycle aging data source] (https://doi.org/10.1016/j.est.2020.101695)
+ - Panasonic NCA-Gr battery aging data
+     - [Calendar aging data source](https://dx.doi.org/10.1149/2.0411609jes)
+     - [Cycle aging data source](https://doi.org/10.1149/1945-7111/abae37)
+ - [Commercial NMC-LTO battery aging data source](https://doi.org/10.1016/j.jpowsour.2020.228566)
+ - [Kokam NMC111-Gr battery aging data source](https://ieeexplore.ieee.org/iel7/7951530/7962914/07963578.pdf)
+ - [Sanyo NMC111-Gr battery aging model source](http://dx.doi.org/10.1016/j.jpowsour.2014.02.012)
+ - [LG MJ1 NMC811-GrSi battery aging data source](https://everlasting-project.eu/wp-content/uploads/2020/03/EVERLASTING_D2.3_final_20200228.pdf)
+ - Stationary storage battery use profiles are [provided open-source](https://dataserv.ub.tum.de/index.php/s/m1510254) by [Kucevic et al](https://www.sciencedirect.com/science/article/pii/S2352152X19309016)
+  - Electric vehicle battery use profiles were generated using NREL's [FASTSim tool](https://www.nrel.gov/transportation/fastsim.html).
+
+## Authors
+Paul Gasper, Kandler Smith
+
+NREL SWR-22-69

license.txt

@@ -0,0 +1,27 @@
+Copyright (c) 2023, Alliance for Sustainable Energy, LLC
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

matlab/Application profiles/sources.txt

@@ -0,0 +1,11 @@
+Electric vehicle battery application profiles:
+EV battery profiles are 1 week of simulated battery use.
+'personal_ev_largebatt.csv': generated using NREL's FASTSim tool. Large EV battery (300+ miles range) only requires 1 partial charge and 1 full charge per week, assuming daily driving and a long weekend trip.
+'personal_ev_smallbatt.csv': generated using NREL's FASTSim tool. Smaller EV battery has less range, so it is discharged more deeply and requires more charging events than a large battery.
+
+Stationary storage battery application profiles:
+All stationary storage application profiles are simulated battery load from real-world power demands. FCR and PV_BESS are 1 month, peak shaving is 1 year.
+Profiles were provided open-source by Kucevic et al, https://doi.org/10.1016/j.est.2019.101077.
+'FCR_1PE_LFP.mat': Frequency containment reserve application - constant low depth-of-discharge use.
+'PS_3Cluster_1.mat': Peak-shaving application - only shaves high demand peaks in warmer months.
+'PVBESS_FeedInDamp.mat': Residential solar + battery - battery charges and discharges ~1 time every single day, with varying depth of discharge.

matlab/Functions/func_LifeMdl_LoadParameters.m

@@ -0,0 +1,79 @@
+function lifeMdl = func_LifeMdl_LoadParameters(battery_chemistry)
+switch battery_chemistry
+    case 'NMC|Gr'
+        lifeMdl.model = 'NMC|Gr';
+        % battery life model parameters
+        lifeMdl.ref          = 'NREL/K. Smith et al. American Control Conference, 2017';
+        lifeMdl.description  = 'Model fit to 13 aging conditions for Kokam 75Ah cell. Data collected by NREL';
+        lifeMdl.chemistry    = 'Graphite negative electrode, Ni-Mn-Co (NMC) positive electrode';
+        lifeMdl.names_states = {...
+            'dq_LLI_t:   Delta relative capacity-Li loss due to SEI growth, time dependent';...
+            'dq_LLI_EFC: Delta relative capacity-Li loss due to SEI growth, charge-throughput dependent';...
+            'dq_LAM:     Delta relative capacity-Electrode damage and lithium plating, charge-throughput dependent';...
+            'dr_LLI_t:   Delta relative resistance-SEI growth, time dependent';...
+            'dr_LLI_EFC: Delta relative resistance-SEI growth, charge-throughput dependent';...
+            'dr_LAM:     Delta relative resistance-Electrode damage and lithium plating, charge-throughput dependent'};
+        lifeMdl.names_outputs= {...
+            'q:                   relative discharge capacity at C/5 rate (5 hour discharge)';...
+            'qLoss_LLI:           relative capacity loss due to SEI growth';...
+            'qLoss_LLI_t:         relative capacity loss due to SEI growth, time dependent';...
+            'qLoss_LLI_t_LLI_EFC: relative capacity loss due to SEI growth, charge-throughput dependent';...
+            'qLoss_LLI_t_LAM:     relative capacity loss due to electrode damage and lithium plating, charge-throughput dependent';...
+            'r:                   relative DC pulse resistance';...
+            'rGain_LLI:           relative resistance gain due to SEI growth';...
+            'rGain_LLI_t:         relative resistance gain due to SEI growth, time dependent';...
+            'rGain_LLI_EFC:       relative resistance gain due to SEI growth, charge-throughput dependent';...
+            'rGain_LAM:           relative resistance gain due to electrode damage and lithium plating, charge-throughput dependent'};
+
+        lifeMdl.n_states = 6;      % number of states in model that are integrated with time or cycles
+        lifeMdl.n_outputs= 10;      % number of outputs in model
+
+        % Model parameters
+        lifeMdl.q1_0 = 2.66e7;
+        lifeMdl.q1_1 = -17.8;
+        lifeMdl.q1_2 = -5.21;
+        lifeMdl.q2 = 0.357;
+        lifeMdl.q3_0 = 3.80e3;
+        lifeMdl.q3_1 = -18.4;
+        lifeMdl.q3_2 = 1.04;
+        lifeMdl.q4 = 0.778;
+        lifeMdl.q5_0 = 1e4;
+        lifeMdl.q5_1 = 153;
+        lifeMdl.p_LAM = 10;
+        lifeMdl.r1 = 0.0570;
+        lifeMdl.r2 = 1.25;
+        lifeMdl.r3 = 4.87;
+        lifeMdl.r4 = 0.712;
+        lifeMdl.r5 = -0.08;
+        lifeMdl.r6 = 1.09;
+    case 'LFP|Gr'
+        % Life model information
+        lifeMdl.model = 'LFP|Gr';
+        lifeMdl.names_states = {...
+            'dqCal: time dependent capacity loss',...
+            'dqCyc: energy-throughput dependent capacity loss',...
+            'drCal: time dependent resistance growth',...
+            'drCyc: energy-throughput dependent resistance growth'};
+        lifeMdl.names_outputs= {...
+            'q: Relative capacity',...
+            'qLossCal: Time dependent capacity loss',...
+            'qLossCyc: Energy-throughput dependent capacity loss',...
+            'r: Relative resistance',...
+            'rGainCal: Time dependent resistance growth',...
+            'rGainCyc: Energy-throughput dependent resistance growth',...
+            };
+        lifeMdl.n_states = 4;      % number of states in model that are integrated with time or cycles
+        lifeMdl.n_outputs= 6;      % number of outputs in model
+
+        % Instantiate parameter table
+        pVars = {'p1','p2','p3','pcal','p4','p5','p6','pcyc'};
+        p = [798.352102069967,-3781.20591728016,-2.04881543271387,0.494361078068684,0.000170691091574341,0.000476715369722458,0.00217277383391531,0.439274000189976];
+        p = array2table(p, 'VariableNames', pVars);
+        lifeMdl.p = p;
+
+        %         pvars_rdc = {'r1','r2','r3','rcal','r4','r5','r6','rcyc'};
+        p_rdc = [537306.726059759,-6747.41305124106,-1.66987638954611,0.994115565719101,4.64136140815756e-08,-5.79557082332511e-11,7.28339782547118e-06,0.931662358013465];
+        p_rdc = array2table(p_rdc, 'VariableNames', pVars);
+        lifeMdl.p_rdc = p_rdc;
+end
+end

matlab/Functions/func_LifeMdl_StressStatistics.m

@@ -0,0 +1,70 @@
+function cycle = func_LifeMdl_StressStatistics(cycle)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Adds to "cycle" structure additional aging stress statistics.
+%  Scalars:
+%   cycle.Crate_rms,Crate_disrms,Crate_chgrms
+%  Vectors: (using rainflow algorithm)
+%   cycle.dsoc_i,ncyc_i,Crate_dis_i,Crate_chg_i,TavgK_i,socavg_i...
+%
+% Copyright Kandler Smith & Ying Shi, NREL 6/2013
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%{
+Input timeseries:
+    cycle.t:     time in days
+    cycle.soc:   soc in normal units, 0 to 1
+    cycle.TdegC: temperature in celsius
+Outputs:
+    cycle.dt:    Time difference in days
+    cycle.dEFC:  Equvialent full cycles during this period
+    cycle.TdegK: Time weighted average of temperature
+    cycle.soc:   Time weighted average of soc
+    cycle.Ua:    Time weighted average of Ua
+    cycle.dod:   Difference b/w max and min SOC
+    cycle.Crate: Time weighted average of Crate in non-resting periods
+%}
+% dt
+dt = cycle.t(end) - cycle.t(1);
+
+% TdegK
+TdegK = trapz(cycle.t, cycle.TdegC + 273.15)/dt;
+
+% soc
+soc = trapz(cycle.t, cycle.soc)/dt;
+
+% Ua:
+x_a_eq = @(SOC) 8.5e-3 + SOC.*(7.8e-1 - 8.5e-3);
+Ua_eq = @(x_a) 0.6379 + 0.5416.*exp(-305.5309.*x_a) + 0.044.*tanh(-1.*(x_a-0.1958)./0.1088) - 0.1978.*tanh((x_a-1.0571)./0.0854) - 0.6875.*tanh((x_a+0.0117)./0.0529) - 0.0175.*tanh((x_a-0.5692)./0.0875);
+Ua = Ua_eq(x_a_eq(soc));
+
+% DOD
+dod = max(cycle.soc) - min(cycle.soc);
+
+% Crate
+Crate = diff(cycle.soc)./diff(cycle.t .* 24);
+% Crate is the average of the non-zero rates (Absolute value)
+% Any rate smaller than a 1000 hour chg/dischg rate is
+% considered to be negligible.
+Crate = abs(Crate);
+maskNonResting = Crate > 1e-3;
+if any(maskNonResting)
+    Crate = Crate(maskNonResting);
+    difft = diff(cycle.t);
+    difft = difft(maskNonResting);
+    tNonResting = cumsum(difft);
+    Crate = trapz(tNonResting, Crate)/(tNonResting(end)-tNonResting(1));
+else
+    Crate = 0;
+end
+% dEFC
+dEFC = sum(abs(diff(cycle.soc)))/2;
+
+% output the extracted stress statistics
+cycle = struct();
+cycle.dt =      dt;
+cycle.dEFC =    dEFC;
+cycle.TdegK =   TdegK;
+cycle.soc =     soc;
+cycle.Ua =      Ua;
+cycle.dod =     dod;
+cycle.Crate =   Crate;
+end

matlab/Functions/func_LifeMdl_UpdateOutput.m

@@ -0,0 +1,50 @@
+function yy = func_LifeMdl_UpdateOutput(lifeMdl, xx)
+% Calculates capacity fade & resistance growth at a given state-of-life
+%
+%   - Author: Kandler Smith, NREL, Copyright 2010
+%   - Based on paper by EVS-24 paper by Smith,Markel,Pesaran, 2009
+%   - Modified for new NCA model, KAS 6/2010
+%   - Added FeP model based on A123 2.3Ah 26650 cell (KAS, 4/28/15)
+%   - Adapted FeP model for Kokam 75Ah NMC cell (KAS,8/1/16)
+%
+% Inputs:
+%   perfMdl - structure containing battery performance model parameters
+%
+%   lifeMdl - structure containing battery life model parameters
+%           * see func_LoadLifeParameters.m for definition of life-parameters
+%
+%   xx      - vector of life model states at present state-of-life
+%
+% Outputs:
+%   yy      - vector of life model capacity fade and resistance growth
+
+switch lifeMdl.model
+    case 'NMC|Gr'
+        q_LLI = 1 - xx(1) - xx(2);
+        q_LLI_t = 1 - xx(1);
+        q_LLI_EFC = 1 - xx(2);
+        q_LAM = 1.01 - xx(3);
+        q = min([q_LLI, q_LAM]);
+
+        % Resistance
+        r_LLI = 1 + xx(4) + xx(5);
+        r_LLI_t = 1 + xx(4);
+        r_LLI_EFC = 1 + xx(5);
+        r_LAM = lifeMdl.r5 + lifeMdl.r6 * (1 / q_LAM);
+        r = max([r_LLI, r_LAM]);
+
+        % Assemble output
+        yy = [q, q_LLI, q_LLI_t, q_LLI_EFC, q_LAM, r, r_LLI, r_LLI_t, r_LLI_EFC, r_LAM];
+    case 'LFP|Gr'
+        qLossCal = xx(1);
+        % Scale cycling loss to assume cells aren't so long lived as the
+        % Sony Murata cells are (10-15 thousand EFCs).
+        scalingfactor = 4;
+        qLossCyc = scalingfactor*xx(2);
+        q = 1 - qLossCal - qLossCyc;
+        rGainCal = xx(3);
+        rGainCyc = scalingfactor*xx(4);
+        r = 1 + rGainCal + rGainCyc;
+        yy = [q qLossCal qLossCyc r rGainCal rGainCyc];
+end
+end

matlab/Functions/func_LifeMdl_UpdateStates.m

@@ -0,0 +1,133 @@
+function xx = func_LifeMdl_UpdateStates(lifeMdl,cycle,xx0)
+% Calculates capacity fade & resistance growth rates & integrates rates
+%   forward in time to update states for a complex cycling profile
+%   at a given state-of-life.
+%
+%   - Author: Kandler Smith, NREL, Copyright 2010
+%   - Based on paper by EVS-24 paper by Smith,Markel,Pesaran, 2009
+%   - Modified for new NCA model, KAS 6/2010
+%   - Added FeP model based on A123 2.3Ah 26650 cell (KAS, 4/28/15)
+%   - Adapted FeP model for Kokam 75Ah NMC cell (KAS,8/1/16)
+%   - Adapted for Shell Battery Life Models (Paul Gasper, 2021)
+%
+% Inputs:
+%   perfMdl - structure containing battery performance model parameters
+%
+%   lifeMdl - structure containing battery life model parameters
+%           * see func_LoadLifeParameters.m for definition of life-parameters
+%
+%   cycle   - structure containing duty-cycle & temperature variables, e.g.
+%           cycle.tsec  - time stamp array
+%           cycle.soc   - state-of-charge array corresponding to time stamps (0=empty, 1=full charge)
+%           cycle.TdegC - temperature array corresponding to time stamps (or use a
+%                         single value to represent constant temperature) (degrees Celcius)
+%           cycle.dsoc_i- array of delta-state-of-charge swings encountered by
+%                         the battery during cycling calulated by Rainflow algorithm
+%           cycle.ncycle_i-array of number of cycles corresponding to each
+%                         delta-state-of-charge calculated by Rainflow algorithm
+%
+%   dtdays_life - timestep between previous and present battery
+%            state-of-life (days)
+%
+%   xx0  - vector of life model states at previous state-of-life
+%
+% Outputs:
+%   xx    - vector of life model states integrated to next state-of-life
+%
+%   dxxdt - vector of life model states time rate of change
+
+% Rate equations are derived from trajectory equations in this manner:
+%   Traj. eq: y = f(x)
+%   Invert to solve for x: x_hat = g(y)
+%   Rate eq: dy/dx = df(x_hat)/dx
+%   The change of the state in each timestep is then dy/dx * dx (x = time
+%       or energy throughput or whatever)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Normalized stressor values
+cycle.TdegKN = cycle.TdegK ./ (273.15 + 35);
+cycle.TdegC = cycle.TdegK - 273.15;
+cycle.UaN = cycle.Ua ./ 0.123;
+
+switch lifeMdl.model
+    case 'NMC|Gr'
+        % Calculate degradation rates
+        q1 = lifeMdl.q1_0 .* exp(lifeMdl.q1_1 .* (1 ./ cycle.TdegKN)) .* exp(lifeMdl.q1_2 .* (cycle.UaN ./ cycle.TdegKN));
+        q3 = lifeMdl.q3_0 .* exp(lifeMdl.q3_1 .* (1 ./ cycle.TdegKN)) .* exp(lifeMdl.q3_2 .* exp(cycle.dod.^2));
+        q5 = lifeMdl.q5_0 + lifeMdl.q5_1 .* (cycle.TdegC - 55) .* cycle.dod;
+
+        % Calculate incremental state changes
+        dq_LLI_t = dynamicPowerLaw(xx0(1), cycle.dt, 2*q1, lifeMdl.q2);
+        dq_LLI_EFC = dynamicPowerLaw(xx0(2), cycle.dEFC, q3, lifeMdl.q4);
+        dq_LAM = dynamicSigmoid(xx0(3), cycle.dEFC, 1, 1/q5, lifeMdl.p_LAM);
+        dr_LLI_t = dynamicPowerLaw(xx0(4), cycle.dt, lifeMdl.r1 * q1, lifeMdl.r2);
+        dr_LLI_EFC = dynamicPowerLaw(xx0(5), cycle.dEFC, lifeMdl.r3 * q3, lifeMdl.r4);
+
+        % Pack up
+        dxx = [dq_LLI_t; dq_LLI_EFC; dq_LAM; dr_LLI_t; dr_LLI_EFC; xx0(6)];
+        xx = xx0 + dxx;
+
+    case 'LFP|Gr'
+        % capacity
+        p = lifeMdl.p;
+        
+        % Degradation rate equations
+        kcal = @(S,p) abs(p.p1) * exp(p.p2 / S.TdegK) * exp(p.p3 * S.Ua);
+        kcyc = @(S,p) (abs(p.p4) + abs(p.p5) * S.dod + abs(p.p6) * S.Crate);
+        
+        % Calendar loss
+        dqLossCal = dynamicPowerLaw(xx0(1), cycle.dt, kcal(cycle,p), p.pcal);
+        % Cycling loss
+        dqLossCyc = dynamicPowerLaw(xx0(2), cycle.dEFC, kcyc(cycle,p), p.pcyc);
+        
+        % resistance
+        p = lifeMdl.p_rdc;
+        drGainCal = dynamicPowerLaw(xx0(3), cycle.dt, kcal(cycle,p), p.pcal);
+        drGainCyc = dynamicPowerLaw(xx0(4), cycle.dEFC, kcyc(cycle,p), p.pcyc);
+        
+        % Pack up & accumulate state vector forward in time
+        dxx = [dqLossCal; dqLossCyc; drGainCal; drGainCyc];
+        xx = xx0 + dxx;
+end
+
+    function dy = dynamicPowerLaw(y0, dx, k, p)
+        % DYNAMICPOWERLAW calculates the change of state for states modeled by a
+        % power law equation, y = k*x^p. The output is instantaneous slope of dy
+        % with respect to dx.
+        if y0 == 0
+            if dx == 0
+                dydx = 0;
+            else
+                y0   = k * dx^p;
+                dydx = y0 / dx;
+            end
+        else
+            if dx == 0
+                dydx = 0;
+            else
+                dydx = k * p * (y0 / k)^((p-1) / p);
+            end
+        end
+        dy = dydx * dx;
+    end
+
+    function dy = dynamicSigmoid(y0, dx, y_inf, k, p)
+        if y0 == 0
+            if dx == 0
+                dydx = 0;
+            else
+                dy = 2 * y_inf * (1/2 - 1 / (1 + exp((k * dx) ^ p)));
+                dydx = dy / dx;
+            end
+        else
+            if dx == 0
+                dydx = 0;
+            else
+                x_inv = (1 / k) * ((log(-(2 * y_inf/(y0-y_inf)) - 1)) ^ (1 / p) );
+                z = (k * x_inv) ^ p;
+                dydx = (2 * y_inf * p * exp(z) * z) / (x_inv * (exp(z) + 1) ^ 2);
+            end
+        end
+        dy = dydx * dx;
+    end
+end

matlab/Functions/func_ThrowOutBadData.m

@@ -0,0 +1,19 @@
+function cycle1 = func_ThrowOutBadData(cycle1)
+% gets rid of NaN and out-of-range data
+% "keep" is indices of good data
+% keep = find(  cycle1.I>-200    & cycle1.I<200 ...
+%             & cycle1.V>0       & cycle1.V<100 ...
+%             & cycle1.TdegC>-40 & cycle1.TdegC<100 ...
+%             & cycle1.soc>=0    & cycle1.soc<=1 ...
+%             & [1; diff(cycle1.tsec)]>0        );
+keep = find(  cycle1.TdegC>-40 & cycle1.TdegC<100 ...
+            & cycle1.soc>=0    & cycle1.soc<=1 ...
+            & [1; diff(cycle1.tsec)]>0        );
+  
+cycle1.I    = cycle1.I(keep);
+cycle1.soc  = cycle1.soc(keep);
+cycle1.TdegC= cycle1.TdegC(keep);
+cycle1.V    = cycle1.V(keep);
+cycle1.tsec = cycle1.tsec(keep);
+
+end

matlab/Functions/rainflow.m

@@ -0,0 +1,30 @@
+% RAINFLOW cycle counting.
+%   RAINFLOW counting function allows you to extract 
+%   cycle from random loading.
+%
+% SYNTAX
+%   rf = RAINFLOW(ext)
+%   rf = RAINFLOW(ext, dt)
+%   rf = RAINFLOW(ext, extt)
+%
+% OUTPUT
+%   rf - rainflow cycles: matrix 3xn or 5xn dependend on input,
+%     rf(1,:) Cycles amplitude,
+%     rf(2,:) Cycles mean value,
+%     rf(3,:) Number of cycles (0.5 or 1.0),
+%     rf(4,:) Begining time (when input includes dt or extt data),
+%     rf(5,:) Cycle period (when input includes dt or extt data),
+%
+% INPUT
+%   ext  - signal points, vector nx1, ONLY TURNING POINTS!,
+%   dt   - sampling time, positive number, when the turning points
+%          spaced equally,
+%   extt - signal time, vector nx1, exact time of occurrence of turning points.
+%
+%
+%   See also SIG2EXT, RFHIST, RFMATRIX, RFPDF3D.
+
+%   RAINFLOW
+%   Copyright (c) 1999-2002 by Adam Nieslony,
+%   MEX function.
+

matlab/Functions/sig2ext.m

@@ -0,0 +1,119 @@
+function [ext, exttime] = sig2ext(sig, dt, clsn)
+% SIG2EXT - search for local extrema in the time history (signal),
+%
+% function [ext, exttime] = sig2ext(sig, dt, clsn)
+%
+% SYNTAX
+%   sig2ext(sig)
+%   [ext]=sig2ext(sig)
+%   [ext,exttime]=sig2ext(sig)
+%   [ext,exttime]=sig2ext(sig, dt)
+%   [ext,exttime]=sig2ext(sig, dt, clsn)
+%
+% OUTPUT
+%   EXT     - found extrema (turning points of the min and max type)
+%             in time history SIG,
+%   EXTTIME - option, time of extremum occurrence counted from
+%             sampling time DT (in seconds) or time vector DT.
+%             If no sampling time present, DT = 1 is assumed.
+%
+% INPUT
+%   SIG     - required, time history of loading,
+%   DT      - option, descripion as above, scalar or vector of
+%             the same length as SIG,
+%   CLSN    - option, a number of classes of SIG (division is performed
+%             before searching of extrema), no CLSN means no division
+%             into classes.
+%
+% The function caused without an output draws a course graph with
+% the searched extrema.
+%
+
+% By Adam Nies³ony
+% Revised, 10-Nov-2009
+% Visit the MATLAB Central File Exchange for latest version
+
+error(nargchk(1,3,nargin))
+
+% Is the time analysed?
+TimeAnalize=(nargout==0)|(nargout==2);
+
+% Sprawdzam czy przyrost dt jest podany prawid³owo
+if nargin==1,
+    dt=1;
+else
+    dt=dt(:);
+end
+
+% Zamieniam dane sig na jedn¹ kolumnê
+sig=sig(:);
+
+% Dzielimy na klasy je¿eli jest podane CLSN
+if nargin==3,
+    if nargout==0,
+        oldsig=sig;
+    end
+    clsn=round(clsn);
+    smax=max(sig);
+    smin=min(sig);
+    sig=clsn*((sig-smin)./(smax-smin));
+    sig=fix(sig);
+    sig(sig==clsn)=clsn-1;
+    sig=(smax-smin)/(clsn-1)*sig+smin;
+end
+
+% Tworzê wektor binarny w gdzie 1 oznacza ekstremum lub równoœæ,
+% Uznajê ¿e pierwszy i ostatni punkt to ekstremum
+w1=diff(sig);
+w=logical([1;(w1(1:end-1).*w1(2:end))<=0;1]);
+ext=sig(w);
+if TimeAnalize,
+    if length(dt)==1,
+        exttime=(find(w==1)-1).*dt;
+    else
+        exttime=dt(w);
+    end
+end
+
+% Usuwam potrójne wartoœci
+w1=diff(ext);
+w=~logical([0; w1(1:end-1)==0 & w1(2:end)==0; 0]);
+ext=ext(w);
+if TimeAnalize,
+    exttime=exttime(w);
+end
+
+% Usuwam podwójne wartoœci i przesuwam czas na œrodek
+w=~logical([0; ext(1:end-1)==ext(2:end)]);
+ext=ext(w);
+if TimeAnalize,
+    w1=(exttime(2:end)-exttime(1:end-1))./2;
+    exttime=[exttime(1:end-1)+w1.*~w(2:end); exttime(end)];
+    exttime=exttime(w);
+end
+
+% Jeszcze raz sprawdzam ekstrema
+if length(ext)>2,  % warunek: w tym momencie mo¿e ju¿ byæ ma³o punktów
+    w1=diff(ext);
+    w=logical([1; w1(1:end-1).*w1(2:end)<0; 1]);
+    ext=ext(w);
+    if TimeAnalize,
+        exttime=exttime(w);
+    end
+end
+
+if nargout==0,
+    if length(dt)==1,
+        dt=(0:length(sig)-1).*dt;
+    end
+    if nargin==3,
+        plot(dt,oldsig,'b-',dt,sig,'g-',exttime,ext,'ro')
+        legend('signal','singal divided in classes','extrema')
+    else
+        plot(dt,sig,'b-',exttime,ext,'ro')
+        legend('signal','extrema')
+    end
+    xlabel('time')
+    ylabel('signal & extrema')
+    clear ext exttime
+end

matlab/README.txt

@@ -0,0 +1,5 @@
+Run the script 'simulate_battery_aging.m' in MATLAB r2020b (or later version) to simulate battery lifetime. Battery life has three inputs: battery chemistry, climate, and the application. The stress on the battery is due to both climate and application.
+Climate is picked from annual temperature data recorded for the 100 most populous U.S. regional airports.
+Application profiles are uploaded by the user into the 'Application profiles' folder. Some preloaded stationary storage and mobile application profiles are provided. Users can upload any arbitrary application profiles as a '.csv' file. Application profiles should have two columns: 'tsecs', the time in seconds, and 'soc', the state-of-charge of the battery. Ensure the start and end state-of-charge of the battery are within 1% for an accurate life simulation.
+Temperature management of the battery system can be turned on or off. The simplistic temperature management system assumes a maximum and minimum battery temperature when the battery is being used (charged or discharged). These maximum and minimum limits are defined by the user. Battery temperature during rest is assumed to be equal to ambient.
+Simulation results are plotted as figures and automatically written to an excel file, saved in the 'Simulation results' folder.

matlab/simulate_battery_aging.m

@@ -0,0 +1,387 @@
+clear; clc;
+path(path,'Functions')
+addpath(genpath('Aging scenarios'))
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 0) Select model and load parameters for battery life (degradation),
+%    performance (ocv relationships), and load aging test data
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+validModels = {...
+    'NMC|Gr',...
+    'LFP|Gr'...
+    };
+
+idxModels = listdlg('ListString', validModels,...
+    'PromptString', "Specify battery models.",...
+    'Name', 'Specify model',...
+    'ListSize', [300 300]);
+assert(~isempty(idxModels), "Must select a model.")
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 1) Select application profile to simulate
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+applicationProfiles =  dir('Application profiles/*.mat');
+driveCycles = dir('Application profiles\*.csv');
+applicationProfiles = {applicationProfiles.name, driveCycles.name};
+idx = listdlg('ListString', applicationProfiles,...
+    'PromptString', "Pick an application profile.",...
+    'ListSize', [300 300]);
+assert(~isempty(idx), "Must select an application profile.")
+
+% Ask for user options
+inputs = {'Ambient temperature in Celsius (-Inf to select a climate):', 'Simulation length in years:'};
+dims = [1 50];
+defaults = {'-Inf', '20'};
+out = inputdlg(inputs, 'Specify simulation options.', dims, defaults);
+TdegC = str2double(out{1});
+dtdays = 1;
+tYears = str2double(out{2});
+if TdegC == -Inf
+    load('Climate data\Top100MSA_hourlyClimateData.mat', 'TMY_data_locations','TMY_ambTemp_degC')
+    idxClimate = listdlg('ListString', TMY_data_locations',...
+        'PromptString', "Specify climate(s) for battery life simulation:",...
+        'ListSize', [300 300]);
+    TdegC = cell(length(idxClimate), 1);
+    for i = 1:length(TdegC)
+        TdegC{i} = TMY_ambTemp_degC(:, idxClimate(i));
+    end
+    climates = TMY_data_locations(idxClimate);
+end
+% Thermal management
+thermalManagement = questdlg('Specify thermal management strategy:',...
+    'Thermal management',...
+    'None', 'During charge/discharge', 'Always',...
+    'None');
+switch thermalManagement
+    case {'During charge/discharge'; 'Always'}
+        out = inputdlg(...
+            {'Minimum allowed temperature (C)', 'Maximum allowed temperature (C)'},...
+            'Specify temperature limits.', ...
+            [1 50],...
+            {'10', '40'});
+        thermalManagement_MinT = str2double(out{1});
+        thermalManagement_MaxT = str2double(out{2});
+        assert(thermalManagement_MinT < thermalManagement_MaxT, "Max temperature must be greater than min temperature.")
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 2) Loop through selected data
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% struct to store outputs:
+simulations = struct(); simIdx = 1;
+for scenarioIdx = idx
+    scenario = applicationProfiles{scenarioIdx};
+    fprintf("Application profile: %s\n", scenario)
+    disp("Loading application profile...")
+    if contains(scenario, '.mat')
+        switch scenario
+            case 'FCR_1PE_LFP.mat'
+                load(['Application profiles/', scenario], 'reference_soc');
+                soc = reference_soc./100; clearvars reference_soc
+                tsec = [0:1:length(soc)-1]';
+                cycle = table(tsec, soc);
+                clearvars t tsec soc
+                % Downsample the data from 1 second resolution
+                cycle = cycle(1:30:height(cycle), :);
+            case 'PS_3Cluster_1.mat'
+                load(['Application profiles/', scenario], 'reference_soc');
+                soc = reference_soc./100; clearvars reference_soc
+                tsec = [0:1:length(soc)-1]';
+                cycle = table(tsec, soc);
+                clearvars t tsec soc
+                % Downsample the data from 1 second resolution
+                cycle = cycle(1:30:height(cycle), :);
+            case 'PVBESS_FeedInDamp.mat'
+                load(['Application profiles/', scenario], 'reference_soc');
+                soc = reference_soc./100; clearvars reference_soc
+                tsec = [0:1:length(soc)-1]';
+                cycle = table(tsec, soc);
+                clearvars t tsec soc
+                % Downsample the data from 1 second resolution
+                cycle = cycle(1:30:height(cycle), :);
+                % For the BESS system, apply a 5% minimum SOC
+                cycle.soc = cycle.soc.*0.95 + 0.05;
+        end
+    elseif contains(scenario, '.csv')
+        % FastSim input profile
+        opts = delimitedTextImportOptions("NumVariables", 2);
+        % Specify range and delimiter
+        opts.DataLines = [2, Inf];
+        opts.Delimiter = ",";
+        % Specify column names and types
+        opts.VariableNames = ["tsec", "soc"];
+        opts.VariableTypes = ["double", "double"];
+        % Specify file level properties
+        opts.ExtraColumnsRule = "ignore";
+        opts.EmptyLineRule = "read";
+        % Import the data
+        cycle = readtable(['Application profiles/', scenario], opts);
+        % Downsample the data from (assumed) 1 second resolution
+        cycle = cycle(1:30:height(cycle), :);
+        % Ensure SOC beginning and end points are identical
+        assert(cycle.soc(1) - cycle.soc(end) < 1e-2, "Start and end SOC of application profile must be within 1% for realistic simulation.")
+        % Repeat the profile for 1 year
+        cycle = repmat(cycle, ceil((3600*24*365)/cycle.tsec(end)), 1);
+        cycle.tsec(:) = 0:30:(30*(height(cycle)-1));
+        cycle = cycle(cycle.tsec <= 3600*24*365, :);
+    else
+        error('Application file-type not recognized.')
+    end
+
+    % Extract stress statistics for this scenario
+    cycleLengthDays = round(cycle.tsec(end)./(24*3600));
+    for iClimate = 1:length(TdegC)
+        % 7 output stressors x number of days
+        stressors = zeros(ceil(cycleLengthDays/dtdays), 7);
+        % temperature
+        if iscell(TdegC)
+            TdegC_cycle = TdegC{iClimate};
+            t_sec_TdegC = [1:8760].*3600;
+            cycle.TdegC = makima(t_sec_TdegC, TdegC_cycle, cycle.tsec);
+        else
+            cycle.TdegC = repmat(TdegC, length(cycle.tsec), 1);
+        end
+        disp("Processing application data...")
+        if height(cycle) > 1e5
+            flagwb = 1;
+            wb = waitbar(0, "Processing application data..."); pause(1e-9);
+            tic;
+        else
+            flagwb = 0;
+        end
+
+        for iDay=1:dtdays:cycleLengthDays
+            if flagwb
+                secs_elapsed = toc;
+                waitbar(iDay/cycleLengthDays, wb, "Calculating daily battery stress...")
+            end
+            % Break the long supercycle up into a smaller cycle for this timestep
+            startDaySubcycle = mod(iDay-1, cycleLengthDays);
+            [~, ind_start] = min(abs(cycle.tsec - startDaySubcycle*24*3600));
+            [~, ind_end] = min(abs(cycle.tsec - (startDaySubcycle + dtdays)*24*3600));
+            subcycle = struct();
+            subcycle.tsec  = cycle.tsec(ind_start:ind_end) - cycle.tsec(ind_start);
+            subcycle.t = subcycle.tsec ./ (24*3600);
+            subcycle.soc   = cycle.soc(ind_start:ind_end);
+            subcycle.TdegC   = cycle.TdegC(ind_start:ind_end);
+            if length(subcycle.t) == 1
+                error('help!')
+            end
+            % Thermal management
+            switch thermalManagement
+                case 'During charge/discharge'
+                    mask = abs([0; diff(subcycle.soc)]) > 1e-6;
+                    mask_lowT = mask & subcycle.TdegC < thermalManagement_MinT;
+                    mask_highT = mask & subcycle.TdegC > thermalManagement_MaxT;
+                    subcycle.TdegC(mask_lowT) = thermalManagement_MinT;
+                    subcycle.TdegC(mask_highT) = thermalManagement_MaxT;
+                case 'Always'
+                    subcycle.TdegC(subcycle.TdegC < thermalManagement_MinT) = thermalManagement_MinT;
+                    subcycle.TdegC(subcycle.TdegC > thermalManagement_MaxT) = thermalManagement_MaxT;
+            end
+            % Extract stressors for this timestep
+            subcycle = struct2table(subcycle);
+            subcycle = func_LifeMdl_StressStatistics(subcycle);
+            % input stressors into the matrix
+            stressors(ceil(iDay/dtdays), :) = [subcycle.dt, subcycle.dEFC, subcycle.TdegK, subcycle.soc, subcycle.Ua, subcycle.dod, subcycle.Crate];
+        end
+        % Fix last day (dt != dtdays, if cycle length isn't perfect)
+        if iDay == cycleLengthDays
+            subcycle.dt = dtdays;
+        end
+        stressors = array2table(stressors, 'VariableNames', {'dt','dEFC','TdegK','soc','Ua','dod','Crate'});
+        if flagwb
+            % close the wait bar
+            close(wb);
+        end
+            
+        % Step through battery types
+        for iModel = idxModels
+            model = validModels{iModel};
+            disp("Running " + model + " battery life simulation...")
+            lifeMdl = func_LifeMdl_LoadParameters(model);
+
+            % Run life simulation
+            lifeSim = runLifeSim(lifeMdl, stressors, tYears);
+    
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            % 5) Plot outputs
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            figure; tl = tiledlayout(2,1);
+            if iscell(TdegC)
+                title(tl, sprintf("Model: %s\nApplication: %s\nClimate: %s", model, scenario, climates{iClimate}), 'Interpreter', 'none')
+            else
+                title(tl, sprintf("Model: %s\nApplication: %s\nClimate: %d C", model, scenario, TdegC), 'Interpreter', 'none')
+            end
+            switch model
+                case 'NMC|Gr'
+                    nexttile; box on; hold on;
+                    plot(lifeSim.t./365, lifeSim.q, '-k', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, lifeSim.q_LLI, ':b', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, lifeSim.q_LAM, ':m', 'LineWidth', 1.5)
+                    xlabel('Time (years)'); ylabel('Relative discharge capacity');  axis([0 tYears 0.7 1.02])
+                    legend('Overall', 'Lithium Inventory', 'Active Material', 'Location', 'southwest')
+                    nexttile; box on; hold on;
+                    plot(lifeSim.t./365, lifeSim.r, '-k', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, lifeSim.r_LLI, ':b', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, lifeSim.r_LAM, ':m', 'LineWidth', 1.5)
+                    xlabel('Time (years)'); ylabel('Relative DC resistance');  axis([0 tYears 0.98 4])
+                    legend('Overall', 'Lithium Inventory', 'Active Material', 'Location', 'northwest')
+                case 'LFP|Gr'
+                    nexttile; box on; hold on;
+                    plot(lifeSim.t./365, lifeSim.q, '-k', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, 1-lifeSim.qLossCal, ':b', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, 1-lifeSim.qLossCyc, ':m', 'LineWidth', 1.5)
+                    xlabel('Time (years)'); ylabel('Relative discharge capacity');  axis([0 tYears 0.7 1.02])
+                    legend('Overall', 'Calendar degradation', 'Cycling degradation', 'Location', 'southwest')
+                    nexttile; box on; hold on;
+                    plot(lifeSim.t./365, lifeSim.r, '-k', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, 1+lifeSim.rGainCal, ':b', 'LineWidth', 1.5)
+                    plot(lifeSim.t./365, 1+lifeSim.rGainCyc, ':m', 'LineWidth', 1.5)
+                    xlabel('Time (years)'); ylabel('Relative DC resistance');  axis([0 tYears 0.98 2])
+                    legend('Overall', 'Calendar degradation', 'Cycling degradation', 'Location', 'northwest')
+            end
+    
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            % 5) Store outputs in struct
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            simulations(simIdx).model = model;
+            simulations(simIdx).application = scenario;
+            if iscell(TdegC)
+                simulations(simIdx).climate = climates{iClimate};
+            else
+                simulations(simIdx).climate = sprintf("%d Celcius", TdegC);
+            end
+            simulations(simIdx).thermalManagement = thermalManagement;
+            switch thermalManagement
+                case {'During charge/discharge'; 'Always'}
+                    simulations(simIdx).MaxT = thermalManagement_MaxT;
+                    simulations(simIdx).MinT = thermalManagement_MinT;
+            end
+            simulations(simIdx).results = lifeSim;
+            simIdx = simIdx + 1;
+        end
+    end
+end
+
+disp(" ")
+disp('Exporting simulation results to file.')
+fnameout = "Simulation results\battery life simulation " + string(datetime("now", "Format", "yyyy-MM-dd-HH-mm-ss")) + ".xls";
+warning('off')
+for i = 1:length(simulations)
+    c = {...
+        'Model:', simulations(i).model;...
+        'Application:', simulations(i).application;...
+        'Climate:', simulations(i).climate;...
+        'Thermal Management:', simulations(i).thermalManagement};
+    switch thermalManagement
+        case {'During charge/discharge', 'Always'}
+            c = [c; {'Min Temp.:', simulations(i).MinT;...
+                     'Max Temp.:', simulations(i).MaxT}];
+    end
+    if i == 1
+        writecell(c, fnameout, 'Sheet', i)
+    else
+        writecell(c, fnameout, 'Sheet', i, 'WriteMode', 'append')
+    end
+    switch simulations(i).model
+        case 'NMC|Gr'
+            columnheaders = {'t','EFC','q','q_LLI','q_LAM','r','r_LLI','r_LAM'};
+            writecell(columnheaders, fnameout, 'Sheet', i, 'WriteMode', 'append')
+            columnvalues = [...
+                simulations(i).results.t';...
+                simulations(i).results.EFC';...
+                simulations(i).results.q';...
+                simulations(i).results.q_LLI';...
+                simulations(i).results.q_LAM';...
+                simulations(i).results.r';...
+                simulations(i).results.r_LLI';...
+                simulations(i).results.r_LAM';...
+                ]';
+            writematrix(columnvalues, fnameout, 'Sheet', i, 'WriteMode', 'append')
+        case 'LFP|Gr'
+            columnheaders = {'t','EFC','q','qLossCal','qLossCyc','r','rGainCal','rGainCyc'};
+            writecell(columnheaders, fnameout, 'Sheet', i, 'WriteMode', 'append')
+            columnvalues = [...
+                simulations(i).results.t';...
+                simulations(i).results.EFC';...
+                simulations(i).results.q';...
+                simulations(i).results.qLossCal';...
+                simulations(i).results.qLossCyc';...
+                simulations(i).results.r';...
+                simulations(i).results.rGainCal';...
+                simulations(i).results.rGainCyc';...
+                ]';
+            writematrix(columnvalues, fnameout, 'Sheet', i, 'WriteMode', 'append')
+    end
+end
+
+%% Helper functions:
+function lifeSim = runLifeSim(lifeMdl, stressors, t_years)
+lifeSim.t = [0:stressors.dt(1):(t_years*365)]';
+lifeSim.EFC = zeros(size(lifeSim.t));
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 5) Life simulation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+xx0 = zeros(lifeMdl.n_states, 1);                             % initialize life model states to zero at beginning of life
+xx0(3) = 1e-8;
+xx  = zeros(lifeMdl.n_states, length(lifeSim.t)); % initialize memory for storing states throughout life
+yy = ones(lifeMdl.n_outputs, length(lifeSim.t)); % initialize memory for storing outputs throughout life
+% customize initial states/outputs for each model
+switch lifeMdl.model
+    case 'NMC|Gr'
+        yy(5, :) = 1.01; % Initialize LAM for NMC|Gr model
+    case 'LFP|Gr'
+        yy([2,3,5,6], :) = 0;
+end
+% wait bar
+wb = waitbar(0, sprintf('Running %s simulation', lifeMdl.model)); pause(1e-9);
+tic;
+% time marching loop
+steps_life_sim = length(lifeSim.t)-1;
+for j=1:steps_life_sim
+    secs_elapsed = toc;
+    waitbar(j/steps_life_sim, wb,...
+        sprintf('Running simulation (step %d/%d, %d seconds left)', j, steps_life_sim, floor(((secs_elapsed/j)*steps_life_sim)-secs_elapsed))); 
+    pause(1e-9);
+    
+    % Extract stressors for this timestep
+    idxRow = mod(j-1, size(stressors, 1)) + 1;
+    
+    % Cumulatively sum equivalent full cycles
+    lifeSim.EFC(j+1) = lifeSim.EFC(j) + stressors.dEFC(idxRow);
+    
+    % integrate state equations to get updated states, xx
+    xx(:,j+1) = func_LifeMdl_UpdateStates(lifeMdl, stressors(idxRow, :), xx0(:, 1));
+    
+    % calculate output equations to get life model outputs, yy
+    yy(:,j+1) = func_LifeMdl_UpdateOutput(lifeMdl, xx(:, j+1));
+    
+    % store states for next timestep
+    xx0(:,1) = xx(:,j+1); 
+end
+% Store results
+switch lifeMdl.model
+    case 'NMC|Gr'
+        lifeSim.q           = yy(1,:)';
+        lifeSim.q_LLI       = yy(2,:)';
+        lifeSim.q_LLI_t     = yy(3,:)';
+        lifeSim.q_LLI_EFC   = yy(4,:)';
+        lifeSim.q_LAM       = yy(5,:)';
+        lifeSim.r           = yy(6,:)';
+        lifeSim.r_LLI       = yy(7,:)';
+        lifeSim.r_LLI_t     = yy(8,:)';
+        lifeSim.r_LLI_EFC   = yy(9,:)';
+        lifeSim.r_LAM       = yy(10,:)';
+    case 'LFP|Gr'
+        lifeSim.q           = yy(1,:)';
+        lifeSim.qLossCal    = yy(2,:)';
+        lifeSim.qLossCyc    = yy(3,:)';
+        lifeSim.r           = yy(4,:)';
+        lifeSim.rGainCal    = yy(5,:)';
+        lifeSim.rGainCyc    = yy(6,:)';
+end
+% close the wait bar
+close(wb);
+end

python/battery_profile.py

@@ -0,0 +1,298 @@
+import pandas as pd
+import numpy as np
+import math
+import datetime
+import random
+
+class BatProfile:
+    # attributes
+    intervals = None
+    df = None
+    soc_sequence_list = []
+    soc_time_sequence_list = []
+    temp_sequence_list = []
+    temp_time_sequence_list = []
+
+    # initialize dataframe
+    def __init__(self):
+        # process dataframe
+        # file_path = "/home/richard/Projects/06_research/battery_degradation_study/battery-anomaly-detection/ISS_data/EP_Battery.Thing_HMD8310.csv"
+        file_path = "/home/richard/Projects/06_research/battery_degradation_study/BLAST-Lite/data/EP_Battery.Thing_HMD8310.csv"
+        fields = ['PACK1_CRIDATA_SOC', 'time']
+        df = pd.read_csv(file_path, skipinitialspace=True, usecols=fields)
+        df['time'] = pd.to_datetime(df['time'])
+        # filter only 2023 data
+        threshold_date = pd.to_datetime('2023-01-01').tz_localize('UTC')
+        df = df[df['time'] >= threshold_date].reset_index(drop = True) 
+        df[fields[0]] = df[fields[0]].replace(86, 85)
+        self.df = df
+
+
+    # methods
+    # obtain clean intervals from data
+    def init_intervals(self):
+        # this dataset contains the full data for all 7 packs
+        df = self.df
+        fields = ['PACK1_CRIDATA_SOC', 'time']
+        
+        def find_intervals_below_threshold(data, threshold):
+            below_threshold = data < threshold
+            # shift series by one, then take bitwise AND, only start points will be 1
+            starts = np.where(below_threshold & ~np.roll(below_threshold, 1))[0]
+            ends = np.where(below_threshold & ~np.roll(below_threshold, -1))[0]
+
+            if below_threshold[0]: # case when first value already is in interval
+                starts = np.insert(starts, 0, 0)
+            if below_threshold[len(below_threshold)-1]: # case when last value is also in interval
+                ends = np.append(ends, len(data) - 1)
+
+            intervals = list(zip(starts, ends))
+
+            return intervals
+
+        # intervals time length
+        def filter_time_length(intervals):
+            return [(x,y) for x,y in intervals if ((df[fields[1]][y] - df[fields[1]][x]) > datetime.timedelta(minutes=10) and 
+                                                (df[fields[1]][y] - df[fields[1]][x]) < datetime.timedelta(hours=3))]
+
+        # intervals depth
+        def interval_depth(intervals):
+            # Find the minimum value within the specified range
+            return [(x,y) for x,y in intervals if (np.min(df[fields[0]][x:y+1]) < 75) and (np.min(df[fields[0]][x:y+1]) > 5)]
+        
+        def has_only_horizontal_line(series, start_index, end_index):
+            interval = series[start_index+1:end_index]
+            gradient_series = np.gradient(interval)
+            return all(value == 0 for value in gradient_series)
+
+        def filter_only_horizontal(intervals):
+            return [(x,y) for x,y in intervals if not has_only_horizontal_line(df['PACK1_CRIDATA_SOC'], x, y)]
+
+        def is_valley(time_series):
+            gradient = np.gradient(time_series)
+            has_negative_value = np.any(gradient < 0)
+            has_positive_value = np.any(gradient > 0)
+            return has_negative_value and has_positive_value
+
+        def filter_valley(intervals):
+            return [(x,y) for x,y in intervals if is_valley(df['PACK1_CRIDATA_SOC'][x:y])]
+
+        bounding_threshold = 80
+        intervals = find_intervals_below_threshold(df[fields[0]], bounding_threshold)
+        intervals = filter_time_length(intervals)
+        intervals = interval_depth(intervals)
+        intervals = filter_only_horizontal(intervals)
+        intervals = filter_valley(intervals)
+        self.intervals = intervals
+
+    # method to ensure that each interval soc begins and ends at 85
+    # this means extending both the soc and time sequences
+    def preprocess_soc_intervals(self):
+        intervals = self.intervals
+        df = self.df
+        raw_time_list = [df["time"][start:end].reset_index(drop = True) for start,end in intervals]
+        raw_soc_list = [df['PACK1_CRIDATA_SOC'][start:end].reset_index(drop = True) for start,end in intervals]
+
+        def extend_segment(discharge_soc_list, discharge_time_list, index):
+            soc_sequence = discharge_soc_list[index]
+            discharge_time_sequence_datetime = discharge_time_list[index]
+
+            start_time = discharge_time_sequence_datetime[0]
+            discharge_time_sequence = [ int((time - start_time).total_seconds()) for time in discharge_time_sequence_datetime ]
+            num_points = 5
+            num_extrapolate = 100
+
+            # extend the series in the beginning to 85%
+            coefficients = np.polyfit(np.arange(num_points), soc_sequence[:num_points], 1)
+            extended_points = np.polyval(coefficients, np.arange(-1, -num_extrapolate-1, -1))
+            extended_points = np.clip(extended_points, None, 85) # ensure values reach 85, but not over
+            # Find the index where differences start repeating
+            repeating_index = np.where(np.diff(extended_points) != 0)[0][-1] + 2
+            # Truncate the array to remove repeating values at the end
+            extended_points = extended_points[:repeating_index]
+            extended_len_1 = len(extended_points)
+            extended_soc_sequence = np.concatenate((extended_points[::-1], soc_sequence))
+            discharge_time_sequence = [ 60 * time for time in range(-extended_len_1, 0)] + discharge_time_sequence
+
+            # extend the series after the end to 85%
+            coefficients = np.polyfit(np.arange(-num_points, 0), soc_sequence[-num_points:], 1)
+            extended_points = np.polyval(coefficients, np.arange(1, num_extrapolate+1))
+            extended_points = np.clip(extended_points, None, 85) # ensure values reach 85, but not over
+            # Find the index where differences start repeating
+            repeating_index = np.where(np.diff(extended_points) != 0)[0][-1] + 2
+            # Truncate the array to remove repeating values at the end
+            extended_points = extended_points[:repeating_index]
+            extended_len = len(extended_points)
+            extended_soc_sequence = np.concatenate((extended_soc_sequence, extended_points,))
+            end_time = discharge_time_sequence[-1]
+            discharge_time_sequence = discharge_time_sequence + [ end_time + 60 * time for time in range(1,extended_len+1)] 
+            # reset index to start at 0 
+            discharge_time_sequence = [ time + 60 * extended_len_1 for time in discharge_time_sequence ]
+
+            # return the modified soc and time series
+            return extended_soc_sequence, discharge_time_sequence
+        
+        # process intervals to start and end at 85
+        soc_time_sequence_list = []
+        soc_sequence_list = []
+        for index in range(len(intervals)):
+            soc_sequence, time_sequence = extend_segment(raw_soc_list, raw_time_list, index)
+            soc_sequence_list.append(soc_sequence)
+            soc_time_sequence_list.append(time_sequence)
+        
+        # save into class variable
+        self.soc_sequence_list = soc_sequence_list
+        self.soc_time_sequence_list = soc_time_sequence_list
+
+    # method to generate temperature sequence for each of the intervals
+    def temp_sequence_generation(self):
+        intervals = self.intervals
+        soc_time_sequence_list = self.soc_time_sequence_list
+
+        def gen_temp_sequence(soc_time_sequence, start_temp):
+            baseline_temp = 24.0
+            max_temp = 30.0
+            temp_rate = 4 / 60
+
+            # we will increase the temperature at a rate of 4/60 degrees per minutes
+            # we will then clip at 30
+            # we use the whole soc discharge+charge interval as the entire warm-up period
+            num_gen = len(soc_time_sequence)
+            x = np.linspace(0, 10, 20)
+            known_gradient = temp_rate
+            y = baseline_temp + known_gradient * x
+            # Use polyfit with deg=1 (linear polynomial) and known gradient as the weight for the first coefficient
+            coefficients = np.polyfit(x, y, deg=1, w=[known_gradient] * len(x))
+            warmup_temp = np.polyval(coefficients, np.arange(1, num_gen))
+            warmup_temp = np.clip(warmup_temp, None, max_temp)
+
+            # we then use the last temperature as the start of the cooldown phase
+            start_temp = warmup_temp[-1]
+            x = np.linspace(0, 10, 20)
+            known_gradient = temp_rate  # Replace this with your known gradient
+            y = start_temp - known_gradient * x
+            coefficients = np.polyfit(x, y, deg=1, w=[known_gradient] * len(x))
+            cooldown_temp = np.polyval(coefficients, np.arange(1, 100))
+            cooldown_temp = np.clip(cooldown_temp, baseline_temp, None)
+            # find where there is no change
+            # take difference in consecutive elements
+            # note that numpy returns a tuple of elements, we only need the first
+            # take the last element of the array
+            repeating_index = np.where(np.diff(cooldown_temp) != 0)[0][-1]
+            # Truncate the array to remove repeating values at the end
+            cooldown_temp = cooldown_temp[:repeating_index]
+
+
+            temp_sequence = np.concatenate((warmup_temp, cooldown_temp))
+            temp_time_sequence = [ 60 * time for time in range(0,len(temp_sequence))] 
+
+            return temp_sequence, temp_time_sequence
+
+        temp_time_sequence_list = []
+        temp_sequence_list = []
+        baseline_temp = 24
+        for index in range(len(intervals)):
+            temp_sequence, time_sequence = gen_temp_sequence(soc_time_sequence_list[index], baseline_temp)
+            temp_sequence_list.append(temp_sequence)
+            temp_time_sequence_list.append(time_sequence)
+
+        self.temp_sequence_list = temp_sequence_list
+        self.temp_time_sequence_list = temp_time_sequence_list
+                
+    # there is a mismatch in number of values between soc and temp
+    # we will pad soc sequence to match that of temp
+    def process_soc_time(self):
+        # we will pad the end soc values with 85 
+        # ensure that the number of values matches that of temperature sequence
+        def extend_soc_time(soc_sequence_list, temp_time_sequence_list, index):
+            previous_soc_count = len(soc_sequence_list[index]) 
+            new_soc_count = len(temp_time_sequence_list[index])
+            num_to_generate = new_soc_count - previous_soc_count
+            padding = np.repeat(85, num_to_generate)
+            extended_soc_sequence = np.concatenate((soc_sequence_list[index], padding))
+            return extended_soc_sequence
+
+        temp_time_sequence_list = self.temp_time_sequence_list
+        soc_sequence_list = self.soc_sequence_list 
+        soc_sequence_list = [ extend_soc_time(soc_sequence_list, temp_time_sequence_list, i) for i in range(len(soc_sequence_list))]
+        self.soc_sequence_list = soc_sequence_list
+
+
+    # generate day values
+    def generate_day_values(self, num_discharges):
+        soc_sequence_list = self.soc_sequence_list
+        temp_sequence_list = self.temp_sequence_list
+        temp_time_sequence_list = self.temp_time_sequence_list
+        # function to give which intervals to include
+        # and where in the day to insert these intervals
+        def sample_intervals(time_sequence_list, num_discharges):
+            # sample with repeats from the list of discharge samples
+            selections = np.random.choice(range(len(soc_sequence_list)), num_discharges, replace=True)
+            # create soc, temp and time lists
+            time_list = []
+            for index in selections:
+                time_list.append(time_sequence_list[index])
+            total_day_time = 60 * 60 * 24 # in seconds
+            # function to check for overlap
+            def is_overlap(range1, range2):
+                a, b = range1
+                c, d = range2
+                return not (b <= c or d <= a)
+            event_duration_list = [ time_sequence[-1] - time_sequence[0] for time_sequence in time_list]
+            time_intervals = []
+            iterations = 0
+            max_iterations = 1000 # to ensure that it ends even if candidate cannot be found
+            for event_duration in event_duration_list:
+                while iterations < max_iterations:
+                    iterations += 1
+                    random_start_time = random.randint(0, total_day_time - event_duration)
+                    proposed_range = (random_start_time, random_start_time + event_duration)
+                    if any(is_overlap(proposed_range, time_interval) for time_interval in time_intervals):
+                        continue
+                    else:
+                        time_intervals.append(proposed_range)
+                        break
+            sorted_order = sorted(range(len(time_intervals)), key=lambda i: time_intervals[i][0])
+            selections = [ selections[i] for i in sorted_order ]
+            time_intervals = [ time_intervals[i] for i in sorted_order ]
+            return selections, time_intervals
+
+        # generate day soc values
+        def gen_day_soc(selections, time_intervals, soc_sequence_list, time_sequence_list):
+            # prepare the start of each sequence
+            soc_day_sequence = np.array([85])
+            time_day_sequence = [0]
+            # add each segment of interest
+            for index in range(len(selections)):
+                soc_day_sequence = np.concatenate((soc_day_sequence, soc_sequence_list[selections[index]]))
+                start_time = time_intervals[index][0]
+                time_day_sequence = time_day_sequence + [ start_time + time for time in time_sequence_list[selections[index]]]
+            # finishing touch
+            soc_day_sequence = np.concatenate((soc_day_sequence, np.array([85])))
+            total_day_time = 60 * 60 * 24
+            time_day_sequence = time_day_sequence + [total_day_time]
+            return soc_day_sequence, time_day_sequence
+
+        def gen_day_temp(selections, time_intervals, temp_sequence_list, time_sequence_list):
+            baseline_temp = 24
+            # prepare the start of each sequence
+            temp_day_sequence = np.array([baseline_temp])
+            time_day_sequence = [0]
+            # add each segment of interest
+            for index in range(len(selections)):
+                temp_day_sequence = np.concatenate((temp_day_sequence, temp_sequence_list[selections[index]]))
+                start_time = time_intervals[index][0]
+                time_day_sequence = time_day_sequence + [ start_time + time for time in time_sequence_list[selections[index]]]
+            # finishing touch
+            temp_day_sequence = np.concatenate((temp_day_sequence, np.array([baseline_temp])))
+            total_day_time = 60 * 60 * 24
+            time_day_sequence = time_day_sequence + [total_day_time]
+            return temp_day_sequence, time_day_sequence
+
+        selections, time_intervals = sample_intervals(temp_time_sequence_list, num_discharges)
+        soc_day_sequence, _ = gen_day_soc(selections, time_intervals, soc_sequence_list, temp_time_sequence_list)
+        temp_day_sequence, time_day_sequence = gen_day_temp(selections, time_intervals, temp_sequence_list, temp_time_sequence_list)
+        return soc_day_sequence, temp_day_sequence, time_day_sequence
+
+            
+