通过PSO实现不同函数的目标值计算和搜索

Posted 2022-01-26 fpga&matlab

 

 

 

clc
clear
% 

time = cputime;

% computering run time
w = 1;        % inertia weight
c1 = 2;        % acceleration
c2 = 2;        % acceleration
eta =0.9;    %  updatingthe value of inertia weight in every iteration
%global  Npop ;    % particles number
Npop  = 25   ;
%Vmax = 25;    % the max value of velocity
%Vmin = 0;    % the min value of velocity
%threshold = 1;    % 
itera = 1;    % the times of iteration

rand('state',sum(100*clock));    % generate variable random numbers with time

% initialization
x1 = rand(Npop, 1)       + 2.6;             %  x1 = [2.6, 3.6]
x2 = rand(Npop, 1)*0.1   + 0.7;             %  x2 = [0.7, 0.8]
x3 = round(rand(Npop, 1) * 11 + 17);        %  x3 = [17, 28] integer value
x4 = rand(Npop, 1)       + 7.3;             %  x4 = [7.3, 8.3]
x5 = rand(Npop, 1)       + 7.3;             %  x5 = [7.3, 8.3]
x6 = rand(Npop, 1)       + 2.9;             %  x6 = [2.9, 3.9]
x7 = rand(Npop, 1)*0.5   + 5.0;             %  x7 = [5.0, 5.5]

v1max = 0.3*( max(x1) - min(x1));
v2max = 0.3*( max(x2) - min(x2));
v3max = 0.3*( max(x3) - min(x3));
v4max = 0.3*( max(x4) - min(x4));
v5max = 0.3*( max(x5) - min(x5));
v6max = 0.3*( max(x6) - min(x6));
v7max = 0.3*( max(x7) - min(x7));

v1 = rand(Npop, 1) * v1max;    % 
v2 = rand(Npop, 1) * v2max;    % 
v3 = rand(Npop, 1) * v3max;    % 
v4 = rand(Npop, 1) * v4max;    % 
v5 = rand(Npop, 1) * v5max;    % 
v6 = rand(Npop, 1) * v6max;    % 
v7 = rand(Npop, 1) * v7max;    % 
% computing the fitness value of every particles
fit = fitness(x1, x2, x3, x4, x5, x6, x7,itera);

% 
%pbest = [x1 x2 x3 x4 x5 fit];                % pbest
pworst = [x1 x2 x3 x4 x5 x6 x7 fit];

% 
%[maxfit, indfit] = max(fit);                % find the particle whose fitness value is the max
%gbest = [x1(indfit) x2(indfit) x3(indfit) x4(indfit) x5(indfit) maxfit];     % take that particle as the gbest particle

% 
[minfit, indfit] = min(fit);                % find the particle whose fitness value is the min
gworst = [x1(indfit) x2(indfit) x3(indfit) x4(indfit) x5(indfit) x6(indfit) x7(indfit) minfit];    % take that particle as the gworst

%the error between the gbest and gworst
%error = (maxfit - minfit)/maxfit;           % 

% -------------------------------------------------------------------------
% --◆initialization complete
G_BEST = zeros(1,401);
ITERA_N =1000;
itera = 0;
Diversity = zeros(ITERA_N+1,1);
dir       =    1;
DivH      =   0.001;
DivL      =   0.00005;

while itera <= ITERA_N,                        % start iteration

    itera = itera + 1;        %  the times of iteration
    
    Diversity(itera)  =  (var(x1) + var(x2) + var(x3) + var(x4) + var(x5) + var(x6) + var(x7))/Npop ;
 

    
    h1 = rand;    h2 = rand;    % random number: h1 and  h2
   % if Diversity(itera) > DivH  & dir < 0
   %     dir = 1; 
   % elseif Diversity(itera) < DivL & dir > 0
    %    dir = -1;
  %  end
        
        
        v1 = w * v1 + dir*(c1 * h1 * (pworst(:,1) - x1) + c2 * h2 * (gworst(1) - x1));    % updatingv1
        v2 = w * v2 + dir*(c1 * h1 * (pworst(:,2) - x2) + c2 * h2 * (gworst(2) - x2));    %   updating  v2
        v3 = w * v3 + dir*(c1 * h1 * (pworst(:,3) - x3) + c2 * h2 * (gworst(3) - x3));    %   updating  v3
        v4 = w * v4 + dir*(c1 * h1 * (pworst(:,4) - x4) + c2 * h2 * (gworst(4) - x4));    %   updating  v4
        v5 = w * v5 + dir*(c1 * h1 * (pworst(:,5) - x5) + c2 * h2 * (gworst(5) - x5));  %   updating  v5
        v6 = w * v6 + dir*(c1 * h1 * (pworst(:,6) - x6) + c2 * h2 * (gworst(6) - x6));  %   updating  v6
        v7 = w * v7 + dir*(c1 * h1 * (pworst(:,7) - x7) + c2 * h2 * (gworst(7) - x7));  %   updating  v7

    
    v1(find( v1 > v1max)) = v1max;  v1(find( v1 < -v1max)) = -v1max;
    v2(find( v2 > v2max)) = v2max;  v2(find( v2 < -v2max)) = -v2max;
    v3(find( v3 > v3max)) = v3max;  v3(find( v3 < -v3max)) = -v3max;
    v4(find( v4 > v4max)) = v4max;  v4(find( v4 < -v4max)) = -v4max;
    v5(find( v5 > v5max)) = v5max;  v5(find( v5 < -v5max)) = -v5max;
    v6(find( v6 > v6max)) = v6max;  v6(find( v6 < -v6max)) = -v6max;
    v7(find( v7 > v7max)) = v7max;  v7(find( v7 < -v7max)) = -v7max;
    
    
    
    x1 = x1 + v1;            % updating x1
    x2 = x2 + v2;            % updating x2
    x3 = round(x3 + v3);    % updating x3,integer value
    x4 = x4 + v4;            % updating x4
    x5 = x5 + v5;            % updating x5
    x6 = x6 + v6;            % updating x6
    x7 = x7 + v7;            % updating x7
    
    x1(find( x1 > 3.6 ))  = 3.6; x1(find( x1 < 2.6 )) = 2.6;   % 
    x2(find( x2 > 0.8 ))  = 0.8; x2(find( x2 < 0.7 )) = 0.7;   % make the values of x1, x2, x3, x4, x5 limited in its region
    x3(find( x3 > 28  ))  = 28 ; x3(find( x3 < 17 ))  = 17;
    x4(find( x4 > 8.3 ))  = 8.3; x4(find( x4 < 7.3 )) = 7.3;
    x5(find( x5 > 8.3 ))  = 8.3; x5(find( x5 < 7.3 )) = 7.3;
    x6(find( x6 > 3.9 ))  = 3.9; x6(find( x6 < 2.9 )) = 2.9;   
    x7(find( x7 > 5.5 ))  = 5.5; x7(find( x7 < 5.0 )) = 5.0;
    % updating pbest and gbest 
    [new_fit ]= fitness(x1, x2, x3, x4, x5, x6 ,x7,itera);                      % computing every particles' fitness value

    %for i = 1 : Npop,   % 1:40
    %    if  new_fit(i) > pbest(i,6),
    %        pbest(i,:) = [x1(i) x2(i) x3(i) x4(i) x5(i) new_fit(i)];    % updating pbest 
    %    end
    %end
    
    for i = 1 : Npop,   % 1:40
        if  new_fit(i) < pworst(i,8),
            pworst(i,:) = [x1(i) x2(i) x3(i) x4(i) x5(i) x6(i) x7(i) new_fit(i)];    % updating pbest 
        end
    end
    
    
% 
% [maxfit, indfit] = max(fit);                %
% gbest = [x1(indfit) x2(indfit) maxfit];

%    [maxfit(itera), indfit] = max(new_fit);             % find the max fitness value in this iteration 
%    if maxfit(itera) > gbest(3),
%        gbest = [x1(indfit) x2(indfit) maxfit(itera)];    % updating gbest     
%    end;

    
    [minfit(itera), indfit] = min(new_fit);             % find the min fitness value in this iteration 
    if minfit(itera) < gworst(8)
        gworst = [x1(indfit) x2(indfit) x3(indfit) x4(indfit) x5(indfit) x6(indfit) x7(indfit) minfit(itera)];    % updating gworst
    end;
 
    G_BEST(itera) = gworst(8);
    
    %  
    avgfit(itera) = mean(new_fit);    %computing the average fitness value of all particles 

    %  error
%    error = (maxfit(itera)- minfit(itera))/maxfit(itera);            % 
    
    %  w 
    w = w * eta;

% plotting result

%subplot(2,1,1)
%plot(x1, x2, '.',gbest(1),gbest(2),'r.')
%axis([-500 500 -500 500])
%title('The distribution of particles')
%xlabel('x position');
%ylabel('y position');

%subplot(2,1,2)

end;
%--------------------------------------◆
figure(1)
t = 1 : ITERA_N+1;        % 
plot(t, minfit,'b', t, avgfit,'g')
hold on
plot(t,G_BEST,'r','LineWidth',3);

axis([0 ITERA_N+10  2500 5000]);
title('fitness values');
xlabel('itera');
ylabel('fitness values');
legend('MinFit','AvgFit','Gbest')
x1      =      gworst(1)*ones(40,1);
x2      =      gworst(2)*ones(40,1);
x3      =      gworst(3)*ones(40,1);
x4      =      gworst(4)*ones(40,1);
x5      =      gworst(5)*ones(40,1);
x6      =      gworst(6)*ones(40,1);
x7      =      gworst(7)*ones(40,1);

g1      =      gcon1(x1, x2, x3);
g2      =      gcon2(x1, x2, x3);
g3      =      gcon3(x2, x3, x4, x6);
g4      =      gcon4(x2, x3, x5, x7);
g5      =      gcon5(x2, x3, x4, x6);
g6      =      gcon6(x2, x3, x5, x7);
g7      =      gcon7(x2, x4);
g8      =      gcon8(x1, x2);
g9      =      gcon9(x1, x2);
g10     =      gcon10(x4, x6);
g11     =      gcon11(x5, x7);

time = cputime - time
D131