Paramter optimization using fmincon

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

[Torsten]

You forgot to include "Target_SAC.xlsx".

[Pramodya]

Hi Torsten

Thank you for the reply. Here is the file.

Answer 1

[Torsten]

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)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

x_opt = 1x5

1.0e+04 * 0.0008 0.0015 8.4840 1.9942 0.0163

fval = 130.5470

exitflag = 1

output = struct with fields:
         iterations: 47
          funcCount: 313
    constrviolation: 0
           stepsize: 0.0015
          algorithm: 'interior-point'
      firstorderopt: 5.1329e-07
       cgiterations: 21
            message: 'Local minimum found that satisfies the constraints....'
       bestfeasible: [1x1 struct]

% Display optimization process information
disp(output);

         iterations: 47
          funcCount: 313
    constrviolation: 0
           stepsize: 0.0015
          algorithm: 'interior-point'
      firstorderopt: 5.1329e-07
       cgiterations: 21
            message: 'Local minimum found that satisfies the constraints....'
       bestfeasible: [1x1 struct]

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

Optimized Parameter Values:

disp(['VL: ', num2str(x_opt(1))]);

VL: 8.0173

disp(['TL: ', num2str(x_opt(2))]);

TL: 14.8636

disp(['T: ', num2str(x_opt(3))]);

T: 84840.3064

disp(['FR: ', num2str(x_opt(4))]);

FR: 19942.1905

disp(['P: ', num2str(x_opt(5))]);

P: 163.2276

disp(['Optimized Objective Value: ', num2str(fval)]);

Optimized Objective Value: 130.547

disp(['Exit Flag: ', num2str(exitflag)]);

Exit Flag: 1

disp(['Number of Function Evaluations: ', num2str(output.funcCount)]);

Number of Function Evaluations: 313

% 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

[Pramodya]

Hi Torsten,

Thank you very much for your kind feedback. However, I didnt understand for not mentioning lower and upper bounds so that is has given abnormal values. But as per the given code, when bounds are also mention, the output is something like this.

Iter Func-count Fval Feasibility Step Length Norm of First-order

step optimality

0 6 3.434672e+03 0.000e+00 1.000e+00 0.000e+00 8.933e+01

1 12 3.430209e+03 0.000e+00 1.000e+00 5.002e-02 6.363e-02

2 18 3.430208e+03 0.000e+00 1.000e+00 1.850e-03 1.784e-03

3 24 3.430208e+03 0.000e+00 1.000e+00 9.413e-03 1.784e-03

4 30 3.430208e+03 0.000e+00 1.000e+00 4.705e-02 1.784e-03

5 36 3.430207e+03 0.000e+00 1.000e+00 9.220e-02 1.784e-03

6 42 3.430207e+03 0.000e+00 1.000e+00 3.066e-02 1.784e-03

7 48 3.430207e+03 0.000e+00 1.000e+00 1.533e-01 1.784e-03

8 54 3.430206e+03 0.000e+00 1.000e+00 7.663e-01 2.152e-03

9 60 3.430199e+03 0.000e+00 1.000e+00 3.832e+00 6.387e-03

10 66 3.430165e+03 0.000e+00 1.000e+00 1.916e+01 7.145e-03

11 72 3.429993e+03 0.000e+00 1.000e+00 9.596e+01 7.167e-03

12 78 3.429106e+03 0.000e+00 1.000e+00 4.832e+02 2.437e-02

13 84 3.423920e+03 0.000e+00 1.000e+00 2.506e+03 7.462e-01

14 90 3.400767e+03 0.000e+00 1.000e+00 6.890e+03 1.148e+01

15 96 3.400767e+03 0.000e+00 1.000e+00 3.056e-05 9.153e-05

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in

feasible directions, to within the value of the optimality tolerance,

and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

iterations: 15

funcCount: 96

algorithm: 'sqp'

message: 'Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 7.316143e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'

constrviolation: 0

stepsize: 3.0561e-05

lssteplength: 1

firstorderopt: 9.1534e-05

bestfeasible: [1×1 struct]

Optimized Parameter Values:

VL: 0.00014944

TL: 0.0004

T: 1.06

FR: 10000

P: 0.8

Optimized Objective Value: 3400.7674

Exit Flag: 1

Number of Function Evaluations: 96

And as you mentioned, The curves look qualitatively different. It is really appreciated if you could provide with a solution or any other paramter optimization method to use in MATLAB for the same issue?

[Torsten]

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.

[Pramodya]

Yes Torsten. Thank you for the feedback.