Cure fitting - lsqnonlin

Hi all,

I have been attempting to develop a code to determine optimized values for four unknown parameters. I utilized experimental data for curve fitting, aiming for the modeled data (SAC) to best match the experimental data (target_SAC) using optimized values. However, I consistently observe differences between the modeled and experimental values. To investigate this issue, I initially inputted my guess, lower bound, and upper bound with the true optimized values obtained from a published paper. Surprisingly, the program correctly produced the optimized values (as expected, given the correct inputs). Nonetheless, the modeled values still differ from the experimental ones. When plotted, the two graphs exhibit discrepancies, even though the optimized values should theoretically yield identical modeled values. I am unable to pinpoint the source of this issue. Could you please assist me with this? I will provide the three files, all of which run without any errors.

%MAIN FILE
% Additional constants
% 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.03;  % 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 = 'Target_SAC.xlsx';
data = readmatrix(filename);

Error using readmatrix
Unable to find or open 'Target_SAC.xlsx'. Check the path and filename or file permissions.

% 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) calculate_residuals(params, frequencies, target_SAC, AD, DV, p0, SH, PN, c0,P,d);
% Initial guess for parameters
initial_guess = [0.0000845, 0.0001023, 1.1, 25543]; % Initial values for VL, TL, T, FR
% Set lower and upper bounds for parameters
lb = [0.0000845, 0.0001023, 1.1, 25543];
ub = [0.0000845, 0.0001023, 1.1, 25543];
% Optimize using lsqnonlin
optimized_params = lsqnonlin(objective_function, initial_guess, lb, ub);
% 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))]);
% Calculate SAC with the optimized parameters
optimized_SAC = calculate_SAC(optimized_params, frequencies, AD, DV, p0, SH, PN, c0,P,d); 
% Create figure for target SAC
figure;
subplot(2,1,1); % Create a 2x1 grid of plots and select the first subplot
plot(frequencies, target_SAC, 'o-', 'DisplayName', 'Target SAC');
xlabel('Frequency (Hz)');
ylabel('SAC');
legend('Location', 'best');
title('Target SAC');
grid on;
% Create figure for optimized SAC
subplot(2,1,2); % Select the second subplot
plot(frequencies, optimized_SAC, 'x-', 'DisplayName', 'Optimized SAC');
xlabel('Frequency (Hz)');
ylabel('SAC');
legend('Location', 'best');
title('Optimized SAC');
grid on;


%Function file
function SAC = calculate_SAC(params, frequencies, AD, DV, p0, SH, PN, c0, d, P)
% 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 frequency
    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

%residual file
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]

I initially inputted my guess, lower bound, and upper bound with the true optimized values obtained from a published paper. Surprisingly, the program correctly produced the optimized values (as expected, given the correct inputs).

That's not surprising. If lb=ub there is only 1 feasible point, which of course lsqnonlin will return. You should instead use whatever bounds the earlier paper used.

If the previously published values are not fitting the data, then the possible reasons are,

1. You have a bug in your implementation of the model, so that it disagrees with the previous author's implementation.
2. The previous paper's work has an error, either in the reporting of the model, or in its experimental implementation.
3. The model you've implemented agrees with the earlier paper, but you have a bug in the way your data was acquired so that it does not suit the model.

[Pramodya]

Thank you matt!!Will check the three options and come up with the outcome. Thank you so much!

[Matt J]

You're welcome, but if your question is resolved please Accept-click the answer.