add matlab code
This commit is contained in:
pauljgasper
parent
7da74274f4
commit
7da049f0db
Binary file not shown.
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -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.
|
||||
Binary file not shown.
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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.
|
||||
|
||||
Binary file not shown.
Binary file not shown.
|
|
@ -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
|
||||
|
|
@ -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.
|
||||
Binary file not shown.
|
|
@ -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
|
||||
Loading…
Reference in New Issue