nlinfit seems to find only local minimum

2 次查看（过去 30 天）
Silke
Silke 2018-1-30
评论： Torsten 2018-2-1
Hello!
I am busy fitting my experimental data with this code:
function dCCconv = DiffEqSolverALLSEandSO(param, t, k_1n, k_1p, mu_n, phi, FA, WL, lf, thickness, T);
RC = 18e-9;
k_2 = param(1);
k_3 = param(2);
filename1 = 'pulse300.dat'
pathname1 = 'D:\My Documents\Experiments\TRMC\PulseShapes\';
pulse = importfile1([pathname1, filename1]);
assignin('base', 't',t);
assignin('base', 'RC',RC);
[P,I] = max(pulse(:,1)); %find the coordinates of the maximum value of t
[P1, I1] = find(pulse(:,2)>=0);
pulse_begin_time = pulse(I1(1),1);
pulse_begin = min(find(t>=pulse_begin_time));
pulse_end = max(find(t<=P));
[P3, I3] = max(t);
dx = t(2)-t(1);
t_x = t(pulse_begin : pulse_end);
pulses(:,2)=spline(pulse(:,1), pulse(:,2), t_x); % make the pulse and the experimental data have the same time scale
pulses(end+1:I3,2)= 0;
pulses(:,1) = t;
assignin('base', 'pulses',pulses);
pmax = trapz(pulses(:,1), pulses(:,2));
NA = lf*1e-3*WL*1e-9*FA/(6.626e-34*3e8*2.1); % number of absorbed photons per cm^2 (2.1 is in cm^2)
NP = lf*1e-3*WL*1e-9/(6.626e-34*3e8*2.1); % number of incident photons per cm^2
G0 = NA/(pmax*thickness*1e2);
function dx = myode(t, x, pulses, G0, k_1n, k_1p, k_2,k_3, phi, T) %k_3,phi_n,k_1p,k_2,
dx = zeros(2,1);
pulse_actual = interp1(pulses(:,1),pulses(:,2),t );
SE_n = (((T/497))*real(abs(k_1n*1e7))* (t*real(abs(k_1n*1e7)))^((T/497)-1) );
SE_p = ((T/497)*real(abs(k_1p*1e6))* (t*real(abs(k_1p*1e6)))^((T/497)-1) );
dx(1) = real((G0 * pulse_actual*real(abs(phi))) - SE_n * real(x(1)) - abs(k_2*1e-10) * x(1) * x(2))- real(abs(k_3)) * 1e-28 * (x(1)^2 * x(2) + x(1) * x(2)^2) ;
dx(2) = real((G0 * pulse_actual*real(abs(phi))) - SE_p * real(x(2)) - abs(k_2*1e-10) * x(1) * x(2))- real(abs(k_3)) * 1e-28 * (x(1)^2 * x(2) + x(1) * x(2)^2);
end
tspan = [t(2): dx : t(end)];
opts = odeset('RelTol',1e-5,'AbsTol',1e-6);
ic = [0;0];
[t,x] = ode45(@(t,x) myode(t, x, pulses, G0, k_1n, k_1p, k_2, k_3, phi, T), tspan, ic, opts);
CCx = x;
f1x = (real(abs(mu_n)) * CCx(:,1)+ (1/2.06)*real(abs(mu_n))*CCx(:,2))*thickness*1e2 / NP;
f2x = conv(exp(-t/RC),f1x); %convolution to take into acount of the time response of the cavity cell
tr_cav = sum(exp(-t/RC));
f(1:I3) = 100*f2x(1:I3)/tr_cav;
dCCconv = f;
end
using nlinfit:
[xfit,resnorm, Jacob, CovB, MSE] = nlinfit( handles.timecorr, 100*handles.datacorr',@(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature) , handles.x0);
Doing this finds usually a local minimum and I vary the guess parameters until I get a minimum for a chi2 (by hand). Now I tried to fit using lsqcurvefit (just to see if this finds the absolute minimum) using this:
xfit,resnorm] = lsqcurvefit( @(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature), handles.x0, handles.timecorr, 100*handles.datacorr',handles.lb, handles.ub, fitopt ); %
But lsqcurvefit runs infinitely long without returning any result. Is there a better fitting function for my case, or is there something in general wrong? I am not sure how I could improve my fit.
Thanks for your help!

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

回答（2 个）

Matt J
Matt J 2018-1-30
Make sure you follow guidelines for Optimizing a Simulation or Ordinary Differential Equation. If you continue to have difficulty finding global solutions, you should perhaps consider MultiStart or other tools in the Global Optimization Toolbox.
  13 个评论
显示 11更早的评论
Silke
Silke 2018-2-1
编辑：Silke 2018-2-1
But how close does my initial guess have to be to find a good solution? If I would already know the parameters for all my experiments, I would not need to perform the fit. I am quite frustrated now that the fitting does not work better at the moment. Even if I start off with 3 and 2, the fit does not find the good parameters.
Matt J
Matt J 2018-2-1
编辑：Matt J 2018-2-1
It is never apriori clear how close the initial guess has to be, but that is why I recommended the Global Optimization Toolbox for cases (like yours possibly) where it seems overly difficult.
Alternatively, because you only have two unknown parameters, it might be computationally affordable for you to do a sweep of the parameters and plot the objective function as a 2D surface. From this you could see, approximately, where the global minimum lies and get a more informed initial guess.

请先登录，再进行评论。

Silke
Silke 2018-1-31
If I triple the data in the peak, but not in the tail, I will not have a continous dataset anymore, won't I? And in how far will this change my fitting parameters? I am not sure if I fully understand your suggestion.
  10 个评论
显示 8更早的评论
Silke
Silke 2018-2-1
I did it like that:
W = zeros(3600,1);
W(1:49) = 5;
W(50:64) = 1000;
W(65) = 10000000;
W(66:80) = 1000;
W(81:500) = 5;
W(501:3600) = 0.000001;
assignin('base', 'W',W);
opts = statset('nlinfit');
opts.RobustWgtFun = [];
[xfit,resnorm, Jacob, CovB, MSE] = nlinfit( handles.timecorr, 100*handles.datacorr',@(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature, sel2) , handles.x0, opts, 'Weights', W');
But it seems not to change anything on the fit. Is this not correct?
Torsten
Torsten 2018-2-1
If the estimated parameters remain the same, you can check whether a weight matrix is used at all by looking at the output variable "resnorm". resnorm(i) should be sqrt(w(i))*res(i) where "res" is the output without using the weighting option.
Best wishes
Torsten.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Choose a Solver 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by