Parameters identification providing derivative

Hello, I'd need to perform parameters identification of a battery model. Let's say I have a simplified model that can be described by a single RC branch:

where the parameters I want to identify are R0, R and C (or the time constant tau instead of C). I have the data of V_measured, V_oc and I. The RC branch voltage needs to be calculated using the provided derivative.

Now, my question is: how can I set up the problem to estimate the parameters using the two equations that I have provided? I don't really need the full code, just some workflow would be fine, because I have no idea how to setup a problem like this.

I know that there is a toolbox that performs this kind of parameters estimation, but I'd like to have a little more control on how the algorithm performs.

[Torsten]

For the theoretical equation

C*dV_1/dt = I(t) - V_1(t)/R

do you have a separate analytical equation for I(t), or do you have to use your measurement data here ?

[Daniel Lotano]

I have measurement data for I.

[Torsten]

I know, but I wanted to know if there is a theoretical equation for I as a function of V_1 that could be used in the differential equation or if it is necessary to use your measurement data in the theoretical equation for V_1.

The background is that theoretical model and measurements should be as independent from each other as possible.

[Daniel Lotano]

The model input is the current, hence, V_1 is a function of I, not the opposite.

[Torsten]

So there is no separate theoretical equation for I such that you could solve the two equations

dI/dt = ? or f(t,I(t),V_1(t)) = 0
C*dV_1/dt = I(t) - V_1(t)/R

for I and V_1 ?

[Daniel Lotano]

No, I is a free parameter. I could set any given current in the model, i.e. I can charge the battery with any profile current or I could connect any load that could absorb any current with any profile. There is no real equation to describe that.

Answer 1

[Star Strider]

The differential equation needs to be integrated in order to use . This does not appear to be as straightforward as I would have hoped, however it may do what you want. (It has the virtue of running without error using random data, however that is all I can say for it.) Note thet the time vector needs to be an input as well as the variables acquired at those times.

Try this —

syms Vmeas V_oc I R_0 V_1(t) C R V_10 
sympref('AbbreviateOutput',false);
eqn1 = Vmeas == V_oc - I * R_0 - V_1

eqn1(t) = 

eqn2 = C*diff(V_1) == I - V_1/R

eqn2(t) = 

V_1(t) = dsolve(eqn2, V_1(0) == V_10)

V_1(t) = 

eqn1 = subs(eqn1)

eqn1(t) = 

eqn1 = simplify(eqn1, 500)

eqn1(t) = 

clearvars
% % % R_0 = b(1),  R = b(2),  C = b(3)
fcn = @(b,V_oc,Vmeas,V_10,I,t) Vmeas + exp(t/(b(3)*b(2))).*(V_10 - I*b(2)) + I*b(2) + I*b(1) - V_oc;

t = 0:10;
Vmeas =  randn(size(t));
V_10 = 0;
I =  randn(size(t));
V_oc =  randn(size(t));

B0 = rand(3,1)

B0 = 3x1

0.7564 0.9149 0.8188

B = fminunc(@(b)norm(fcn(b,V_oc,Vmeas,V_10,I,t)), B0)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

B = 3x1

0.7961 2.8837 3.4689

I used fminunc here, it also works with fminsearch, (although producing differnt parameter estimates) and would probably work as well with fmincon if you wanted to constrain the parameters.

.

[Torsten]

@Star Strider

I must be defined as a function of t.

[Star Strider]

All the variables are implicitly defined as functions of t, and I assume they are vectors, each collected at specific values of t.

[Torsten]

The solution of the differential equation is not correct if I is assumed as a constant (as is done in your code).

syms I

instead of

syms I(t)

[Star Strider]

I initially had all the independent variables as funcitons of t and the integrated result was impossible to code.

[Daniel Lotano]

@Star Strider thanks for you answer! Yes, I forgot to mentiont that I obviously have the sampling time, too, so t can easily be an input. Your solution is something that I was thinking about, but I couldn't manage to write down.

The code unfortunately won't run with fminunc as I get the error on fminusub "Objective function is undefined at initial point. Fminunc cannot continue.", but it will run using fminsearch, even though the parameters are far off from the values that I expect. I'm sure that fmincon could solve this issue with some bounds.

It is anyways a great starting point. Huge thank you!

[Torsten]

@Star Strider

I initially had all the independent variables as funcitons of t and the integrated result was impossible to code.

I think it's necessary to use a numerical integrator then with the measurement I-values interpolated to the time when the integrator requests the time derivatives for V_1. Or there is another theoretical equation for I(t) - algebraic or differential.

[Star Strider]

As always, my pleasure!

I might be able to help with fminunc error if I have your data and your initial code. The only possible division-by-zero that I can see is if ‘R’ (‘b(2)’) or ‘C’ (‘b(3)’) are initially defined to be zero. For this reason, I always suggest using rand (or a multiple or added offset value of it) to set the initial esitimates, unless you know in advance what the actual estimated value is likely to be, and never zeros.

The only other possibility are NaN elements in your data. Using either fillmissing or rmmissing can deal with that problem.

[Daniel Lotano]

The data is cleared from NaNs as I use it for several other algorithms. I initialized the parameters as [5e-3, 3e-3, 2], as they are the mean value i get from the toolbox above mentioned.

A minus sign is missing in the function handle in the exp function:

fcn = @(b,V_oc,Vmeas,V_10,I,t) ...
    Vmeas + exp(-t/(b(3)*b(2))).*(V_10 - I*b(2)) + I*b(2) + I*b(1) - V_oc;

I will try to implement the fmincon to avoid negative values, or I'll even try lsqnonlin.

I will update soon.

[Star Strider]

@Daniel Lotano — The missing minus sign could be the reason for the ‘undefined at initial point’ problem. I did not clearly see it. when I hand-coded ‘fcn’.

@Torsten — I gave some thought to numeircal integration, however I was uncertain how to implement it (integral or trapz) and then code it.

.

[Torsten]

But you are the author of the standard example.

Let me just say: Monod Kinetics :-)

[Star Strider]

Thank you!

I considered that, however I could not figure out how to make that work with these equations. I decided to use the analytic solution to the differential equation, and then rearranged that result to get the final function expression.

[Daniel Lotano]

I'm sorry for the late reply, but I had to work quite a bit on the code!

As said, huge thank you to @Star Strider that provided a great starting point, but the implementation didnt work as expected. In fact the equation couldn't catch the dynamics of the system as V_10 was not chaging over time (probably my fault?).

First of all, I decided to implement a 2nd order model, with 2 RC nets, just to have a little better estimation. Then, I calculated the derivatives as finite differences:

for k = 1:length(I)-1

V1(k+1,1) = (I(k) - V1(k)/R1) * dt(k)/tau1*R1 + V1(k);

V2(k+1,1) = (I(k) - V2(k)/R2) * dt(k)/tau2*R2 + V2(k);

end

Hence, I could simply calculate the residual battery voltage to use as objective function as:

res = Voc - I.*R0 - V1 - V2 - Vmeas;

This methodology allowed to capture the dynamics of the charge and the discharge of the imaginary capacitor in the model. Previously I was just getting sharp corners without charge/discharge dynamics.

I attached the files needed to run the example and a subset of the data that I have sampled.

Again, huge thank you for helping out with this!

[Star Strider]

As always, my pleasure!

The ‘V_10’ (actually ) is the initial condition for the integrated solution of so it should not change. My apologies for not clarifying that. It could also be a parameter to be estimated, and I changed my function to do that here.

It is going to be very difficult to fit these data.

I intended to run your code here, and while I can type all the files, I cannot run the main file, even if I copy-paste its full, complete name to the run argument. I have no idea what that problem could be. (That should work, since it has worked successfully when I have used it in other Answers.)

mfiles = dir('*.m');
hyphenstring = repmat('-',1,75);
for k = 1:numel(mfiles)
    filename = mfiles(k).name
    type(filename)
    fprintf('\n\n%s\n',hyphenstring)
end

filename = 'batteryObjectiveFunction.m'

function res = batteryObjectiveFunction(par, Voc, Vmeas, I, t, V_x0)
%BATTERYOBJECTIVEFUNCTION Objective function for the battery
%   
%   Vettore dei parametri da stimare
%   par = [R0 R1 R2 tau1 tau2]'
%
%   Vettore delle condizioni iniziali
%   V_x0 = [V_10 V_20]'
%

R0 = par(1);
R1 = par(2);
R2 = par(3);
tau1 = par(4);
tau2 = par(5);
V_10 = V_x0(1);
V_20 = V_x0(2);

V1 = NaN(size(I));
V2 = NaN(size(I));

V1(1) = V_10;
V2(1) = V_20;

dt = diff(t);

for k = 1:length(I)-1
    V1(k+1,1) = (I(k) - V1(k)/R1) * dt(k)/tau1*R1 + V1(k);
    V2(k+1,1) = (I(k) - V2(k)/R2) * dt(k)/tau2*R2 + V2(k);
end

res = Voc - I.*R0 - V1 - V2 - Vmeas;
end

---------------------------------------------------------------------------

filename = 'estimateVt.m'

function [Vt_est] = estimateVt(par, Voc, I, t, V_x0)
%ESTIMATEVT Stima della tensione ai terminali
%   

R0 = par(1);
R1 = par(2);
R2 = par(3);
tau1 = par(4);
tau2 = par(5);
V_10 = V_x0(1);
V_20 = V_x0(2);

V1 = NaN(size(I));
V2 = NaN(size(I));

V1(1) = V_10;
V2(1) = V_20;

dt = diff(t);

for k = 1:length(I)-1
    V1(k+1,1) = (I(k) - V1(k)/R1) * dt(k)/tau1*R1 + V1(k);
    V2(k+1,1) = (I(k) - V2(k)/R2) * dt(k)/tau2*R2 + V2(k);
end

Vt_est = Voc - I.*R0 - V1 - V2;
end

---------------------------------------------------------------------------

filename = 'working_parame...ptimization.m'

load('kth_data.mat')

% Initial guess
par0 = [5e-3, 3e-3, 3e-3, 15, 600]';

% Capacitors voltage initial conditions
V_x0 = [0, 0]';

% Physical limits
lb = [1e-3, 1e-3, 1e-3, 2, 40];
ub = [10e-3, 40e-3, 40e-3, 100, 3500];

% Optimizer options
opts = optimset('TolFun', 1e-15, 'Display', 'off');

% Function handler
objFun = @(par) batteryObjectiveFunction(par, V_oc_kth, V_kth, I_kth, t_kth, V_x0);

cycle_parameters = lsqnonlin(objFun, par0, lb, ub, opts)

% Plot results
Vestim = estimateVt(cycle_parameters, V_oc_kth, I_kth, t_kth, V_x0);
figure
plot(t_kth, V_kth, 'o'), hold on
plot(t_kth, Vestim, LineWidth=2)
xlabel('Time (s)'), ylabel('Voltage (V)')
legend({'True value', 'Estimation'})

---------------------------------------------------------------------------

load('kth_data.mat')
% whos
% return
I = I_kth;
Vmeas = V_kth;
V_oc = V_oc_kth;
t = t_kth;

figure
yyaxis left
plot(t, I)
ylabel('I(t)')
yyaxis right
plot(t, [Vmeas V_oc])
ylabel('V_{meas}(t),  V_{oc}(t)')
grid
xlabel('t')
title('kth\_data')
legend('I','V_{meas}','V_{oc}', 'Location','best')

run(mfiles(3).name)             % It Cannot Find It To Run It, However It Can Find It To Type It

Error using run
working_parame...ptimization.m not found.

return

% ----- RUNNING MY ORIGINAL CODE — 

load('kth_data.mat')
% whos
% return
I = I_kth;
Vmeas = V_kth;
V_oc = V_oc_kth;
t = t_kth;

figure
yyaxis left
plot(t, I)
ylabel('I(t)')
yyaxis right
plot(t, [Vmeas V_oc])
ylabel('V_{meas}(t),  V_{oc}(t)')
grid
xlabel('t')
legend('I','V_{meas}','V_{oc}', 'Location','best')

% fcn = @(b,V_oc,Vmeas,V_10,I,t) Vmeas + exp(-t/(b(3)*b(2))).*(V_10 - I*b(2)) + I*b(2) + I*b(1) - V_oc;
fcn = @(b,V_oc,Vmeas,I,t) Vmeas + exp(-t/(b(3)*b(2))).*(b(4) - I*b(2)) + I*b(2) + I*b(1) - V_oc;

B0 = rand(4,1)
opts = optimoptions('fminunc', 'MaxFunctionEvaluations',1E4);
[B,fv] = fminunc(@(b)norm(fcn(b,V_oc,Vmeas,I,t)), B0,  opts)
fprintf('%s = %.6E\n', [["R0" "R" "C" "V_10"].' B].')

I will run it on MATLAB Online later, to see how it works.

.

[Daniel Lotano]

@Star Strider I don't really know what's wrong running the files. I have run them on 2 different machines and there are no issues. I just open working_parameters_optimization.m and press Run.

[Star Strider]

The problem is not with the files.

The problem is the Answers site. I have no idea wwhy they will not run here using the run function. They should, and this has worked in other Answers.

I am reporting it to MathWorks as a bug. It needs to be fixed.

EDIT — (7 Mar 2024 at 17:21)

Report a bug Case 06860945 created successfully

Our offices are closed March 7 to March 12 for a staff development program. If this is an urgent matter that cannot wait for our return, call your local office to leave a message. We will respond as soon as possible.

.