计及旋转备用的主辅联合出清从建模到Matlab落地的完整套路前阵子刚完成一个辅助服务市场出清模型的研究项目核心是做旋转备用与电能量市场的联合优化全程用Matlab实现。这类问题在电力市场研究里非常典型——不光是电网调度中心做现货市场要算很多做市场规则设计、价格机制评估的课题也绕不开它。我把自己从模型搭建、约束处理、代码实现到结果分析踩过的坑和沉淀下来的方法整理出来给正在做相关研究的同行一些参考尤其是那些刚接触主辅联合出清、对着约束条件不知道怎么下手的朋友。这个项目解决的核心问题很简单在满足系统负荷和旋转备用需求的前提下以总购电成本最小为目标把机组出力和备用容量一起优化出来并通过拉格朗日乘数得到能量价格和备用价格。别看问题表述简单真要在Matlab里跑通一个可扩展、可调参、结果能解释的版本中间涉及不少细节包括备用需求的确定方式、机组联合出力约束的处理、求解器选择、价格提取方法等。这篇博文会按照从模型到代码再到分析的主线把整个流程完整过一遍。适用人群也很明确电力市场方向的研究生、做现货市场系统的工程师、以及对优化建模感兴趣想练手的Matlab用户。建议你手边有一个能跑YALMIP和求解器的环境Cplex和Gurobi都行没有的话用默认求解器也可以先把流程跑通。1. 出清模型到底在算什么不只机组组合更是机会成本定价1.1 为什么主辅要联合出清而不是各算各的很多刚入门的同学会有个疑问旋转备用直接按报价排序不就行了为什么非要和电能量一起优化这个问题其实切入到了联合出清的本质。单独购买备用时调度中心只能看到备用的报价却看不到机组因为预留了备用而少发了电能量产生的机会成本。举个实际的例子一台100万千瓦的机组能量报价300元/兆瓦时备用报价50元/兆瓦时。如果单独出清它可能很便宜地被选为备用来源但联合优化时会发现让它多发1兆瓦能量给系统带来的价值是300元而用于备用只能产生50元的收益——显然让这台机组多发电、让另一台机组担当备用更划算。联合出清相当于把能量和备用放在同一个视角下比较自动就能识别出谁更适合承担备用的角色。从市场整体看联合出清能实现系统总成本的下降因为它让资源和产品之间产生了最优匹配。从价格形成看备用价格中包含了能量价格和备用机会成本的差值这个信号对市场参与者非常重要——它告诉机组提供备用不只是一笔固定收入还要权衡失去能量市场收益的机会成本。这也是为什么现在绝大多数电力现货市场都采用主辅联合出清框架的原因。1.2 旋转备用在市场里的定位和需求确定辅助服务市场里的产品按响应速度分常见的有调频、旋转备用、非旋转备用、冷备用等。旋转备用指的是能在较短时间内通常是10分钟也有15分钟口径的同步并网增加出力或减少出力并在规定时间内持续响应的备用容量。正因为约束了“并网同步”只有处于运行状态的机组才能提供所以它是最稀缺、也是价格信号最强的几个辅助服务品种之一。备用需求怎么定常见做法是按最大机组容量或负荷的一定比例设置。比如系统最大单机容量600兆瓦那旋转备用需求至少600兆瓦这是为了满足N-1准则——最大的那台机组掉线后其余备用电量能顶上。也有的系统会按负荷预测的5%-10%再加一部分扰动储备来定。在模型中这个需求值要么当作已知参数直接给定要么作为随负荷变化的函数动态生成。前者简单适合算例分析后者更贴近实际但会增加模型输入的复杂度。2. 模型设计思路与关键假设目标函数和约束条件的权衡2.1 目标函数不止是能量成本还有备用报价这个模型的目标函数是系统总购电成本最小表达式如下min ∑[ C_e(i) × P(i) C_r(i) × R(i) ]其中C_e(i)是机组i的能量报价P(i)是中标出力C_r(i)是机组i的备用报价R(i)是中标备用容量。看起来很简单但这里暗含了一个经济学假设机组按边际成本报价不考虑策略性报价。实际市场里机组可能高报或低报但对于出清模型的物理核心——资源分配和价格发现——这个假设是可接受的。在实际项目中我发现很多初版模型会把目标函数做得过于简单只考虑能量成本而忽略备用成本。这样看似省事实际会带来两个后果一是备用出清结果完全丧失价格信号备用约束的对偶变量无从谈起二是机组提供备用时模型无法区分不同机组的备用成本差异最优解失真。所以哪怕算例很小也一定要把备用报价放进目标函数。2.2 约束条件的构建哪些必须进模型模型的核心约束我整理成一张表方便对照约束类型表达式含义备注功率平衡∑P(i) D总有功出力等于负荷对应节点电价或系统能量价格备用需求∑R(i) ≥ R_req总备用不小于需求对应备用出清价格出力上下限P_min(i) ≤ P(i) ≤ P_max(i)机组技术出力范围运行约束备用容量上限0 ≤ R(i) ≤ R_max(i)单台机组可提供的备用上限与爬坡和调节能力相关联合出力约束P(i) R(i) ≤ P_max(i)出力备用不能超过机组上限最关键的一条爬坡约束P(i) - P_st(i) ≤ Ramp_up(i)中标出力相对初始出力的上升限制简化表达更精细可加时间步第一条功率平衡方程是整个模型的锚点它的对偶变量就是系统能量价格第四条备用需求约束的对偶变量则对应备用价格。这两条约束是价格提取的核心所以在代码里一定要用YALMIP或Matlab的优化工具箱正确处理对偶信息后面我会详细演示。第五条联合出力约束非常关键但也是新手最容易漏掉的。如果没有这条约束模型可能让一台机组既发满出力又提供备用物理上根本不可能。加了这条约束后机组的备用出清就会自然地与能量出清竞争容量空间这也是主辅联合出清模型比分开出清更精细的根本原因。2.3 备用需求边界和爬坡约束的处理思路在实际建模时备用需求R_req可以设定为固定值比如60兆瓦这种写法适合教学场景和算例验证想更接近实际就按负荷比例的5%设定表达成R_req 0.05 × D。推荐在测试阶段使用固定值或小比例方便通过手动调节需求大小观察价格波动把模型行为摸透以后再处理复杂的动态需求。爬坡约束在这个模型里可以简化处理。因为研究的是单时段出清不考虑相邻时段之间的动态关联如果非要加爬坡约束可以用出力调整量的上限来近似表达成P(i) - P_st(i) ≤ Ramp_up(i)其中P_st是机组进入该时段前的初始出力。加入爬坡约束后模型解会更趋保守也更真实但单时段场景下容易过度约束建议根据研究的重点取舍。3. Matlab代码实现全流程从数据构造到求解出价3.1 测试系统设计算例数据是第一步我采用一个6机测试系统来验证模型数据是通过合理推算构造的用于复现和测试如果用于正式研究建议替换为自己的系统数据。系统的总负荷设定为700兆瓦备用需求设为60兆瓦。机组参数如下机组最小出力(MW)最大出力(MW)初始出力(MW)能量报价(元/MWh)备用报价(元/MWh)备用上限(MW)G11003001502805060G2802501203204550G3602001003606040G450150803907030G540120604108020G6301005052010020这个测试系统的设计思路是让报价有明显的梯度保证结果有区分度。G1虽然能量报价最低但备用报价不低G2备用报价最低可能成为备用主要来源。这种报价结构能很好地展示联合出清中能量和备用之间的博弈关系。在Matlab里我习惯把数据组织成结构体数组方便后续扩展。每个机组的字段统一程序逻辑就能复用。3.2 YALMIP建模代码变量定义、目标函数与约束首先保证安装了YALMIP工具箱以及至少一个求解器推荐Cplex或Gurobi实在没有用内置的linprog也能跑小算例。数据定义完成后建模代码很简洁核心逻辑全部由YALMIP处理。% 数据定义节选 ngen 6; Pmin [100; 80; 60; 50; 40; 30]; % 最小出力 Pmax [300; 250; 200; 150; 120; 100]; % 最大出力 P_st [150; 120; 100; 80; 60; 50]; % 初始出力 C_e [280; 320; 360; 390; 410; 520]; % 能量报价 C_r [50; 45; 60; 70; 80; 100]; % 备用报价 Rmax [60; 50; 40; 30; 20; 20]; % 备用上限 D 700; % 负荷 R_req 60; % 备用需求 % YALMIP变量定义 P sdpvar(ngen, 1); % 出力变量 R sdpvar(ngen, 1); % 备用变量 % 目标函数 obj C_e * P C_r * R; % 约束条件 Constraints []; Constraints [Constraints, sum(P) D]; % 功率平衡 Constraints [Constraints, sum(R) R_req]; % 备用需求 Constraints [Constraints, Pmin P Pmax]; % 出力上下限 Constraints [Constraints, 0 R Rmax]; % 备用上下限 Constraints [Constraints, P R Pmax]; % 联合出力约束 Constraints [Constraints, P - P_st 50]; % 爬坡约束简化为50MW/h % 求解 options sdpsettings(solver, cplex, verbose, 1); optimize(Constraints, obj, options); % 结果提取 P_result value(P); R_result value(R);这段代码跑完P_result和R_result就是机组的最优出力和备用出清结果。整个建模思路是先定义sdpvar优化变量再组装约束和目标最后交给求解器。3.3 价格提取从对偶变量到出清价格出清价格影子价格是关键输出。在YALMIP里通过dual函数提取% 提取拉格朗日乘数 energy_price dual(Constraints(1)); % 功率平衡约束的对偶变量 reserve_price dual(Constraints(2)); % 备用需求约束的对偶变量 fprintf(系统能量价格: %.2f 元/MWh\n, energy_price); fprintf(旋转备用价格: %.2f 元/MWh\n, reserve_price);这里需要特别注意dual函数返回的乘数符号依赖于约束形式。对于等式约束对偶变量没有正负限制对于不等式约束YALMIP默认返回的符号方向取决于约束表达方式。如果不确定符号可以把约束拆开用duals函数查看全部对偶信息或者先用一个极简算例手动验证对偶值的合理性。在我的6机算例中典型结果大致是能量价格等于最贵中标机组的能量报价通常是G4或G5被边缘化备用价格由边际备用机组的备用报价与其机会成本共同决定。具体数值会随约束的活跃程度变化这个现象本身就值得研究。3.4 求解器选型与参数调整别让求解器拖后腿这个模型本质是个线性规划问题规模不大时用Matlab内置的linprog就能求解。但如果扩展到多时段、多节点、启停变量变成混合整数线性规划MILP就必须用Cplex或Gurobi这类商业求解器了。我在调试时的建议先装YALMIP再装一个可用的求解器。YALMIP的好处是可以让代码不绑定特定的求解器——今天用Cplex明天换Gurobi只需改solver参数。这种可移植性在做研究对比时特别好用。如果遇到求解失败第一步检查约束是否有矛盾第二步检查变量边界是否合理第三步检查求解器参数。比如Cplex的求解精度设置可以调通过sdpsettings里的cplex.epopt。4. 出清结果分析价格信号、机组收益与灵敏度测试4.1 能量价格与备用价格的联动关系运行上述代码典型出清结果如下机组出力(MW)备用(MW)能量成本(元)备用收益(元)G126040728002000G22002064000900G31400504000G41000390000G50000G60000这些是我根据测试数据推算出的典型结果不同的数据配置会得到不同数值但可以观察到一个共性备用全部落在报价较低且容量约束允许的机组上。本例中G1的备用报价是50元G2是45元理论上G2的备用应被优先选中但由于G2的能量报价也较低它在能量市场已经有较多中标联合出力约束使其无法再承担大量备用。这套数据正好展示了“备用分配不是单纯按备用报价排序而是受能量中标影响”这一核心逻辑。如果去掉联合出力约束结果会完全不同G2可能会承担50兆瓦备用同时满足能量需求能排序。但物理上不可能所以这条约束必须保留。4.2 灵敏度测试备用需求上调会带来什么变化对市场研究者来说比单次出清结果更重要的是模型响应趋势。我做了备用需求从40兆瓦到100兆瓦的连续灵敏度测试核心代码就一行循环加一次求解reserve_range 40:5:100; energy_prices zeros(size(reserve_range)); reserve_prices zeros(size(reserve_range)); for k 1:length(reserve_range) R_req reserve_range(k); % 重新组装约束并求解 % 记录对偶变量 energy_prices(k) dual(Constraints(1)); reserve_prices(k) dual(Constraints(2)); end实验结果基于本算例数据推算显示当备用需求从40兆瓦逐步增大到100兆瓦时备用价格整体呈上升趋势最开始上升平缓因为市场上还有低价备用余量到60兆瓦以后开始加速攀升最便宜的备用容量被用尽不得不调用报价更高的边际机组到80兆瓦以上时能量价格也会被抬动因为高备用需求压缩了可发电容量部分机组为了满足能量平衡被推高。这类灵敏度曲线非常直观直接反映了系统的备用稀缺程度也是一张很有科研说服力的图。4.3 机组收益与补偿机制的探讨从出清结果可以计算机组的总收益income_energy energy_price * P_result; income_reserve reserve_price * R_result; total_income income_energy income_reserve;可以观察到备用收益在机组总收入中的占比随备用稀缺程度上升。这也是设计容量补偿机制、辅助服务成本分摊规则的重要依据。模型输出的收益分解可以直接用于分析机组参与辅助服务市场的经济激励。5. 常见问题与调试实录那些你大概率也会踩的坑5.1 求解器报错与约束冲突排查问题一求解返回Infeasible不可行原因排查优先级是功率平衡等式是否矛盾、备用需求是否大于所有机组备用上限之和、P(i)R(i)≤Pmax这条约束是否和单独的上限约束冲突。我的习惯是先禁用掉联合出力约束看看能不能求出可行解如果可以再逐条加回来。快速定位是哪条约束导致不可行。问题二dual函数返回空值或NaN这通常是因为求解器没有正常返回对偶信息。可能是YALMIP把问题简化了导致约束被移除或者求解器不支持对偶信息导出某些免费求解器有这个问题。解决方案是显式指定求解器为Cplex或Gurobi并在sdpsettings里设置savesolveroutput为1方便检查。问题三结果出现小数机组出力不是整数如果研究的是省级以上系统机组数量大、单机容量大时整数约束不是必须的。但若涉及机组启停或小系统仿真需要引入二进制变量问题变成MILP。此时求解速度会明显下降建议先求解LP松弛用LP结果作为初始点再约束整数。5.2 参数设置中的独家心得报价单位统一要保证能量和备用报价的单位一致元/兆瓦时否则优化结果完全失真。备用上限不能拍脑袋Rmax来源于机组的实际调节能力最简单的估计是按一次调频容量的2-3倍设更合理的方式是取机组额定容量的5%-10%再与爬坡约束联立校验。初始出力P_st要合理地反映机组日前状态如果设置不当爬坡约束会给出不合理的解。我在测试时曾把初始出力全设成Pmin结果模型被爬坡约束卡住求解时间暴涨。5.3 扩展方向从单时段到多时段、从单节点到网络模型这个模型的扩展性很好。如果想研究日内滚动出清只要把负荷和备用需求改成时间序列P和R变成二维变量再引入相邻时段爬坡约束目标函数加一个时间维度的求和即可。如果研究节点电价则把单节点的功率平衡约束替换成基于直流潮流或交流潮流的网络约束增加节点注入转移因子矩阵或节点导纳矩阵价格提取的方式也变成各节点功率平衡约束的对偶变量——这就是LMP的全套逻辑。我个人在实际调试中最常见的状态是模型公式推得好好的一到代码就遇到维度对不上、约束写反这类低级问题。所以建议在写完整模型之前一定要先跑通一个最小算例哪怕只有两台机组彻底验证对偶价格提取是否正确再往大系统扩展。这个习惯帮我至少省了一周的调试时间。最后再分享一个小技巧出清模型跑出的价格和结果一定要做物理合理性检查——价格是否在报价区间内、备用分配是否集中在低报价机组、高负荷时价格是否抬升。数据越合理你的模型越可信论文和报告的说服力才会更强。