做电力市场优化这些年我最怕看到的题目关键词就是“两级市场”“风险”“强对偶”凑在一起——不是怕难是怕绕。省间交易商要在省间、省内两级市场里倒腾电量价格随场景波动购电成本不仅看期望值还要看尾部风险而市场出清环节又是一个内生优化问题交易商的决策会反过来影响出清价格。这套耦合关系放到MATLABCplex里实现最关键的坎就是怎么把下层市场出清问题用强对偶转成单层约束。最近我把《两级电力市场环境下计及风险的省间交易商最优购电模型》这套思路完整跑了一遍从建模、推导到代码落地踩了不少坑这里把全过程拆开讲透。这篇文章适合三类人看一是做电力市场、综合能源交易相关课题的研究生需要复现论文或自己搭双层优化模型二是刚接触CplexYalmip、对强对偶转化只会背公式不知道怎么下手的初学者三是已经跑通基础LP模型想进一步把风险度量CVaR和策略性购电加进自己代码里的工程师。已经熟悉双层规划基本概念的同学可以直接跳到第3节看强对偶的落地方案还在起步阶段的建议从头读每一步我都标了“为什么这么做”。1. 先说清楚这个模型到底在优化什么1.1 两级市场里交易商的生存逻辑先捋一下业务背景。我国电力市场正在从省间、省内两级市场协同运作的方向演进省间市场负责跨省资源优化配置省内市场负责省内的电力电量平衡。省间交易商就是夹在中间的“倒爷”——不过这个“倒爷”要遵循严格的电力市场规则。交易商的典型动作是在省间市场比如省间日前交易买入电量然后在省内市场通过合约、现货或增量挂牌等方式卖出或自用。买入价格取决于省间市场的出清结果卖出价格由省内市场决定两边都有不确定性。两级市场的价格信号互相传导交易商一旦在省间买了高价电省内卖不出去就要自己消化成本。更麻烦的是省间市场的出清价格不是外生给定的——如果交易商申报的购电量大可能推动节点边际电价上升这就是策略性交易行为也正是论文里“最优”二字的含义不是被动接受价格而是主动选择购电量来影响出清实现自身利益最大化。这里有个容易混淆的点普通用户的购电模型里价格是外生参数优化问题只有一个决策变量就是买多少。但省间交易商面对的是市场出清机制价格由全市场的供需决定交易商的申报量会进入出清模型的约束条件中。于是问题天然变成了双层结构上层是交易商决策下层是市场出清。这个结构在数学上不好直接求解强对偶就是用来拆解它的工具。1.2 风险度量为什么绕不开CVaR购电成本包含两部分一部分是确定性成本比如固定合约的分解曲线另一部分是现货市场购电成本由场景化的现货价格决定。场景一多成本的分布就出来了——不是一个确定数而是一条概率分布曲线。如果只优化期望成本模型会倾向于在低价场景多买在高价场景少买但实际运行时你不知道会落到哪个场景。万一落到高价场景购电成本可能比期望值高出几个量级这就是尾部风险。VaR风险价值能告诉你“最坏情况下有5%的概率损失超过某个值”但它不关心超过之后到底亏多少。CVaR条件风险价值更进一步它度量的是“超过VaR那5%场景下的平均损失”数学性质也比VaR好得多——CVaR是凸的、相合的coherent在离散场景下还能用线性约束精确表达天然适合放进LP/MILP框架里求解。具体到本模型决策变量里的购电量为连续变量目标函数是期望购电成本和CVaR的加权和权重由交易商的风险偏好系数决定。风险偏好系数越大模型越保守越倾向于牺牲期望成本来压低尾部损失。1.3 两层优化结构从哪来前面提到下层市场出清问题是一个标准的线性规划在满足电力平衡、线路潮流限额如果模型里考虑电网拓扑、机组出力上下限等约束下最小化系统购电成本。交易商的上层决策变量会出现在下层模型的约束条件中比如购电申报量进入省间市场的负荷需求项。数学上可以写成如下形式上层min购电成本期望 λ × CVaR决策变量是购电申报量下层min 系统出清成本决策变量是机组出力、节点电价等且下层问题中包含上层决策变量上下层变量耦合直接套用商业求解器无法求解。常规做法有两种一是用KKT条件替换下层问题二是用强对偶替换下层问题。KKT条件需要写出下层的拉格朗日函数并补充互补松弛条件代码量稍大强对偶的路线则是利用线性规划原对偶最优值相等的性质把下层目标函数用对偶目标代换再加上原对偶可行约束和强对偶等式就能把整体变成单层MILP。论文标题里点明“强对偶”走的就是第二条路线也是我认为工程上最容易实现、最不容易出错的一条。2. 数学建模目标函数与约束的完整推导2.1 场景生成与削减计及风险的前提是知道现货价格可能的分布形态现实中我们拿不到连续分布只能通过历史数据生成有限个典型场景。常用的方法有两类一类是直接对历史现货价格做统计拟合抽样生成场景另一类是利用电价预测模型输出每个时段的预测区间再结合蒙特卡洛或拉丁超立方抽样生成场景。我这次用的是拉丁超立方抽样LHS加同步回代消除。LHS相比普通蒙特卡洛抽样的优势在于分层采样能以更少的样本覆盖整个输入空间避免随机抽样带来的样本聚集问题。具体过程分三步对每个时段的价格随机变量将其累积分布函数等分为N个区间在每个区间内随机取一个代表点对所有时段的代表点做随机排列组合生成初始场景集用同步回代消除算法fast forward selection合并相似场景把场景数从几百个削减到十几个同时记录每个场景的概率。场景削减这一步不能省。Cplex求解MILP的时间跟二进制变量数量和约束规模强相关场景每多一倍求解时间可能翻几倍。把200个场景削减到20个目标函数值变化通常在5%以内但求解速度能提升一到两个数量级性价比极高。削减后的每个场景包含一组完整的省间现货价格时间序列同时保留原概率p_s。这样就把随机优化问题转化成了一个确定性的等价问题——用多个确定性场景去逼近随机过程是电力市场优化里最主流的处理方式。2.2 目标函数与CVaR线性化模型的上层目标函数写成min Σs ps × (省间购电成本_s 省内售电收益_s的负值) β × CVaR_α这里的β是风险权重α是置信水平一般取0.95或0.99。省间购电成本_s由场景s的出清价格和购电量相乘得到——注意出清价格本身受购电量影响这正是下层出清模型的输出所以不能简单把价格当作常量。CVaR的线性化是这套模型里最经典也最容易被忽略的一步。引入辅助变量ρ对应VaR值和辅助变量u_s每个场景下成本超过ρ的溢出量CVaR可以写成CVaR ρ (1 / ((1-α) × Σs ps)) × Σs ps × u_s约束为u_s ≥ C_s - ρ u_s ≥ 0其中C_s是场景s下的总购电成本。这样处理之后CVaR就从“排序取分位数”这种不可导的操作变成了一个线性表达式。Cplex求解器完全不需要做任何特殊处理直接把ρ和u_s当普通连续变量放进模型即可。有个细节值得注意ρ的取值范围不需要限制因为它会被u_s和C_s的约束自然拉到一个合理区间但如果想帮求解器加速可以给ρ设一个宽泛的上下界比如[min(C_s), max(C_s)]能明显减少分支定界的搜索范围。2.3 约束条件详细拆解模型的约束从功能上分四块第一块是交易商的电量平衡约束。任意时段t、任意场景s下省间购电量等于省内售电量加自身净负荷这是交易商的“物理守恒”必须严格满足。如果模型允许弃电或回购还需要引入对应的松弛变量并给松弛变量一个惩罚系数否则模型可能通过“凭空消失电量”来“优化”成本。第二块是省间市场的出清约束这部分来自下层模型。包括系统电力平衡约束、机组出力上下限、线路潮流限额等所有约束都带场景下标。因为下层是线性规划这些约束在对偶转化后会以对偶可行约束的形式重新出现。第三块是交易上限约束。省间交易商在某个时点的购电量有市场规则限制比如占省间通道输送能力的比例上限、月度交易电量上限等这部分约束通常写成简单的box constraint但注意要加场景下标否则等于假设所有场景下交易商的运行约束完全一致对风险场景会严重失真。第四块是CVaR相关约束即上节提到的u_s和ρ的关系约束。这部分虽然简短但它是整个模型从“期望优化”升级为“风险优化”的关键。建模时我习惯先把约束分类每类写一个函数最后统一拼装Yalmip的约束对象。这样后面排查模型错误时直接注释掉某一类约束就能定位问题——是电量平衡崩了、市场出清崩了还是风险约束崩了一目了然。3. 强对偶转化把双层问题变成单层问题3.1 下层市场出清问题的原始形式下层省间市场出清问题是一个标准LP目标是最小化系统购电成本决策变量是机组出力g、切负荷量l、以及各节点的相角θ如果考虑直流潮流。为了说清楚对偶转化这里给一个简化的原始形式min Σ (C_g × g V_OLL × l)s.t.系统功率平衡Σg l D_total其中D_total包含交易商的申报购电量机组出力上下限g_min ≤ g ≤ g_max切负荷上限0 ≤ l ≤ l_max线路潮流约束如果简化可以省略把这个LP规范化成min c^T x, s.t. Ax ≤ b的形式然后写出对偶问题。因为原问题包含等式约束功率平衡对偶变量是自由变量λ不等式约束对应的对偶变量要求非负。对偶问题的形式大概是max λ × D_total - μ_max × g_max μ_min × g_min - ...s.t. 对偶可行约束由原问题的列生成这里D_total里包含了交易商的购电量所以上层决策变量会出现在对偶目标函数中——这正是强对偶转化能“穿透”双层结构的关键点。3.2 对偶问题与原-对偶约束线性规划强对偶定理告诉我们如果原问题有最优解那么对偶问题也有最优解且两个最优目标函数值相等。反过来如果原问题可行、对偶问题可行、且原目标等于对偶目标那么这两个解分别是各自问题的最优解。这个性质的价值在于我们不需要显式求解下层问题只需要保证存在一组原变量和对偶变量满足以下三个条件原问题约束原可行对偶问题约束对偶可行原目标函数值 对偶目标函数值强对偶等式三个条件全部写进上层模型下层优化就“退化”为一组约束。上层模型的解自然会让这组约束成立因为交易商不会选择一个导致市场无法出清的决策。强对偶等式是这一步的核心但也是麻烦的来源——它包含原变量与对偶变量的乘积项比如λ × D_total其中D_total是交易商购电量这是双线性项不是线性的。直接交给Cplex如果用的是Yalmip会得到“检测到双线性项”的报错。3.3 互补松弛条件的线性化与大M法处理双线性项的标准思路是放弃强对偶等式改用KKT互补松弛条件。对原问题中的每个约束写出其对应的对偶变量然后补充以下条件对偶变量 × (约束松弛量) 0这个等式本身还是非线性的乘积为0但可以用大M法线性化。以g_min ≤ g ≤ g_max为例引入二进制变量δ_min和δ_maxg - g_min ≤ M × (1 - δ_min)μ_min ≤ M × δ_min这样当δ_min 1时μ_min被迫为0此时约束g ≥ g_min不起作用g可以大于g_min当δ_min 0时g被固定在g_minμ_min可以取正值。注意这里的MM要选得足够大但也不能太大否则数值稳定性会出问题——这是整个模型里最容易翻车的地方我在第5节会详细讲。用互补松弛条件替换强对偶等式之后模型变成MILP包含少量二进制变量。Cplex对MILP的求解效率比对MINLP高得多这也是论文标题里点出“强对偶”背后的工程含义——它把难以处理的均衡约束问题变成了可求解的混合整数问题。转化完成后还要验证一件事原问题的约束是“紧凑”的即不存在冗余约束导致对偶变量不唯一。否则互补条件可能出现多个解上层模型的最优值会不稳定。实操中可以给模型加一个很小的正则项来消除这种退化或者直接检查约束矩阵是否满秩。4. MATLABCplex代码实现全流程4.1 代码框架与模块划分我用的求解环境是MATLAB Yalmip Cplex。Yalmip是一个建模层能把优化问题翻译成Cplex能吃的LP/MILP格式代码可读性比直接调Cplex API高得多。完整代码我建议拆成五个模块main.m主程序定义参数、调用建模函数、求解、输出结果scenario_generation.m场景生成与削减build_upper_model.m上层目标函数与约束build_lower_model.m下层出清模型及对偶约束post_process.m结果分析输出购电计划、成本分布、CVaR值这种模块划分的好处是改场景生成策略、改风险参数、改市场规则约束都只需要动对应函数不用把几百行代码翻来覆去。参数部分有几个值需要提前定好时段数比如24、场景数削减后取20、置信水平0.95、风险权重从0到1扫描、机组参数、负荷数据、省间通道容量上限。建议把所有参数集中写在main.m开头的结构体里方便批量实验。4.2 核心代码实现与讲解下面给出关键建模代码基于Yalmip语法。上层购电量的定义和CVaR约束可以这样写%% 上层变量 q sdpvar(T, S, full); % 交易商购电量T时段时间S个场景 rho sdpvar(1, 1); % VaR辅助变量 u sdpvar(T, S, full); % 尾部溢出变量 delta binvar(...); % 互补松弛用二进制变量按需定义 %% 购电成本省间市场出清价格由下层对偶变量lambda给出 cost sum(sum(lambda .* q)); % lambda为场景化出清电价 %% CVaR约束 Constraints [Constraints, u cost - rho]; Constraints [Constraints, u 0]; %% 目标函数期望成本 beta * CVaR obj sum(ps .* sum(cost, 1)) beta * (rho 1/((1-alpha)*S) * sum(ps .* sum(u, 1)));注意这里的lambda是下层对偶变量它和q是同时被优化的——这就是“价格被内生化”的体现。如果直接把lambda定义为常数模型就退化成普通的价格接受者模型丢失了策略性投标的核心行为。对偶约束部分以系统功率平衡约束为例。原约束是sum(g) l demand q对应自由对偶变量lambda其对偶可行约束由机组成本列生成。代码里可以用Yalmip的dual函数和constraint对象配合但更可控的方式是手写对偶约束%% 下层原问题约束g为机组出力变量 Constraints [Constraints, sum(g, 1) l demand q]; % 功率平衡 Constraints [Constraints, g_min g g_max]; Constraints [Constraints, 0 l l_max]; %% 下层对偶问题约束 Constraints [Constraints, lambda 0]; % 按对偶变量性质定义 % 对每一列原变量写出对应的对偶可行不等式 % 例如机组g对应的对偶约束C_g - lambda mu_max - mu_min 0为了代码整洁我建议把“原问题约束”和“对偶问题约束”写成两个子函数分别返回Yalmip约束对象最后再拼接。这样当模型规模扩大时调试起来非常方便。4.3 求解设置与结果解读模型拼接完成后调用Cplex求解ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.mip.tolerances.integrality 1e-5; ops.cplex.timelimit 3600; optimize(Constraints, obj, ops);这几个参数很关键。mipgap设到1e-4已经足够工程精度设太严会让Cplex陷入长尾搜索integrality是二进制变量的整数容忍度默认1e-5可行时间限制建议一定设上否则一个难解的算例会把你一个下午都吃掉。结果输出部分我关注四个指标购电量计划q的期望值、CVaR值、每个场景下的成本分位数、以及省间出清电价水平。把这组数据画出来能看到风险中性beta0和风险规避beta较大情况下购电曲线的差异——通常风险规避会让交易商在高价时段的购电量更保守整体成本分布更集中。还有一个容易被忽略的指标对偶变量的值。功率平衡约束的对偶变量就是节点边际电价它的经济含义是“在该节点增加1MW负荷时系统成本的变化量”。模型跑完后检查这个值是否落在合理区间能快速发现约束写错或大M值过大的问题。5. 踩坑记录与高效调试技巧5.1 大M取值带来的数值灾难大M法是线性化互补松弛条件最常用的手段但M取值需要极小心。M太小互补条件无法正确表达模型结果错误M太大Cplex的数值稳定性会崩——典型症状是求解器报“infeasible or unbounded”或者连续变量的解在一串无意义的尾数上震荡。我的经验是M不是全局统一值而是按约束逐条取。对于机组出力约束M取“机组最大出力与最小出力的差值”就足够了对于线路潮流约束M取“线路容量的两倍”对于爬坡约束M取“爬坡速率的若干倍”。这样做比全局设一个巨大值要稳得多。另外一个实用技巧是把互补松弛条件优化成“大M值越小越好”的形式。比如原约束是g - g_min ≥ 0松弛量最大只有g_max - g_min那么M就设成g_max - g_min 0.1既能覆盖全部可能取值又不会因为过大导致求解器数值问题。5.2 场景削减与求解时间权衡场景数是最直接影响求解时间的超参数。我做过一组对比20个场景时求解时间约2分钟50个场景时约15分钟100个场景时直接超过1小时——时间呈指数上升主要原因是二进制变量数量跟着场景数一起涨。但场景太少也不行CVaR估计的准确性会大幅下降。一个比较稳妥的做法是先用快速前向选择削减到20个场景做机制验证确定模型逻辑没有错误后再逐步放大到30~50个场景看结果稳定性。如果结果在不同场景数下变化不大购电成本差异小于3%就说明场景数量已经足够收敛。再分享一个实践技巧两层不确定性是可以解耦的。上层交易商的风险决策只需要对价格场景敏感而价格波动的主要来源是负荷和新能源出力。所以场景生成时可以先用关键因子分析法降维再对低维因子做LHS抽样生成速度能快一倍以上。5.3 Cplex求解器设置与Yalmip常见报错YalmipCplex的组合我用了很多年报错大多集中在三个位置。第一调用optimize时提示“No suitable solver for bilinear terms”。这表示模型里还有双线性项没有被消掉。排查办法用Yalmip的expand命令检查目标函数和约束里是否存在乘积项找到后用前面讲的互补松弛大M法处理。第二Cplex报“License Error”或者“No available license”。这不是代码问题是许可证没配置好。确认环境变量CPLEX_HOME指向正确并且MATLAB中Cplex类能正常实例化。第三求解结果里对偶变量出现异常值比如电价达到1e8。这通常是M值过大或模型存在重复约束导致对偶不唯一。先检查约束矩阵的秩和重叠约束再有针对性地调M。调试这类模型的通用方法论我总结成一句话先缩小算例、再放大规模。遇到问题先把时段数降到2、场景数降到2、机组数降到1跑通了再逐步加回来。如果在最小算例上都报错问题一定出在约束逻辑本身而不是计算资源不足。我在实际跑这套模型时最深的一个体会是强对偶转化的难点不在“背出对偶问题”而在“原问题要写规范”——变量分清楚谁是原变量、谁是对偶变量约束分清楚谁对应哪个对偶变量互补松弛条件逐条对上号模型基本一遍过。反过来如果一开始就急着堆代码原问题约束抄漏了一条后面的对偶约束全是错的排查一晚上都不一定找得到根因。最后分享一个提升效率的小技巧给模型加一个“冷启动初始解”。先用风险中性模型跑一个解作为热启动初值传给风险规避模型Cplex在MILP分支定界时能明显加快下界收敛。具体做法是在优化前给二进制变量赋一个合理的试探值assign(delta, 0.5); % 给二进制变量一个初始赋值帮助求解器剪枝注意Yalmip的assign是建模层面的初始值不改变模型本身。配合热启动我实测过整体求解时间能缩短30%左右。这套代码的后续扩展方向其实很多可以加入多时段耦合的抽水蓄能、加入绿电交易环境价值、或者把省内市场出清过程的网络约束细化——只要下层模型还是线性规划强对偶的整个推导框架都能直接复用不需要从零开始重写建模逻辑。