尧图网站建设 
开题那阵子我琢磨微网电源容量配置满脑子都是“怎么让这套系统在现实里扛得住”。风光出力的随机性、负荷的波动任何一个没照顾好之前算的“最优方案”到实际运行就可能变成“最贵方案”甚至“失稳方案”。后来我把目光锁定在两阶段鲁棒优化算法上配合MATLAB、YALMIP和CPLEX这条经典工具链才真正把这个问题从理论推到了可仿真的层面。这篇东西我想把我从建模到代码实现、再到踩坑排错的完整过程拆开讲清楚给正在做相关研究的同学一条能直接上手的路。1. 为什么容量配置必须考虑不确定性确定性模型的局限1.1 传统确定性优化的问题描述与短板在微网电源容量优化里传统做法通常是把风光出力、负荷大小当作已知量来处理。设定几个典型日场景比如“夏季典型日”“冬季典型日”把每个时段的风速、光照、负荷都取一个固定值然后建立优化模型求解出光伏、风电、储能、燃气轮机这些设备的最优安装容量。这类确定性模型的数学形式很简洁[ \min_{C, P} ; C_{inv}(C) C_{ope}(P) ]其中 (C_{inv}) 是设备投资成本(C_{ope}) 是运行成本约束条件包括功率平衡、设备出力上下限、储能充放电约束等。所有参数都是确定值模型直接扔给求解器就能出结果。但这个思路有个致命的隐含假设未来风光出力和负荷曲线就是那几个典型日场景。实际情况呢光照可能连续三天低于预期风速可能整月都处于“要死不活”的状态负荷也可能因为极端天气而飙升。一旦实际运行与预设场景偏差过大之前算好的容量配置就会出现两类问题供电可靠性下降切负荷风险增加运行成本大幅偏离预期甚至出现必须高价购电的局面。我最早用确定性模型做仿真时结果很漂亮总成本最低、新能源渗透率最高。但一换到真实历史数据里回放最恶劣的那几天系统几乎撑不住。这个问题不解决论文写出去、方案落地都站不住脚。1.2 从“预估场景”到“不确定集”两类常用处理思路对比面对不确定性学术界和工程界的主流做法大致分成两类随机规划和鲁棒优化。随机规划的核心思路是给不确定性参数指定概率分布然后通过各种采样方法如蒙特卡洛模拟生成大量场景以期望成本最小化为目标进行优化。它的优点是结果相对“经济”不会过度保守缺点也很明显——需要知道准确的概率分布而现实中这个分布很难精确获得且计算量往往非常大。鲁棒优化则走了另一条路不再假设概率分布而是构造一个不确定性集合uncertainty set要求优化方案在这个集合中的任意一种情况下都可行且目标值尽量优。它不追求期望最优而是追求“最坏情况下的最优”也就是让系统在环境最恶劣时也不至于崩溃。两阶段鲁棒优化更是把这种思路往前推了一步第一阶段先做“现在必须定下来”的决策比如设备容量第二阶段等不确定参数实际暴露后再做“可以事后调整”的决策比如运行调度。这种“先决策、后调整”的架构非常契合微网电源容量配置加经济调度的实际流程。1.3 两阶段鲁棒优化为什么适合这个场景微网电源容量配置本质上是一个“投资决策 运行决策”耦合的问题。设备买多大、装多少属于长期投资决策一旦确定很难更改而系统怎么运行、储能怎么充放、燃气轮机出多少力属于短期运行决策可以根据实时风光出力动态调整。两阶段鲁棒优化天然匹配这个结构第一阶段确定各类电源的安装容量决策变量包括光伏装机容量、风电装机容量、储能额定容量与额定功率、燃气轮机装机容量等第二阶段在给定容量配置和实际风光出力的条件下最小化运行成本决策变量包括各机组出力、储能充放电功率、与主网交互功率等。同时两阶段结构允许不确定性在第二阶段“暴露”后再做运行调度这样求出来的容量配置既保证了系统在最恶劣风光出力下依然能安全运行又不会像单阶段鲁棒优化那样把所有决策都绑死导致结果过分保守。2. 两阶段鲁棒优化模型建模目标函数与约束条件拆解2.1 第一阶段决策变量设备选型与容量两阶段鲁棒优化的第一阶段核心任务是为微网选择合理的电源组成和容量配置。在我的仿真模型里候选设备包括光伏机组PV、风力发电机组WT、储能系统BESS和微型燃气轮机MT。第一阶段决策变量定义如下变量含义单位(C_{PV})光伏安装容量kW(C_{WT})风电安装容量kW(C_{BESS}^{E})储能额定容量kWh(C_{BESS}^{P})储能额定功率kW(C_{MT})燃气轮机安装容量kW这些变量的取值范围受限于微网建设场地的物理条件比如屋顶面积限制了光伏的最大安装容量风速资源决定了风电的上限这些会在模型里以容量上限约束的形式体现。第一阶段目标函数是所有设备的年化投资成本[ F_1 \sum_{i} C_{inv,i} \cdot C_i \cdot \frac{r(1r)^Y}{(1r)^Y - 1} ]其中 (r) 是折现率(Y) 是设备使用寿命(C_{inv,i}) 是单位容量投资成本。这个年化处理非常重要它把一次性投资摊到每一年让投资成本和运行成本可以在同一个时间尺度上相加。2.2 第二阶段决策变量运行调度第二阶段是在给定第一阶段容量配置方案 (C)、且实际风光出力 (P_{PV}^{actual}, P_{WT}^{actual}) 和负荷 (P_L) 已知的条件下制定系统的最优运行方案。第二阶段决策变量包括微型燃气轮机出力 (P_{MT,t})储能充电功率 (P_{ch,t})、放电功率 (P_{dis,t})储能荷电状态 (SOC_t)向主网购电功率 (P_{buy,t})、售电功率 (P_{sell,t})切负荷功率 (P_{curt,t})在极端场景下允许少量切负荷弃风弃光功率 (P_{abandon,t})。第二阶段目标函数是系统在一个典型日内的运行成本最小化[ F_2 \sum_{t} \left( c_{gas} \cdot \frac{P_{MT,t}}{\eta_{MT}} c_{buy} \cdot P_{buy,t} - c_{sell} \cdot P_{sell,t} c_{curt} \cdot P_{curt,t} \right) ]其中 (c_{gas}) 是天然气单价(\eta_{MT}) 是燃气轮机的发电效率(c_{buy}) 和 (c_{sell}) 分别是购电和售电电价(c_{curt}) 是切负荷惩罚系数。2.3 目标函数与约束的数学表达整体两阶段鲁棒优化模型可以写成如下紧凑形式[ \min_{C} ; \left( F_1(C) \max_{u \in U} \min_{y} ; F_2(y, u, C) \right) ]其中外层 (\min) 对第一阶段容量变量 (C) 进行优化内层 (\max) 寻找最恶劣的不确定性场景 (u)最内层 (\min) 在给定容量和场景下求最优运行成本。约束条件包括功率平衡约束24小时逐时满足[ P_{PV,t} P_{WT,t} P_{MT,t} P_{dis,t} P_{buy,t} P_{L,t} P_{ch,t} P_{sell,t} ]储能约束[ SOC_{t1} SOC_t \eta_{ch} \cdot P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}} ][ 0 \le SOC_t \le C_{BESS}^{E}, \quad 0 \le P_{ch,t}, P_{dis,t} \le C_{BESS}^{P} ]燃气轮机出力约束[ 0 \le P_{MT,t} \le C_{MT} ]主网交互功率约束[ 0 \le P_{buy,t} \le P_{buy}^{max}, \quad 0 \le P_{sell,t} \le P_{sell}^{max} ]风光出力约束[ 0 \le P_{PV,t} \le P_{PV,t}^{available}, \quad 0 \le P_{WT,t} \le P_{WT,t}^{available} ]2.4 不确定集的设计箱式与预算约束两阶段鲁棒优化的核心在于不确定性集合 (U) 的构造。为了兼顾计算复杂度和解的保守程度我采用箱式不确定集加预算约束的形式。设想我们要处理光伏出力和风电出力两方面的不确定性。设某时段风光出力的预测值为 (P_{PV,t}^{forecast}) 和 (P_{WT,t}^{forecast})实际取值范围为[ P_{PV,t}^{actual} \in \left[ (1 - \alpha_{PV}) P_{PV,t}^{forecast}, ; (1 \alpha_{PV}) P_{PV,t}^{forecast} \right] ][ P_{WT,t}^{actual} \in \left[ (1 - \alpha_{WT}) P_{WT,t}^{forecast}, ; (1 \alpha_{WT}) P_{WT,t}^{forecast} \right] ]其中 (\alpha_{PV}) 和 (\alpha_{WT}) 表示预测偏差比例。但这里有个问题如果允许每一个时段的不确定量都取最恶劣值求解出来的方案会过于保守因为现实中不可能全天候24小时都处于“最差风光”状态。为了控制保守度引入预算约束[ \sum_{t} z_{PV,t} \sum_{t} z_{WT,t} \le \Gamma ]其中 (z_{PV,t}) 和 (z_{WT,t}) 是0-1变量或连续辅助变量表示该时段是否达到最恶劣偏差。(\Gamma) 是预算参数物理含义是“最坏情况最多在多少个时段同时发生”。(\Gamma) 越大鲁棒性越强经济性也越差(\Gamma 0) 时退化为确定性模型。我仿真时常用的做法是(\Gamma) 取调度周期内总时段的1/3左右既保证鲁棒性又不至于让成本涨得太离谱。3. CCG求解算法从max-min结构到可计算形式3.1 模型整体结构MP与SP两阶段鲁棒优化模型直接求解是不现实的因为内层有一个“(\max \min)”嵌套结构。处理这种问题的主流方法是Benders对偶切平面法和CCG列与约束生成Column-and-Constraint Generation算法。我实际使用的是CCG它在处理混合整数鲁棒优化问题时比Benders法收敛更快尤其在场景数较多的情况下优势明显。CCG将原问题分解为主问题Master ProblemMP和子问题SubproblemSP主问题MP在第一阶段变量 (C) 的基础上加入一系列已生成的最恶劣场景对应的第二阶段变量和约束求解一个确定性的优化问题子问题SP固定第一阶段变量 (C)在不确定集 (U) 中寻找使运行成本最大化的场景 (u^*)这个场景称为“最恶劣场景”。两者迭代求解子问题找到新的最恶劣场景后将其对应的第二阶段变量和约束加入主问题重新求解主问题获得新的容量配置如此循环直到主问题目标值和子问题目标值之间的间隙小于预设容忍度。3.2 子问题的对偶转化子问题本质上是[ SP: \quad Q(C^) \max_{u \in U} \min_{y} ; F_2(y, u, C^) ]由于第二阶段是线性规划假设所有变量连续根据强对偶定理可以把内层 (\min_y) 转化为对偶最大化问题从而把子问题改写成单层最大化问题[ SP: \quad Q(C^) \max_{u \in U, \lambda} ; \lambda^T \cdot g(u, C^) ]这里的 (\lambda) 是对偶变量(g(u, C^*)) 是包含不确定参数 (u) 的约束右端项。转化后的单层问题是一个带有双线性项对偶变量乘不确定变量的优化问题处理方法是采用大M法引入辅助变量进行线性化。这在实际建模中非常繁琐但 YALMIP 工具箱提供了便捷方式可以直接用dualize函数处理对偶问题大大降低了建模难度。我实测下来YALMIP 的dualize在处理中小规模对偶问题时比较稳定但一旦约束规模比较大生成的对偶模型会有冗余变量偶尔也会出问题。因此我更推荐一种变通手段手动推导对偶问题。虽然麻烦但可控性强调试起来也容易定位错误。3.3 CCG迭代流程与收敛判据CCG 的完整迭代步骤如下初始化设置下界 (LB -\infty)上界 (UB \infty)迭代次数 (k 1)求解主问题此时包含已生成的最恶劣场景集合 ({u_1^, u_2^, ..., u_{k-1}^*})得到最优解 ((C_k, \eta_k))更新下界 (LB \eta_k)固定 (C C_k)求解子问题得到当前容量下的最恶劣场景 (u_k^*) 和最恶劣运行成本 (Q(C_k))更新上界 (UB \min(UB, F_1(C_k) Q(C_k)))若 ((UB - LB) / UB \le \epsilon)(\epsilon) 通常取0.01或0.001则算法收敛输出当前最优容量配置否则将新找到的最恶劣场景 (u_k^*) 对应的第二阶段变量 (y_k) 和约束加入主问题令 (k k 1)回到第2步。收敛判据这里有个细节容易被忽略主问题更新下界时取的是主问题的目标值包括投资成本和已生成场景的运行成本而上界是“当前投资成本 最恶劣场景运行成本”的较小者。这个关系必须严格对应否则可能出现上下界不收敛、循环无法退出的情况。我在仿真中发现CCG算法通常在5到8次迭代内就能收敛到1%的间隙以内速度相当快。个别情况下如果子问题求解不精确会出现震荡这时候需要检查子问题求解器的容忍度设置。4. MATLABYALMIPCPLEX环境搭建与工具选型原因4.1 为什么是YALMIPCPLEX而不是其他组合做优化仿真工具链的选择直接影响开发效率和求解性能。我最终选定 MATLAB YALMIP CPLEX是经过反复对比后确定的组合。YALMIP是一个MATLAB下的建模工具箱它最大的价值是让你用接近自然语言的方式描述优化问题而不用关心底层的求解器调用细节。比如你要定义一个变量、一组约束、一个目标函数YALMIP的语法几乎是“所见即所得”。CPLEX则是目前公认的线性规划、混合整数线性规划求解器里性能最稳的之一。对于两阶段鲁棒优化这种需要反复迭代求解的模型每一轮迭代都要调用CPLEX求解主问题和子问题求解速度直接决定你能否在可接受的时间内完成仿真。我之前用开源求解器做过对比同样的模型规模CPLEX的求解速度能快3到5倍。这一点在迭代几十轮的场景下非常关键。也有人会考虑用Gurobi替代CPLEX两者都是顶级求解器。我选CPLEX主要是学术圈内的代码资料更多遇到问题更容易找到参考。4.2 环境配置步骤与常见坑环境搭建本身不复杂但有几个坑非常折磨人。第一步安装MATLAB并确认版本建议使用R2018b及以上版本。YALMIP持续更新旧版MATLAB在部分新语法支持上会出问题。第二步下载YALMIP工具箱从YALMIP官网或GitHub仓库下载最新版压缩包解压后将整个文件夹路径添加到MATLAB路径中。操作步骤是主页 - 设置路径 - 添加文件夹 - 选择YALMIP根目录。第三步安装CPLEX并配置接口CPLEX需要去IBM官网注册下载选择与MATLAB版本匹配的版本。安装完成后在MATLAB里运行addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab\x64_win64); savepath;路径根据实际安装位置调整。我用的是CPLEX Studio 12.10或22.1版和R2020a、R2021b搭配都没问题。第四步验证YALMIP能否正确调用CPLEX运行以下测试代码sdpvar x y; Constraints [x y 1, x 0, y 0]; Objective -x - y; options sdpsettings(solver, cplex); optimize(Constraints, Objective, options);如果返回值problem为0说明CPLEX配置成功。这里有一个非常常见的坑YALMIP在调用CPLEX时如果cplex的许可证书未激活或者环境变量没配置好MATLAB会报“could not find solver”的错误。解决办法是检查CPLEX安装目录下的许可文件以及MATLAB的setenv环境变量设置。第五步配置CPLEX的参数在求解较大规模模型时CPLEX默认参数未必是最优的。我在仿真里经常调整的参数有options sdpsettings(solver, cplex, ... cplex.mip.tolerances.mipgap, 0.001, ... % 混合整数规划的相对间隙 cplex.mip.tolerances.integrality, 1e-9, ...% 整数变量的整性容忍度 cplex.mip.strategy.search, 1, ... % 分支策略 cplex.parallel, 4); % 并行核心数这些参数的合理设置能显著减少单次求解时间。我之前用默认参数跑一个包含数千个约束的主问题单次求解要70多秒调整参数后降到20秒左右。5. 核心代码实现主问题、子问题与迭代框架5.1 数据与参数初始化先定义算例的基础参数。以一个包含光伏、风电、储能、燃气轮机的小型微网为例%% 基础参数 T 24; % 调度周期24小时 dt 1; % 时间步长1小时 % 负荷数据kW P_load [420, 400, 380, 350, 360, 380, 450, 520, 600, 680, ... 720, 750, 730, 700, 680, 690, 720, 760, 820, 860, ... 800, 720, 620, 500]; % 光伏预测出力p.u.基于安装容量的归一化值 PV_forecast [0, 0, 0, 0, 0.05, 0.15, 0.35, 0.55, 0.75, 0.88, ... 0.95, 0.98, 0.96, 0.90, 0.78, 0.60, 0.40, 0.22, ... 0.08, 0, 0, 0, 0, 0]; % 风电预测出力p.u. WT_forecast [0.52, 0.48, 0.42, 0.38, 0.35, 0.32, 0.30, 0.28, ... 0.25, 0.27, 0.30, 0.34, 0.38, 0.42, 0.45, 0.48, ... 0.52, 0.55, 0.58, 0.60, 0.58, 0.55, 0.53, 0.50]; % 经济参数 c_buy 0.8; % 购电价元/kWh c_sell 0.4; % 售电价元/kWh c_gas 0.35; % 天然气单价元/kWh eta_MT 0.38; % 燃气轮机发电效率 eta_ch 0.95; % 储能充电效率 eta_dis 0.95; % 储能放电效率 % 不确定度参数后续会以这些预测值为中心生成不确定集 alpha_PV 0.2; % 光伏出力偏差比例 alpha_WT 0.15; % 风电出力偏差比例 Gamma 8; % 预算参数这里有一点要注意光伏归一化出力值是根据典型日光照曲线得到的实际使用时建议用本地气象数据统计生成这样仿真结果更有说服力。5.2 主问题的YALMIP建模主问题是CCG迭代框架的核心部分。初始迭代时主问题里只有一个名义场景随着迭代进行新的最恶劣场景会不断加入。%% 主问题建模函数 function [C_opt, LB, details] solve_master_problem(worst_scenarios, params) % 解包参数 T params.T; C_PV_max params.C_PV_max; C_WT_max params.C_WT_max; ... % 第一阶段变量容量决策 C_PV sdpvar(1, 1); C_WT sdpvar(1, 1); C_bess_E sdpvar(1, 1); C_bess_P sdpvar(1, 1); C_MT sdpvar(1, 1); % 引入辅助变量 eta表示当前主问题的总成本上界 eta sdpvar(1, 1); % 目标函数投资成本 运行成本通过eta传递 InvestCost inv_cost_PV * C_PV inv_cost_WT * C_WT ... inv_cost_BESS_E * C_bess_E inv_cost_BESS_P * C_bess_P ... inv_cost_MT * C_MT; Objective InvestCost eta; % 容量上限约束 Constraints [C_PV_max C_PV 0, ... C_WT_max C_WT 0, ... C_bess_E_max C_bess_E 0, ... C_bess_P_max C_bess_P 0, ... C_MT_max C_MT 0]; % 遍历所有已生成的最恶劣场景加入对应的第二阶段约束 scenario_count size(worst_scenarios, 2); for k 1:scenario_count P_pv_worst worst_scenarios(k).P_pv; P_wt_worst worst_scenarios(k).P_wt; % 第二阶段变量当前场景 P_MT sdpvar(T, 1); P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); SOC sdpvar(T1, 1); P_buy sdpvar(T, 1); P_sell sdpvar(T, 1); P_curt sdpvar(T, 1); % 功率平衡约束 Constraints [Constraints, ... P_pv_worst P_wt_worst P_MT P_dis P_buy ... P_load P_ch P_sell P_curt]; % 储能动态约束 Constraints [Constraints, ... SOC(2:T1) SOC(1:T) eta_ch * P_ch - P_dis / eta_dis]; Constraints [Constraints, ... SOC(1) params.SOC_init, SOC(T1) params.SOC_init]; Constraints [Constraints, ... 0 SOC(1:T1) C_bess_E, ... 0 P_ch C_bess_P, ... 0 P_dis C_bess_P]; % 燃气轮机出力约束 Constraints [Constraints, 0 P_MT C_MT]; % 主网交互约束 Constraints [Constraints, ... 0 P_buy params.P_buy_max, ... 0 P_sell params.P_sell_max]; % 运行成本约束通过eta夹逼 RunCost sum(c_gas * P_MT / eta_MT c_buy * P_buy - ... c_sell * P_sell c_curt * P_curt); Constraints [Constraints, eta RunCost]; end % 求解 options sdpsettings(solver, cplex, verbose, 0, ... cplex.mip.tolerances.mipgap, 0.001); optimize(Constraints, Objective, options); % 输出结果 C_opt.C_PV value(C_PV); C_opt.C_WT value(C_WT); C_opt.C_bess_E value(C_bess_E); C_opt.C_bess_P value(C_bess_P); C_opt.C_MT value(C_MT); LB value(Objective); details.solver_output ... end这里有个YALMIP使用上的细节沿着场景循环添加变量和约束时每次循环都要新建sdpvar变量不能把不同场景的第二阶段变量共用同一个变量名否则YALMIP会将其视为同一个变量导致场景之间的约束耦合在一起模型就完全错了。我在第一次写这个主问题时就在这里栽过跟头当时所有场景共用一组P_MT、P_ch变量结果求解出的容量配置明显偏大检查半天才发现是变量名作用域出了问题。5.3 子问题对偶建模与最坏场景求解子问题的目标是在给定容量配置 (C^*) 的条件下找到使运行成本最大的风光出力场景。%% 子问题求解函数寻找最恶劣场景 function [Q_val, worst_scenario] solve_subproblem(C_current, params) % 解包当前容量配置 C_PV C_current.C_PV; C_WT C_current.C_WT; C_bess_E C_current.C_bess_E; C_bess_P C_current.C_bess_P; C_MT C_current.C_MT; T params.T; % 不确定变量各时段实际风光出力 P_pv_actual sdpvar(T, 1); P_wt_actual sdpvar(T, 1); % 第二阶段运行变量 P_MT sdpvar(T, 1); P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); SOC sdpvar(T1, 1); P_buy sdpvar(T, 1); P_sell sdpvar(T, 1); P_curt sdpvar(T, 1); % 不确定集约束箱式约束加预算约束 z_pv binvar(T, 1); % 光伏偏差指示变量 z_wt binvar(T, 1); % 风电偏差指示变量 Constraints []; for t 1:T % 光伏出力不确定范围 Constraints [Constraints, ... P_pv_actual(t) params.PV_forecast(t) * C_PV * (1 - alpha_PV 2 * alpha_PV * z_pv(t))]; % 风电出力不确定范围 Constraints [Constraints, ... P_wt_actual(t) params.WT_forecast(t) * C_WT * (1 - alpha_WT 2 * alpha_WT * z_wt(t))]; end % 预算约束 Constraints [Constraints, sum(z_pv) sum(z_wt) Gamma]; % 第二阶段运行约束与主问题类似 ... % 子问题目标最大化运行成本 RunCost sum(c_gas * P_MT / eta_MT c_buy * P_buy - ... c_sell * P_sell c_curt * P_curt); Objective -RunCost; % 注意max 转 min % 求解 options sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, Objective, options); Q_val value(RunCost); worst_scenario.P_pv value(P_pv_actual); worst_scenario.P_wt value(P_wt_actual); end这里采用的是“直接构造原始子问题”的写法利用YALMIP可以对整数变量直接建模的能力来处理不确定集中的预算约束。这种做法的好处是直观、不容易出错代价是求解速度比纯线性对偶问题要慢一些。如果模型规模特别大建议还是走对偶加线性化的路。5.4 CCG循环控制主程序和子程序都准备好后把它们串起来的迭代框架如下%% CCG主循环 function result CCG_main(params) % 初始化 UB inf; LB -inf; iter 0; max_iter 20; epsilon 0.01; % 初始场景名义预测场景 worst_scenarios struct([]); worst_scenarios(1).P_pv params.PV_forecast .* params.C_PV_init; worst_scenarios(1).P_wt params.WT_forecast .* params.C_WT_init; while iter max_iter iter iter 1; fprintf(迭代次数: %d\n, iter); % 1. 求解主问题得到容量配置和新的下界 [C_current, LB_new, mp_detail] solve_master_problem(worst_scenarios, params); LB LB_new; % 2. 固定当前容量配置求解子问题得到最恶劣场景 [Q_val, worst_scenario_new] solve_subproblem(C_current, params); % 3. 计算当前上界投资成本 最恶劣场景运行成本 InvestCost inv_cost_PV * C_current.C_PV inv_cost_WT * C_current.C_WT ... inv_cost_BESS_E * C_current.C_bess_E inv_cost_BESS_P * C_current.C_bess_P ... inv_cost_MT * C_current.C_MT; UB min(UB, InvestCost Q_val); fprintf(LB %.4f, UB %.4f, gap %.4f%%\n, LB, UB, abs((UB - LB) / UB * 100)); % 4. 判断收敛 if abs((UB - LB) / UB) epsilon fprintf(算法收敛于第 %d 次迭代\n, iter); result.C C_current; result.UB UB; result.LB LB; result.iter iter; result.worst_scenario worst_scenario_new; return; end % 5. 将新场景加入主问题场景集合 if isempty(worst_scenarios) worst_scenarios(1) worst_scenario_new; else worst_scenarios(end1) worst_scenario_new; end end error(达到最大迭代次数算法未收敛); end这里要特别强调一个细节每次迭代后保存的新场景 (u_k^*) 是子问题在当前容量配置 (C_k) 下求出的“最恶劣场景”。把这个场景加入主问题本质上是告诉主问题“你算出来的这套容量配置在遇到这个场景时可能要付出很高的运行成本你重新优化一下。”这个过程反复进行容量配置方案会逐步趋向于“在所有可能的最恶劣场景下都可接受”的鲁棒方案。6. 仿真算例与结果分析配置方案与算法性能6.1 算例参数为了验证模型和算法的有效性我设计了一个典型日微网算例。微网包含光伏、风电、储能、燃气轮机并与配电网保持连接。主要参数如下表参数数值单位微网峰值负荷860kW光伏最大安装容量1200kW风电最大安装容量600kW储能最大额定容量1000kWh储能最大额定功率500kW燃气轮机最大安装容量800kW光伏单位投资成本3500元/kW风电单位投资成本6000元/kW储能单位容量成本1500元/kWh燃气轮机单位投资成本3000元/kW折现率0.06-设备寿命20年6.2 迭代收敛与最坏场景运行CCG主循环后典型的收敛过程如下迭代次数下界LB万元/年上界UB万元/年间隙%1452.3489.67.62468.5485.23.43474.8481.31.44476.2479.50.75477.0478.40.3算法在5次迭代内收敛到0.3%的间隙效果相当理想。这说明CCG处理这个规模的问题收敛速度很快不需要担心计算负担。最恶劣场景揭示了一个重要现象出现最恶劣组合的时段集中在傍晚17:00到21:00——此时光伏出力变为0风电出力同时处于低谷负荷却达到全天最高。这个时段系统几乎完全依赖储能放电和燃气轮机出力加上从主网购电才能维持供电。这个结论很有价值它告诉我们在容量配置时要把“晚高峰时段”作为最关键的校核场景。6.3 配置结果对比鲁棒解vs确定性解为了展示两阶段鲁棒优化的效果我把鲁棒模型的结果与确定性模型进行对比。配置方案光伏容量kW风电容量kW储能容量kWh储能功率kW燃气轮机kW年总成本万元确定性解900420650280300458.6鲁棒解Γ8820380780320420478.3可以看到鲁棒解相比确定性解光伏和风电的安装容量有所降低储能和燃气轮机的容量有所增加年总成本上升了约4.3%。这个结果从工程角度非常好理解光伏和风电是不确定电源过多依赖它们会在恶劣场景下造成供电缺口因此鲁棒解主动降低了这两个的配置比例储能和燃气轮机是可控电源它们才是应对不确定性的“压舱石”因此鲁棒解增加了这两个的容量。多花的4.3%投资换来的是系统在极端天气条件下的可靠供电能力这笔账在经济性分析中需要结合停电损失一起评估。我在论文中通常会加一组对比如果把切负荷惩罚成本拉高模拟重要负荷场景鲁棒解的经济优势会立刻体现出来。6.4 保守度参数对结果的影响预算参数 (\Gamma) 直接控制鲁棒模型的保守程度。我扫描了不同的 (\Gamma) 值观察它对成本的影响。Γ值年总成本万元相比Γ0增幅%最恶劣场景运行成本万元/天0确定性458.6-0.854466.21.70.928478.34.31.0412489.56.71.1516497.88.51.2424503.29.71.30可以看到随着 (\Gamma) 增大总成本单调上升但增幅逐渐放缓。这说明鲁棒性带来的成本增加呈现边际递减效应当 (\Gamma) 达到一定值后继续增加对系统鲁棒性的提升已不明显但对经济性的损害却在累积。实际工程中(\Gamma) 不应一味取大。我建议的做法是基于历史气象数据统计一个观察周期内风光出力的最大偏差时长占比以这个统计值作为 (\Gamma) 的选取参考这样既避免了拍脑袋定参数又保证了模型的鲁棒性贴合实际需求。7. 实操中容易踩的坑与调试建议7.1 YALMIP求解器配置问题YALMIP 在使用中最大的坑是求解器调用不明确。当同时安装了多个求解器时比如既装了CPLEX又装了GurobiYALMIP有可能会自动选择别的求解器导致结果异常或者性能骤降。我在仿真中遇到过一个问题明明安装好了CPLEX求解时却提示“No suitable solver”。排查后发现是因为YALMIP版本过旧无法识别新版本CPLEX的接口。解决办法是升级YALMIP到最新版并确保CPLEX的MATLAB接口路径已正确添加。建议在每次调用前显式指定求解器options sdpsettings(solver, cplex, verbose, 1);不要依赖YALMIP的默认求解器选择机制。7.2 CPLEX整数规划求解缓慢在子问题中预算约束使用了0-1整数变量 (z_{PV,t}) 和 (z_{WT,t})这意味着子问题是一个混合整数规划。当 (T) 增大到96即15分钟一个间隔或者不确定变量的时段数变多时求解时间会指数级增长。我最初的仿真中T24时子问题求解只需1到2秒但当我把时间分辨率提高到15分钟T96后单次子问题求解时间飙升到接近60秒整个CCG流程跑完需要近10分钟。解决思路有两个方向一是对偶线性化。将子问题转化为纯线性规划去掉整数变量这样求解速度会快很多。具体做法是引入不确定变量的连续辅助变量用KKT条件或强对偶理论处理内层min问题。二是场景缩减。对历史数据进行聚类分析将相似的时段合并减少需要单独建模的时段数。我实际采用的是第一种方法通过对偶转化将子问题中的整数变量替换为连续变量求解时间降到了5秒以内。7.3 对偶问题的正确性验证手动推导对偶问题时最容易犯的错误是对偶方向写反、约束符号搞错、对偶变量漏项。这些问题不会导致模型无解但会给出一个错误的最恶劣场景进而污染整个迭代过程。验证对偶问题是否正确有一个非常实用的方法取一组相同的参数分别求解原始子问题和对偶问题对比两者的目标函数值。由于强对偶定理保证两者相等如果数值不一致说明对偶推导或实现有问题。我在代码里写了这样一个自检函数function verify_dual(C_current, params) % 求解原始子问题 Q_primal solve_subproblem_primal(C_current, params); % 求解对偶子问题 Q_dual solve_subproblem_dual(C_current, params); % 对比 fprintf(原始问题目标值: %.4f\n, Q_primal); fprintf(对偶问题目标值: %.4f\n, Q_dual); if abs(Q_primal - Q_dual) 1e-6 error(对偶验证失败); end end每次都跑一遍自检发现问题即刻排查能节省大量后期调试时间。7.4 大M法与数值稳定性在使用大M法处理双线性项时M值的选取直接影响求解稳定性。M值过小会导致不可行解M值过大又可能引起数值问题导致CPLEX出现严重病态矩阵警告。我在调试时发现M值并非越大越好。CPLEX官方文档里建议M的取值应该比模型中的其他系数大一到两个数量级即可不要盲目取1e9这类超大值。一般我会根据模型参数估算M的范围功率平衡约束中的M值取最大可能的功率缺额如 (M \ge \max(P_{load}) \max(P_{ch}) \max(P_{sell}))容量上限约束中的M值取对应容量的上限值即可。此外在YALMIP中使用大M法时尽量把M写成具体的数值而不是一个很大的1e6常量这样CPLEX在做预求解时能够更有效地收紧边界提升求解速度。8. 从仿真到论文结果呈现与延伸方向8.1 仿真结果的可视化呈现仿真做完后出图是必不可少的一步。我在论文中一般会准备四类图第一类是迭代收敛曲线图展示CCG算法中上下界随迭代次数的变化过程证明算法的收敛性。可以生成这样的图figure; plot(1:iter, LB_history, -o, LineWidth, 1.5); hold on; plot(1:iter, UB_history, -s, LineWidth, 1.5); xlabel(迭代次数); ylabel(成本万元); legend(下界LB, 上界UB, Location, northeast); grid on; saveas(gcf, convergence_curve.png);第二类是容量配置对比图用堆叠柱状图展示确定性模型和鲁棒模型在不同电源上的容量配置差异。第三类是最恶劣场景下的功率平衡图展示在最恶劣场景下各电源的出力分配情况和系统整体功率平衡关系。这张图最能直观展示鲁棒方案在极端工况下的应对能力。第四类是用折线图展示不同 (\Gamma) 值下年总成本和切负荷量的变化趋势直观呈现鲁棒性与经济性的折中关系。8.2 两阶段鲁棒优化模型的扩展方向两阶段鲁棒优化在微网电源容量配置中的应用远不止我上面实现的这些内容。如果你做完基础模型想往更深的方向走有几个我认为很有价值的延伸方向第一多微网互联场景。多个微网之间通过联络线互济可以大幅提高系统整体的鲁棒性降低成本。但这种模型会引入更多的一阶段决策变量计算复杂度急剧上升需要配合分解算法或分布式求解框架。第二考虑需求响应的鲁棒配置。把可调负荷纳入优化框架让负荷侧也参与削峰填谷。需求响应本质上是给系统增加了柔性资源在不确定环境下价值更为明显。第三多阶段鲁棒优化。两阶段模型适合“先配置、后运行”的场景但实际运行中调度决策往往是滚动时域的。多阶段鲁棒模型能更精确地刻画这种时序决策关系但建模难度和求解复杂度都会显著上升。第四将鲁棒优化与机器学习结合。用LSTM或Transformer对风光出力进行概率预测将预测区间作为不确定集的边界这样不确定集不再是人为设定而是由数据驱动生成模型的实用性会更强。8.3 一些个人心得做这个课题前前后后大概花了两个月。最大的感受是两阶段鲁棒优化本身并不难理解难的是把模型转成代码、并让代码在各种边界条件下都能稳定运行。数学推导和代码实现之间存在一条需要反复跳转的鸿沟跳过去了前面的路就顺了。如果你刚开始接触这个方向我的建议是先不管YALMIP手推一个最简单的两阶段鲁棒模型比如单时段、单设备然后手动用对偶转换加CCG求解一遍完全吃透整个流程。再回到MATLAB里用YALMIP写代码你会觉得一切都是顺理成章的。工具是辅助对模型本质的理解才是决定你能否做深做透的关键。 两阶段鲁棒优化在微网电源容量配置中的应用与实现
开题那阵子我琢磨微网电源容量配置满脑子都是“怎么让这套系统在现实里扛得住”。风光出力的随机性、负荷的波动任何一个没照顾好之前算的“最优方案”到实际运行就可能变成“最贵方案”甚至“失稳方案”。后来我把目光锁定在两阶段鲁棒优化算法上配合MATLAB、YALMIP和CPLEX这条经典工具链才真正把这个问题从理论推到了可仿真的层面。这篇东西我想把我从建模到代码实现、再到踩坑排错的完整过程拆开讲清楚给正在做相关研究的同学一条能直接上手的路。1. 为什么容量配置必须考虑不确定性确定性模型的局限1.1 传统确定性优化的问题描述与短板在微网电源容量优化里传统做法通常是把风光出力、负荷大小当作已知量来处理。设定几个典型日场景比如“夏季典型日”“冬季典型日”把每个时段的风速、光照、负荷都取一个固定值然后建立优化模型求解出光伏、风电、储能、燃气轮机这些设备的最优安装容量。这类确定性模型的数学形式很简洁[ \min_{C, P} ; C_{inv}(C) C_{ope}(P) ]其中 (C_{inv}) 是设备投资成本(C_{ope}) 是运行成本约束条件包括功率平衡、设备出力上下限、储能充放电约束等。所有参数都是确定值模型直接扔给求解器就能出结果。但这个思路有个致命的隐含假设未来风光出力和负荷曲线就是那几个典型日场景。实际情况呢光照可能连续三天低于预期风速可能整月都处于“要死不活”的状态负荷也可能因为极端天气而飙升。一旦实际运行与预设场景偏差过大之前算好的容量配置就会出现两类问题供电可靠性下降切负荷风险增加运行成本大幅偏离预期甚至出现必须高价购电的局面。我最早用确定性模型做仿真时结果很漂亮总成本最低、新能源渗透率最高。但一换到真实历史数据里回放最恶劣的那几天系统几乎撑不住。这个问题不解决论文写出去、方案落地都站不住脚。1.2 从“预估场景”到“不确定集”两类常用处理思路对比面对不确定性学术界和工程界的主流做法大致分成两类随机规划和鲁棒优化。随机规划的核心思路是给不确定性参数指定概率分布然后通过各种采样方法如蒙特卡洛模拟生成大量场景以期望成本最小化为目标进行优化。它的优点是结果相对“经济”不会过度保守缺点也很明显——需要知道准确的概率分布而现实中这个分布很难精确获得且计算量往往非常大。鲁棒优化则走了另一条路不再假设概率分布而是构造一个不确定性集合uncertainty set要求优化方案在这个集合中的任意一种情况下都可行且目标值尽量优。它不追求期望最优而是追求“最坏情况下的最优”也就是让系统在环境最恶劣时也不至于崩溃。两阶段鲁棒优化更是把这种思路往前推了一步第一阶段先做“现在必须定下来”的决策比如设备容量第二阶段等不确定参数实际暴露后再做“可以事后调整”的决策比如运行调度。这种“先决策、后调整”的架构非常契合微网电源容量配置加经济调度的实际流程。1.3 两阶段鲁棒优化为什么适合这个场景微网电源容量配置本质上是一个“投资决策 运行决策”耦合的问题。设备买多大、装多少属于长期投资决策一旦确定很难更改而系统怎么运行、储能怎么充放、燃气轮机出多少力属于短期运行决策可以根据实时风光出力动态调整。两阶段鲁棒优化天然匹配这个结构第一阶段确定各类电源的安装容量决策变量包括光伏装机容量、风电装机容量、储能额定容量与额定功率、燃气轮机装机容量等第二阶段在给定容量配置和实际风光出力的条件下最小化运行成本决策变量包括各机组出力、储能充放电功率、与主网交互功率等。同时两阶段结构允许不确定性在第二阶段“暴露”后再做运行调度这样求出来的容量配置既保证了系统在最恶劣风光出力下依然能安全运行又不会像单阶段鲁棒优化那样把所有决策都绑死导致结果过分保守。2. 两阶段鲁棒优化模型建模目标函数与约束条件拆解2.1 第一阶段决策变量设备选型与容量两阶段鲁棒优化的第一阶段核心任务是为微网选择合理的电源组成和容量配置。在我的仿真模型里候选设备包括光伏机组PV、风力发电机组WT、储能系统BESS和微型燃气轮机MT。第一阶段决策变量定义如下变量含义单位(C_{PV})光伏安装容量kW(C_{WT})风电安装容量kW(C_{BESS}^{E})储能额定容量kWh(C_{BESS}^{P})储能额定功率kW(C_{MT})燃气轮机安装容量kW这些变量的取值范围受限于微网建设场地的物理条件比如屋顶面积限制了光伏的最大安装容量风速资源决定了风电的上限这些会在模型里以容量上限约束的形式体现。第一阶段目标函数是所有设备的年化投资成本[ F_1 \sum_{i} C_{inv,i} \cdot C_i \cdot \frac{r(1r)^Y}{(1r)^Y - 1} ]其中 (r) 是折现率(Y) 是设备使用寿命(C_{inv,i}) 是单位容量投资成本。这个年化处理非常重要它把一次性投资摊到每一年让投资成本和运行成本可以在同一个时间尺度上相加。2.2 第二阶段决策变量运行调度第二阶段是在给定第一阶段容量配置方案 (C)、且实际风光出力 (P_{PV}^{actual}, P_{WT}^{actual}) 和负荷 (P_L) 已知的条件下制定系统的最优运行方案。第二阶段决策变量包括微型燃气轮机出力 (P_{MT,t})储能充电功率 (P_{ch,t})、放电功率 (P_{dis,t})储能荷电状态 (SOC_t)向主网购电功率 (P_{buy,t})、售电功率 (P_{sell,t})切负荷功率 (P_{curt,t})在极端场景下允许少量切负荷弃风弃光功率 (P_{abandon,t})。第二阶段目标函数是系统在一个典型日内的运行成本最小化[ F_2 \sum_{t} \left( c_{gas} \cdot \frac{P_{MT,t}}{\eta_{MT}} c_{buy} \cdot P_{buy,t} - c_{sell} \cdot P_{sell,t} c_{curt} \cdot P_{curt,t} \right) ]其中 (c_{gas}) 是天然气单价(\eta_{MT}) 是燃气轮机的发电效率(c_{buy}) 和 (c_{sell}) 分别是购电和售电电价(c_{curt}) 是切负荷惩罚系数。2.3 目标函数与约束的数学表达整体两阶段鲁棒优化模型可以写成如下紧凑形式[ \min_{C} ; \left( F_1(C) \max_{u \in U} \min_{y} ; F_2(y, u, C) \right) ]其中外层 (\min) 对第一阶段容量变量 (C) 进行优化内层 (\max) 寻找最恶劣的不确定性场景 (u)最内层 (\min) 在给定容量和场景下求最优运行成本。约束条件包括功率平衡约束24小时逐时满足[ P_{PV,t} P_{WT,t} P_{MT,t} P_{dis,t} P_{buy,t} P_{L,t} P_{ch,t} P_{sell,t} ]储能约束[ SOC_{t1} SOC_t \eta_{ch} \cdot P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}} ][ 0 \le SOC_t \le C_{BESS}^{E}, \quad 0 \le P_{ch,t}, P_{dis,t} \le C_{BESS}^{P} ]燃气轮机出力约束[ 0 \le P_{MT,t} \le C_{MT} ]主网交互功率约束[ 0 \le P_{buy,t} \le P_{buy}^{max}, \quad 0 \le P_{sell,t} \le P_{sell}^{max} ]风光出力约束[ 0 \le P_{PV,t} \le P_{PV,t}^{available}, \quad 0 \le P_{WT,t} \le P_{WT,t}^{available} ]2.4 不确定集的设计箱式与预算约束两阶段鲁棒优化的核心在于不确定性集合 (U) 的构造。为了兼顾计算复杂度和解的保守程度我采用箱式不确定集加预算约束的形式。设想我们要处理光伏出力和风电出力两方面的不确定性。设某时段风光出力的预测值为 (P_{PV,t}^{forecast}) 和 (P_{WT,t}^{forecast})实际取值范围为[ P_{PV,t}^{actual} \in \left[ (1 - \alpha_{PV}) P_{PV,t}^{forecast}, ; (1 \alpha_{PV}) P_{PV,t}^{forecast} \right] ][ P_{WT,t}^{actual} \in \left[ (1 - \alpha_{WT}) P_{WT,t}^{forecast}, ; (1 \alpha_{WT}) P_{WT,t}^{forecast} \right] ]其中 (\alpha_{PV}) 和 (\alpha_{WT}) 表示预测偏差比例。但这里有个问题如果允许每一个时段的不确定量都取最恶劣值求解出来的方案会过于保守因为现实中不可能全天候24小时都处于“最差风光”状态。为了控制保守度引入预算约束[ \sum_{t} z_{PV,t} \sum_{t} z_{WT,t} \le \Gamma ]其中 (z_{PV,t}) 和 (z_{WT,t}) 是0-1变量或连续辅助变量表示该时段是否达到最恶劣偏差。(\Gamma) 是预算参数物理含义是“最坏情况最多在多少个时段同时发生”。(\Gamma) 越大鲁棒性越强经济性也越差(\Gamma 0) 时退化为确定性模型。我仿真时常用的做法是(\Gamma) 取调度周期内总时段的1/3左右既保证鲁棒性又不至于让成本涨得太离谱。3. CCG求解算法从max-min结构到可计算形式3.1 模型整体结构MP与SP两阶段鲁棒优化模型直接求解是不现实的因为内层有一个“(\max \min)”嵌套结构。处理这种问题的主流方法是Benders对偶切平面法和CCG列与约束生成Column-and-Constraint Generation算法。我实际使用的是CCG它在处理混合整数鲁棒优化问题时比Benders法收敛更快尤其在场景数较多的情况下优势明显。CCG将原问题分解为主问题Master ProblemMP和子问题SubproblemSP主问题MP在第一阶段变量 (C) 的基础上加入一系列已生成的最恶劣场景对应的第二阶段变量和约束求解一个确定性的优化问题子问题SP固定第一阶段变量 (C)在不确定集 (U) 中寻找使运行成本最大化的场景 (u^*)这个场景称为“最恶劣场景”。两者迭代求解子问题找到新的最恶劣场景后将其对应的第二阶段变量和约束加入主问题重新求解主问题获得新的容量配置如此循环直到主问题目标值和子问题目标值之间的间隙小于预设容忍度。3.2 子问题的对偶转化子问题本质上是[ SP: \quad Q(C^) \max_{u \in U} \min_{y} ; F_2(y, u, C^) ]由于第二阶段是线性规划假设所有变量连续根据强对偶定理可以把内层 (\min_y) 转化为对偶最大化问题从而把子问题改写成单层最大化问题[ SP: \quad Q(C^) \max_{u \in U, \lambda} ; \lambda^T \cdot g(u, C^) ]这里的 (\lambda) 是对偶变量(g(u, C^*)) 是包含不确定参数 (u) 的约束右端项。转化后的单层问题是一个带有双线性项对偶变量乘不确定变量的优化问题处理方法是采用大M法引入辅助变量进行线性化。这在实际建模中非常繁琐但 YALMIP 工具箱提供了便捷方式可以直接用dualize函数处理对偶问题大大降低了建模难度。我实测下来YALMIP 的dualize在处理中小规模对偶问题时比较稳定但一旦约束规模比较大生成的对偶模型会有冗余变量偶尔也会出问题。因此我更推荐一种变通手段手动推导对偶问题。虽然麻烦但可控性强调试起来也容易定位错误。3.3 CCG迭代流程与收敛判据CCG 的完整迭代步骤如下初始化设置下界 (LB -\infty)上界 (UB \infty)迭代次数 (k 1)求解主问题此时包含已生成的最恶劣场景集合 ({u_1^, u_2^, ..., u_{k-1}^*})得到最优解 ((C_k, \eta_k))更新下界 (LB \eta_k)固定 (C C_k)求解子问题得到当前容量下的最恶劣场景 (u_k^*) 和最恶劣运行成本 (Q(C_k))更新上界 (UB \min(UB, F_1(C_k) Q(C_k)))若 ((UB - LB) / UB \le \epsilon)(\epsilon) 通常取0.01或0.001则算法收敛输出当前最优容量配置否则将新找到的最恶劣场景 (u_k^*) 对应的第二阶段变量 (y_k) 和约束加入主问题令 (k k 1)回到第2步。收敛判据这里有个细节容易被忽略主问题更新下界时取的是主问题的目标值包括投资成本和已生成场景的运行成本而上界是“当前投资成本 最恶劣场景运行成本”的较小者。这个关系必须严格对应否则可能出现上下界不收敛、循环无法退出的情况。我在仿真中发现CCG算法通常在5到8次迭代内就能收敛到1%的间隙以内速度相当快。个别情况下如果子问题求解不精确会出现震荡这时候需要检查子问题求解器的容忍度设置。4. MATLABYALMIPCPLEX环境搭建与工具选型原因4.1 为什么是YALMIPCPLEX而不是其他组合做优化仿真工具链的选择直接影响开发效率和求解性能。我最终选定 MATLAB YALMIP CPLEX是经过反复对比后确定的组合。YALMIP是一个MATLAB下的建模工具箱它最大的价值是让你用接近自然语言的方式描述优化问题而不用关心底层的求解器调用细节。比如你要定义一个变量、一组约束、一个目标函数YALMIP的语法几乎是“所见即所得”。CPLEX则是目前公认的线性规划、混合整数线性规划求解器里性能最稳的之一。对于两阶段鲁棒优化这种需要反复迭代求解的模型每一轮迭代都要调用CPLEX求解主问题和子问题求解速度直接决定你能否在可接受的时间内完成仿真。我之前用开源求解器做过对比同样的模型规模CPLEX的求解速度能快3到5倍。这一点在迭代几十轮的场景下非常关键。也有人会考虑用Gurobi替代CPLEX两者都是顶级求解器。我选CPLEX主要是学术圈内的代码资料更多遇到问题更容易找到参考。4.2 环境配置步骤与常见坑环境搭建本身不复杂但有几个坑非常折磨人。第一步安装MATLAB并确认版本建议使用R2018b及以上版本。YALMIP持续更新旧版MATLAB在部分新语法支持上会出问题。第二步下载YALMIP工具箱从YALMIP官网或GitHub仓库下载最新版压缩包解压后将整个文件夹路径添加到MATLAB路径中。操作步骤是主页 - 设置路径 - 添加文件夹 - 选择YALMIP根目录。第三步安装CPLEX并配置接口CPLEX需要去IBM官网注册下载选择与MATLAB版本匹配的版本。安装完成后在MATLAB里运行addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab\x64_win64); savepath;路径根据实际安装位置调整。我用的是CPLEX Studio 12.10或22.1版和R2020a、R2021b搭配都没问题。第四步验证YALMIP能否正确调用CPLEX运行以下测试代码sdpvar x y; Constraints [x y 1, x 0, y 0]; Objective -x - y; options sdpsettings(solver, cplex); optimize(Constraints, Objective, options);如果返回值problem为0说明CPLEX配置成功。这里有一个非常常见的坑YALMIP在调用CPLEX时如果cplex的许可证书未激活或者环境变量没配置好MATLAB会报“could not find solver”的错误。解决办法是检查CPLEX安装目录下的许可文件以及MATLAB的setenv环境变量设置。第五步配置CPLEX的参数在求解较大规模模型时CPLEX默认参数未必是最优的。我在仿真里经常调整的参数有options sdpsettings(solver, cplex, ... cplex.mip.tolerances.mipgap, 0.001, ... % 混合整数规划的相对间隙 cplex.mip.tolerances.integrality, 1e-9, ...% 整数变量的整性容忍度 cplex.mip.strategy.search, 1, ... % 分支策略 cplex.parallel, 4); % 并行核心数这些参数的合理设置能显著减少单次求解时间。我之前用默认参数跑一个包含数千个约束的主问题单次求解要70多秒调整参数后降到20秒左右。5. 核心代码实现主问题、子问题与迭代框架5.1 数据与参数初始化先定义算例的基础参数。以一个包含光伏、风电、储能、燃气轮机的小型微网为例%% 基础参数 T 24; % 调度周期24小时 dt 1; % 时间步长1小时 % 负荷数据kW P_load [420, 400, 380, 350, 360, 380, 450, 520, 600, 680, ... 720, 750, 730, 700, 680, 690, 720, 760, 820, 860, ... 800, 720, 620, 500]; % 光伏预测出力p.u.基于安装容量的归一化值 PV_forecast [0, 0, 0, 0, 0.05, 0.15, 0.35, 0.55, 0.75, 0.88, ... 0.95, 0.98, 0.96, 0.90, 0.78, 0.60, 0.40, 0.22, ... 0.08, 0, 0, 0, 0, 0]; % 风电预测出力p.u. WT_forecast [0.52, 0.48, 0.42, 0.38, 0.35, 0.32, 0.30, 0.28, ... 0.25, 0.27, 0.30, 0.34, 0.38, 0.42, 0.45, 0.48, ... 0.52, 0.55, 0.58, 0.60, 0.58, 0.55, 0.53, 0.50]; % 经济参数 c_buy 0.8; % 购电价元/kWh c_sell 0.4; % 售电价元/kWh c_gas 0.35; % 天然气单价元/kWh eta_MT 0.38; % 燃气轮机发电效率 eta_ch 0.95; % 储能充电效率 eta_dis 0.95; % 储能放电效率 % 不确定度参数后续会以这些预测值为中心生成不确定集 alpha_PV 0.2; % 光伏出力偏差比例 alpha_WT 0.15; % 风电出力偏差比例 Gamma 8; % 预算参数这里有一点要注意光伏归一化出力值是根据典型日光照曲线得到的实际使用时建议用本地气象数据统计生成这样仿真结果更有说服力。5.2 主问题的YALMIP建模主问题是CCG迭代框架的核心部分。初始迭代时主问题里只有一个名义场景随着迭代进行新的最恶劣场景会不断加入。%% 主问题建模函数 function [C_opt, LB, details] solve_master_problem(worst_scenarios, params) % 解包参数 T params.T; C_PV_max params.C_PV_max; C_WT_max params.C_WT_max; ... % 第一阶段变量容量决策 C_PV sdpvar(1, 1); C_WT sdpvar(1, 1); C_bess_E sdpvar(1, 1); C_bess_P sdpvar(1, 1); C_MT sdpvar(1, 1); % 引入辅助变量 eta表示当前主问题的总成本上界 eta sdpvar(1, 1); % 目标函数投资成本 运行成本通过eta传递 InvestCost inv_cost_PV * C_PV inv_cost_WT * C_WT ... inv_cost_BESS_E * C_bess_E inv_cost_BESS_P * C_bess_P ... inv_cost_MT * C_MT; Objective InvestCost eta; % 容量上限约束 Constraints [C_PV_max C_PV 0, ... C_WT_max C_WT 0, ... C_bess_E_max C_bess_E 0, ... C_bess_P_max C_bess_P 0, ... C_MT_max C_MT 0]; % 遍历所有已生成的最恶劣场景加入对应的第二阶段约束 scenario_count size(worst_scenarios, 2); for k 1:scenario_count P_pv_worst worst_scenarios(k).P_pv; P_wt_worst worst_scenarios(k).P_wt; % 第二阶段变量当前场景 P_MT sdpvar(T, 1); P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); SOC sdpvar(T1, 1); P_buy sdpvar(T, 1); P_sell sdpvar(T, 1); P_curt sdpvar(T, 1); % 功率平衡约束 Constraints [Constraints, ... P_pv_worst P_wt_worst P_MT P_dis P_buy ... P_load P_ch P_sell P_curt]; % 储能动态约束 Constraints [Constraints, ... SOC(2:T1) SOC(1:T) eta_ch * P_ch - P_dis / eta_dis]; Constraints [Constraints, ... SOC(1) params.SOC_init, SOC(T1) params.SOC_init]; Constraints [Constraints, ... 0 SOC(1:T1) C_bess_E, ... 0 P_ch C_bess_P, ... 0 P_dis C_bess_P]; % 燃气轮机出力约束 Constraints [Constraints, 0 P_MT C_MT]; % 主网交互约束 Constraints [Constraints, ... 0 P_buy params.P_buy_max, ... 0 P_sell params.P_sell_max]; % 运行成本约束通过eta夹逼 RunCost sum(c_gas * P_MT / eta_MT c_buy * P_buy - ... c_sell * P_sell c_curt * P_curt); Constraints [Constraints, eta RunCost]; end % 求解 options sdpsettings(solver, cplex, verbose, 0, ... cplex.mip.tolerances.mipgap, 0.001); optimize(Constraints, Objective, options); % 输出结果 C_opt.C_PV value(C_PV); C_opt.C_WT value(C_WT); C_opt.C_bess_E value(C_bess_E); C_opt.C_bess_P value(C_bess_P); C_opt.C_MT value(C_MT); LB value(Objective); details.solver_output ... end这里有个YALMIP使用上的细节沿着场景循环添加变量和约束时每次循环都要新建sdpvar变量不能把不同场景的第二阶段变量共用同一个变量名否则YALMIP会将其视为同一个变量导致场景之间的约束耦合在一起模型就完全错了。我在第一次写这个主问题时就在这里栽过跟头当时所有场景共用一组P_MT、P_ch变量结果求解出的容量配置明显偏大检查半天才发现是变量名作用域出了问题。5.3 子问题对偶建模与最坏场景求解子问题的目标是在给定容量配置 (C^*) 的条件下找到使运行成本最大的风光出力场景。%% 子问题求解函数寻找最恶劣场景 function [Q_val, worst_scenario] solve_subproblem(C_current, params) % 解包当前容量配置 C_PV C_current.C_PV; C_WT C_current.C_WT; C_bess_E C_current.C_bess_E; C_bess_P C_current.C_bess_P; C_MT C_current.C_MT; T params.T; % 不确定变量各时段实际风光出力 P_pv_actual sdpvar(T, 1); P_wt_actual sdpvar(T, 1); % 第二阶段运行变量 P_MT sdpvar(T, 1); P_ch sdpvar(T, 1); P_dis sdpvar(T, 1); SOC sdpvar(T1, 1); P_buy sdpvar(T, 1); P_sell sdpvar(T, 1); P_curt sdpvar(T, 1); % 不确定集约束箱式约束加预算约束 z_pv binvar(T, 1); % 光伏偏差指示变量 z_wt binvar(T, 1); % 风电偏差指示变量 Constraints []; for t 1:T % 光伏出力不确定范围 Constraints [Constraints, ... P_pv_actual(t) params.PV_forecast(t) * C_PV * (1 - alpha_PV 2 * alpha_PV * z_pv(t))]; % 风电出力不确定范围 Constraints [Constraints, ... P_wt_actual(t) params.WT_forecast(t) * C_WT * (1 - alpha_WT 2 * alpha_WT * z_wt(t))]; end % 预算约束 Constraints [Constraints, sum(z_pv) sum(z_wt) Gamma]; % 第二阶段运行约束与主问题类似 ... % 子问题目标最大化运行成本 RunCost sum(c_gas * P_MT / eta_MT c_buy * P_buy - ... c_sell * P_sell c_curt * P_curt); Objective -RunCost; % 注意max 转 min % 求解 options sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, Objective, options); Q_val value(RunCost); worst_scenario.P_pv value(P_pv_actual); worst_scenario.P_wt value(P_wt_actual); end这里采用的是“直接构造原始子问题”的写法利用YALMIP可以对整数变量直接建模的能力来处理不确定集中的预算约束。这种做法的好处是直观、不容易出错代价是求解速度比纯线性对偶问题要慢一些。如果模型规模特别大建议还是走对偶加线性化的路。5.4 CCG循环控制主程序和子程序都准备好后把它们串起来的迭代框架如下%% CCG主循环 function result CCG_main(params) % 初始化 UB inf; LB -inf; iter 0; max_iter 20; epsilon 0.01; % 初始场景名义预测场景 worst_scenarios struct([]); worst_scenarios(1).P_pv params.PV_forecast .* params.C_PV_init; worst_scenarios(1).P_wt params.WT_forecast .* params.C_WT_init; while iter max_iter iter iter 1; fprintf(迭代次数: %d\n, iter); % 1. 求解主问题得到容量配置和新的下界 [C_current, LB_new, mp_detail] solve_master_problem(worst_scenarios, params); LB LB_new; % 2. 固定当前容量配置求解子问题得到最恶劣场景 [Q_val, worst_scenario_new] solve_subproblem(C_current, params); % 3. 计算当前上界投资成本 最恶劣场景运行成本 InvestCost inv_cost_PV * C_current.C_PV inv_cost_WT * C_current.C_WT ... inv_cost_BESS_E * C_current.C_bess_E inv_cost_BESS_P * C_current.C_bess_P ... inv_cost_MT * C_current.C_MT; UB min(UB, InvestCost Q_val); fprintf(LB %.4f, UB %.4f, gap %.4f%%\n, LB, UB, abs((UB - LB) / UB * 100)); % 4. 判断收敛 if abs((UB - LB) / UB) epsilon fprintf(算法收敛于第 %d 次迭代\n, iter); result.C C_current; result.UB UB; result.LB LB; result.iter iter; result.worst_scenario worst_scenario_new; return; end % 5. 将新场景加入主问题场景集合 if isempty(worst_scenarios) worst_scenarios(1) worst_scenario_new; else worst_scenarios(end1) worst_scenario_new; end end error(达到最大迭代次数算法未收敛); end这里要特别强调一个细节每次迭代后保存的新场景 (u_k^*) 是子问题在当前容量配置 (C_k) 下求出的“最恶劣场景”。把这个场景加入主问题本质上是告诉主问题“你算出来的这套容量配置在遇到这个场景时可能要付出很高的运行成本你重新优化一下。”这个过程反复进行容量配置方案会逐步趋向于“在所有可能的最恶劣场景下都可接受”的鲁棒方案。6. 仿真算例与结果分析配置方案与算法性能6.1 算例参数为了验证模型和算法的有效性我设计了一个典型日微网算例。微网包含光伏、风电、储能、燃气轮机并与配电网保持连接。主要参数如下表参数数值单位微网峰值负荷860kW光伏最大安装容量1200kW风电最大安装容量600kW储能最大额定容量1000kWh储能最大额定功率500kW燃气轮机最大安装容量800kW光伏单位投资成本3500元/kW风电单位投资成本6000元/kW储能单位容量成本1500元/kWh燃气轮机单位投资成本3000元/kW折现率0.06-设备寿命20年6.2 迭代收敛与最坏场景运行CCG主循环后典型的收敛过程如下迭代次数下界LB万元/年上界UB万元/年间隙%1452.3489.67.62468.5485.23.43474.8481.31.44476.2479.50.75477.0478.40.3算法在5次迭代内收敛到0.3%的间隙效果相当理想。这说明CCG处理这个规模的问题收敛速度很快不需要担心计算负担。最恶劣场景揭示了一个重要现象出现最恶劣组合的时段集中在傍晚17:00到21:00——此时光伏出力变为0风电出力同时处于低谷负荷却达到全天最高。这个时段系统几乎完全依赖储能放电和燃气轮机出力加上从主网购电才能维持供电。这个结论很有价值它告诉我们在容量配置时要把“晚高峰时段”作为最关键的校核场景。6.3 配置结果对比鲁棒解vs确定性解为了展示两阶段鲁棒优化的效果我把鲁棒模型的结果与确定性模型进行对比。配置方案光伏容量kW风电容量kW储能容量kWh储能功率kW燃气轮机kW年总成本万元确定性解900420650280300458.6鲁棒解Γ8820380780320420478.3可以看到鲁棒解相比确定性解光伏和风电的安装容量有所降低储能和燃气轮机的容量有所增加年总成本上升了约4.3%。这个结果从工程角度非常好理解光伏和风电是不确定电源过多依赖它们会在恶劣场景下造成供电缺口因此鲁棒解主动降低了这两个的配置比例储能和燃气轮机是可控电源它们才是应对不确定性的“压舱石”因此鲁棒解增加了这两个的容量。多花的4.3%投资换来的是系统在极端天气条件下的可靠供电能力这笔账在经济性分析中需要结合停电损失一起评估。我在论文中通常会加一组对比如果把切负荷惩罚成本拉高模拟重要负荷场景鲁棒解的经济优势会立刻体现出来。6.4 保守度参数对结果的影响预算参数 (\Gamma) 直接控制鲁棒模型的保守程度。我扫描了不同的 (\Gamma) 值观察它对成本的影响。Γ值年总成本万元相比Γ0增幅%最恶劣场景运行成本万元/天0确定性458.6-0.854466.21.70.928478.34.31.0412489.56.71.1516497.88.51.2424503.29.71.30可以看到随着 (\Gamma) 增大总成本单调上升但增幅逐渐放缓。这说明鲁棒性带来的成本增加呈现边际递减效应当 (\Gamma) 达到一定值后继续增加对系统鲁棒性的提升已不明显但对经济性的损害却在累积。实际工程中(\Gamma) 不应一味取大。我建议的做法是基于历史气象数据统计一个观察周期内风光出力的最大偏差时长占比以这个统计值作为 (\Gamma) 的选取参考这样既避免了拍脑袋定参数又保证了模型的鲁棒性贴合实际需求。7. 实操中容易踩的坑与调试建议7.1 YALMIP求解器配置问题YALMIP 在使用中最大的坑是求解器调用不明确。当同时安装了多个求解器时比如既装了CPLEX又装了GurobiYALMIP有可能会自动选择别的求解器导致结果异常或者性能骤降。我在仿真中遇到过一个问题明明安装好了CPLEX求解时却提示“No suitable solver”。排查后发现是因为YALMIP版本过旧无法识别新版本CPLEX的接口。解决办法是升级YALMIP到最新版并确保CPLEX的MATLAB接口路径已正确添加。建议在每次调用前显式指定求解器options sdpsettings(solver, cplex, verbose, 1);不要依赖YALMIP的默认求解器选择机制。7.2 CPLEX整数规划求解缓慢在子问题中预算约束使用了0-1整数变量 (z_{PV,t}) 和 (z_{WT,t})这意味着子问题是一个混合整数规划。当 (T) 增大到96即15分钟一个间隔或者不确定变量的时段数变多时求解时间会指数级增长。我最初的仿真中T24时子问题求解只需1到2秒但当我把时间分辨率提高到15分钟T96后单次子问题求解时间飙升到接近60秒整个CCG流程跑完需要近10分钟。解决思路有两个方向一是对偶线性化。将子问题转化为纯线性规划去掉整数变量这样求解速度会快很多。具体做法是引入不确定变量的连续辅助变量用KKT条件或强对偶理论处理内层min问题。二是场景缩减。对历史数据进行聚类分析将相似的时段合并减少需要单独建模的时段数。我实际采用的是第一种方法通过对偶转化将子问题中的整数变量替换为连续变量求解时间降到了5秒以内。7.3 对偶问题的正确性验证手动推导对偶问题时最容易犯的错误是对偶方向写反、约束符号搞错、对偶变量漏项。这些问题不会导致模型无解但会给出一个错误的最恶劣场景进而污染整个迭代过程。验证对偶问题是否正确有一个非常实用的方法取一组相同的参数分别求解原始子问题和对偶问题对比两者的目标函数值。由于强对偶定理保证两者相等如果数值不一致说明对偶推导或实现有问题。我在代码里写了这样一个自检函数function verify_dual(C_current, params) % 求解原始子问题 Q_primal solve_subproblem_primal(C_current, params); % 求解对偶子问题 Q_dual solve_subproblem_dual(C_current, params); % 对比 fprintf(原始问题目标值: %.4f\n, Q_primal); fprintf(对偶问题目标值: %.4f\n, Q_dual); if abs(Q_primal - Q_dual) 1e-6 error(对偶验证失败); end end每次都跑一遍自检发现问题即刻排查能节省大量后期调试时间。7.4 大M法与数值稳定性在使用大M法处理双线性项时M值的选取直接影响求解稳定性。M值过小会导致不可行解M值过大又可能引起数值问题导致CPLEX出现严重病态矩阵警告。我在调试时发现M值并非越大越好。CPLEX官方文档里建议M的取值应该比模型中的其他系数大一到两个数量级即可不要盲目取1e9这类超大值。一般我会根据模型参数估算M的范围功率平衡约束中的M值取最大可能的功率缺额如 (M \ge \max(P_{load}) \max(P_{ch}) \max(P_{sell}))容量上限约束中的M值取对应容量的上限值即可。此外在YALMIP中使用大M法时尽量把M写成具体的数值而不是一个很大的1e6常量这样CPLEX在做预求解时能够更有效地收紧边界提升求解速度。8. 从仿真到论文结果呈现与延伸方向8.1 仿真结果的可视化呈现仿真做完后出图是必不可少的一步。我在论文中一般会准备四类图第一类是迭代收敛曲线图展示CCG算法中上下界随迭代次数的变化过程证明算法的收敛性。可以生成这样的图figure; plot(1:iter, LB_history, -o, LineWidth, 1.5); hold on; plot(1:iter, UB_history, -s, LineWidth, 1.5); xlabel(迭代次数); ylabel(成本万元); legend(下界LB, 上界UB, Location, northeast); grid on; saveas(gcf, convergence_curve.png);第二类是容量配置对比图用堆叠柱状图展示确定性模型和鲁棒模型在不同电源上的容量配置差异。第三类是最恶劣场景下的功率平衡图展示在最恶劣场景下各电源的出力分配情况和系统整体功率平衡关系。这张图最能直观展示鲁棒方案在极端工况下的应对能力。第四类是用折线图展示不同 (\Gamma) 值下年总成本和切负荷量的变化趋势直观呈现鲁棒性与经济性的折中关系。8.2 两阶段鲁棒优化模型的扩展方向两阶段鲁棒优化在微网电源容量配置中的应用远不止我上面实现的这些内容。如果你做完基础模型想往更深的方向走有几个我认为很有价值的延伸方向第一多微网互联场景。多个微网之间通过联络线互济可以大幅提高系统整体的鲁棒性降低成本。但这种模型会引入更多的一阶段决策变量计算复杂度急剧上升需要配合分解算法或分布式求解框架。第二考虑需求响应的鲁棒配置。把可调负荷纳入优化框架让负荷侧也参与削峰填谷。需求响应本质上是给系统增加了柔性资源在不确定环境下价值更为明显。第三多阶段鲁棒优化。两阶段模型适合“先配置、后运行”的场景但实际运行中调度决策往往是滚动时域的。多阶段鲁棒模型能更精确地刻画这种时序决策关系但建模难度和求解复杂度都会显著上升。第四将鲁棒优化与机器学习结合。用LSTM或Transformer对风光出力进行概率预测将预测区间作为不确定集的边界这样不确定集不再是人为设定而是由数据驱动生成模型的实用性会更强。8.3 一些个人心得做这个课题前前后后大概花了两个月。最大的感受是两阶段鲁棒优化本身并不难理解难的是把模型转成代码、并让代码在各种边界条件下都能稳定运行。数学推导和代码实现之间存在一条需要反复跳转的鸿沟跳过去了前面的路就顺了。如果你刚开始接触这个方向我的建议是先不管YALMIP手推一个最简单的两阶段鲁棒模型比如单时段、单设备然后手动用对偶转换加CCG求解一遍完全吃透整个流程。再回到MATLAB里用YALMIP写代码你会觉得一切都是顺理成章的。工具是辅助对模型本质的理解才是决定你能否做深做透的关键。