企业官网建设流程全解析

做电力市场优化这些年，我最怕看到的题目关键词就是“两级市场”“风险”“强对偶”凑在一起——不是怕难，是怕绕。省间交易商要在省间、省内两级市场里倒腾电量，价格随场景波动，购电成本不仅看期望值，还要看尾部风险；而市场出清环节又是一个内生优化问题，交易商的决策会反过来影响出清价格。这套耦合关系放到MATLAB+Cplex里实现，最关键的坎就是怎么把下层市场出清问题用强对偶转成单层约束。最近我把《两级电力市场环境下计及风险的省间交易商最优购电模型》这套思路完整跑了一遍，从建模、推导到代码落地踩了不少坑，这里把全过程拆开讲透。

这篇文章适合三类人看：一是做电力市场、综合能源交易相关课题的研究生，需要复现论文或自己搭双层优化模型；二是刚接触Cplex+Yalmip、对强对偶转化只会背公式不知道怎么下手的初学者；三是已经跑通基础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相比普通蒙特卡洛抽样的优势在于分层采样，能以更少的样本覆盖整个输入空间，避免随机抽样带来的样本聚集问题。具体过程分三步：

1. 对每个时段的价格随机变量，将其累积分布函数等分为N个区间，在每个区间内随机取一个代表点；
2. 对所有时段的代表点做随机排列组合，生成初始场景集；
3. 用同步回代消除算法（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 对偶问题与原-对偶约束

线性规划强对偶定理告诉我们：如果原问题有最优解，那么对偶问题也有最优解，且两个最优目标函数值相等。反过来，如果原问题可行、对偶问题可行、且原目标等于对偶目标，那么这两个解分别是各自问题的最优解。

这个性质的价值在于：我们不需要显式求解下层问题，只需要保证存在一组原变量和对偶变量满足以下三个条件：

1. 原问题约束（原可行）
2. 对偶问题约束（对偶可行）
3. 原目标函数值 = 对偶目标函数值（强对偶等式）

三个条件全部写进上层模型，下层优化就“退化”为一组约束。上层模型的解自然会让这组约束成立，因为交易商不会选择一个导致市场无法出清的决策。

强对偶等式是这一步的核心，但也是麻烦的来源——它包含原变量与对偶变量的乘积项（比如λ × D_total，其中D_total是交易商购电量），这是双线性项，不是线性的。直接交给Cplex，如果用的是Yalmip，会得到“检测到双线性项”的报错。

3.3 互补松弛条件的线性化与大M法

处理双线性项的标准思路是放弃强对偶等式，改用KKT互补松弛条件。对原问题中的每个约束，写出其对应的对偶变量，然后补充以下条件：

对偶变量 × (约束松弛量) = 0

这个等式本身还是非线性的（乘积为0），但可以用大M法线性化。以g_min ≤ g ≤ g_max为例，引入二进制变量δ_min和δ_max：

• g - 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. MATLAB+Cplex代码实现全流程

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值、每个场景下的成本分位数、以及省间出清电价水平。把这组数据画出来，能看到风险中性（beta=0）和风险规避（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常见报错

Yalmip+Cplex的组合我用了很多年，报错大多集中在三个位置。

第一，调用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%左右。这套代码的后续扩展方向其实很多：可以加入多时段耦合的抽水蓄能、加入绿电交易环境价值、或者把省内市场出清过程的网络约束细化——只要下层模型还是线性规划，强对偶的整个推导框架都能直接复用，不需要从零开始重写建模逻辑。