Using Multistart with constrainted fmincon

Hey,

i want to optimize a constrainted function using fmincon. To find the global optimum, i found out that the use of MultiStart is helpful. My function looks like this

Here is my code for the use of fmincon

k = [5,7];
p = [1:4];


a = [pi/11, pi/7, pi/6, pi/3]; %Initial guesses
weights=(1./k.^4);
weights=weights/sum(weights);
fun = @(a) sqrt(sum( weights(:).*( 1 + 2*sum( (-1).^p.*cos(k(:).*a(p)) ,2) ).^2) );

A = [1, -1,0,0;0, 1, -1, 0; 0,0,1,-1]; % a1 < a2 < a3 ..
B = [0,0,0];
Aeq = [];
Beq = [];
lb = [0,0,0,0];
ub = [pi/4,pi/4,pi/4,pi/4];

function [c,ceq] = modindex(a,p) % saved as a seperate script
c = [];
ceq = 4/pi.*( 1 + 2*sum( (-1).^p.*cos(a(p)) ,2 ))-0.6;
end

nlcon = @(a) modindex(a,p); % this is a constraint, implemented as a function
x = fmincon(fun, a, A, B, Aeq, Beq, lb, ub,nlcon)

For the use of fmincon this code works, but now i want to run it with MultiStart. My question is how to respect all of my constraints?

My first try was this

problem = createOptimProblem('fmincon', 'objective', fun, 'x0', [pi/11, pi/7, pi/6, pi/3],'A', [1, -1,0,0;0, 1, -1, 0; 0,0,1,-1],...
    'b',[0,0,0], 'Aeq', [], 'beq', [],  'lb', [0,0,0,0], 'ub', [pi/4,pi/4,pi/4,pi/4], 'nonlcon', @(a) modindex(a,p));
ms = MultiStart( 'UseParallel', 'allways', 'StartPointstoRun', 'bounds');
[x,f] = run(ms, problem, 10)

This doesnt work, "No field A exists for PROBLEM structure."

What does my code must look like, that all of the constraints are respected?

Kind regards

Answer 1

[Matt J]

problem = createOptimProblem('fmincon', 'objective', fun, 'x0', [pi/11, pi/7, pi/6, pi/3],'Aineq', [1, -1,0,0;0, 1, -1, 0; 0,0,1,-1],...
    'bineq',[0,0,0], 'Aeq', [], 'beq', [],  'lb', [0,0,0,0], 'ub', [pi/4,pi/4,pi/4,pi/4], 'nonlcon', @(a) modindex(a,p));

[Meikel Vollmers]

Now it works, thanks.

But for this constraint

function [c,ceq] = modindex(a,p)
c = [];
ceq = 4/pi.*( 1 + 2*sum( (-1).^p.*cos(a(p)) ,2 ))-m;
end

with m = 0.41 i get these results

x =

    0.0000    0.7854    0.7429    0.7854
    

So the constraints are not respected because x(2) is higher than x(3) and equal to x(4).

[Matt J]

The real question is not whether MultiStart succeeded, but whether it thinks it succeeded. You should be running with an exitflag and other diagnostic output to check.

[Meikel Vollmers]

Okay, the exitflag output says -2, "No feasible local minimum found."

In such case, what do i have to do now? Do i have to increase the amount of starting points?

[Matt J]

How do we know the problem is feasible?

[Meikel Vollmers]

The goal is to minimize the Total harmonic distortion of current, i read some paper where this is done for 0<m<1.2

[Matt J]

Well, it does look like there are feasible points. Perhaps you can initialize with those found by the combinatoric search below.

p=1:4;

[A{1:4}]=ndgrid(linspace(0,pi/4,50));

A=cell2mat( cellfun(@(x) x(:), A,'uni',0));

ceq = abs( 4/pi.*( 1 + 2*sum( (-1).^p.*cos(A) ,2 ))-0.6 );

aFeasible = A( ceq/max(ceq)<1e-6 , :)

aFeasible = 8×4

0.2244 0.7854 0.2244 0.2084 0.3206 0.7373 0.0321 0.3366 0.0321 0.7373 0.3206 0.3366 0.3206 0.3366 0.0321 0.7373 0.6732 0.7373 0.2725 0.7373 0.0321 0.3366 0.3206 0.7373 0.2725 0.7373 0.6732 0.7373 0.2244 0.2084 0.2244 0.7854

[Meikel Vollmers]

When i use that combinatoric search for m = 0.4, aFeasible is a double empty matrix.

Maybe the implementation of my function is not correct, but in the paper it works out, i'm confused

[Matt J]

I don't think you've mentioned what m is. I don't see it anywhere in your code.

[Meikel Vollmers]

function [c,ceq] = modindex(a,p) % saved as a seperate script
c = [];
ceq = 4/pi.*( 1 + 2*sum( (-1).^p.*cos(a(p)) ,2 ))-m;
end

m ist mentioned in the function modindex. This is the right code, when i calculate the optimum i set the desired value in place of m manually. And the problem occurs when the value of m is lower than 0.5

[Matt J]

If you slacken the search tolerance a little bit, m=0.4 does give some hits

p=1:4;

[A{1:4}]=ndgrid(linspace(0,pi/4,50));

A=cell2mat( cellfun(@(x) x(:), A,'uni',0));

ceq = abs( 4/pi.*( 1 + 2*sum( (-1).^p.*cos(A) ,2 ))-0.4);

aFeasible = A( ceq/max(ceq)<2e-6 , :)

aFeasible = 8×4

0.3526 0.7373 0.0160 0.5450 0.0160 0.7373 0.3526 0.5450 0.5129 0.7373 0.1122 0.6732 0.1122 0.7373 0.5129 0.6732 0.3526 0.5450 0.0160 0.7373 0.5129 0.6732 0.1122 0.7373 0.0160 0.5450 0.3526 0.7373 0.1122 0.6732 0.5129 0.7373

[Matt J]

Also if you discretize the search region more finely,

[A{1:4}]=ndgrid(linspace(0,pi/4,100));

[Meikel Vollmers]

Okay, but what does that mean to the use of MultiStart? Do i only have to set the tolerances in Multistart to other values?

[Matt J]

It means you can do,

[x,fval,exitflag,output] = run(ms,problem, CustomStartPointSet(aFeasible))

[Meikel Vollmers]

The exitflag stays -2 and the output says:

output =

struct with fields:

funcCount: 2010

localSolverTotal: 8

localSolverSuccess: 0

localSolverIncomplete: 0

localSolverNoSolution: 8

message: 'No feasible solution found.↵↵MultiStart called the local solver 8 times and did not find a point that satisfies the constraints↵within the local solver constraint tolerance (problem.options.ConstraintTolerance).'

So it does not work out

[Matt J]

We forgot to include your linear inequality constraints in the combinatoric search. When we do, it does indeed appear that the problem is infeasible, even with very lax tolerances on the constraints:

tol=1e-3;

p=1:4;

[AA{1:4}]=ndgrid(linspace(0,pi/4,50));

AA=cell2mat( cellfun(@(x) x(:), AA,'uni',0));

ceq = abs( 4/pi.*( 1 + 2*sum( (-1).^p.*cos(AA) ,2 ))-0.4 );

a = AA( ceq/max(ceq)<tol , :).';

A = [1, -1,0,0;0, 1, -1, 0; 0,0,1,-1]; % a1 < a2 < a3 ..
B = [0,0,0].';

idx=all(A*a<=B+tol,1);

aFeasible=a(:,idx).'

aFeasible =

  0×4 empty double matrix

[Meikel Vollmers]

Could it be that there is a mistake in the code of fun?

[Matt J]

In the objective function? The objective function has no bearing on whether the problem is feasible or not. Only the constraints do. There could be a mistake in your constraints, but only you can know how the problem should be defined.

[Meikel Vollmers]

i got 2 constraints:

this includes the m, which i tried to realize using the modindex(a,p) function.

and the second is realized in the 'Aineq' and 'bineq' constraints

[Torsten]

And the bounds give

0 <= a1 <= a2 <= ... <= ap <= π/4

[Matt J]

We can prove mathematically that the smallest value m can have is about 0.527393087579050.

First, because the sequence of are monotonic in [0,],

Therefore,

So, you cannot consider any m lower than this value.

Also, this is a tight lower bound, achieved for example by choosing and all other .

[Meikel Vollmers]

Yeah, got it, youre right. Thats a mistake made in the paper i read. Matt, youre the best, thanks a lot!

I will replace the upper bound with pi/2 instead of pi/4. That will also guarantee the quater and halfwave symmetrie of my pulses.