企业官网建设流程全解析

复现一篇“计及需求响应的区域综合能源系统双层优化调度策略研究”的核心期刊论文，最难的不是抄公式，而是把论文里的数学模型真正转成能跑的Matlab代码。我这两年断断续续做过几次类似的复现，说实话，双层优化这层窗户纸一旦捅破，代码本身并不复杂，卡住人的往往是建模思路和数学转化的细节。这篇权当是把我自己的完整复现路径摊开来讲一遍，里面既有模型拆解，也有Matlab实现骨架，还有几个我踩过不止一次的坑，希望能让正在做区域综合能源系统、双层优化调度、需求响应相关复现的同学少走点弯路。

先说个总判断：这类论文的复现，本质上是一个“三层翻译”工作。第一层，把文字和公式翻译成模型结构；第二层，把模型结构翻译成数学问题；第三层，把数学问题翻译成Matlab代码。大部分人一上来就直接进第三层，结果代码写了一大堆，模型边界却糊里糊涂，后面怎么调都不对。所以我下面的顺序也是按这个翻译流程来的。

1. 复现前先读懂“双层”：这篇论文到底在博弈谁和谁

1.1 上层在优化什么，下层在优化什么

区域综合能源系统（Regional Integrated Energy System，RIES）里的“双层”不是一个泛泛而谈的说法，它有明确的博弈结构。我在复现之前花了大半天时间只干了一件事：搞清楚上下层各自的决策主体、决策变量和目标函数。

按照这类论文最常见的设定：

• 上层决策者是综合能源系统运营商，它决定的是系统内部的设备出力，比如燃气轮机、燃气锅炉、电储能、电转气这些设备的运行状态，以及从外部电网买电、从气网购气的计划。目标函数一般是系统运行成本最小化，或者考虑碳交易、可再生能源消纳之后的综合成本最小化。
• 下层决策者通常是用能用户或者负荷聚合商，它在收到上层给出的能源价格信号后，通过调整自身用能行为来最小化用能成本或最大化用能效用，也就是我们常说的需求响应。
• 上下层之间的耦合变量就是价格信号和响应后的负荷曲线。上层先给出价格，下层基于价格调整负荷，调整后的负荷又会影响上层的运行成本和设备出力。这其实是一个Stackelberg博弈：上层是领导者，下层是跟随者。

这里有一个很关键的认知：双层模型要解决的不是“分别求解两个独立优化问题”，而是求一个“上层决策使得下层达到最优响应”的均衡解。你在复现时的所有数学处理，归根结底都是为了把这个“下层最优”翻译成上层问题里的一组约束。

1.2 上下层之间的耦合变量：价格、响应量与负荷曲线

复现的时候最容易模糊的地方就是耦合变量到底怎么定义。我建议在动手写代码前，先用一张表把耦合关系列清楚，比如下面这样：

为什么这一步重要？因为很多复现代码跑不出结果，就是耦合变量方向搞反了。比如有些论文把下层用户的负荷都当成固定参数，那下层就不存在决策了，整个模型其实是单层；有些代码把价格当成常数，那价格型需求响应也等于没建。读论文时如果发现作者在某个公式里同时出现上下层变量相乘，比如“价格乘以需求响应量”，那就得格外小心，这说明后面大概率有双线性项要处理。

1.3 复现前从论文中提取的公式清单

我的习惯是，拿到论文先做一张公式清单，把每个公式编号和它的作用列出来。比如：

• 设备约束：燃气轮机出力区间、爬坡约束、热电联产的电热耦合可行域；
• 能量平衡约束：电、热、气、冷各自的供需平衡；
• 需求响应约束：可平移负荷范围、可削减负荷比例、响应量上下限；
• 目标函数：购电购气成本、设备运行维护成本、需求响应补偿成本、碳排放成本等；
• 双层转化后的KKT条件或对偶约束。

这张清单的价值在于，它能让你在写代码时不用反复翻论文，每个约束对应一段代码，一段代码对应一组变量。复现失败的人大多是“边看论文边写代码”，写完发现漏了一堆约束，结果模型能跑但不收敛，或者结果和论文差很远。

2. 设备模型与需求响应的数学化：只有写得出公式才写得出代码

2.1 供能设备建模：CHP、燃气锅炉、储能与多能互补

综合能源系统的设备建模其实比单一电力系统的要繁琐，因为要同时处理电、热、气三种能量的转换。但好在多数论文用的都是线性或分段线性模型，Matlab配合Yalmip完全可以搞定。

以热电联产（CHP）机组为例，最常见的建模方式是用热电比耦合电出力与热出力：

P_gt = η_gt × F_gt × LHV

Q_chp = c_hr × P_gt

其中P_gt是燃气轮机发电功率，F_gt是天然气消耗速率，η_gt是发电效率，LHV是天然气低位热值，Q_chp是余热回收的供热功率，c_hr是热电比。

注意，大多数论文不会只用一条线性关系，而是给一个电热运行可行域，因为CHP的电热出力会互相制约。复现时最稳妥的做法是把可行域写成多边形约束，例如：

P_gt_min ≤ P_gt ≤ P_gt_max

Q_chp_min ≤ Q_chp ≤ Q_chp_max

Q_chp ≤ a×P_gt + b

这里的系数a和b由论文中的运行特性曲线决定，如果论文没给全，就得参考同类设备的典型参数，或者从图上的运行点反推。这一步我一般会在代码注释里标明“该参数为根据XX文献典型数据补充”，避免后面自己都忘了参数来源。

燃气锅炉的模型相对简单：

Q_gb = η_gb × F_gb × LHV

相当于把天然气化学能转换成热能，效率一般在0.85到0.95之间。

储能设备是另一个复现时特别容易写错的地方。电储能的荷电状态递推方程要写成：

SOC_{t+1} = SOC_t + (η_ch × P_ch_t - P_dis_t / η_dis) × Δt / E_cap

同时还要限制充放电功率上下限、储能SOC上下限，以及为了避免“既充又放”添加的互补约束，通常用0-1变量来做：

P_ch_t ≤ U_ch_t × P_ch_max

P_dis_t ≤ (1 - U_ch_t) × P_dis_max

我看到很多复现代码省略了这个0-1变量，结果就是储能一边充电一边放电，系统成本莫名其妙变低。这个问题在第5章会专门讲。

2.2 价格型与激励型需求响应的模型表述

需求响应在这个系统里不是简单地把负荷乘以一个系数，它分两种比较典型的建模方式。

价格型需求响应（Price-based DR，PDR）的经典做法是用价格弹性系数来描述负荷变化率与价格变化率的关系：

ΔL_t / L_base_t = e_tt × (Δρ_t / ρ_base_t) + Σ_{j≠t} e_tj × (Δρ_j / ρ_base_j)

其中e_tt是自弹性，e_tj是交叉弹性，ρ_t是时段t的能源价格。这里要特别注意：弹性系数不是越大越好，一般自弹性取-0.1到-0.3，交叉弹性取0.05到0.2，取值太夸张会让负荷跑到负值去。

激励型需求响应（Incentive-based DR，IDR）则是用户响应上层给出的削减激励，在约束范围内调整负荷。可以写成约束加成本的形式：

0 ≤ ΔP_cut_t ≤ U_cut_t × ΔP_cut_max

C_DR_t = a × (ΔP_cut_t)^2 + b × ΔP_cut_t

成本函数里的二次项会让问题变成二次约束规划（MIQP），很多论文为了稳定求解会把它做分段线性化。我建议复现时也这么做，一方面求解更快，另一方面分段线性化在和KKT条件结合时更方便处理。

2.3 多能负荷的耦合约束与热/电替代逻辑

综合能源系统里需求响应有个特点，就是“多能互补”。比如用户收到了高价电信号，不一定只能削减用电，还可以通过电锅炉、空气源热泵等设备在谷时段多用电、峰时段少用电。这种电热替代关系反映在模型里，通常是热负荷平衡方程里增加了一个“电制热”变量：

Q_load_t = Q_chp_t + Q_gb_t + P_eb_t × η_eb - Q_cut_t

其中P_eb_t是电锅炉消耗的电功率，η_eb是电转热效率，Q_cut_t是可削减热负荷。这种耦合约束让需求响应从“纯电负荷变动”升级成“综合用能替代”，也是这类核心期刊论文和普通单能系统优化文章的重要差别。复现时，多能替代的约束往往定义了模型的核心创新点，所以宁可写得保守一点，也不能写成两个解耦的模型。

3. 从主从博弈到可求解MILP：KKT条件与线性化全流程

3.1 写下层问题的KKT系统

双层优化模型没法直接用求解器跑，除非是形式非常特殊的问题。通用的处理思路是把下层优化问题用它的KKT最优性条件替换，从而把双层问题转化成带互补约束的单层问题。

假设下层用户的决策变量是x_L，目标函数是f_L(x_L; u)，其中u是上层传给下层的参数，约束是A_L x_L ≤ b_L。那KKT系统包含四部分：

• 驻点条件：∇f_L + A_L^T λ = 0
• 原始可行性：A_L x_L ≤ b_L
• 对偶可行性：λ ≥ 0
• 互补松弛条件：λ^T (b_L - A_L x_L) = 0

前三条都是线性不等式或等式，唯一的麻烦是互补松弛条件。因为它含有乘积项，是非线性的。

3.2 互补松弛条件的Big-M线性化

互补松弛条件的典型处理是引入0-1变量z和足够大的常数M，把每条互补条件拆成两组约束。例如对于第i条约束：

λ_i × (b_i - A_i x_L) = 0

可以等价转化为：

λ_i ≤ M × z_i

b_i - A_i x_L ≤ M × (1 - z_i)

其中z_i ∈ {0,1}

这个转化的逻辑是：如果z_i为0，则λ_i必须为0；如果z_i为1，则松弛量b_i - A_i x_L必须为0。这样就把乘积为0的条件拆成了线性不等式加整数变量。

这里Big-M的取值是个大学问。M取小了，会错误地砍掉可行解；M取大了，求解器会产生严重的数值病态，结果里出现大量0.9999和0.0001。我的经验是M不要统一拍脑袋取10^6，而是根据每个约束的物理量级单独设置。比如价格约束的M取50（因为电价一般不超过2元/kWh），功率约束的M取该设备最大容量的1.5倍。如果论文里数据和物理量级对不上，优先检查单位是不是kW和MW混用了。

3.3 双线性项的处理与目标函数最终形态

上层目标函数里经常出现上层决策变量乘以下层决策变量的情况，尤其是下层选择的是“用户用电量”这种和收益直接相关的变量，而上层决定的是“价格”时。价格乘以电量，天然是一个双线性项。

处理双线性项有三种常见手段：

第一种，如果价格是分时固定阶梯值，不是上层连续优化出来的，那这项其实是常数，直接预计算就行，常见于简化模型；

第二种，如果价格是上层可以连续调整的变量，那就要看双线性项是否可以通过变量代换变成凸二次项，有些论文会把双线性项写成用户效用函数的一部分，从而消掉乘法结构；

第三种，也是最常用的，就是设定上层在价格内部让利、下层响应后负荷变化只作用于上层平衡方程，而最终目标函数里只保留线性项。很多论文在数学推导的时候已经提前规避了双线性项，复现时要仔细看原公式到底有没有这个乘积项。如果确认存在并且没有被处理，那这个模型本质上是双层非线性优化，求解难度会上一个台阶，需要非常小心地设计求解策略。

3.4 单层MILP的求解复杂度与求解器选择

单层化之后，模型就变成了一个混合整数线性规划（MILP）或者混合整数二次规划（MIQP），规模取决于时段数、设备数、需求响应约束条数。

一个典型的24时段调度模型，时段数为24，如果CHP、锅炉、储能、需求响应各几个约束，再加上KKT转化带来的0-1变量，变量总数很容易到几百上千，整数变量也有上百个。好在现在的商业求解器和开源求解器都吃得下这个规模。我在Matlab里一般用Yalmip做建模层，求解器选Gurobi或Cplex，如果只有Matlab自带求解器，可以考虑用“intlinprog”，但遇到大规模问题会比较吃力。

这里特别提醒一下：不要轻易在复现中尝试用fmincon处理这种双层单层化模型。因为KKT互补约束让问题非光滑，fmincon对初值极度敏感，经常卡在局部最优甚至不收敛。MILP求解器至少能保证在给定MIP Gap下找到全局近似最优解，这才是主流论文的做法。

4. Matlab代码实现：Yalmip+Gurobi搭建双层调度完整框架

4.1 环境与工具箱配置

我的复现环境是Matlab R2021b + Yalmip + Gurobi 9.5，操作系统Windows和Ubuntu都跑过。Gurobi的学术许可证申请很快，装好后在Matlab里运行gurobi_setup即可。Yalmip是纯Matlab工具箱，解压后加路径就行。

代码组织上我喜欢按模块拆分，而不是几百行塞在一个脚本里。一个典型的文件结构如下：

IES_DOUBLE_LAYER/ ├── main_upper.m % 主程序入口，运行上层调度 ├── model_upper.m % 上层模型：设备约束、目标函数 ├── model_lower_kkt.m % 下层模型KKT条件 ├── load_data.m % 读取负荷、价格、设备参数 ├── plot_results.m % 绘图和结果分析 └── data/ ├── load_curve.xlsx % 典型日电负荷、热负荷 ├── price_data.xlsx % 分时电价、气价 └── device_params.m % 设备参数

4.2 上层调度模型的代码骨架

上层模型核心代码的大致结构能写成这样，注意这是简化示意，参数要根据你的论文调整：

%% load_data.m 中读取参数后 % 24时段典型日 T = 24; dt = 1; % 时间间隔，单位h % 设备参数 P_gt_max = 2000; % CHP最大电出力，kW P_gt_min = 200; % CHP最小电出力 eta_gt = 0.35; % 发电效率 c_hr = 1.2; % 热电比 Q_gb_max = 1500; % 燃气锅炉最大热出力 Soc_max = 1000; % 储能容量 kWh P_ch_max = 300; % 最大充电功率 P_dis_max = 300; % 最大放电功率 eta_ch = 0.95; eta_dis = 0.95; %% 定义决策变量 P_gt = sdpvar(1, T); % CHP电出力 Q_chp = sdpvar(1, T); % CHP热出力 Q_gb = sdpvar(1, T); % 燃气锅炉热出力 P_grid = sdpvar(1, T); % 外购电功率 Soc = sdpvar(1, T); % 储能SOC P_ch = sdpvar(1, T); % 充电功率 P_dis = sdpvar(1, T); % 放电功率 Constraints = []; % CHP约束 Constraints = [Constraints, P_gt_min <= P_gt <= P_gt_max]; Constraints = [Constraints, Q_chp >= 0, Q_chp <= c_hr * P_gt]; % 这里要留意：热电耦合可行域是从论文中提取的，不一定是简单不等号 % 燃气锅炉 Constraints = [Constraints, 0 <= Q_gb <= Q_gb_max]; % 储能 Constraints = [Constraints, Soc(1) == 0.2 * Soc_max]; Constraints = [Constraints, Soc(2:T) == Soc(1:T-1) ... + (eta_ch * P_ch(2:T) - P_dis(2:T) / eta_dis) * dt / Soc_max]; Constraints = [Constraints, 0 <= Soc <= 1]; Constraints = [Constraints, 0 <= P_ch <= P_ch_max]; Constraints = [Constraints, 0 <= P_dis <= P_dis_max]; Constraints = [Constraints, P_ch + P_dis <= P_ch_max]; % 简化同时充放限制 % 电功率平衡 P_load = data.P_load; % 原始电负荷 Constraints = [Constraints, P_gt + P_wt + P_pv + P_grid + P_dis ... == P_load + P_ch + P_heat_pump]; % 目标函数：购电成本 + 购气成本 + 运维成本 + 需求响应补偿 Cost_buy_e = sum(price_e .* P_grid); Cost_buy_g = sum(price_g .* (P_gt / eta_gt + Q_gb / eta_gb) / LHV); ... Objective = Cost_buy_e + Cost_buy_g + Cost_om + Cost_dr; ops = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(Constraints, Objective, ops);

4.3 下层KKT模型在代码里的实现方式

下层KKT条件翻译成代码有两种写法。一种是显式写出对偶变量和Big-M约束，但这要求你手工推导出所有互补条件，非常繁琐而且容易写错符号。另一种是直接用Yalmip的“implies”或“iff”关系来表示逻辑约束，Yalmip会自行转化为混合整数约束。

举个例子，下层用户有一个决策变量是削减量delta_L，约束是0 ≤ delta_L ≤ delta_L_max，且削减会带来负效用，下层目标可以简化为min lambda * delta_L + c_dr * delta_L^2。其KKT条件翻译时，可以用Yalmip写：

delta_L = sdpvar(1, T); lambda_L = sdpvar(1, T); % 对偶变量 mu_L = sdpvar(1, T); % 对偶变量 % 下层驻点条件，简化示意 Constraints = [Constraints, lambda_L >= 0, mu_L >= 0]; Constraints = [Constraints, lambda_L + cost_coef == mu_L]; % 互补约束：lambda_L 和 (delta_L_max - delta_L) 互补 Constraints = [Constraints, implies(lambda_L > 0, delta_L == delta_L_max)]; Constraints = [Constraints, implies(delta_L < delta_L_max, lambda_L == 0)];

注意，implies里比较表达式严格大于0的写法在数值上可能不太干净，实际使用时更稳妥的是引入0-1变量，把互补条件写成：

z1 = binvar(1, T); z2 = binvar(1, T); Constraints = [Constraints, lambda_L <= M * z1]; Constraints = [Constraints, delta_L_max - delta_L <= M * (1 - z1)]; Constraints = [Constraints, delta_L <= M * z2]; Constraints = [Constraints, lambda_L <= M * (1 - z2)];

哪种自己看着舒服就用哪种。我实际更喜欢显式Big-M写法，因为M的取值可控，而implies在某些Yalmip版本里生成的中间变量让人难以排查数值问题。

4.4 数据组织与运行框架

数据部分的坑比想象中多。典型的输入数据包括：

• 典型日负荷曲线：电负荷、热负荷、冷负荷，通常取春夏秋冬四个典型日或夏季/冬季两个典型日；
• 分时电价：峰平谷时段划分和各时段电价；
• 天然气价格：按单位体积或单位能量计价均有可能，注意统一换算；
• 新能源出力曲线：风电、光伏24小时出力标幺值；
• 设备参数表：容量、效率、爬坡率、运维成本系数。

我习惯把所有数据放到一个Excel里，load_data.m统一读取，然后转成数组。这样后续换一组数据做灵敏度分析非常方便。表格的设计尽量和论文中Table对应，比如设备参数表可以这样组织：

主程序的运行逻辑是：先load_data，然后建立上层模型，再把下层KKT条件加进去形成单层MILP，求解，最后把结果画出来。由于有0-1变量，一次完整求解可能从几秒到几分钟不等，我一般会把gurobi的输出级别设低，避免刷屏。

5. 复现踩坑实录：五个最容易翻车的细节与修正

5.1 模型无界或不可行的排查思路

这类问题在双层优化复现里出现频率最高。模型无界，先检查目标函数里有没有变量没写成本项但约束又没限制它，比如漏了SOC初值，储能就会无限放电来“赚钱”。模型不可行，大概率是互补松弛条件的方向反了，或者KKT符号写错。

我的排查方法是分步验证：先只跑上层模型，把下层响应量全部固定为常数，看看模型能不能正常求解；如果能，再加下层KKT约束；加了之后报不可行，那问题基本100%在KKT转化那一段。然后再逐条注释互补约束，缩小范围。

5.2 Big-M数值病态问题

Big-M取太大直接导致求解器报“Numerical trouble”或者结果出现大量小数。一个看起来很吓人的报错就像下面这样：

Warning: Model may be infeasible or unbounded.

遇到这种，别急着改模型，先检查所有M值。我的经验是M要贴合具体约束的量级，比如负荷变量的M可以用峰值负荷的2倍，价格变量的M用最高价格的2倍。如果模型中已经用了很多0-1变量，还可以把求解器参数里的NumericalFocus设成2，让Gurobi花更多精力处理数值条件。

5.3 需求响应量脱离物理约束

有些复现结果里，需求响应后的负荷曲线会出现明显的“过度削峰填谷”，峰时段的负荷被砍到负值，或者谷时段负荷暴涨，这说明需求响应模型少了物理约束。要记住三条兜底约束：

第一，负荷在任意时刻都不能为负；第二，可削减负荷总量有上限，不能一天之内把工业负荷全砍掉；第三，可平移负荷要满足总量守恒，即平移前后总电量不变。如果论文里没有明确写这三条，复现时也应该主动补上，因为不加这些约束的模型生产出来的结果在工程上是没有意义的。

我在代码里通常这么加：

Constraints = [Constraints, 0 <= P_load_after <= 1.2 * P_load_base]; Constraints = [Constraints, sum(P_load_after) >= 0.95 * sum(P_load_base)]; Constraints = [Constraints, sum(P_load_after) <= 1.05 * sum(P_load_base)];

5.4 主从迭代不收敛

如果你没有走KKT单层化路线，而是用“上层优化-下层优化-更新价格-再上层优化”这种启发式迭代，那大概率会碰上不收敛。原因很简单：上下层之间是强耦合的博弈问题，交替迭代本质上是一个不动点迭代，如果没有松弛因子，很容易在两点之间震荡。

如果确实要用迭代法，建议加入权重更新，例如下一次迭代的上层负荷取“本次下层响应结果”和“上次结果”的加权平均：

P_load_{k+1} = α × P_response_{k} + (1 - α) × P_load_k

其中α取0.3到0.5比较稳。但归根结底，我建议优先做单层化，KKT方法至少收敛性有保障。

5.5 复现结果与论文对不上时的检查顺序

有时候模型能跑，结果趋势也对，但数值和论文差不少。我的检查顺序是：

1. 单位是否统一：论文里的热负荷是MW还是kW，电价是元/kWh还是元/MWh，气价是按体积还是按热值；
2. 典型日数据是否一致：有些论文用的是春季典型日，你用了夏季典型日，结果当然不一样；
3. 是否考虑了网络损耗或线路容量约束：一些论文在综合能源系统里加入了天然气网络和热网管道约束，如果漏掉会直接影响购电购气量；
4. 目标函数里是否含碳交易成本：这个项经常在正文公式里出现，但在摘要或图表里不明显，很容易漏；
5. 下层用户效用函数的参数：弹性系数、补偿系数这些论文可能只给了一部分，复现时需要用敏感性分析来匹配，找到一组和论文结果最接近的参数值。

最后聊两句我自己的习惯

复现这类双层优化调度，我会把论文里的每个公式编号都写在代码注释里，哪怕代码已经精简过，也会保留原始公式编号，比如% Eq.(12): CHP feasible region。这个习惯救过我很多次，因为来回翻论文查公式的效率实在太低。另外就是每调完一个版本，把生成的图表和参数存一份，记录下当时的求解器设置和M值，这样对比不同参数下的结果时能快速定位变动来源。

这个方向后续如果想深入，可以加的东西还有很多，比如把单层的MILP扩展成考虑多场景随机优化的两阶段模型，或者在目标函数里引入碳捕集与电转氢的耦合。但前提都是先把双层调度这个框架搞清楚、跑通，它才是所有扩展的底盘。