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!
回答（1 个）
Alan Weiss
2018-3-14
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
类别
在 帮助中心 和 File Exchange 中查找有关 Surrogate Optimization 的更多信息
2018-3-13
2018-3-14
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
