1, Introduction

1 concept of particle swarm optimization algorithm
Particle swarm optimization (PSO) is an evolutionary computing technology. From the study of predation behavior of birds. The basic idea of particle swarm optimization algorithm is to find the optimal solution through the cooperation and information sharing among individuals in the group
The advantage of PSO is that it is simple and easy to implement, and there is no adjustment of many parameters. At present, it has been widely used in function optimization, neural network training, fuzzy system control and other application fields of genetic algorithm.

2 particle swarm optimization analysis
2.1 basic ideas
Particle swarm optimization (PSO) simulates birds in a flock of birds by designing a massless particle. The particle has only two attributes: speed and position. Speed represents the speed of movement and position represents the direction of movement. Each particle searches for the optimal solution separately in the search space, records it as the current individual extreme value, shares the individual extreme value with other particles in the whole particle swarm, and finds the optimal individual extreme value as the current global optimal solution of the whole particle swarm, All particles in the particle swarm adjust their speed and position according to the current individual extreme value found by themselves and the current global optimal solution shared by the whole particle swarm. The following dynamic diagram vividly shows the process of PSO algorithm:

2 update rules
PSO is initialized as a group of random particles (random solutions). Then the optimal solution is found through iteration. In each iteration, the particle updates itself by tracking two "extreme values" (pbest, gbest). After finding these two optimal values, the particle updates its speed and position through the following formula.

The first part of formula (1) is called [memory item], which represents the influence of the last speed and direction; The second part of formula (1) is called [self cognitive item], which is a vector from the current point to the best point of the particle itself, indicating that the action of the particle comes from its own experience; The third part of formula (1) is called [group cognition item], which is a vector from the current point to the best point of the population, reflecting the cooperation and knowledge sharing among particles. Particles determine the next movement through their own experience and the best experience of their companions. Based on the above two formulas, the standard form of PSO is formed.

Formulas (2) and (3) are considered as standard PSO algorithms.
3 PSO algorithm flow and pseudo code

2, Source code

clear
clc
close all
%% Parameter initialization
c1=2.05;
c2=2.05;
maxgen=5000;
sizepop=30;
k=0.6;
% wV=1.1;
wP=1.1;
v=5;
popmax=30;
popmin=-30;

%      pso_option = struct('c1',1.5,'c2',1.7,'maxgen',200,'sizepop',20, ...
%                    'k',0.6,'wV',1,'wP',1,'v',5, ...
%                     'popcmax',100,'popcmin',0.1,'popgmax',10^3,'popgmin',10^(-2));
% c1:The initial value is 1.5,pso Parameter local search capability
% c2:The initial value is 1.7,pso Parameter global search capability
% maxgen:Initially 200,Maximum evolutionary quantity
% sizepop:Initial 20,Maximum population
% k:Initial 0.6(k belongs to [0.1,1.0]),Rate and x Relationship between(V = kX)
% wV:The initial value is 1(wV best belongs to [0.8,1.2]),Elastic coefficient in front of velocity in rate update formula
% wP:The initial value is 1,Elasticity coefficient in front of velocity in population renewal formula
% v:Initially 5,SVM Cross Validation parameter
% popcmax:Initially 100,SVM parameter c Maximum value of change.
% popcmin:Initial 0.1,SVM parameter c Minimum value of change.
% popgmax:Initially 1000,SVM parameter g Maximum value of change.
% popgmin:Initial 0.01,SVM parameter c Minimum value of change.
D=10;  %%%dimension
Vmax =k*popmax;
Vmin = -Vmax ;

eps =1E-5;

%% Generate initial particles and velocities
pop=zeros(sizepop,D);
V=zeros(sizepop,D);
fitness=zeros(sizepop,1);
for i=1:sizepop
    
    % Randomly generated population and speed
    pop(i,:) = (popmax-popmin)*rand(1,D)+popmin;
    
    V(i,:)=Vmax*rands(1,D);
    
    % Calculate initial fitness
    fitness(i)=myfunc_fit1(pop(i,:));
    
end
Xd_ave0=repmat(sum(pop)/sizepop,sizepop,1);
D_t0=sum((sum((pop-Xd_ave0).^2,2)).^0.5)/sizepop/(popmax-popmin);
wV=1/(1+exp(-12*(D_t0-0.5)));

D_min=D_t0*0.2;

% Find extreme value and extreme point
[global_fitness bestindex]=min(fitness); % Global extremum
local_fitness=fitness;   % Individual extremum initialization

global_x=pop(bestindex,:);   % Global extreme point
local_x=pop;    % Individual extreme point initialization

% Average fitness of each generation population
avgfitness_gen = zeros(maxgen,1);
fit_gen=zeros(maxgen,1);

%% Iterative optimization
for i=1:maxgen
    
    for j=1:sizepop
        
        %Speed update
        V(j,:) = wV*V(j,:) + c1*rand*(local_x(j,:) - pop(j,:)) + c2*rand*(global_x - pop(j,:));
        if find(V(j,:) > Vmax)
            V_maxflag=find(V(j,:) > Vmax);
            V(j,V_maxflag) = Vmax;
        end
        if find(V(j,1) < Vmin)
            V_minflag=find(V(j,1) < Vmin);
            V(j,V_minflag) = Vmin;
        end
        
        %Population regeneration
        pop(j,:)=pop(j,:) + wP*V(j,:);
        if find(pop(j,:) > popmax)
            pop_maxflag=find(pop(j,:) > popmax);
            pop(j,pop_maxflag) = popmax;
        end
        if find(pop(j,:) < popmin)
            pop_minflag=find(pop(j,:) < popmin);
            pop(j,pop_minflag) = popmin;
        end
        
        % Adaptive particle mutation
        if rand>0.5
            k=ceil(2*rand);
            pop(j,k) = (popmax-popmin)*rand + popmin;
        end
        
        %Fitness value
        fitness(j)=myfunc_fit1(pop(j,:));
        
        %Individual optimal update
        if fitness(j) < local_fitness(j)
            local_x(j,:) = pop(j,:);
            local_fitness(j) = fitness(j);
        end
        
%         if fitness(j) == local_fitness(j) && length(pop(j,:) < local_x(j,:))
%             local_flag=find(pop(j,:) < local_x(j,:));
%             local_x(j,local_flag) = pop(j,local_flag);
%             local_fitness(j) = fitness(j);
%         end
        
        %Group optimal update
        if fitness(j) < global_fitness
            global_x = pop(j,:);
            global_fitness = fitness(j);
        end
        function C=PSO_FUNC(X)
global G_AC T_C hour_num Wind_V
C_W=110;   %%%Wind power generation
P_W=100;
u_PW=6;
m_WG=20;
r0=0.06;
v_ci=3; %Cut in wind speed
v_r=12; %Rated wind speed
v_co=25; %Cut off wind speed
P_r=P_W;
V_t=repmat(Wind_V(:,3),52,1);

C_S=0.7;   %%%%Photovoltaic power generation
P_S=0.2;
u_PS=0.009;
m_PV=25;

C_B=0.5;   %%%Battery
u_WB=0.0014;
m_B=5;

sigam_bat=1e-4;   %%%Self discharge rate
N_B=2000;
W_B=0.64;
Wbat_0=0.5*N_B*W_B;
Pbat_max=0.2*N_B*W_B;
Pbat_min=-0.2*N_B*W_B;
Pbat_maxt=Pbat_max;
Pbat_mint=Pbat_min;
Wbat_t=zeros(hour_num,1);
Wbat_t(1)=Wbat_0; %%%Remaining electricity

Pbat_t=zeros(hour_num,1);
Pbat_t(1)=Pbat_max;


C_d=10;     %%%Diesel generator
u_Pd=0.95;
P=4.62;
Q_d0=0.22;
m_die=10;

P_STC=0.2;
G_STC=1;
K=-0.47;
Tr=298.15;
T_C=T_C+273.15;

P1_t=300;     %%%Peak residential load  KW
Pdes_t=200;   %%%Desalination load   KW
P_des=25;      %%%Rated power of single seawater desalination unit   KW
N_des=8;     %%%Total number of seawater desalination units
G_des=100/24;    %%%The hourly water yield of a single unit is 100 t/d
Rwater_t=500/24;  %%%All day water demand on the island   t
Rdes_min=0;
Rdes_max=8*100/24;
eta_c=0.97;

Rdes_t=zeros(hour_num,1);
Rdes_t(1)=Rdes_max;   %%Initial storage 

P_PV=zeros(hour_num,1);
P_WG=zeros(hour_num,1);
P_PVM_t=zeros(hour_num,1);
P_WGM_t=zeros(hour_num,1);
P_net_t=zeros(hour_num,1);
Pdie_t=zeros(hour_num,1);
C_f=0;   %%%Annual cost of diesel
yeushu1=0;
yeushu2=0;
yeushu3=0;
yeushu4=0;
yeushu5=0;
yeushu6=0;
yeushu7=0;
yeushu8=0;
for i=1:hour_num
    P_PV(i)=P_STC*G_AC(i).*(1+K*(T_C(i)-Tr))/G_STC;
    a=P_r/(v_r^3-v_ci^3);
    b=v_ci^3/(v_r^3-v_ci^3);
    if (V_t(i)<v_ci)
        P_WG(i)=0;
    elseif (v_ci<V_t(i)<v_r)
        P_WG(i)=a*V_t(i)^3-b*P_r;
    elseif (v_r<V_t(i)<v_co)
        P_WG(i)=P_r;
    else
        P_WG(i)=0;
    end
    
    if(Rdes_t(i)-Rdes_max>=Rwater_t)
        Ndes_mint=0;
    else
       Ndes_mint=(Rwater_t-(Rdes_t(i)-Rdes_max))/G_des;
    end
    if(Rdes_t(i)+N_des*P_des-Rwater_t<=Rdes_max)
        Ndes_maxt=N_des*P_des;
    else
        Ndes_maxt=(Rdes_max+Rwater_t-Rdes_t(i))/P_des;
    end
    Pdes_mint=Ndes_mint*P_des;
    Pdes_maxt=Ndes_maxt*P_des;
    
    P_PVM_t(i)=P_PV(i);
    P_WGM_t(i)=P_WG(i);
    P_net_t(i)=P_WGM_t(i)+P_PVM_t(i)-P1_t;
    
   if(i>=2)
           Pbat_maxt=min([Pbat_max ((1-sigam_bat)*Wbat_t(i)-((1-sigam_bat)*Wbat_t(i)+Pbat_t(i)))*eta_c]);
           Pbat_mint=max([Pbat_min ((1-sigam_bat)*Wbat_t(i)-((1-sigam_bat)*Wbat_t(i)+Pbat_t(i)))/eta_c]);
   end   

3, Operation results

4, Remarks

Complete code or write on behalf of QQ 1564658423