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

05-10 1871阅读

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为微网运行经济成本，主要包括：与主网交换功率成本，光伏、风机、蓄电池的维护成本，燃气轮机运行成本，污染气体治理成本，燃气轮机开机成本，可转移负荷补偿成本。

考虑用户因参与需求响应改变用电习惯所带来用电感受的变化，另第二个目标函数值为用户不满意度，优化调度的目标是在保证微网经济性的同时也使用户保留较好的用电感受。用户用电不满意度从用电舒适度和用电经济度来考量，式中 $基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）$ 为舒适度， $基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）$ 为经济度， $基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）$ 为t时刻优化前后负荷变化量， $基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）$ 为t时刻的分时电价， $基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）$ 和 $基于免疫粒子群算法的考虑负荷需求相应的热电联供系统优化调度（MATLAB实现）$ 分别为优化前后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= mdt/stp )%MT关闭时间大于下降时间约束
       %       statusMT(i) = 1;%标志MT打开
       %   end
        %  if( statusMT(i) == 0 )
          %    ton=0;%开启时间记为0
         %     toff=toff+1;%关闭时间增加
          %    if(i0
			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