Paramter optimization using fmincon
1 view (last 30 days)
Show older comments
Dear all,
Pls kindly help me with this! I have now engaged this for three months and still couldnt get a solution. I have the following code which is working. I use that to optimize five paramters. The thing is when I run, it gives some optimized paramters but the plot does not show the modeld data curve fit to the experimental data. When I seperatly checked the model with the given optimzied paramters, it gives a closer plot. But in the code, the graph does show almost zero for all frequeices for the SAC_model. Pls help me to correct it.
clc;
clear;
close all;
% Load experimental data from Excel file
exp_data = readmatrix('Target_SAC.xlsx'); % Assuming first column contains frequencies and second column contains corresponding SAC values
f_exp = exp_data(:, 1); % frequency data
SAC_exp = exp_data(:, 2); % experimental SAC data
% Define SAC function
SAC_fun = @(x, f) your_SAC_function(x, f);
% Initial guess for parameters
initial_guess = [0.0001, 0.00015, 1.02, 20000, 0.75];
% Lower and upper bounds for parameters
lb = [0.00006, 0.0001, 1.0, 10000, 0.65];
ub = [0.00018, 0.0004, 1.06, 50000, 0.8];
% Optimization using fmincon
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp');
[x_opt, fval, exitflag, output] = fmincon(@(x) objective_fun(x, SAC_fun, SAC_exp, f_exp), initial_guess, [], [], [], [], lb, ub, [], options);
% Display optimization process information
disp(output);
% Display optimized parameter values in command window
disp('Optimized Parameter Values:');
disp(['VL: ', num2str(x_opt(1))]);
disp(['TL: ', num2str(x_opt(2))]);
disp(['T: ', num2str(x_opt(3))]);
disp(['FR: ', num2str(x_opt(4))]);
disp(['P: ', num2str(x_opt(5))]);
disp(['Optimized Objective Value: ', num2str(fval)]);
disp(['Exit Flag: ', num2str(exitflag)]);
disp(['Number of Function Evaluations: ', num2str(output.funcCount)]);
% Plot optimized SAC curve and experimental data
f = 1:1:5500; % Generate frequencies from 1 to 5500 with a step size of 1
SAC_opt = your_SAC_function(x_opt, f); % Compute optimized SAC values
figure;
plot(f, SAC_opt, 'r-', 'LineWidth', 2);
hold on;
plot(f_exp, SAC_exp, 'bo'); % Plot experimental data points
title('Optimized SAC vs. Experimental Data');
xlabel('Frequency');
ylabel('SAC');
legend('Optimized SAC', 'Experimental Data');
grid on;
function SAC = your_SAC_function(x, f)
% Parameters
VL = x(1);
TL = x(2);
T = x(3);
FR = x(4);
P = x(5);
% Constants
AD = 1.2; % Air density
DV = 1.84e-5; % Dynamic viscosity
p0 = 101320; % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.03; % Assuming a constant value for d
% SAC calculation based on provided equations
E3 = sqrt(1 + (4 * 1i * T * T * DV * 2 * pi * f * AD) / (FR * FR * VL * VL * p0 * p0));
E4 = FR * p0 ./ (1i * 2 * pi * f * AD * T);
E1 = T * AD .* (1 + E4 .* E3);
E6 = 1i * 8 * 2 * pi * f ./ (PN * AD * 2 * pi * f * TL * TL);
E7 = sqrt(1 + (1i * PN * AD * 2 * pi * f * TL * TL / 16 * DV));
E8 = 1./(1 - (E6 .* E7));
E2 = (SH * p0) ./ (SH - ((SH - 1).* E8));
KC = 2 * pi * f .* sqrt(E1 ./ E2);
ZC = sqrt(E1 .* E2);
ZS = -1i * ZC .* cot(KC .* d) / P;
R = (ZS - AD * c0) ./ (ZS + AD * c0);
SAC = 1 - abs(R).^2;
end
function error = objective_fun(x, SAC_fun, SAC_exp, f_exp)
% Compute SAC values using the model
SAC_model = SAC_fun(x, f_exp);
% Compute the squared errors at each frequency
squared_errors = (SAC_exp - SAC_model).^2;
% Compute the average squared error
error = mean(squared_errors);
end
2 Comments
Accepted Answer
Torsten
on 13 Mar 2024
Moved: Torsten
on 13 Mar 2024
It seems that the model is not able to reproduce your measurement data. The curves look qualitatively different.
clc;
clear;
close all;
% Load experimental data from Excel file
exp_data = readmatrix('Target_SAC.xlsx'); % Assuming first column contains frequencies and second column contains corresponding SAC values
f_exp = exp_data(:, 1); % frequency data
SAC_exp = exp_data(:, 2); % experimental SAC data
% Define SAC function
SAC_fun = @(x, f) your_SAC_function(x, f);
% Initial guess for parameters
initial_guess = [0.0001, 0.00015, 1.02, 20000, 0.75];
% Lower and upper bounds for parameters
%lb = [0.00006, 0.0001, 1.0, 10000, 0.65];
%ub = [0.00018, 0.0004, 1.06, 50000, 0.8];
% Optimization using fmincon
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp');
[x_opt, fval, exitflag, output] = fmincon(@(x) objective_fun(x, SAC_fun, SAC_exp, f_exp), initial_guess)%, [], [], [], [], [], [], [], options)
% Display optimization process information
disp(output);
% Display optimized parameter values in command window
disp('Optimized Parameter Values:');
disp(['VL: ', num2str(x_opt(1))]);
disp(['TL: ', num2str(x_opt(2))]);
disp(['T: ', num2str(x_opt(3))]);
disp(['FR: ', num2str(x_opt(4))]);
disp(['P: ', num2str(x_opt(5))]);
disp(['Optimized Objective Value: ', num2str(fval)]);
disp(['Exit Flag: ', num2str(exitflag)]);
disp(['Number of Function Evaluations: ', num2str(output.funcCount)]);
% Plot optimized SAC curve and experimental data
f = 1:1:5500; % Generate frequencies from 1 to 5500 with a step size of 1
SAC_opt = your_SAC_function(x_opt, f); % Compute optimized SAC values
figure;
plot(f, SAC_opt, 'r-', 'LineWidth', 2);
hold on;
plot(f_exp, SAC_exp, 'bo'); % Plot experimental data points
title('Optimized SAC vs. Experimental Data');
xlabel('Frequency');
ylabel('SAC');
legend('Optimized SAC', 'Experimental Data');
grid on;
function SAC = your_SAC_function(x, f)
% Parameters
VL = x(1);
TL = x(2);
T = x(3);
FR = x(4);
P = x(5);
% Constants
AD = 1.2; % Air density
DV = 1.84e-5; % Dynamic viscosity
p0 = 101320; % Air pressure
SH = 1.4; % Specific heat of air
PN = 0.7; % Prandtl constant
c0 = 343; % Air velocity
d = 0.03; % Assuming a constant value for d
% SAC calculation based on provided equations
E3 = sqrt(1 + (4 * 1i * T * T * DV * 2 * pi * f * AD) / (FR * FR * VL * VL * p0 * p0));
E4 = FR * p0 ./ (1i * 2 * pi * f * AD * T);
E1 = T * AD .* (1 + E4 .* E3);
E6 = 1i * 8 * 2 * pi * f ./ (PN * AD * 2 * pi * f * TL * TL);
E7 = sqrt(1 + (1i * PN * AD * 2 * pi * f * TL * TL / 16 * DV));
E8 = 1./(1 - (E6 .* E7));
E2 = (SH * p0) ./ (SH - ((SH - 1).* E8));
KC = 2 * pi * f .* sqrt(E1 ./ E2);
ZC = sqrt(E1 .* E2);
ZS = -1i * ZC .* cot(KC .* d) / P;
R = (ZS - AD * c0) ./ (ZS + AD * c0);
SAC = 1 - abs(R).^2;
end
function error = objective_fun(x, SAC_fun, SAC_exp, f_exp)
% Compute SAC values using the model
SAC_model = SAC_fun(x, f_exp);
% Compute the squared errors at each frequency
squared_errors = (SAC_exp - SAC_model).^2;
% Compute the average squared error
error = sum(squared_errors);
end
3 Comments
Torsten
on 14 Mar 2024
That seems to be the best your can get if the parameter bounds have to be respected. That's why I wrote that you should check whether your model is really appropriate to reflect your simulation results.
More Answers (0)
See Also
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!