-
Notifications
You must be signed in to change notification settings - Fork 1
/
EstECMPara.m
224 lines (175 loc) · 6.61 KB
/
EstECMPara.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
function [thetaOpt, paraInfo, vFit] = EstECMPara(t,u,v,varargin)
% Estimate the parameters of an ECM. Postive current is assumed as
% discharge
%
% W.D. Widanage 07/01/2016 (Fade to Black)
% W.D. Widanage 21/12/2018 (Aint my ****)
p = inputParser; % Create an input parse object to handle positional and property-value arguments
% Create variable names and assign default values after checking the value
addRequired(p,'t', @isnumeric);
addRequired(p,'u', @isnumeric);
addRequired(p,'v', @isnumeric);
% Optional parameters
addParameter(p,'Jacobian','on')
addParameter(p,'order',2)
addParameter(p,'plotFit',0,@isnumeric)
addParameter(p,'ecmFitSeriesCap','off') % Set this to on to fit a ECM with a series capacitor
addParameter(p,'dispMsg','on')
% options = optimset('TolFun',1e-5,'MaxFunEvals',2000,'Jacobian','on','Algorithm','trust-region-reflective');
% Re-parse parObj
parse(p,t, u, v, varargin{:})
t = p.Results.t;
u = p.Results.u;
v = p.Results.v;
order = p.Results.order;
plotFit = p.Results.plotFit;
% Ensure inputs are column vectors
t = t(:);
u = u(:);
v = v(:);
% De-trend OCV from the voltage
ocvStart = v(1); tStart = t(1);
ocvEnd = v(end); tEnd = t(end);
ocvRef = interp1([tStart;tEnd],[ocvStart;ocvEnd],t);
volOCVRemoved = v - ocvRef + ocvStart;
OCV = v(1);
% Fit equivalent circuit with series capacitior on time data
if ismember(p.Results.ecmFitSeriesCap,{'on','On'})
dtTmp = diff(t);
dt = [dtTmp(1);dtTmp];
fh = @(theta,u)EquivalentCircuitModel(theta,u,OCV,dt,order); % Model function handle
% Set initial values for ECM parameters
[~, idx] = max(abs(u));
if sign(u(idx)) < 0 % For a charge pulse
Ro = (max(v) - OCV)/abs(u(idx));
OCVacc0 = (v(end) - OCV)/(abs(u(idx))*10); % OCV', Steady state voltage = OCVacc*Area under applied current + OCV
else % For a discharge pulse
Ro = (OCV - min(v))/abs(u(idx));
OCVacc0 = (OCV - v(end))/(abs(u(idx))*10); % OCV', Steady state voltage = OCVacc*Area under applied current + OCV
end
for mm = 1: order
Rp(mm) = 1e-3;
tau(mm) = 10^(mm-1);
end
theta0 = [Ro OCVacc0 Rp tau]';
options.Jacobian = p.Results.Jacobian;
options.dispMsg = p.Results.dispMsg;
[thetaOpt,paraInfo] = LMAlgorithm(fh, theta0, u, v, options);
% [thetaOpt,resNorm] = lsqcurvefit(fh,theta0,u,volOCVRemoved);
% Simulate ECM with estimation data set and estimated optimum parameteres
vFit = EquivalentCircuitModel(thetaOpt,u,OCV,dt,order);
else
% Set initial values for ECM parameters
[~, idx] = max(abs(u));
if sign(u(idx)) < 0 % For a charge pulse
Ro = (max(v) - OCV)/abs(u(idx));
else % For a discharge pulse
Ro = (OCV - min(v))/abs(u(idx));
end
for mm = 1: order
Rp(mm) = 1e-3;
tau(mm) = 5*(mm);
end
theta0 = [Ro Rp tau]';
fh = @(theta,u)ECMStable(theta, u, OCV, t, order); % Model function handle
options.Jacobian = p.Results.Jacobian;
options.dispMsg = p.Results.dispMsg;
[thetaOpt,paraInfo] = LMAlgorithm(fh, theta0, u, volOCVRemoved, options);
% [thetaOpt,resNorm] = lsqcurvefit(fh,theta0,u,volOCVRemoved);
yECMTmp = ECMStable(thetaOpt, u, OCV, t, order);
vFit = yECMTmp + ocvRef - ocvStart;
end
errorECM = v - vFit;
paraInfo.resNorm = rms(errorECM);
paraInfo.pkErrECM = max(abs(errorECM));
if plotFit == 1
figure();
ax(1) = subplot(2,1,1);
t = t - t(1);
plot(t,u,'. -');
ylabel('Current (A)');
ax(2) = subplot(2,1,2);
plot(t,v,'- .',t,vFit,'-')
xlabel('Time (s)'); ylabel('Voltage (V)'); legend('Measured','ECM sim')
title(['RMSE ECM: ',num2str(paraInfo.resNorm*1000),'mV Pk Err ECM: ', num2str(paraInfo.pkErrECM*1000),'mV']);
linkaxes(ax,'x')
end
function [Vl, J]= EquivalentCircuitModel(theta, Il, OCV, del, order)
%
%
% Parameters are arranged as:
% theta = Ro, OCV', Rp1,..,Rpn, tau1,...,taun
%
% W. D. Widanage 19/10/2013 (Regain)
%
% Change parameter vector to a column vector
theta = theta(:);
% Extract parameters
Ro = theta(1);
OCVacc = theta(2);
RpAll = theta(3:2+order);
tauAll = theta(3+order:2+2*order);
dataPts = length(Il); % Number of data points
% Initialise
Ip = zeros(order,1);
Vl = zeros(dataPts,1);
Integral = 0;
% Initialise for Jacobian
J = zeros(dataPts,length(theta));
dIpdtau = zeros(order,1);
Vl(1) = OCV - Ro*Il(1) - OCVacc*Integral - RpAll'*Ip;
J(1,:) = [-Il(1), -Integral, -Ip', -dIpdtau'];
for ii = 2:dataPts
Ts = del(ii-1); % Sampling interval
expTau = exp(-Ts./tauAll);
for jj = 1:order
Ip(jj,1) = expTau(jj)*(Ip(jj,1) - Il(ii)) + Il(ii);
end
Integral = Integral + (Il(ii) + Il(ii-1))*Ts/2;
Vl(ii) = OCV - Ro*Il(ii) - OCVacc*Integral - RpAll'*Ip;
% Create Jacobian matrix
for jj = 1:order
dIpdtau(jj,1) = expTau(jj)*(dIpdtau(jj,1) + Ip(jj,1)*Ts/tauAll(jj)^2 - Il(ii)*Ts/tauAll(jj)^2);
end
RdIpdtau = RpAll.*dIpdtau;
J(ii,:) = [-Il(ii), -Integral, -Ip', -RdIpdtau'];
end
function [Vl, J] = ECMStable(theta, Il, OCV, timeVec, order)
% Stable equivalent cirucit model, without series capacitor
%
% Parameters are arranged as:
% theta = Ro, Rp1,..,Rpn, tau1,...,taun
%
% W. D. Widanage 04/07/2014 (Serene)
%
dataPts = length(Il); % Number of data points
% Check if weighting function is provided. else initialise to ones
% Change parameter vector to a column vector
theta = theta(:);
% Extract parameters
Ro = theta(1);
RpAll = theta(2:1+order);
tauAll = theta(2+order:1+2*order);
% Initialise
Ip = zeros(order,1);
Vl = zeros(dataPts,1);
% Initialise for Jacobian
J = zeros(dataPts,length(theta));
dIpdtau = zeros(order,1);
Vl(1) = OCV - Ro*Il(1) - RpAll'*Ip;
J(1,:) = [-Il(1), -Ip', -dIpdtau'];
for ii = 2:dataPts
Ts = timeVec(ii)-timeVec(ii-1); % Sampling interval
expTau = exp(-Ts./tauAll);
for jj = 1:order
Ip(jj,1) = expTau(jj)*(Ip(jj,1) - Il(ii)) + Il(ii);
end
Vl(ii) = OCV - Ro*Il(ii) - RpAll'*Ip;
% Create Jacobian matrix
for jj = 1:order
dIpdtau(jj,1) = expTau(jj)*(dIpdtau(jj,1) + Ip(jj,1)*Ts/tauAll(jj)^2 - Il(ii)*Ts/tauAll(jj)^2);
end
RdIpdtau = RpAll.*dIpdtau;
J(ii,:) = [-Il(ii), -Ip', -RdIpdtau'];
end
Vl = Vl(:); % Change to a column vector