Clear Filters
Clear Filters

The issue of optimization (minimization) of the average relative error between experimental and calculated data

126 views (last 30 days)
hello
I want to share the difficulties that I faced. Can someone help
problem statement:
there is a 'x' column where the values ​​of the independent variable are written and there is a 'y' column where the experimental values ​​of the dependent variable are written.
approximation model is considered:
y_calculate=A*x^B+C,
and based on this model, an objective function is created, which is equal to the average value of the relative deviations between y and y_calculate:
error_function = mean(abs(y - y_calculate)) / y)=mean(abs(y - =A*x^B+C)) / y);
Our goal is to select parameters A,B,C in such a way that 'error_function' takes the value of the global minimum.
I calculated the optimal values ​​of A, B, C and got:
A = 85.5880, B = -0.0460, C = 4.8824,
at which error function value for optimized parameters: 0.0285.
but I know in advance the specific values ​​of A, B, C:
A = 1005.6335852931, B = -1.59745963925582, C = 73.54149744754400,
at which error function value for specific parameters: 0.002680472178434,
which is much better than with optimization
Below is the code with visualization, which confirms the above.
clear
close all
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guesses for parameters A, B, C
initial_guess = [1, 1, 1];
% Error function
error_function = @(params) mean(abs(y - (params(1) * x.^params(2) + params(3))) ./ y);
% Optimization of parameters
optimized_params = fminsearch(error_function, initial_guess);
% Results of optimization
A_optimized = optimized_params(1);
B_optimized = optimized_params(2);
C_optimized = optimized_params(3);
% Calculation of the fitted function for optimized parameters
y_calculate_optimized = A_optimized * x.^B_optimized + C_optimized;
% Calculate and display the error function value for optimized parameters
value_error_optimized = error_function(optimized_params);
fprintf('Optimized parameters:\nA = %.4f\nB = %.4f\nC = %.4f\n', A_optimized, B_optimized, C_optimized);
fprintf(' error function value for optimized parameters: %.4f\n', value_error_optimized);
% Other specific parameters A, B, C
A_specific = 1005.63358529310;
B_specific = -1.59745963925582;
C_specific = 73.541497447544;
% Calculation of the fitted function for specific parameters
y_calculate_specific = A_specific * x.^B_specific + C_specific;
% Calculate and display the error function value for specific parameters
value_error_specific = error_function([A_specific, B_specific, C_specific]);
fprintf('Specific parameters:\nA = %.10f\nB = %.14f\nC = %.14f\n', A_specific, B_specific, C_specific);
fprintf(' error function value for specific parameters: %.4f\n', value_error_specific);
% Visualization
figure;
plot(x, y, 'bo-', 'DisplayName', 'Experimental data');
hold on;
plot(x, y_calculate_optimized, 'r--', 'DisplayName', 'Fitted model (Optimized)');
plot(x, y_calculate_specific, 'g-.', 'DisplayName', 'Fitted model (Specific)');
xlabel('x');
ylabel('y');
legend('Location', 'best');
title('Approximation of experimental data');
grid on;
Obviously, my optimization code does not lead to a global minimum of the objective function, since there is a better approximation for specific values ​​of A,B,C. Maybe this is caused by a random selection of the initial values ​​of the parameters A=1, B=1, c=1 and therefore my code is stuck in a local minimum?
who can write a code that will select the A,B,C parameters so as to achieve the global minimum of the target function 'error_function', for any initial iteration data of the variables A,B,C. Thoughts for testing: the value of the target function 'error_function' should not be worse (that is, more) than 0.002680472178434, which is obtained with the specific value of A,B,C: A = 1005.6335852931, B = -1.59745963925582, C = 73.54149744754400

Accepted Answer

Matt J
Matt J on 22 Aug 2024
Edited: Matt J on 5 Sep 2024
If you download minL1lin from
then you can do,
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter B
initial_guess = 1;
% Error function
error_function = @(B) objFun(B, x,y);
% Optimization of parameters
B = fminsearch(error_function, initial_guess);
[fval,ABC]=error_function(B);
[A,B,C]=deal(ABC{:})
A = 1.0460e+03
B = -1.6184
C = 73.5535
fval
fval = 0.0025
function [fval,p]=objFun(theta, x,y)
arguments
theta; x (:,1); y(:,1);
end
B=theta;
d=ones(size(y));
Q=[x.^B, d]./y;
[AC,fval]=minL1lin(Q,d,[],[],[],[],[],[],[],optimoptions('linprog','Display','off'));
fval=fval/numel(y);
A=AC(1); C=AC(2);
p={A, B, C};
end
  39 Comments
roborrr
roborrr on 1 Nov 2024 at 14:08
thank you dear Torshen for correcting the error in my second script - instead of model function y=A.*(1+(lambda*x).^m)+C, there should indeed be function y=A./(1+(lambda*x).^m)+C:
clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter theta
initial_guess = [1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m
lambda0 =0;
lambda1=inf;
m0 =0;
m1 =inf;
A0 = 0;
A1 = inf;
C0 = 0.2*min(y);
C1 = 1.1*min(y);
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0],[lambda1,m1]);
fval=error_function(theta);
fval
function fval=objFun(theta, x,y)
x=x(:);y=y(:); theta=theta(:);
lambda=theta(1); m=theta(2);
% Set A and B to appropriate values !
A = 131.5161945547369;
C = 73.9206126555175;
fval = sum(abs((y - (A./(1 + (lambda * x).^m)+ C))./y));
fval= fval/numel(y);
end
This construct does give the best result - Fval = 0.015, when the best coefficients A and C are known in advance. But I tried to use this construct to optimize all coefficients A, C, lambda and m and got the worst result - Fval = 0.0240. Here is the script I wrote:
clear
% Data
x = [7.3392, 14.6784, 22.0176, 29.3436, 36.6828, 44.0088, 51.3348, 58.674, 66, 73.3392, 80.6652, 88.0044, 95.3304, 102.6696, 109.9956, 117.3348, 124.6608, 132];
y = [115.1079, 87.7698, 80.5755, 78.1611, 76.5743, 75.7074, 74.9375, 74.9453, 74.59, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990, 74.2990];
% Initial guess for parameter theta
initial_guess = [1,1,1,1];
% Error function
error_function = @(theta) objFun(theta, x,y);
% Optimization of parameters
%Bounds on lambda and m
lambda0 =0;
lambda1=inf;
m0 =0;
m1 =inf;
A0 = 0;
A1 = inf;
C0 = 0.2*min(y);
C1 = 1.1*min(y);
%Call the optimizer
theta = fminsearchbnd(error_function, initial_guess,[lambda0,m0,A0,C0],[lambda1,m1,A1,C1]);
fval=error_function(theta);
fval
function fval=objFun(theta, x,y)
x=x(:);y=y(:); theta=theta(:);
lambda=theta(1); m=theta(2);
A=theta(3);
C=theta(4);
% the best values ​​of A and C that should be obtained in this script
% A = 131.5161945547369;
% C = 73.9206126555175;
fval = sum(abs((y - (A./(1 + (lambda * x).^m)+ C))./y));
fval= fval/numel(y);
end
why does the value - Fval = 0.015, obtained from the first script not match the value obtained from the second script - Fval = 0.0240? where am I making a mistake in the second script?
how universal is this construct? is it possible to use it for models of any complexity? or is the class of models limited?
Torsten
Torsten on 1 Nov 2024 at 14:22
Edited: Torsten on 1 Nov 2024 at 14:30
You tried two solution strategies: Solving for all 4 parameters simultaneously (first strategy) and solving for the 2 linear and the 2 nonlinear coefficients separately with result coupling (second strategy). The second strategy turned out to be better and more robust - especially because no initial guesses for A and C are needed (your values A = 1 and C = 1 used for the first strategy most probably are too far off).
That's life in numerical analysis: usually, you have to go the way of trial and error for each problem you want to solve, and generalizations are not possible.

Sign in to comment.

More Answers (0)

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!