Clear Filters
Clear Filters

Paramter optimization by curve fitting

14 views (last 30 days)
Pramodya
Pramodya on 19 May 2024
Commented: Star Strider on 20 May 2024
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.
  2 Comments
Torsten
Torsten on 19 May 2024
We cannot test your code without the data you are trying to fit.
Pramodya
Pramodya on 19 May 2024
Edited: Pramodya on 19 May 2024
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.

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 19 May 2024
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
.
  4 Comments
Pramodya
Pramodya on 20 May 2024
Thank you very much for the feedback. Yes I will check with the function, those ub paramters and let you know the output.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!