1 主要内容

该程序复现《A Modified Binary PSO to solve the Thermal Unit Commitment Problem》第五章内容，主要做的是一个考虑需求响应的机组组合问题，首先构建了机组组合问题的基本模型，在此基础上，进一步考虑负荷侧管理，也就是需求响应，在调控过程中通过补偿引导负荷侧积极进行需求响应，在模型的求解上，采用了一种基于改进二进制粒子群算法的求解方法，相较于传统的粒子群算法，更加创新，而且求解的效果更好，代码出图效果非常好。该程序函数比较多，主函数为Swarm_generator，运行结果已经保存在Graphs文件夹内部，可以通过运行Graphs.m直接得到出图结果。程序采用matlab编程，注释为英文，适合具有编程经验的同学下载学习！

2 部分代码

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_BEE0); %%%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)