参考文献：
[1]张浩鹏,李泽宁,薛屹洵,等.基于共享储能服务的智能楼宇双层优化配置[J/OL].中国电机工程学报,1-12[2024-05-22].
1.摘要
为降低城市化进程中楼宇储能投资成本，提出一种基于共享储能服务的智能楼宇（Intelligent Buildings，IBs）双层优化配置方法。首先，建立计及寿命周期的共享储能电站（Shared Energy Storage Station，SESS）模型。其次，基于楼宇建筑物热惯性，构建含空调系统的 IBs 数学模型。然后，综合考虑 SESS 与 IBs 的差异化利益诉求，建立基于 SESS 的 IBs 双层优化模型。上层模型目标函数旨在降低 SESS 的规划成本，下层模型目标函数旨在降低IBs 的 年 运 行 成 本 ， 并 采 用 卡罗需 - 库 恩 - 塔 克（Karush-Kuhn-Tucher，KKT）条件将原双层优化问题转换为单层混合整数线性规划问题进行求解。最后，以三个IBs 社区的四季典型日为例，对比分析了不同优化配置方法对于 IBs 运行和 SESS 配置结果的影响。结果表明，在满足 IBs 用户温度舒适性的同时，所提双层优化配置方法可充分满足 SESS 运营商与 IBs 的差异化利益诉求，实现双方的共赢。
2.原理介绍
SESS-IBs 双层优化框架如图 1 所示，包括SESS、N 个 IBs 和配电系统。每个 IBs 社区可由本地光伏与风电、SESS 和电网供电，其能耗主要来自于常规电负荷和空调负荷。考虑到 IBs 用户受到相关政策和技术的限制，设定 IBs 用户不能售电给电网。IBs 通过向 SESS 缴纳服务费参与共享储能服务，满足自身的储能需求，同时节省了巨大的储能投资成本。
SESS 内设置储能电站调控中心，根据每个时段 IBs 的充放电需求下达调度指令，在同一时段内接入 SESS 的所有 IBs 总充放电需求若为充电，则储能电站调控中心使用储能装置充电来存储 IBs 的多余电能；若为放电，则储能电站调控中心使用储能装置放电来满足 IBs 的用能需求，IBs 用能以购售电的形式进行结算。储能电站调控中心根据一个周期内所有 IBs 的用电行为来配置相应的容量和最大充放电功率。SESS 根据分时购售电电价向用户支付或收取相应的费用，主要以不同时段购售电的价格差以及向 IBs 收取服务费实现盈利，其中服务费均以单位功率收取。
2.1 上层共享储能电站规划模型
2.1.1 上层模型目标函数
基于共享储能服务的 IBs 双层优化模型以降低 SESS 储能投资成本和提高 SESS 运行经济效益为上层目标，其中上层目标函数包括投资和运维成本、SESS 和 IBs 交易电量的费用和从 IBs 收取的共享储能服务费。
2.1.2 上层模型约束条件
SESS 需要满足的约束条件包括荷电状态约束、充放电功率约束、储能电站倍率约束和储能电站充放电功率平衡约束。
1）荷电状态约束
2.2 下层智能楼宇优化运行模型
2.2.1 下层模型目标函数
基于共享储能服务的 IBs 双层优化模型以降低 IBs 的年运行成本为下层目标，其中下层目标函数包括从电网购电成本、IBs 和 SESS 交易电量的费用和向 SESS 缴纳的服务费。
2.2.2 下层模型约束条件
下层 IBs 需要满足的约束条件包括墙体热平衡约束、室内热平衡约束、电功率平衡约束、空调运行约束、爬坡约束、温度舒适区间约束、从电网购电功率约束和 IBs 与 SESS 间购售电功率约束。具体约束如下：
1）墙体热平衡约束
基于热容-热阻网络模型可构建考虑建筑物热惯性的 IBs 详细热动态模型。以 IBs 单个室内区域为例，楼宇建筑物的热动态模型如图 2 所示。
热动态模型描述的是单个制热/冷区域，楼宇模型则是由多个类似构造的区域聚合而成。本文近似每栋楼宇内的制热/冷区域构造一致，因此在同种场景及相同环境参数下，楼宇群中每个制热/冷区域的空调系统功率是一致的。IBs 所有室内区域的围护结构均相同，楼宇围护结构的热平衡约束可以描述为：
2）室内热平衡约束
室内温度受到多种因素的影响，包括外界环境温度、光照强度、墙体温度和空调系统制冷/热等。室内热平衡约束如下所示：
3）电功率平衡约束
4）空调运行约束
空调系统主要通过制冷/热来满足 IBs 用户的温度舒适性，其运行功率如下所示：
5）爬坡约束
空调运行功率变化应在一定的区间范围内，其爬坡约束具体如下：
6）温度舒适区间约束
满足 IBs 用户的温度舒适性应在一定的区间范围内，该温度舒适性约束如下：
7）从电网购电功率约束
8）IBs 与 SESS 购售电功率约束
2.3 求解方法
本文构建的双层模型目标函数中含有绝对值，为非线性项，并且上层模型和下层模型存在耦合关系，难以直接进行求解。
因此，先将目标函数中的绝对值项线性化。再通过构建下层 IBs 模型的拉格朗日函数，利用卡罗需-库恩-塔克（Karush-Kuhn-Tucher，KKT）条件代替下层模型，同时作为上层模型的附加约束，将原双层优化问题转化为单层非线性优化问题。最后采用Big-M 法将非线性约束进行线性化，问题最终变为单层混合整数线性规划问题，求解过程具体如图 3 所示。
3.编程思路
3.1参数和变量定义
3.2编程思路
原文献中对下层优化的KKT方程进行了详细的推导。但实际上yalmip工具箱中已经存在kkt命令可以直接求解得到优化问题的kkt方程。
此外，代码中采用了yalmip工具箱+gurobi求解器进行编程，需要安装求解器才能得到正确的结果。如果没有gurobi但是有其他求解器，也可以在设置中进行修改。
4.部分代码
%% 清除内存空间
clc
clear
close all
warning off%% 系统参数
NT = 24; % 控制时段数
Nw = 4; % 典型日数目
Tw = 91; % 每个典型日的天数
NIB = 3; % IBs社区数
tau = 0; % 自放电效率
n_abs = 0.98; % SESS 的充电效率
n_relea = 0.98; % SESS 的放电效率
dt = 1; % 单位时间间隔
SOC_min = 0.1; % SOC最小值
SOC_max = 0.9; % SOC最大值
SOC0 = 0.2; % SOC初始值
beta = 2.8; % SESS 能量倍率
Cwall = [7.9e5*ones(3,1);2.6e7]; % 墙体热容
Rwall = [0.06*ones(3,1);0.08]; % 墙体热阻
vij = [0.45;0.35]; % 墙体吸热率
Awall = [12;15]; % 墙体面积
Qrad = [1000;1500]; % 墙体对应方向的光照强度
Croom = 2.5e5; % 房间内热容
wij = 0.9; % 窗体透射率
Awin = 6; % 窗体面积
E_EER = 3; % 空调能效比
PAC_max = 1.6; % 空调系统运行功率的上限
Ppb_max = 0.2; % 空调运行功率变化的上限
Ppb_min = -0.2; % 空调运行功率变化的下限
Tin_max = 24; % IBs 用户舒适温度区间的上限
Tin_min = 20; % IBs 用户舒适温度区间的下限
Pgrid_max = 4000; % IBs从电网购电的最大功率
Pess_mg_max = 4000; % IBs 和 SESS 进行电量交易的最大功率
Ppv = zeros(NT,Nw,NIB); % 第w个典型日第n个IBs社区在t时刻的光伏功率
Pwt = zeros(NT,Nw,NIB); % 第w个典型日第n个IBs社区在t时刻的风电功率
Pload = zeros(NT,Nw,NIB); % 第w个典型日第n个IBs社区在t时刻的负荷需求
data1 = xlsread('风光负荷数据.xlsx', 'IB1');
data2 = xlsread('风光负荷数据.xlsx', 'IB2');
data3 = xlsread('风光负荷数据.xlsx', 'IB3');
Pload(:,:,1) = data1(:,[2,5,8,11]);
Pload(:,:,2) = data2(:,[2,4,6,8]);
Pload(:,:,3) = data3(:,[2,5,8,11]);
Ppv(:,:,1) = data1(:,[3,6,9,12]);
Ppv(:,:,2) = data2(:,[3,5,7,9]);
Ppv(:,:,3) = data3(:,[3,6,9,12]);
Pwt(:,:,1) = data1(:,[4,7,10,13]);
Pwt(:,:,3) = data3(:,[4,7,10,13]);Tout = ... % 室外温度
[26.50 14.71 8.17 -8.15
26.36 13.57 7.32 -8.31
26.22 11.58 6.18 -8.47
26.07 9.45 4.48 -8.48
25.93 8.31 4.33 -8.79
25.22 7.60 3.05 -9.26
24.80 8.45 1.07 -9.42
24.94 10.02 3.05 -9.58
26.36 12.15 7.32 -8.35
26.50 14.28 10.44 -7.12
28.35 18.54 12.43 -6.20
28.35 20.25 13.29 -5.59
30.20 21.95 13.57 -5.59
31.47 22.24 13.71 -5.60
32.18 23.09 13.85 -6.07
33.89 23.66 14.14 -6.38
33.61 22.24 13.85 -6.85
33.89 22.10 13.14 -7.48
32.47 21.67 12.15 -7.33
31.47 20.11 11.01 -8.11
30.76 18.69 10.16 -8.27
29.91 17.69 9.45 -8.74
29.20 16.84 9.02 -8.44
29.20 16.27 8.88 -8.91];%% 上层决策变量
Pess_max = sdpvar(1); % SESS的最大充放电功率
Eess_max = sdpvar(1); % SESS的最大容量
Eess_t = sdpvar(NT + 1,Nw); % t时刻SESS的电能储量
Pess_t = sdpvar(NT,Nw); % t时刻SESS和IBs交易电量的总功率
Pess_abs = sdpvar(NT,Nw); % t时刻SESS的充电功率
Pess_relea = sdpvar(NT,Nw); % t时刻SESS的放电功率%% 下层决策变量
Pess_bs = sdpvar(NT,Nw,NIB); % 第w个典型日第n个IBs社区在t时刻与SESS交易电量的功率
Pgrid = sdpvar(NT,Nw,NIB); % 第w个典型日第n个IBs社区在t时刻从电网购电的功率
Twall_12 = sdpvar(NT,Nw,NIB); % 墙体温度
Twall_13 = sdpvar(NT,Nw,NIB); % 墙体温度
Twall_14 = sdpvar(NT,Nw,NIB); % 墙体温度
Twall_15 = sdpvar(NT,Nw,NIB); % 墙体温度
Tin = sdpvar(NT,Nw,NIB); % 房间内温度
Qint = sdpvar(NT,Nw,NIB); % 房间内部热源
PAC = sdpvar(NT,Nw,NIB); % 第w个典型日第n个IBs社区的房间在t时刻的空调功率%% 上层约束条件%% 目标函数[KKTsystem , details] = kkt(Constraints_down,objective_down,[Pess_bs;Pgrid;Twall_12;Twall_13;Twall_14;Twall_15;Tin;Qint;PAC]);%% 设求解器
% gurobi求解器
ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1);
ops.gurobi.TimeLimit = 600; % 运行时间限制为10min
ops.gurobi.MIPGap = 0.01; % 收敛精度限制为0.01% cplex求解器
% ops = sdpsettings('verbose', 3, 'solver', 'cplex','showprogress',1,'debug',1);
% ops.cplex.timelimit = 600; % 运行时间限制为10min
% ops.cplex.mip.tolerances.mipgap = 0.01; % 收敛精度限制为0.01% mosek求解器
% ops=sdpsettings('verbose', 3, 'solver', 'MOSEK','cachesolvers',1);
% ops.mosek.MSK_DPAR_OPTIMIZER_MAX_TIME=600;% 运行时间限制为10min
% ops.mosek.MSK_DPAR_MIO_TOL_REL_GAP=0.01; % 收敛精度限制为0.01
sol = optimize([KKTsystem,Constraints_up,Constraints_down], objective_up, ops);%% 分析错误标志
if sol.problem == 0disp('求解成功');
elsedisp('运行出错');yalmiperror(sol.problem)
end%% 显示结果
show_result;完整代码获取方式：
(文章复现)基于共享储能服务的智能楼宇双层优化配置matlab代码资源-CSDN文库
5.代码运行结果
