% function [xm,fv] = PSO(N,c1,c2,w,M,D)
format long
clear;
clc;
N=500;c1=1.2;c2=1.2;wmax=1;wmin=0.01;M=200;D=16;
UB=64;LB=1;
lc=[0.111 0.141 0.196 0.25 0.307 0.391 0.442 0.563 0.602 0.766 0.785 0.994 1 1.228 1.266 1.457 1.563 1.62 1.8 1.99 2.13 2.38 2.62 2.63 2.88 2.93 3.09 3.13 3.38 3.47 3.55 3.63 3.84 3.87 3.88 4.18 4.22 4.49 4.59 4.8 4.97 5.12 5.74 7.22 7.97 8.53 9.3 10.85 11.5 13.5 13.9 14.2 15.5 16 16.9 18.8 19.9 22 22.9 24.5 26.5 28 30 33.5];
x=zeros(N,D);v=zeros(N,D);y=zeros(N,D);p=zeros(1,D);a=zeros(1,N);best=zeros(1,M);count=0;nouse=zeros(D,0);nouse1=zeros(D,0);
x1=zeros(N,D);y1=zeros(N,D);
yy=c1+c2;xx=2/abs(2-yy-sqrt(yy^2-4*yy));
%---------------------------------------------------------------------
% 参数初始的x、v，得到符合约束的N个初始x,写i=1:100000原因是100000个里总有N个符合约束的初始点；
%---------------------------------------------------------------------
for i=1:50000
for j=1:D
nouse(i,j)=randperm(64,1);
nouse1(i,j)=lc(nouse(i,j));
end
%---------------------------------------------------------------------
%强度约束
%---------------------------------------------------------------------
c = opt_intensity(nouse1(i,:));
if c(1)>0||c(2)>0||c(3)>0||c(4)>0
continue
end
count=count+1;
x(count,:)=nouse(i,:);
x1(count,:)=nouse1(i,:);
if count>=N
break
end
end
for i=1:N %参数初始的x、v
for j=1:D
v(i,j)=randperm(20,1)-5;
end
y(i,:)=x(i,:);
y1(i,:)=x1(i,:);
p(i)=opt_mass(x1(i,:));
end
[~,i]=min(p);
pg=x(i,:);
pg1=x1(i,:);
%---------------------------------------------------------------------
%开始逐代筛选
%---------------------------------------------------------------------
for t=1:M
for i=1:N
for j=1:D
w(j)=((wmax-wmin)*(M-t)/M+wmin);
v(i,j)=xx*(w(j)*v(i,j)+c1*rand*(y(i,j)-x(i,j))+c2*rand*(pg(j)-x(i,j)));
v(i,j)=round(v(i,j));
x(i,j)=x(i,j)+v(i,j);
end
%---------------------------------------------------------------------
%上下限约束
%---------------------------------------------------------------------
for j=1:D
if x(i,j)>UB
x(i,j)= UB;
end
if x(i,j)<LB
x(i,j)=LB;
end
x1(i,j)=lc(x(i,j));
end
a(i)=opt_mass(x1(i,:));
%---------------------------------------------------------------------
%强度约束
%---------------------------------------------------------------------
c = opt_intensity(x1(i,:));
if c(1)>0||c(2)>0||c(3)>0||c(4)>0
a(i)=a(i)+10^7;
end
%---------------------------------------------------------------------
%筛选
%---------------------------------------------------------------------
if a(i)<p(i)
p(i)=a(i);
y1(i,:)=x1(i,:);
y(i,:)=x(i,:);
end
if opt_mass(pg1)>p(i)
pg=y(i,:);
pg1=y1(i,:);
end
end
best(t)=opt_mass(pg1);
end
xm=pg1;
fv=min(best);
plot(1:M,best)
%
%
% end
