minimizing parameterized function with MultiStart (ideally in a loop)

Hi, I have an objective function in three variables (x(1),x(2),x(3)) and two parameters q1,q2. The variables (x(1),x(2),x(3)) must lie between 0 and 1 and must satisfy a non-linear constraint, x(1)*x(2) + (1 - x(2))*x(3) = p0, where p0 is another parameter in [0,1]. I would like to compute the maximum of the objective(minimum of -(objective)) as the parameter p0 varies, store the values in a vector v and plot it against p0.

As shown below, if I do the maximization for a single fixed value of p0, everything works out well. However, as soon as define p0 as a vector (0.01:0.01:1) and embed the maximization in a loop, everything breaks down and Matlab returns many errors. There's surely something wrong in my code. Hope someone can spot it!

The following is the code for the objective function, which I have in a script called Objective11.m:

   % code
function s = Objective11(x,q1,q2)
if (x(1)<=0.25) && (x(3)<=0.25)
s = -((4.*x(1)).*(x(2).* (q1 + q2 - 1) + (1 - q2)) + (4.*x(3)).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.25<x(1)) && (x(1)<=0.375) && (x(3)<=0.25)
s = -((-8.*x(1) + 3).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (4.*x(3)).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.25<x(1)) && (x(1)<=0.375) && (0.25<x(3)) && (x(3)<=0.375)
s = -((-8.*x(1) + 3).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.375<x(1)) && (x(1)<=0.5) && (x(3)<=0.25)
s = -((12.*x(1) - 9/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (4.*x(3)).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.375<x(1)) && (x(1)<=0.5) && (0.25<x(3)) && (x(3)<=0.375)
s = -((12.*x(1) - 9/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));   
elseif (0.375<x(1)) && (x(1)<=0.5) && (0.375<x(3)) && (x(3)<=0.5)
s = -((12.*x(1) - 9/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));  
elseif (0.5<x(1)) && (x(1)<=0.625) && (x(3)<=0.25) 
s = -((-12.*x(1) + 15/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (4.*x(3)).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.5<x(1)) && (x(1)<=0.625) && (0.25<x(3)) && (x(3)<=0.375)
s = -((-12.*x(1) + 15/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1))); 
elseif (0.5<x(1)) && (x(1)<=0.625) && (0.375<x(3)) && (x(3)<=0.5)
s= -((-12.*x(1) + 15/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.5<x(1)) && (x(1)<=0.625) && (0.5<x(3)) && (x(3)<=0.625)
s = -((-12.*x(1) + 15/2).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-12.*x(3) + 15/2).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.625<x(1)) && (x(1)<=0.75) && (x(3)<=0.25)
s = -((14.*x(1) - 35/4).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (4.*x(3)).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));    
elseif (0.625<x(1)) && (x(1)<=0.75) && (0.25<x(3)) && (x(3)<=0.375)
s = -((14.*x(1) - 35/4).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));    
elseif (0.625<x(1)) && (x(1)<=0.75) && (0.375<x(3)) && (x(3)<=0.5)
s = -((14.*x(1) - 35/4).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.625<x(1)) && (x(1)<=0.75) && (0.5<x(3)) && (x(3)<=0.625)
s = -((14.*x(1) - 35/4).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-12.*x(3) + 15/2).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.625<x(1)) && (x(1)<=0.75) && (0.625<x(3)) && (x(3)<=0.75)
s = -((14.*x(1) - 35/4).* (x(2) .* (q1 + q2 - 1) + (1 - q2)) + (14.*x(3) - 35/4).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.75<x(1)) && (x(1)<=1) && (x(3)<=0.25) 
s = -((-7.*x(1) + 7).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (4.*x(3)).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.75<x(1)) && (x(1)<=1) && (0.25<x(3)) && (x(3)<=0.375)
s = -((-7.*x(1) + 7).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.75<x(1)) && (x(1)<=1) && (0.375<x(3)) && (x(3)<=0.5)
s = -((-7.*x(1) + 7).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-8.*x(3) + 3).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.75<x(1)) && (x(1)<=1) && (0.5<x(3)) && (x(3)<=0.625)
s = -((-7.*x(1) + 7).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-12.*x(3) + 15/2).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));    
elseif (0.75<x(1)) && (x(1)<=1) && (0.625<x(3)) && (x(3)<=0.75)
s = -((-7.*x(1) + 7).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (14.*x(3) - 35/4).*((1 - x(2)).*(q2 + q1 - 1) + (1 - q1)));
elseif (0.75<x(1)) && (x(1)<=1) && (0.75<x(3)) && (x(3)<=1)    
s = -((-7.*x(1) + 7).*(x(2) .* (q1 + q2 - 1) + (1 - q2)) + (-7.*x(3)) + 7).*((1 - x(2))*(q2 + q1 - 1) + (1 - q1));   
end
end

The reason why it is so long is that it was constructed from a piecewise function with 6 branches. Anyway, for some input, it gives the right output so it should be correct.

Next, to minimize, I use MultiStart with fmincon. Using ony fmincon gives the local minima, which I don't need.

First, I've defined a function for the constraint, which I saved in a script called plausibility.m

    % code
function [c,ceq]=plausibility(x,p0)
c =[];
ceq=x(1)*x(2) + (1 - x(2))*x(3) - p0;
end

Then, I defined the problem and run the minimization in a third script, which is the following:

    % code
p0 = 3/10;
q1=1;
q2=1/2;
nonlcon=@(x)plausibility(x,p0);
objective = @(x)Objective11(x,q1,q2);
problem = createOptimProblem('fmincon','objective',objective,'x0',[p0,0,0],'lb',[p0,0,0],'ub',[1,1,p0],'nonlcon',nonlcon);
ms = MultiStart;
[x,fval] = run(ms,problem,30);

Thus, I assigned values to the three parameters that enter the objective function (q1,q2) and the constraint (p0). Then, I defined a function handle to the non-linear constraint in plausibility.m and defined the Optimization problem. There, I have a starting value x0, upper bound 'ub' and lower bound 'lb' that are correct. The code works and I can find the maximum as well as the maximizers.

    % code
MultiStart completed some of the runs from the start points.
29 out of 30 local solver runs converged with a positive local solver exit flag.
x =
    0.7500    0.1000    0.2500
fval =
   -1.4125
v =
    1.4125

which is good. I tried with different values for p0 and the output seems correct.

Now, i need the program to automatically compute the maximum and the maximizers of Objective11 for a series of p0 = 0.01,0.02,0.03,...etc and store them in a matrix. I wrote the following code to at least store in a vector v the maximum of the function at each p0.

    % code
p0 = (0.01:0.01:1);
q1 = 1;
q2=1/2;
n = length(p0);
v = zeros(n,1);
x0 = zeros(n,2);
lb = zeros(n,2);
fval = zeros(n);
ms = MultiStart;
objective = @(x)Objective11(x,q1,q2);
for i = 1:100
    x0(i,1:2) = [p0(i), 0];
    lb(i,1:2) = [p0(i), 0];
    ub = [1,1];
    nonlcon=@(x)plausibility(x,p0(i));
    problem(i) = createOptimProblem('fmincon','objective',objective,'x0',[p0(i),0],'lb',[p0(i),0],'ub',[1,1]);
    [fval(i)] = run(ms,problem(i),30);
    v(i) = -fval(i);
end

The non-linear equality constraint, as well as the starting value x0, the lower and the upper bound for the variables in x, vary with p0, so I inserted them in the loop. I'm not very sure about

    % code
nonlcon=@(x)plausibility(x,p0(i));

though.

When I run up to the line problem = .... I obtain no errors, and the program creates 100 different problems with the right x0,lb and ub (it doesn't seem to get the 'nonlcon' though). However, when I run also the last two lines, therefore the actual minimization, I get loads of errors:

    % code
Index exceeds matrix dimensions.
Error in Objective11 (line 17)
elseif (0.375<x(1)) && (x(1)<=0.5) && (x(3)<=0.25) %%%Third Constraint
Error in @(x)Objective11(x,q1,q2)
Error in fmincon (line 536)
      initVals.f = feval(funfcn{3},X,varargin{:});
Error in fmultistart
Error in MultiStart/run (line 260)
               fmultistart(problem,startPointSets,msoptions);
Caused by:
    Failure in initial objective function evaluation. FMINCON cannot continue.
    Failure in call to the local solver with supplied problem structure.

I can figure out the first. I don't understand the others. Especially 'Error in @(x)Objective11(x,q1,q2)', since I do the same thing when I maximize without the loop.

Please help me out!

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Answer 1

[Alan Weiss]

I suggest that you remove your first line of code:

x = [x(1),x(2)];

And, in addition to the link that Walter gave, you might want to consult this link, which is a slightly different description of how to pass extra parameters.

Your objective function looks nonsmooth. For all I know, it may be discontinuous. As such, I suggest that you use a solver that can handle nonsmoothness, such as patternsearch or, if you don't have a Global Optimization Toolbox license, try fmincon from a variety of start points, and maybe try the 'sqp' algorithm.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation