1 主要内容

本程序参考文献《考虑可再生能源消纳的建筑综合能源系统日前经济调度模型》模型，建立了电热综合能源系统优化调度模型，包括燃气轮机、燃气锅炉、余热锅炉、蓄电池、储热、风电和光伏等设备模型，注意：不含文献中的地源热泵和空调设备，可以作为比较经典的电热综合能源优化调度学习资料。

模型示意图

注意：本程序不含地源热泵模型。

目标函数

目标函数:经济性=燃料成本+设备运维成本+环境保护

程序亮点

程序对于蓄电池和储热设备的约束形式较常规方式有所区别，见部分程序59-66行，虽然这种约束方式略显笨重，但提供了另外一种程序编写方式，顺便了解implies命令应用方式。
程序对于蓄电池和储热采用自放电/热系数方式，因此可以合理解释结果图2中储热满足24小时充放电功率一致，但是热能却一直下降的原因。

2 部分程序

%% 程序初始化
clc;
clear;
close all;

%% 定义变量
%电设备
P_PV=sdpvar(1,24,'full');%光伏发电
P_WT=sdpvar(1,24,'full');%风力发电
P_GT=sdpvar(1,24,'full');%燃气轮机电功率出力
Eb=sdpvar(1,24,'full');% 电储能的容量
P_ES=sdpvar(1,24,'full');%电储能输出，有正有负
P_ESc=sdpvar(1,24,'full');% 电储能充电功率
P_ESd=sdpvar(1,24,'full');% 电储能放电功率
%热设备
Q_GT=sdpvar(1,24,'full');%燃气轮机输出的热功率
Q_RB=sdpvar(1,24,'full');%余热锅炉输出的热功率
Q_GB=sdpvar(1,24,'full');%燃气锅炉输出的热功率
Eh=sdpvar(1,24,'full');%热储能的容量
P_HS=sdpvar(1,24,'full');%热储能输出，有正有负
P_HSc=sdpvar(1,24,'full');%热储能储热功率
P_HSd=sdpvar(1,24,'full');%热储能放热功率
%辅助变量
Ubc=binvar(1,24,'full');%电储能充电状态，1表示充电
Ubd=binvar(1,24,'full');%电储能放电状态，1表示放电
Uhc=binvar(1,24,'full');%热储能储热状态，1表示储热
Uhd=binvar(1,24,'full');%热储能放热状态，1表示放热
F_GT=sdpvar(1,24,'full');%燃气轮机消耗的天然气量
F_GB=sdpvar(1,24,'full');%燃气锅炉消耗的天然气量
%% 常量数据,包括电热冷负荷和风光出力
Pwt=[66.9  68.2  71.9  72  78.8  94.8  114.3  145.1  155.5  142.1  115.9  127.1  141.8  145.6   145.3  150  206.9  225.5  236.1  210.8  198.6  177.9  147.2  58.7];
Ppv=[0,0,0,  0,  0.06,  6.54,  20.19,  39.61,  49.64,  88.62,  101.59,  66.78,  110.46,  67.41,  31.53, 50.76,20.6,22.08,2.07,0,0,0,0,0];
load_e=[88.24   83.01   80.15   79.01   76.07   78.39   89.95   128.85   155.45   176.35   193.71   182.57   179.64   166.31   164.61   164.61   174.48   203.93   218.99   238.11   216.14   173.87   131.07   94.04];
load_h=[122 125 128 125 128 119 113 109 100,90 88,79 66,55 63,71 98 109 139 150 142 139 136 134];

%% 目标函数:经济性=燃料成本+设备运维成本+环境保护
z=0;
for t=1:24
    z=z+0.025*P_GT(t)+0.0097*Q_GB(t)+0.025*Q_RB(t)...%设备投资运行费用+0.043*P_WT(t)+0.0096*P_PV(t)
    +0.03*(abs(P_ESc(t))+abs(P_ESd(t)))+0.02*(abs(P_HSc(t))+abs(P_HSd(t)))...%储能设备
    +3.15*(F_GT(t)+F_GB(t))...%买天然气的价格
    +0.088*(P_GT(t)*0.242+Q_GB(t)*0.085)+0.624*(P_GT(t)*0.0036+Q_GB(t)*0.764)+1.95*(P_GT(t)*0.2+Q_GB(t)*0.54);%污染物处理
end
      F=z
      
%% 约束条件

C=[];
%机组约束,决策变量
C=[C,0.8*Ppv<=P_PV<=Ppv];%光伏发电上下限约束
C=[C,0.8*Ppv<=P_WT<=Pwt];%风力发电上下限约束
for t=1:24
    C=[C,0<=P_GT(t)<=550];%燃气轮机上下限约束√
    C=[C,0<=Q_RB(t)<=750];%余热锅炉上下限约束√
    C=[C,0<=Q_GB(t)<=700];%燃气锅炉上下限约束√
    C=[C,0<=Q_GT(t)<=200];%燃气轮机输出的热功率
%设备模型

    C=[C,Ubc(t)+Ubd(t)<=1];
    C=[C,Uhc(t)+Uhd(t)<=1];%确定充放电状态不可能同时发生
    C=[C,-80<=P_ES(1,t)<=80,0<=P_ESc(1,t)<=80,-80<=P_ESd(1,t)<=0]; 
    C=[C,-80<=P_HS(1,t)<=80,0<=P_HSc(1,t)<=80,-80<=P_HSd(1,t)<=0];   
    C=[C,implies(Ubc(1,t),[P_ES(1,t)>=0,P_ESc(1,t)==P_ES(1,t),P_ESd(1,t)==0])];%储电=负荷
    C=[C,implies(Ubd(1,t),[P_ES(1,t)<=0,P_ESd(1,t)==P_ES(1,t),P_ESc(1,t)==0])];
    C=[C,implies(Uhc(1,t),[P_HS(1,t)>=0,P_HSc(1,t)==P_HS(1,t),P_HSd(1,t)==0])];%储热=负荷
    C=[C,implies(Uhd(1,t),[P_HS(1,t)<=0,P_HSd(1,t)==P_HS(1,t),P_HSc(1,t)==0])];