粒子群优化G-SFR辨识附matlab代码.zip

function [pbest best_record gbest]=PSO_Step2() %打开仿真模型 ModelName='SFRRR'; ModelFile=[ModelName '.slx']; %打开模型 open(ModelFile) %最大迭代次数 20180630 POPSIZE=30; MAXGEN=20; %模型名称 global MN; MN=ModelName; %计算步长 global Ts; Ts=0.02; global ST; ST=0; global ET; ET=100/Ts; global T; T=0:Ts:100; global dPValue; dPValue=-0.04; global A3_ok; A3_ok=-10; load('dfreq_real.mat'); global dFreqA; % dFreqA=dfreq; dFreqA=data_analysis_1(:,3); % global dFreqL; % dFreqL=tcxs120(1:find(tcxs120==min(tcxs120))); % global dFreqY; % dFreqY=tcxs120(find(fhtp150==min(fhtp150)):end); % global dFreqG; % dFreqG=tcxs120(1:length(dFreqL)+find(dFreqY==max(fhtp150(10000:end)))); %% 粒子群求最小值 format long; warning off; tstart=tic; %记录时间 global region; %% ////////////////////////////////////////////////////////////////////以下为需要修改部分（模型的一些初始化工作）////////////////////////////////////////////////// %% 待辨识参数 %% ////////////////////////////////////////////////////////////////////以上为需要修改部分（模型的一些初始化工作）////////////////////////////////////////////////// disp(['*************PSO算法开始*****************']); %% ////////////////////////////////////////////////////////////////////以下为需要修改部分（PSO算法的一些初始化工作）////////////////////////////////////////////////// %% PSO的参数设置 pop_size = POPSIZE; % pop_size 种群大小（包含有十五个粒子） max_gen = MAXGEN; % max_gen 最大迭代次数20 part_size = 5; % part_size 需要辨识参数的数量 region=[ 10 30 1 20 40 90 1 60 -300 300 ]; %% ////////////////////////////////////////////////////////////////////以上为需要修改部分（PSO算法的一些初始化工作）////////////////////////////////////////////////// gbest=zeros(1,part_size+1); % gbest当前搜索到的最小的值 rng('shuffle'); % 重新“洗牌”，产生随机数初始化 %% 初始化 arr_present = ini_pos(pop_size,part_size); % 当前位置,随机初始化,rand()的范围为0~1 v=ini_v(pop_size,part_size); % 初始化当前速度 pbest = zeros(pop_size,part_size+1); % pbest 粒子以前搜索到的最优值，最后一列包括这些值的适应度 w_max = 0.9; % w_max 权系数最大值 w_min = 0.45; c1 = 2; % 学习因子 c2 = 2; % 学习因子 best_record = zeros(max_gen,part_size+1); % best_record记录每次迭代最好的粒子的适应度。 arr_present(:,end)=ini_fit(arr_present,pop_size,part_size); % 初始化适应度。（目标函数值） pbest = arr_present; % pbest初始化各个粒子最优值 [best_value best_index] = min(arr_present(:,end)); % 初始化全局最优（此时取最小值为最优值） gbest = arr_present(best_index,:); % gbest是当前的全局最优，1*4维矩阵，前三列维参数，最后一列为目标函数值 %% 迭代 station=0; %% 变异时机初始化 %% 迭代 for i=1:max_gen gengxin=1; disp(['第' num2str(i) '轮迭代.']); gbest %************************************ try %当前的最优曲线 SFR_N=[ 1]; param=['[' num2str(SFR_N) ']']; set_param('SFRRR/SFR','Numerator',param); SFR_D=[gbest(1) gbest(2) ]; param=['[' num2str(SFR_D) ']']; set_param('SFRRR/SFR','Denominator',param); SFR1_N=[ gbest(5) A3_ok-gbest(2)]; param=['[' num2str(SFR1_N) ']']; set_param('SFRRR/SFR1','Numerator',param); SFR1_D=[gbest(3) gbest(4) 1 ]; param=['[' num2str(SFR1_D) ']']; set_param('SFRRR/SFR1','Denominator',param); %设置功率突变率 % param=['[' num2str(dPValue) ']']; % set_param('SFRRR/dP','Value',param); %运行仿真程序 sim('SFRRR',27.98); figure(3); plot(dFreqA,'b'); hold on; plot(F_SFR,'r'); hold off; catch disp('仿真出错！'); end %************************************ w = w_max-(w_max-w_min)*i/max_gen; %当前迭代的w %% 求取速度并限制速度在一定范围内 for j=1:pop_size v(j,:) = w.*v(j,:)+c1.*rand*(pbest(j,1:part_size)-arr_present(j,1:part_size))+c2.*rand*(gbest(1:part_size)-arr_present(j,1:part_size)); % 粒子速度更新 (a) %% 速度大小限制在运动范围的20% [mm,nn]=size(region); for aa=1:mm if abs(v(j,aa)) > (region(aa,2)-region(aa,1))*0.2 v(j,aa)=sign(v(j,aa))*(region(aa,2)-region(aa,1))*0.2; end end %% 更新迭代 arr_present(j,1:part_size) = arr_present(j,1:part_size)+v(j,1:part_size); % 粒子位置更新 (b) %% 变异 if station>=8 %适应度结果不好误差过大且多次没有更新最优值 if j~=best_index&&rand<1/(1+sqrt(i)) arr_present(j,1:part_size)= arr_present(j,1:part_size)+randn(1,part_size).*(region(:,2)-region(:,1))'*0.1; disp('变异'); end end %% 判断是否越界更新适应度 [ arr_present(j,1:part_size) v(j,:)]=Region_in(arr_present(j,1:part_size),v(j,:)); arr_present(j,end) = fitness(arr_present(j,1:part_size)); % 粒子对应的目标函数更新 if arr_present(j,end)<pbest(j,end) % 根据条件更新pbest,找到更新的最小值且没有越界的参数值，就更新最小值 pbest(j,:) = arr_present(j,:); end end [best best_index] = min(arr_present(:,end)); % 求取当前最小值，如果是最小的值为min if best<gbest(end) % 如果当前最好的结果比以前的好，则更新最优值gbest,如果是最小的值为小于号,相反则为大于号 gbest = arr_present(best_index,:); gengxin=0; end %% 计算变异时机 if gengxin==1 %没更新 station=station+1; else station=0; % 变异时机清零 end best_record(i,:) = gbest; %记录每一次迭代的最优值（目标函数值） end %% 输出历史最优值、最后一次迭代最优位置以及目标函数值、每次迭代过程的最优目标值。 %显示最终结果 disp('最终结果') gbest %************************************ try %当前的最优曲线 SFR_N=[ 1]; param=['[' num2str(SFR_N) ']']; set_param('SFRRR/SFR','Numerator',param); SFR_D=[gbest(1) gbest(2) ]; param=['[' num2str(SFR_D) ']']; set_param('SFRRR/SFR','Denominator',param); SFR1_N=[ gbest(5) A3_ok-gbest(2)]; param=['[' num2str(SFR1_N) ']']; set_param('SFRRR/SFR1','Numerator',param); SFR1_D=[gbest(3) gbest(4) 1 ]; param=['[' num2str(SFR1_D) ']']; set_param('SFRRR/SFR1','Denominator',param); %运行仿真程序 sim('SFRRR',27.98); figure(4); plot(dFreqA,'b'); hold on; plot(F_SFR,'r'); hold off; catch disp('仿真出错！'); end %************************************ %% 子函数定义 %% ////////////////////////////////////////////////////////////////////以下为需要修改部分////////////////////////////////////////////////// %% 计算适应度 function fit = fitness(present) %模型名称 global Ts; global ST; global ET; global t; t=0:Ts:100; global A3_ok; global dFreqA; global dFreqL; global dPValue; global dFreqG; global dFreqY; % present(1): B0 % present(2): B1 % present(3): A0 % present(4): A1 % present(5): A2 SFR_N=[ 1]; param=['[' num2str(SFR_N) ']']; set_param('SFRRR/SFR','Numerator',param); SFR_D=[present(1) present(2) ]; param=['[' num2str(SFR_D) ']']; set_param('SFRRR/SFR','Denomina