Paramter optimization by curve fitting

I tried paramter optimization by curve fitting. I have five paramters to optimize. If possible, I would kindly like to know the below things.

1. I kindly would like to know, how opts could be chnaged to have good fit? I have this code. It says options are not excecuted but doesn't give an error.
2. How to decide lower and upper bounds?Is it soley based on literature?(In my case, I have only one literature similar to my material, but it didnt give a good fit)

Main file

%rng default % For reproducibility

% Constants

AD = 1.2; % Air density

p0 = 101320; % Air pressure

SH = 1.4; % Specific heat of air

PN = 0.7; % Prandtl constant

c0 = 343; % Air velocity

d = 0.05; % Assuming a constant value for d

% Define the model equation

SAC_fun = @(x, f) your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d);

% Load experimental data from Excel file

exp_data = readmatrix('/Users/pramodya/Desktop/PhD/Modelling/MATLAB Final codes/My research/Exp.results/Done/C-1.xlsx'); % Assuming first column contains frequencies and second column contains corresponding SAC values

f = exp_data(:, 1); % Frequency data

ydata = exp_data(:, 2); % SAC experimental data

% Plot the experimental data

figure

plot(f, ydata, 'r*')

xlabel 'Frequency'

ylabel 'SAC'

title('Experimental Data')

% Define the objective function for GA

objectiveFunction = @(x) sum((SAC_fun(x, f) - ydata).^2);

% Define bounds

lb = [0, 0, 0, 0.1]; % Lower bounds for parameters

ub = [2,20, 1e8, 1]; % Upper bounds for parameters

% Define GA options with a specified number of generations

opts = optimoptions(@ga,'SelectionFcn',@selectiontournament, ...

'FitnessScalingFcn',@fitscalingprop);

% Run GA to fit the model to the data

nvars = 4; % Number of variables to optimize

[sol, fval] = ga(objectiveFunction, nvars, [], [], [], [], lb, ub, [], options);

% Plot the fitted curve

figure

responsedata = SAC_fun(sol, f);

plot(f, ydata, 'y*', f, responsedata, 'g-')

legend('Experimental Data', 'Fitted Curve')

xlabel 'Frequency'

ylabel 'SAC'

title('Fitted Response')

% Display optimized parameter values

disp('Optimized Parameter Values:');

disp(['b: ', num2str(sol(1))]);

disp(['T: ', num2str(sol(2))]);

disp(['FR: ', num2str(sol(3))]);

disp(['P: ', num2str(sol(4))]);

% Calculate R-squared value

SS_total = sum((ydata - mean(ydata)).^2);

SS_residual = sum((ydata - responsedata).^2);

R_squared = 1 - (SS_residual / SS_total);

% Display R-squared value

disp(['R-squared Value: ', num2str(R_squared)]);

Function file:

function SAC = your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d)

% Parameters

%VL = x(1);

b = x(1);%shape factor

T = x(2);

FR = x(3);

P = x(4);

% SAC calculation based on provided equations

s = (1/2*b)*sqrt((8*T*AD*2*pi*f)/(FR*P));

T_bes = besselj(1, s*sqrt(-1i)) ./ besselj(0, s*sqrt(-1i));

p = (T*AD/P) * (1 - (2./(s*sqrt(-1i))) .* T_bes).^-1;

T_bes2 = besselj(1, PN^0.5*s*sqrt(-1i)) ./ besselj(0, PN^0.5*s*sqrt(-1i));

K_w = (SH*p0/P) * (SH - ((2*(SH-1))./(PN^0.5*s*sqrt(-1i))) .* T_bes).^-1;

%Characteristic Impedance

Zc = sqrt(K_w.* p);

%Complex wave number

k_c = 2 * pi * f .* sqrt(p./K_w);

%Surface impedance of sample

Zs = -1i * Zc .* cot(k_c * d);

z0 = AD * c0;

%Normalised impedance

Zn = Zs / z0;

%Reflection coefficient

R = (Zs - z0) ./ (Zs + z0);

%Absorption coefficient

SAC = 1 - abs(R).^2;

end

Thank you in advance.

[Torsten]

We cannot test your code without the data you are trying to fit.

[Pramodya]

Thank you for leeting me know. I just attached it.Moreover, the code worked perfectly well for the past literature paper. It gave optimized values as they have mentioned and the fit also was good. But those initial, lower and upper bounds didnt work for me. So I want to know regarding my experimental cruve, how to improve the fit.

Answer 1

[Star Strider]

The ‘options’ problem was easy to fix. Just use the name of the options structure you created in your ga call.

There are a few ways to improve this, however the first would be to check to be sure that the code in ‘your_ATTENBOROUGH_function’ is correct. I can’t check that because I have no idea what you are coding in it. Providing that information would be helpful.

Beyond that, I usually provide an 'InitialPopulationMatrix' and set a 'FunctionTolerance' however thesse are just my preferences.

exp_data = readtable('Control.xlsx')

exp_data = 451x2 table
    Var1      Var2  
    ____    ________

    200     0.009049
    204     0.024959
    208     0.040138
    212     0.027335
    216     0.020957
    220     0.024509
    224     0.023416
    228     0.021707
    232     0.031569
    236     0.035655
    240     0.029074
    244     0.022684
    248     0.025787
    252     0.026313
    256      0.02873
    260     0.027512

f = exp_data.Var1;
ydata = exp_data.Var2;

% figure
% plot(f, ydata)
% grid

% Constants
AD = 1.2;  % Air density
p0 = 101320;  % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.05;  % Assuming a constant value for d
% Define the model equation
SAC_fun = @(x, f) your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d);
% Load experimental data from Excel file
% % exp_data = readmatrix('/Users/pramodya/Desktop/PhD/Modelling/MATLAB Final codes/My research/Exp.results/Done/C-1.xlsx'); % Assuming first column contains frequencies and second column contains corresponding SAC values
% % f = exp_data(:, 1); % Frequency data
% % ydata = exp_data(:, 2); % SAC experimental data
% Plot the experimental data
figure
plot(f, ydata, 'r*')
xlabel 'Frequency'
ylabel 'SAC'
title('Experimental Data')

% Define the objective function for GA
objectiveFunction = @(x) sum((SAC_fun(x, f) - ydata).^2);
% Define bounds
lb = [0, 0, 0, 0.1]; % Lower bounds for parameters
ub = [2,20, 1e8, 1]; % Upper bounds for parameters
% Define GA options with a specified number of generations
opts = optimoptions(@ga,'SelectionFcn',@selectiontournament, ...
                        'FitnessScalingFcn',@fitscalingprop);
% Run GA to fit the model to the data
nvars = 4; % Number of variables to optimize
[sol, fval, ~, output] = ga(objectiveFunction, nvars, [], [], [], [], lb, ub, [], opts);

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

NrGen = output.generations

NrGen = 75

% Plot the fitted curve
figure
responsedata = SAC_fun(sol, f);
plot(f, ydata, 'y*', f, responsedata, 'g-')
legend('Experimental Data', 'Fitted Curve')
xlabel 'Frequency'
ylabel 'SAC'
title('Fitted Response')

% Display optimized parameter values
disp('Optimized Parameter Values:');

Optimized Parameter Values:

disp(['b: ', num2str(sol(1))]);

b: 1.8165

disp(['T: ', num2str(sol(2))]);

T: 0.15516

disp(['FR: ', num2str(sol(3))]);

FR: 4181963.7861

disp(['P: ', num2str(sol(4))]);

P: 0.92243

% Calculate R-squared value
SS_total = sum((ydata - mean(ydata)).^2);
SS_residual = sum((ydata - responsedata).^2);
R_squared = 1 - (SS_residual / SS_total);
% Display R-squared value
disp(['R-squared Value: ', num2str(R_squared)]);

R-squared Value: -0.1502

% % Function file:
function SAC = your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d)
% Parameters
%VL = x(1);
b = x(1);%shape factor
T = x(2);
FR = x(3);
P = x(4);
% SAC calculation based on provided equations
s = (1/2*b)*sqrt((8*T*AD*2*pi*f)/(FR*P));
T_bes = besselj(1, s*sqrt(-1i)) ./ besselj(0, s*sqrt(-1i));
p = (T*AD/P) * (1 - (2./(s*sqrt(-1i))) .* T_bes).^-1;
T_bes2 = besselj(1, PN^0.5*s*sqrt(-1i)) ./ besselj(0, PN^0.5*s*sqrt(-1i));
K_w = (SH*p0/P) * (SH - ((2*(SH-1))./(PN^0.5*s*sqrt(-1i))) .* T_bes).^-1;
%Characteristic Impedance
Zc = sqrt(K_w.* p);
%Complex wave number
k_c = 2 * pi * f .* sqrt(p./K_w);
%Surface impedance of sample
Zs = -1i * Zc .* cot(k_c * d);
z0 = AD * c0;
%Normalised impedance
Zn = Zs / z0;
%Reflection coefficient
R = (Zs - z0) ./ (Zs + z0);
%Absorption coefficient
SAC = 1 - abs(R).^2;
end 

.

[Pramodya]

Thank you so much for the response. In the code, I have used an empirical model to calculate sound absorption coefficnet of composite. So the fucntion has that empirical model. I checked that using previous literature and found that it gives a good fit for previous literature.

And at the same time,as you have given how to use 'InitialPopulationMatrix' and 'FunctionTolerance' to improve the fit?

[Star Strider]

With respect to those, I would do something like this —

exp_data = readtable('Control.xlsx')

exp_data = 451x2 table
    Var1      Var2  
    ____    ________

    200     0.009049
    204     0.024959
    208     0.040138
    212     0.027335
    216     0.020957
    220     0.024509
    224     0.023416
    228     0.021707
    232     0.031569
    236     0.035655
    240     0.029074
    244     0.022684
    248     0.025787
    252     0.026313
    256      0.02873
    260     0.027512

f = exp_data.Var1;
ydata = exp_data.Var2;

% figure
% plot(f, ydata)
% grid

% Constants
AD = 1.2;  % Air density
p0 = 101320;  % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.05;  % Assuming a constant value for d
% Define the model equation
SAC_fun = @(x, f) your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d);
% Load experimental data from Excel file
% % exp_data = readmatrix('/Users/pramodya/Desktop/PhD/Modelling/MATLAB Final codes/My research/Exp.results/Done/C-1.xlsx'); % Assuming first column contains frequencies and second column contains corresponding SAC values
% % f = exp_data(:, 1); % Frequency data
% % ydata = exp_data(:, 2); % SAC experimental data
% Plot the experimental data
figure
plot(f, ydata, 'r*')
xlabel 'Frequency'
ylabel 'SAC'
title('Experimental Data')

% Define the objective function for GA
objectiveFunction = @(x) sum((SAC_fun(x, f) - ydata).^2);
% Define bounds
lb = [0, 0, 0, 0.1]; % Lower bounds for parameters
ub = [2,20, 1e8, 1]; % Upper bounds for parameters
% Define GA options with a specified number of generations
PopSz = 500;
Parms = 4;
opts = optimoptions(@ga,'SelectionFcn',@selectiontournament, ...
                        'FitnessScalingFcn',@fitscalingprop, 'PopulationSize',PopSz, ... 
                        'InitialPopulationMatrix', randi(1E4,PopSz,Parms).*[1E-3 1E-2 1E+4 1E-3], ... 
                        'FunctionTolerance',1E-10, 'MaxGenerations',5E+3);
% Run GA to fit the model to the data
nvars = 4; % Number of variables to optimize
tic
[sol, fval, ~, output] = ga(objectiveFunction, nvars, [], [], [], [], lb, ub, [], opts);

ga stopped because it exceeded options.MaxGenerations.

toc

Elapsed time is 82.371415 seconds.

NrGen = output.generations

NrGen = 1000

% Plot the fitted curve
figure
responsedata = SAC_fun(sol, f);
plot(f, ydata, 'y*', f, responsedata, 'g-')
legend('Experimental Data', 'Fitted Curve')
xlabel 'Frequency'
ylabel 'SAC'
title('Fitted Response')

% Display optimized parameter values
disp('Optimized Parameter Values:');

Optimized Parameter Values:

disp(['b: ', num2str(sol(1))]);

b: 1.7701

disp(['T: ', num2str(sol(2))]);

T: 0.015362

disp(['FR: ', num2str(sol(3))]);

FR: 4222619.0327

disp(['P: ', num2str(sol(4))]);

P: 0.98243

% Calculate R-squared value
SS_total = sum((ydata - mean(ydata)).^2);
SS_residual = sum((ydata - responsedata).^2);
R_squared = 1 - (SS_residual / SS_total);
% Display R-squared value
disp(['R-squared Value: ', num2str(R_squared)]);

R-squared Value: -0.15005

% % Function file:
function SAC = your_ATTENBOROUGH_function(x, f, AD, p0, SH, PN, c0, d)
% Parameters
%VL = x(1);
b = x(1);%shape factor
T = x(2);
FR = x(3);
P = x(4);
% SAC calculation based on provided equations
s = (1/2*b)*sqrt((8*T*AD*2*pi*f)/(FR*P));
T_bes = besselj(1, s*sqrt(-1i)) ./ besselj(0, s*sqrt(-1i));
p = (T*AD/P) * (1 - (2./(s*sqrt(-1i))) .* T_bes).^-1;
T_bes2 = besselj(1, PN^0.5*s*sqrt(-1i)) ./ besselj(0, PN^0.5*s*sqrt(-1i));
K_w = (SH*p0/P) * (SH - ((2*(SH-1))./(PN^0.5*s*sqrt(-1i))) .* T_bes).^-1;
%Characteristic Impedance
Zc = sqrt(K_w.* p);
%Complex wave number
k_c = 2 * pi * f .* sqrt(p./K_w);
%Surface impedance of sample
Zs = -1i * Zc .* cot(k_c * d);
z0 = AD * c0;
%Normalised impedance
Zn = Zs / z0;
%Reflection coefficient
R = (Zs - z0) ./ (Zs + z0);
%Absorption coefficient
SAC = 1 - abs(R).^2;
end 

This runs, however I still believe that you need to check the code in ‘your_ATTENBOROUGH_function’ as well as the upper parameter limits (and the scaling vector ‘[1E-3 1E-2 1E+4 1E-3]’ used with the initial populatioon matrix if the ‘ub’ vector changes), since it does not even come close to fitting the data.

.

[Pramodya]

Thank you very much for the feedback. Yes I will check with the function, those ub paramters and let you know the output.

[Star Strider]

My pleasure!