psi     = ac1 + ac2;     chi_den = abs(2-psi-sqrt(psi^2 - 4*psi));     chi_num = 2*kappa;     chi     = chi_num/chi_den; endend% INITIALIZE END INITIALIZE END INITIALIZE END INITIALIZE ENDrstflg = 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 <  posmaskmax;             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) & (gbestval>=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%% clear temp outputs% evalin('base','clear temp_pso_out temp_te temp_tr;');% output & return OUT=[gbest';gbestval]; varargout{1}=[1:te]; varargout{2}=[tr(find(~isnan(tr)))];  return