paretosearch does not satisfy nonlinear constraints

8 次查看（过去 30 天）
ray
ray 2024-8-25
编辑： Torsten 2024-8-26
Ran in:
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
  3 个评论
显示 1更早的评论
ray
ray 2024-8-26
编辑：ray 2024-8-26
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 2024-8-26
Ran in:
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-8-26
编辑：Torsten 2024-8-26
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Write Constraints 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by