paretosearch does not satisfy nonlinear constraints
3 views (last 30 days)
Show older comments
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)
% 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);
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
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)
format long g
mask2 = abs(x * Aeq_.' - beq_.') <= 0.0001;
x_match = x(mask2, :)
fval_match = fval(mask2)
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
Accepted Answer
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);
x0 = [0.7 0.1 0.2 50 60];
[sol2,fval] = fmincon(@obj,x0,A,b,Aeq,beq,lb,ub);
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)
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 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!