【微电网以后】基于改进二进制粒子群优化算法解决热电联产机组组合问题附matlab代码
✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。
🔥 内容介绍
在全球能源转型的浪潮中，微电网正崭露头角，成为能源领域的焦点。它作为一种小型、模块化的电力系统，集成了分布式电源、储能系统、电力电子装置以及负荷，宛如一个微型的能源生态系统，具备 “并网” 和 “离网” 两种运行模式，可根据实际需求灵活切换 ，实现自我控制与能源自治。
微电网在能源体系中的重要地位愈发凸显。一方面，它是促进可再生能源消纳的关键一环。太阳能、风能等可再生能源虽清洁环保，但其间歇性和不稳定性给传统电网的接纳带来挑战。微电网通过搭配储能系统和智能控制系统，能够有效平滑可再生能源的出力波动，实现高比例的本地消纳，减少对主电网的冲击，推动能源向绿色低碳转型，堪称能源转型的 “理想试验田”。另一方面，微电网显著提升了供电可靠性。当极端天气、自然灾害或大电网故障导致大面积停电时，微电网可迅速切换至离网模式，成为保障区域内关键负荷不间断供电的 “生命线能源孤岛”。此外，它还能实现能源的梯级利用，提高能源利用效率，带来可观的经济效益，如在电费高的时段，优先使用自发清洁电力或储能放电，降低用电成本，甚至可将富余电力卖给电网获取收益。
在微电网的诸多组成部分与运行环节中，热电联产机组扮演着重要角色。热电联产机组能够同时生产电能和热能，极大地提高了能源的综合利用效率，相比传统的发电和供热分离模式，优势明显。然而，要充分发挥热电联产机组的效能，实现其最优运行，机组组合问题亟待解决。所谓机组组合问题，就是要在满足电力和热力负荷需求、设备运行约束等条件下，确定各个热电联产机组在不同时段的启停状态和出力水平，以达到最小化运行成本、最大化能源利用效率等目标。这一问题涉及众多变量和复杂约束，求解难度较大。
为攻克热电联产机组组合问题，各种优化算法不断涌现，其中改进二进制粒子群优化算法脱颖而出。它源于对标准粒子群优化算法的改进，采用二进制编码来表示搜索空间中的解，巧妙地解决了实数编码中的精度问题，在组合优化等领域展现出独特优势，为解决热电联产机组组合问题提供了新的思路与有效途径。
热电联产机组组合问题剖析
（一）热电联产的运作原理
热电联产，英文缩写为 CHP（Combined Heat and Power），是一种先进的能源综合利用技术，它打破了传统的发电与供热相互分离的模式，实现了在同一生产过程中同时产出电能和热能 。
从原理上讲，热电联产主要基于能量的梯级利用。以常见的蒸汽轮机热电联产系统为例，燃料（如煤炭、天然气、生物质等）在锅炉中充分燃烧，释放出大量的热能，将水加热转化为高温高压的蒸汽。这些蒸汽首先推动汽轮机高速旋转，进而带动发电机运转，完成将热能转化为电能的关键步骤。在蒸汽做功发电后，其压力和温度虽有所降低，但仍蕴含着可观的能量。通过科学合理的系统设计，这些低品位的蒸汽或热水会被回收利用，用于满足周边区域的供暖需求，比如为居民住宅、商业建筑提供暖气；或是在工业生产中，为一些对温度要求不太高的工艺流程提供热能，像食品加工、纺织印染等行业，需要利用蒸汽进行加热、烘干等操作。这种对发电余热的二次利用，避免了能源的白白浪费，大大提高了能源的综合利用效率。据统计，传统的火力发电方式，能源利用效率通常仅在 35% - 40%，而热电联产系统凭借余热回收机制，综合能源利用效率能够达到 70% - 90%，实现了能源利用效率的大幅提升 。
从能量流向的角度，热电联产可分为先发电式和后发电式。先发电式中，蒸汽率先用于发电，之后再将剩余的蒸汽热能投入工业制程，适合食品加工、化工等中等温度制程需求的产业；后发电式则是蒸汽先服务于工业制程，完成制程后的余热再用于发电，契合造纸、冶金等高温需求场景 。
（二）机组组合问题的挑战
在微电网中，确定热电联产机组的最优组合并非易事，机组组合问题存在诸多挑战。
首先是经济成本方面。在满足热、电负荷需求的过程中，不同类型的热电联产机组有着各异的运行成本，包括燃料费用、设备维护成本以及启停成本等。燃料成本会随着燃料价格的波动而变化，且不同燃料的价格和燃烧效率不同，如天然气价格相对稳定但成本较高，煤炭价格波动大但燃烧效率较低，选择何种燃料以及如何搭配使用，都对总成本产生影响。设备维护成本与机组的运行时长、使用频率和技术复杂度相关，复杂先进的机组维护成本通常较高。启停成本也不容忽视，频繁启停机组会增加设备损耗和能源消耗，产生额外费用。以某工业园区的热电联产系统为例，在夏季制冷需求高、冬季供暖需求高的不同时段，若不合理安排机组的启停和出力，每月的运行成本可能会增加数十万元。
能源效率也是一个重要挑战。热电联产机组的能源转换效率并非固定不变，而是会受到多种因素干扰。机组的运行工况对效率影响显著，在低负荷运行时，能源转换效率会大幅下降。例如，当热电联产机组的负荷率低于 50% 时，其能源效率可能会降低 10% - 20%。此外，设备老化也会导致性能衰退，使得能源利用率逐渐降低。像运行年限超过 10 年的热电联产机组，与新机组相比，能源转换效率可能会下降 5% - 10% 。
运行稳定性是机组组合问题面临的又一挑战。热、电负荷需求具有动态变化的特性，难以精准预测。在实际运行中，可能会出现突发情况，如极端天气导致热负荷或电负荷急剧增加，或者部分机组突发故障。当遭遇这些情况时，要保障系统的稳定运行，就需要合理安排机组的组合，确保有足够的备用容量来应对负荷的波动和机组的故障。但在实际操作中，由于热电联产机组之间存在复杂的耦合关系，调整一台机组的运行状态，可能会对其他机组乃至整个系统的稳定性产生连锁反应，使得运行稳定性的保障难度加大。
⛳️ 运行结果
📣 部分代码
%%%%%%Unit commitment schedule generator for swarm based metaheuristics%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MAIN%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Coded by Vick @NSYSU%%%
clc;
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%TEST CASE I-II-V-VI%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%IEEE 10 UNIT-SYSTEM DATA
% %Pmax %%Pmin %%ai %%bi %%ci %MUT %MDT %Shot %Scold %CSH %I.S %%D/URR %%%MR %%%UN
I=[ 455 150 0.00048 16.19 1000 8 8 4500 9000 5 8 91 0 0 ; %U1
455 150 0.00031 17.26 970 8 8 5000 10000 5 8 91 0 0 ; %U2
130 20 0.002 16.60 700 5 5 550 1100 4 -5 26 0 0 ; %U3
130 20 0.00211 16.50 680 5 5 560 1120 4 -5 26 0 0 ; %U4
162 25 0.00398 19.70 450 6 6 900 1800 4 -6 33 0 0 ; %U5
80 20 0.00712 22.26 370 3 3 170 340 2 -3 21 0 0 ; %U6
85 25 0.00079 27.74 480 3 3 260 520 2 -3 20 0 0 ; %U7
55 10 0.00413 25.92 660 1 1 30 60 0 -1 11 0 0 ; %U8
55 10 0.00222 27.27 665 1 1 30 60 0 -1 11 0 0 ; %U9
55 10 0.00173 27.79 670 1 1 30 60 0 -1 11 0 0 ]; %U10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%TEST CASE III%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%IEEE 10 UNIT-SYSTEM DATA CONSIDERING EMISSION COEFFICIENTS (ton/h)
% %Pmax %Pmin %alpha %beta %gamma %MUT %MDT %Shot %Scold %CSH %I.S %%D/URR %%%MR %%%UN
% I=[ 455 150 0.00312 -0.24444 10.33908 8 8 4500 9000 5 8 0 0 0 ; %U1
% 455 150 0.00312 -0.24444 10.33908 8 8 5000 10000 5 8 0 0 0 ; %U2
% 130 20 0.00509 -0.40695 30.03910 5 5 550 1100 4 -5 0 0 0 ; %U3
% 130 20 0.00509 -0.40695 30.03910 5 5 560 1120 4 -5 0 0 0 ; %U4
% 162 25 0.00344 -0.38132 32.00006 6 6 900 1800 4 -6 0 0 0 ; %U5
% 80 20 0.00344 -0.38132 32.00006 3 3 170 340 2 -3 0 0 0 ; %U6
% 85 25 0.00465 -0.39023 33.00056 3 3 260 520 2 -3 0 0 0 ; %U7
% 55 10 0.00465 -0.39023 33.00056 1 1 30 60 0 -1 0 0 0 ; %U8
% 55 10 0.00465 -0.39524 35.00056 1 1 30 60 0 -1 0 0 0 ; %U9
% 55 10 0.00470 -0.39864 36.00012 1 1 30 60 0 -1 0 0 0 ]; %U10
% % % %%%%%%%%%%%%%%%%%%%DEMAND DATA FOR 10 UNIT SYSTEM%%%%%%%%%%%%%%%%%%%%
% % Demand & Spinning reserve requirement per hour (CASES I,III,IV,V)
P_D=[700;750;850;950;1000;1100;1150;1200;1300;1400;1450;1500;1400;1300;1200;1050;1000;1100;1200;1400;1300;1100;900;800];
% %
% % % Demand & Spinning reserve requirement per hour (CASE VI-DRUC)
% % P_D=[700;750;850;950;1000;1100;1150;1200;1040;1120;1160;1200;1120;1040;1200;1050;1000;1100;1200;1120;1040;1100;900;800];
%
% %%%Common for all cases
SRREQ=0.10*P_D;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%TEST CASE II%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%IEEE 26 UNIT-SYSTEM DATA (ROBUSTNESS TEST)
% %Pmax %%Pmin %%ai %%bi %%ci %MUT %MDT %Shot %Scold %CSH %I.S %%D/URR %%%MR %%%UN
% I=[ 400 100 0.0019 7.5031 311.9102 8 5 500 500 10 10 0 0 0 ; %U1
% 400 100 0.0019 7.4921 310.0021 8 5 500 500 10 10 0 0 0 ; %U2
% 350 140 0.0015 10.8616 177.0575 8 5 300 200 8 10 0 0 0 ; %U3
% 197 68.95 0.0026 23.2000 260.1760 5 4 200 200 8 -4 0 0 0 ; %U4
% 197 68.95 0.0026 23.1000 259.6490 5 4 200 200 8 -4 0 0 0 ; %U5
% 197 68.95 0.0026 23.0000 259.1310 5 4 200 200 8 -4 0 0 0 ; %U6
% 155 54.25 0.0049 10.7583 143.5972 5 3 150 150 6 5 0 0 0 ; %U7
% 155 54.25 0.0048 10.7367 143.3719 5 3 150 150 6 5 0 0 0 ; %U8
% 155 54.25 0.0047 10.7154 143.0288 5 3 150 150 6 5 0 0 0 ; %U9
% 155 54.25 0.0046 10.6940 142.7348 5 3 150 150 6 5 0 0 0 ; %U10
% 100 25.00 0.0060 18.2000 218.7752 4 2 70 70 4 -3 0 0 0 ; %U11
% 100 25.00 0.0061 18.1000 218.3350 4 2 70 70 4 -3 0 0 0 ; %U12
% 100 25.00 0.0062 18.0000 217.8952 4 2 70 70 4 -3 0 0 0 ; %U13
% 76 15.20 0.0093 13.4073 81.6259 3 2 50 50 3 3 0 0 0 ; %U14
% 76 15.20 0.0091 13.3805 81.4641 3 2 50 50 3 3 0 0 0 ; %U15
% 76 15.20 0.0089 13.3538 81.2980 3 2 50 50 3 3 0 0 0 ; %U16
% 76 15.20 0.0088 13.3272 81.1364 3 2 50 50 3 3 0 0 0 ; %U17
% 20 4.000 0.0143 37.8896 118.8206 1 1 20 20 2 -1 0 0 0 ; %U18
% 20 4.000 0.0136 37.7770 118.4576 1 1 20 20 2 -1 0 0 0 ; %U19
% 20 4.000 0.0126 37.6637 118.1083 1 1 20 20 2 -1 0 0 0 ; %U20
% 20 4.000 0.0120 37.5510 117.7551 1 1 20 20 2 -1 0 0 0 ; %U21
% 12 2.400 0.0285 26.0611 24.8882 1 1 0 0 1 -1 0 0 0 ; %U22
% 12 2.400 0.0284 25.9318 24.7605 1 1 0 0 1 -1 0 0 0 ; %U23
% 12 2.400 0.0280 25.8027 24.6382 1 1 0 0 1 -1 0 0 0 ; %U24
% 12 2.400 0.0265 25.6753 24.4110 1 1 0 0 1 -1 0 0 0 ; %U25
% 12 2.400 0.0253 25.5472 24.3891 1 1 0 0 1 -1 0 0 0 ]; %U26
% %
% % % %%%%%%%%%%%%%%%DEMAND DATA FOR 26 UNIT SYSTEM%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %
% % % %Demand & Spinning reserve requirement per hour
% P_D=[2223;2052;1938;1881;1824;1825.5;1881;1995;2280;2508;2565;2593.5;2565;2508;2479.5;2479.5;2593.5;2850;2821.5;2764.5;2679;2662;2479.5;2308.5];
%
% SRREQ=0.05*P_D;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%INPUT DATA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=size(I,1); %%%Number of TGU's
T=24; %%%Study period duration
%%%%fmincon.m solver options for F_LIM_ED
OPTS = optimoptions('fmincon','Algorithm','sqp','display','off','ConstraintTolerance',1e-2,'OptimalityTolerance',1.0000e-02,'MaxIterations',100, 'StepTolerance',1e-2,'MaxFunctionEvaluations',100);
%%%%%%%%%%%%%%%%%%%%%%%%%%UNIT DATA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Generator limits
PGI_MAX=I(:,1); %%%Upper generation limit for each generator
PGI_MIN=I(:,2); %%%Lower generation limit for each generator
%%%I/O curve for each generator is modeled by a smooth quadratic function
%%%i.e C(Pgi)=ai*Pgi^2+bi*Pgi+ci
ai=I(:,3);
bi=I(:,4);
ci=I(:,5);
[I_C_SORT_EXS,IDX_EXS,I_C_SORT_INS,IDX_INS] = AFLC(ai,bi,ci,PGI_MAX); %%%Unit priority list assuming same fuel cost for every unit
%%MINIMUM TIME DATA
MUT=I(:,6); %%%Minimum up-time for each generator
MDT=I(:,7); %%%Minimum down-time for each generator
%%START-UP AND SHUT-DOWN COSTS
SU_H=I(:,8);
SU_C=I(:,9);
CSH=I(:,10);
init_status=I(:,11);
%%%%OTHER CONSTRAINTS
M_R=I(:,13); %%MR units
U_N=I(:,14); %%Unavailable units
%%%% MR & UN hours & units
%%Hr begin Hr end
MR=[ 10 15 ];
UN=[ ];
%%%% Prohibitive operative zones (POZ)
%%%Unit %%%POZ
POZ=[ [1 150 165; 1 448,453]
[2 90 110; 2 240 250]
[8 20 30; 8 40 45]
[10 12 17; 10 35 45]];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%MAIN LOOP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%Initiaize solution (NxT) schedule solution as an empty array
Init_SOL=zeros(N,T);
complete=ones(1,T); %%%Nth unit operation schedule is ON for all t
sum_Pgi_max=sum(PGI_MAX);
%%%%%Define swarm size (Either particle or bee)
SWARM_SZ=20;
PART_BEE=0;
YES_CNT=0;
NO_CNT=0;
population=zeros(SWARM_SZ,N*T);
total_COST=zeros(SWARM_SZ,1);
while PART_BEE<SWARM_SZ
PART_BEE=PART_BEE+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Repairing initial schedule according to fixed constraints%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Determine the Must-Run and Unavailable units (Rules B&C)
for j=1:N
if M_R(j)==0 && U_N(j)==0
continue
elseif M_R(j)==1
Init_SOL(j,MR(1):MR(end))=1;
continue
elseif U_N(j)==1
Init_SOL(j,UN(1):UN(end))=0;
continue
end
end
MR_idx=find(M_R>0); %%%Find the exception indices for must run units
UN_idx=find(U_N>0); %%%Find the exception indices for unavailable units
%%%%%%%%%%%%Determine initial status of units %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tdiff_UP=zeros(N,3);%%%Exception indices for coupling constraints Tdiff_UP
Tdiff_DW=zeros(N,3);%%%Exception indices for coupling constraints Tdiff_DW
for k=1:N
Tdiff_UP(k,1)=k;
Tdiff_DW(k,1)=k;
if ismember(k,MR_idx) || ismember(k,UN_idx) %%%If must run or unavailable unit, continue schedule repair without modifying INIT_SOL
continue
else
P=1+(PGI_MAX(k)/sum_Pgi_max);
if init_status(k)>0 %%Unit Online for ti=init_status hours
Tdiff_U=init_status(k)-MUT(k);
if Tdiff_U<0
Tdiff_UP(k,2)=abs(Tdiff_U);
Tdiff_UP(k,3)=1;
Init_SOL(k,1:Tdiff_UP(k,2))=1;
else
Tdiff_UP(k,2)=0;
Tdiff_UP(k,3)=0;
end
for j=Tdiff_UP(k,2)+1:T
X=-1+(P-(-1))*rand;
if X>=0
Init_SOL(k,j)=1;
else
Init_SOL(k,j)=0;
end
end
else %%Unit is offline for ti=init_status hours
Tdiff_D=MDT(k)+init_status(k);
if Tdiff_D>0
Tdiff_DW(k,2)=abs(Tdiff_D);
Init_SOL(k,1:Tdiff_D)=0;
Tdiff_DW(k,3)=1;
else
Tdiff_DW(k,2)=0;
Tdiff_DW(k,3)=0;
end
for j=Tdiff_DW(k,2)+1:T
X=-1+(P-(-1))*rand;
if X>=0
Init_SOL(k,j)=1;
else
Init_SOL(k,j)=0;
end
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%Check & Repair Coupling constraints %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Decompose commitment schedule into T single-hour unit combinations
Exs_On_cap=zeros(3,T); %%Excess online capacity
Ins_On_cap=zeros(3,T); %%Insufficient online capacity
f_unit_ins_all=zeros(N,T); %%%This unit is not part of the exceptions MR,UN,Initialization
free_unit_ins=zeros(N,T); %%%These units are offline and available to be turned off
%%Check Excess online capacity first for full commitment
check_full_exs=ones(1,N)*PGI_MIN; %%%Check if system has excess online capacity for a complete commitment
check_1=find(check_full_exs>P_D);
if isempty(check_1)
clear excp_exs f_unit_exs_all free_unit_exs Exs_On_cap
else
%%%Excess online capacity schedule repair
for t=1:T
u_t=Init_SOL(:,t); %%Decompose into T single hour unit combinations
Exs_On_cap(1,t)=u_t'*PGI_MIN;
Exs_On_cap(2,t)=P_D(t);
if Exs_On_cap(1,t)>P_D(t)
Exs_On_cap(3,t)=1;
format_exs='Excess capacity at hour %u. Repairing schedule...\n';
fprintf(format_ins,t);
end
end
end
status_e=0;
status_f=0;
status_g=0;
status=[status_e,status_f,status_g];
trials=0;
max_trials=15;
while ~all(status)
trials=trials+1;
if trials>max_trials
break
end
if trials==1
%%%Insufficient online capacity initial schedule repair
for t=1:T
u_t=Init_SOL(:,t); %%Decompose into T single hour unit combinations
Ins_On_cap(1,t)=u_t'*PGI_MAX;
Ins_On_cap(2,t)=P_D(t)+SRREQ(t);
if Ins_On_cap(1,t)<P_D(t)+SRREQ(t)
Ins_On_cap(3,t)=1;
for k=1:N
%%%Check Exceptions
if t<=Tdiff_UP(k,2)|| t<=Tdiff_DW(k,2) %%%If the unit has been initialized (MUT/MDT) continue with next unit
continue
else
%%%Obtain available offline generators
f_unit_ins_all(k,t)=k; %%Available units (Both offline and online)
if u_t(k)==0
free_unit_ins(k,t)=k; %%%Available offline units
end
end
end
%%%Repair initial schedule such that online capacity is sufficient
free_unit_ins_tmp=find(free_unit_ins(:,t));
[bool, idx]=ismember(free_unit_ins_tmp,IDX_INS);
idx=sort(idx);
s_idx_off=size(idx,1);
while Ins_On_cap(3,t)==1
for i=1:s_idx_off
u_t(IDX_INS(idx(i)))=1;
Ins_On_cap(1,t)=u_t'*PGI_MAX;
if Ins_On_cap(1,t)<P_D(t)+SRREQ(t)
continue
else
Ins_On_cap(3,t)=0;
Init_SOL(:,t)=u_t;
break
end
end
if Ins_On_cap(3,t)==1
format_ins='Insufficent capacity at hour %u. Check input data\n';
fprintf(format_ins,t);
end
break
end
end
end
[check_F_G,~,~]=check_MUT_MDT(N,T,Init_SOL,init_status,MUT,MDT);
else
prev_SOL=Init_SOL;
v=find(check_E(1,:));
sz_v=length(v);
for t=1:sz_v
v_t=v(t);
for k=1:N
if check_E(1,v_t)==0
break
end
commit_first=IDX_INS(k); %%%Commit the cheapeast generator
u_t=Init_SOL(commit_first,:);
if u_t==complete
continue
else
if u_t(v_t)==0
u_t(v_t)=1;
Init_SOL(commit_first,:)=u_t;
check_E=check_SR_PD(T,P_D,SRREQ,Init_SOL,PGI_MAX);
end
end
end
end
[check_F_G,~,~]=check_MUT_MDT(N,T,Init_SOL,init_status,MUT,MDT);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%Check & Repair Time constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if trials==1
[Init_SOL] = sf_repair_MDT_MUT(N,T,MUT,MDT,Init_SOL,init_status,check_F_G);
[~,status_f,status_g]=check_MUT_MDT(N,T,Init_SOL,init_status,MUT,MDT);
[check_E,status_e]=check_SR_PD(T,P_D,SRREQ,Init_SOL,PGI_MAX);
else %%%Pivot heuristic for MUT/MDT correction
if prev_SOL==Init_SOL
%%%%Repair directly since no pivots where found
vio_u=zeros(N,1);
for j=1:N
if ~all(check_F_G(j,:)==0)
vio_u(j)=j;
end
end
vio_u=vio_u((find(vio_u)));
else
pivot_idx=abs(prev_SOL-Init_SOL);
vio_u=zeros(N,1);
for j=1:N
if ~all(check_F_G(j,:)==0)
vio_u(j)=j;
end
end
vio_u=vio_u((find(vio_u)));
sz_vio_u=length(vio_u);
pvt_u=zeros(sz_vio_u,1);
for j=1:sz_vio_u
k=vio_u(j);
if all(pivot_idx(k,:)==0)
%%%%Repair directly schedule strings with no pivots
x=1;
else
pvt_u(j)=k;
end
end
pvt_u=pvt_u((find(pvt_u)));
sz_pvt_u=length(pvt_u);
%%%Initialize pivot structure
p_s=struct('unit',zeros(sz_pvt_u,1),'pivot_no',zeros(sz_pvt_u,1),'begin_end',zeros(sz_pvt_u,6));
p_s_type=cell(sz_pvt_u,3);
%%%Count number of pivots and define type(Store in a structure)
for j=1:sz_pvt_u
k=pvt_u(j);
p_s.unit(j)=k;
pivot=find(pivot_idx(k,:));
sz_pivot=length(pivot);
diff_pivot=diff(pivot);
if sz_pivot>1
if all(diff_pivot==1)
no_pivot=1;
p_s.begin_end(j,1:2)=[pivot(1),pivot(end)];
p_s_type{j,1}='cluster';
else
cluster=pivot(end)-pivot(1)+1;
if cluster<=MUT(k)
no_pivot=1;
p_s.begin_end(j,1:2)=[pivot(1),pivot(end)];
p_s_type{j,1}='cluster';
else %%%Multiple pivots in a string
if sz_pivot==2
no_pivot=2;
p_s.begin_end(j,1:4)=[pivot(1),pivot(1),pivot(end),pivot(end)];
p_s_type{j,1}='single';
p_s_type{j,2}='single';
else
no_pivot=0;
while ~isempty(pivot)
no_pivot=no_pivot+1;
add=pivot(1)+MUT(k);
sub_pivot=pivot(pivot<=add);
p_s.begin_end(j,2*no_pivot-1:2*no_pivot)=[sub_pivot(1),sub_pivot(end)];
if sub_pivot(1)==sub_pivot(end)
p_s_type{j,no_pivot}='single';
else
p_s_type{j,no_pivot}='cluster';
end
pivot(sub_pivot>0)=0;
pivot=pivot(pivot>0);
end
end
end
end
else
no_pivot=1;
p_s.begin_end(j,1)=pivot(1);
p_s_type{j,1}='single';
end
p_s.pivot_no(j)=no_pivot;
no_pivot=0;
end
%%%%%Begin repair strategy for MUT/MDT violated hours
p_s_complete.pvt=p_s;
p_s_complete.typ=p_s_type;
search_l_r=zeros(sz_pvt_u,1);
for k=1:sz_pvt_u
u=p_s_complete.pvt.unit(k);
n_sch=Init_SOL(u,:);
pivot_no=p_s_complete.pvt.pivot_no(k);
%%%%Single pivot
if pivot_no==1;
t=pivot_no;
typ=p_s_complete.typ{k,t};
switch typ
case 'single'
idx_begin=p_s_complete.pvt.begin_end(k);
idx_end=idx_begin;
search=abs(MUT(u)-1);
case 'cluster'
idx_begin=p_s_complete.pvt.begin_end(k,2*t-1);
idx_end=p_s_complete.pvt.begin_end(k,2*t);
D=abs(idx_end-idx_begin+1);
if ~all(n_sch(idx_begin:idx_end)==1)
n_sch(idx_begin:idx_end)=1;
end
if D~=MUT(u) %%%%Checks left and right to make decision
search=abs(MUT(u)-D);
else
search=3;
end
end
search_l_r(k,1)=search;
search_right=idx_end+search;
search_left= idx_begin-search;
if search_right<=T && search_left>=1
idx_begin_mone=idx_begin-1;
idx_end_pone=idx_end+1;
all_right_on=all(n_sch(idx_end_pone:search_right)==1);
all_left_on=all(n_sch(idx_begin_mone:-1:search_left)==1);
if all_right_on==1 && all_left_on==1
dec='r';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,0,dec,idx_end,T,1);
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
elseif all_right_on==1 && all_left_on==0
l=randi([0,1]);
if l==1
z=until_zero(n_sch,idx_begin,0);
n_sch(idx_begin_mone:-1:z)=1;
else
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
end
dec='r';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,0,dec,idx_end,T,1);
elseif all_right_on==0 && all_left_on==1
r=randi([0,1]);
if r==1
z=until_zero(n_sch,idx_end,1);
n_sch(idx_end_pone:1:z)=1;
else
dec='r';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,0,dec,idx_end,T,1);
end
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
else
l_r=randi([0,1]);
if l_r==1
z=until_zero(n_sch,idx_end,1);
n_sch(idx_end_pone:1:z)=1;
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
else
z=until_zero(n_sch,idx_begin,0);
n_sch(idx_begin_mone:-1:z)=1;
dec='r';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,0,dec,idx_end,T,1);
end
end
else
if search_right>T
n_sch(idx_end:1:T)=1;
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
elseif search_left<1
if init_status(u)<0 && idx_begin<MUT(u)
n_sch(1:MUT(u))=1;
elseif init_status(u)>0 && idx_begin<MDT(u)
n_sch(1:MUT(u))=1;
end
dec='r';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,0,dec,idx_end,T,1);
end
end
%%%%%%Multiple pivots
else
last=pivot_no;
n_sch=Init_SOL(u,:);
decomp=zeros(pivot_no,T);
for j=1:last %%%Decompose n_sch into no_pivot strings
typ=p_s_complete.typ{k,j};
if j<last
else
idx_begin= p_s_complete.pvt.begin_end(k,2*j-1);
idx_end= p_s_complete.pvt.begin_end(k,2*j);
decomp(j,idx_begin:T)=n_sch(idx_begin:T);
if ~all(n_sch(idx_begin:idx_end)==1)
n_sch(idx_begin:idx_end)=1;
end
if idx_end<T
idx_begin_mone=idx_begin-1;
l_r=randi([0,1]);
if l_r==1
z=until_zero(n_sch,idx_begin,0);
n_sch(idx_begin_mone:-1:z)=1;
n_sch(idx_end+1:T)=0;
else
n_sch(idx_end+1:T)=1;
end
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
else
dec='l';
n_sch=mod_repair_MDT_MUT(u,MUT,MDT,n_sch,init_status,dec,idx_begin,0,1);
end
end
end
end
Init_SOL(u,:)=n_sch;
end
[~,status_f,status_g]=check_MUT_MDT(N,T,Init_SOL,init_status,MUT,MDT);
[check_E,status_e]=check_SR_PD(T,P_D,SRREQ,Init_SOL,PGI_MAX);
end
end
status=[status_e,status_f,status_g];
end
if ~all(status==1)
PART_BEE=PART_BEE-1;
continue
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%De-committ excessive units%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%De-committment process
for t=1:T
SRR_PL=check_E(3,t);
Rxs=check_E(2,t);
uj_t=Init_SOL(:,t);
u_on=find(uj_t);
sz_u_on=numel(u_on);
still_xs=true;
while still_xs
Rxs=0;
for k=sz_u_on:-1:1
uj_t(u_on(k))=0;
Rxs=uj_t'*PGI_MAX;
if Rxs>SRR_PL
continue
else
uj_t(u_on(k))=1; %%%When uj_t'*Pgi_max<P_L+SRR re-commit last unit
still_xs=false;
break
end
end
end
Init_SOL(:,t)=uj_t;
end
[check_F_G]=check_MUT_MDT(N,T,Init_SOL,init_status,MUT,MDT);
%%%%%%%%%Repair the schedule in case of any further violations
temp_SOL=Init_SOL;
trial=0;
max_trial=10;
while ~all(all(check_F_G==0))
if trial>=max_trial
break
end
trial=trial+1;
for k=1:N
test_F_G=all(check_F_G(k,:)==0);
if test_F_G==true;
continue
else
if trial==1
n_sch=Init_SOL(k,:);
else
n_sch=temp_SOL(k,:);
end
[intr_on,~] = count_intervals(n_sch);
%%%%MUT violation correction (Correct first)
if intr_on==0;
temp_SOL(k,:)=n_sch;
continue
else
int_vio_on=find(intr_on(:,3)<MUT(k));
sz_ivon=numel(int_vio_on);
for j=1:sz_ivon
u=int_vio_on(j);
idx_begin_on=intr_on(u,1);
idx_end_on=intr_on(u,2);
crr_on=MUT(k)-intr_on(u,3);
dec=randi([0,1]);
if dec==1
if idx_end_on+crr_on>=T
n_sch(idx_end_on:T)=1;
else
n_sch(idx_end_on:idx_end_on+crr_on)=1;
end
else
if idx_begin_on-crr_on<=1
n_sch(idx_begin_on:-1:1)=1;
else
n_sch(idx_begin_on:-1:idx_begin_on-crr_on)=1;
end
end
end
end
[intr_on,intr_off] = count_intervals(n_sch);
%%%%MDT violation (Decommitt units only if possible)
if intr_off==0;
temp_SOL(k,:)=n_sch;
continue
else
int_vio_off=find(intr_off(:,3)<MDT(k));
if isempty(int_vio_off)
temp_SOL(k,:)=n_sch;
continue
else
sz_ivoff=numel(int_vio_off);
for j=1:sz_ivoff
u=int_vio_off(j);
idx_begin_off=intr_off(u,1);
idx_end_off=intr_off(u,2);
if init_status(k)<1 && idx_begin_off==1
continue
elseif idx_end_off==24
continue
else
if idx_end_off==idx_begin_off
n_sch(idx_begin_off)=1;
else
n_sch(idx_begin_off:idx_end_off)=1;
end
end
end
end
end
temp_SOL(k,:)=n_sch;
end
end
[check_F_G,status_f,status_g]=check_MUT_MDT(N,T,temp_SOL,init_status,MUT,MDT);
end
Init_SOL=temp_SOL;
[check_E,status_e]=check_SR_PD(T,P_D,SRREQ,Init_SOL,PGI_MAX);
status=[status_e,status_f,status_g];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%Economic-Dispatch%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Ramping Rate constraints not included
[P_SOL,P_srv,P_COST,tot_gen_COST,itt] = F_LIM_ED(Init_SOL,T,I,P_D,OPTS);
%%%Ramping Rate constraints included
% [P_SOL,P_srv,P_COST,tot_gen_COST,itt,Init_SOL] = F_LIM_ED_RR(Init_SOL,T,I,P_D,IDX_INS,OPTS,POZ);
dff=abs(P_srv-P_D');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%Start-up-cost computation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[su_COST,tot_su_COST] = SU_COST(N,Init_SOL,init_status,SU_H,SU_C,CSH,MUT,MDT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%Display results%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if all((dff)<1e-2)
YES_CNT=YES_CNT+1;
else
NO_CNT=NO_CNT+1;
disp(PART_BEE);
end
population(PART_BEE,:)=reshape(Init_SOL,1,N*T);
total_COST(PART_BEE)=tot_gen_COST+tot_su_COST;
end
disp('%%%%%%%%%%%%%%%%%%%%%SWARM GENERATOR MODULE%%%%%%%%%%%%%%%%%%%%%%%%');
disp('%%%%%%%%%%%%%%%%Coded by Vick @2018 NSYSU EMS-LAB%%%%%%%%%%%%%%%%%%');
fprintf('\n');
format_swarm='%u initial particles created for swarm.\n';
fprintf(format_swarm,SWARM_SZ);
fprintf('\n');
format_particle='%u particles succesfully repaired satisfying all constraints.\n';
fprintf(format_particle,YES_CNT);
fprintf('\n');
disp('Total cost for each particle (Generation+Startup):')
fprintf('\n');
disp(total_COST);
fprintf('\n');
toc
🔗 参考文献
🎈 部分理论引用网络文献，若有侵权联系博主删除
🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：
🌟 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位、冷链、时间窗、多车场等、选址优化、港口岸桥调度优化、交通阻抗、重分配、停机位分配、机场航班调度、通信上传下载分配优化
🌟 机器学习和深度学习时序、回归、分类、聚类和降维
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类
2.10 DBN深度置信网络时序、回归预测和分类
2.11 FNN模糊神经网络时序、回归预测
2.12 RF随机森林时序、回归预测和分类
2.13 BLS宽度学习时序、回归预测和分类
2.14 PNN脉冲神经网络分类
2.15 模糊小波神经网络预测和分类
2.16 时序、回归预测和分类
2.17 时序、回归预测预测和分类
2.18 XGBOOST集成学习时序、回归预测预测和分类
2.19 Transform各类组合时序、回归预测预测和分类
方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断
🌟图像处理方面
图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知
🌟 路径规划方面
旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻、公交车时间调度、水库调度优化、多式联运优化
🌟 无人机应用方面
无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划、
🌟 通信方面
传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化、水声通信、通信上传下载分配
🌟 信号处理方面
信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化、心电信号、DOA估计、编码译码、变分模态分解、管道泄漏、滤波器、数字信号处理+传输+分析+去噪、数字信号调制、误码率、信号估计、DTMF、信号检测
🌟电力系统方面
微电网优化、无功优化、配电网重构、储能配置、有序充电、MPPT优化、家庭用电、电/冷/热负荷预测、电力设备故障诊断、电池管理系统（BMS）SOC/SOH估算（粒子滤波/卡尔曼滤波）、 多目标优化在电力系统调度中的应用、光伏MPPT控制算法改进（扰动观察法/电导增量法）、电动汽车充放电优化、微电网日前日内优化、储能优化、家庭用电优化、供应链优化\智能电网分布式能源经济优化调度，虚拟电厂，能源消纳，风光出力，控制策略，多目标优化，博弈能源调度，鲁棒优化
电力系统核心问题经济调度：机组组合、最优潮流、安全约束优化。新能源消纳：风光储协同规划、弃风弃光率量化、爬坡速率约束建模多能耦合系统：电-气-热联合调度、P2G与储能容量配置新型电力系统关键技术灵活性资源：虚拟电厂、需求响应、V2G车网互动、分布式储能优化稳定与控制：惯量支撑策略、低频振荡抑制、黑启动预案设计低碳转型：碳捕集电厂建模、绿氢制备经济性分析、LCOE度电成本核算风光出力预测：LSTM/Transformer时序预测、预测误差场景生成（GAN/蒙特卡洛）不确定性优化：鲁棒优化、随机规划、机会约束建模能源流分析、PSASP复杂电网建模，经济调度，算法优化改进，模型优化，潮流分析，鲁棒优化，创新点，文献复现微电网配电网规划，运行调度，综合能源，混合储能容量配置，平抑风电波动，多目标优化，静态交通流量分配，阶梯碳交易，分段线性化，光伏混合储能VSG并网运行，构网型变流器， 虚拟同步机等包括混合储能HESS：蓄电池+超级电容器，电压补偿,削峰填谷，一次调频，功率指令跟随，光伏储能参与一次调频，功率平抑，直流母线电压控制；MPPT最大功率跟踪控制，构网型储能，光伏，微电网调度优化，新能源，虚拟同同步机，VSG并网，小信号模型
🌟 元胞自动机方面
交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀
🌟 雷达方面
卡尔曼滤波跟踪、航迹关联、航迹融合、SOC估计、阵列优化、NLOS识别
🌟 车间调度
零等待流水车间调度问题NWFSP 、 置换流水车间调度问题PFSP、 混合流水车间调度问题HFSP 、零空闲流水车间调度问题NIFSP、分布式置换流水车间调度问题 DPFSP、阻塞流水车间调度问题BFSP
👇
5 往期回顾扫扫下方二维码
更多推荐
21
0- 0
已为社区贡献8条内容
所有评论(0)