function [sys,x0,str,ts] = spso(t,x,u,flag)
persistent first
% if isempty(first)
% first = 0;
% end
persistent stop
% if isempty(stop)
% stop = 0;
% end
persistent i
% if isempty(i)
% i = 0;
% end
persistent mg
% if isempty(mg)
% mg = 1;
% end
%Initialize the parameters
popsize = 10;
maxgen = 30;
Wmin = 0.4; Wmax = 0.9;
Pmax = 4; Pmin = 0;
Vmax = 4; Vmin = -4; %注意min 必须是负数，速度应该可加可减
c1 = 2; c2 = 2;
persistent w ;
persistent Gbest;
persistent Gmaxval;
persistent P;
persistent V;
persistent Pbest;
persistent Pbestval;
% count = 0;
% w = getweight(Wmin, Wmax);
switch flag,
case 0,
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 1;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0=[];
str=[];
ts=[0.01 ,0];
disp(ts);
% initialize the static variables
first = 0;
stop = 0
i = 0;
mg = 1;
w = getweight(Wmin, Wmax)
Gbest = 0;
Gmaxval = 0;
P = zeros(1, popsize); %particles' position，matrix 1X10
V = zeros(1, popsize); %particles' velocity，matrix 1X10
% fitnessval = zeros(1, popsize);
Pbest = zeros(1, popsize);
Pbestval= zeros(1, popsize);
%initialize the particles' position and velocity, and limit the values
P(:) = rand(1, popsize) / (1 / (Pmax - Pmin)) %position
V(:) = rand(1, popsize) / (1 / (Vmax - Vmin)) %velocity
case 3,
%判断是否迭代完成
if mg > maxgen
stop = 1;
end
if stop == 1;
sys = Gbest;
return;
end
if first ==0
% disp(P);
first = 1;
sys = P(1);
return;
end
i = i + 1;
% count = count+1;
%First, calculate the fitness for each single particle
fitnessval =adaptfunc(P(i), u);
%Find the best position for every single particle
if i <= popsize
if(fitnessval > Pbestval(i))
Pbestval(i) = fitnessval;
Pbest(i) = P(i);
end
end
if i == popsize %种群迭代一次结束，则更新全局最优值，更新速度和位置
%Find the global best position in all particles
[best, best_index]= max( Pbestval(:));
if (best > Gmaxval)
Gmaxval = best;
Gbest = Pbest(best_index);
disp(mg);
disp(Gbest);
disp(Gmaxval);
end
%Update the velocity and limit it
for i = 1: popsize
newv(i) = w * V(i) + c1 * unifrnd(0, 1) * (Pbest(i) - P(i)) + c2 * unifrnd(0, 1) * (Gbest - P(i));
if newv(i) > Vmax
newv(i) = Vmax;
elseif newv(i) < Vmin
newv(i) = Vmin;
end
end
V = newv;
%end if
%Update the position and limit it
for i = 1: popsize
newp(i) = P(i) + V(i);
if newp(i) > Pmax
newp(i) = Pmax;
elseif newp(i) < Pmin
newp(i) = Pmin;
end
end
P = newp
% i = 0; %为下一次迭代做准备
% disp(Gmaxval);
first = 0;
i = 0;
mg = mg +1;
% stop = 1;
end
if i == popsize
sys = P(1);
else
sys = P(i+1);
end
% disp(i);
case {1,2,4,9},
sys = [];
end
function w=getweight(x, y)
w=x+rand(1)*(y - x);
function ret=adaptfunc(d, u)
% ret = 1- cos(3*x)*exp(-x);
ret = u;
function val=AdaptFuncForSinglePartical(size, p)
for i = 1: size
val(i) = adaptfunc(p(i));
end
