fmincon with and without parallel computing yields different results

Hi all, I have a question concerning the use of fmincon (nonlinear constraints) with the parallelisation toolbox. I use the following code:

options=optimset('Display','iter','Algorithm','active-set','UseParallel','always');

[man_var,f,eflag,outpt] = fmincon(@calculate_obj,man_var_0,[],[],[],[],man_var_lb,man_var_ub,@calculate_ineqc,options);

In the functions calculate_obj and calculate_ineqc I evalute a process model and use its results to calculate the nonlinear constraints and objective value. My problem is that when I add 'UseParallel','always'in the options, I get a different result from fmincon and it seems that fmincon with 'UseParallel','always' solves the problem differntly compared to the fmincon without parallelisation. Here, the problem appers that in the case with parallel computing fmincon starts to oszillate and does not find a solution. As fare as I understood the parallelisation toolbox, it should only distribute the computation to different workes, but the final solution should remain the same. I'm currently using version 2013b, because the model was written in this version of Matlab.

It would be nice if someone could give me a hint how to solve this problem. Thanks a lot.

Best regards Carsten

Answer 1

[Alan Weiss]

The only difference between using parallel computing and serial computing is that, in parallel, finite difference approximations to the gradient are done in parallel.

So I believe that your calculate_obj function or calculate_ineqc function use random numbers in some way, or write to files that get overwritten in parallel, or some such thing. I can guarantee that the underlying fmincon algorithm is 100% identical in parallel and serial.

For caveats and similar issues, see Improving Performance with Parallel Computing, particularly the section on Factors That Affect Results.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Carsten Asmanoglo]

Hi Alan, thank you for your fast response. I exactly thouhgt the same, but I really don't know where this could happen because I do not use random numbers anywhere in my underlying process model and I also don't safe /overwrite anything in the current Folder. As calculate_obj(man_var) and calculate_ineqc(man_var) are seperate functions, they should work in their own Workspace and only depend on the input variable(man_var) and therefore should not exchange/change/overwrite any varible (I do not use global variables). Do you have any suggenstions due to what effect this "overwrite" might could have happened?

These are my functions:

function perform_optimization
clear all
close all
clc
%Remakrs
% 1) define man_var_0 (this function)
% 1a)define boundaries of man_var  (this function)
% 2) grey out param.... which will be applied in optimization problem in calculate_obj_con(man_var)
% 3) define: param.X=man_variable(i) in calculate_obj_con(man_var)
% 4) Check constraints and obj function in calculate_obj_con

% Additional for embed function

%%Definition of optimization problem and call of algorithsm
man_var_0=[11];
man_var_lb=[9];
man_var_ub=[14];

options=optimset('Display','iter','Algorithm','active-set');%,'UseParallel','always');

[man_var,f,eflag,outpt] = fmincon(@calculate_obj,man_var_0,[],[],[],[],man_var_lb,man_var_ub,@calculate_ineqc,options);

save('optiresults_test')
end

Calculate_obj

if true
  function obj=calculate_obj(man_var)       

            [obj,ineqc]=calculate_obj_con(man_var);

    end
  end

calculate_ineqc

      function [ineqc,eqc]=calculate_ineqc(man_var)

                  [obj,ineqc]=calculate_obj_con(man_var);

          eqc=[];

      end

calculate_obj_con

function [obj,ineqc]=calculate_obj_con(man_var)

%%basic setting
    param.L_B_x_ad=5;                                          % /m length of adsorption region
    param.L_B_x_heat=6;                                        % /m length of heating region
    param.L_B_x_q_ex=6;                                        % /m length of q_ext region
    param.L_B_x_de=6;                                          % /m length of q_ext region
    param.L_B_z=5*10^-3;                                       % /m thickness of band in normal direction
    param.L_G_z_lb=2*10^-2;                                    % /m thickness of gas region normal to the band (heating region)
    param.u_G_in=-0.004;                                       % /m s^-1      Gas velocity at hydrogen inlet
    param.u_B_value=0.001;                                     % /m s^-1
%     param.q_heat_ex_heating_value=11;                          % W m^-2 heat input in heating section

%%definition of manipulated variables
    param.q_heat_ex_heating_value=man_var(1);

%%Calculation of obj function and constraint via the model

      opti_results=solve_problem(param);

      obj=-opti_results.mass_Product_CO2_out_year;
      ineqc=[0.10-opti_results.y_Product_CO2_out;...
             opti_results.T_max-373.15;...
             1-opti_results.q_heat_ratio]; 
  end

In solve_problem the process model is solved by means of method of lines (ode15s a for integration over time and a finite volume scheme for the spatial discretization. As the functions do not interfere with each other (at least I have know idea where this should come from ;) ) I really do not know where to search for an error.

Best regards Casten

[Carsten Asmanoglo]

Hi Alan,

I have again checked my code and tried to find out when/where the problem appears. It seems that the use of ode15s leads to slight deviations in the solution. At time step 27 a deviation in the order of 10^-30 appears which seems to self-accelerate in following steps of the odefunction. I doublechecked the input parameters to the odefunction to make sure that they are exactly the same. As I'm very sure that "my" functions called during the run of ode15s do not load any variables, I have really no idea why this deviation appers.

If you have any idea further ideas, please let me know.

Best regards,

Carsten

[Alan Weiss]

It looks like your problem has just one control variable, man_var, and a lower bound of 9 and an upper bound of 14. Have you tried plotting the objective and max of the nonlinear constraints over this region, even with a sparse set of points?

You would undoubtedly benefit from using the technique in Objective and Nonlinear Constraints in the Same Function. And you might want to set the finite differences to larger than default; see Optimizing a Simulation or ODE.

As for your problem running in parallel, I don't understand it any better than before, sorry. Something is going on with your function, but I don't know what.

Alan Weiss

MATLAB mathematical toolbox documentation

[Carsten Asmanoglo]

Hi Alan,

thanks again for your fast answer. Yes the problem has (currently) only one manipulated variable, due to debugging purpose. However, I would like to do an optimization with about 7-10 manipulated variables.

I know that I can embed the functions to fasten the code, but during debugging I'm not using this option.

I have now clearly located the problem. It is caused during the use of ode15s. I compared the exchange of data between ode15s and my function (output to ode 15s of my function which calculates the differential for ode15s and input state variables form ode15s to my function) for the case with and without parallelisation. Do you have any idea, why with the same output of my function in timestep x ode15s at timestep x+1 sets another input to my function? (The input output behavior before timestep x is excalty the same).

Best regards Carsten

[Alan Weiss]

Sorry, I really don't know what is going on. The only idea I have is that you might be running in parallel on a network of computers, rather than a multicore machine, and the miniscule error comes from an answer being off by machine eps, which is to be expected whenever you run software on different machines or software versions. These kinds of errors can quickly expand on some ODEs, as the well-known chaos examples show.

If I am right about this, then the answer in parallel is not wrong, any more than the answer in serial is wrong. It just means that your problem is sensitive to miniscule differences in computations, which is not a fault of the solution computations, but is an inherent feature of your problem.

Sorry, that's all I know.

Alan Weiss

MATLAB mathematical toolbox documentation

[Lipe]

I would like to report a similar problem (R2024b), and I am not on a network of computers.

The optimization halts at an infeasible point after iteration 0 or 1 (see comments next to the algorithms).

The problem is persistent regardless of optimization algorithm.

The output of fmincon is not used to construct the problem and the order with which I run with or without UseParallel does not change the outcome. It is always solved when false, and it always returns infeasible when true.

clc

obj.optimizer_options.Algorithm = 'interior-point'; % Func-count per iter: 263
obj.optimizer_options.Algorithm = 'sqp';            % Func-count per iter: 263 264
obj.optimizer_options.Algorithm = 'sqp-legacy';     % Func-count per iter: 263 264
obj.optimizer_options.Algorithm = 'active-set';     % Func-count per iter: 263
% Obs: There are 263 decision variables.

%%
disp('-----------------------------------------------------')
disp('TRUE')
disp('-----------------------------------------------------')
obj.optimizer_options.UseParallel = true;

[wopt_0,fval,exitflag,output] = fmincon(@(w) obj.objfun(w, u_init),w0,...
                    [],[],[],[],obj.wL,obj.wU,...
                    @(w) obj.confun(w, x_init),...
                    obj.optimizer_options);

%%
disp('-----------------------------------------------------')
disp('FALSE')
disp('-----------------------------------------------------')
obj.optimizer_options.UseParallel = false;


[wopt_1,fval,exitflag,output] = fmincon(@(w) obj.objfun(w, u_init),w0,...
                    [],[],[],[],obj.wL,obj.wU,...
                    @(w) obj.confun(w, x_init),...
                    obj.optimizer_options);

%%
disp('-----------------------------------------------------')
disp('TRUE')
disp('-----------------------------------------------------')
obj.optimizer_options.UseParallel = true;

[wopt_2,fval,exitflag,output] = fmincon(@(w) obj.objfun(w, u_init),w0,...
                    [],[],[],[],obj.wL,obj.wU,...
                    @(w) obj.confun(w, x_init),...
                    obj.optimizer_options);

%%
disp('-----------------------------------------------------')
disp('FALSE')
disp('-----------------------------------------------------')
obj.optimizer_options.UseParallel = false;

[wopt_3,fval,exitflag,output] = fmincon(@(w) obj.objfun(w, u_init),w0,...
                    [],[],[],[],obj.wL,obj.wU,...
                    @(w) obj.confun(w, x_init),...
                    obj.optimizer_options);