Clear Filters
Clear Filters

paretosearch does not satisfy nonlinear constraints

3 views (last 30 days)
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。
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
  3 Comments
ray
ray on 26 Aug 2024
Edited: ray on 26 Aug 2024
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
Walter Roberson on 26 Aug 2024
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

Sign in to comment.

Accepted Answer

Torsten
Torsten on 26 Aug 2024
Edited: Torsten on 26 Aug 2024
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

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!