Arrhenius model fitting using fminsearch

I am trying to use the Arrhenius model σ(T) = σ_ref * exp(-Ea/(kT)) to fit a set of data using fminsearch. The data describes the behaviour of the electrical conductivity with respect to temperature.

T (C) = [30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60, 62.5, 65, 67.5, 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95, 97.5, 100];

conductivity (s/m) = [0.36, 0.38, 0.4, 0.41, 0.43, 0.45, 0.47, 0.49, 0.52, 0.54, 0.56, 0.58, 0.61, 0.63, 0.66, 0.68, 0.71, 0.74, 0.77, 0.78, 0.79, 0.79, 0.78, 0.76, 0.75, 0.71, 0.7, 0.66, 0.65];

The σ_ref (Electrical conductivity reference) = 0.36 s/m, Ea = 1.326 ev and k (Boltzmann constant) = 8.617333262145e-5 ev/K are all known values. I generated the code where I used the objective function, created the function of the Arrhenius model and then called fminsearch. I'm getting a fitting equal to 0 at all the temperature range (Figure below).

The code also provide a value prediction of the electrical conductivity for the temperature range = [100,120].

N.B: I converted the Temperature to kelvin in the Arrhenius model.

The code is as below:

T (C) = [30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60, 62.5, 65, 67.5, 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95, 97.5, 100];

conductivity (s/m) = [0.36, 0.38, 0.4, 0.41, 0.43, 0.45, 0.47, 0.49, 0.52, 0.54, 0.56, 0.58, 0.61, 0.63, 0.66, 0.68, 0.71, 0.74, 0.77, 0.78, 0.79, 0.79, 0.78, 0.76, 0.75, 0.71, 0.7, 0.66, 0.65];

x_data = T;

y_data = conductivity;

% Plot

plot(x_data, y_data, 'r+', 'linewidth', 1.5, 'Markersize', 9);

title('Arrhenius Model Measurement');

xlabel('Temperature (°C)');

ylabel('Electrical Conductivity (S/m)');

grid on;

% Initial guess

initial_guess = [1.326, 8.617333262145e-5, 0.36];

%fminsearch

x = fminsearch(@(x) obj_func_Arrhenius(x, x_data, y_data_k), initial_guess);

%parameters

Ea = x(1);

k = x(2);

ElectricalCond_ref = x(3);

%model predictions

x_model = linspace(min(x_data), max(x_data), 100);

y_model = Arrhenius_model(sigma_ref, Ea, x_model, k);

% Plot the measured data and model predictions

hold on;

plot(x_model, y_model, 'b-', 'linewidth', 1.5);

% Predict for temperatures from 90 to 120

new_temperatures = 100:120;

new_predictions = Arrhenius_model(sigma_ref, Ea, new_temperatures, k);

plot(new_temperatures, new_predictions, '--', 'linewidth', 1.5);

legend('Measured Data', 'Model Prediction', 'Predictions for 100-120', 'Location','northwest');

hold off;

% Arrhenius Model function

function ElectCond = Arrhenius_model(sigma_ref, Ea, Temp, k)

ElectCond = sigma_ref .* exp(-Ea./(k .* Temp));

end

% Objective function

function SSE = obj_func_Arrhenius(x, Temp_r, ElectricalCond_r)

Ea = x(1);

k = x(2);

ElectricalCond_ref = x(3);

ElectricalCond_m = Arrhenius_model(ElectricalCond_ref, Ea, Temp_r, k);

SSE = sum((ElectricalCond_m - ElectricalCond_r).^2);

end

I'm not able to figure out the problem here.

Answer 1

[VBBV]

Give a suitable values for intial guess as inputs to objective function, and set values for physical constants in correct units ,

T = [30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60, 62.5, 65, 67.5, 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95, 97.5, 100];
conductivity = [0.36, 0.38, 0.4, 0.41, 0.43, 0.45, 0.47, 0.49, 0.52, 0.54, 0.56, 0.58, 0.61, 0.63, 0.66, 0.68, 0.71, 0.74, 0.77, 0.78, 0.79, 0.79, 0.78, 0.76, 0.75, 0.71, 0.7, 0.66, 0.65];
x_data = T;
y_data_k = conductivity;
% Initial guess
initial_guess = [10.26, 8.617333262145e-3, 1.06];


%fminsearch 
x = fminsearch(@(x) obj_func_Arrhenius(x, x_data, y_data_k), initial_guess)

x = 1×3

0.6401 0.0179 1.1127

%parameters
Ea = x(1);
k = x(2);
ElectricalCond_ref = x(3);
sigma_ref = 1.1;
%model predictions
x_model = linspace(min(x_data), max(x_data), 100);

y_model = Arrhenius_model(sigma_ref, Ea, x_model, k);
% Plot
plot(x_data, y_data_k, 'r+', 'linewidth', 1.5, 'Markersize', 9);
title('Arrhenius Model Measurement');
xlabel('Temperature (°C)');
ylabel('Electrical Conductivity (S/m)');
grid on;

% Plot the measured data and model predictions
hold on;
plot(x_model, y_model, 'b-', 'linewidth', 1.5);
% Predict for temperatures from 90 to 120
new_temperatures = 100:120;
new_predictions = Arrhenius_model(sigma_ref, Ea, new_temperatures, k);
plot(new_temperatures, new_predictions, '--', 'linewidth', 1.5);
legend('Measured Data', 'Model Prediction', 'Predictions for 100-120', 'Location','northwest');
hold off;

% Objective function
function SSE = obj_func_Arrhenius(x, Temp_r, ElectricalCond_r)
    Ea = x(1);
    k = x(2);
    ElectricalCond_ref = x(3);
    
    ElectricalCond_m = Arrhenius_model(ElectricalCond_ref, Ea, Temp_r, k);
    
    SSE = sum((ElectricalCond_m - ElectricalCond_r).^2);
end
% Arrhenius Model function
function ElectCond = Arrhenius_model(sigma_ref, Ea, Temp, k)
    ElectCond = sigma_ref .* exp(-Ea./(k .* Temp));
end

Answer 2

[Torsten]

Imagine you get Ea = 3 and k = 6 in one run and Ea = 12 and k = 24 in a second run. Both solutions give the same quality of fitting, but which one should be preferred ?

We call this "overfitting". The Arrhenius model cannot be fitted with three free parameters, but only two. Those parameters are p1 = sigma_ref and the quotient of Ea and k: p2 = Ea/k (or Ea alone if you fix the value of k to be the Boltzmann constant).

[Ahmad]

Hello,

I really appreciate the help !

I took your suggestion into consideration and i tried to fit the Arrhenius model with 2 parameters, sigma_ref and Ea and by giving the value of k = 8.617 333 262 x 10-5 eV/K. Still the plot shows 0 fitting at all temperature.

Answer 3

[Torsten]

Your model function σ(T) = σ_ref * exp(-Ea/(kT)) does not have a maximum. Thus it is not suited to reflect your measurement data.

That's the best you can get:

T= [30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60, 62.5, 65, 67.5, 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95, 97.5, 100];
conductivity  = [0.36, 0.38, 0.4, 0.41, 0.43, 0.45, 0.47, 0.49, 0.52, 0.54, 0.56, 0.58, 0.61, 0.63, 0.66, 0.68, 0.71, 0.74, 0.77, 0.78, 0.79, 0.79, 0.78, 0.76, 0.75, 0.71, 0.7, 0.66, 0.65];
T_ref = 273.15;
k = 8.617333262145e-5;
fun = @(x) x(1)*exp(-x(2)./(k*(T+T_ref)))-conductivity;
x0 = [1 1e-2];
sol = lsqnonlin(fun,x0)

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol = 1×2

13.1054 0.0894

hold on
plot(T,conductivity,'o')
plot(T,fun(sol)+conductivity)
hold off
grid on