1 内容介绍
经济负载调度问题是在满足约束的同时以最小的成本将负载分配给工厂。它被制定为一个优化问题，即在满足需求和损失的同时，最大限度地降低所有承诺工厂的总燃料成本。这些问题的变体很多，它们以不同的方式模拟目标和约束。
基本经济调度问题在数学上可以描述为最小化受约束的所有承诺工厂的总燃料成本的最小化问题。
2 部分代码
% pso_Trelea_vectorized.m
% 通用粒子群优化器
% 以查找任何项的最小值或最大值
% 味精矩阵函数
%
% 使用率：
% [选择选项]=PSO（函数名，D）
% 或：
% [选择输出，tr，te]=...
% PSO（functname，D，mv，VarRange，minmax，PSOparams，plotfcn，PSOseedValue）
%
% 输入：
% 函子名 - 要优化的 matlab 函数字符串
% D - 函数输入的 # （问题的维度）
%
% 可选输入：
% mv - 最大粒子速度，标量或长度为 D 的向量
%（这允许每个组件具有自己的最大速度），
% 默认值 = 4，如果未输入，则设置为 NaN 或输入
%
% 变量范围 - 每个输入变量的范围矩阵，
% 默认值 -100 到 100，格式为：
% [ 最小值1 最大值1
% 最小值2 最大值2
% ...
% 最小值 最大值 ]
%
% 最小值 = 0，函数最小化（默认值）
% = 1，函数最大化
% = 2，函数以 P（12） 为目标（最小化到遍历的距离）
% PS蛋白酶 - PSO 参数
% P（1） - 更新显示之间的时间，默认值 = 100。如果为 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,
默认值 % = 150（如果 SSE 在 250 个纪元内未更改）
% 然后退出
% P（11）- 误差目标，如果 NaN 则无约束最小值或最大值，则默认值为 NaN
% P（12）- 类型标志（使用哪种 PSO）
% 0 = 带间质的普通 PSO（默认）
% 1，2 = 特雷利类型 1，2
% 3 = 克莱尔克收缩 PSO，类型 1”
% P（13）- PS 已确定，默认值为 0
% = 0 表示初始仓位全部随机
% = 1 表示作为用户输入的初始粒子
%
% 绘图fcn - 绘图函数的可选名称，默认的“主题”，
%自己做，放在这里
%
% PS已摄取值 - 初始粒子位置，取决于 P（13），必须为
如果 P（13） 是 1 或 2，不用于 P（13）=0，则需要
% 为 nXm，其中 n<= ps，并且 m<= D
% 如果 n<ps 和/或 m<D，则剩余值是随机设置的
% 在瓦朗日上
% 输出：
% optOUT - 最佳输入和关联的最小/最大输出功能，形式为：
% [ bestin1
% bestin2
如果存在（“PSOSE值”）
tmpsz=size（PSOSEValue）;
如果 D <（2）
错误（“PSOSE值列大小必须为D或更小”）;
结束
如果 ps < tmpsz（1）
错误（“PSOSE值行长度必须为粒子数或更少”）;
结束
结束
% 设置绘图标志
如果 （P（1））~=0
情节flg=1;
还
情节流 =0;
结束
% 预分配变量以加快速度
tr = 一（1，me）*纳纳;
% 负责在此处设置最大速度和位置参数
如果长度（mv）==1
velmaskmin = -mv*一% 最小值，psXD 矩阵
velmaskmax = mv*一（ps，D）;% 最大垂直度
elseif 长度（mv）==D
if （迭代贝斯特-错误）^2 <= （gbestval-ergoal）^2
格贝斯特瓦尔 = 伊特贝斯特瓦尔;
gbest = pbest（idx1，:);
% 与训练方案一起使用，用于神经网络训练
% 在每次迭代时将 gbest 分配给净值，这些临时分配
% 主要用于绘图
if strcmp（functname，'pso_neteval（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（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（(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（gbestoffset,ps,1);
% pbest = pbest - repmat(gbestoffset,ps,1);
% gbest = gbest + gbestoffset;
end % end epoch loop
%% clear temp outputs
% eval（find(~isnan(tr)))];
return
3 运行结果
4 参考文献
[1]芮钧, 陈守伦. MATLAB粒子群算法工具箱求解水电站优化调度问题[J]. 中国农村水利水电, 2009(1):3.
[2]解玖霞. 基于粒子群算法的主动配电网经济优化调度[J]. 电气技术, 2017(9):4.