电力市场的出清计算说白了就是在一个巨大的经济调度问题里找平衡点。这几年我做过不少相关的优化项目最头疼的往往不是连续量的经济调度而是那些带 0-1 整数变量的均衡问题。尤其当场景里再加上补偿费用、机组上升爬坡约束问题立刻从普通的 LP 变成混合整数均衡问题常规求解器直接罢工手写对偶条件又到处都是坑。今天就把这一类“带补偿和电力市场上升问题的二元平衡问题精确求解”的方法从头到尾拆开讲一遍。后面所有内容都基于我实际跑过的项目经验模型和参数都尽量给全方便你直接参考。这类问题表面上是数学规划实际上是一个包含“市场出清”和“机组决策”两层关系的博弈。所谓的“二元”一方面是指模型里有离散决策变量比如机组启停0 或者 1另一方面也指这个平衡发生在“价格机制”和“物理约束”两个维度之间。加上补偿项以后目标函数不再是简单的成本最小化而要考虑额外的激励费用加上“上升问题”本质上是机组向上爬坡速率限制与时间耦合。几件事叠在一起精确求解的难度不是加法而是乘法。1. 先把这个“二元平衡问题”掰开揉碎1.1 从电力市场出清说起为什么会有 0-1 变量和连续变量共存很多刚接触电力市场建模的工程师第一反应是市场出清不就是线性规划吗把负荷、机组报价放进去目标函数最小化总购电成本约束满足功率平衡和潮流完事了。没错如果所有机组都已经确定开机或停机并且不考虑启停成本那确实是一个标准的连续线性规划问题。但真实市场不是这样。机组在日前市场中必须决定明天的开机状态这个状态是离散的开就是 1关就是 0。而开机的机组还要决定每时段的出力这是连续变量。于是整个问题出现了两类变量混扫的情况0-1 变量决定结构性选择连续变量决定运行点选择。这本身就是一个混合整数二次规划或混合整数线性规划的雏形。更麻烦的是在均衡视角下市场出清方系统运营商和发电商之间存在策略互动。发电商知道自己会影响出清价格因此会在报价里加入策略性抬价或压价而系统运营商必须按既定的市场规则解出清。这个博弈的解就是一个均衡点。因为策略空间里包含离散状态这就不再是普通纳什均衡而是“二元均衡”或“混合整数均衡”。1.2 “补偿”在模型里到底补什么“带补偿”这三个字说起来容易建模时却要非常小心。补偿从哪儿来最常见的场景有两类第一类是机组因提供调频、备用等辅助服务而获得的容量补偿第二类是系统为了激励机组在低电价时段不关机、高电价时段快速响应给出的“可用性补偿”。从数学上讲补偿通常是一个与机组启停状态或承诺容量有关的线性或二次项。比如一台机组如果被要求必须保持开机状态那么即使它实际出力很低也应该获得一笔固定补偿用来覆盖其空载成本或机会成本。这笔钱在目标函数里表现为补偿系数乘以一个 0-1 变量。补偿项最微妙的地方在于它改变了目标函数的凸性结构。本来单纯购电成本最小化时目标函数是凸的加上一个与整数变量相乘的补偿项后目标函数不再保证凸性。这种情况下直接用 KKT 条件全局求解就会失效。这也是为什么需要“精确求解方法”——也就是能够在有限步内证明找到的解确实是最优均衡解的方法而不是靠运气找一个局部均衡。1.3 “上升问题”不是小问题向上爬坡约束的建模细节标题里的“上升问题”我猜不少同行也会脑门一亮。不是疫情的那种上升也不是股票指数上升而是机组出力在一个时间段内向上调节的能力限制。专业上叫“爬坡约束”或“ramping constraint”。它描述的是t 时刻的出力不能比 t-1 时刻高出超过某一个速率。建模的时候很多人会把它简单写成P_{t} - P_{t-1} RU但加上启停变量以后问题就来了。如果一台机组在 t-1 时刻是关机的P_{t-1}0t 时刻开机那它必须在一个小时内从零升到额定出力甚至更高这显然不合理。所以实际的模型必须区分“运行状态下爬坡”和“启动过程”两类约束。更细一点说还要引入启动变量和停机变量。通常用三个 0-1 变量状态变量 u_t1 表示运行、启动变量 v_t1 表示本时段刚开机、停机变量 w_t1 表示本时段刚停机。于是上升约束往往写成P_t - P_{t-1} RU * u_{t-1} SU * v_t这里的 SU 是启动速率RU 是正常运行时的向上爬坡速率。两个速率可能差好几倍。这个约束虽然只是一个小式子但它把从前时段的运行状态和历史出力全部耦合在一起导致整个时间维度的变量之间形成带状稀疏结构。求解时矩阵性质变得非常复杂也让局限的启发式算法很容易顾此失彼。1.4 平衡问题在数学上是什么从单层优化到均衡约束先别笑这个名字拗口。所谓“平衡问题”在数学优化里往往指代一类变分不等式或互补问题。最典型的形式是0 x ⊥ F(x) 0意思是 x 满足非负性F(x) 满足非负性且两者不能同时严格大于零至少有一个是零。这个互补关系在电力市场里到处都是。比如节点电价等于发电出力与报价之间的边际成本互补约束。把市场出清的 KKT 条件与发电商的策略性问题组合到一起就形成了均衡约束数学规划MPEC 或 EPEC。当其中含有离散变量时问题从“非凸连续问题”进一步恶化成“混合整数非凸问题”。这也就是我们标题里所谓的“二元平衡问题”。很多初学者会误以为这只是一个更大的 MILP实际上它的可行域可能是非凸、非连通、甚至包含孤立的可行点。精确求解这类问题的核心就是把非凸结构转化为可处理的混合整数线性或二次规划。2. 精确求解的核心难点和思路拆解2.1 混合整数带来的组合爆炸0-1 变量带来的最大问题是组合爆炸。假设系统里有 100 台机组、24 个时段那么状态变量 u_{i,t} 就有 2400 个。如果直接枚举所有启停组合那规模是 2^{2400} 量级宇宙寿命耗尽也算不完。所以必须依赖分支定界、割平面这些算法。但更本质的难点在于平衡约束的加入让原本的 MILP 松弛结构变差。如果你把整数变量松弛成 [0,1] 连续变量KKT 条件可能会出现劣质解或奇异点。精确求解算法必须保证在整数变量取离散值的同时连续子问题达到全局最优。这就涉及到求解器内部的“强分支”和“割平面循环”不是一股脑扔给求解器就行。2.2 平衡约束为什么不能用常规求解器直接解常规求解器比如 Gurobi、Cplex、COPT都能解 MILP 和 MIQP但它们默认的目标函数和约束都是显式的。而平衡约束里有一个极其讨厌的互补条件比如0 λ ⊥ (a x - b) 0这里 λ 是乘子a x - b 是某个不等式约束。Gurobi 并不直接支持这种描述。你需要把它转化为一组“大 M 线性约束”λ 0a x - b 0λ M za x - b M(1-z)其中 z 是新的 0-1 变量。这个方法说起来简单做起来全是细节。M 取值太小可能截断可行域导致误判无解M 取值太大求解器数值稳定性崩溃割平面收敛极慢。我后面会详细讲怎么调 M这里先提个醒大 M 的选择本身就是一门手艺。2.3 精确求解的意义不是“差不多”而是要证明全局最优有人会问启发式算法比如遗传算法、粒子群也能找到不错的解为什么非要精确求解因为电力市场出清是涉及真金白银的系统运行决策。你给调度员一个次优解可能意味着多付了几十万的购电成本更严重的是可能违反安全约束。更重要的是在博弈场景下启发式算法无法提供“均衡性证明”——你没法说明为什么对方没有动机偏离这个策略。精确求解的意义在于它能在有限步内给出一个带有最优性间隙上界的解并且当间隙为 0 时你可以用数学上严格的逻辑证明这个解是全局最优的。在市场结算、投资规划等场景下这个证明是必要的合规要求。2.4 可行的算法路线概览大M法、强对偶/KKT、分段线性化、割平面法把这么多难点放在一起总得给几条能落地的路。我实际验证过四条路线各有适用场景这里先总结一下。大 M 法把互补约束展开为混合整数线性约束加入二进制辅助变量。适用于小规模或中等规模问题思路简单调试方便但对 M 敏感。强对偶转换利用线性规划强对偶定理把下层出清问题的最优性条件替换为对偶可行性和强对偶等式从而把双层问题转成单层约束。这种方式在连续下层问题里非常好用但要求下层问题没有整数变量。分段线性化把非线性补偿项或二次成本项在离散点上线性化配合 SOS2 约束或增量线性化方法得到精确的 MILP 近似。注意是近似不是精确需要控制分段数量以平衡精度与速度。Benders 分解/割平面法把整数主问题和连续子问题解耦通过子问题的线性化拉格朗日乘子生成割平面。这种方式适合大规模机组组合问题但实现复杂度高对问题结构依赖强。3. 实操过程一个带补偿的机组组合-市场出清联合模型求解示例3.1 模型假设与参数设计我设计一个简化但不失真实性的算例给你看看完整的建模和求解过程。假设系统有 3 台机组4 个时段负荷曲线为 [400, 550, 650, 500] MW。机组参数如下机组最大出力(MW)最小出力(MW)启动成本(元)空载成本(元/h)边际成本(元/MWh)向上爬坡速率(MW/h)G1300605005080120G22504040040100100G3200303003012080补偿机制设为如果系统要求 G1 在任何时段保持开机但出力低于 120 MW则每个低出力时段支付 200 元补偿。这个补偿会进入目标函数你会发现它直接改变了 G1 的开机决策。3.2 写出完整的 MICP/MIQCP 模型我们建立如下模型。设 i 为机组集合t 为时段。决策变量u_{i,t}0-1表示机组 i 在时段 t 的运行状态v_{i,t}0-1表示机组 i 在时段 t 是否启动w_{i,t}0-1表示机组 i 在时段 t 是否停机p_{i,t}连续变量机组 i 在时段 t 的出力s_{i,t}连续变量表示补偿状态对应的出力或者我们可以用一个辅助二元变量 c_{i,t} 表示是否触发补偿。为方便这里直接引入补偿项为线性项如果 u_{i,t}1 且 p_{i,t} 120则支付 200 元。这个条件可以用二进制辅助变量 y_{i,t} 表示并为 1 表示满足补偿条件。约束条件功率平衡约束每个时段所有机组出力之和等于负荷。机组出力上下限u_{i,t} * P_i_min p_{i,t} u_{i,t} * P_i_max。启动/停机逻辑u_{i,t} - u_{i,t-1} v_{i,t} - w_{i,t}且 v_{i,t} w_{i,t} 1。向上爬坡约束p_{i,t} - p_{i,t-1} RU_i * u_{i,t-1} SU_i * v_{i,t}。启动速率 SU_i 这里简化取 200 MW/h大于爬坡速率表示启动过程更灵活但也不能瞬间满发。补偿触发约束p_{i,t} 120 M8 * (1 - y_{i,t})p_{i,t} 120 - M9 * (1 - y_{i,t})其实不需要大于等于。我们只需要在 p120 时让 y1。可以用 p_{i,t} - 120 M_y * (1 - y_{i,t})y_{i,t} u_{i,t}y_{i,t} 为 0-1。 同时目标函数中加上 200 * y_{i,t}注意这里是补偿成本计入最小化。或者换一种更简单的补偿建模不引入 y直接令补偿成本为 200 * u_{i,t} * (p_{i,t} 120)但这是非线性不可导的需要整数变量所以还是用 y 好。目标函数 min Σ_{i,t} (启停成本 空载成本 边际发电成本 补偿成本)其中启动成本 SU_i * v_{i,t}停机成本可以设为零或一个值空载成本 50 * u_{i,t}边际发电成本 80 * p_{i,t}以 G1 为例补偿成本 200 * y_{i,t}。这里要特别说明如果补偿项是支付给机组那么市场运营者的总成本应包含补偿所以是加在目标函数里。如果你是站在发电商角度建模目标可能变成收益最大化那另当别论。我们这里按照系统运营商最小化总社会成本的角度来建模。3.3 关键参数计算与工况选择为什么选这么小的规模因为要演示精确求解而不是靠蒙。小规模能够枚举验证和大规模共享同一个数学结构。你可以先用这个小模型测试你的大 M 参数再放大规模。我们先手动看一眼负荷曲线400、550、650、500。如果只有这三台机组G1 的最大出力 300G2 250G3 200总最大出力 750足够应对 650。最小出力 604030130远低于 400所以可行域是存在的。接下来考验的是上升约束。我们从 1 时段开始。假设初始状态所有机组在 0 时段都是停机u_{i,0}0。因此 1 时段如果开机需要启动。启动速率假设为 200 MW/h这意味着机组在启动后一个小时内最多出力从 0 升到 200但还有最小出力要求呢实际上启动时段出力通常必须大于等于最小出力但速率上限可能不允许从0直接跳到最小出力。这可能是不可行的。所以实际中我们会要求启动前如果状态为 0 且之前一段时间未运行则有一个最小停机时间约束。这里为了简化假设 0 时段机组已经运行且出力为各自最小出力但那样负荷400又不能满足所以通常会在第一时段引入初始状态参数。我们这里设初始出力为G160, G240, G330总出力130负荷400差值270必须靠启动机组增加出力。G1 爬坡上限120G2 100G3 80最大增加和是300270可行。但如果全部开机上升约束依然可行。所以我们选择这个初始状态是有意让解有趣。3.4 求解器配置与代码级细节我用 Gurobi 作为示例因为它的 MIQP 和 MICP 支持很成熟。但在用之前必须把互补约束自己展开。下面是核心伪码不是完整程序但足够说明逻辑for i in I: for t in T: u[i,t] model.addVar(vtypeGRB.BINARY, namefu_{i}_{t}) v[i,t] model.addVar(vtypeGRB.BINARY, namefv_{i}_{t}) w[i,t] model.addVar(vtypeGRB.BINARY, namefw_{i}_{t}) p[i,t] model.addVar(lb-GRB.INFINITY, vtypeGRB.CONTINUOUS, namefp_{i}_{t}) # 逻辑关系 model.addConstr(u[i,t] - u[i,t-1] v[i,t] - w[i,t]) model.addConstr(v[i,t] w[i,t] 1) # 爬坡约束注意 t 从1开始 model.addConstr(p[i,t] - p[i,t-1] RU[i] * u[i,t-1] SU[i] * v[i,t]) # 补偿触发变量与 G1 的 p120 绑定 model.addConstr(p[G1,t] - 120 BIG_M * (1 - y[G1,t])) model.addConstr(y[G1,t] u[G1,t]) model.addConstr(y[G1,t] 0) # y 是二进制BIG_M 的取值很关键。对于补偿触发约束p 上限是 300所以 p - 120 最大 180BIG_M 取 200 就够了。别取 10000。更关键的是如果你的下层是市场出清你要写出对偶约束。此时需要使用强对偶等式例如Σ 边际成本 * p - Σ 负荷 * λ 0这是把双层转化为单层的关键。我建议先在小模型上单独测试这个强对偶等式是否成立——如果对偶约束写错了解出来会发现目标值和你手工算的对不上。3.5 结果解读怎么判断求出来的均衡“合不合理”求解结束后不要直接采纳结果。先检查几个东西每个时段的机组启停状态是不是和负荷波动匹配例如负荷最高峰时段时段3650 MW应该保证可开机的高效率机组G1、G2开着。补偿项有没有被滥用如果 G1 一直低出力补偿触发 y1同时机组状态 u1目标函数里的补偿成本增加。系统会理性地选择不要无谓地开 G1 低出力除非 G2、G3 爬坡不够必须 G1 顶上。上升约束是否被触发检查每台机组相邻时段的出力差是否达到爬坡上限。如果某条约束乘子大于零说明该约束是紧的市场价格会被这类约束抬高这正是“上升问题”对市场的真实影响。最后也是最重要的验证这个解是不是均衡。固定解中的整数变量然后只优化连续变量看原问题目标值是否一致。如果一致说明整数解是稳定的再尝试翻转任意一个整数变量的值查看目标函数是否变大最小化问题如果是说明这个整数解是局部最优。如果在全局枚举所有小规模组合后都能确认那就是精确解。4. 常见问题与排查技巧实录4.1 大M取值不当一不小心就“太松”或“太紧”大 M 太紧会误删可行解。比如上面的补偿触发约束如果 M 取 100而 p-120 可能最多 180那么当 p300 时左边 180右边最终会变成 180 100 ... 其实表达式是 p-120 M*(1-y)当 y0 时右边 M必须是 180 100 才能让 y0 可行。这显然不成立因此 y 会被强制等于 1导致即使 G1 高出力也要求补偿——目标函数多了钱解偏保守。你可能会发现结果居然还有解但成本虚高。大 M 太松比如取 100000问题数值稳定性奇差Gurobi 的分支定界里会出现非常多“零整数”判定收敛很慢。经验做法每个互补约束独立取最小的安全上界。先写个小脚本扫描所有可行解统计该表达式的最大绝对值再加一个 10% 的余量。4.2 补偿项导致目标函数非凸求解器报“infeasible”有一次我加补偿项时用了分段线性函数其中某段是凹的情况目标函数变成非凸最小化结果 Gurobi 直接报“Model is infeasible”。其实不是真无可做解而是在预求解阶段检测到非凸性后拒绝求解。解决办法如果补偿项与 0-1 变量相乘那么本质上是一个双线性项需要线性化。比如 y u * f(p) 时可以用标准 big-M 引入辅助变量 z 替换添加约束 z u * M、z 0、z f(p)、z f(p) - M(1-u)。这样可以保持 MILP。不要试图直接使用乘积形式。4.3 上升约束和启停变量耦合时数值病态怎么处理上升约束写 p_t - p_{t-1} RUu_{t-1} SUv_t。当 u 和 v 都是 0 的时候右边是 0这意味着 p_t - p_{t-1} 0。但此时 p 应该为 0机组停机没问题。但如果机组处于运行状态且 v0、u_{t-1}1右边是 RU没问题。真正的病态来自 SU 远大于 RU。比如 SU1000RU100两系数相差 10 倍在预求解器里引发的尺度问题容易让求解器把整数变量当成连续变量处理。我的办法是把启动过程拆分得更细引入“启动出力轨迹”变量或者将启动速率约束直接写成 p_t SU * v_t P_i_max * u_{t-1}实际上常见公式是 p_t P_i_max * u_t再补充另一个约束 p_t - p_{t-1} RU (SU - RU) * v_t。这样可以避免右边出现大 M 的松弛。实测下来数值稳定性好很多。4.4 精确解与启发式解差距调查有同行问我“我启发式找到的解跟精确解只差 0.01%是不是可以不用精确求解了”我的回答是看场景。如果是运行决策0.01% 可能价值几十万且违反约束的可能性大。如果只是规划前期评估可以接受。不过要注意启发式解的间隙 0.01% 是相对于它自己能搜到的范围而不是全局最优的间隙。要验证真正间隙还是要求解松弛后的下界再比较。我自己做过一次测试一个 3 机 4 时段的例子遗传算法在某次运行里找到了目标值 26800而精确解是 26500。遗传法看起来只差 1%但仔细检查发现遗传法的启停状态在时段2和3之间多开了一台 G3导致空载成本高却因为补偿机制而掩盖了。这种“假均衡”在人工检查时很难分辨。4.5 速查表几类典型症状和对应处理症状可能原因排查/处理方法求解器报 infeasible 但人工判断有解大M太小截断可行域二元变量约束写反检查每个 big-M 约束调大并验证可行性解出后补偿变量 y 全是 1 或全是 0 且目标异常补偿触发约束逻辑反了验证 y 与 p 的关系打印未补贴时的解爬坡约束乘子巨大上升约束与启停耦合松弛太大采用 SU-RU 拆分的改进形式整数变量解出来了但连续子问题无界下层对偶约束缺少某些边界条件加入强对偶等式检查对偶变量范围分支定界卡死下界一直不涨大M太大导致 LP 松弛太弱收紧每个互补约束的 M改为有效上界解在数学上最优但实际不满足市场规则缺少最小开机/停机时间约束导致频繁启停添加最小运行/停机时间约束5. 最后再分享一个小的经验技巧很多人在建模时会把补偿放到约束里而不是目标函数里比如“机组出力低于 120 时不可以申请补偿”或“必须支付补偿”这种写法容易造成可行域变形。准确做法是将补偿视为一个可选触发项在目标函数里把它作为成本纳入最小化用二元变量表示触发条件。这样求解器会在成本和物理约束之间自动权衡而不是死板地满足某个规则。这比我见过的一些论文模型要干净得多。还有所有均衡模型写完以后我都建议做一次“纯枚举验证”小规模问题比如 3 机 4 时段所有启停状态最多 2^{12}4096 个遍历每个状态解连续 LP找到全局最优。把这个穷举结果和你的 MILP 精确解对比只要一致你才能放心部署到更大规模场景。这个步骤看起来笨但却是防止模型写错最管用的手段。我在项目里被坑过太多次因为一个小符号的错整个对偶条件悄悄失守如果不是穷举验证根本发现不了。哪怕规模放大到 10 台机组你也可以用随机抽样的方式验证不要嫌麻烦。
