【优化调度】基于粒子群算法解决经济调度附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代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。