风光荷不确定性是这几年电力系统调度里绕不开的话题做优化调度研究的同学基本都会碰到。CVaR条件风险价值这个从金融领域来的风险度量指标被引入电力系统后和鲁棒优化结合处理风光出力和负荷预测误差已经成了一种很主流的研究范式。标题里这个“基于 CVaR 的风光负荷不确定性鲁棒调度优化研究Matlab代码实现”实际上就是要解决一个很具体的问题在风光和负荷都有预测误差的情况下怎么做一个既经济又抗风险的机组调度方案并且用 Matlab 把整套模型跑通。这篇博文我就从模型思路、数学推导、代码实现到调试经验完整拆一遍这个课题。这套东西适合谁呢如果你的研究方向是电力系统优化调度、微电网经济调度、新能源并网运行或者正在做涉及不确定性的优化建模那这篇文章正好能帮你把 CVaR 和鲁棒优化这条技术路线从理论落到代码。就算你只是刚接触 Yalmip 和求解器跟着这篇文章走一遍也能把整套框架搭起来。1. 问题拆解风光荷不确定性下的调度为什么要用 CVaR 加鲁棒优化1.1 风光荷不确定性的本质与调度困境风电和光伏的出力本质上是气象驱动的风速的间歇性、云层遮挡的随机性直接导致风电场和光伏电站的出力无法被精确预测。负荷侧也没好到哪里去虽然负荷预测技术发展了很多年但居民的用电行为、极端天气带来的制冷制热需求都会让实际负荷偏离预测曲线。这三个不确定性叠加在一起调度员面临的就不是一个确定的优化问题而是一个带随机参数的决策问题。传统的确定性调度只考虑预测场景把风电、光伏、负荷都当成已知的固定值求解一个标准的线性规划或混合整数规划。这种方式在新能源渗透率低的时候问题不大预测误差带来的偏差可以由系统备用容量兜底。但随着风光装机容量快速增加预测误差的绝对值越来越大确定性调度方案的鲁棒性就出问题了——实际运行中可能因为风电出力比预测低、负荷比预测高导致系统需要紧急调用高成本机组甚至切负荷。所以调度模型必须显式地考虑不确定性这也是为什么随机规划、鲁棒优化、分布鲁棒优化这些方法在电力系统调度里越来越受重视。1.2 CVaR 的引入从金融风控到电力调度的迁移逻辑CVaR 全称 Conditional Value at Risk中文一般叫条件风险价值。它最早用在金融投资组合里度量的是投资组合在给定置信水平下发生尾部损失时的平均损失程度。和它经常一起出现的 VaRValue at Risk风险价值定义是分位数比如 95% 置信水平下的 VaR 意味着有 95% 的概率损失不会超过这个值。但 VaR 有个众所周知的缺点——它只关注分位点本身完全忽略超过分位点之后尾部损失的严重程度。CVaR 恰好补上了这个短板它计算的是损失超过 VaR 之后的条件期望也就是尾部损失的平均水平。在电力系统调度里我们可以把“损失”理解成调度方案在不确定性场景下偏离预期成本的那部分。风光出力不足、负荷攀升时系统需要调用高成本的备用机组或者因为切负荷产生惩罚费用这些都属于尾部风险。用 CVaR 度量这种风险就可以在目标函数里显式地纳入对极端场景的考量决策者通过调整风险偏好系数在“期望成本最低”和“最坏情况可控”之间做权衡。这正是 CVaR 相比传统期望值模型的核心优势——它不追求所有场景下都最优而是保证了风险场景下不出现灾难性损失。1.3 鲁棒优化与 CVaR 的互补关系鲁棒优化的核心思想是在不确定参数属于某个集合的前提下寻找一个对所有集合内场景都可行的调度方案目标是最小化最坏情况下的成本。它的优势是计算效率高模型通常能转化为等价的确定性优化问题尤其是不确定集用盒式、椭球式或多面体式描述时可以通过对偶理论或线性化手段直接求解。但它也有一个天然的问题——鲁棒优化只关心最坏场景可能导致方案过于保守经济性变差。CVaR 和鲁棒优化结合恰好是互补的。鲁棒优化的保守性来自它对所有不确定集合内的场景都一视同仁但实际上极端场景出现的概率很低为之付出过高的成本并不划算。CVaR 的引入可以看作是在“鲁棒性”和“经济性”之间插了一个旋钮通过调节风险偏好系数让模型自动决定要在多大程度上为极端场景买单。从数学上看如果风险偏好系数取到足够大CVaR 模型会退化为纯鲁棒模型取到 0则退化为期望值模型。所以这个框架实际上是统一了确定性与鲁棒性之间的整个谱系。2. 数学建模不确定集、CVaR 表达式与完整约束推导2.1 风光荷不确定集的构建方式要做鲁棒优化第一步是定义不确定参数的取值范围。最常用的是盒式不确定集形式很直观每个不确定参数的真实值在预测值附近的一个区间内波动。比如风电出力可以写成W_t ∈ [W_t^pre - ΔW_t, W_t^pre ΔW_t]其中 W_t^pre 是 t 时段的预测出力ΔW_t 是波动偏差上限。光伏和负荷也做同样的处理。盒式集合简单、线性但边界太“硬”把所有参数同时取到边界的最坏情况往往是不现实的所以工程上通常会引入一个预算参数 Γbudget of uncertainty来限制波动参数偏离预测值的总个数或总程度。这个思想来源是 Bertsimas 和 Sim 的经典鲁棒优化框架它能在不确定集的保守程度和计算复杂度之间取得很好的平衡。不确定集还有椭球式和多面体式等选择。椭球式集合能更好地刻画参数之间的相关性但会引入二阶锥约束增加求解难度多面体式介于盒式和椭球式之间。我的经验是对初学者来说先不要贪多把盒式或盒式预算约束的不确定集吃透能调通代码再考虑更复杂的集合结构。实际研究里最常用也最容易写进论文的是带预算参数的盒式集合既严谨又方便做保守性分析。2.2 CVaR 的数学定义与场景化表达CVaR 的标准数学定义如下。设 f(x, ξ) 是决策变量 x 在随机参数 ξ 下的损失函数给定置信水平 β比如 0.95VaR 定义为VaR_β(x) min{ α | P(f(x, ξ) ≤ α) ≥ β }CVaR 则是损失超过 VaR 之后的条件期望CVaR_β(x) E[ f(x, ξ) | f(x, ξ) ≥ VaR_β(x) ]直接算这个条件期望很麻烦但 Rockafellar 和 Uryasev 在 2000 年提出了一个关键定理CVaR 可以等价表示为一个凸优化的最小值问题CVaR_β(x) min_{α} { α (1/(1-β)) · E[ max(f(x, ξ) - α, 0) ] }这个形式太重要了它把 CVaR 转换成了一个线性规划可处理的结构配合蒙特卡洛抽样的场景集可以很容易地离散化。假设我们有 N 个抽样场景 ξ_1, ξ_2, ..., ξ_N那么 CVaR 的样本估计是CVaR_β(x) ≈ min_{α} { α (1/((1-β)N)) · Σ_{k1}^N max(f(x, ξ_k) - α, 0) }引入辅助变量 z_k max(f(x, ξ_k) - α, 0)可以线性化为z_k ≥ f(x, ξ_k) - αz_k ≥ 0这是标准套路后面 Matlab 建模时直接照这个写就行。2.3 目标函数期望成本与 CVaR 风险成本的加权有了上面的数学工具调度模型的目标函数就可以写成min (1 - λ) · E[Cost] λ · CVaR_β(Cost)其中 λ 是风险偏好权重从 0 到 1 变化。λ 0 时是纯期望值模型λ 1 时是最小化尾部风险也就是接近最坏情况决策。Cost 是调度总成本包括常规机组发电成本、启停成本、弃风和弃光惩罚、以及失负荷惩罚。在场景化建模中期望成本就是所有场景成本的平均值。用脚本记住一个要点CVaR 模型里 α 是一个决策变量它自己也有经济含义——近似对应置信水平下的 VaR 值。有些论文里会把 α 直接解释成“风险成本阈值”方便做结果分析。实际实现中我建议把 α 初始值设成所有场景成本的均值附近有助于求解器更快收敛。2.4 系统约束功率平衡、机组出力与备用约束调度模型的约束分几类。第一类是功率平衡约束每个时段系统总发电等于总负荷加网损在场景化模型里每个场景、每个时段都要满足。第二类是机组技术约束火电机组有功出力要在一个区间内爬坡速率要受限启停状态和最小启停时间也要满足。这些属于混合整数约束是模型求解难度的主要来源。第三类是备用约束这是处理不确定性最直接的手段。传统做法是要求系统上备用和下调备用分别大于风光预测误差的某个置信区间比如上备用容量大于风电和光伏出力下降、负荷上升的最坏组合。引入 CVaR 之后备用约束可以和风险度量结合起来——不是保证所有情况下都有足够备用而是保证在置信水平 β 覆盖的场景集合内有足够备用用风险指标来刻画备用的充裕度。这样做的好处是避免为了极低概率事件配置过度冗余的备用容量经济性明显提升。3. Matlab 实现从数学模型到可运行代码3.1 代码架构与文件组织一套完整的 Matlab 实现我建议按下面这个结构组织文件├── main.m % 主程序参数设置、场景生成、建模求解、结果输出 ├── data/ │ ├── wind_data.m % 风电预测与误差参数 │ ├── pv_data.m % 光伏预测与误差参数 │ └── load_data.m % 负荷预测与误差参数 ├── model/ │ ├── build_model.m % 构建 Yalmip 优化模型 │ ├── add_cvar.m % 添加 CVaR 相关变量和约束 │ └── add_uc_ed.m % 添加机组组合与经济调度约束 ├── scenarios/ │ └── gen_scenarios.m % 蒙特卡洛场景生成与缩减 └── utils/ ├── plot_results.m % 结果可视化 └── print_table.m % 生成结果表格这种模块化的好处是很明显的调试的时候只需要关注单个文件改参数不用翻整个模型写论文做灵敏度分析也方便。很多新手喜欢把所有代码堆在一个 main 脚本里图省事但到了要改模型加约束的时候就会非常痛苦。数据、模型、求解、可视化分开是一种值得从一开始就养成的习惯。主程序 main.m 的执行顺序大致是加载数据、生成场景、构建模型、调用求解器、提取结果、可视化。每一步之间要留好接口特别是场景生成的输出格式要和模型构建函数能对上这是最容易出 bug 的地方。3.2 场景生成蒙特卡洛抽样与场景缩减场景生成是整个不确定性建模的基础环节。简单做法是假设风电预测误差服从正态分布或 t 分布用 Matlab 自带的随机数生成函数抽样。但这里有个工程细节值得注意风电预测误差严格来说并不是均值为零的正态分布偏度和峰度在不同风速区间差异很大。如果只是做方法验证用正态分布也够但如果要发高质量论文建议用基于历史数据的预测误差分布或者用核密度估计方法拟合误差分布再进行采样。我提供一个比较实用的场景生成框架。先对每个时段的风电、光伏、负荷分别抽样 N 个独立误差值然后叠加到预测曲线上得到 N 个完整的不确定场景。为了避免场景数量过多导致计算量爆炸通常需要做场景缩减常用的算法是快速前向选择法Fast Forward Selection在 Matlab File Exchange 上有现成代码可以借鉴。缩减后保留的场景能覆盖原始场景集的主要概率特征同时把求解规模降下来。具体来说蒙特卡洛抽样的代码逻辑大致是% 抽样得到 N 个场景 N 500; wind_scenarios zeros(N, T); for t 1:T error normrnd(0, sigma_w(t), N, 1); wind_scenarios(:, t) wind_forecast(t) error; end % 同理生成光伏和负荷场景场景缩减之后每个场景要配一个概率权重所有场景权重之和为 1。在目标函数的期望成本计算里用权重乘以对应场景成本即可。注意概率权重的归一化处理否则会出现数值上的偏差。3.3 Yalmip 建模CVaR 约束的标准写法Matlab 里做优化建模最省力的工具就是 Yalmip。它是一个建模层把优化模型的描述和具体求解器解耦你可以先写模型最后再指定用 Cplex 还是 Gurobi 求解。下面这段是 CVaR 相关变量和约束在 Yalmip 里的标准写法% 决策变量定义 alpha sdpvar(1, 1); % VaR 阈值变量 z sdpvar(N_scen, 1); % 每个场景的超过量 % 每个场景的总成本 cost_by_scen 是 sdpvar 表达式 % CVaR 约束 Constraints [Constraints, z 0]; Constraints [Constraints, z cost_by_scen - alpha]; % 目标函数期望成本 风险权重 * CVaR cvar_term alpha (1 / ((1 - beta) * N_scen)) * sum(z); Objective (1 - lambda) * mean(cost_by_scen) lambda * cvar_term;这里 cost_by_scen 是每个场景下的系统总成本它在每个场景里都是关于机组出力、启停状态等决策变量的线性表达式。把 CVaR 相关项加到目标函数之后整个模型就是一个混合整数线性规划或混合整数二次规划取决于目标函数里是否包含机组的二次成本函数。Yalmip 建模的坑主要在维度匹配上。sdpvar 定义的变量维度和约束矩阵维度必须严格一致否则会报维度不匹配错误或者更隐蔽地给出一个错误的模型。遇到这种情况先用 size 命令检查所有变量维度再用 assign 命令给变量赋初值排查哪个约束构造不对。3.4 求解器配置与参数调优Matlab 的 linprog 和 intlinprog 内置求解器可以处理规模适中的混合整数线性规划但大型调度问题的求解效率非常差。建议用 Cplex 或 Gurobi这两个求解器都支持 Yalmip 接口学术版免费。以 Gurobi 为例在求解前设置好参数能很大程度提升效率options sdpsettings(solver, gurobi, ... verbose, 2, ... gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 3600, ... gurobi.NumThreads, 4);MIPGap 设为 0.01 意味着求解器可以在 1% 的相对最优间隙内停止这在实际调度中完全够用。TimeLimit 设为 3600 秒是防止模型太复杂导致求解时间不可控。NumThreads 设成你机器的逻辑核心数多线程求解单机就能明显加速。很多论文习惯把 MIPGap 设成 0即求得最优整数解才停这在中小规模算例上没问题但规模一大等待时间会非常感人。我的实际经验是先跑一个宽松的 MIPGap比如 0.05验证模型正确性再逐步收紧到 0.005 或 0.001 做最终计算。这样能在保证结果可靠的前提下大幅节省调试时间。4. 调试经验、结果分析与避坑指南4.1 常见报错与排查思路这个课题的代码实现过程中有几类报错出现频率特别高。第一个是 Yalmip 维度不匹配问题这通常是因为场景矩阵和决策变量矩阵的形状没对齐用 size 逐个变量查一遍就能定位。第二个是求解器返回“infeasible”不可行原因多半是约束条件之间自相矛盾——比如备用约束要求上备用大于某个值但机组最大出力不满足这个要求。排查方法是先去掉 CVaR 相关约束跑一遍看基础模型是否可行再加回来锁定是哪些约束组合导致了不可行。第三个是求解器报“numerical issues”数值问题在混合整数模型里很常见。解决思路是检查模型中的系数量级是否统一。如果风光出力是兆瓦级、成本是元级量级差出几个数量级就容易触发数值问题。一个实际经验是把所有物理量统一到标幺值或者按同一个数量级归一化数值稳定性会明显改善。还有一个容易被忽视的问题Yalmip 在建模二次成本函数如机组成本为 a·P² b·P c时如果二次项里混入了整数变量模型会变成 MIQP求解难度大幅上升。工程上常用的方法是分段线性化piecewise linearization把二次成本函数近似成若干段线性函数把模型退化成 MILP求解速度能提高一到两个数量级。4.2 结果分析怎么从数值结果里提取有价值的结论模型跑通之后不能只是得到几个数字就完事要有逻辑地对结果进行分析这也是论文审稿人最看重的内容。首先是 λ 的参数灵敏度分析让风险偏好系数从 0 变化到 1画出期望成本和 CVaR 之间的权衡曲线。这条曲线通常叫什么有效前沿efficient frontier。在曲线上可以清楚看到随着风险偏好增加期望成本抬升但尾部风险下降决策者可以根据对风险的承受能力选择合适的工作点。其次是置信水平 β 的影响分析。β 增大意味着关注的尾部概率变小CVaR 度量的是更极端情况下的平均损失此时系统的备用需求会增加运行成本上升。在论文里这张不同 β 下的运行成本对比表几乎是标配它证明了 CVaR 模型对风险定义参数的一致性。最后是将所提出的模型和两个基准模型做对比一个是纯确定性模型无视不确定性一个是纯鲁棒模型只优化最坏情况。确定性模型的成本最低但可能在极端场景下失负荷纯鲁棒模型最保守成本也最高本文模型居中且通过调节 λ 可以灵活过渡。这个对比实验是证明方法有效性的关键证据一定要做扎实。4.3 几个提升模型性能的实用技巧调试和优化过程中我积累了几个很实用的技巧。第一个是给决策变量赋初值。混合整数规划求解器通常内置了启发式算法来找初始可行解但如果你自己定义了一个合理的初值特别是机组启停状态变量可以大幅缩短求解时间。Yalmip 里用 assign 给初始变量赋值再给 sdpsettings 里的 usex0 设为 1 就能启用。第二个技巧是尽量合并同类约束。比如功率平衡约束如果写成逐场景、逐时段的数百条约束可以用矩阵运算一次性构造Yalmip 支持向量化约束这一点对减少建模开销和有约束爆炸的问题都很有帮助。对 loop 循环逐条添加约束会让构建模型的时间变得很长模型规模大的时候甚至会成为瓶颈。第三个技巧是做场景缩减后再求解。500 个场景和 50 个场景的求解时间差距不是线性的而是呈指数级增长。很多方法验证只需要证明模型和算法的有效性不一定需要特别多的场景场景缩减后跑出来的结果差异通常很小但时间节省非常可观。还有一点要提醒的是代码实现过程中一定要加注释。Matlab 本身是脚本语言重新运行成本低但如果你中间隔了一周再回来看代码没有注释的时候真的会一头雾水。关键公式对应的代码行旁边最好把论文里的公式标号也写上这样调试和写论文时一目了然。4.4 成果扩展方向从单目标到多目标、从离线到在线这个模型做完之后后续扩展空间非常大。一个方向是引入多目标优化把运行成本和碳排放同时作为优化目标用 epsilon-约束法或 NSGA-II 求解 Pareto 前沿分析碳交易机制对调度决策的影响。另一个方向是结合分布鲁棒优化DRO用 Wasserstein 距离构造模糊集在不需要精确概率分布的情况下进行调度这是目前优化领域非常活跃的方向。再往应用层面走可以把这个静态调度模型扩展到多时间尺度比如日前调度加日内滚动修正的两阶段结构。日前阶段用本文的 CVaR 鲁棒模型求出机组组合日内阶段随着风光和负荷的超短期预测不断刷新用滚动优化修正出力计划。这种结构更贴近调度生产的实际流程做出来后续发论文和应用都更有说服力。我自己做这个课题时最大的体会是风光不确定性的调度优化重点不全在模型有多高级而在于你对不确定性建模的假设是否合理、对风险度量的含义是否吃得透。CVaR 和鲁棒优化结合的好处是可以把“风险”这个抽象概念量化到可优化的目标函数里让调度方案在成本和安全之间得到真正可控的平衡。你可以从最简单的盒式不确定集加 CVaR 目标起步一步步把模型扩展到自己论文需要的样子。Matlab 加 Yalmip 这套组合现在仍然是做算法验证和仿真实验最高效的工具链没有之一。 基于CVaR的风光负荷不确定性鲁棒调度优化Matlab实现
风光荷不确定性是这几年电力系统调度里绕不开的话题做优化调度研究的同学基本都会碰到。CVaR条件风险价值这个从金融领域来的风险度量指标被引入电力系统后和鲁棒优化结合处理风光出力和负荷预测误差已经成了一种很主流的研究范式。标题里这个“基于 CVaR 的风光负荷不确定性鲁棒调度优化研究Matlab代码实现”实际上就是要解决一个很具体的问题在风光和负荷都有预测误差的情况下怎么做一个既经济又抗风险的机组调度方案并且用 Matlab 把整套模型跑通。这篇博文我就从模型思路、数学推导、代码实现到调试经验完整拆一遍这个课题。这套东西适合谁呢如果你的研究方向是电力系统优化调度、微电网经济调度、新能源并网运行或者正在做涉及不确定性的优化建模那这篇文章正好能帮你把 CVaR 和鲁棒优化这条技术路线从理论落到代码。就算你只是刚接触 Yalmip 和求解器跟着这篇文章走一遍也能把整套框架搭起来。1. 问题拆解风光荷不确定性下的调度为什么要用 CVaR 加鲁棒优化1.1 风光荷不确定性的本质与调度困境风电和光伏的出力本质上是气象驱动的风速的间歇性、云层遮挡的随机性直接导致风电场和光伏电站的出力无法被精确预测。负荷侧也没好到哪里去虽然负荷预测技术发展了很多年但居民的用电行为、极端天气带来的制冷制热需求都会让实际负荷偏离预测曲线。这三个不确定性叠加在一起调度员面临的就不是一个确定的优化问题而是一个带随机参数的决策问题。传统的确定性调度只考虑预测场景把风电、光伏、负荷都当成已知的固定值求解一个标准的线性规划或混合整数规划。这种方式在新能源渗透率低的时候问题不大预测误差带来的偏差可以由系统备用容量兜底。但随着风光装机容量快速增加预测误差的绝对值越来越大确定性调度方案的鲁棒性就出问题了——实际运行中可能因为风电出力比预测低、负荷比预测高导致系统需要紧急调用高成本机组甚至切负荷。所以调度模型必须显式地考虑不确定性这也是为什么随机规划、鲁棒优化、分布鲁棒优化这些方法在电力系统调度里越来越受重视。1.2 CVaR 的引入从金融风控到电力调度的迁移逻辑CVaR 全称 Conditional Value at Risk中文一般叫条件风险价值。它最早用在金融投资组合里度量的是投资组合在给定置信水平下发生尾部损失时的平均损失程度。和它经常一起出现的 VaRValue at Risk风险价值定义是分位数比如 95% 置信水平下的 VaR 意味着有 95% 的概率损失不会超过这个值。但 VaR 有个众所周知的缺点——它只关注分位点本身完全忽略超过分位点之后尾部损失的严重程度。CVaR 恰好补上了这个短板它计算的是损失超过 VaR 之后的条件期望也就是尾部损失的平均水平。在电力系统调度里我们可以把“损失”理解成调度方案在不确定性场景下偏离预期成本的那部分。风光出力不足、负荷攀升时系统需要调用高成本的备用机组或者因为切负荷产生惩罚费用这些都属于尾部风险。用 CVaR 度量这种风险就可以在目标函数里显式地纳入对极端场景的考量决策者通过调整风险偏好系数在“期望成本最低”和“最坏情况可控”之间做权衡。这正是 CVaR 相比传统期望值模型的核心优势——它不追求所有场景下都最优而是保证了风险场景下不出现灾难性损失。1.3 鲁棒优化与 CVaR 的互补关系鲁棒优化的核心思想是在不确定参数属于某个集合的前提下寻找一个对所有集合内场景都可行的调度方案目标是最小化最坏情况下的成本。它的优势是计算效率高模型通常能转化为等价的确定性优化问题尤其是不确定集用盒式、椭球式或多面体式描述时可以通过对偶理论或线性化手段直接求解。但它也有一个天然的问题——鲁棒优化只关心最坏场景可能导致方案过于保守经济性变差。CVaR 和鲁棒优化结合恰好是互补的。鲁棒优化的保守性来自它对所有不确定集合内的场景都一视同仁但实际上极端场景出现的概率很低为之付出过高的成本并不划算。CVaR 的引入可以看作是在“鲁棒性”和“经济性”之间插了一个旋钮通过调节风险偏好系数让模型自动决定要在多大程度上为极端场景买单。从数学上看如果风险偏好系数取到足够大CVaR 模型会退化为纯鲁棒模型取到 0则退化为期望值模型。所以这个框架实际上是统一了确定性与鲁棒性之间的整个谱系。2. 数学建模不确定集、CVaR 表达式与完整约束推导2.1 风光荷不确定集的构建方式要做鲁棒优化第一步是定义不确定参数的取值范围。最常用的是盒式不确定集形式很直观每个不确定参数的真实值在预测值附近的一个区间内波动。比如风电出力可以写成W_t ∈ [W_t^pre - ΔW_t, W_t^pre ΔW_t]其中 W_t^pre 是 t 时段的预测出力ΔW_t 是波动偏差上限。光伏和负荷也做同样的处理。盒式集合简单、线性但边界太“硬”把所有参数同时取到边界的最坏情况往往是不现实的所以工程上通常会引入一个预算参数 Γbudget of uncertainty来限制波动参数偏离预测值的总个数或总程度。这个思想来源是 Bertsimas 和 Sim 的经典鲁棒优化框架它能在不确定集的保守程度和计算复杂度之间取得很好的平衡。不确定集还有椭球式和多面体式等选择。椭球式集合能更好地刻画参数之间的相关性但会引入二阶锥约束增加求解难度多面体式介于盒式和椭球式之间。我的经验是对初学者来说先不要贪多把盒式或盒式预算约束的不确定集吃透能调通代码再考虑更复杂的集合结构。实际研究里最常用也最容易写进论文的是带预算参数的盒式集合既严谨又方便做保守性分析。2.2 CVaR 的数学定义与场景化表达CVaR 的标准数学定义如下。设 f(x, ξ) 是决策变量 x 在随机参数 ξ 下的损失函数给定置信水平 β比如 0.95VaR 定义为VaR_β(x) min{ α | P(f(x, ξ) ≤ α) ≥ β }CVaR 则是损失超过 VaR 之后的条件期望CVaR_β(x) E[ f(x, ξ) | f(x, ξ) ≥ VaR_β(x) ]直接算这个条件期望很麻烦但 Rockafellar 和 Uryasev 在 2000 年提出了一个关键定理CVaR 可以等价表示为一个凸优化的最小值问题CVaR_β(x) min_{α} { α (1/(1-β)) · E[ max(f(x, ξ) - α, 0) ] }这个形式太重要了它把 CVaR 转换成了一个线性规划可处理的结构配合蒙特卡洛抽样的场景集可以很容易地离散化。假设我们有 N 个抽样场景 ξ_1, ξ_2, ..., ξ_N那么 CVaR 的样本估计是CVaR_β(x) ≈ min_{α} { α (1/((1-β)N)) · Σ_{k1}^N max(f(x, ξ_k) - α, 0) }引入辅助变量 z_k max(f(x, ξ_k) - α, 0)可以线性化为z_k ≥ f(x, ξ_k) - αz_k ≥ 0这是标准套路后面 Matlab 建模时直接照这个写就行。2.3 目标函数期望成本与 CVaR 风险成本的加权有了上面的数学工具调度模型的目标函数就可以写成min (1 - λ) · E[Cost] λ · CVaR_β(Cost)其中 λ 是风险偏好权重从 0 到 1 变化。λ 0 时是纯期望值模型λ 1 时是最小化尾部风险也就是接近最坏情况决策。Cost 是调度总成本包括常规机组发电成本、启停成本、弃风和弃光惩罚、以及失负荷惩罚。在场景化建模中期望成本就是所有场景成本的平均值。用脚本记住一个要点CVaR 模型里 α 是一个决策变量它自己也有经济含义——近似对应置信水平下的 VaR 值。有些论文里会把 α 直接解释成“风险成本阈值”方便做结果分析。实际实现中我建议把 α 初始值设成所有场景成本的均值附近有助于求解器更快收敛。2.4 系统约束功率平衡、机组出力与备用约束调度模型的约束分几类。第一类是功率平衡约束每个时段系统总发电等于总负荷加网损在场景化模型里每个场景、每个时段都要满足。第二类是机组技术约束火电机组有功出力要在一个区间内爬坡速率要受限启停状态和最小启停时间也要满足。这些属于混合整数约束是模型求解难度的主要来源。第三类是备用约束这是处理不确定性最直接的手段。传统做法是要求系统上备用和下调备用分别大于风光预测误差的某个置信区间比如上备用容量大于风电和光伏出力下降、负荷上升的最坏组合。引入 CVaR 之后备用约束可以和风险度量结合起来——不是保证所有情况下都有足够备用而是保证在置信水平 β 覆盖的场景集合内有足够备用用风险指标来刻画备用的充裕度。这样做的好处是避免为了极低概率事件配置过度冗余的备用容量经济性明显提升。3. Matlab 实现从数学模型到可运行代码3.1 代码架构与文件组织一套完整的 Matlab 实现我建议按下面这个结构组织文件├── main.m % 主程序参数设置、场景生成、建模求解、结果输出 ├── data/ │ ├── wind_data.m % 风电预测与误差参数 │ ├── pv_data.m % 光伏预测与误差参数 │ └── load_data.m % 负荷预测与误差参数 ├── model/ │ ├── build_model.m % 构建 Yalmip 优化模型 │ ├── add_cvar.m % 添加 CVaR 相关变量和约束 │ └── add_uc_ed.m % 添加机组组合与经济调度约束 ├── scenarios/ │ └── gen_scenarios.m % 蒙特卡洛场景生成与缩减 └── utils/ ├── plot_results.m % 结果可视化 └── print_table.m % 生成结果表格这种模块化的好处是很明显的调试的时候只需要关注单个文件改参数不用翻整个模型写论文做灵敏度分析也方便。很多新手喜欢把所有代码堆在一个 main 脚本里图省事但到了要改模型加约束的时候就会非常痛苦。数据、模型、求解、可视化分开是一种值得从一开始就养成的习惯。主程序 main.m 的执行顺序大致是加载数据、生成场景、构建模型、调用求解器、提取结果、可视化。每一步之间要留好接口特别是场景生成的输出格式要和模型构建函数能对上这是最容易出 bug 的地方。3.2 场景生成蒙特卡洛抽样与场景缩减场景生成是整个不确定性建模的基础环节。简单做法是假设风电预测误差服从正态分布或 t 分布用 Matlab 自带的随机数生成函数抽样。但这里有个工程细节值得注意风电预测误差严格来说并不是均值为零的正态分布偏度和峰度在不同风速区间差异很大。如果只是做方法验证用正态分布也够但如果要发高质量论文建议用基于历史数据的预测误差分布或者用核密度估计方法拟合误差分布再进行采样。我提供一个比较实用的场景生成框架。先对每个时段的风电、光伏、负荷分别抽样 N 个独立误差值然后叠加到预测曲线上得到 N 个完整的不确定场景。为了避免场景数量过多导致计算量爆炸通常需要做场景缩减常用的算法是快速前向选择法Fast Forward Selection在 Matlab File Exchange 上有现成代码可以借鉴。缩减后保留的场景能覆盖原始场景集的主要概率特征同时把求解规模降下来。具体来说蒙特卡洛抽样的代码逻辑大致是% 抽样得到 N 个场景 N 500; wind_scenarios zeros(N, T); for t 1:T error normrnd(0, sigma_w(t), N, 1); wind_scenarios(:, t) wind_forecast(t) error; end % 同理生成光伏和负荷场景场景缩减之后每个场景要配一个概率权重所有场景权重之和为 1。在目标函数的期望成本计算里用权重乘以对应场景成本即可。注意概率权重的归一化处理否则会出现数值上的偏差。3.3 Yalmip 建模CVaR 约束的标准写法Matlab 里做优化建模最省力的工具就是 Yalmip。它是一个建模层把优化模型的描述和具体求解器解耦你可以先写模型最后再指定用 Cplex 还是 Gurobi 求解。下面这段是 CVaR 相关变量和约束在 Yalmip 里的标准写法% 决策变量定义 alpha sdpvar(1, 1); % VaR 阈值变量 z sdpvar(N_scen, 1); % 每个场景的超过量 % 每个场景的总成本 cost_by_scen 是 sdpvar 表达式 % CVaR 约束 Constraints [Constraints, z 0]; Constraints [Constraints, z cost_by_scen - alpha]; % 目标函数期望成本 风险权重 * CVaR cvar_term alpha (1 / ((1 - beta) * N_scen)) * sum(z); Objective (1 - lambda) * mean(cost_by_scen) lambda * cvar_term;这里 cost_by_scen 是每个场景下的系统总成本它在每个场景里都是关于机组出力、启停状态等决策变量的线性表达式。把 CVaR 相关项加到目标函数之后整个模型就是一个混合整数线性规划或混合整数二次规划取决于目标函数里是否包含机组的二次成本函数。Yalmip 建模的坑主要在维度匹配上。sdpvar 定义的变量维度和约束矩阵维度必须严格一致否则会报维度不匹配错误或者更隐蔽地给出一个错误的模型。遇到这种情况先用 size 命令检查所有变量维度再用 assign 命令给变量赋初值排查哪个约束构造不对。3.4 求解器配置与参数调优Matlab 的 linprog 和 intlinprog 内置求解器可以处理规模适中的混合整数线性规划但大型调度问题的求解效率非常差。建议用 Cplex 或 Gurobi这两个求解器都支持 Yalmip 接口学术版免费。以 Gurobi 为例在求解前设置好参数能很大程度提升效率options sdpsettings(solver, gurobi, ... verbose, 2, ... gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 3600, ... gurobi.NumThreads, 4);MIPGap 设为 0.01 意味着求解器可以在 1% 的相对最优间隙内停止这在实际调度中完全够用。TimeLimit 设为 3600 秒是防止模型太复杂导致求解时间不可控。NumThreads 设成你机器的逻辑核心数多线程求解单机就能明显加速。很多论文习惯把 MIPGap 设成 0即求得最优整数解才停这在中小规模算例上没问题但规模一大等待时间会非常感人。我的实际经验是先跑一个宽松的 MIPGap比如 0.05验证模型正确性再逐步收紧到 0.005 或 0.001 做最终计算。这样能在保证结果可靠的前提下大幅节省调试时间。4. 调试经验、结果分析与避坑指南4.1 常见报错与排查思路这个课题的代码实现过程中有几类报错出现频率特别高。第一个是 Yalmip 维度不匹配问题这通常是因为场景矩阵和决策变量矩阵的形状没对齐用 size 逐个变量查一遍就能定位。第二个是求解器返回“infeasible”不可行原因多半是约束条件之间自相矛盾——比如备用约束要求上备用大于某个值但机组最大出力不满足这个要求。排查方法是先去掉 CVaR 相关约束跑一遍看基础模型是否可行再加回来锁定是哪些约束组合导致了不可行。第三个是求解器报“numerical issues”数值问题在混合整数模型里很常见。解决思路是检查模型中的系数量级是否统一。如果风光出力是兆瓦级、成本是元级量级差出几个数量级就容易触发数值问题。一个实际经验是把所有物理量统一到标幺值或者按同一个数量级归一化数值稳定性会明显改善。还有一个容易被忽视的问题Yalmip 在建模二次成本函数如机组成本为 a·P² b·P c时如果二次项里混入了整数变量模型会变成 MIQP求解难度大幅上升。工程上常用的方法是分段线性化piecewise linearization把二次成本函数近似成若干段线性函数把模型退化成 MILP求解速度能提高一到两个数量级。4.2 结果分析怎么从数值结果里提取有价值的结论模型跑通之后不能只是得到几个数字就完事要有逻辑地对结果进行分析这也是论文审稿人最看重的内容。首先是 λ 的参数灵敏度分析让风险偏好系数从 0 变化到 1画出期望成本和 CVaR 之间的权衡曲线。这条曲线通常叫什么有效前沿efficient frontier。在曲线上可以清楚看到随着风险偏好增加期望成本抬升但尾部风险下降决策者可以根据对风险的承受能力选择合适的工作点。其次是置信水平 β 的影响分析。β 增大意味着关注的尾部概率变小CVaR 度量的是更极端情况下的平均损失此时系统的备用需求会增加运行成本上升。在论文里这张不同 β 下的运行成本对比表几乎是标配它证明了 CVaR 模型对风险定义参数的一致性。最后是将所提出的模型和两个基准模型做对比一个是纯确定性模型无视不确定性一个是纯鲁棒模型只优化最坏情况。确定性模型的成本最低但可能在极端场景下失负荷纯鲁棒模型最保守成本也最高本文模型居中且通过调节 λ 可以灵活过渡。这个对比实验是证明方法有效性的关键证据一定要做扎实。4.3 几个提升模型性能的实用技巧调试和优化过程中我积累了几个很实用的技巧。第一个是给决策变量赋初值。混合整数规划求解器通常内置了启发式算法来找初始可行解但如果你自己定义了一个合理的初值特别是机组启停状态变量可以大幅缩短求解时间。Yalmip 里用 assign 给初始变量赋值再给 sdpsettings 里的 usex0 设为 1 就能启用。第二个技巧是尽量合并同类约束。比如功率平衡约束如果写成逐场景、逐时段的数百条约束可以用矩阵运算一次性构造Yalmip 支持向量化约束这一点对减少建模开销和有约束爆炸的问题都很有帮助。对 loop 循环逐条添加约束会让构建模型的时间变得很长模型规模大的时候甚至会成为瓶颈。第三个技巧是做场景缩减后再求解。500 个场景和 50 个场景的求解时间差距不是线性的而是呈指数级增长。很多方法验证只需要证明模型和算法的有效性不一定需要特别多的场景场景缩减后跑出来的结果差异通常很小但时间节省非常可观。还有一点要提醒的是代码实现过程中一定要加注释。Matlab 本身是脚本语言重新运行成本低但如果你中间隔了一周再回来看代码没有注释的时候真的会一头雾水。关键公式对应的代码行旁边最好把论文里的公式标号也写上这样调试和写论文时一目了然。4.4 成果扩展方向从单目标到多目标、从离线到在线这个模型做完之后后续扩展空间非常大。一个方向是引入多目标优化把运行成本和碳排放同时作为优化目标用 epsilon-约束法或 NSGA-II 求解 Pareto 前沿分析碳交易机制对调度决策的影响。另一个方向是结合分布鲁棒优化DRO用 Wasserstein 距离构造模糊集在不需要精确概率分布的情况下进行调度这是目前优化领域非常活跃的方向。再往应用层面走可以把这个静态调度模型扩展到多时间尺度比如日前调度加日内滚动修正的两阶段结构。日前阶段用本文的 CVaR 鲁棒模型求出机组组合日内阶段随着风光和负荷的超短期预测不断刷新用滚动优化修正出力计划。这种结构更贴近调度生产的实际流程做出来后续发论文和应用都更有说服力。我自己做这个课题时最大的体会是风光不确定性的调度优化重点不全在模型有多高级而在于你对不确定性建模的假设是否合理、对风险度量的含义是否吃得透。CVaR 和鲁棒优化结合的好处是可以把“风险”这个抽象概念量化到可优化的目标函数里让调度方案在成本和安全之间得到真正可控的平衡。你可以从最简单的盒式不确定集加 CVaR 目标起步一步步把模型扩展到自己论文需要的样子。Matlab 加 Yalmip 这套组合现在仍然是做算法验证和仿真实验最高效的工具链没有之一。