Parameters renormalization in lsqcurvefit to estimate parameters in ODE15s

Hello!

I am trying to solve a parametric ODE system. I woul like to estimate the parameters by fitting the calculated data with experimental data.

The problem is that these parameters differ by many orders of magnitude.From a physical point of view I have already an idea of the order of magnitude the estimated value should get:

n_abs = 3.3e15;  % (cm-3)
p0 = 9e17; % (cm-3)
Br = 2e-10; % (cm3/s)
kte = 7e-8; %(cm3/s) 
Nte = 1e15; % (cm-3)
gamma = 9.3e-30; %(cm6/s)

Now.. the lsqcurvefit combined with the ode solved can't deal with such values so i was thinking to renormalize the parameters to let the ode solver see only values around unity. Like this I renormalize the initial values of the parameters:

%normalization
t_n_abs = n_abs * 1e-16;
t_p0 = p0 * 1e-18;
t_Br = Br * 1e09;
t_kte = kte * 1e7;
t_Nte = Nte * 1e-16;
t_gamma = gamma * 1e29;

I don't know how to normalize and go back to the right values the parameters and initial conditions within the "kinetics" function.

Here my code that does not work. Can you please tell me how would you solve it?

T1 = readtable('data.txt');
t = T1.t*1e9; %in ns seconds
c = T1.c;
%parameters
n_abs = 3.3e15;  % (cm-3)
p0 = 9e17; % (cm-3)
Br = 2e-10; % (cm3/s)
kte = 7e-8; %(cm3/s) 
Nte = 1e15; % (cm-3)
gamma = 9.3e-30; %(cm6/s)

%normalization
t_n_abs = n_abs * 1e-16;
t_p0 = p0 * 1e-18;
t_Br = Br * 1e09;
t_kte = kte * 1e7;
t_Nte = Nte * 1e-16;
t_gamma = gamma * 1e29;
%initial values
theta0 = [t_p0 , t_Br , t_kte , t_Nte, t_gamma, t_n_abs , t_n_abs + t_p0 , 0 ,0];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@(theta,t) kinetics(theta,t),theta0,t,c,zeros(1,9));
[Cfit,T] = kinetics(theta,t);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %d\n', k1, theta(k1))
end
fprintf('\nNOTE —  The last four ‘theta’ values are the initial conditions.\n\n')

figure(1)
plot(t, c, 'p')
hold on
plot(t, Cfit);
ylim([1e-4 1])
hold off
grid
xlabel('Time (ns)')
ylabel('normalized units')
legend('data', 'fitting', 'Location','N')
set(gca,'yscale','log')

function [C,T] = kinetics(theta,t)
    y0 = theta(6:end);
    [T,y] = ode15s(@odefcn, t, y0);
    
    % Now i have to go back to the normal units
    theta(1) = theta(1) * 1e18;
    theta(2) = theta(2) * 1e-09;
    theta(3) = theta(3) * 1e-7;
    theta(4) = theta(4) * 1e16;
    theta(5) = theta(5) * 1e-29;
    theta(6) = theta(6) * 1e16;
    theta(7) = theta(7) * 1e18;
    theta(8) = theta(8) * 1e16;
    theta(9) = theta(9);
    
    function dydt = odefcn(t,y)
               
        dydt = zeros(4,1);
        dydt(1) = -theta(2).*y(1).*y(2)-theta(3).*y(1).*(theta(4)-y(3))-theta(5)*y(1).*y(2).^2;
        dydt(2) = -theta(2).*y(1).*y(2)-theta(5)*y(1).*y(2).^2;
        dydt(3) = theta(3).*y(1).*(theta(4)-y(3));
        dydt(4) = 0;
        
    end
    C=theta(2)*y(:,1).*y(:,2)/max(theta(2)*y(:,1).*y(:,2));
    %Normalization
    theta(1) = theta(1) * 1e-18;
    theta(2) = theta(2) * 1e09;
    theta(3) = theta(3) * 1e7;
    theta(4) = theta(4) * 1e-16;
    theta(5) = theta(5) * 1e29;
    theta(6) = theta(6) * 1e-16;
    theta(7) = theta(7) * 1e-18;
    theta(8) = theta(8) * 1e-16;
    theta(9) = theta(9);
end

[Star Strider]

Parameters such as the ones described can be difficult to estimate.

I characteristically use the ga (genetic algorithm) to deal with such problems, since it searches the entire parameter space for acceptable parameters. The approach in Is it possible to use Simulated Annealing or Genetic Algorithm for parameter estimation of a parametric System of ODEs? is one such. I have had success with it.

If normalisation is the approach you want, then I would err on normalising the data first, fitting the function, and then using the re-scaled parameters on the original data.

[Antonella Treglia]

Thank you for the tip!

I looked into the code but I have a question. When you define the stiffness function as:

ftns = @(theta) norm(c-kinetics(theta,t));

you use the set of solutions c that in my case corresponds to y. However in my case I do not want to fit directly the y solution but

C=theta(2)*y(:,1).*y(:,2)/max(theta(2)*y(:,1).*y(:,2));

How would you define the ftns function in this case?

Thank you!

[Star Strider]

That would appear to be the output of the ‘kinetics’ function. It should produce a column vector equal in size to the second column in the data file.

I caution against re-defining the parameters by re-scaling them in the ‘kinetics’ function. Scale the data, then scale the parameters (if necessary) after the optimisation converges.

To use the ga function correctly, since you have 9 parameters that you want to scale, the 'InitialPopulationMatrix' in my code should be:

PopSz = 100;
Parms =   9;
randi(1E+4,PopSz,Parms).*[1E-18 1E+9 1E+7 1E-16 1E29 1E-16 1E-18 1E-16 1]

ans = 100×9

1.0e+32 * 0.0000 0.0000 0.0000 0.0000 3.1350 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.4250 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.6760 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.2770 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 6.0790 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.1590 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 9.8760 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.1870 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0690 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 5.4080 0.0000 0.0000 0.0000 0.0000

That will scale it correctly, however it may be extremely difficult to get a reeliable convergence, since the parameter magnitudes are less than eps (about 1E-16) or greater than 1/eps, so those parameters would be either 0 or Inf. The problem are the relative magnitudes, not the magnitudes themselves.

It would lieklly be best to experiment with scaling the differential equation terms and let the parameterrs have reasonable values.

I do not understand the system you are simulating or estimating, or the reason such disparate parameter magnitudes are necessary.

.

[Antonella Treglia]

It is a system of rate equation describing carrier dynamics in semiconductors. The orders of magnitude are well known and accepted in the scientific community. Moreover I can not rescale the equations themselves because the initial conditions are also in the order of 1e14/1e17.

[Star Strider]

I would still see if the ga function can fit them acceptably. The range of the parameters is going to be a problem, regardless. I doubt anything other than some sort of efficient extended-precision (more than 64-bit floating point) implementation of ga is going to be able to deal with them.