my_pso_zyk_cgaus_rldc.m

给予粒子群算法的优点,将其用于电力系统中进行参数优化有明显效果

%PROGRAM: Modified PARTICLE SWARM OPTIMIZATION
% Particle Swarm Optimization - PSO
% Haupt & Haupt
% 2003
clear all
close all
%ff = 'testfunction'; % Objective Function
% Initializing variables
trnum=3;    %???????????????????????????????????????????????????????????????????????????????????????????????????????
Size=8;
G=12;                   %????????????????????????????????????????????????????????????????????????????????????????
Variables=3;              %????????????????????????????????????????????????????????????????????????????????????????
%************************************************  reconstraction;
data1=load('rldc.txt');
ttt=data1(:,1);
aaa=data1(:,2);

% data11=load('fenduan2000.txt');
% tttt=data11(1,:)';
% aaaa=data11(2,:)';

data=ttt;
%datat=tttt;
tau = 3;
m=6;
tsmax=max(data);
tsmin=min(data);
data=(data-tsmin)/(tsmax-tsmin);
%datat=(datat-tsmin)/(tsmax-tsmin);

n1=data(1:4:188);
x1=aaa(1:4:188)';
n2=data(1:3:188);
x2=aaa(1:3:188)';

% data11=data(850:1056);
% N=length(data11);
% [n1 x1]=reconstitution(data11,N,m,tau);%reconstitute state space
% 
% data11=data(1041:1104);
% N=length(data11);
% [n2 x2]=reconstitution(data11,N,m,tau);%reconstitute state space
%-------------------my program

Variablesmax=[10  1e-003   1e+009];%[20.133 0.001766 8.6698e+005][40.472 0.01 3.0316e+005][10 0.011705 5.8847e+005]%[2.0926 0.0089522 2729.7 0.010215 1][1.6323 0.0005266 2929.4 0.010799 1.3103][0.5 0.0001 3869.1 0.0051837 3.0148][1.75 2901.5 100 0.025];%[7  10000   100];
Variablesmin=[0.01 0   1e+003];%[1.75 2901.5 0.01 0.01];%[0.01  10 1 ];1.8939 0.001[2.7579 0.001097 3748.2 0.022802]

thresholdFun=0.006; 

failnum=0;  %记录失败次数
average_Fun=zeros(1,G);
average_k=0;  %平均代数（达到精度后所需代数）
average_time=0; 


fFun=zeros(trnum,G);

%npar = 2; % Dimension of the problem  问题维数，即变量个数
%maxit = 100; % Maximum number of iterations 最大代数
sum_c12=4;
%c1 = 1; % cognitive parameter 认知参数
c1 = 2;
%c2 = 4-c1; % social parameter  社会参数
c2 = sum_c12-c1; % social parameter  社会参数

CF=1; % constriction factor 收缩因子

ratio_par_vel=0.3;% 0.3;

velhi=ratio_par_vel*(Variablesmax-Variablesmin);
vello=-velhi;

popsize=Size;
npar=Variables;
maxit=G;

fBestFun=inf*ones(1,trnum);
fBestpar=ones(trnum,npar);

for tr=1:1:trnum,
    time1=cputime;
    flag=0;
% Initializing swarm and velocities
%par=rand(popsize,npar); % random population of continuous values 范围[0，1]
par=ones(Size,1)*Variablesmin+(ones(Size,1)*(Variablesmax-Variablesmin)).*rand(Size,Variables); % random   初始化个体
par(1,:)=[0.049628	4.718e-005	2.8568e+008];%[2.0926 0.0089522 2729.7 0.010215 1];%[2.2344 0.0071012 2206.3 0.021967 1];%[1.2067 546.94];
%vel = rand(popsize,npar); % random velocities

%vel=(velhi-vello)*rand(popsize,npar)+vello; % random   初始化个体
vel=zeros(popsize,npar);
% Evaluate initial population
%cost=feval(ff,par); % calculates population cost using ff
for s=1:1:popsize,
    %cost(s)=Test_function(par(s,:),funnum);
    GS=[0,s] 
    cost(s)=Reg_yun_cgaus_fenduan(par(s,:),n1,x1,n2,x2,tsmax,tsmin);
end
minc(1)=min(cost); % min cost
meanc(1)=mean(cost); % mean cost
globalmin(1)=minc(1); % initialize global minimum
% Initialize local minimum for each particle
localpar = par; % location of local minima
localcost = cost; % cost of local minima
% Finding best particle in initial population
[globalcost,indx] = min(cost);
globalpar=par(indx,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start iterations
w_max=0.9;
w_min=0.4;


iter = 0; % counter
while iter < maxit
iter = iter + 1;
w_flag=2;
if w_flag==1
%%%%%  w 自适应算法 1 %%%%%%%%
w=(maxit-iter)/maxit; %inertia weiindxht  
elseif w_flag==2
%%%%%  w 自适应算法 3 （《IIR数字滤波器设计的粒子群优化算法》） %%%%%%%%
w=w_max-iter*(w_max-w_min)/maxit;
end

%%%%%  c1 的自适应算法 %%%%%%%%
%c1=c1_max-iter*(c1_max-c1_min)/maxit;
%c2=c2_min+iter*(c2_max-c2_min)/maxit;
%c1=c1_min+iter*(c1_max-c1_min)/maxit;
%sum_c12=c12_min+iter*(c12_max-c12_min)/maxit;
%sum_c12=c12_max-iter*(c12_max-c12_min)/maxit;
%c2=sum_c12-c1;

r1 = rand(popsize,npar); % random numbers
r2 = rand(popsize,npar); % random numbers
vel = CF*(w*vel + c1 *r1.*(localpar-par) + c2*r2.*(ones(popsize,1)*globalpar-par)); %par为npar维的当前个体
for j=1:popsize
    vel(j,1:npar)=(vel(j,1:npar)<vello).*(vello)+(vel(j,1:npar)>=vello).*vel(j,1:npar);
    vel(j,1:npar)=(vel(j,1:npar)>velhi).*velhi+(vel(j,1:npar)<=velhi).*vel(j,1:npar);
end
% update particle positions
par = par + vel; % updates particle position
%overlimit=par<=1; %如果变量小于等于1，则overlimit＝1，
%underlimit=par>=0;%如果变量大于0，则overlimit＝1 
%par=par.*overlimit+not(overlimit); %如果par小于等于1，则par不变，否则赋值1
%par=par.*underlimit; %如果par小于0，则使par＝0
for j=1:popsize
    par(j,1:npar)=(par(j,1:npar)<Variablesmin).*(Variablesmin)+(par(j,1:npar)>=Variablesmin).*par(j,1:npar);
    par(j,1:npar)=(par(j,1:npar)>Variablesmax).*Variablesmax+(par(j,1:npar)<=Variablesmax).*par(j,1:npar);
end
% Evaluate the new swarm
%cost = feval(ff,par); % evaluates cost of swarm
for s=1:1:popsize
    GS=[iter,s] 
    cost(s)=Reg_yun_cgaus_fenduan(par(s,:),n1,x1,n2,x2,tsmax,tsmin);
end
% Updating the best local position for each particle
bettercost = cost < localcost;
localcost = localcost.*not(bettercost) +cost.*bettercost; %局部最佳个体
localpar(find(bettercost),:) =par(find(bettercost),:);
% Updating index g  更新全局最佳个体
[temp, t] = min(localcost); %从局部最佳个体找全局最佳个体，好主意,temp为所有局部最佳个体中的最小代价函数值，t为具有最小代价函数值的最佳局部个体所在的下标
if temp<globalcost
  globalpar=par(t,:); 
  indx=t; 
  globalcost=temp;
end
%[iter globalpar globalcost]; % print output each iteration
minc(iter+1)=min(cost); % min for this iteration 存放每代的最小代价，相当于GA中最大种群适应度值的倒数，也相当于GA中求最小问题的函数值
globalmin(iter)=globalcost; % best min so far
meanc(iter+1)=mean(cost); % avg. cost for this iteration
    if globalcost<=thresholdFun & flag==0   %检查是否满足   求最小值问题
        flag=1;
        average_k=average_k+ iter+1;
    elseif globalcost>thresholdFun & iter==maxit-1   %  检查是否失败   
        average_k=average_k+maxit;
        failnum=failnum+1;
    end     
end% while
fBestFun(tr)=globalcost;
fBestpar(tr,:)=globalpar;
average_Fun=average_Fun+globalmin;
fFun(tr,:)=globalmin;
average_time=average_time+cputime-time1;

tr
end % for
average_Fun=average_Fun/trnum;
average_k=average_k/trnum;
average_time=average_time/trnum;

average_fBestFun=mean(fBestFun);
Std_fBestFun=std(fBestFun);
index_Best=find(fBestFun==min(fBestFun));
index_Worst=find(fBestFun==max(fBestFun));
Best_fBestFun=fBestFun(index_Best(1));
Worst_fBestFun=fBestFun(index_Worst(1));

save outPSO_new fFun fBestpar globalpar failnum average_Fun average_k average_time fBestFun average_fBestFun Std_fBestFun Best_fBestFun Worst_fBestFun;
% outPSO_n_1 没有其他程序干扰的运行结果
figure(1)
iters=1:length(average_Fun);
%plot(iters,minc,'–',iters,globalmin,':');
plot(iters,average_Fun);
xlabel('generation');ylabel('cost');
%text(0,minc(1),'best');text(1,minc(2),'population average')