Run Nested FOR-loop Parallelly for Multivariable Function Optimization

F_max = 0; % Temp var for max of F
F_curr = 0; % Temp var for current F

for x = -0.02:0.001:0.02
    for x_1 = 20:100
        for x_2 = 20:100
            for x_3 = 20:100
                for x_4 = 20:100
                    for x_5 = -15:15
                        for x_6 = -15:15
                            for x_7 = -15:15
                                for x_8 = -15:15
                                    
                                    F_curr = double(F(x,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8));
                                    if F_curr>F_max
                                        F_max = F_curr;
                                    end
                                    
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

F is a (symbolic) function of 9 (symbolic) variables x, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8. The formula of F is given in the following Link1.

I have to find maximum value F takes for bounded-constraints of the variables given in Link2.

How do I edit this code so that it can run parallelly? The reason I can't apply the usual solution is because the variables aren't looping by integer iterations starting from 0 (e.g. i = 1:n (some integer)).

[Sam Marshalik]

Hey Joshua, I do not currently have an opportunity to play around with the code, but you will want to employ parfor to speed this up. The issue is that F_max is a temporary variable on all of the workers and they will not be able to exchange information to determine who has the highest F_max (parfor workers are not able to communicate with one another).

I think something you can try is making F_max a sliced output variable (https://www.mathworks.com/help/parallel-computing/sliced-variable.html#bq_tiga) - this will give you a large array with all of the values from your loops. You can then determine the highest value (max(F_max)) from the entire list.

You can also do the following to deal with the outermost parfor-loop, since the non-integers will be an issue:

x = 0:0.1:1;

parfor idx = 1:length(x)

disp(x(idx))

end

Answer 1

[Matt J]

I don't think a loop over all 9 variables is going to be practical (10^15 combinations).

An important observation, though, is that your function F() is linear with respect to x5,x6,x7,x8,x9. This means that the maximum will be achieved at one of the extreme values of these variables. Consequently, you don't have to search -15:15. You only have to search the two end points -15 and 15 for each of these variables for a total of 16 combinations. For each of the 16 combinations, you need to do a grid search over X, X_1,X_2,X_3,X_4 but the dimension of that search (EDIT:) can be done with increased vectorization as follows.

Xgrid ={  -0.02:.001:0.02;  
          [-15,15]};

 Xgrid=  Xgrid([1,2,2,2,2]);
 
 
 sizeT=repelem(numel(20:100),1,4);
      
[X_1,X_2,X_3,X_4]=ndgridVecs(20:100); % ndgridVecs available at https://www.mathworks.com/matlabcentral/fileexchange/74956-ndgridvecs

[X_0, X_5,X_6,X_7,X_8]=ndgrid(Xgrid{:});

Fdouble=@(x1,x2,x3,x4,x5,x6,x7,x8,x9) -(x1.^2+x2.^2+x3.^2+x4.^2+x5.^2+x6.^2+x7.^2+x8.^2+x9.^2); %example function
%Fdouble=matlabFunction(F); %true function

N=numel(X_0);


F_max=nan(1,N);
loc=cell(1,N);

tic
parfor n=1:numel(X_0)
                
    [x_0, x_5,x_6,x_7,x_8] = deal(X_0(n), X_5(n),X_6(n),X_7(n), X_8(n));
    
    
                T=Fdouble(x_0,X_1,X_2,X_3,X_4,   x_5,x_6,x_7,x_8);
                   
                [F_max(n),I]=max(T,[],'all','linear');
            
                loc{n}=[x_0,I,   x_5,x_6,x_7,x_8]; %location of maximum


end

[F_max,nmax]=max(F_max);

loc=loc{nmax};

[j,k,l,m]=ind2sub(sizeT,loc(2));

loc=[loc(1), X_1(j),X_2(k),X_3(l),X_4(m) ,loc(3:end)];

toc%Elapsed time is 46.687547 seconds.

[Matt J]

As you told, maximum occurs at [x5,x6,x7,x8] = [-15,15,15,-15].

That's stronger than what I said. I said there were 16 possible combinations of [x5,x6,x7,x8] that needed to be considered.

[Matt J]

The reason for stating that is because, the coefficients of [x5,x6,x7,x8] in the function is [-1,1,1,-1]

The coefficients that I see are complicated functions of x,x1..x4. I don't see how you can anticipate their signs.

I've modified my answer above, however, and tested that it runs with more modest memory consumption.

[Matt J]

And coefficient of 1st fraction is +1

No, it isn't. The coefficient of x7 in the first fraction is,

which means that if is negative at the optimum, then the whole coefficeint will be negative and the value of x7 that will maximize the function is -15, not +15.

Answer 2

[Matt J]

Xgrid ={ -0.02:0.001:0.02;
    20:100;
    20:100;
    20:100;
    20:100;
    -15:15;
    -15:15;
    -15:15;
    -15:15};

sz=cellfun(@numel,Xgrid);
N=prod(sz);
J=numel(sz);

Fdouble=matlabFunction(F);

F_max=-inf;

parfor n=1:N

   sub=cell(J,1);
   
   [sub{1:J}]=ind2sub(sz,n); %convert to subscripts

   X=cellfun(@(A,B) A(B), Xgrid,sub,'uni',0);  %lookup grid values

   F_max=max(F_max, Fdouble(X{:}) );  %reduction

end
                  

[Matt J]

The formula of F is given in the following Link1.

The formula shown at your link has a number of recurring intermediate quantities. It would therefore improve performance if you used matlabFuntion's code optimization options, like in this example.

https://www.mathworks.com/help/symbolic/matlabfunction.html#mw_d01438db-f7c5-471a-b3d7-b0aaa7eb67ba

This would mean converting F to an mfile instead of to an anonymous function.

[Joshua Roy Palathinkal]

Seems like I did something similar. However, I could only extend it to perform local optimization. Each time I change the starting point, I am getting different results.

I am attaching the codes for reference:

Optimization_OPTIM.m:

clc, clear,

prob = optimproblem('ObjectiveSense','max');

x = optimvar("x","LowerBound",-0.02,"UpperBound",0.02);
x_1 = optimvar("x_1","LowerBound",20,"UpperBound",100);
x_2 = optimvar("x_2","LowerBound",20,"UpperBound",100);
x_3 = optimvar("x_3","LowerBound",20,"UpperBound",100);
x_4 = optimvar("x_4","LowerBound",20,"UpperBound",100);
x_5 = optimvar("x_5","LowerBound",-15,"UpperBound",15);
x_6 = optimvar("x_6","LowerBound",-15,"UpperBound",15);
x_7 = optimvar("x_7","LowerBound",-15,"UpperBound",15);
x_8 = optimvar("x_8","LowerBound",-15,"UpperBound",15);

obj = fcn2optimexpr(@objfun,x,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8);

prob.Objective = obj;

x0.x = 0;
x0.x_1 = 20;
x0.x_2 = 20;
x0.x_3 = 20;
x0.x_4 = 20;
x0.x_5 = -15;
x0.x_6 = -15;
x0.x_7 = -15;
x0.x_8 = -15;

[sol,fval,exitflag] = solve(prob,x0)

objfun.m

function f = objfun(x,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8)

    f = (x_1*x_2*x_7*(x_2/2 + 1/50)^2*(x_1 - 2*x + 2*x_2 + 1/25))/((x_2/2 + 1/50)^2 + (x_1/2 - x + x_2 + 1/50)^2)^(5/2) - (x_3*x_4*x_5*(x_4/2 + 1/50)^2*(2*x + 2*x_1 + 2*x_2 + x_3 + 1/25))/((x_4/2 + 1/50)^2 + (x + x_1 + x_2 + x_3/2 + 1/50)^2)^(5/2) - (x_3*x_4*x_8*(x_4/2 + 1/50)^2*(2*x_1 - 2*x + 2*x_2 + x_3 + 1/25))/((x_1 - x + x_2 + x_3/2 + 1/50)^2 + (x_4/2 + 1/50)^2)^(5/2) + (x_1*x_2*x_6*(x_2/2 + 1/50)^2*(2*x + x_1 + 2*x_2 + 1/25))/((x + x_1/2 + x_2 + 1/50)^2 + (x_2/2 + 1/50)^2)^(5/2);

end

Thank you soo much for your inputs. It is very helpful. Much appreciated.

[Joshua Roy Palathinkal]

Could you also help me in figuring/storing the variable values at which the function is mazimum? That is the objective of this problem - "Where does the function attain the global maximum?" Issue is that even when we use parfor, the maximum would be different for different workers. Probably one way to go around this is that I can store variable values for each worker (where maximum is obtained for each workers); and then take the maximum. But how do I do that?

Sorry for the trouble.