PSO.rar_PSO c++_粒子群算法

clear all, clc % pso example iter = 1000; % number of algorithm iterations np = 2; % number of model parameters ns = 10; % number of sets of model parameters Wmax = 0.9; % maximum inertial weight Wmin = 0.4; % minimum inertial weight c1 = 2.0; % parameter in PSO methodology c2 = 2.0; % parameter in PSO methodology Pmax = [10 10]; % maximum model parameter value Pmin = [-10 -10]; % minimum model parameter value Vmax = [1 1]; % maximum change in model parameter Vmin = [-1 -1]; % minimum change in model parameter modelparameters(1:np,1:ns) = 0; % set all model parameter estimates for all model parameter sets to zero modelparameterchanges(1:np,1:ns) = 0; % set all change in model parameter estimates for all model parameter sets to zero bestmodelparameters(1:np,1:ns) = 0; % set best model parameter estimates for all model parameter sets to zero setbestcostfunction(1:ns) = 1e6; % set best cost function of each model parameter set to a large number globalbestparameters(1:np) = 0; % set best model parameter values for all model parameter sets to zero bestparameters = globalbestparameters'; % best model parameter values for all model parameter sets (to plot) globalbestcostfunction = 1e6; % set best cost function for all model parameter sets to a large number i = 0; % indicates ith algorithm iteration j = 0; % indicates jth set of model parameters k = 0; % indicates kth model parameter for k = 1:np % initialization for j = 1:ns modelparameters(k,j) = (Pmax(k)-Pmin(k))*rand(1) + Pmin(k); % randomly distribute model parameters modelparameterchanges(k,j) = (Vmax(k)-Vmin(k))*rand(1) + Vmin(k); % randomly distribute change in model parameters end end for i = 2:iter for j = 1:ns x = modelparameters(:,j); % calculate cost function costfunction = 105*(x(2)-x(1)^2)^2 + (1-x(1))^2; if costfunction < setbestcostfunction(j) % best cost function for jth set of model parameters bestmodelparameters(:,j) = modelparameters(:,j); setbestcostfunction(j) = costfunction; end if costfunction < globalbestcostfunction % best cost function for all sets of model parameters globalbestparameters = modelparameters(:,j); bestparameters(:,i) = globalbestparameters; globalbestcostfunction(i) = costfunction; else bestparameters(:,i) = bestparameters(:,i-1); globalbestcostfunction(i) = globalbestcostfunction(i-1); end end W = Wmax - i*(Wmax-Wmin)/iter; % compute inertial weight for j = 1:ns % update change in model parameters and model parameters for k = 1:np modelparameterchanges(k,j) = W*modelparameterchanges(k,j) + c1*rand(1)*(bestmodelparameters(k,j)-modelparameters(k,j))... + c2*rand(1)*(globalbestparameters(k) - modelparameters(k,j)); if modelparameterchanges(k,j) < -Vmax(k), modelparameters(k,j) = modelparameters(k,j) - Vmax(k); end if modelparameterchanges(k,j) > Vmax(k), modelparameters(k,j) = modelparameters(k,j) + Vmax(k); end if modelparameterchanges(k,j) > -Vmax(k) & modelparameterchanges(k,j) < Vmax(k), modelparameters(k,j) = modelparameters(k,j) + modelparameterchanges(k,j); end if modelparameters(k,j) < Pmin(k), modelparameters(k,j) = Pmin(k); end if modelparameters(k,j) > Pmax(k), modelparameters(k,j) = Pmax(k); end end end i end bp = bestparameters; index = linspace(1,iter,iter); figure; semilogy(globalbestcostfunction,'k'); set(gca,'FontName','Arial','Fontsize',14); axis tight; xlabel('iteration'); ylabel('cost function'); figure; q = plot(index,bp(1,,'k-',index,bp(2,,'k:'); set(gca,'FontName','Arial','Fontsize',14); axis tight; legend(q,'x_1','x_2'); xlabel('iteration'); ylabel('parameter')