Using lsqcurvefit for parameter optimisation

I am replicating the data from a published paper to comprehend the utilization of lsqcurvefit for parameter optimization using best fit to predicted data to target data. The paper provides initial guesses as well as lower and upper bound values. Although the paper used the 'ga' function for parameter optimization, I am attempting to achieve the same results using lsqcurvefit. However, my implementation consistently yields predicted values close to zero. I am seeking guidance to identify any potential issues related to using the lsqcurvefit function in this context. Here I am providng the scripts.

% main_script.m
% Set actual constant parameters
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.04;  % Assuming a constant value for d
P = 0.701;  % Assuming a constant value for P
% Read frequencies and target_SAC values from Excel file
filename = 'MATLABSample_1_SAC.xlsx';
data = readmatrix(filename);
% Extract frequencies and target_SAC from the data matrix
frequencies = data(:, 1);
target_SAC = data(:, 2);
% Define the function to minimize (residuals)
objective_function = @(params,frequencies) calculate_residuals(params, frequencies, target_SAC, AD, DV, p0, SH, PN, c0);
% Initial guess for parameters
initial_guess = [0.0001, 0.00015, 1.02, 20000]; % Initial values for VL, TL, T, FR
% Set lower and upper bounds for parameters
lb = [0.00006, 0.0001, 1, 10000];
ub = [0.00018, 0.0004, 1.06, 50000];
% Optimize using lsqnonlin
optimized_params = lsqcurvefit(objective_function, initial_guess, frequencies, target_SAC, lb, ub);
% Calculate SAC with the optimized parameters
optimized_SAC = calculate_SAC(optimized_params, frequencies, AD, DV, p0, SH, PN, c0);
% Display optimized values
disp('Optimized Values:');
disp(['   VL: ', num2str(optimized_params(1))]);
disp(['   TL: ', num2str(optimized_params(2))]);
disp(['   T: ', num2str(optimized_params(3))]);
disp(['   FR: ', num2str(optimized_params(4))]);
% Plotting
figure;
plot(frequencies, target_SAC, 'o-', 'DisplayName', 'Target SAC');
hold on;
plot(frequencies, optimized_SAC, 'x-', 'DisplayName', 'Optimized SAC');
xlabel('Frequency (Hz)');
ylabel('SAC');
legend('Location', 'best');
title('Target vs Optimized SAC');
grid on;
% Save SAC values to an Excel file in a specific directory
directory = '/Users/pramodya/Documents/MATLAB';
filename = sprintf('%sSample_%d_wrong.xlsx', directory, i);
writematrix([f', SAC'], filename);



% calculate_SAC.m
function SAC = calculate_SAC(params, frequencies, AD, DV, p0, SH, PN, c0, P, d)
% Extract parameters
VL = params(1);
TL = params(2);
T = params(3);
FR = params(4);
% Initialize SAC array
SAC = zeros(size(frequencies));
% Loop through different frequencies
for i = 1:numel(frequencies)
    f = frequencies(i);
    % Evaluate equations for the current set of parameters
    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(i) = 1 - abs(R).^2;
end
end



% calculate_residuals.m
function residuals = calculate_residuals(params, frequencies, target_SAC, AD, DV, p0, SH, PN, c0, P, d)
% Calculate SAC with given parameters
calculated_SAC = calculate_SAC(params, frequencies, AD, DV, p0, SH, PN, c0, P, d);
% Calculate residuals (difference between calculated and target SAC)
residuals = calculated_SAC - target_SAC;
end

Answer 1

[Matt J]

Since you are just trying to verify the results of a previous paper, initialize with the solution parameters that the paper provides, not the initial guess. No reason to make it harder on the optimizer than it has to be.

If nothing else, this will be an important sanity check on the objective function you've provided. If the paper's solution doesn't cause lsqcurvefit to terminate within 1 iteration, you know that your calculate_SAC isn't the same objective as what the paper considered.

[Pramodya]

Thank you for your response. Indeed, I attempted the approach you suggested, and it provided accurate results. This allowed me to confirm the correctness of my SAC model. Essentially, when I utilized the obtained solution parameters and applied them to the calculate_SAC model, the outcomes aligned accurately. However, my current endeavor involves developing a program for parameter optimization in the context of my research. Therefore, I utilized this reference to validate and verify my optimization program. But it gives almost zero optimized SAC results.

[Matt J]

Then clearly your initial guess is inadequate. The test you've done proves it.

[Alan Weiss]

It is possible that there are multiple local optima. You might want to try MultiStart on your problem, as in this example.

Alan Weiss

MATLAB mathematical toolbox documentation

[Pramodya]

Dear Matt and Weiss,

Thank you for your comments. Will try as instructed and let you know the feedback.