简介：面向工业园区负荷聚合商开展日内需求响应的场景，提供一份基于Matlab与Yalmip实现的需求响应资源聚合优化配置代码，可复现《电网技术》2022年相关文献的聚合建模与求解流程，适合电力系统优化方向的研究生、工程师借鉴。压缩包共4个文件，包含两个.m脚本、一个.p文件与一份PDF代码说明，整体仅689KB，结构精简、便于快速部署。已有395人学习浏览，常用于课题复现、算法对比与课程设计参考。代码涵盖日前聚合优化与日内调用两阶段逻辑，借助Yalmip完成建模与求解，读者可据此掌握大规模分散资源的聚合方法、目标约束构建思路，并快速迁移到类似园区级需求响应问题中展开扩展实验。 如果你也和我一样，面对过这样一篇文章——题目写着“工业园区需求响应资源聚合优化配置方法”，点进去一看，满屏的都是双层优化、混合整数规划、需求响应潜力评估这些让人头皮发麻的术语，而你的任务偏偏是把它复现成一套能在Matlab里跑起来的代码，那这篇文章就是写给你看的。
我最早接到这个类型的复现需求时，犯了一个典型错误：打开论文附件就开始对着公式敲代码。结果三天过去，代码写了八百行，运行报错二十几次，连论文里的模型到底决策哪些变量都没理清楚。后来我把整套复现流程推翻重来，换了个思路——先翻译模型，再搭工程框架，最后才是填代码。这个过程走顺之后，类似的优化配置类论文基本都能在几天之内跑通。
这篇文章不会去复读论文本身，因为每篇具体论文的目标函数和约束条件多少有差异。我以工业园区里最常见的需求响应资源配置为例，把从“看论文公式”到“Matlab跑出优化结果”的完整链路拆开讲清楚，包括代码架构怎么搭、Yalmip建模怎么处理双线性项和整数变量、求解器怎么选、论文图表怎么复刻，以及在复现过程中我踩过的那些坑。
1. 复现论文先别急着写代码：把“公式-模型-代码”三层翻译顺
很多刚接触论文复现的人最大的问题，不是不会写代码，而是跳过模型理解，直接试图把LaTeX公式一行行翻译成Matlab语句。这种做法在简单线性模型上勉强能成，一旦遇到需求响应这种多设备、多时间尺度、含整数决策变量的优化问题，一定会被卡死。
1.1 先搞清楚论文在优化什么：一个运行逻辑的梳理模板
拿到论文后，我建议先花两个小时做一件事：把论文里的“优化主体”和“决策变量”单独抽出来，写在一张纸上。以这篇题目一般涉及的内容来说，园区里的需求响应资源通常包括以下几类：
- 可削减负荷（空调、照明等柔性负荷，切掉一部分不影响生产的主流程）
- 可转移负荷（某些工业流程可以在不同时段之间平移，但总用电量不变）
- 储能系统（充放电功率和SOC状态都在优化范围内）
- 光伏或分布式电源（出力曲线可以作为边界条件，也可以作为决策变量）
怎么判断论文到底把哪些设备作为决策变量？最快的办法是看摘要里的关键词是“最小化运行成本”还是“最大化新能源消纳”，然后去结论部分找优化后的结果图。结论图里出现了哪几条曲线，论文的优化对象就有几类。
以我复现过的常见框架为例，典型的目标通常不是单一的经济目标，而是“经济性+削峰填谷”的组合目标。经济性一般是园区从电网购电的成本，削峰填谷则是让负荷曲线尽量平缓，通常体现为峰时负荷最小化或者负荷方差最小化。两个目标加权组合，权重系数论文里一般会给出几组对比值。
1.2 从结论图表里反推算例构成
这是我觉得最实用的一招：翻到论文倒数几页的算例分析，看结果图里到底画了什么，再反推代码需要准备什么数据。比如结果图里有一张“优化前后负荷曲线对比”，那说明你需要准备基础负荷数据；如果有一张“储能SOC变化曲线”，那说明储能模型一定在约束里；如果有一张“各时段可削减负荷调用量”的堆积图，那就是可削减负荷参与了优化。
我当时照着这个方法，从论文图表里反推出来的算例构成包括：24小时的园区基础负荷曲线、分时电价、光伏出力、储能参数，以及可削减负荷的容量和补偿成本系数。把这些东西整理成一个Excel数据文件，后面的建模工作就有了具体的对象。
2. 优化模型的解构与数学化：目标函数和约束怎么写进代码
模型翻译是复现工作的核心，这一步做好了，后面的代码只是体力活。我在这一步采用的方式是：先把论文里的每个公式编号，然后逐个确定它在代码里对应的是目标函数、约束条件还是参数定义，最后再做维度统一。
2.1 目标函数：经济性主导时怎么处理非线性项
以某篇典型的工业园区需求响应论文为例，目标函数可以表述为：
[ \min \sum_{t=1}^{T} \left[ C_t^{buy} P_t^{grid} \Delta t + \sum_{i} (a_i P_{i,t}^{cut} + b_i) \right] ]
其中第一项是购电成本，第二项是调用可削减负荷的补偿成本。这个函数本身是线性的，但如果你遇到的目标函数里有[P_{cut}^2]这类平方项，直接用Yalmip加进去就会导致模型变成二次规划。二次规划对中小规模问题还能解，但论文里如果真的用二次项表示设备损耗或惩罚，复现时建议先确认原文是用了二次规划求解器还是做了一次线性化处理。
从可复现性的角度，我通常优先按论文原文的形式建模，然后用Yalmip自带的求解器识别问题类型，再决定配置哪类求解器。
2.2 约束条件的类型与矩阵化表达
对需求响应聚合优化问题而言，约束条件通常可以分成四组：
功率平衡约束：任意时段的电网购电+光伏出力+储能放电 = 基础负荷+可转移负荷+储能充电+可削减负荷的削减前功率。
设备物理约束：储能的充放电功率上限、SOC上下限、充放电状态互斥；可削减负荷的单次削减比例上限、削减次数上限。
需求响应约束：被削减的总量不能超过合同约定的响应容量，或者响应后的负荷不能低于某个安全值。
系统运行约束：与电网交互的功率不能越限。
把约束条件全部列出来之后，下一步是检查每个约束里的变量是否都有明确的上下界。这是新手最容易漏掉的地方——变量没有边界，求解器会直接告诉你“模型无界”或者给出一个离谱的结果。
2.3 单位、时段、上下界的统一约定
复现过程中最伤害代码质量的不是模型复杂，而是单位不统一。很多论文把功率单位写成MW，但电量单位用MWh，时间步长可能是一小时也可能是一刻钟。写进Matlab的时候，必须统一成一套体系：功率用kW、电量用kWh、时间步长用小时，价格用元/kWh。
另外，时段数的确定也很关键。大多数论文用24个时段，每段一小时；如果是15分钟一个点，那T就是96。我通常建议代码里的T做成参数而不是写死，这样复现不同论文的不同分辨率时，只改一个数就行。
我在代码里维护了一个结构体params，里面存了时段数T、设备参数和电价序列，凡是后续函数要用到的常量一律通过params传入，而不是在多个函数里重复硬编码。这样的好处是调数据时只需要动一个文件。
3. Matlab工程架构：一套能跑通、能改参、能复现结果的代码组织方式
写优化代码和写一般仿真代码不一样，它的特点是模型和数据耦合度高，求解失败时排查链路长。如果你把一千行代码写在同一个main脚本里，报错的时候你根本分不清是数据问题还是约束写错。
3.1 模块化文件组织与执行流程
我复现这类论文时使用的文件结构如下：
industrial_park_project/ ├── main.m ├── config/ │ └── config_params.m ├── data/ │ ├── load_data.xlsx │ ├── price_data.xlsx │ └── pv_data.xlsx ├── model/ │ ├── build_optimization_model.m │ ├── objective.m │ └── constraints.m ├── solver/ │ └── solve_model.m └── plot/ └── plot_results.mmain.m的角色只是按顺序调用各模块，保证整个目录下任何一个文件单独打开都能看懂它负责什么。
config_params.m负责读数据和定义全局参数。model下的三个文件分别构建目标函数和约束，solver里调用Yalmip的optimize函数，plot负责把结果画成论文里的样子。
这种架构最大的好处是：论文里的每个公式都能在代码里找到对应位置，导师问你这个约束怎么实现的，你直接打开constraints.m指给他看就行。
3.2 数据结构设计：用结构体统一管理时变参数
我习惯把时段性参数（电价、光伏、基础负荷）和静态参数（储能容量、效率、上下限）分开放到两个结构体里。时段性参数用列向量表示，长度与T相同；静态参数用标量。这样定义约束的时候，写P_grid >= 0和P_grid <= P_grid_max都特别直观，不用到处查下标。
举个典型的配置例子：
% config_params.m T = 24; % 时段数，1小时一个点 params.dt = 1; % 时间步长（小时） % 分时电价，元/kWh params.price = [ 0.33 * ones(8,1); % 谷段 00:00-08:00 0.68 * ones(6,1); % 平段 08:00-14:00 1.10 * ones(6,1); % 峰段 14:00-20:00 0.68 * ones(4,1); % 平段 20:00-24:00 ]; % 储能参数 params.battery_capacity = 4000; % kWh params.battery_power_max = 1000; % kW params.soc_min = 0.2; params.soc_max = 0.9; params.eta_ch = 0.95; % 充电效率 params.eta_dis = 0.95; % 放电效率 % 可削减负荷参数 params.load_cut_max = 0.15; % 最多削减基础负荷的15% params.cut_cost = 0.45; % 补偿单价，元/kWh这些参数我没有完全照搬某一篇论文，而是按行业常见量级给出的参考值，具体复现时应该换成目标论文的数据。
3.3 为什么我优先用Yalmip而不是手写优化求解
Matlab里做优化有三条路：直接调linprog/intlinprog、手写内点法或单纯形法、用Yalmip建模再调外部求解器。
我的建议非常明确：用Yalmip。原因很简单——需求响应聚合优化本质上是混合整数线性规划或混合整数二次规划，手写求解器既不稳定也不现实；linprog只能处理纯线性且退化能力有限，而Yalmip帮你把模型描述和求解器解耦了，建模写的是数学表达式，求解时才决定用哪个求解器。
% main.m 中初始化Yalmip环境 yalmip('clear') ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'showprogress', 1);后面详细写建模的时候还会展开说明怎么用Yalmip表示约束。
4. 核心代码片段拆解：优化求解部分是怎么写出来的
这一节直接上干货。以“含储能+可削减负荷的工业园区日前优化调度”为蓝本，我把整个优化模型的关键代码片段过一遍。如果你的目标论文里还有可转移负荷或电锅炉，在这个框架上加约束就行。
4.1 决策变量定义与维度检查
% model/build_optimization_model.m P_grid = sdpvar(T, 1); % 电网购电功率，kW P_ch = sdpvar(T, 1); % 储能充电功率，kW P_dis = sdpvar(T, 1); % 储能放电功率，kW SOC = sdpvar(T+1, 1); % 储能SOC状态，保留T+1个点（含初始） P_cut = sdpvar(T, 1); % 可削减负荷功率，kW u_ch = binvar(T, 1); % 充电状态0/1变量 u_dis = binvar(T, 1); % 放电状态0/1变量定义完变量后先别急着写约束，用size检查一下每个变量的维度是否和预期的(T,1)一致。这一行检查能避免后面大量因维度问题引发的报错。
4.2 功率平衡约束的写法
功率平衡是每篇优化论文都会有的核心约束，它描述的物理含义是：园区里所有“进”的功率等于所有“出”的功率。在Matlab代码里它不是一个等式，而是一组按时间展开的等式或不等式。
% 功率平衡：购电+光伏+放电 = 基础负荷-削减+充电 P_load = params.load_basic; % 基础负荷序列，kW P_pv = params.pv_output; % 光伏出力序列，kW Constraints = []; Constraints = [Constraints, P_grid + P_pv + P_dis == P_load - P_cut + P_ch];这里有个细节要提醒：P_cut表示的是”削减掉的功率“，还是“削减后依然要供应的功率”？不同论文的符号定义不一样。我的习惯是P_cut表示被削减掉的那部分功率，所以等式右边基础负荷要减去它。写约束前先在注释里明确自己的符号约定，否则后面画图的时候会乱。
4.3 储能SOC的状态转移和逻辑约束
储能是需求响应里最灵活的设备，也是最容易写出问题的部分。SOC的状态转移方程虽然简单，但要配合充放电状态互斥约束、功率上下限约束一起写，才算完整。
% 储能SOC递推公式 SOC0 = 0.5; % 初始SOC，常见取0.5或论文给定值 % SOC(t+1) = SOC(t) - P_ch*eta_ch*dt/Cap + P_dis*dt/(eta_dis*Cap) Constraints = [Constraints, SOC(2:T+1) == SOC(1:T) ... - P_ch * params.eta_ch * params.dt / params.battery_capacity ... + P_dis * params.dt / (params.eta_dis * params.battery_capacity)]; % SOC上下限 Constraints = [Constraints, params.soc_min <= SOC(2:T+1) <= params.soc_max]; % 充放电状态互斥：同一时刻只能充或只能放 Constraints = [Constraints, P_ch <= params.battery_power_max * u_ch]; Constraints = [Constraints, P_dis <= params.battery_power_max * u_dis]; Constraints = [Constraints, u_ch + u_dis <= 1];互斥约束用了一个0-1变量和一个不等式来实现：u_ch和u_dis不能同时为1，这样就从数学上保证了不可能同时充放电。虽然实际接近工程中也存在“同时充放电效率损失”的情况，但绝大多数学术模型不会这么做。
还有一点是SOC向量的长度。我用的是T+1而不是T，因为SDP变量从1到T+1，其中SOC(1)是初始值，SOC(2)到SOC(T+1)对应每个调度时段的结束状态。不这么设置，递推公式就不方便写成向量化的形式。
4.4 可削减负荷与其他需求响应约束
可削减负荷的约束看起来简单，实际操作时要注意“削减量不得超过合同上限”和“削减量不能为负”两件事，前一个体现对用户舒适度和生产过程的保护，后一个保证优化器不会为了让目标函数更小去“强行增加负荷”。
% 可削减负荷约束：削减比例限值和功率非负 Constraints = [Constraints, 0 <= P_cut <= params.load_cut_max * P_load]; % 全天削减电量上限，防止过度响应 Constraints = [Constraints, sum(P_cut) <= params.total_cut_limit];有的论文还会要求削峰填谷效果不能把峰变成谷，也就是削减后的负荷曲线峰谷差要小于某个阈值，这种约束在Yalmip里表达也非常直接：max(P_load - P_cut) - min(P_load - P_cut) <= 某个值。不过max/min在约束里属于非光滑函数，建议用辅助变量+引入峰谷差上下界的方式处理，或者先不加这个约束，等基础模型跑通之后再补充。
4.5 优化问题的装配与求解
% model/build_optimization_model.m 收尾部分 objective = sum(params.price .* P_grid * params.dt) ... + sum(params.cut_cost * P_cut * params.dt) ... + 1000 * sum(u_ch + u_dis); % 最小化设备动作次数，惩罚项 ops = sdpsettings('solver', 'cplex', 'verbose', 2); optimize(Constraints, objective, ops); % 取出结果，存入结构体 result.P_grid = value(P_grid); result.P_ch = value(P_ch); result.P_dis = value(P_dis); result.SOC = value(SOC); result.P_cut = value(P_cut);那1000乘上0-1变量和的惩罚项是我自己加的。如果不加，系统可能在储能SOC和电价没差异的时段来回切换充放电状态，导致结果里有大量高频动作，这在工程上完全不可接受。加一个足够大的惩罚系数，设备动作次数就会被压缩到合理范围。这个套路在复现论文里很常用，因为很多论文不写这个细节，但不加代码就是会出问题。
5. 求解器选型与数值调优：解不出来、解太慢、解不对怎么办
代码写完只是第一步，真正折磨人的是求解阶段。同一个模型，用不同求解器的表现差异可能非常大。
5.1 混合整数规划为什么是标配，求解器怎么选
需求响应聚合优化里只要涉及设备启停、削负荷决策这样的“有或无”的问题，就一定会有整数变量。模型带整数变量之后，问题类型就从线性规划变成了混合整数线性规划，对求解器的要求高了一个层次。
常见的选择是CPLEX、Gurobi，或者学术免费的SCIP。Yalmip官方文档对每个求解器的支持范围写得很清楚，MILP问题这几个都能解。我个人的经验是，Gurobi在中小规模模型上的求解速度快一些，CPLEX在问题数值稳定性上更省心，SCIP适合预算有限的学生。
ops = sdpsettings('solver', 'gurobi', 'verbose', 2);如果没有商业求解器，也可以用Matlab自带的intlinprog配合Yalmip。不过intlinprog在大规模问题上的表现确实不如专业的商业求解器，如果你的园区设备数量很多、时段分辨率又取到96点，还是装个CPLEX或者Gurobi更稳。
5.2 双线性项怎么处理：三类常见线性化手段
复现过程中最让人头疼的是目标函数或约束里出现变量乘积，比如两个变量相乘：SOC(t)和P_ch(t)乘积表示一种耦合关系，或者某论文里有“负荷削减的0-1状态 × 连续功率”的组合。这类项在数学上叫双线性项，会让问题变成非凸的，没法直接求解。
处理手段按优先级有三种：
第一种，能用逻辑约束表达的就用大M法。一个0-1变量u和连续变量x相乘，如果x的范围已知，就把它替换成辅助变量z，然后用四个约束限定z的取值：z <= Mu、z <= x、z >= x - M(1-u)、z >= 0。这是最普遍的手段。
第二种，如果乘积项是储能充放电状态和SOC相乘，很多时候可以通过约束结构避免乘积——比如限制充电时SOC在某范围内，或者直接用状态转移方程把SOC写成充电功率的函数再化简。
第三种，实在无法线性化，就用Yalmip的非线性建模能力硬算，把问题交给ipopt或者fmincon。但对需求响应这种长时间尺度模型，非线性求解器的稳定性和求解速度都很差，建议只把它当作兜底方案。
5.3 求解时间爆炸和数值异常时的排查清单
模型规模大了以后，最常见的症状是求解器“半天转不出来”，或者提示“ numerical issues（数值问题）”。我按踩过的坑总结了一份排查顺序：
- 检查变量和约束是否出现了数量级差异巨大的系数。比如电价单位是元/千瓦时数值大约在1附近，储能容量却是4000，两个数值放一起，求解器的容差设置会受影响。
- 把所有约束都打印出来检查一遍维度，尤其在矩阵运算中任何一行的维度不匹配Yalmip会直接报错。
- 检查是否存在冗余变量和冗余约束。例如某个0-1变量在目标函数里的权重是0，求解器会耗费大量节点去试探它的取值——这种情况应该直接删掉该变量。
- 如果模型本身没问题但还是慢，可以给求解器设置一个合理的MIP gap（相对最优间隙），比如1%甚至5%。论文复现不追求小数点后五位的精度，一个在可接受范围内的次优解完全够你画出和论文趋势一致的结果图。
6. 结果复现与可视化：把论文里的图还原出来
代码能求解出结果，只完成了一半工作，另一半是把结果图还原成论文里的样式。学术界和工业界看一篇文章，最直观的并不是你的目标函数收敛到多少，而是你的图和原文的趋势是否一致。
6.1 论文图表复刻的基本套路
我复刻论文图表的流程是：先把论文里的图放大仔细看，确认横轴是时刻还是设备编号，纵轴单位是什么，然后拿求解出来的数据用相同坐标范围去画。
以最经典的“优化前后负荷曲线对比图”为例：
% plot/plot_results.m t = 1:24; figure; plot(t, params.load_basic, 'o-', 'LineWidth', 1.2); hold on; plot(t, params.load_basic - result.P_cut, 's-', 'LineWidth', 1.2); stairs(t, result.P_grid, '--', 'LineWidth', 1.2); legend('优化前负荷', '优化后负荷', '电网购电功率'); xlabel('时刻/h'); ylabel('功率/kW'); grid on;如果论文里有储能SOC曲线，我一般用stairs而不是plot，因为SOC在每个时段内是保持不变的，stairs更能反映这个离散时变的属性。很多复现出来的图看起来和论文有差异，不是数据不对，而是绘图方式不对。
6.2 对比实验的设计：怎么体现优化效果
论文里通常会有几组对比实验：基准方案（无需求响应）、只优化储能、储能+削负荷联合优化等。复现时，建议不要在一个模型里来回改配置，而是为每个场景准备一个独立的config文件。
% config_scenario1.m - 无需求响应（基准） params.enable_battery = 0; params.enable_cut = 0; % config_scenario2.m - 仅储能 params.enable_battery = 1; params.enable_cut = 0; % config_scenario3.m - 储能+可削减负荷联合优化 params.enable_battery = 1; params.enable_cut = 1;然后在构建约束时用这些开关控制是否加入对应的约束。这样做的好处是场景之间的切换极其平滑，不会因为注释代码而导致版本混乱。
6.3 让结果经得起追问：保存中间变量
关于结果可视化，还有一条经验：务必把每次求解的P_grid、P_ch、P_dis、SOC、P_cut全部保存到mat文件里，文件名带上场景编号和时间戳。因为复现一个完整的研究往往需要多次迭代，你可能改了一个参数重跑，结果发现新的结果更差了，想对比旧的又找不到原始数据。把每一步中间结果都备份，能让你的复盘省出大量时间。
7. 复现过程中的典型卡点与我的调试清单
最后是复现这类论文最容易踩的五个坑，每一条我都真实遇到过，希望你看完能绕开。
7.1 五大卡点逐条说明
第一个坑是初始SOC设错。论文如果不给初始SOC，很多复现者默认设为0，结果求解器为了满足SOC上下限，被迫一开始就大幅充电，目标函数值偏高。规范做法是初始SOC取中间值0.5，并在论文正文或附录里找依据，找不到就设为0.5并在报告里说明假设。
第二个坑是储能SOC写成了T维而不是T+1维，导致递推公式只能写成for循环，不仅代码冗长而且求解速度慢。用向量化表达，SOC(2:T+1)和SOC(1:T)之间做一个差分，直观又高效。
第三个坑是可削减负荷约束“削减量取整”。部分论文里削负荷是按整数步长进行的（比如只能按100kW为单位切），这就要额外定义整数变量而不是连续变量。如果你忘了加整数约束，复现出来的结果会与论文偏差很大。
第四个坑是充放电效率的位置。充电和放电效率到底是在功率上乘还是在SOC递推里除，不同论文写法不同。建议全程统一用“充电时用电量比充电功率大、放电时供电量比放电功率小”这个物理直觉来检查：SOC(t+1)相比SOC(t)的增量为充电功率×效率，减少量为放电功率÷效率。凡是写成相反的，都会造成储能“凭空多出来电量”的荒唐结果。
第五个坑是单位混乱。我见过有人把光伏出力写成MW，电网购电功率写成kW，功率平衡约束怎么都不满足，最后发现是单位问题。所以我在config里加了一个强制约定：所有与能量、功率、电价相关的参数一律在注释中标注单位，加载数据后先做一轮单位校验，不一致就直接报错。
7.2 调试阶段每天使用的检查清单
调试优化模型时，我每天必做三件事。第一，打印约束数量和变量数量，如果某个场景的约束数量与上一个场景差异不符合预期，一定是约束条件被意外跳过。第二，检查求解状态，Yalmip的optimize返回的problem字段如果非0，必须查清楚是“infeasible（不可行）”还是“unbounded（无界）”。第三，取单时段结果人工验算，把第一个时段的P_grid、P_ch、SOC代入手算功率平衡方程，确认等式成立。
我在实际复现中还有一个习惯：先用一个极小的算例验证模型。比如只取T=4、一台储能设备、一条基础负荷曲线，跑通后再扩大到24时段。这样做能大幅减少排查错误的时间，因为小算例可以手算出预期结果，代码和手算结果对上了，再跑全尺寸算例就有底气。
复现不是终点：把模型改造成“你的版本”
如果你已经完全跑通了这篇论文的代码，接下来我建议你做一件事：把论文里没做但你认为合理的改进加进去。比如给储能寿命加上衰退成本，或者把可削减负荷的补偿模型改成阶梯计价。这不仅是写论文时需要的创新点，也是检验你建模能力的试金石——能改别人的模型，才算真正读懂别人的模型。
我个人体会最深的一点是，复现论文切忌追求“和原文一模一样”。数值完全一致几乎不可能，因为参数取值、求解器设定、初始条件任何一个微小差别都会让结果略有不同，关键是结果的趋势和论文一致，优化曲线在哪个时段抬升、哪个时段削峰，这些和论文对上了，复现就算成功。现在你可以打开Matlab，先跑一个T=24的基础模型，祝顺利。
本文还有配套的精品资源，点击获取
