ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

基于线性决策规则的分布鲁棒优化机组组合Matlab实现

基于线性决策规则的分布鲁棒优化机组组合Matlab实现 搞电力调度的人最怕的就是风电出力预测不准带来的连锁反应。你按照预报值把明天96个时段的机组开机计划排好后半夜风突然小了火电还顶着最小技术出力系统频率往下掉备用被吃掉一大块后半夜风又突然大了火电没办法那么快压下来只能眼睁睁看着弃风。这种场景在新能源渗透率高的电网里几乎是常态。传统确定性机组组合对这种随机性基本没有招架之力于是大家开始把不确定性处理直接塞进调度模型。这篇文章要分享的就是目前学术和工程界都很有分量的做法基于线性准则线性决策规则LDR的分布鲁棒优化机组组合并用Matlab完整落地实现。它解决的核心问题是风电出力到底怎么波动、误差分布到底长什么样这个不确定性源头下的最优决策适合电力系统优化方向的研究生、做调度算法开发的工程师、以及对两阶段鲁棒优化感兴趣的程序员参考。1. 从确定性机组组合到不确定性环境下的建模逻辑1.1 经典机组组合的本质和痛点机组组合Unit Commitment简称UC是电力日前调度里最核心的优化问题本质上是问两件事明天哪些机组要开机以及开着的机组每个时段发多少电。决策变量分两类一类是0/1整数变量代表机组启停状态另一类是连续变量代表每台机组在某个时段的出力。目标函数通常由启动成本、空载成本和燃料成本三部分组成约束条件包括功率平衡、旋转备用、出力上下限、最小开停机时间、爬坡速率等等。标准模型是混合整数线性规划或混合整数二次规划规模一放大到IEEE 118节点这种量级变量数量就奔着十万以上去了对建模质量和求解器性能都非常有考验。我的体会是经典UC最常见的坑反而不是模型本身难建而是参数设置不严谨。比如爬坡约束的单位问题冷启动和热启动成本的区分最小开停机时间的递推初始化这些细节稍微写错求解器就会给出一个看似合理、实则无法落地的调度计划。所以在做不确定性问题之前我强烈建议先把确定性UC在几个标准算例上跑通再往上叠加不确定性建模否则出了问题都没法定位到底是基础模型写错了还是鲁棒化过程引入的bug。1.2 风电不确定性让调度变得多难风电出力的不确定性来自天气过程的随机性。预测模型给一条期望出力曲线但实际出力会在预测值附近波动波动幅度大时能到额定容量的百分之二三十。当风电场数量多、分布广时预测误差之间还带有空间相关性整个系统的总误差不再是单机误差的简单叠加。把这些不确定性写进机组组合模型意味着排计划时必须给系统留出足够余量既要能应对风电出力不足时的缺电风险也要能应对出力盈余时的调峰压力。确定性模型处理这种问题的传统办法是加固定比例备用比如让开机机组的可用容量比负荷预测高一个固定的百分比。这招简单有效但在风电占比高的系统里固定比例备用要么太保守要么太激进。更关键的是备用约束只卡住容量这一个维度根本不关心风电误差的分布形态也就没法回答最有可能发生的坏事是什么、概率有多大这类问题。要做更精细的调度必须把风电的随机特性真正放进优化模型而不是简单加一个备用系数了事。1.3 三种主流处理不确定性的策略对比学术界处理调度不确定性的主流框架有三类随机优化、传统鲁棒优化、分布鲁棒优化。我用一个表格先把它们放在一起对比后面再做详细解释。方法需要的信息核心数学结构优点缺点随机优化精确概率分布或场景集场景期望下的混合整数规划充分利用分布信息场景少失真、场景多计算爆炸对分布假设敏感传统鲁棒优化不确定量的支撑区间盒式约束下的min-max问题模型简单、计算可控不考虑概率信息结果偏保守分布鲁棒优化历史样本模糊集半径min-max-min结构可转化为确定性MILP兼顾分布信息与最坏情况稳健性可调模糊集构造影响结果转化过程较复杂用一个生活化的类比帮助理解。随机优化就像出门前看天气预报说下雨概率70%你带把伞传统鲁棒优化则是不管天气预报怎么报都穿全套雨衣雨靴因为最坏情况就是下大暴雨分布鲁棒优化介于两者之间——参考天气预报但同时又考虑万一这个天气预报本身不太准真实天气状况可能会偏到哪个方向然后针对这种不确定性再做一个更保守的决策。三种方法没有绝对优劣关键看场景、数据条件和保守度接受程度。1.4 引入线性决策规则的动机在风电不确定环境下的机组组合天然是一个两阶段决策问题。第一阶段在日前提前安排机组启停第二阶段等实际风电出力揭晓后再对火电机组出力做实时调整。严格来说第二阶段决策应该是实际风电出力向量到调整量的任意函数关系这个函数空间是无限维的求解器根本没有办法直接处理。线性决策规则在标题里叫线性准则就是把这个任意函数限制成仿射函数。举个例子火电机组i的实时调整量写成dg_i a_i b_i乘以风电预测误差其中a_i和b_i是待求的常数系数。这样一来第二阶段策略就只由有限个系数完全刻画无限维的函数空间自然塌缩成一个有限维的系数空间。从电力系统运行的角度看这也不是一个脱离实际的假设调度中心安排机组参与AGC调节时本来就是按偏差比例分配调节量本质上就是一种线性规则所以线性决策规则这类方法不仅仅是数学上的权宜之计它也有合理的工程背景。2. 分布鲁棒机组组合的模型与线性决策规则2.1 模糊集给概率分布上保险分布鲁棒优化要解决的第一个问题是我该怎么描述真实分布到底在哪里。答案就是构造一个集合把可能的概率分布都圈在集合里这个集合就叫模糊集。工程上最常看到的是两类模糊集矩约束模糊集和Wasserstein距离模糊集。矩约束模糊集限制分布的均值、协方差等矩信息在一定范围内。它的好处是表达式简洁对偶转换相对直接缺点是只刻画低阶矩很多形态完全不同的分布在同一个模糊集里是难以区分的。Wasserstein距离模糊集则是以历史样本构造的经验分布为中心以某个半径ε为半径的Wasserstein球。从几何意义上讲Wasserstein距离度量的是把一个分布搬移成另一个分布的最小运输成本所以这个模糊集可以直观控制真实分布和样本分布差多远。实际项目里我更倾向于用Wasserstein模糊集因为它直接吃历史样本不需要假设分布形态而且半径的大小和鲁棒程度是一一对应的调参方向非常清楚。统计学里还有一个有用的结论样本量越大真实分布越可能落在一个半径跟样本量倒数平方根同阶的Wasserstein球里。这可以作为初始半径的上界参考后面再用交叉验证修正。2.2 两阶段模型与min-max-min结构加入模糊集之后分布鲁棒机组组合的完整结构是一个三层的min-max-min问题。最外层min对应第一阶段的机组启停决策这部分在日前阶段必须提前敲定。最内层min对应第二阶段实时调整即风电出力揭晓后在已知基点出力的前提下微调各机组。夹在中间的max是模糊集内最坏分布的选取我们希望在最不利的现实分布下期望运行成本仍然尽可能低。三层套在一起直观的理解就是启停方案可以先定下来但必须保证无论现实按哪种合理分布发展通过最优实时调整都能把成本控制在可接受范围。单独看这个min-max-min结构直接求解是不可能的。求解器的能力边界在有限维确定性规划所以必须想办法把无穷维的分布变量和函数变量消掉。模糊集的对偶变换对应了max这一层线性决策规则对应了内层min这一半两个工具搭在一起才能把整个问题压成可解的确定性规划。这个点理解透了后面看代码和推导都会顺畅很多。2.3 线性决策规则如何化简模型用线性决策规则替换第二阶段决策后实时调整量变成不确定性变量ξ的仿射函数d D0 D1乘以ξ其中D0是常数向量D1是系数矩阵都是待优化变量。把这种表达式代入功率平衡约束、爬坡约束、出力上下限约束以后每条约束都变成含随机变量ξ的线性不等式。对等式约束来说因为要求对任意ξ都成立所以ξ的常数项和一次项系数必须分别满足等式约束这个处理比较直接。对不等式约束来说则需要做鲁棒化处理把对任意ξ都满足不等式的要求转换成一组确定的约束。在分布鲁棒框架下这个转换通常会引入与模糊集半径相关的惩罚项惩罚项最终会和目标函数的期望成本联合在一起形成新的确定性表达式。如果你不用分布鲁棒而用盒式支撑集那其实就是退化成了传统鲁棒优化所有点都要求满足表达式更简单但结果更保守。2.4 对偶转换后得到的确定性MILP经过模糊集对偶和线性决策规则替换这两个步骤原来的min-max-min三层结构被压成一个确定性的混合整数线性规划。目标函数变成第一阶段启动与空载成本加上第二阶段期望成本的上界约束条件包含原本的机组运行约束、LDR系数相关约束以及对偶变换产生的辅助变量约束。这个确定性MILP可以直接交给Gurobi或CPLEX求解。这里要说明一点对偶转换的数学推导步骤比较多实际建模时可以利用YALMIP这类建模工具减少手算负担但前提是你要清楚模型里哪些变量是决策变量、哪些是不确定参数、哪些是辅助对偶变量。如果只是照抄别人的代码一旦遇到无解或者结果不对很难定位问题。我自己做这个项目时最开始就是手推一遍小规模的6节点算例然后才在Matlab里扩展后续调试才没那么痛苦。3. Matlab实现从数据准备到求解器调用3.1 代码架构与文件模块划分这类优化代码最忌讳所有东西堆在一个脚本里。我建议把项目拆成几个模块main主程序管整体流程case_data负责定义机组参数和负荷风电数据build_deterministic_uc负责构建确定性UC模型build_dro_constraints负责叠加模糊集、LDR系数和对偶约束solve_milp负责调用求解器并检查状态plot_results负责出图。文件划分清晰了后面改参数、换数据、对比不同模糊集半径时能省下大量时间。文件模块职责main.m主流程数据、建模、求解、后处理全链路case_data.m机组参数、负荷曲线、风电历史数据build_deterministic_uc.m构建确定性UC的变量、目标、约束build_dro_constraints.m追加模糊集、线性决策规则与对偶约束solve_milp.m设置求解器参数并执行优化plot_results.m绘制启停状态、出力和后验对比图模块化还有一个隐藏的好处调试时可以先用deterministic版本验证基础模型完全正确再打开DRO开关。出问题时定位范围一下就缩小了——是先确认基础模型能跑通再去研究新增的分布鲁棒约束哪里写得不对。3.2 数据准备风电样本与机组参数数据准备第一件事是统一单位。风电出力、负荷、机组容量尽量都换算成同一个基准值。我一般用100 MVA做标幺化这样目标函数数值不至于因为数量级差异在求解器里出现数值病态。我踩过一次很深的坑负荷数据用MW风电数据也用了MW但忘记了把某个风电场的历史数据减去基准值导致功率平衡约束在数值上总有几百的残差求解结果画出来曲线看不出问题一校验才发现特定时段严重不满足平衡。风电不确定性样本的构造方式取决于你手上的数据。如果有历史实测出力可以直接把预测误差作为样本如果只有预测和实测两条曲线可以按小时构造误差样本并考虑时序相关性。具体到线性决策规则稳妥的做法是把每个时段的风电误差都作为多维随机变量ξ的一个分量样本集就是历史误差向量的集合。用经验分布作为模糊集的中心再结合交叉验证选半径要比随手指定一个半径可靠得多。如果原始数据里有明显的坏数据或者通信中断导致的零值记得先做数据清洗否则模糊集中心会被垃圾样本带偏。3.3 YALMIP建模与求解器配置Matlab里最省事的优化建模路线是YALMIP加Gurobi或CPLEX。YALMIP负责把数学表达式翻译成求解器输入Gurobi负责解混合整数线性规划。安装和配置不展开说了重点讲代码逻辑和容易出错的地方。%% 主流程核心代码片段YALMIP Gurobi nUnits 3; % 机组数 nPeriods 24; % 时段数 nSamples 500; % 风电误差样本数 u binvar(nUnits, nPeriods, full); % 机组启停 p0 sdpvar(nUnits, nPeriods, full); % 基点出力 a0 sdpvar(nUnits, nPeriods, full); % LDR常数项系数 A1 sdpvar(nUnits*nPeriods, nSamples, full); % LDR对风电误差的线性系数实际建模需按系统维度展平 %% 目标函数启动成本 空载成本 燃料成本 分布鲁棒惩罚项 objective startupCost(u) noLoadCost(u) fuelCost(p0) droPenalty(a0, A1, radius, samples); %% 约束功率平衡、机组限值、爬坡、最小启停时间、分布鲁棒约束 constraints [powerBalance, unitLimit, rampLimit, minUpDownTime, droConstraints]; ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.001, gurobi.TimeLimit, 600); optimize(constraints, objective, ops);代码的核心是变量层级要规划清楚。u和p0是第一层变量直接出现在目标函数和确定性约束里a0和A1是线性决策规则层的系数几乎每个与风电误差相关的约束里都要用到。写代码时要特别注意A1的索引顺序别把时段、机组和样本的顺序搞混。我建议每定义一个变量都写清注释否则这类三维以上的数据展开在调试时真的会让人怀疑人生。分布式鲁棒惩罚项在代码里通常不是一个简单的函数而是通过对偶辅助变量构成的线性项放在独立函数里更利于测试。3.4 结果分析与后处理要点模型求解完之后别急着出图按顺序做四件事。第一检查求解器返回状态确认MIPgap是否在可接受范围。第二检查机组启停状态是否满足最小开停时间约束有时候初始状态设置不当会出现机组连续启停的不合理序列。第三把最优解代回原始的两阶段问题做蒙特卡洛回验在随机风电样本下统计系统实际的失负荷量和弃风量。这个回验非常重要因为它才能说明这个分布鲁棒调度方案在模拟真实运行中到底靠不靠谱。第四做不同Wasserstein半径的敏感性分析画出半径—期望成本/最坏成本曲线看结果的鲁棒性如何随保守程度变化。后处理是论文里最有说服力的部分。我的习惯是把确定性UC、传统鲁棒优化、DRO不同半径的三组方案放在同一张图上横轴是风电出力的波动幅度纵轴是系统运行成本或者失负荷量。这样能直观展示DRO相对传统方法的优势和成本代价也能在审稿人面前展现出对模型特性的深入理解。4. 调试经验与常见问题速查4.1 求解慢的排查方案求解慢是这类项目最常见的难题但慢的原因往往不在模型而在于建模质量。我按经验列一个排查顺序。第一看MIP下界是不是长时间不动。如果下界一点改善都没有多半是约束里大M常数太多或者太大尝试把所有大M压缩到恰好够用的量级求解速度通常立刻提升。第二看整数变量规模是不是过大。机组数量多时可以尝试对称性破缺把同类型机组合并成机组簇或者添加对称性破缺约束能有效压缩搜索空间。第三调整求解器参数。Gurobi的MIPFocus可以设为1或2Threads适当增加对中小规模算例一般有正向收益。分布式鲁棒模型求解慢还有一个来源模糊集半径设置过大导致对偶辅助变量数量膨胀线性规划松弛质量变差。如果你只是做原理验证先在3机6节点系统上测试通再逐步放大到IEEE 30节点和118节点别一上来就拿大系统调参那样很容易被漫长的求解时间磨掉耐心。4.2 Wasserstein半径的调节技巧半径的选择本质上是风险偏好的量化。太小模型几乎完全信任样本分布相当于退化成随机优化太大模型倾向于为极端概率分布做准备成本一路走高。我的调节思路是两步走。先用较粗的网格扫描半径比如从0.01到0.5观察成本曲线变化找到曲线斜率变化最明显的拐点然后在拐点附近加密扫描。拐点处的半径通常是性价比最高的选择——再增大半径成本大幅上升但鲁棒性收益有限再缩小半径鲁棒性明显下降但成本节省不多。如果手头数据是真实风电历史出力半径还可以参考样本量做上界估计。样本量越大经验分布越接近真实分布半径应取得越小。这种基于统计直觉的调参方式写进论文里比人为设定半径要显得扎实得多。另外提醒一句每次改完半径都要重新跑一遍蒙特卡洛回验因为纸上谈兵的最坏分布成本和实际随机样本下的表现不一定完全对应。4.3 无可行解问题分析遇到无可行解我的习惯不是盯着报错信息看而是先做松弛定位。具体做法是逐个去掉约束组看看哪一组约束拿掉之后问题重新可行。根据经验常见原因有这么几个。最小开停机时间和爬坡约束在小系统中特别容易打架尤其当负荷曲线峰谷差大、机组最小启停时间又长时备用约束和机组最大出力之间的配合不好导致某个时段总可用容量低于负荷加备用还有一个容易被忽略的情况模糊集半径设得太大、样本里存在极端点会让对偶辅助变量的可行域被冲散。另外线性决策规则系数A1如果没有设置对应的物理上下限优化器可能会利用变量自由度产生技术可行但实际反物理的调度方案。做可行性测试的时候把机组出力值打印出来和上下限肉眼比对一下很多问题一眼就能看出来。如果发现某台机组出力长期贴在边界上同时启停状态还在频繁切换大概率是初始状态或者爬坡约束出了问题。4.4 隐藏在实现细节中的坑整理几个我实测下来最容易踩的坑。第一YALMIP中二元变量和连续变量相乘时需要用big-M做线性化这个M的值必须基于物理边界设置得尽量紧。如果M太大线性松弛后的可行域过宽求解过程会异常缓慢M太小又会切掉真实可行解导致结果不对。第二Gurobi许可证和Matlab版本之间的兼容性装不上就换Mosek或者用YALMIP自带求解器但默认求解器对大规模MILP非常吃力不建议用在大算例上。第三做蒙特卡洛后验时抽样样本数量一定要大于模糊集中心的样本数量不然你验证的稳健性可能只是对少量样本的重放统计意义不足。第四所有风电相关数据进入模型前务必统一为标幺值。我因为这个吃过一次大亏某节点的功率平衡残差显示500多MW排查了半天最后发现只是一个风电场的数据漏了基准值换算。这些坑单个看都不大但叠加起来足以让一个好端端的模型变成看起来很对、跑起来就废的鸡肋。做项目时最好把数据校验和结果校验写成独立脚本每次改动模型后自动跑一遍在早期就能拦住绝大部分问题。4.5 可以继续扩展的方向这个模型还有不少值得扩展的空间。实际系统往往同时存在负荷不确定性、风电不确定性甚至电价和检修计划也是随机的把多源不确定性都纳入模糊集是更贴近现实的建模方式。另一个有价值的方向是把线性决策规则升级为分段线性决策规则模型会变大但保守性会明显降低。还有些团队把分布鲁棒机组组合和储能规划、需求响应联合优化在计算复杂度可控的前提下有效降低系统运行成本。如果你打算在这个方向发论文这几个角度比单纯换算例改参数更有含金量。我个人做完这个模型的体会是数学推导再漂亮最终还是要落到代码细节上才见真章。分布鲁棒、线性决策规则这些名词听起来很前沿但实现过程中的很多问题其实是基础性的单位换算了没有变量索引对没对大M松不松MIPgap设得合不合理。把这些基础层面的东西打磨扎实模型自然就能跑得动、也跑得稳。希望这篇文章能帮你少走一点弯路把风电不确定性调度这个方向的研究和工程落地做得更顺手。
返回列表