Paramter optimization using fmincon
3 views (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
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!