Parameters identification providing derivative

2024 3 6
1 个回答
回答已采纳
更新时间：2024 3 7
8 次查看（30 天）

0 个投票

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.

6 个评论
显示 4更早的评论 隐藏 4更早的评论

Torsten
Torsten 2024-3-6
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
Daniel Lotano 2024-3-6
I have measurement data for I.
Torsten
Torsten 2024-3-6
编辑：Torsten 2024-3-6
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
Daniel Lotano 2024-3-6
The model input is the current, hence, V_1 is a function of I, not the opposite.
Torsten
Torsten 2024-3-6
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
Daniel Lotano 2024-3-6
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.

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Star Strider
Star Strider 2024-3-6
Ran in:

1 个投票

Ran in:
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.
.

15 个评论
显示 13更早的评论 隐藏 13更早的评论

Torsten
Torsten 2024-3-6
@Star Strider
I must be defined as a function of t.
Star Strider
Star Strider 2024-3-6
All the variables are implicitly defined as functions of t, and I assume they are vectors, each collected at specific values of t.
Torsten
Torsten 2024-3-6
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
Star Strider 2024-3-6
I initially had all the independent variables as funcitons of t and the integrated result was impossible to code.
Daniel Lotano
Daniel Lotano 2024-3-6
@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
Torsten 2024-3-6
编辑：Torsten 2024-3-6
@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
Star Strider 2024-3-6
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
Daniel Lotano 2024-3-6
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
Star Strider 2024-3-6
@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
Torsten 2024-3-6
But you are the author of the standard example.
Let me just say: Monod Kinetics :-)
Star Strider
Star Strider 2024-3-6
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
Daniel Lotano 2024-3-7
编辑：Daniel Lotano 2024-3-7
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
Star Strider 2024-3-7
Ran in:
Ran in:
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
Daniel Lotano 2024-3-7
@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
Star Strider 2024-3-7
编辑：Star Strider 2024-3-7
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.
.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by