SQP fmincon optimization problem

Hi there,

I am trying to set up a code that has an aim to minimze an error function between an experimental dataset and canonical dataset (that is relied on equations and formulas).

I have tried to compile a Sequential Quadratic Programming code by using the fmincon MATLAB function. The main concept is to minimize as said earlier an error function that connects the experimental with the canonical dataset by figuring ou the values of 5 parameters (u_tau, thelta_y, k, pi (Π) and thelta.The MATLAB code is the following:

clear all
clc
%%

% Load experimental data
% data = xlsread('file name.xlsx', 'sheet name', 'Cells');
data = xlsread('Spalding profile_NTUA+Oesterlund Data_yo .xlsx', 'Spald._PD_VP_5cm_4spires.', 'H2:I37');
y_exp = data(:,1);
u_exp = data(:,2);
hold on
plot(y_exp,u_exp)

% Define function to minimize
Cost_Function = @(x) mean(abs((u_plus_Canonical(x,y_exp,u_exp)-u_exp)./u_plus_Canonical(x,y_exp,u_exp)));

% Set initial guess for parameters
x0 = [0.7, -.05, .4, .5, 1.4];

% Set lower and upper bounds for parameters
lb = [0,   -2  ,  .2, 0, 1.2];
ub = [1, -.037 ,  .8, 5, 12];


% Use fmincon with SQP algorithm to optimize the function
options = optimoptions('fmincon','Algorithm','sqp','MaxFunctionEvaluations',10000,'MaxIterations',5000,'OptimalityTolerance',1e-6,'StepTolerance',1e-6, 'Display','iter-detailed');%,'PlotFcn',{@optimplotx,@optimplotfval});
[x_opt, fval] = fmincon(Cost_Function, x0, [], [], [], [], lb, ub, [], options);

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           6    8.650001e-01     0.000e+00     1.000e+00     0.000e+00     1.646e+01  
    1          12    8.252488e-01     0.000e+00     1.000e+00     1.408e-01     1.645e+01  
    2          23    8.251611e-01     0.000e+00     1.681e-01     5.543e-02     1.482e+01  
    3          31    8.248298e-01     0.000e+00     4.900e-01     2.170e-02     6.545e+00  
    4          38    8.243266e-01     0.000e+00     7.000e-01     4.266e-02     2.525e+00  
    5          44    8.240541e-01     0.000e+00     1.000e+00     1.828e-02     1.457e-01  
    6          50    8.228004e-01     0.000e+00     1.000e+00     8.951e-02     1.529e-02  
    7          56    8.211383e-01     0.000e+00     1.000e+00     1.526e-01     1.192e-02  
    8          62    8.210327e-01     0.000e+00     1.000e+00     8.731e-03     3.757e-02  
    9          71    8.210131e-01     0.000e+00     3.430e-01     1.545e-02     7.555e-02  
   10          77    8.209802e-01     0.000e+00     1.000e+00     1.123e-02     3.609e-02  
   11          83    8.209412e-01     0.000e+00     1.000e+00     3.107e-03     3.477e-02  
   12          92    8.209243e-01     0.000e+00     3.430e-01     5.389e-03     8.523e-02  
   13          98    8.209222e-01     0.000e+00     1.000e+00     3.879e-03     3.469e-02  
   14         104    8.209083e-01     0.000e+00     1.000e+00     1.094e-03     3.483e-02  
   15         113    8.208947e-01     0.000e+00     3.430e-01     1.883e-03     8.772e-02  
   16         121    8.208927e-01     0.000e+00     4.900e-01     6.602e-04     3.501e-02  
   17         127    8.208903e-01     0.000e+00     1.000e+00     2.068e-04     3.506e-02  
   18         136    8.208897e-01     0.000e+00     3.430e-01     3.548e-04     8.863e-02  
   19         142    8.208892e-01     0.000e+00     1.000e+00     2.308e-04     3.508e-02  
   20         148    8.208883e-01     0.000e+00     1.000e+00     1.123e-04     3.509e-02  
   21         157    8.208880e-01     0.000e+00     3.430e-01     1.926e-04     8.875e-02  
   22         163    8.208879e-01     0.000e+00     1.000e+00     8.403e-05     3.510e-02  
   23         169    8.208876e-01     0.000e+00     1.000e+00     9.592e-05     3.510e-02  
   24         178    8.208874e-01     0.000e+00     3.430e-01     1.645e-04     8.880e-02  
   25         185    8.208873e-01     0.000e+00     7.000e-01     2.631e-05     3.510e-02  
   26         191    8.208873e-01     0.000e+00     1.000e+00     9.370e-05     3.510e-02  
   27         200    8.208872e-01     0.000e+00     3.430e-01     1.607e-04     8.881e-02  
   28         206    8.208872e-01     0.000e+00     1.000e+00     1.378e-05     3.510e-02  
   29         212    8.208872e-01     0.000e+00     1.000e+00     1.052e-04     3.510e-02  
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
   30         221    8.208872e-01     0.000e+00     3.430e-01     1.803e-04     8.882e-02  
   31         227    8.208872e-01     0.000e+00     1.000e+00     2.771e-05     3.510e-02  
   32         233    8.208871e-01     0.000e+00     1.000e+00     1.605e-04     3.510e-02  
   33         239    8.208871e-01     0.000e+00     1.000e+00     8.019e-04     8.882e-02  
   34         245    8.208871e-01     0.000e+00     1.000e+00     4.768e-04     8.882e-02  
   35         251    8.208868e-01     0.000e+00     1.000e+00     2.375e-03     3.511e-02  
   36         257    8.208851e-01     0.000e+00     1.000e+00     2.235e-02     3.510e-02  
   37         263    8.208782e-01     0.000e+00     1.000e+00     1.046e-01     3.511e-02  
   38         269    8.208663e-01     0.000e+00     1.000e+00     3.044e-01     3.513e-02  
   39         275    8.208618e-01     0.000e+00     1.000e+00     2.714e-01     8.879e-02  
   40         281    8.208606e-01     0.000e+00     1.000e+00     1.575e-01     8.878e-02  
   41         287    8.208599e-01     0.000e+00     1.000e+00     2.021e-01     8.875e-02  
   42         293    8.208597e-01     0.000e+00     1.000e+00     1.153e-01     8.874e-02  
   43         299    8.208596e-01     0.000e+00     1.000e+00     6.958e-02     8.874e-02  
   44         305    8.208596e-01     0.000e+00     1.000e+00     5.699e-02     8.873e-02  
   45         311    8.208594e-01     0.000e+00     1.000e+00     9.363e-02     8.873e-02  
   46         317    8.208591e-01     0.000e+00     1.000e+00     1.544e-01     8.872e-02  
   47         323    8.208582e-01     0.000e+00     1.000e+00     2.920e-01     8.873e-02  
Objective function returned NaN; trying a new point...
   48         330    8.208572e-01     0.000e+00     7.000e-01     1.606e-01     8.874e-02  
Objective function returned NaN; trying a new point...
   49         337    8.208548e-01     0.000e+00     7.000e-01     4.756e-02     3.511e-02  
Objective function returned NaN; trying a new point...
   50         344    8.208545e-01     0.000e+00     7.000e-01     1.462e-02     3.512e-02  
Objective function returned NaN; trying a new point...
   51         360    8.208544e-01     0.000e+00     2.825e-02     2.050e-04     8.883e-02  
Objective function returned NaN; trying a new point...
   52         368    8.208544e-01     0.000e+00     4.900e-01     2.952e-03     3.512e-02  
Objective function returned NaN; trying a new point...
   53         375    8.208544e-01     0.000e+00     7.000e-01     2.165e-03     3.513e-02  
Objective function returned NaN; trying a new point...
   54         383    8.208544e-01     0.000e+00     4.900e-01     4.995e-04     8.883e-02  
Objective function returned NaN; trying a new point...
   55         390    8.208544e-01     0.000e+00     7.000e-01     3.190e-04     3.513e-02  
Objective function returned NaN; trying a new point...
   56         397    8.208544e-01     0.000e+00     7.000e-01     1.146e-04     3.513e-02  
Objective function returned NaN; trying a new point...
   57         406    8.208544e-01     0.000e+00     3.430e-01     5.399e-05     1.270e-02  
Objective function returned NaN; trying a new point...
   58         416    8.208544e-01     0.000e+00     2.401e-01     7.093e-06     2.969e-02  
Objective function returned NaN; trying a new point...
   59         423    8.208544e-01     0.000e+00     7.000e-01     2.961e-05     1.765e-02  
Objective function returned NaN; trying a new point...
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
   60         434    8.208544e-01     0.000e+00     1.681e-01     1.904e-05     1.123e-02  
Objective function returned NaN; trying a new point...
   61         476    8.208544e-01     0.000e+00     3.120e-07     2.677e-06     1.123e-02  

Optimization stopped because the relative changes in all elements of x are
less than options.StepTolerance = 1.000000e-06, and the relative maximum constraint
violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.

plot(y_exp,u_plus_Canonical(x_opt,y_exp,u_exp))

Warning: Imaginary parts of complex X and/or Y arguments ignored.

hold off

%% Optimization Results

% Display the results
disp('Optimal solution:');

Optimal solution:

disp(x_opt);

    1.0000   -0.0370    0.5393    0.0000    3.4201

disp('Optimal cost:');

Optimal cost:

disp(fval);

    0.8209

%%

% Define function for u_plus_Canonical
function u_plus = u_plus_Canonical(x,y,u_exp)
    % Define constants
    v = 1.458*1e-5;
    %k = 0.41;
    B = 5.0;
    s = 0.001093;
    M1 = 30;
    M2 = 2.85;
    a2 = 132.8410;
    a3 = -166.2041;
    a4 = 71.9114;
    
    % Extract variables
    u_tau = x(1);
    thelta_y = x(2);
    k = x(3);
    pi = x(4);
    thelta = x(5);
    
    % Calculate yPlus
    yPlus = (y + thelta_y)*u_tau/v;
    
    % Calculate u_plus_musker
    [yPlus_Musker, uPlus_Musker] = ode45(@(yPlus, uPlus_Musker) (((yPlus.^2)/k) + 1/s) ./ ((yPlus.^3) + ((yPlus.^2)/k) + 1/s), [0, 5000], 0);
    u_plus_musker = interp1(yPlus_Musker, uPlus_Musker, yPlus, 'linear', 'extrap');
    
    % Calculate u_plus_Bump
    uPlusBump = exp(-(log(yPlus/M1).^2)/M2);
    
    % compute ita
    ita = y ./ thelta;

    % compute the numerator and denominator
    numerator = 1 - exp((-1/4)*(5*a2 + 6*a3 + 7*a4)*ita.^4 + (a2*ita.^5) + (a3*ita.^6) + (a4*ita.^7));
    denominator = 1 - exp(-(1/4)*(a2 + 2*a3 + 3*a4));

    % compute w
    W = numerator ./ denominator .* (1 - (log(ita) ./ (2*pi)));

    %
    %data = xlsread('Spalding profile_NTUA+Oesterlund Data_yo .xlsx', 'Spald._PD_VP_5cm_4spires.', 'H2:I37');
    %u_exp = data(:,2);
    % u_plus_e
    u_plus_e = u_exp./u_tau;

    % Calculate u_plus_Canonical
    u_plus = (u_plus_musker + uPlusBump + (2*pi/k)*W).*(y>=0 & y<=thelta) + u_plus_e.*(y>=thelta);
    
end 

The problem lies in the different resutls I get for different lower and upper bounds, that is not supposed to happen because based on the scientific paper my project relies it is clearly stated the following: "Different bounds have been tested without effect either on the result of the method or on the convergence rate."

For the experimental data the code reads only two columns of my Excel ("Spalding profile_NTUA+Oesterlund Data_yo .xlsx") y-the height in m and u-velocity in m/s, nothing more complex than that.

I am thinking that maybe the minimum that is presented through my error function is local and not global?

I know that the results I get are incorrect as the right friction velocity that is reported on the scientific paper for the same data is u_tau=0.7458

Could someone spread some light in the problem because I feel completely stucked?

Thanks

[Pavlos Dionysopoulos]

Did you edit the code?

[Torsten]

Did you edit the code?

Yes, I facilitated the comparison between experimental and fitted curve, I guess.

And the code took too long to run under MATLAB online when you read the Excel file each time your model function was called. So I passed u_exp to the function in order to run the code here.

I wonder why u_plus depends on u_exp:

u_plus = (u_plus_musker + uPlusBump + (2*pi/k)*W).*(y>=0 & y<=thelta) + u_exp./u_tau.*(y>=thelta);

[Pavlos Dionysopoulos]

ignore this term because it is always y<=thelta

[Torsten]

How does it happen that your function becomes complex-valued ? Maybe the log(ita) ? You must prevent this.

[Pavlos Dionysopoulos]

Is this the problem of my optimization?

[Torsten]

It's at least one problem of your code that must be solved. There might be other issues - I can't tell in advance.

If you check the output starting from iteration 42 to the end, it says that your objective function returned NaN values. You will have to find the reason why this happens.