基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）

1.研究背景

        随着“双碳”战略的提出，各种分布式能源的开发和利用收到越来越多的重视。冷热电联供(Combined Cooling Heating and Power, CCHP)系统在发电的同时可以将燃气轮机产生的废热用于制热或制冷，实现能量的梯级利用，并减少系统的污染气体排放量，具有良好的社会和经济效益。
        需求侧管理(Demand response, DR)可以改变微电网的负荷特性。对用户侧来说，不仅提高了用户需求侧灵活响应能力，且节约了用户需求侧的电能；对供给侧来说，需求侧管理则可以提高能源的利用率和可再生能源消纳水平，还提高了微电网的经济性、运行可靠性。
        单纯的追求微电网运行的经济性，不仅会影响用户的用电成本，还会改变用户用电的舒适度。因此对于微电网的优化调度不能只考虑经济最优，还应将用户用电满意度纳入寻优目标。

2.系统结构

地源热泵作为一种高效的可再生能源利用设备，近几年来备受关注。因此，将地源热泵与CCHP结合，建立一种包含地源热泵、风机、光伏、微型燃气轮机、吸收式制冷机组、蓄电装置和蓄冷/热装置的冷热电联供微电网系统，同时，该系统与大电网相连，可以在电价较低时从电网购电，电价较高时向电网售电。该系统所需的天然气由燃气公司提供。结构图如下所示：

2.1 吸收式制冷机模型

制冷机一般选用溴化锂吸收式制冷机，将微型燃气轮机排出的废热制冷。参考文献[1]其模型表示为：

式中，P_MT 为微型燃气轮机输出电功率大小；η_MT为微型燃气轮机的发电效率；Q_MT (t)、Q_AM (t)分别为t时刻微燃机的废热排放量和吸收式制冷机组的制冷功率。

2.2 地源热泵机组模型

参考文献[1]，地源热泵是一种利用浅层地热资源，既能制热又能制冷的高效节能空调技术，比常规中央空调节能40%以上，在大型公共建筑节能中潜力巨大。其模型为：

式中，Q_HP (t)、P_HP (t)分别为t时刻地源热泵产生的冷/热功率和地源热泵消耗的电功率；C_OP为地源热泵的制冷/制热系数，本模型中取C_OP=4。

2.3 储热/冷装置模型：

冷热储能装置采用水蓄冷设备，与电储能装置具有相似的运行特征，其模型可表示为：

式中，S_ES (t)、S_ES (t-1)分别为储能装置t和t-1时刻的容量；𝝉为储能损失系数；P_(ES.ch) (t)、P_(ES.dis) (t)分别为t时刻储能装置充放能功率大小； η_ch 、η_dis分别为充放能效率。

3.优化模型

3.1 目标函数

        目标函数1为微网运行经济成本，主要包括：与主网交换功率成本，光伏、风机、蓄电池的维护成本，燃气轮机运行成本，污染气体治理成本，燃气轮机开机成本，可转移负荷补偿成本。
        考虑用户因参与需求响应改变用电习惯所带来用电感受的变化，另第二个目标函数值为用户不满意度，优化调度的目标是在保证微网经济性的同时也使用户保留较好的用电感受。用户用电不满意度从用电舒适度和用电经济度来考量，式中 $F_{com}$ 为舒适度， $F_{eco}$ 为经济度， $\Delta P_{L}(t)$ 为t时刻优化前后负荷变化量， $\lambda _{t}$ 为t时刻的分时电价， $P_{La}(t)$ 和 $P_{Lb}(t)$ 分别为优化前后t时刻的负荷量，α为用户用电舒适度占用户满意度的权重系数，β为用户用电经济度占用户满意度的权重系数。
        可以看出，经济度表示负荷转移后电价变化对于用户用电感受的影响，而舒适度表示负荷转移后因用电习惯改变而带来的用电感受的影响。本文将用户定义为三类：第一类为用电舒适度占用电满意度权重较大的用户，代表生产时间以及生产流程较为严格的工厂车间等用电负荷；第二类代表用户对于用电舒适度与用电满意度的权重相同，适用于住宅居民区等用电负荷；第三类代表电费占生产成本较大且对用电时间没有严格要求的商业化用户，具体的满意度权重取值如下表所示：

3.2 约束条件

4.求解算法

本文采用基于疫苗接种的免疫粒子群算法对所建立的优化模型进行求解。

4.1 粒子群（PSO）算法

用一种粒子来表示一个个体，每个粒子可视为N维搜索空间中的一个搜索个体，粒子的当前位置即为对应优化问题的一个候选解，粒子的飞行过程即为该个体的搜索过程．粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整．粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值，粒子群中最优的个体极值作为当前全局最优解。不断迭代，更新速度和位置。最终得到满足终止条件的最优解。

算法流程如下：

① 初始化

首先，我们设置最大迭代次数，目标函数的自变量个数，粒子的最大速度，位置信息为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模，每个粒子随机初始化一个速度。

② 个体极值与全局最优解

定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。

③ 更新粒子的速度和位置

④ 终止条件

（1）达到设定迭代次数；（2）代数之间的差值满足最小界限。

4.2 基于疫苗接种的免疫粒子群算法

首先将传统固定惯性权重系数更改为根据全局最优点自适应调整惯性权重系数，即权重根据粒子的位置不同而动态变化。采用的非线性动态惯性权重系数公式为：

此外，引入基于疫苗接种的免疫优化来提升算法的综合性能，面对待求解问题时，相对于面对各种抗原，可以提前注射“疫苗”来抑制退化问题，从而更加保持优胜劣汰的特点，使算法一直优化下去，达到免疫的目的。

算法流程图：

5. MATLAB代码

以下给出main函数代码，具体完整版代码请前往博主资源里下载。

clear all;
close all;
clc;
tic
%% INPUTS
pso_flag = 1;
swarm_size = 100;
BATTpmax = 200000;%电池出力上限200kw
BATTprated=400000;%电池容量400kWh
BATTpmin = -150000;%电池出力下限-150KW
Chureqmax=100000;%储热罐出力上限100kW
Chureqrated=300000;%储热罐基准容量300kW
Chureqmin=-150000;%储热罐出力下限-150kW
PVpmax =100000;%光伏出力上限100KW
WTpmax =80000;%风机出力上限80KW
MTpmax = 50000;%燃气轮机出力上限50KW
MTpmin=5000;%不停机运行的出力下限5kw
Putilmax=200000;%与主网交换功率最大值200kW
    %% INIALIZE TIME STEPS AND AC BUS（初始化时间步长以及交流母线）
    stp = 3600;
    pertime = 24;
    totalTime = pertime*1; % total hours simulation
    totalUnits = 5;%总共五个单元
    p_bus = zeros(totalUnits,totalTime);%规格：（5，24）
    index_GenMT = 1;
    index_GenPV = 2;
    index_GenWT = 3;
    index_GenBTR = 4;
    index_Load = 5;
    index_Util = 6;
    zeros24(1:pertime) = 0.5;
    %% MICROURBINE PARAMETERS FOR EMS （微型燃气轮机参数）
    %pMTmax = 20*10^3;
    cmt = 0.3571; % ($/m^3)
    fmt = 0.0085; % (m^3/Wh)
    dt = 1; % time stp, the optimazation is updated every time stp (h)仿真步长为1h
    koc = 0.0587; % ($/Wh)Wh:瓦特小时
    hot_startup = 30; %sec
    cold_startup = 200; %sec
    cooling_time = 520; %sec
    mut = 600;% minimum up time (sec)
    mdt = 300;% minimum down time (sec)
    ton=0;% counting the times turned on
	tonmin=2;
    toff=0;% counting the times turned off
	toffmin=2;
    statusMT=ones(1,pertime);%先假设MT一直开机
	Rup=75000;%爬坡上升速率限制 75kWh
	Rdown=75000;%爬坡下降速率限制 75kWh
	
        %% WIND TURBINE PARAMETERS FOR EMS（风机参数）
    Cwt = 0.0001;% $0.1/Wh
    statusWT=ones(1,totalTime);
        %% BATTTERY PARAMETERS FOR EMS（蓄电池参数）
    soc = zeros(1,totalTime);
    soc(1,:) = 0.5; %initialize battery's sate of charge to 50%
    statusBTR=ones(1,totalTime);
	
	 %% 储热罐参数
	socH=zeros(1,totalTime);
	socH(1,:) = 0.5; %initialize sate of charge to 50%
	%% 地源热泵参数
	Cop=4;%地源热泵的制冷制热系数
	
    minchargpwr = 1000;%最小充电功率
    maxchargpwr = 16000;%最大充电功率
    effBTR = 0.9;%电池效率
        %% PV array PARAMETERS FOR EMS（光伏系统参数）
    Cpv = 0.2;% $0.2/kWh
        %% LOAD PARAMETERS FOR EMS（负载侧参数）
    powerDemandmin = 150*10^3;%负载需求最小值
    powerDemandmax = 300*10^3;%负载需求最大值
       %% UTILITY
    cutil = 1.05;
    counterDisp=0;
    util = zeros(1,totalTime);
    stp_time =0;
    irrad = zeros(totalTime);
    tempr = zeros(totalTime);
    vwind = zeros(totalTime);
    p_min=zeros(3,pertime);
    gb_mt=zeros(1,pertime);
    gb_batt=zeros(1,pertime);
    gb_util=zeros(1,pertime);
	irrad=[0 0 0 0 0 0 0 0.2 0.3 0.4 0.6 0.8 0.8 0.8 0.8 0.7 0.6 0.4 0.3 0.1 0 0 0 0];%定义光照强度
	tempr=[10 10 9 9 8 9 10 12 15 18 20 22 25 25 25 25 24 21 18 16 14 12 11 10];%定义温度
	%vwind=[10 10 11 11 10 9 8 8 7 9 10 10 8 7 7 8 9 9 9 10 11 10 9 8];%定义风速，3-13
	vwind=[10,11,12,12,12,12,12,13,15,16,17,17,18,19,12,13,21,22,23,24,25,26,27,28];%每小时的风速
	Pwind=[0,0,0,0,0,0,5000,25000,80000,90000,100000,800000,90000,100000,90000,40000,42000,45000,90000,100000,100000,30000,5000,0];%直接定义风机出力值
	Qload=[95000,96000,97000,94000,95000,92000,90000,88000,88000,80000,75000,73000,65000,45000,60000,72000,80000,82000,85000,90000,91000,90000,94000,94500]*1;%定义冬季典型日热负荷
	%Qload=[60000,58000,62000,64000,70000,68000,70000,74000,79000,82000,83000,100000,120000,125000,130000,94000,92000,90000,88000,90000,80000,62000,60000,58000]*1;%定义夏季典型日的冷负荷
	
	
     for i=1:totalTime
        %% LOAD DEMAND
        stp_time=stp_time+1;

        if(stp_time<=pertime)% stp_time<=24，分时负荷需求赋值
            if (stp_time==1 || stp_time==8 || stp_time==9 || stp_time==16 || stp_time==17 || stp_time==24||stp_time==23||stp_time==22)
                p_bus(index_Load,i) = 50000;
            elseif (stp_time==11 || stp_time==12 || stp_time==15 ||stp_time==19)
                p_bus(index_Load,i) = 120000;
            elseif (stp_time==10||stp_time==18)
                p_bus(index_Load,i) = 100000;
            elseif (stp_time==21)
                p_bus(index_Load,i) = 220000;
            elseif (stp_time==14 || stp_time==20)
                p_bus(index_Load,i) = 180000;
            elseif (stp_time==13)
                p_bus(index_Load,i) = 280000;
            else
                p_bus(index_Load,i) = 10000;
            end
        else %当运行时间大于24h
              if (mod(stp_time,pertime)==1 || mod(stp_time,pertime)==8 || mod(stp_time,pertime)==9 || mod(stp_time,pertime)==16 || mod(stp_time,pertime)==17 || mod(stp_time,pertime)==24)
                p_bus(index_Load,i) = 30000;
            elseif (mod(stp_time,pertime)==11 || mod(stp_time,pertime)==12 || mod(stp_time,pertime)==15 || mod(stp_time,pertime)==18 || mod(stp_time,pertime)==19 || mod(stp_time,pertime)==23)
                p_bus(index_Load,i) = 60000;
            elseif (mod(stp_time,pertime)==10)
                p_bus(index_Load,i) = 90000;
            elseif (mod(stp_time,pertime)==21 || mod(stp_time,pertime)==22)
                p_bus(index_Load,i) = 150000;
            elseif (mod(stp_time,pertime)==14 || mod(stp_time,pertime)==20)
                p_bus(index_Load,i) = 210000;
            elseif (mod(stp_time,pertime)==13)
                p_bus(index_Load,i) = 300000;
            else
                p_bus(index_Load,i) = 3000;
            end
        end
              %% CONSTRAINTS MT（微型燃气轮机的约束）***
       % if( ton >= mut/stp )%600/3600，将单位由s转换为h；MT打开时间大于上升时间约束
        %      statusMT(i) = 0;%标志MT关闭
        %  end

       %   if( toff >= mdt/stp )%MT关闭时间大于下降时间约束
       %       statusMT(i) = 1;%标志MT打开
       %   end

        %  if( statusMT(i) == 0 )
          %    ton=0;%开启时间记为0
         %     toff=toff+1;%关闭时间增加
          %    if(i<length(statusMT))
         %         statusMT(i+1)=0;
         %     end
        %  else
         %     toff = 0;
         %     ton = ton+1;
         %     if(i<length(statusMT))
          %        statusMT(i+1)=1;%当时间小于24h，标记下一个小时MT的状态为打开，顺延上一时刻的MT的状态
          %    end
         % end
                %% POWER GENERATORS (RES) 发出电能
%         irrad(i) = 0.4;
%         tempr(i) = 25 ;
%         vwind(i) = 23;
      %  irrad(i) = rand();%当前光照强度，随机选取0-1
      %  tempr(i) = 25 + (75-25).*rand();%当前温度***
       % vwind(i) = 10 + (28-10).*rand();%当前风速
       
	   p_bus(index_GenWT,i) = windturbine(vwind(i),statusWT(i),WTpmax);%P_bus第三行记录风机各个小时的发电量
	%	p_bus(index_GenWT,i) =Pwind(i);
		
        %if(mod(stp_time,pertime)<=19 && mod(stp_time,pertime)>=7)%一天中的7：00到19：00光伏出力
         %   p_bus(index_GenPV,i) = pv_array(irrad(i),tempr(i),PVpmax);%Pbus第二行记录光伏出力
        %end
         p_bus(index_GenPV,i) = pv_array(irrad(i),tempr(i),PVpmax);%Pbus第二行记录光伏出力
		 
         if(mod(stp_time,pertime)==0)%到达第24h
            if(pso_flag==1)
                soc_tmp=soc;
                %% INITIALIZE POPULATION
                [ DV, putil, soc_ ,status_,socH_] = init_swarms(p_bus(index_GenWT,i-pertime+1:i),p_bus(index_GenPV,i-pertime+1:i),p_bus(index_Load,i-pertime+1:i),soc(i-pertime+1:i),statusMT(i-pertime+1:i), BATTpmax, BATTpmin, MTpmax,MTpmin,Putilmax,swarm_size,Qload,Chureqmax,Chureqmin,socH,BATTprated,Chureqrated);%初始化种群，因为此时i是24的整数倍，故输入参数是p_bus每行的第1到24列。
                %输出DV：蓄电池出力、燃气轮机出力；putil与电网交换功率；蓄电池soc值
                %% PSO
                [final_global_best,ind, gbest, val, swarm, bestval, gb_mt, gb_batt, gb_util, soc,soc_best,Pload_af_sheding,flj,gb_chure,socH_global,socH_best,Php_,QAM_] = pso(DV, putil, p_bus(index_Load,i-pertime+1:i), p_bus(index_GenPV,i-pertime+1:i), p_bus(index_GenWT,i-pertime+1:i), soc_, status_,BATTpmax, BATTpmin, MTpmax,MTpmin,Putilmax,swarm_size,Qload,Chureqmax,Chureqmin,socH_,BATTprated,Chureqrated);
                p_bus(index_GenMT,i-pertime+1:i) = gb_mt(1:pertime);
                p_bus(index_GenBTR,i-pertime+1:i) = gb_batt(1:pertime);
                p_bus(index_Util,i-pertime+1:i) = gb_util(1:pertime);
				Qhp_=Php_*Cop;%最优解对应地源热泵输出热功率

                soc_tmp(i-pertime+1:i) = soc;
                soc = soc_tmp;

                h=figure;
                set(gcf,'Visible','on');
                plot(gbest)
                xlabel('Generations')
                ylabel('Cost($)')
                grid on
                title('最优解对应的目标函数值')

               % baseFileName =  sprintf('%d_%d_%d_%d.jpg', BATTpmax, PVpmax, WTpmax, MTpmax);
               % saveas(h, baseFileName ); 
			 end
		end
	end
	toc
	%% 画图
	t=1:24;
	h2=figure;
	set(gcf,'Visible','on');
	plot(t,p_bus(index_GenMT,1:24),'-o',t,p_bus(index_GenPV,1:24),'-*',t,p_bus(index_GenWT,1:24),'-+',t,p_bus(index_GenBTR,1:24),'-h',t,Pload_af_sheding(:),'-s',t,p_bus(index_Util,1:24),'-d');
	xlabel('Time(Hour of day)')
    ylabel('Power Output (Watt)')                
    grid on;
	title('最优解对应的各单元出力情况')
    legend('MT','PV','WT','BTR','Load','Utility');
    
    h3=figure;
    set(gcf,'Visible','on');
    plot(t,soc,'-*');
    xlabel('Time(Hour of day)');
    ylabel('SOC');
    title('最优解对应的蓄电池SOC值变化情况');
    
	h4=figure;
    set(gcf,'Visible','on');
     C_buy=[0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.8	0.8	0.8	1.1	1.1...
    1.1	1.1	1.1	0.8	0.8	0.8	1.1	1.1	1.1	0.8	0.8	0.8].*10^(-3);%购电价格
	C_sell=[0.35 0.35 0.35 0.35 0.35 0.35 0.35 0.56 0.56 0.56 0.76 0.76 0.76 0.76 0.76 0.56 0.56 0.56 0.76 0.76 0.76 0.56 0.56 0.56].*10^(-3);%售电价格
    plot(t,C_buy,'-o',t,C_sell,'-*');
	legend('购电价格','售电价格');
    xlabel('Time(Hour of day)');
    ylabel('每小时电价（￥/W）');
    title('分时电价');
	
	Psell=zeros(1,24);
	Pbuy=zeros(1,24);
	Pcha=zeros(1,24);
	Pdis=zeros(1,24);
	for	i=1:24
		if	p_bus(index_Util,i)>0
			Pbuy(i)=p_bus(index_Util,i);
		else 
			Psell(i)=p_bus(index_Util,i);
		end
		if	p_bus(index_GenBTR,i)>0
			Pdis(i)=p_bus(index_GenBTR,i);
		else
			Pcha(i)=p_bus(index_GenBTR,i);
		end
	end
	h5=figure;
	PP=[p_bus(index_GenMT,1:24);p_bus(index_GenPV,1:24);p_bus(index_GenWT,1:24);Pdis(1:24);Pbuy(1:24)];
	PP_neg=[Psell(1:24);Pcha(1:24)];
	bar(PP','stack');
	hold on
	bar(PP_neg','stack');
	plot(t,Pload_af_sheding(:),'r','linewidth',2);
	h=legend('燃气轮机出力','光伏出力','风机出力','蓄电池放电','从主网购电','向主网售电','蓄电池充电','预测负荷曲线');
	set(h,'Orientation','horizon')
	xlabel('时段');ylabel('功率/W');
	hold off
	
	h6=figure;
	set(gcf,'Visible','on');
	p_bus(index_Load,1:24)=p_bus(index_Load,1:24)+Php_;
	plot(t,p_bus(index_Load,1:24),'-s',t,Pload_af_sheding(:),'-h');
	xlabel('时间/h');
	ylabel('需求响应前日负荷曲线/W');
	title('日负荷曲线');
	legend('DR前的日负荷曲线','DR后的日负荷曲线');
	
	h7=figure;
	set(gcf,'Visible','on');
	plot(flj.obj1,flj.obj2,'*');
	xlabel('经济成本Cost(￥)');
	ylabel('用户参与DR的不满意度')
	title('多目标寻优的Pareto解集');
	
	h8=figure;
	set(gcf,'Visible','on');
	plot(t,socH_global,'-h');
    xlabel('Time(Hour of day)');
    ylabel('SOCH');
    title('最优解对应的储热装置SOCH值变化情况');
	
	h9=figure;
	set(gcf,'Visible','on');
	plot(t,QAM_,'-o',t,Qhp_,'-*',t,gb_chure,'-+',t,Qload,'-s');
	xlabel('Time(Hour of day)')
    ylabel('Power Output (Watt)')                
    grid on;
	title('最优解对应的各热源单元出力情况')
    legend('吸收式制冷机','地源热泵','储热装置','热负荷');
	
	Qcha=zeros(1,24);
	Qdis=zeros(1,24);
	for	i=1:24
		if	gb_chure(i)>0
			Qdis(i)=gb_chure(i);
		else 
			Qcha(i)=gb_chure(i);
		end
	end
	h10=figure;
	PP2=[QAM_;Qhp_;Qdis];
	PP2_neg=[Qcha];
	bar(PP2','stack');
	hold on
	bar(PP2_neg','stack');
	plot(t,Qload,'r','linewidth',2);
	h=legend('吸收式制冷机出力','地源热泵出力','储热装置热出力','储热装置储热','热负荷曲线');
	set(h,'Orientation','horizon')
	xlabel('时段');ylabel('功率/W');
	hold off
	
	
	
	flj

6.算例仿真结果

6.1 分时电价情况

6.2 目标函数1的收敛情况

6.3 系统电部分各单元出力情况

6.4 蓄电池SOC变化情况

6.5 负荷需求响应情况

6.6 系统热部分各单元出力情况

6.7 储热装置SOCH变化情况

参考文献

[1]李坚,吴亮红,张红强,王维,贾睿.基于排序交叉优化算法的冷热电联供微电网经济调度[J].电力系统保护与控制,2021,49(18):137-145.DOI:10.19783/j.cnki.pspc.201556.
[2]杨欢红,唐芃芃,黄文焘.考虑用户满意度的基于需求侧管理的微电网多目标优化调度[J/OL].电测与仪表:1-8.
[3]毛晓明,陈深,吴杰康,郭壮志.分时电价机制下含蓄电池微网的优化调度[J].电网技术,2015,39(05):1192-1197.DOI:10.13335/j.1000-3673.pst.2015.05.004.
[4]《MATLAB智能算法》