【优化调度】基于粒子群算法解决经济调度附Matlab代码
1 内容介绍
Economic load dispatch problem is allocating loads to plants for minimum cost while meeting the constraints. It is formulated as an optimization problem of minimizing the total fuel cost of all committed plant while meeting the demand and losses .The variants of the problems are numerous which model the objective and the constraints in different ways.
The basic economic dispatch problem can described mathematically as a minimization of problem of minimizing the total fuel cost of all committed plants subject to the constraints.
2 部分代码
% pso_Trelea_vectorized.m
% a generic particle swarm optimizer
% to find the minimum or maximum of any
% MISO matlab function
%
% Usage:
% [optOUT]=PSO(functname,D)
% or:
% [optOUT,tr,te]=...
% PSO(functname,D,mv,VarRange,minmax,PSOparams,plotfcn,PSOseedValue)
%
% Inputs:
% functname - string of matlab function to optimize
% D - # of inputs to the function (dimension of problem)
%
% Optional Inputs:
% mv - max particle velocity, either a scalar or a vector of length D
% (this allows each component to have it's own max velocity),
% default = 4, set if not input or input as NaN
%
% VarRange - matrix of ranges for each input variable,
% default -100 to 100, of form:
% [ min1 max1
% min2 max2
% ...
% minD maxD ]
%
% minmax = 0, funct minimized (default)
% = 1, funct maximized
% = 2, funct is targeted to P(12) (minimizes distance to errgoal)
% PSOparams - PSO parameters
% P(1) - Epochs between updating display, default = 100. if 0,
% no display
% P(2) - Maximum number of iterations (epochs) to train, default = 2000.
% P(3) - population size, default = 24
%
% P(4) - acceleration const 1 (local best influence), default = 2
% P(5) - acceleration const 2 (global best influence), default = 2
% P(6) - Initial inertia weight, default = 0.9
% P(7) - Final inertia weight, default = 0.4
% P(8) - Epoch when inertial weight at final value, default = 1500
% P(9)- minimum global error gradient,
% if abs(Gbest(i+1)-Gbest(i)) < gradient over
% certain length of epochs, terminate run, default = 1e-25
% P(10)- epochs before error gradient criterion terminates run,
% default = 150, if the SSE does not change over 250 epochs
% then exit
% P(11)- error goal, if NaN then unconstrained min or max, default=NaN
% P(12)- type flag (which kind of PSO to use)
% 0 = Common PSO w/intertia (default)
% 1,2 = Trelea types 1,2
% 3 = Clerc's Constricted PSO, Type 1"
% P(13)- PSOseed, default=0
% = 0 for initial positions all random
% = 1 for initial particles as user input
%
% plotfcn - optional name of plotting function, default 'goplotpso',
% make your own and put here
%
% PSOseedValue - initial particle position, depends on P(13), must be
% set if P(13) is 1 or 2, not used for P(13)=0, needs to
% be nXm where n<=ps, and m<=D
% If n<ps and/or m<D then remaining values are set random
% on Varrange
% Outputs:
% optOUT - optimal inputs and associated min/max output of function, of form:
% [ bestin1
% bestin2
if exist('PSOseedValue')
tmpsz=size(PSOseedValue);
if D < tmpsz(2)
error('PSOseedValue column size must be D or less');
end
if ps < tmpsz(1)
error('PSOseedValue row length must be # of particles or less');
end
end
% set plotting flag
if (P(1))~=0
plotflg=1;
else
plotflg=0;
end
% preallocate variables for speed up
tr = ones(1,me)*NaN;
% take care of setting max velocity and position params here
if length(mv)==1
velmaskmin = -mv*ones(ps,D); % min vel, psXD matrix
velmaskmax = mv*ones(ps,D); % max vel
elseif length(mv)==D
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
3 运行结果
4 参考文献
[1]芮钧, 陈守伦. MATLAB粒子群算法工具箱求解水电站优化调度问题[J]. 中国农村水利水电, 2009(1):3.
[2]解玖霞. 基于粒子群算法的主动配电网经济优化调度[J]. 电气技术, 2017(9):4.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机、雷达通信、无线传感器等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。