EI论文复现：考虑冷热运行特性的综合能源系统多时间尺度优化调度程序代码！

适用平台/参考文献：Matlab+Yalmip+Cplex；

参考文献：电力系统自动化《含冰蓄冷空调的冷热电联供型微网多时间尺度优化调度》

提出考虑冷热特性的综合能源系统多时间尺度优化调度模型，日前计划中通过多场景描述可再生能源的不确定性，侧重于一个运行优化周期内系统的经济运行；日内调度基于日前计划方案，根据冷热电在不同时间尺度上的相关性和互补性，提出考虑冷热负荷变化的双层滚动优化平抑模型。程序中算例丰富、注释清晰、干货满满，创新性很高！下面对文章和程序做简要介绍！

程序创新点：

日前调度中，考虑新能源处理的不确定性，利用拉丁超立方抽样生成新能源典型出力场景；

日内调度中，建立考虑不同能源时间尺度差异的双层滚动优化调度模型，使冷热能在“时”时间尺度上进行调度，使电能在“分”时间尺度上进行调度，及时纠正预测误差引起的功率波动，实现对系统日前方案的修正，保证综合能源系统的稳定运行。

主要工作：

建立考虑冷热特性的多时间尺度优化调度模型，在日前计划模型，研究冷-热特性的运行方式对综合能源系统综合效益的影响，并考虑可再生能源不确定性，以运行费用最小为优化目标。在日内调度，考虑不同能源的调节时间尺度差异，建立双层滚动优化调度模型，基于滚动时长的差异，使冷热能在“时”时间尺度上进行调度，使电能在“分”时间尺度上进行调度，及时有效纠正预测误差引起的功率波动，实现对系统日前方案的修正，实现系统的稳定运行。

拉丁抽样方法：

拉丁超级方抽样（Latin Hypercube Sampling，简称LHS）是一种用于在多维空间中进行抽样的统计方法。它是在拉丁方抽样的基础上进行改进的一种技术，常用于数值模拟、实验设计和优化问题中。以下是拉丁超级方抽样的基本步骤：

确定抽样维度：首先，确定需要抽样的变量和其相应的分布。这些变量可以是模型中的输入参数，例如在数值模拟中的各种参数。

划分区间：将每个变量的分布范围划分为相等的子区间。这些区间的数量通常是相对较少的，而且是平均分布的。

生成拉丁超级方：创建一个拉丁方，确保在每个维度上都有且仅有一个样本。这可以通过在每个行和列上使用每个值一次的方式来实现。这个步骤确保了在每个维度上的均匀抽样。

填充超级方格：将每个区间的中心与相应的拉丁方中的一个样本关联起来。这样就得到了一个拉丁超级方抽样。

应用于模型：使用生成的抽样点作为模型的输入，执行模拟、实验或优化。这有助于在多维参数空间中更全面地探索。

拉丁超级方抽样的优势在于它能够提供相对均匀的抽样，而且相较于传统的随机抽样，它更有效地覆盖了参数空间。在有限的样本数量下，对模型或系统的全面性进行评估成为可能。

实时调度框架：

程序结果：

（1）上图是文中设备出力柱状图，下图是复现程序中设备出力柱状图。

（2）日内实时调度设备出力折线图

（3）日前调度主要热冷设备出力折线图

（4）风光典型场景出力

（5）日前调度设备出力

部分程序：

%电转热
Peh_max=;n_Peh=;%转换设备/转换率%热储能
H_storage_max=; h_n=;h_charge=;h_discharge=;%热储能容量/自损/充热/放热；%电储能
E_storage_max=; e_n=;e_charge=;e_discharge=;%电储能容量/自损/充电/放电;
%冷储能C_storage_max=;Ncool=;ncool=;kcool=%容量/自损量/充冷系数%PGrid
Pgrid_max=;%PCC交互容量；%%储冰空调+普通空调Qa_max=[];%空调模式Qa_min=];
Qc_max=[];%制冰模式Qc_min=[0;Qd_max=[120];%融冰模式Qd_min=[0];%设置变量
P=intvar(3,24,'full');
UP=binvar(3,24,'full'); %三台机组Pc=intvar(3,24,'full');%分配转冷部分Ph=intvar(3,24,'full');%分配转热部分
buy_Pgrid=intvar(1,24,'full');Ubuy_Pgrid=binvar(1,24,'full');%Pgrid购
Pcharge=intvar(1,24,'full');UPcharge=binvar(1,24,'full');%  充电  
Pdischarge=intvar(1,24,'full');UPdischarge=binvar(1,24,'full');%  放电  
H=intvar(1,24,'full');%锅炉Hh=intvar(1,24,'full');%转热Hc=intvar(1,24,'full');%转冷
Peh=intvar(1,24,'full');%电转热
Hti=intvar(1,24,'full');UHti=binvar(1,24,'full');%充热
Hto=intvar(1,24,'full');UHto=binvar(1,24,'full');%放热
Qa=intvar(1,24,'full');Uqa=binvar(1,24,'full');%
Qc=intvar(1,24,'full');Uqc=binvar(1,24,'full');%
Qd=intvar(1,24,'full');Uqd=binvar(1,24,'full');%
%设置目标函数buy_Price=[0.45 0.73 1.21];%谷/峰 购电价
for i=1:24
    H_h(1,i)=aH(1)*Ph(1,i)+bH(2)*UP(2,i)+Hh(1,i)*0.98;%转热功率
    H_c(1,i)=(aH(1)*Pc(1,i)+bH(2)*UP(2,i)+Hc(1,i))*0.95;%转冷功率
for i=1:24 %蓄冰模式下的耗电量
    Paa(1,i)=Qa(1,i)*0.45+3.6*Uqa(1,i);
    Pcc(1,i)=Qc(1,i)*0.44;
    Pdd(1,i)=8*Uqd(1,i); 
K1=;%热效益K2=;%冷效益aaa=0;
for i=1:24
    if i<=6%%谷
       objective=aaa+Pcharge(1,i)*0.00362+Pdischarge(1,i)*0.00182+Cgas*(aF(1)*P(1,i)+bF(1)*UP(1,i)+aF(2)*P(2,i)+bF(2)*UP(2,i))+aF(3)*P(3,i)+bF(3)*UP(3,i)...
           +Cgas*H(1,i)/H_n+buy_Price(1).*buy_Pgrid(1,i);
       aaa=objective;
       mPaa(1,i)=Paa(1,i)*buy_Price(1);
       mPcc(1,i)=Pcc(1,i)*buy_Price(1);
       mPdd(1,i)=Pdd(1,i)*buy_Price(1);
       mPgrid(1,i)=buy_Price(1).*buy_Pgrid(1,i);
    elseif i<=11&&i>=7%峰
       objective=aaa+Pcharge(1,i)*0.00362+Pdischarge(1,i)*0.00182+Cgas*(aF(1)*P(1,i)+bF(1)*UP(1,i)+aF(2)*P(2,i)+bF(2)*UP(2,i))+aF(3)*P(3,i)+bF(3)*UP(3,i)...
           +Cgas*H(1,i)/H_n+buy_Price(3).*buy_Pgrid(1,i);
       aaa=objective;
       mPaa(1,i)=Paa(1,i)*buy_Price(3);
       mPcc(1,i)=Pcc(1,i)*buy_Price(3);
       mPdd(1,i)=Pdd(1,i)*buy_Price(3);
       mPgrid(1,i)=buy_Price(3).*buy_Pgrid(1,i);
    elseif i>=12&&i<=17%平
       objective=aaa+Pcharge(1,i)*0.00362+Pdischarge(1,i)*0.00182+Cgas*(aF(1)*P(1,i)+bF(1)*UP(1,i)+aF(2)*P(2,i)+bF(2)*UP(2,i))+aF(3)*P(3,i)+bF(3)*UP(3,i)...
           +Cgas*H(1,i)/H_n+buy_Price(2).*buy_Pgrid(1,i);
       aaa=objective;
       mPaa(1,i)=Paa(1,i)*buy_Price(2);
       mPcc(1,i)=Pcc(1,i)*buy_Price(2);
       mPdd(1,i)=Pdd(1,i)*buy_Price(2);
       mPgrid(1,i)=buy_Price(2).*buy_Pgrid(1,i);
    elseif i>=18&&i<=22%峰
       objective=aaa+Pcharge(1,i)*0.00362+Pdischarge(1,i)*0.00182+Cgas*(aF(1)*P(1,i)+bF(1)*UP(1,i)+aF(2)*P(2,i)+bF(2)*UP(2,i))+aF(3)*P(3,i)+bF(3)*UP(3,i)...
           +Cgas*H(1,i)/H_n+buy_Price(3).*buy_Pgrid(1,i);
    aaa=objective;
    mPaa(1,i)=Paa(1,i)*buy_Price(3);
    mPcc(1,i)=Pcc(1,i)*buy_Price(3);
    mPdd(1,i)=Pdd(1,i)*buy_Price(3);
    mPgrid(1,i)=buy_Price(3).*buy_Pgrid(1,i);
    else i>=23&&i<=24;%谷
       objective=aaa+Pcharge(1,i)*0.00362+Pdischarge(1,i)*0.00182+Cgas*(aF(1)*P(1,i)+bF(1)*UP(1,i)+aF(2)*P(2,i)+bF(2)*UP(2,i))+aF(3)*P(3,i)+bF(3)*UP(3,i)...
           +Cgas*H(1,i)/H_n+buy_Price(1).*buy_Pgrid(1,i);
       aaa=objective;
  mPaa(1,i)=Paa(1,i)*buy_Price(1);mPcc(1,i)=Pcc(1,i)*buy_Price(1);mPdd(1,i)=Pdd(1,i)*buy_Price(1);
       mPgrid(1,i)=buy_Price(1).*buy_Pgrid(1,i);

欢迎感兴趣的小伙伴关注哦，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！