paretosearch does not satisfy nonlinear constraints

I want to make V<=0.001 at a position where Wn=1, but the non-linear constraint cannot do it. I am sure there is this answer. The following is my program。

https://drive.mathworks.com/sharing/df23a873-b7f8-420c-b085-dedccef7860c

clear;
fun = @object1;
nonlcon = @unitdisk52;
nvar=5;
lb = [0,0,0,0,0];
ub = [1,1,1,2*pi,2*pi];
A = [0 0 0 -1 0;0 0 0 1 -1];
b = [0;0];
Aeq = [1 1 1 0 0];
beq = [1];
opts = optimoptions(@paretosearch,'PlotFcn','psplotparetof','ParetoSetChangeTolerance',1e-3);
[x,fval,exitflag,output] = paretosearch(fun,nvar,A,b,Aeq,beq,lb,ub,nonlcon,opts)

Unable to find a feasible point.


x =

     []


fval =

     []

exitflag = -2

output = struct with fields:
         iterations: 22
          funccount: 0
             volume: []
    averagedistance: []
             spread: []
      maxconstraint: []
            message: 'Unable to find a feasible point.'
           rngstate: [1x1 struct]
        maxmeshsize: []

% max_iter = 100;
% for iter = 1:max_iter
% opts = optimoptions("gamultiobj","PlotFcn","gaplotpareto",...
%     "PopulationSize",100, 'ParetoFraction', 0.60,"ConstraintTolerance",1e-2);
% [x,fval,exitflag,output] = gamultiobj(fun,nvar,A,b,Aeq,beq,lb,ub,nonlcon,opts)               %mo   
% [c,ceq] =unitdisk52(x)
% d=pi;
% e=2*pi;
% Wn=x(1,4)+(x(1,5)-x(1,4))*0.5;
% V=sqrt((x(1,1)+(x(1,2)*cos(Wn*d))+(x(1,3)*...
%     cos(Wn*e))).^2+(x(1,2)*sin(Wn*d)+...
%     x(1,3)*sin(Wn*e)).^2)
% if (exitflag ==1) && V<0.05
%     break;
% end
% end
% opts = optimoptions("gamultiobj","PlotFcn","gaplotpareto",...
%     "PopulationSize",100, 'ParetoFraction', 0.60,"ConstraintTolerance",1e-2);
% [x,fval,exitflag,output] = gamultiobj(fun,nvar,A,b,Aeq,beq,lb,ub,nonlcon,opts)               %mo
%% curve test
figure;
% length(x)
gg=1:1:length(x); 
Wn=0:0.001:2;
A0=x(gg,1);

Index in position 2 exceeds array bounds.

A1=x(gg,2);
A2=x(gg,3);
t1=x(gg,4);
t2=x(gg,5);
% t1=x(gg,6);
% t2=x(gg,7);
C=A0+(A1.*cos(Wn.*t1))+(A2.*cos(Wn.*t2));
S=A1.*sin(Wn.*t1)+A2.*sin(Wn.*t2);
V=sqrt(C.^2+S.^2);
plot(Wn,V)
yline(0.05,'-.k');
array1=zeros(1,length(gg));
% for gg=1:60
% WL=find(V(gg,:)<=0.05,1,'first');
% WH=find(Wn(1,:)>=1&V(gg,:)>=0.05001,1,'first')-1;
% array1(1,gg)=Wn(1,WH)-Wn(1,WL);
% end
% trapz(V);
% plot(Wn,V)
% [row1,col1]=find(V<=0.05,1,'first');
% [row2,col2]=find(Wn>=1&V>0.0500000000001,1,'first');
% I=Wn(row2,col2-1)-Wn(row1,col1)
% yline(0.05,'-.k');
hold on

function f=object1(x)
%Gaussian+面積
Wn=0:0.001:2; 
V=sqrt((x(1)+(x(2)*cos(Wn*x(4)))+(x(3)*...
    cos(Wn*x(5)))).^2+(x(2)*sin(Wn*x(4))+...
    x(3)*sin(Wn*x(5))).^2);
% guassian = normpdf(Wn,1,0.2);

%頻率寬度倒數
% [col1]=find(V<=0.05,1,'first');
% [col2]=find(Wn>1&V>0.05000001,1,'first');
% I=Wn(1,col2-1)-Wn(1,col1);

f1=x(5);
% f2=1./I;
f2=trapz(V)/2000;

f=[f1 f2];
end

function [c,ceq] =unitdisk52(x)
Wn=1;
Cn=((x(1))+(x(2)*cos(Wn.*x(4)))+(x(3)*cos(Wn.*x(5)))).^2;
Sn=((x(2)*sin(Wn.*x(4)))+(x(3)*sin(Wn.*x(5)))).^2;
c=[(sqrt(Cn+Sn)-0.001)];
ceq=[];
x;
end

[Walter Roberson]

A = [0 0 0 -1 0;0 0 0 1 -1];
b = [0;0];

The first row of A is enforcing that -x(4) < = 0

lb = [0,0,0,0,0];
ub = [1,1,1,2*pi,2*pi];

but x is non-negative, so -x(4) is always <= 0

So the first row of A (together with the first row of b) is unnecessary

[ray]

Thank you for your answer. I have completed the changes, but the same problem still occurs. Are there any other errors? Thank you.

clear;
fun = @object1;
nonlcon = @unitdisk52;
nvar=5;
lb = [0,0,0,0,0];
ub = [1,1,1,2*pi,2*pi];
A = [0 0 0 1 -1];
b = [0];
Aeq = [1 1 1 0 0];
beq = [1];
opts = optimoptions(@paretosearch,'PlotFcn','psplotparetof','ParetoSetChangeTolerance',1e-3);
[x,fval,exitflag,output] = paretosearch(fun,nvar,A,b,Aeq,beq,lb,ub,nonlcon,opts)

[Walter Roberson]

Testing by deactivating Aeq and beq, and later select the locations that fit Aeq and beq:

fun = @object1;
nonlcon = @unitdisk52;
nonlcon = [];
nvar=5;
lb = [0,0,0,0,0];
ub = [1,1,1,2*pi,2*pi];
A = [0 0 0 1 -1];
b = [0];
Aeq_ = [1 1 1 0 0];
beq_ = [1];
Aeq = [];
beq = [];
opts = optimoptions(@paretosearch,'PlotFcn','psplotparetof','ParetoSetChangeTolerance',1e-3);
[x,fval,exitflag,output] = paretosearch(fun,nvar,A,b,Aeq,beq,lb,ub,nonlcon,opts)

Paretosearch stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

x = 55x5

1.0e+00 * 0 0 0 0.0003 0.0003 -0.0000 0.0000 -0.0000 0.0003 0.0003 -0.0000 0.0000 -0.0000 0.0003 0.0003 -0.0000 -0.0000 0.0000 0.0003 0.0003 -0.0000 0.0000 -0.0000 0.0003 0.0003 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0003 0.0003 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0003 0.0003

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 55x2

1.0e+00 * 0.0003 0 0.0003 0.0000 0.0003 0.0000 0.0003 0.0000 0.0003 0.0000 -0.0000 0.0000 0.0003 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0003 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

exitflag = 0

output = struct with fields:
         iterations: 150
          funccount: 21000
             volume: 3.0006
    averagedistance: 0.0997
             spread: 2.6841
      maxconstraint: 0
            message: 'Paretosearch stopped because it exceeded the function evaluation limit set by ...'
           rngstate: [1x1 struct]

format long g

mask2 = abs(x * Aeq_.' - beq_.') <= 0.0001;
x_match = x(mask2, :)

x_match = 1x5

3.8637781108584e-10 0 1 -8.82470907060096e-09 -8.90473036352407e-09

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval_match = fval(mask2)

fval_match = 
     -8.90473036352407e-09

This one match is pretty much at the boundary conditions, except that it has negative components for x_match while lb for x_match should be strictly zero.

function f=object1(x)
%Gaussian+面積
Wn=0:0.001:2; 
V=sqrt((x(1)+(x(2)*cos(Wn*x(4)))+(x(3)*...
    cos(Wn*x(5)))).^2+(x(2)*sin(Wn*x(4))+...
    x(3)*sin(Wn*x(5))).^2);
% guassian = normpdf(Wn,1,0.2);
%頻率寬度倒數
% [col1]=find(V<=0.05,1,'first');
% [col2]=find(Wn>1&V>0.05000001,1,'first');
% I=Wn(1,col2-1)-Wn(1,col1);
f1=x(5);
% f2=1./I;
f2=trapz(V)/2000;
f=[f1 f2];
end
function [c,ceq] =unitdisk52(x)
Wn=1;
Cn=((x(1))+(x(2)*cos(Wn.*x(4)))+(x(3)*cos(Wn.*x(5)))).^2;
Sn=((x(2)*sin(Wn.*x(4)))+(x(3)*sin(Wn.*x(5)))).^2;
c=[(sqrt(Cn+Sn)-0.001)];
ceq=[];
x;
end

Answer 1

[Torsten]

You may get feasible initial points if you run this code before calling "paretosearch" (maybe for different x0 to get several feasible points). The solution vectors can then be used in the call to "paretosearch" in the options structure as "InitialPoints".

lb = [0,0,0,0,0];
ub = [1,1,1,2*pi,2*pi];
A = [0 0 0 1 -1];
b = [0];
Aeq = [1 1 1 0 0];
beq = [1];
x0 = [1/3,1/3,1/3,1,2];
[sol1,fval] = fmincon(@obj,x0,A,b,Aeq,beq,lb,ub);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

x0 = [0.7 0.1 0.2 50 60];
[sol2,fval] = fmincon(@obj,x0,A,b,Aeq,beq,lb,ub);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

fun = @object1;
nonlcon = @unitdisk52;
nvar=5;
opts = optimoptions(@paretosearch,'PlotFcn','psplotparetof','ParetoSetChangeTolerance',1e-3,'InitialPoints',[sol1;sol2]);
[x,fval,exitflag,output] = paretosearch(fun,nvar,A,b,Aeq,beq,lb,ub,nonlcon,opts)

Pareto set found that satisfies the constraints. 

Mesh size of all incumbents less than 'options.MeshTolerance' and constraints are 
satisfied to within 'options.ConstraintTolerance'.

x = 6x5

0.2634 0.4920 0.2446 2.2571 4.5243 0.2813 0.4619 0.2569 2.4339 4.9032 0.2813 0.4619 0.2569 2.4340 4.9033 0.2813 0.4619 0.2569 2.4341 4.9033 0.2813 0.4619 0.2569 2.4341 4.9033 0.2813 0.4619 0.2569 2.4340 4.9033

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 6x2

4.5243 0.3847 4.9032 0.3776 4.9033 0.3776 4.9033 0.3776 4.9033 0.3776 4.9033 0.3776

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

exitflag = 1

output = struct with fields:
         iterations: 70
          funccount: 139
             volume: 1.3861
    averagedistance: 0.0300
             spread: 0.8656
      maxconstraint: 0
            message: 'Pareto set found that satisfies the constraints. ...'
           rngstate: [1x1 struct]

figure;
% length(x)
gg=1:1:length(x); 
Wn=0:0.001:2;
A0=x(gg,1);
A1=x(gg,2);
A2=x(gg,3);
t1=x(gg,4);
t2=x(gg,5);
% t1=x(gg,6);
% t2=x(gg,7);
C=A0+(A1.*cos(Wn.*t1))+(A2.*cos(Wn.*t2));
S=A1.*sin(Wn.*t1)+A2.*sin(Wn.*t2);
V=sqrt(C.^2+S.^2);
plot(Wn,V)
yline(0.05,'-.k');

function value = obj(x)
  Wn=1.5;
  Cn=((x(1))+(x(2)*cos(Wn.*x(4)))+(x(3)*cos(Wn.*x(5)))).^2;
  Sn=((x(2)*sin(Wn.*x(4)))+(x(3)*sin(Wn.*x(5)))).^2;
  value=sqrt(Cn+Sn);
end

function f=object1(x)
%Gaussian+面積
Wn=0:0.001:2; 
V=sqrt((x(1)+(x(2)*cos(Wn*x(4)))+(x(3)*...
    cos(Wn*x(5)))).^2+(x(2)*sin(Wn*x(4))+...
    x(3)*sin(Wn*x(5))).^2);
% guassian = normpdf(Wn,1,0.2);
%頻率寬度倒數
% [col1]=find(V<=0.05,1,'first');
% [col2]=find(Wn>1&V>0.05000001,1,'first');
% I=Wn(1,col2-1)-Wn(1,col1);
f1=x(5);
% f2=1./I;
f2=trapz(V)/2000;
f=[f1 f2];
end

function [c,ceq] =unitdisk52(x)
Wn=1.5;
Cn=((x(1))+(x(2)*cos(Wn.*x(4)))+(x(3)*cos(Wn.*x(5)))).^2;
Sn=((x(2)*sin(Wn.*x(4)))+(x(3)*sin(Wn.*x(5)))).^2;
c=[(sqrt(Cn+Sn)-0.001)];
ceq=[];
x;
end