MyPsoClerc.rar_pso identification_pso辨识_粒子群 辨识_粒子群辨识_系统辨识

function [OutObj,OutCon] = MyPsoClerc(AdaptFunc,SwarmSize,ParticleSize,ParticlePosScope,ParticleVelScope,LoopCount) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %功能描述：向量化的粒子群算法，返回最后全局粒子的历史最优解及对应的适应度值。 % Clerc(1999)年提出的带收敛因子（constrictor factor）的PSO算法（该程序中的全局历史最优值采用在个体历史极值中选择的方法，不是在个体当前值中选择，实质一样） %输入参数：AdaptFunc：适应度函数（以min为最优），采用向量化方法 % [F] = Func(ParticleSwarm) 其中，F = [F1;F2;...;FN], x1 = ParticleSwarm(:,1) % 则这样使用 AdaptFunc = @MyRastrigin % 例如：[OutObj,OutCon] = MyPsoClerc(@MyRastrigin,100,3,[-5 5],[-4 4],150); %输入参数：SwarmSize：种群大小个数 %输入参数：ParticleSize：一个粒子的维数（说明中用D表示） %输入参数：ParticlePosScope：粒子各维的位置范围 % default [-100 100], of form: % [ min1 max1 % min2 max2 % ... % minD maxD ] %输入参数：ParticleVelScope：粒子各维的速度范围 % default [-4 4], of form: % [ min1 max1 % min2 max2 % ... % minD maxD ] %输入参数：LoopCount：最大迭代次数 default = 100 %输出参数：OutObj：经过迭代后得到的最优解及其适应度值 %输出参数：OutCon：输出约束函数的值以判断处理约束能力 %编制人：赵志斌（重庆大学自动化学院）Email: zhaozb1983@gmail.com %编制时间：2012.3.27 %参考文献：Particle swarm optimization: developments, applications and resources %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %输入变量检测及默认值设置 if nargin<3 error('至少输入3个变量：AdaptFunc,SwarmSize,ParticleSize'); elseif nargin>6 error('最多输入6个变量'); end if nargin==3 ParticlePosScope = [-100 100]; ParticleVelScope = [-4 4]; LoopCount = 100; elseif nargin==4 ParticleVelScope = [-4 4]; LoopCount = 100; elseif nargin==5 LoopCount = 100; end %输入变量维数检测及修正 if numel(SwarmSize)~=1 error('种群大小SwarmSize应该为1X1维'); end if numel(ParticleSize)~=1 error('粒子维数ParticleSize应该为1X1维'); end [row,colum] = size(ParticlePosScope); if ~((row==1||row==ParticleSize)&&(colum==2)) error('粒子各维位置范围ParticlePosScope应该为1X2维或ParticleSizeX2维'); elseif row==1 ParticlePosScope = repmat(ParticlePosScope,ParticleSize,1); %粒子位置范围扩展为ParticleSizeX2维 end [row,colum] = size(ParticleVelScope); if ~((row==1||row==ParticleSize)&&(colum==2)) error('粒子各维速度范围ParticleVelScope应该为1X2维或ParticleSizeX2维'); elseif row==1 ParticleVelScope = repmat(ParticleVelScope,ParticleSize,1); %粒子速度范围扩展为ParticleSizeX2维 end if numel(LoopCount)~=1 error('迭代次数LoopCount应该为1X1维'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %初始化粒子群，限定粒子的位置及速度在指定的范围内 %用矩阵ParSwarm表示当前粒子的位置、速度及当前适应度值。我们用W来表示位置，用V来代表速度， %用F来代表当前的适应度值,在这里我们假设粒子个数为N，每个粒子的维数为D。 % [ W11,W12,...,W1D,V11,V12,...,V1D,F1 % W21,W22,...,W2D,V21,V22,...,V2D,F2 % ...,...,...,...,...,...,...,...,... % WN1,WN2,...,WND,VN1,VN2,...,VND,FN ] %用矩阵OptSwarm记录每个粒子的历史最优解、全部粒子搜索到的全局历史最优解及其对应的适应度值。 %用Wg代表全局最优解，W.,1代表每个粒子的历史最优解。 % [ Wj1,Wj2,...,WjD,Fj % Wk1,Wk2,...,WkD,Fk % ...,...,...,...,... % Wm1,Wm2,...,WmD,Fm % Wg1,Wg2,...,WgD,Fg ] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %初始化粒子群矩阵及历史最优解矩阵，全部设为[0-1]随机数 %rng('shuffle'); %确保随机值每次都不同 %rng(0); %确保随机值每次都相同 ParSwarm = rand(SwarmSize,2*ParticleSize+1); %在粒子位置空间随机产生初始位置，在速度空间产生随机初始速度（对位置、速度范围进行调整） for k = 1:ParticleSize ParSwarm(:,k) = ParticlePosScope(k,1)+(ParticlePosScope(k,2)-ParticlePosScope(k,1))*ParSwarm(:,k); ParSwarm(:,ParticleSize+k) = ParticleVelScope(k,1)+(ParticleVelScope(k,2)-ParticleVelScope(k,1))*ParSwarm(:,ParticleSize+k); end %计算粒子适应度函数值，采用向量化的计算方法 [ParSwarm(:,2*ParticleSize+1),~,~] = AdaptFunc(ParSwarm(:,1:ParticleSize)); %初始化粒子群历史最优解矩阵,最优解矩阵全部设为零 OptSwarm = zeros(SwarmSize+1,ParticleSize+1); %历史位置最优解暂设置为初始值 OptSwarm(1:SwarmSize,1:ParticleSize) = ParSwarm(:,1:ParticleSize); %计算粒子适应度函数值，采用向量化的计算方法 [OptSwarm(1:SwarmSize,ParticleSize+1),~,~] = AdaptFunc(OptSwarm(1:SwarmSize,1:ParticleSize)); %寻找目标函数值最小的解在个体极值矩阵中的位置(行数) [~,row] = min(OptSwarm(1:SwarmSize,ParticleSize+1)); %历史全局最优值为适应度最大的解(目标值最小)，在个体极值矩阵中找 OptSwarm(SwarmSize+1,:) = OptSwarm(row,:); %带有适应度值的 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %初始化一些备用变量 %清除图形为画图做准备，这样不需要重复地关闭图形。 figure(1); clf(1); %figure(2); clf(2); %如果IsPlot为1 %figure(3); clf(3); %如果IsPlot为1 %初始化每一次迭代的最好适应度(最小目标值) MinObjV = zeros(1,LoopCount); %初始化每一次迭代的平均适应度(目标值) MeanObjV = zeros(1,LoopCount); %初始化每一次迭代的全局历史最好适应度 GlobalAdapt = zeros(1,LoopCount); %初始化每一次迭代后的离线与在线性性能 OnLine = zeros(1,LoopCount); TempOffLine = zeros(1,LoopCount); OffLine = zeros(1,LoopCount); %初始化离线及在线曲线画图标志 IsPlot = 0; %绘图为1，不绘图为0 %此处也可将向量形式的速度、位置范围变成矩阵形式，放在此处减少每次迭代的计算开销，注意不要再clear这些变量了。 %TempMinVel = repmat((ParticleVelScope(:,1))',[SwarmSize,1]); %SwarmSizeXParticleSize %TempMaxVel = repmat((ParticleVelScope(:,2))',[SwarmSize,1]); %SwarmSizeXParticleSize %TempMinPos = repmat((ParticlePosScope(:,1))',[SwarmSize,1]); %SwarmSizeXParticleSize %TempMaxPos = repmat((ParticlePosScope(:,2))',[SwarmSize,1]); %SwarmSizeXParticleSize %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %开始更新迭代算法 for CurCount = 1:LoopCount %*********************************************************************% %功能描述：Clerc(1999)年提出的带收敛因子（constrictor factor）的PSO算法 %参考文献：改进的粒子群优化算法及应用 %---------------------------------------------------------------------% %不带惯性权重因子：w w = 1; %有关位置约束因子：a a = 1; %有关认知系数c1、社会系数c2，满足(c1+c2)>4 c1 = 2.05; c2 = 2.05; %收敛因子：K psi = c1+c2; K = 2/abs(2-psi-sqrt(psi^2-4*psi)); %---------------------------------------------------------------------% %SubTract1为个体历史最优位置减当前最优位置，矩阵SwarmSizeXParticleSize SubTract1 = OptSwarm(1:SwarmSize,1:ParticleSize)-ParSwarm(:,1:ParticleSize); %扩展全局历史最优向量1XParticleSize为矩阵SwarmSizeXParticleSize以与当前最优位置矩阵同维 TempOptGlobal = repmat(OptSwarm(SwarmSize+1,1:ParticleSize),SwarmSize,1); %SubTract2为全局历史最优位置减当前最优位置,矩阵SwarmSizeXParticleSize SubTract2 = TempOptGlobal-ParSwarm(:,1:ParticleSize); clear TempOptGlobal; %更新速度矩阵，位于ParSwarm %首先产生一随机分布矩阵SwarmSizeXParticleSize且随机数对每个粒子的每一维均相同，不同粒子的随机数不同 TempRand1 = repmat(unifrnd(0,1,SwarmSize,1),[1,ParticleSize]); TempRand2 = repmat(unifrnd(0,1,SwarmSize,1),[1,ParticleSize]); TempVel = K*(w*ParSwarm(:,ParticleSize+1:2*ParticleSize)+c1*TempRand1.*SubTract1+c2*TempRand2.*SubTract2); clear TempRand1 TempRand2; %限制速度用ParticleVelScope，用最后的TempVel更新ParSwarm的速度空间 TempMinVel = repmat((ParticleVelScope(:,1))',SwarmSize,1); %SwarmSizeXParticleSize TempMaxVel = repmat((ParticleVelScope(:,2))',SwarmSize,1); %SwarmSizeXParticleSize TempVel = max(TempMinVel,TempVel); %先对每维速度最小限制 TempVel = min(TempMaxVel,TempVel); %再对每维速度最大限制 ParSwarm(:,ParticleSize+1:2*P