基于改进二进制粒子群算法的含需求响应机组组合问题研究（matlab代码）

目录

1 主要内容

2 部分代码

3 程序结果

 4 下载链接

---

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)

3 程序结果

 4 下载链接

详见文后联系方式-->程序目录