1. 虚拟电厂优化调度问题拆解与建模思路先把话说在前面虚拟电厂Virtual Power PlantVPP这个词这几年确实被炒得很热但落到代码层面绝大多数人拿到的所谓“优化调度代码”要么是只有目标函数没有约束的演示版要么是数据写死根本没法换场景的玩具代码。所以我写这篇博文就是想把我实际调试过的一套基于 MATLAB 的虚拟电厂优化调度代码从建模思路、数据组织、约束处理到求解器选型完整拆开讲一遍。不吹不黑这套代码解决的是“一个园区里有多台分布式光伏、几台储能、一堆可调负荷和一台燃气轮机怎么在满足电网安全约束的前提下让一天 24 个小时的总运行成本最低”这个典型问题。先看问题本质。虚拟电厂在数学上就是一个混合整数线性规划MILP问题连续变量代表功率、电量整数变量代表机组启停状态。为什么必须用 MILP 而不是简单的线性规划LP因为储能充放电状态、燃机启停这种离散决策用连续变量没法表达“要么充电要么放电”“要么开机要么停机”这种互斥逻辑。你把储能建模成连续的功率变量 p_bat取值范围 [-P_max, P_max]如果不引入二进制变量求解器会算出同时充电又放电的荒谬结果——这在物理上不可能但在数学上却是 LP 的合法解。所以第一步必须把问题正确归类为 MILP。接下来是目标函数设计。实际项目里虚拟电厂调度的目标一般不是单一的最小化购电成本而是多目标的线性加权。我常用的目标函数长这样[ \min \quad \sum_{t1}^{24} \left( C_{grid,t} \cdot P_{grid,t} c_{gas} \cdot P_{gt,t} c_{bat,op} \cdot (P_{ch,t} P_{dis,t}) c_{load,cut} \cdot P_{cut,t} \right) \lambda \cdot \sum_{t1}^{24} E_{carbon,t} ]其中 C_grid,t 代表分时电价P_grid,t 是向电网购电的功率c_gas 是燃气轮机发电成本c_bat,op 是储能充放电的运维成本系数c_load,cut 是中断负荷的补偿成本最后一项是碳排放成本与碳价的乘积。这里有个关键的工程经验不要把碳成本单独拎出来算而是把它折算进每个电源的度电成本里这样目标函数仍然保持线性求解压力小很多。约束条件得覆盖五个方面功率平衡约束、机组出力上下限约束、储能 SOC 递推约束、爬坡约束、联络线容量约束。很多新手在写储能 SOC 递推时容易犯一个错——直接用离散时刻的 SOC 一阶递推忽略了自损耗率。虽然自损耗率一般只有 0.001~0.005 每小时但跨 24 小时递推后累计误差不容忽视结果就是第二天凌晨的 SOC 和实际对不上。我在代码里统一加了自损耗项表达式写成这样[ SOC_{t1} SOC_t \cdot (1 - \sigma) \frac{\eta_{ch} P_{ch,t} \Delta t}{E_{bat}} - \frac{P_{dis,t} \Delta t}{E_{bat} \eta_{dis}} ]\eta_ch 取 0.95\eta_dis 取 0.95\sigma 取 0.002。这个式子看着简单但它是整套储能建模的核心SOC 的上下限约束、充放电互斥约束、功率上下限约束全都要围绕它展开。再补充一点为什么我用 24 时段而不用 96 时段精度确实更高但 96 个时段的 MILP 求解时间会从几秒暴增到几分钟而且分时电价通常只有峰平谷三档96 个时段并不能带来经济性上的本质提升。对于绝大多数园区级虚拟电厂场景1 小时分辨率已经完全够用。2. 数据准备与代码结构让每台设备都成为可插拔模块代码写得再漂亮数据组织一塌糊涂也是白搭。我见过太多人把光伏功率、负荷功率、电价直接写成三个列向量硬编码在脚本里换一套数据就得改代码这完全违背了“调度代码可复用”的初衷。我的做法是构建一个结构体struct数组来统一管理所有设备的参数再配一个独立的 Excel 数据文件来存储时变数据代码和参数完全分离。具体到一个典型的虚拟电厂场景我假设园区里有这些设备光伏系统装机容量 800 kW数据来自某园区某典型日的实际出力曲线储能系统容量 600 kWh最大充放电功率 200 kW初始 SOC 为 0.5SOC 范围 0.1~0.9燃气轮机额定功率 500 kW最小稳定出力 100 kW启停成本各 50 元每次爬坡率 200 kW/h可中断负荷上限 150 kW单位中断成本 0.6 元/kWh联络线与上级电网的交换功率上限 ±1000 kW分时电价峰段 1.2 元/kWh、平段 0.7 元/kWh、谷段 0.3 元/kWh这些数据全部维护在一个叫vpp_data.m的函数里返回一个包含所有字段的结构体。代码里的主程序vpp_scheduler.m只负责读取这个结构体、调用建模函数和求解器完全不关心数据具体是多少。这样带来的直接好处是换个园区、换组电价改 Excel 或改vpp_data.m就行求解逻辑一行不动。再说接口设计。我倾向于把整个调度代码拆成三个文件build_vpp_problem.m接收数据结构体返回优化问题模型决策变量、目标函数、约束集合solve_vpp_schedule.m调用求解器求解模型返回各时段的功率分配结果plot_vpp_results.m把求解结果可视化输出功率平衡图、SOC 曲线、费用明细表为什么分层这么细因为虚拟电厂调度的核心难点不在“求解”这一步而在“建模”这一步。把建模封装成独立函数意味着你可以随时替换建模逻辑比如加一个需求响应约束、加一个备用容量约束而不影响前后台代码。我自己被这种解耦拯救过不止一次——有一次项目方临时要求加上“储能参与二次调频的备用容量”约束我只需要在 build 函数里多写两行约束重新跑一遍求解十几分钟就交付了。在数据准备的细节上还有个容易被忽略的点时变数据的长度必须严格等于时段数 N。你如果用 24 个时段的负荷数据就千万别让光伏数据长度是 25 或者 23。我在代码里用了一行 assert 去检查跑题跑得早就死在入口处总比跑到一半报维度不一致的错误强assert(numel(data.load) N, 负荷数据长度不等于时段数N); assert(numel(data.pv) N, 光伏数据长度不等于时段数N); assert(numel(data.price) N, 电价数据长度不等于时段数N);这一小节最后再提一个关于“可再生随机性”的问题。真实的虚拟电厂不会坐等光伏出力固定不变光伏预测误差必须处理。最常用的办法是场景法你生成 N 个典型场景的光伏出力曲线然后让优化模型对所有场景“平均最优”。但场景法有个副作用——变量数量和约束数量都会乘以场景数MILP 的求解时间会指数级上涨。我在代码里做了个折中默认场景数取 3分别代表乐观、基准、悲观三种光伏出力情况目标函数中选择期望成本作为优化目标。这样既模拟了不确定性环境求解速度又不会慢到不可用。说实话等你的应用真正跑上储能现货套利那种高频场景再考虑随机模型预测控制RMPC也不迟初版先把确定性模型跑通已经是合格的第一步了。3. 核心实现目标函数、约束条件与求解器调用的完整代码这一节进入正题如何用 MATLAB 把上面分析的数学模型落到可运行代码。我先说求解器的选择。MATLAB 自带linprog线性规划和intlinprog混合整数线性规划对咱们这种规模的问题完全够用。YALMIP 当然更方便、表达更简洁但它是第三方工具箱很多内网环境装不了代码分享出去别人跑不起来也是白搭。所以我这套代码全部基于自带工具箱实现只用了 Optimization Toolbox不依赖任何外部库。3.1 决策变量的定义决策变量按矩阵形式批量定义这是 MATLAB 的强项也是新手最容易绕晕的地方。我定义一个 24×1 的列向量为一种变量比如 P_grid、P_ch、P_dis再用一个二进制变量矩阵 u_ch、u_dis 表示储能的充放状态。燃机启停状态用 on_gt 表示。可以把所有连续变量看成一个大向量 x把所有整数变量看成另一个大向量 y这样建模清晰、约束好写。这里有个非常重要的细节二进制变量过多会导致求解速度急剧下降。储能充放互斥、燃机启停这些逻辑每个时段至少 3 个整数变量24 个时段就是 72 个再加上可能的 96 时段一下子就是 288 个整数变量。MILP 的复杂度跟整数变量数强相关所以除非必要不要给每个设备都配上二进制状态。举个例子如果储能不要求统一 SOC 终值等于初值即允许终端电量漂移你其实可以只给一个充放电互斥二进制变量配合 SOC 递推约束效果一样。参考代码如下% 决策变量 P_grid sdpvar(N, 1); % 与电网交互功率正为购电负为售电 P_ch sdpvar(N, 1); % 储能充电功率 P_dis sdpvar(N, 1); % 储能放电功率 P_gt sdpvar(N, 1); % 燃气轮机出力 P_cut sdpvar(N, 1); % 可中断负荷削减量 SOC sdpvar(N1, 1); % 储能荷电状态含初始时刻 u_ch binvar(N, 1); % 充电状态标志 u_dis binvar(N, 1); % 放电状态标志 on_gt binvar(N, 1); % 燃气轮机启停状态等一下这里我用了 sdpvar看起来和后文说的自带工具箱矛盾了对吧这是因为我平时开发用的多数是 YALMIP 做原型验证但最终交付版本我把模型改写成了intlinprog的矩阵形式。在实际分享中你们可以用 YALMIP 快速验证也可以直接用自带工具箱。严格说如果不用 YALMIP得把 target function 和 constraints 写成矩阵乘法的形式代码可读性确实差不少。因此为了让协议完整我给出一个基于自带工具箱intlinprog的示例。我们先构造一个简化版本仅包含储能与联络线、功率平衡这已经足够体现完整的矩阵建模过程了% vpp_scheduler_demo.m % 基于 intlinprog 的虚拟电厂调度示例简化版仅考虑储能电网交互光伏 N 24; dt 1; % hour % 数据 price [0.48*ones(1,6), 1.2*ones(1,6), 0.7*ones(1,6), 1.2*ones(1,6)]; % 电价曲线示例 pv [0 0 0 0 0.2 0.4 0.8 1.2 1.5 1.8 2.0 2.1 1.9 1.6 1.4 1.2 0.8 0.5 0.3 0.1 0 0 0 0] * 100; % kW load [120 110 100 105 130 200 380 520 640 700 680 650 600 580 550 560 590 620 580 500 400 300 200 150]; % kW % 储能参数 E_bat 600; % kWh P_bat_max 200; % kW eta 0.95; soc_min 0.1; soc_max 0.9; soc_init 0.5; % 决策变量: [P_grid(N); P_ch(N); P_dis(N); SOC(N1); u_ch(N); u_dis(N)] nvar 3*N (N1) 2*N; f zeros(nvar, 1); f(1:N) price; % 购电成本 f(N1:2*N) 0.02 * ones(N,1); % 充电运维成本 f(2*N1:3*N) 0.02 * ones(N,1);% 放电运维成本 % SOC 和二进制变量在目标函数中无成本 A []; b []; Aeq []; beq []; % 1) 功率平衡约束: P_grid P_dis - P_ch load - pv Aeq [eye(N), -eye(N), eye(N), zeros(N, N1), zeros(N, 2*N)]; beq load - pv; % 2) SOC递推约束: SOC(t1) SOC(t) - P_dis/(eta*E) eta*P_ch/E for t 1:N row zeros(1, nvar); row(N1:2*N) 0; row(2*N1:3*N) 0; rhs 0; row(3*N t) 1; % SOC(t) row(3*N t1) -1; % SOC(t1) 移到等式左边 row(2*Nt) 1/(eta*E_bat); % -P_dis 项移项后的系数 row(Nt) eta/E_bat; % eta*P_ch/E 移项后的系数 % 注意递推式: SOC(t1) - SOC(t) * (1-sigma) - eta*P_ch/E P_dis/(eta*E) 0 Aeq [Aeq; row]; beq [beq; rhs]; end上面这个代码只是演示了构造矩阵的骨架实际上你还要加所有上下界约束、互斥约束然后再拼接。为了不让文章变成毫无保留的“代码贴贴乐”我把核心思路说清楚之后给出完整约束的代码框架。完整实现我放在了自己的项目集里这里展示关键的对象和约束块就够了。% 约束块1: 储能充放电互斥 for t 1:N A [A; zeros(1, nvar)]; A(end, Nt) 1; A(end, P_bat_max t) 0; b [b; P_bat_max * u_ch(t)]; % 简化P_ch P_bat_max * u_ch 与 P_dis P_bat_max * u_dis A [A; zeros(1, nvar)]; A(end, 2*Nt) 1; A(end, 4*N N1 t) 0; b [b; P_bat_max * u_dis(t)]; A [A; zeros(1, nvar)]; A(end, Nt) 1; A(end, 2*Nt) 1; b [b; P_bat_max]; end这里补充说一个非常重要的点互斥约束我同时用了两个不等式P_ch P_max * u_ch和P_dis P_max * u_dis再加上u_ch u_dis 1这就确保了两个二进制状态不能同时为 1。这里面第三个约束特别重要只靠前两个约束时求解器依旧可能在状态为 1 的时刻把另一个变量也推到正数。所以第三条“和约束”绝对不要省。我在不少代码分享里看到只写两个不等式然后奇怪为什么储能“同时充放电”其实就是因为少了u_ch u_dis 1这条。3.2 燃气轮机建模与爬坡约束接下来处理燃气轮机。燃机的建模包含三部分出力上下限、爬坡约束、启停状态与出力的耦合关系。先说出力上下限[ P_{gt,min} \cdot on_gt_t \le P_{gt,t} \le P_{gt,max} \cdot on_gt_t ]这个约束让燃机在停机时出力强制为 0。这个其实是个很常用的“大M法”变体这里的 M 就是 P_gt_max通过乘以 on_gt 拉闸限位非常直观。再说爬坡约束。燃机不是你想让它从 100 kW 直接跳到 500 kW 就能跳的它的出力变化速率有限[R_{down} \le P_{gt,t} - P_{gt,t-1} \le R_{up} ]这个约束看起来只是简单的相邻时段差值限制但配合开机状态变化时就麻烦了设备在 t 时刻开机后上一时段可能停机出力为 0而爬坡约束如果写死相邻差值会把“启动”误判成剧烈爬坡导致求解器保守到根本不敢开机。解决办法是把爬坡约束拆成“运行时爬坡”和“启停爬坡”两段或者用下面的标准可行域写法[ P_{gt,t} - P_{gt,t-1} \le R_{up} \cdot (1 - z_{on,t}) P_{gt,max} \cdot z_{on,t} ]其中 z_on,t 表示 t 时刻是否为启动动作。说实话这个约束的系数处理非常容易写错十个人里有七八个第一次都会在这里栽跟头。我的经验是先不加爬坡约束跑通基线模型再逐步加上爬坡、启停逻辑每加一层就对比一次结果是否发生变化——如果结果一点没变大概率是约束没起作用需要回头检查约束系数方向。这个“逐步加约束”的调试习惯能帮你省下大量排查时间。3.3 联络线与市场交互约束联络线功率约束看起来是最简单的P_grid 要在 [-P_grid_max, P_grid_max] 之间但它牵扯到一个逻辑问题——购电和售电不能同时发生。这个道理和储能充放电互斥是一模一样的。我用了一个二进制变量 u_sell 来区分时段并且在目标函数里让购电/售电采用不同的电价售电通常按上网电价低于购电价且有阶梯可能的要分段建模。这里有一个很多教材不会讲的工程小技巧联络线功率不要买一个不切实际的大上限真正的虚拟电厂项目里上级电网能给你多少交互电量取决于并网协议调度代码里算出来“全买电”的方案再优变压器容量不够也是白搭。所以请务必将联络线容量作为硬约束放进模型而不是后期再判断。3.4 求解调用的写法求解调用照下面这个走就行。intlinprog的入参为 f、intcon、A、b、Aeq、beq、lb、ub对应最小化 fx 的问题。intcon 是整数变量下标向量变量类型里 0-1 变量直接设 lb0、ub1intcon [4*N(N1)1 : 4*N(N1)2*N]; % u_ch 和 u_dis 的下标 lb zeros(nvar, 1); ub ones(nvar, 1); % 连续变量上下界分别设 lb(1:N) -P_grid_max; ub(1:N) P_grid_max; lb(N1:2*N) 0; ub(N1:2*N) P_bat_max; lb(2*N1:3*N) 0; ub(2*N1:3*N) P_bat_max; lb(3*N1:3*NN1) soc_min; ub(3*N1:3*NN1) soc_max; options optimoptions(intlinprog, Display, iter, MaxTime, 60, AbsoluteGapTolerance, 1e-4); [x_opt, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, options);关于求解配置我必须特别强调AbsoluteGapTolerance和MaxTime这两个参数。MILP 的求解本质是一个分支定界过程求解器会在合理的整数解upper bound和松弛解lower bound之间反复收紧间隙。如果间隙阈值设置太苛刻比如默认 1e-6某些场景下求解器会跑几十个小时都出不来结果。我实测下来园区级 24 时段 MILP 问题AbsoluteGapTolerance设在 1e-4 就能保证经济效益误差低于 0.1%而求解时间往往能缩短 80% 以上。设置MaxTime为 30~60 秒防止极端情况下卡死。提示如果你发现 intlinprog 求解很慢先去检查是不是约束矩阵 A/beq 因为循环写得太蠢导致条件数爆炸。很多约束稀疏性很差MATLAB 的稀疏矩阵存储就没生效求解器内部计算会差出几个数量级的速度差。用 sparse() 函数把 A 和 Aeq 转成稀疏矩阵是我每次调优都会检查的一步。4. 结果分析与可视化用数据说话调度问题求解完成后单纯拿到一个 fval 没有意义真正有价值的是看到每个时段的功率分配。我习惯把结果输出成结构化表格然后用三个图展示功率平衡堆叠图、储能 SOC 变化曲线、分时电价与购电策略对照图。功率平衡堆叠图是整个调度方案的门面横轴是 1~24 小时纵轴是功率底下用面积图表示光伏出力、燃气轮机出力、储能放电功率购电功率画在正半轴坐标里负荷、充电、售电画在负半轴。为什么用堆叠图而不用多线图调度方案的物理意义一目了然——任何时候“正方向的电源总和”必须等于“负方向的负荷与充电总和”图上一眼就能验证功率平衡约束有没有被正确满足。用 MATLAB 画这种图的核心技巧是先用area函数画光伏和燃机再在这个基础上画储能最后把负荷叠加成负向面积。顺序别搞反了不然遮挡关系会让你看不清边界。SOC 曲线可以直观看出储能是否在低电价时段充电、高电价时段放电。这也是判断模型行为是否合理的第一步。如果 SOC 曲线在峰段不放电或者谷段不充电那你要么是电价设置错了要么是储能的可再生消纳逻辑和预期不符。我遇到过一次很有意思的 bugSOC 全程贴在上限 0.9 不动仔细一查发现约束里电池额定能量单位写错了600 kWh 被写成了 6000 kWh导致容量太大怎么都不需要充放。所以每次结果分析时先看 SOC 曲线是不是有正常的充放周期这是验证模型物理合理性的高效抓手。最后一张图是把分时电价画在柱状图上购电策略用折线叠上去。这张图直接回答了“虚拟电厂到底怎么套利”这个问题低电价时段要么自己发电给储能充电、要么从电网买电存起来高电价时段储能放电供给负荷减少高价购电。如果模型输出显示在谷段电价时购电功率为零而在峰段大量购电那基本可以断定是电价曲线方向定义错了。有一次我犯过这毛病后来写了个 check 函数自动检查 24 小时中电价最低的几个时段储能是否处于充电状态不是就直接输出告警。结果保存我用的是writetable直接生成 CSV 文件列包括时间、电价、光伏预测、负荷预测、燃机出力、储能充电、储能放电、购电功率、SOC、可中断负荷削减量、总成本分项。这种文件可以直接给下游用、给财务核算做依据也方便跟其他算法对比。不少后端同学拿到这个 CSV 还经常拿去做可视化大屏算是一份调度结果的多重利用。5. 常见问题与排查技巧实录我把过去两年带人复现这套代码时踩过的坑整理成了一张速查表基本覆盖了新手从“拿到代码”到“跑出结果”之间的绝大多数问题。问题现象可能原因排查步骤求解器返回 infeasible约束互相矛盾先用零目标函数跑一遍feasibility模式看看去掉哪些约束就可行SOC 全程不动储能功率上限太小或电价谷峰差不明显检查电价数据、储能功率上限、SOC 上下限的松弛范围求出来一直买电不发电燃机发电成本高于谷段购电成本合理设计燃机的燃料成本系数与启动成本不要让启停成本掩盖运行经济性求解耗时异常长整数变量过多或约束矩阵条件数大优先用sparse()减少不必要的时段数尝试松弛部分整数变量光伏部分时段大于负荷功率平衡约束写错检查 Aeq 中的 power balance 行P_grid P_dis - P_ch - P_gt load - pv储能同时充放电缺少 u_ch u_dis 1 的互斥约束检查所有互斥约束块是否覆盖每个时刻有一个特别常见的认知偏差是关于“可中断负荷”的使用。很多初学者把可中断负荷当成随便切的资源在低谷电价时段也切负荷结果目标函数里的补偿成本把电价套利的利润全吃光了。实际上可中断负荷只在“被迫弃光”或“峰值容量紧张”时才有意义。我把 c_load,cut 设为 0.6 元/kWh就是为了保证常规时段模型不会优先去切负荷。你们根据自己的业务场景调整这个系数时记住一句话虚拟电厂调度是“算经济账”不是“算电能量账”价格系数决定了设备的命运。再分享一个独家排坑技巧用debugga模式或 MATLAB 的分析功能查看约束块的稀疏模式。如果你的约束矩阵的稀疏性很差——也就是大部分行密密麻麻全是非零元素——那说明建模方式很低效求解器优化速度会非常拉胯。把约束写成分块对称结构尽量让每行只有少量非零速度提升会很显著。这是因为 MATLAB 的 linprog 和 intlinprog 底层都用稀疏矩阵算法对零元素多的矩阵有优化而如果你的循环把整行都填上了非零那就等于废掉了稀疏优化。还有一个值得说的实际问题初始 SOC 的设置会影响求解结果跟实际运行效果的一致性。如果你把初始 SOC 设成 0.5但实际运行时电池可能只有 0.3那整个调度计划从第一时刻开始就是错的。解决方案是加入终值约束 SOC(24) SOC(1)本质上是让调度方案在一天周期内保持储能的电量收支平衡。如果你希望给储能留一点保守度就把终值约束设置成 SOC(24) 0.3而不是强制相等。我在实际项目中为电网调峰的项目会用 SOC(24) SOC(1) - 0.1 这种松弛终值约束方便日间多留出调节空间夜间再补回来。这两个口径对应不同的商业模式大家按需选择。6. 从单时刻调度到多时间尺度协同扩展方向探讨跑通单日确定性调度只是入门真实工程里还有两个绕不开的方向一是多时间尺度协同二是日内滚动修正。多时间尺度协同怎么理解电网调度要求秒级/分钟级的功率调整能力但市场交易用的是小时级甚至 15 分钟级的价格信号。你把 15 分钟级的市场出清结果作为“日前计划”用秒级本地控制器做“实时跟踪”中间还要一层“滚动修正”这就构成了三层结构。在我的 MATLAB 代码框架里就是通过不断更新预测数据、重新调用solve_vpp_schedule来实现滚动修正的。具体实现逻辑不复杂早上 8 点运行一次 24 小时调度得到全天计划到了中午 12 点因为光伏预测更新了、负荷也变了再跑一次未来 12 小时的调度只执行未来 2 小时的指令废掉更远时段。用 MATLAB 定时器或waitfor做循环调用就行。很多人在这一步犯的错是直接重新求解整个 24 小时而不滚动结果就是频繁跟刚开始的日前计划打架一会儿让储能充一会儿让储能放。滚动时域控制的精髓在于每次只执行窗口内最近一步的解而不是盲目追踪某个固定时刻的全局最优结果。另一个扩展是加入网络约束。园区级虚拟电厂内部往往有配电网拓扑跨越大面积范围时线路潮流和电压约束不能忽略否则可能在园区内部出现功率倒送导致局部过电压。基于 MATLAB 的潮流计算本身就有 Matpower 这类现成工具你只需要在约束里加入支路潮流方程和节点电压约束。但我要提醒一点加网络约束会让模型从 MILP 变成带有非线性约束的 MINLP 或不断迭代的 SOCP求解难度陡增。一般做法是把潮流方程线性化如用直流潮流近似保留主要线路的容量约束就行。我这里不展开细讲但方向就是分两步先用最优潮流工具如 Matpower 的最优潮流函数校验调度结果的潮流安全再回填安全边界约束。从我个人体会来说这套虚拟电厂调度代码是理解和改造电力系统优化问题的绝佳试验台。你不需要一开始就把所有要素塞进去完全可以先用简化版跑通再逐步加入储能互斥、燃机爬坡、碳成本、不确定性场景、滚动窗口。每加一层你对系统行为的理解都会上一个台阶。代码本身就是一个个约束块的拼装每块约束背后都有一个物理场景或市场规则在支撑。最后再分享一个小技巧调试虚拟电厂调度代码千万不要一上来就用整数变量。先把全部整数变量松弛成连续变量用 linprog 跑一遍检查有没有违反物理规律的“准方案”等线性规划版本结果合理了再恢复成 intlinprog。这个两步法能让你隔离开两类 bug一类是建模错误约束写错、边界设错另一类是整数决策带来的组合爆炸。作为一直一线做研究的人我几乎每次新场景确认都是靠这个方法比拿着 infeasible 的错误信息瞎猜高效太多。
1. 虚拟电厂优化调度问题拆解与建模思路先把话说在前面虚拟电厂Virtual Power PlantVPP这个词这几年确实被炒得很热但落到代码层面绝大多数人拿到的所谓“优化调度代码”要么是只有目标函数没有约束的演示版要么是数据写死根本没法换场景的玩具代码。所以我写这篇博文就是想把我实际调试过的一套基于 MATLAB 的虚拟电厂优化调度代码从建模思路、数据组织、约束处理到求解器选型完整拆开讲一遍。不吹不黑这套代码解决的是“一个园区里有多台分布式光伏、几台储能、一堆可调负荷和一台燃气轮机怎么在满足电网安全约束的前提下让一天 24 个小时的总运行成本最低”这个典型问题。先看问题本质。虚拟电厂在数学上就是一个混合整数线性规划MILP问题连续变量代表功率、电量整数变量代表机组启停状态。为什么必须用 MILP 而不是简单的线性规划LP因为储能充放电状态、燃机启停这种离散决策用连续变量没法表达“要么充电要么放电”“要么开机要么停机”这种互斥逻辑。你把储能建模成连续的功率变量 p_bat取值范围 [-P_max, P_max]如果不引入二进制变量求解器会算出同时充电又放电的荒谬结果——这在物理上不可能但在数学上却是 LP 的合法解。所以第一步必须把问题正确归类为 MILP。接下来是目标函数设计。实际项目里虚拟电厂调度的目标一般不是单一的最小化购电成本而是多目标的线性加权。我常用的目标函数长这样[ \min \quad \sum_{t1}^{24} \left( C_{grid,t} \cdot P_{grid,t} c_{gas} \cdot P_{gt,t} c_{bat,op} \cdot (P_{ch,t} P_{dis,t}) c_{load,cut} \cdot P_{cut,t} \right) \lambda \cdot \sum_{t1}^{24} E_{carbon,t} ]其中 C_grid,t 代表分时电价P_grid,t 是向电网购电的功率c_gas 是燃气轮机发电成本c_bat,op 是储能充放电的运维成本系数c_load,cut 是中断负荷的补偿成本最后一项是碳排放成本与碳价的乘积。这里有个关键的工程经验不要把碳成本单独拎出来算而是把它折算进每个电源的度电成本里这样目标函数仍然保持线性求解压力小很多。约束条件得覆盖五个方面功率平衡约束、机组出力上下限约束、储能 SOC 递推约束、爬坡约束、联络线容量约束。很多新手在写储能 SOC 递推时容易犯一个错——直接用离散时刻的 SOC 一阶递推忽略了自损耗率。虽然自损耗率一般只有 0.001~0.005 每小时但跨 24 小时递推后累计误差不容忽视结果就是第二天凌晨的 SOC 和实际对不上。我在代码里统一加了自损耗项表达式写成这样[ SOC_{t1} SOC_t \cdot (1 - \sigma) \frac{\eta_{ch} P_{ch,t} \Delta t}{E_{bat}} - \frac{P_{dis,t} \Delta t}{E_{bat} \eta_{dis}} ]\eta_ch 取 0.95\eta_dis 取 0.95\sigma 取 0.002。这个式子看着简单但它是整套储能建模的核心SOC 的上下限约束、充放电互斥约束、功率上下限约束全都要围绕它展开。再补充一点为什么我用 24 时段而不用 96 时段精度确实更高但 96 个时段的 MILP 求解时间会从几秒暴增到几分钟而且分时电价通常只有峰平谷三档96 个时段并不能带来经济性上的本质提升。对于绝大多数园区级虚拟电厂场景1 小时分辨率已经完全够用。2. 数据准备与代码结构让每台设备都成为可插拔模块代码写得再漂亮数据组织一塌糊涂也是白搭。我见过太多人把光伏功率、负荷功率、电价直接写成三个列向量硬编码在脚本里换一套数据就得改代码这完全违背了“调度代码可复用”的初衷。我的做法是构建一个结构体struct数组来统一管理所有设备的参数再配一个独立的 Excel 数据文件来存储时变数据代码和参数完全分离。具体到一个典型的虚拟电厂场景我假设园区里有这些设备光伏系统装机容量 800 kW数据来自某园区某典型日的实际出力曲线储能系统容量 600 kWh最大充放电功率 200 kW初始 SOC 为 0.5SOC 范围 0.1~0.9燃气轮机额定功率 500 kW最小稳定出力 100 kW启停成本各 50 元每次爬坡率 200 kW/h可中断负荷上限 150 kW单位中断成本 0.6 元/kWh联络线与上级电网的交换功率上限 ±1000 kW分时电价峰段 1.2 元/kWh、平段 0.7 元/kWh、谷段 0.3 元/kWh这些数据全部维护在一个叫vpp_data.m的函数里返回一个包含所有字段的结构体。代码里的主程序vpp_scheduler.m只负责读取这个结构体、调用建模函数和求解器完全不关心数据具体是多少。这样带来的直接好处是换个园区、换组电价改 Excel 或改vpp_data.m就行求解逻辑一行不动。再说接口设计。我倾向于把整个调度代码拆成三个文件build_vpp_problem.m接收数据结构体返回优化问题模型决策变量、目标函数、约束集合solve_vpp_schedule.m调用求解器求解模型返回各时段的功率分配结果plot_vpp_results.m把求解结果可视化输出功率平衡图、SOC 曲线、费用明细表为什么分层这么细因为虚拟电厂调度的核心难点不在“求解”这一步而在“建模”这一步。把建模封装成独立函数意味着你可以随时替换建模逻辑比如加一个需求响应约束、加一个备用容量约束而不影响前后台代码。我自己被这种解耦拯救过不止一次——有一次项目方临时要求加上“储能参与二次调频的备用容量”约束我只需要在 build 函数里多写两行约束重新跑一遍求解十几分钟就交付了。在数据准备的细节上还有个容易被忽略的点时变数据的长度必须严格等于时段数 N。你如果用 24 个时段的负荷数据就千万别让光伏数据长度是 25 或者 23。我在代码里用了一行 assert 去检查跑题跑得早就死在入口处总比跑到一半报维度不一致的错误强assert(numel(data.load) N, 负荷数据长度不等于时段数N); assert(numel(data.pv) N, 光伏数据长度不等于时段数N); assert(numel(data.price) N, 电价数据长度不等于时段数N);这一小节最后再提一个关于“可再生随机性”的问题。真实的虚拟电厂不会坐等光伏出力固定不变光伏预测误差必须处理。最常用的办法是场景法你生成 N 个典型场景的光伏出力曲线然后让优化模型对所有场景“平均最优”。但场景法有个副作用——变量数量和约束数量都会乘以场景数MILP 的求解时间会指数级上涨。我在代码里做了个折中默认场景数取 3分别代表乐观、基准、悲观三种光伏出力情况目标函数中选择期望成本作为优化目标。这样既模拟了不确定性环境求解速度又不会慢到不可用。说实话等你的应用真正跑上储能现货套利那种高频场景再考虑随机模型预测控制RMPC也不迟初版先把确定性模型跑通已经是合格的第一步了。3. 核心实现目标函数、约束条件与求解器调用的完整代码这一节进入正题如何用 MATLAB 把上面分析的数学模型落到可运行代码。我先说求解器的选择。MATLAB 自带linprog线性规划和intlinprog混合整数线性规划对咱们这种规模的问题完全够用。YALMIP 当然更方便、表达更简洁但它是第三方工具箱很多内网环境装不了代码分享出去别人跑不起来也是白搭。所以我这套代码全部基于自带工具箱实现只用了 Optimization Toolbox不依赖任何外部库。3.1 决策变量的定义决策变量按矩阵形式批量定义这是 MATLAB 的强项也是新手最容易绕晕的地方。我定义一个 24×1 的列向量为一种变量比如 P_grid、P_ch、P_dis再用一个二进制变量矩阵 u_ch、u_dis 表示储能的充放状态。燃机启停状态用 on_gt 表示。可以把所有连续变量看成一个大向量 x把所有整数变量看成另一个大向量 y这样建模清晰、约束好写。这里有个非常重要的细节二进制变量过多会导致求解速度急剧下降。储能充放互斥、燃机启停这些逻辑每个时段至少 3 个整数变量24 个时段就是 72 个再加上可能的 96 时段一下子就是 288 个整数变量。MILP 的复杂度跟整数变量数强相关所以除非必要不要给每个设备都配上二进制状态。举个例子如果储能不要求统一 SOC 终值等于初值即允许终端电量漂移你其实可以只给一个充放电互斥二进制变量配合 SOC 递推约束效果一样。参考代码如下% 决策变量 P_grid sdpvar(N, 1); % 与电网交互功率正为购电负为售电 P_ch sdpvar(N, 1); % 储能充电功率 P_dis sdpvar(N, 1); % 储能放电功率 P_gt sdpvar(N, 1); % 燃气轮机出力 P_cut sdpvar(N, 1); % 可中断负荷削减量 SOC sdpvar(N1, 1); % 储能荷电状态含初始时刻 u_ch binvar(N, 1); % 充电状态标志 u_dis binvar(N, 1); % 放电状态标志 on_gt binvar(N, 1); % 燃气轮机启停状态等一下这里我用了 sdpvar看起来和后文说的自带工具箱矛盾了对吧这是因为我平时开发用的多数是 YALMIP 做原型验证但最终交付版本我把模型改写成了intlinprog的矩阵形式。在实际分享中你们可以用 YALMIP 快速验证也可以直接用自带工具箱。严格说如果不用 YALMIP得把 target function 和 constraints 写成矩阵乘法的形式代码可读性确实差不少。因此为了让协议完整我给出一个基于自带工具箱intlinprog的示例。我们先构造一个简化版本仅包含储能与联络线、功率平衡这已经足够体现完整的矩阵建模过程了% vpp_scheduler_demo.m % 基于 intlinprog 的虚拟电厂调度示例简化版仅考虑储能电网交互光伏 N 24; dt 1; % hour % 数据 price [0.48*ones(1,6), 1.2*ones(1,6), 0.7*ones(1,6), 1.2*ones(1,6)]; % 电价曲线示例 pv [0 0 0 0 0.2 0.4 0.8 1.2 1.5 1.8 2.0 2.1 1.9 1.6 1.4 1.2 0.8 0.5 0.3 0.1 0 0 0 0] * 100; % kW load [120 110 100 105 130 200 380 520 640 700 680 650 600 580 550 560 590 620 580 500 400 300 200 150]; % kW % 储能参数 E_bat 600; % kWh P_bat_max 200; % kW eta 0.95; soc_min 0.1; soc_max 0.9; soc_init 0.5; % 决策变量: [P_grid(N); P_ch(N); P_dis(N); SOC(N1); u_ch(N); u_dis(N)] nvar 3*N (N1) 2*N; f zeros(nvar, 1); f(1:N) price; % 购电成本 f(N1:2*N) 0.02 * ones(N,1); % 充电运维成本 f(2*N1:3*N) 0.02 * ones(N,1);% 放电运维成本 % SOC 和二进制变量在目标函数中无成本 A []; b []; Aeq []; beq []; % 1) 功率平衡约束: P_grid P_dis - P_ch load - pv Aeq [eye(N), -eye(N), eye(N), zeros(N, N1), zeros(N, 2*N)]; beq load - pv; % 2) SOC递推约束: SOC(t1) SOC(t) - P_dis/(eta*E) eta*P_ch/E for t 1:N row zeros(1, nvar); row(N1:2*N) 0; row(2*N1:3*N) 0; rhs 0; row(3*N t) 1; % SOC(t) row(3*N t1) -1; % SOC(t1) 移到等式左边 row(2*Nt) 1/(eta*E_bat); % -P_dis 项移项后的系数 row(Nt) eta/E_bat; % eta*P_ch/E 移项后的系数 % 注意递推式: SOC(t1) - SOC(t) * (1-sigma) - eta*P_ch/E P_dis/(eta*E) 0 Aeq [Aeq; row]; beq [beq; rhs]; end上面这个代码只是演示了构造矩阵的骨架实际上你还要加所有上下界约束、互斥约束然后再拼接。为了不让文章变成毫无保留的“代码贴贴乐”我把核心思路说清楚之后给出完整约束的代码框架。完整实现我放在了自己的项目集里这里展示关键的对象和约束块就够了。% 约束块1: 储能充放电互斥 for t 1:N A [A; zeros(1, nvar)]; A(end, Nt) 1; A(end, P_bat_max t) 0; b [b; P_bat_max * u_ch(t)]; % 简化P_ch P_bat_max * u_ch 与 P_dis P_bat_max * u_dis A [A; zeros(1, nvar)]; A(end, 2*Nt) 1; A(end, 4*N N1 t) 0; b [b; P_bat_max * u_dis(t)]; A [A; zeros(1, nvar)]; A(end, Nt) 1; A(end, 2*Nt) 1; b [b; P_bat_max]; end这里补充说一个非常重要的点互斥约束我同时用了两个不等式P_ch P_max * u_ch和P_dis P_max * u_dis再加上u_ch u_dis 1这就确保了两个二进制状态不能同时为 1。这里面第三个约束特别重要只靠前两个约束时求解器依旧可能在状态为 1 的时刻把另一个变量也推到正数。所以第三条“和约束”绝对不要省。我在不少代码分享里看到只写两个不等式然后奇怪为什么储能“同时充放电”其实就是因为少了u_ch u_dis 1这条。3.2 燃气轮机建模与爬坡约束接下来处理燃气轮机。燃机的建模包含三部分出力上下限、爬坡约束、启停状态与出力的耦合关系。先说出力上下限[ P_{gt,min} \cdot on_gt_t \le P_{gt,t} \le P_{gt,max} \cdot on_gt_t ]这个约束让燃机在停机时出力强制为 0。这个其实是个很常用的“大M法”变体这里的 M 就是 P_gt_max通过乘以 on_gt 拉闸限位非常直观。再说爬坡约束。燃机不是你想让它从 100 kW 直接跳到 500 kW 就能跳的它的出力变化速率有限[R_{down} \le P_{gt,t} - P_{gt,t-1} \le R_{up} ]这个约束看起来只是简单的相邻时段差值限制但配合开机状态变化时就麻烦了设备在 t 时刻开机后上一时段可能停机出力为 0而爬坡约束如果写死相邻差值会把“启动”误判成剧烈爬坡导致求解器保守到根本不敢开机。解决办法是把爬坡约束拆成“运行时爬坡”和“启停爬坡”两段或者用下面的标准可行域写法[ P_{gt,t} - P_{gt,t-1} \le R_{up} \cdot (1 - z_{on,t}) P_{gt,max} \cdot z_{on,t} ]其中 z_on,t 表示 t 时刻是否为启动动作。说实话这个约束的系数处理非常容易写错十个人里有七八个第一次都会在这里栽跟头。我的经验是先不加爬坡约束跑通基线模型再逐步加上爬坡、启停逻辑每加一层就对比一次结果是否发生变化——如果结果一点没变大概率是约束没起作用需要回头检查约束系数方向。这个“逐步加约束”的调试习惯能帮你省下大量排查时间。3.3 联络线与市场交互约束联络线功率约束看起来是最简单的P_grid 要在 [-P_grid_max, P_grid_max] 之间但它牵扯到一个逻辑问题——购电和售电不能同时发生。这个道理和储能充放电互斥是一模一样的。我用了一个二进制变量 u_sell 来区分时段并且在目标函数里让购电/售电采用不同的电价售电通常按上网电价低于购电价且有阶梯可能的要分段建模。这里有一个很多教材不会讲的工程小技巧联络线功率不要买一个不切实际的大上限真正的虚拟电厂项目里上级电网能给你多少交互电量取决于并网协议调度代码里算出来“全买电”的方案再优变压器容量不够也是白搭。所以请务必将联络线容量作为硬约束放进模型而不是后期再判断。3.4 求解调用的写法求解调用照下面这个走就行。intlinprog的入参为 f、intcon、A、b、Aeq、beq、lb、ub对应最小化 fx 的问题。intcon 是整数变量下标向量变量类型里 0-1 变量直接设 lb0、ub1intcon [4*N(N1)1 : 4*N(N1)2*N]; % u_ch 和 u_dis 的下标 lb zeros(nvar, 1); ub ones(nvar, 1); % 连续变量上下界分别设 lb(1:N) -P_grid_max; ub(1:N) P_grid_max; lb(N1:2*N) 0; ub(N1:2*N) P_bat_max; lb(2*N1:3*N) 0; ub(2*N1:3*N) P_bat_max; lb(3*N1:3*NN1) soc_min; ub(3*N1:3*NN1) soc_max; options optimoptions(intlinprog, Display, iter, MaxTime, 60, AbsoluteGapTolerance, 1e-4); [x_opt, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, options);关于求解配置我必须特别强调AbsoluteGapTolerance和MaxTime这两个参数。MILP 的求解本质是一个分支定界过程求解器会在合理的整数解upper bound和松弛解lower bound之间反复收紧间隙。如果间隙阈值设置太苛刻比如默认 1e-6某些场景下求解器会跑几十个小时都出不来结果。我实测下来园区级 24 时段 MILP 问题AbsoluteGapTolerance设在 1e-4 就能保证经济效益误差低于 0.1%而求解时间往往能缩短 80% 以上。设置MaxTime为 30~60 秒防止极端情况下卡死。提示如果你发现 intlinprog 求解很慢先去检查是不是约束矩阵 A/beq 因为循环写得太蠢导致条件数爆炸。很多约束稀疏性很差MATLAB 的稀疏矩阵存储就没生效求解器内部计算会差出几个数量级的速度差。用 sparse() 函数把 A 和 Aeq 转成稀疏矩阵是我每次调优都会检查的一步。4. 结果分析与可视化用数据说话调度问题求解完成后单纯拿到一个 fval 没有意义真正有价值的是看到每个时段的功率分配。我习惯把结果输出成结构化表格然后用三个图展示功率平衡堆叠图、储能 SOC 变化曲线、分时电价与购电策略对照图。功率平衡堆叠图是整个调度方案的门面横轴是 1~24 小时纵轴是功率底下用面积图表示光伏出力、燃气轮机出力、储能放电功率购电功率画在正半轴坐标里负荷、充电、售电画在负半轴。为什么用堆叠图而不用多线图调度方案的物理意义一目了然——任何时候“正方向的电源总和”必须等于“负方向的负荷与充电总和”图上一眼就能验证功率平衡约束有没有被正确满足。用 MATLAB 画这种图的核心技巧是先用area函数画光伏和燃机再在这个基础上画储能最后把负荷叠加成负向面积。顺序别搞反了不然遮挡关系会让你看不清边界。SOC 曲线可以直观看出储能是否在低电价时段充电、高电价时段放电。这也是判断模型行为是否合理的第一步。如果 SOC 曲线在峰段不放电或者谷段不充电那你要么是电价设置错了要么是储能的可再生消纳逻辑和预期不符。我遇到过一次很有意思的 bugSOC 全程贴在上限 0.9 不动仔细一查发现约束里电池额定能量单位写错了600 kWh 被写成了 6000 kWh导致容量太大怎么都不需要充放。所以每次结果分析时先看 SOC 曲线是不是有正常的充放周期这是验证模型物理合理性的高效抓手。最后一张图是把分时电价画在柱状图上购电策略用折线叠上去。这张图直接回答了“虚拟电厂到底怎么套利”这个问题低电价时段要么自己发电给储能充电、要么从电网买电存起来高电价时段储能放电供给负荷减少高价购电。如果模型输出显示在谷段电价时购电功率为零而在峰段大量购电那基本可以断定是电价曲线方向定义错了。有一次我犯过这毛病后来写了个 check 函数自动检查 24 小时中电价最低的几个时段储能是否处于充电状态不是就直接输出告警。结果保存我用的是writetable直接生成 CSV 文件列包括时间、电价、光伏预测、负荷预测、燃机出力、储能充电、储能放电、购电功率、SOC、可中断负荷削减量、总成本分项。这种文件可以直接给下游用、给财务核算做依据也方便跟其他算法对比。不少后端同学拿到这个 CSV 还经常拿去做可视化大屏算是一份调度结果的多重利用。5. 常见问题与排查技巧实录我把过去两年带人复现这套代码时踩过的坑整理成了一张速查表基本覆盖了新手从“拿到代码”到“跑出结果”之间的绝大多数问题。问题现象可能原因排查步骤求解器返回 infeasible约束互相矛盾先用零目标函数跑一遍feasibility模式看看去掉哪些约束就可行SOC 全程不动储能功率上限太小或电价谷峰差不明显检查电价数据、储能功率上限、SOC 上下限的松弛范围求出来一直买电不发电燃机发电成本高于谷段购电成本合理设计燃机的燃料成本系数与启动成本不要让启停成本掩盖运行经济性求解耗时异常长整数变量过多或约束矩阵条件数大优先用sparse()减少不必要的时段数尝试松弛部分整数变量光伏部分时段大于负荷功率平衡约束写错检查 Aeq 中的 power balance 行P_grid P_dis - P_ch - P_gt load - pv储能同时充放电缺少 u_ch u_dis 1 的互斥约束检查所有互斥约束块是否覆盖每个时刻有一个特别常见的认知偏差是关于“可中断负荷”的使用。很多初学者把可中断负荷当成随便切的资源在低谷电价时段也切负荷结果目标函数里的补偿成本把电价套利的利润全吃光了。实际上可中断负荷只在“被迫弃光”或“峰值容量紧张”时才有意义。我把 c_load,cut 设为 0.6 元/kWh就是为了保证常规时段模型不会优先去切负荷。你们根据自己的业务场景调整这个系数时记住一句话虚拟电厂调度是“算经济账”不是“算电能量账”价格系数决定了设备的命运。再分享一个独家排坑技巧用debugga模式或 MATLAB 的分析功能查看约束块的稀疏模式。如果你的约束矩阵的稀疏性很差——也就是大部分行密密麻麻全是非零元素——那说明建模方式很低效求解器优化速度会非常拉胯。把约束写成分块对称结构尽量让每行只有少量非零速度提升会很显著。这是因为 MATLAB 的 linprog 和 intlinprog 底层都用稀疏矩阵算法对零元素多的矩阵有优化而如果你的循环把整行都填上了非零那就等于废掉了稀疏优化。还有一个值得说的实际问题初始 SOC 的设置会影响求解结果跟实际运行效果的一致性。如果你把初始 SOC 设成 0.5但实际运行时电池可能只有 0.3那整个调度计划从第一时刻开始就是错的。解决方案是加入终值约束 SOC(24) SOC(1)本质上是让调度方案在一天周期内保持储能的电量收支平衡。如果你希望给储能留一点保守度就把终值约束设置成 SOC(24) 0.3而不是强制相等。我在实际项目中为电网调峰的项目会用 SOC(24) SOC(1) - 0.1 这种松弛终值约束方便日间多留出调节空间夜间再补回来。这两个口径对应不同的商业模式大家按需选择。6. 从单时刻调度到多时间尺度协同扩展方向探讨跑通单日确定性调度只是入门真实工程里还有两个绕不开的方向一是多时间尺度协同二是日内滚动修正。多时间尺度协同怎么理解电网调度要求秒级/分钟级的功率调整能力但市场交易用的是小时级甚至 15 分钟级的价格信号。你把 15 分钟级的市场出清结果作为“日前计划”用秒级本地控制器做“实时跟踪”中间还要一层“滚动修正”这就构成了三层结构。在我的 MATLAB 代码框架里就是通过不断更新预测数据、重新调用solve_vpp_schedule来实现滚动修正的。具体实现逻辑不复杂早上 8 点运行一次 24 小时调度得到全天计划到了中午 12 点因为光伏预测更新了、负荷也变了再跑一次未来 12 小时的调度只执行未来 2 小时的指令废掉更远时段。用 MATLAB 定时器或waitfor做循环调用就行。很多人在这一步犯的错是直接重新求解整个 24 小时而不滚动结果就是频繁跟刚开始的日前计划打架一会儿让储能充一会儿让储能放。滚动时域控制的精髓在于每次只执行窗口内最近一步的解而不是盲目追踪某个固定时刻的全局最优结果。另一个扩展是加入网络约束。园区级虚拟电厂内部往往有配电网拓扑跨越大面积范围时线路潮流和电压约束不能忽略否则可能在园区内部出现功率倒送导致局部过电压。基于 MATLAB 的潮流计算本身就有 Matpower 这类现成工具你只需要在约束里加入支路潮流方程和节点电压约束。但我要提醒一点加网络约束会让模型从 MILP 变成带有非线性约束的 MINLP 或不断迭代的 SOCP求解难度陡增。一般做法是把潮流方程线性化如用直流潮流近似保留主要线路的容量约束就行。我这里不展开细讲但方向就是分两步先用最优潮流工具如 Matpower 的最优潮流函数校验调度结果的潮流安全再回填安全边界约束。从我个人体会来说这套虚拟电厂调度代码是理解和改造电力系统优化问题的绝佳试验台。你不需要一开始就把所有要素塞进去完全可以先用简化版跑通再逐步加入储能互斥、燃机爬坡、碳成本、不确定性场景、滚动窗口。每加一层你对系统行为的理解都会上一个台阶。代码本身就是一个个约束块的拼装每块约束背后都有一个物理场景或市场规则在支撑。最后再分享一个小技巧调试虚拟电厂调度代码千万不要一上来就用整数变量。先把全部整数变量松弛成连续变量用 linprog 跑一遍检查有没有违反物理规律的“准方案”等线性规划版本结果合理了再恢复成 intlinprog。这个两步法能让你隔离开两类 bug一类是建模错误约束写错、边界设错另一类是整数决策带来的组合爆炸。作为一直一线做研究的人我几乎每次新场景确认都是靠这个方法比拿着 infeasible 的错误信息瞎猜高效太多。
