How to get a surface plot for the given function to know how many minima are there?

I have a function whose code is given below. I want to determine whether it has one local minimum or several local minima? Fo this I want to see its surface plot But how? In this function, there are four decision variables namely A1, A2 ,Theta1 and Theta2. All are given in a row vector u. The range of A1 or A2 is from 0 to 10 and the range of Theta1 or Theta2 is from -90 to +90. How can I get its surface plot?

function e=myfun(b)% b is a rand vector of the same size as is u below
u=[2 7 50 85];% u=[A1 A2 Theta Theta2] 
C=numel(b);                 
P=C/2;
M=2*C;
xo=zeros(1,M);
    for k=1:M 
        for i=1:P 
            xo(1,k)=xo(1,k)+u(i)*exp(-1i*(k-1)*pi*cosd(u(P+i))); 
        end 
    end

xe=zeros(1,M);
    for k=1:M
            for i=1:P
                xe(1,k)=xe(1,k)+b(i)*exp(-1i*(k-1)*pi*cosd(b(P+i)));
            end 
    end 
    
    abc=0.0;
    for m1=1:M
      abc=abc+(abs(xo(1,m1)-xe(1,m1))).^2;  
    end 
    e=abc/M; 
end

[Torsten]

You cannot use a surface plot for a function with four independent variables. Only a function of two independent variables can be visualized this way.

[Sadiq Akbar]

@Torsten

Thanks a lot dear Torsten for your kind reply. So what type of plot can we get for this so that I can know the number of local minima with in the search space i.e., 0 to 10 for A and -90 to 90 for Theta?

[Torsten]

No way to visualize an objective function that depends on more than two unknowns.

[Sadiq Akbar]

@Torsten

Ok if we keep A1 and A2 fixed and keep Theta1 and Theta2 as varying from -90 to +90, then how will we get the surface plot?

[Torsten]

theta1 = -90:1:90;
theta2 = -90:1:90;
emin = Inf;
for i = 1:numel(theta1)
    for j = 1:numel(theta2)
        e(i,j) = myfun([2 7 theta1(i) theta2(j)]);
        if e(i,j) < emin
          imin = i;
          jmin = j;
          emin = e(i,j);
        end
    end
end
theta1(imin)

ans = -50

theta2(jmin)

ans = -85

emin

emin = 0

surf(theta1,theta2,e.')

function e=myfun(b)% b is a rand vector of the same size as is u below
u=[2 7 50 85];% u=[A1 A2 Theta Theta2] 
C=numel(b);                 
P=C/2;
M=2*C;
xo=zeros(1,M);
    for k=1:M 
        for i=1:P 
            xo(1,k)=xo(1,k)+u(i)*exp(-1i*(k-1)*pi*cosd(u(P+i))); 
        end 
    end
xe=zeros(1,M);
    for k=1:M
            for i=1:P
                xe(1,k)=xe(1,k)+b(i)*exp(-1i*(k-1)*pi*cosd(b(P+i)));
            end 
    end 
    
    abc=0.0;
    for m1=1:M
      abc=abc+(abs(xo(1,m1)-xe(1,m1))).^2;  
    end 
    e=abc/M; 
end

[Sadiq Akbar]

@Torsten

Thanks a lot for your kind response. Yes, it works now. But what is the role of -50 and -85 and emin-=0 here in the code which is displayed after execution? Further, if I keep different values of Theta1 and Theta2 in vector u, then the shape of the surface plot becomes different, why is it so? Additionally do the peaks and depths tell me about the maxima and minima within the given range i.e. -90 to 90?

[Torsten]

But what is the role of -50 and -85 and emin-=0 here in the code which is displayed after execution?

theta1(imin) and theta2(jmin) are the mesh coordinates for theta1 and theta2 that give the minimum value of the objective function given that A1 = 2 and A2 = 7.

Further, if I keep different values of Theta1 and Theta2 in vector u, then the shape of the surface plot becomes different, why is it so?

Because u is part of the objective function:

xo(1,k)=xo(1,k)+u(i)*exp(-1i*(k-1)*pi*cosd(u(P+i)));

Changes in u lead to changes in e.

Additionally do the peaks and depths tell me about the maxima and minima within the given range i.e. -90 to 90?

Yes, assuming that A1 = 2 and A2 = 7.

Answer 1

[John D'Errico]

You cannot plot 4 independent variables. In fact, you would need 5 dimensions to plot, since you want to plot the result of that function for each of the 4 independent variables. But your monitor can display only 2 dimensions. Surface plots work as projections into the 2-d plane of the monitor, because your brain is wired to understand what it sees as a surface. Even then, your brain can be confused, and there are many optical illusions that prove this fact.

However, you cannot plot a 5 dimensional thing. It does not matter that you want to do so, or that your teacher told you to do so for 1 or 2 unknown variables. They apparently did not explain to you that this idea fails once you go beyond that point.

Wanting is not enough. Perhaps you have heard the old phrase, "If wishes were horses, beggars would ride"? That applies here. What works in a lower number of dimensions does not always extend to higher dimensions.

Is there anything you can do? In general, for some completely arbitrary nonlinear function one cannot know how many extrema there will be. If it is simple enough of course, then you may be able to derive a result mathematically. But there is no mathematical trickery you can employ that will always work for some fully general function.

One approach (often called multi-start) will be to start a nonlinear optimizer at very, many distinct start points. Allow it to operate until convergence for each start point. Then collect all of the results, and cluster all of the solutions, since you will get subtly different results even for the same extrema due to convergence tolerances. Count the number of substantially distinct solutions. That will be at least an approximation to the number of distinct solutions, unless you missed some due to insufficient start points in the multi-start.

The problem even with multi-start is that is is a computer hungry solution. a 4 dimensional search space can be pretty large, depending on how complex is your function. Not too bad usually in only 4-d, but I'm not going to spend the CPU time to see how complicated is your function.

[Sadiq Akbar]

@Torsten

@John D'Errico

Ok if we keep A1 and A2 fixed and keep Theta1 and Theta2 as varying from -90 to +90, then how will we get the surface plot?

[John D'Errico]

That is now trivial. Fix A1 and A2 at some specific numerical values. That is to say, not as symbolic unknown constants, but as numbers. Choose some values. Then use meshgrid to generate a set of values for Theta1 and Theta2. Compute the function for each combination of those 4 parameters, and plot, using surf. What is the problem? Even if you need to use a loop to evaluate all of those combinations, WTP?

Of course, this won't tell you how many minima there are, since the true minima probably do not have those chosen values of A1 and A2. It tells you absolutely nothing about the minima, their locations, or the number of them.

As I suggested, if that information is your goal, then the use of a multi-start method is probably the best approach.

[Sadiq Akbar]

@John D'Errico

What do you mean by a multi-stat method and how will we do it?

[Torsten]

help MultiStart

 MultiStart A multi-start global optimization solver. 
    A MultiStart object is a solver that attempts to locate multiple local
    solutions (possibly the global one among them) to a given problem by
    starting from different start points. Create this object by typing its
    constructor, MultiStart and supplying parameter values (if different
    from default). 
 
    MS = MultiStart constructs a new multi-start optimization solver
    with its properties set to defaults.
 
    MS = MultiStart(PROP, VAL, ...) specifies a set of property-value
    pairs that are applied to the multi-start optimization solver before
    creating it.
 
    MS = MultiStart(OLDMS, PROP, VAL, ...) creates a copy of the
    MultiStart solver OLDMS with the named properties altered with the
    specified values.
 
    MS = MultiStart(GS) constructs a new MultiStart solver and copies
    the common parameter values in the GlobalSearch solver GS into the new
    solver MS.
 
    MultiStart properties:
        Display             - detail of display
        FunctionTolerance   - minimum distance between two separate 
                              objective function values
        XTolerance          - minimum distance between two separate points 
        MaxTime             - time allowed to run the solver
        UseParallel         - perform calculation in parallel mode
        StartPointsToRun    - points that the solver will start from
        OutputFcn           - user-defined output functions
        PlotFcn             - functions to plot solver progress
 
    MultiStart method:
        run                 - run a local solver from multiple start points 
                              for a given optimization problem
 
 
    Typical workflow to run the MultiStart solver:
    ==============================================
    1. Set up the PROBLEM structure
        PROBLEM = createOptimProblem('fmincon','objective',...)
    2. Construct the MultiStart solver
        MS = MultiStart
    3. Run the solver from NUMPOINTS random start points
        run(MS,PROBLEM,NUMPOINTS)
 
    Example:
       Run fmincon from 20 random start points for the optimization problem 
          minimize x(1)^2 + 4*sin(5*x(2)); subject to 
                    (x(1)-1)^2 + (x(2)-1)^2 <= 25, 
                    -5 <= x(1) <= 5 and
                    -5 <= x(2) <= 5.
 
       Specify the first constraint in a MATLAB file function such as
          function [c,ceq] = mycon(x)
          c = (x(1)-1)^2 + (x(2)-1)^2 - 25;
          ceq = [];
 
       Implement the typical workflow
          opts = optimoptions('fmincon','Algorithm','sqp')
          problem = createOptimProblem('fmincon','objective', ...
          @(x) x(1)^2 + 4*sin(5*x(2)),'x0',[3 3],'lb',[-5 -5], ...
          'ub',[5 5],'nonlcon',@mycon,'options',opts)
          ms = MultiStart
          [x,f] = run(ms,problem,20)
 
    See also GLOBALSEARCH

    Documentation for MultiStart
       doc MultiStart

[Sadiq Akbar]

@Torsten

Thanks a lot for your guiadance. Is multistart a MATLAB builtin function? If yes, can you do it for me? I checked your above code, it works very well. But if the multistart can do it for all the required variables, then I think it will be more good. If I want to convert your above code to multistart, how will I do that?

Regards,

[Torsten]

If I want to convert your above code to multistart, how will I do that?

Do you have problems reading English documentations or are you just too lazy ?

[John D'Errico]

No. I won't write the code for you. You need to make some effort.

Multistart is a tool found in the global optimization toolbox. A great toolbox if you do much optimization. But not everybody has it. Actually, I don't have it myself, though I should pick it up one day.

However, it is also trivially easy to write a simple version. But, since I wont write code directly to do your job, here is an example you can read and follow:

% pick a simple function with many local minima
fun = @(x,y) sin(5*x + 2*y).*cos(-3*x + 7*y) + x.^2 + y.^2;
fsurf(fun,[-1,1,-1,1])

Of course, I know this has many local minima. Suppose you want to know how many minima lie in the box where -1<=x<=1 and -1<=y<=1. I'll use fmincon as an optimizer.

nsets = 100;
xyset = NaN(nsets,2);
fval = NaN(nsets,1);
obj = @(X) fun(X(1),X(2));
lb = [-1 -1];
ub = [1 1];
opts = optimset('fmincon');
opts.Display = 'none';
for i = 1:nsets
    % choose a start point sampled randomly in the box 
    % with -1<=x<=1 and -1<=y<=1
    xystart = rand(1,2)*2 - 1;
    % collect all solutinos found
    [xyset(i,:),fval(i)] = fmincon(obj,xystart,[],[],[],[],lb,ub,[],opts);
end

As you can see, there are many essentially duplicate solutions.

table(xyset,fval)

ans = 100×2 table
           xyset              fval  
    ____________________    ________

     -0.5479     0.63274    -0.27284
     0.10215     0.47726    -0.75051
      0.6967          -1     0.54447
     0.39953    -0.26275    -0.76113
      0.4301           1     0.47413
    -0.25093     -0.1073    -0.92037
      0.7521     0.32157    -0.28364
    -0.25093     -0.1073    -0.92037
      0.6967          -1     0.54447
     0.39953    -0.26275    -0.76113
     0.39953    -0.26275    -0.76113
    -0.25093     -0.1073    -0.92037
      0.7521     0.32157    -0.28364
     0.39953    -0.26275    -0.76113
     0.10215     0.47726    -0.75051
    -0.25093     -0.1073    -0.92037

First, do a simplistic clustering. Since I know the convergence tolerance on fmincon will be approximately 1e-8 by default, 1e-4 will be entirely adequate to cluster the solutions as a clustering tolerance.

[xyuniq,IA] = uniquetol(xyset,1e-4,byrows=true);
funiq = fval(IA);
table(xyuniq,funiq)

ans = 11×2 table
           xyuniq            funiq  
    ____________________    ________

    -0.90031    0.048006    -0.13521
    -0.60328    -0.69084    -0.10346
     -0.5479     0.63274    -0.27284
    -0.25093     -0.1073    -0.92037
    0.046774    -0.84697    -0.25164
     0.10215     0.47726    -0.75051
     0.39953    -0.26275    -0.76113
      0.4301           1     0.47413
      0.6967          -1     0.54447
      0.7521     0.32157    -0.28364
           1      0.8784      1.3155

% We found 11 essentially distinct solutions.
% plot the solutions found
plot(xyuniq(:,1),xyuniq(:,2),'o')

As you should expect, the solutions all fall along straight lines in the plane, although near the edges, that may fail. Even so, this did reasonably.

As far as writing the code to solve your specific poblem, you should be able to use my example. If not, it is time for you to to learn.

[Sadiq Akbar]

@Torsten

No no reading is not a problem for me. But as I have requested you that I am not so much expert in Matlab, that's why I request here for help. Actually I don't understand the documentation because Matlab has used very tough terminlogy and the mathematics behind is very difficult, so that's why I request you people for help.

Regards,

[Sadiq Akbar]

@John D'Errico

Thanks a lot. Ok I try myself to convert it for my function. Regards,