% INITIALIZE END INITIALIZE END INITIALIZE END INITIALIZE END
rstflg = 0; % for dynamic environment checking
% start PSO iterative procedures
cnt = 0; % counter used for updating display according to df in the options
cnt2 = 0; % counter used for the stopping subroutine based on error convergence
iwt(1) = iw1;
for i=1:me % start epoch loop (iterations)
out = feval(functname,[pos;gbest]);
outbestval = out(end,:);
out = out(1:end-1,:);
tr(i+1) = gbestval; % keep track of global best val
te = i; % returns epoch number to calling program when done
bestpos(i,1:D+1) = [gbest,gbestval];
%assignin('base','bestpos',bestpos(i,1:D+1));
%------------------------------------------------------------------------
% this section does the plots during iterations
if plotflg==1
if (rem(i,df) == 0 ) | (i==me) | (i==1)
fprintf(message,i,gbestval);
cnt = cnt+1; % count how many times we display (useful for movies)
eval(plotfcn); % defined at top of script
end % end update display every df if statement
end % end plotflg if statement
% check for an error space that changes wrt time/iter
% threshold value that determines dynamic environment
% sees if the value of gbest changes more than some threshold value
% for the same location
chkdyn = 1;
rstflg = 0; % for dynamic environment checking
if chkdyn==1
threshld = 0.05; % percent current best is allowed to change, .05 = 5% etc
letiter = 5; % # of iterations before checking environment, leave at least 3 so PSO has time to converge
outorng = abs( 1- (outbestval/gbestval) ) >= threshld;
samepos = (max( sentry == gbest ));
if (outorng & samepos) & rem(i,letiter)==0
rstflg=1;
% disp('New Environment: reset pbest, gbest, and vel');
%% reset pbest and pbestval if warranted
% outpbestval = feval( functname,[pbest] );
% Poutorng = abs( 1-(outpbestval./pbestval) ) > threshld;
% pbestval = pbestval.*~Poutorng + outpbestval.*Poutorng;
% pbest = pbest.*repmat(~Poutorng,1,D) + pos.*repmat(Poutorng,1,D);
pbest = pos; % reset personal bests to current positions
pbestval = out;
vel = vel*10; % agitate particles a little (or a lot)
% recalculate best vals
if minmax == 1
[gbestval,idx1] = max(pbestval);
elseif minmax==0
[gbestval,idx1] = min(pbestval);
elseif minmax==2 % this section needs work
[temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2);
gbestval = pbestval(idx1);
end
gbest = pbest(idx1,:);
% used with trainpso, for neural net training
% assign gbest to net at each iteration, these interim assignments
% are for plotting mostly
if strcmp(functname,'pso_neteval')
net=setx(net,gbest);
end
end % end if outorng
sentryval = gbestval;
sentry = gbest;
end % end if chkdyn
% find particles where we have new pbest, depending on minmax choice
% then find gbest and gbestval
%[size(out),size(pbestval)]
if rstflg == 0
if minmax == 0
[tempi] = find(pbestval>=out); % new min pbestvals
pbestval(tempi,1) = out(tempi); % update pbestvals
pbest(tempi,:) = pos(tempi,:); % update pbest positions
[iterbestval,idx1] = min(pbestval);
if gbestval >= iterbestval
gbestval = iterbestval;
gbest = pbest(idx1,:);
% used with trainpso, for neural net training
% assign gbest to net at each iteration, these interim assignments
% are for plotting mostly
if strcmp(functname,'pso_neteval')
net=setx(net,gbest);
end
end
elseif minmax == 1
[tempi,dum] = find(pbestval<=out); % new max pbestvals
pbestval(tempi,1) = out(tempi,1); % update pbestvals
pbest(tempi,:) = pos(tempi,:); % update pbest positions
[iterbestval,idx1] = max(pbestval);
if gbestval <= iterbestval
gbestval = iterbestval;
gbest = pbest(idx1,:);
% used with trainpso, for neural net training
% assign gbest to net at each iteration, these interim assignments
% are for plotting mostly
if strcmp(functname,'pso_neteval')
net=setx(net,gbest);
end
end
elseif minmax == 2 % this won't work as it is, fix it later
egones = errgoal*ones(ps,1); % vector of errgoals
sqrerr2 = ((pbestval-egones).^2);
sqrerr1 = ((out-egones).^2);
[tempi,dum] = find(sqerr1 <= sqrerr2); % find particles closest to targ
pbestval(tempi,1) = out(tempi,1); % update pbestvals
pbest(tempi,:) = pos(tempi,:); % update pbest positions
sqrerr = ((pbestval-egones).^2); % need to do this to reflect new pbests
[temp,idx1] = min(sqrerr);
iterbestval = pbestval(idx1);
if (iterbestval-errgoal)^2 <= (gbestval-errgoal)^2
gbestval = iterbestval;
gbest = pbest(idx1,:);
% used with trainpso, for neural net training
% assign gbest to net at each iteration, these interim assignments
% are for plotting mostly
if strcmp(functname,'pso_neteval')
net=setx(net,gbest);
end
end
end
end
% % build a simple predictor 10th order, for gbest trajectory
% if i>500
% for dimcnt=1:D
% pred_coef = polyfit(i-250:i,(bestpos(i-250:i,dimcnt))',20);
% % pred_coef = polyfit(200:i,(bestpos(200:i,dimcnt))',20);
% gbest_pred(i,dimcnt) = polyval(pred_coef,i+1);
% end
% else
% gbest_pred(i,:) = zeros(size(gbest));
% end
%gbest_pred(i,:)=gbest;
%assignin('base','gbest_pred',gbest_pred);
% % convert to non-inertial frame
% gbestoffset = gbest - gbest_pred(i,:);
% gbest = gbest - gbestoffset;
% pos = pos + repmat(gbestoffset,ps,1);
% pbest = pbest + repmat(gbestoffset,ps,1);
%PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO
% get new velocities, positions (this is the heart of the PSO algorithm)
% each epoch get new set of random numbers
rannum1 = rand([ps,D]); % for Trelea and Clerc types
rannum2 = rand([ps,D]);
if trelea == 2
% from Trelea's paper, parameter set 2
vel = 0.729.*vel... % prev vel
+1.494.*rannum1.*(pbest-pos)... % independent
+1.494.*rannum2.*(repmat(gbest,ps,1)-pos); % social
elseif trelea == 1
% from Trelea's paper, parameter set 1
vel = 0.600.*vel... % prev vel
+1.700.*rannum1.*(pbest-pos)... % independent
+1.700.*rannum2.*(repmat(gbest,ps,1)-pos); % social
elseif trelea ==3
% Clerc's Type 1" PSO
vel = chi*(vel... % prev vel
+ac1.*rannum1.*(pbest-pos)... % independent
+ac2.*rannum2.*(repmat(gbest,ps,1)-pos)) ; % social
else
% common PSO algo with inertia wt
% get inertia weight, just a linear funct w.r.t. epoch parameter iwe
if i<=iwe
iwt(i) = ((iw2-iw1)/(iwe-1))*(i-1)+iw1;
else
iwt(i) = iw2;
end
% random number including acceleration constants
ac11 = rannum1.*ac1; % for common PSO w/inertia
ac22 = rannum2.*ac2;
vel = iwt(i).*vel... % prev vel
+ac11.*(pbest-pos)... % independent
+ac22.*(repmat(gbest,ps,1)-pos); % social
end
% limit velocities here using masking
vel = ( (vel <= velmaskmin).*velmaskmin ) + ( (vel > velmaskmin).*vel );
vel = ( (vel >= velmaskmax).*velmaskmax ) + ( (vel < velmaskmax).*vel );
% update new position (PSO algo)
pos = pos + vel;
% position masking, limits positions to desired search space
% method: 0) no position limiting, 1) saturation at limit,
% 2) wraparound at limit , 3) bounce off limit
minposmask_throwaway = pos <= posmaskmin; % these are psXD matrices
minposmask_keep = pos > posmaskmin;
maxposmask_throwaway = pos >= posmaskmax;
maxposmask_keep = pos
if posmaskmeth == 1
% this is the saturation method
pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos );
pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos );
elseif posmaskmeth == 2
% this is the wraparound method
pos = ( minposmask_throwaway.*posmaskmax ) + ( minposmask_keep.*pos );
pos = ( maxposmask_throwaway.*posmaskmin ) + ( maxposmask_keep.*pos );
elseif posmaskmeth == 3
% this is the bounce method, particles bounce off the boundaries with -vel
pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos );
pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos );
vel = (vel.*minposmask_keep) + (-vel.*minposmask_throwaway);
vel = (vel.*maxposmask_keep) + (-vel.*maxposmask_throwaway);
else
% no change, this is the original Eberhart, Kennedy method,
% it lets the particles grow beyond bounds if psoparams (P)
% especially Vmax, aren't set correctly, see the literature
end
%PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO
% check for stopping criterion based on speed of convergence to desired
% error
tmp1 = abs(tr(i) - gbestval);
if tmp1 > ergrd
cnt2 = 0;
elseif tmp1 <= ergrd
cnt2 = cnt2+1;
if cnt2 >= ergrdep
if plotflg == 1
fprintf(message,i,gbestval);
disp(' ');
disp(['--> Solution likely, GBest hasn''t changed by at least ',...
num2str(ergrd),' for ',...
num2str(cnt2),' epochs.']);
eval(plotfcn);
end
break
end
end
% this stops if using constrained optimization and goal is reached
if ~isnan(errgoal)
if ((gbestval<=errgoal) & (minmax==0)) | ((gbestval>=errgoal) & (minmax==1))
if plotflg == 1
fprintf(message,i,gbestval);
disp(' ');
disp(['--> Error Goal reached, successful termination!']);
eval(plotfcn);
end
break
end
% this is stopping criterion for constrained from both sides
if minmax == 2
if ((tr(i)=errgoal)) | ((tr(i)>errgoal) ...
& (gbestval <= errgoal))
if plotflg == 1
fprintf(message,i,gbestval);
disp(' ');
disp(['--> Error Goal reached, successful termination!']);
eval(plotfcn);
end
break
end
end % end if minmax==2
end % end ~isnan if
% % convert back to inertial frame
% pos = pos - repmat(gbestoffset,ps,1);
% pbest = pbest - repmat(gbestoffset,ps,1);
% gbest = gbest + gbestoffset;
end % end epoch loop
扫一扫
1223
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
