两阶段鲁棒优化在微电网多电源容量配置中的应用与Matlab实现

微电网容量配置这件事我是真吃过亏。早年做园区微网的源网荷储规划拿典型日的光伏、风电、负荷曲线做确定性优化算出配置方案时还挺满意。结果项目落地第三个月连续一周阴雨又叠加一次强阵风光伏出力腰斩风电出力剧烈波动储能很快就触底柴发长时间半载运行那一个月的整体运行成本比设计值高了近三成。那段经历让我彻底明白了一个道理微网多电源容量配置的真正难点不在容量怎么算而在于怎么让方案在最不利的情况下仍然扛得住。而这个问题的标准解法就是两阶段鲁棒优化算法。这篇文章我就把两阶段鲁棒优化 微网多电源容量配置 Matlab代码实现这条线完整讲清楚从原理推导到代码骨架再到调参避坑一次讲透。1. 从拍脑袋选容量到两阶段鲁棒优化为什么这个问题值得做1.1 微网多电源容量配置到底在解决什么问题微网通常是指包含风电、光伏、柴油机/燃气轮机、储能以及各类负荷的小型电力系统可以并网运行也可以离网运行。所谓多电源容量配置就是在满足供电可靠性、储能循环寿命、机组运行边界等一堆约束的前提下决定风电机组装几台、光伏组件装多少、储能电池配多大、柴油机选什么规格使得整个微网在生命周期内的总成本最低。这里有个很容易被忽视的点容量决策是一次性的、前置的而运行调度是反复发生的。你把风机装多了在风速平缓的季节风机就是一堆闲置资产你把储能配大了前期投资沉没成本高还可能常年用不满。所以这个问题本质上就带着先投资、后运行的时间先后特征。1.2 传统确定性规划的局限很多人第一步想到的是确定性规划把风速、光照、负荷统一替换成典型日或最恶劣单一场景。我在早期也是这么做的结果就是上面那个翻车案例。典型日方法的问题在于它用平均值或加权平均掩盖了真实波动的相关性。连续阴雨天、夜间大负荷叠加低风速、台风前后的极端风况这些组合根本不是一条典型日曲线能表达的。也有人用最恶劣单一场景来定容量但场景怎么选、多恶劣才算够人工挑场景本质上是把不确定性拍成几个离散点一旦真实场景超出预定义范围配置方案就会失效。而且单一最恶劣场景往往会同时假设所有时段都极端这又太过于保守投资成本高到业主无法接受。1.3 鲁棒优化用不确定集合代替单一场景鲁棒优化做的事情是把所有可能的波动描述成一个不确定集合 U然后寻找一个容量配置方案使得即便不确定性参数落在集合里最坏的位置系统依然能够通过合理的运行调度满足所有约束。换句话说它优化的是最坏情况下的总成本而不是期望情况下的总成本。这样得到的配置结果看起来会比确定性方案略贵但换来的是很强的不后悔特性。两阶段鲁棒优化和单阶段鲁棒优化的区别在于单阶段要求容量和调度必须同时、一次性决定这往往太保守因为现实中调度是可以根据实际风光的观测结果来调整的而两阶段则允许先决策容量、等不确定性参数实现后再做调度决策。这个先决策、后调整的时序恰好和微网容量配置的物理过程完全一致。2. 容量配置问题的数学建模骨架2.1 系统假设与需要决策的东西建模之前先把系统边界说清楚。以最常见的并网/离网混合微网为例假设系统包含风电机组WT、光伏PV、柴油发电机DE、蓄电池储能ESS以及电负荷。第一阶段变量是装机决策包括各类电源的安装台数/容量。为了贴近实际工程装机台数通常建模为整数变量比如 x_w ∈ {0,1,2,...} 表示风电装机台数x_pv 表示光伏装机容量等级x_ess 表示储能容量等级x_de 表示柴油机台数。这样问题就变成一个混合整数规划MIP问题。第二阶段变量是运行调度变量包括每个时段 t1...T 的风电实际出力 P_w,t、光伏出力 P_pv,t、柴油机出力 P_de,t、储能充放电功率 P_ch,t 和 P_dis,t、蓄电池 SOC_t以及可能发生的切负荷量 P_cur,t。第二阶段变量在容量确定后、不确定场景实现时才能决定这正是两阶段结构的来源。2.2 目标函数年综合成本怎么算容量配置类问题的目标函数一般写成年综合成本最小化。公式不粘贴复杂推导直接给出工程里常用的结构min C_inv 365 × (C_fuel C_om C_cur)其中 C_inv 是年化投资成本把风机、光伏、储能、柴油机各自的一次性投资按利率和寿命折成每年费用C_fuel 是柴油机燃料成本C_om 是运维成本C_cur 是切负荷惩罚成本。为什么要把燃料成本和切负荷惩罚都放进来因为鲁棒优化要在所有最坏场景里找成本最高的那个而切负荷惩罚是约束系统可靠性的代表。没有惩罚项优化器可能会选择最小化投资成本、然后疯狂切负荷这种技术上可行、实际上荒诞的方案。常见做法是切负荷惩罚设成远高于发电成本的数量级比如 10 元/kWh 以上。2.3 约束条件功率平衡、储能SOC、柴油机运行边界约束部分至少要有下面几类。功率平衡约束是所有微网模型的地基每个时段风电、光伏、柴油机、储能放电功率、购电并网时之和等于负荷加储能充电功率加切负荷。一个关键细节是风电和光伏在确定性模型里是固定参数在鲁棒模型里要被替换成不确定参数。P_w,t^u P̄_w,t u_w,t × Δ_w,t其中 P̄_w,t 是预测出力Δ_w,t 是最大偏差u_w,t 限制在 [-1,1]表示实际出力可以在预测值附近波动。储能约束包括 SOC 递推关系、容量上限、SOC 上下限以及充放电功率上限。SOC 递推的经典写法是SOC_t1 SOC_t η_ch × P_ch,t − P_dis,t / η_dis柴油机约束包括出力上下限、爬坡约束、最小启停时间。此外还有容量与运行变量的耦合约束比如储能充放电功率不能超过已装容量对应的上限风电实际出力小于等于装机容量对应的预测上限。2.4 为什么天然是两阶段结构我经常被问为什么这里一定要用两阶段不能把容量和调度写在一起同时优化吗可以但那样做就没有兑现调度决策可以根据实际天气来调整这个信息优势。真实微网的调度是这样的容量是年初就定死的而风光的实际出力要等到当天、当时才知道然后运行人员才据此安排柴油机和储能。两阶段模型把知道不确定值之前的决策和知道不确定值之后的决策分开处理正好反映了这个信息结构。如果合成单阶段优化等于强迫决策者在不看天气的情况下把每个时段的运行功率都定死结果只会更保守、成本更高而且违背物理直觉。3. 两阶段鲁棒优化原理min-max-min到底在解什么3.1 不确定变量的建模方式先看最常用的盒式不确定集合。对风电而言实际出力 P_w,t 可以在预测值附近波动U { u ∈ R^(N_T) : ||u||_∞ ≤ 1, ||u||_1 ≤ Γ }这个公式看起来有点抽象翻译成人话就是每个时段的风电出力相对预测值最多偏离一定百分比无穷范数限制同时所有时段累计偏离总量不超过一个预算 Γ1-范数限制。为什么要加 Γ因为如果允许每个时段都同时处于极端偏差那造出来的方案会过于保守、成本高得离谱。预算 Γ 让模型可以描述虽然极个别时段会很恶劣但不是所有时段同时恶劣的现实。这个思想是鲁棒优化的灵魂也直接影响配置结果——后面第 5 节我会给一个不同 Γ 下配置变化的对照。3.2 两层决策的嵌套结构直观解释两阶段鲁棒优化的标准形式是min_x { c^T x max_(u∈U) min_(y∈Ω(x,u)) d^T y }其中 x 是容量决策y 是调度决策u 是风光不确定参数。只看公式容易晕我习惯用买保险来类比。你想买一台备用发电机来应对小区停电。x 是买多大的发电机属于现在就要拍板的事u 是未来停电最严重的程度属于你无法控制的最坏情况y 是停电发生时你怎么调度这台发电机、要不要给电梯优先供电属于事到临头才能做的安排。两阶段鲁棒就是问在最坏的停电情况下我买的发电机和事后的调度策略配合起来总成本还能不能接受如果太贵就换一个容量方案重新评估。这样理解min-max-min 就不是玄学了而是一个很自然的先定方案、再被最坏现实考验、再事中应对的过程。3.3 CCG算法主问题与子问题的拆分两阶段鲁棒问题直接求解非常困难因为 max 和 min 嵌套不是一个标准可解模型。工程和学术界最常用的解法是列与约束生成CCG算法也叫场景割法。思路是把原问题拆成一个主问题MP和一个子问题SP。主问题在一组已知的恶劣场景 u_1, u_2, ..., u_k 下做容量调度联合优化得到下界解子问题固定主问题给出的容量 x去寻找能够让运行成本最大化的新恶劣场景 u_(k1)把该场景下的运行成本作为上界解。如果上下界差距小于阈值迭代停止否则把这个新场景以及对应的调度变量一起加入主问题再重复。注意这个关键点主问题每次不是只加一条割约束而是把新场景对应的整套第二阶段变量和约束都加进去所以叫列与约束生成而不是经典的 Benders 割。这样处理的好处是收敛速度快很多。在容量配置这类问题上通常迭代十几次以内就能收敛。3.4 子问题求解的关键对偶转换与双线性项处理子问题固定 x 之后是关于 y 和 u 的问题内层 min 是 LP外层是 max。直接 max-min 无法用通用求解器处理通常做法是先对 min 用拉格朗日对偶转成 max再把两个 max 合并于是目标变成关于对偶变量 λ 和不确定变量 u 的函数。麻烦的是这里大概率出现 λ 乘 u 的双线性项——这两个都是变量乘积让问题变成非凸。我见过很多 demo 代码在这里含糊带过实际落地时最常用的招数是引入辅助变量和大 M 法把它线性化。比如 λ_i × u_j引入 z_ij 替代并添加形如 z ≤ M×u_indicator、z ≤ λ 的约束把原问题转换成 MILP 后用 Cplex/Gurobi 直接求解。也可以走 KKT 路线把内层 min 的最优性条件直接写进外层 max变成带互补条件的 MPEC 再用求解器算。这个思路在主问题规模不大的时候没问题但互补条件含大 M 参数设置同样很考验经验。所以我个人更推荐对偶线性化这条路线。4. Matlab代码实现从公式到可运行程序4.1 环境准备Yalmip Cplex/GurobiMatlab 实现我建议用 Yalmip Cplex/Gurobi 这套组合。Yalmip 能让你直接写 sdpvar、intvar 和约束避免手写求解器 API 的体力活而 Cplex 的 MIP 求解性能在容量配置这种问题上明显比 Matlab 内置的 intlinprog 稳。如果你的 Matlab 版本比较新Yalmip 安装只需要把解压后的 yalmip 文件夹加入 matlab path 就行。Cplex 或 Gurobi 记得要安装对应的 Matlab 接口安装完成后在命令行输入yalmiptest能跑通就说明环境没问题。这个环节经常有人卡在 solver 找不到多半是路径没配对或者 license 没激活建议先在自带 yalmip demo 上验证再往下走。4.2 变量定义与参数初始化用 Yalmip 写这个问题的第一步是定义变量。我习惯把第一阶段变量和第二阶段变量分开方便 CCG 迭代时做场景扩展。T 24; % 调度周期小时 nW 10; % 风电候选最大台数 nPV 20; % 光伏候选组数 nESS 5; % 储能候选容量等级 nDE 3; % 柴油机候选台数 % 第一阶段变量容量决策整数 x_w intvar(1, nW, full); x_pv intvar(1, nPV, full); x_ess intvar(1, nESS, full); x_de intvar(1, nDE, full); % 第二阶段变量调度决策连续 P_w sdpvar(1, T, full); P_pv sdpvar(1, T, full); P_de sdpvar(1, T, full); P_ch sdpvar(1, T, full); P_dis sdpvar(1, T, full); SOC sdpvar(1, T, full); P_cur sdpvar(1, T, full);注意装机台数用intvar而不是binvar因为允许装多台同型号机组不是简单的装或不装。如果候选容量是连续变量也可以换成sdpvar但工程上设备都有标准规格整数更真实。4.3 主问题MP的编写CCG 主问题需要给每个已发现的恶劣场景 k 单独建立一套第二阶段变量。我写成 cell 数组每轮迭代新增一个场景就新增一组变量% 假设已经找到场景集 u_set尺寸为 T×K每列是一个场景 P_w_k cell(1, K); P_pv_k cell(1, K); P_de_k cell(1, K); % ... 其他调度变量同理 for k 1:K P_w_k{k} sdpvar(1, T, full); % 对应约束当前场景下的风电出力上限 Constraints [Constraints, P_w_k{k} u_set(:,k) .* P_w_max_pred]; % 功率平衡P_w_k P_pv_k P_de_k P_dis_k P_buy_k % P_load P_ch_k P_cur_k Constraints [Constraints, P_w_k{k} P_pv_k{k} P_de_k{k} ... P_dis_k{k} P_load P_ch_k{k} P_cur_k{k}]; % ... 储能SOC、柴油机上下限等 end Objective inv_cost(x_w, x_pv, x_ess, x_de) ... 365 * sum(op_cost_k); % 各场景运行成本加权主问题相对好写但它有很强的多场景扩展属性。每次迭代新增的不是一行约束而是一整块即新场景对应的全部调度变量与约束。这是 CCG 和 Benders 在编码上的最大差异写的时候别图省事只加约束不加变量否则新约束引用不存在的变量会直接报错。4.4 子问题SP的编写子问题固定 x 后求解恶劣场景。这里按第 3.4 节的对偶大M写法给出核心骨架注意目标方向和约束符号要按你自己的模型推导来调整。%% 子问题SP固定x找最恶劣u lambda sdpvar(1, size(A_dual,1), full); % 对偶变量 u sdpvar(1, T, full); % 不确定变量 % 双线性项 lambda_i * u_j 用大M线性化 M 1e3; % 根据实际变量取值范围设置 z binvar(size(lambda,2), T, full); w sdpvar(size(lambda,2), T, full); for i 1:size(lambda,2) for j 1:T % w(i,j) 替代 lambda(i)*u(j) Constraints [Constraints, 0 w(i,j) M*z(i,j)]; Constraints [Constraints, lambda(i)-M*(1-z(i,j)) w(i,j)]; Constraints [Constraints, w(i,j) lambda(i)M*(1-z(i,j))]; Constraints [Constraints, -M*(1-z(i,j)) w(i,j)-u(j); w(i,j)-u(j) M*(1-z(i,j))]; end end Objective_SP -(lambda*b0 sum(sum(w .* c_dual)) ...); % 注意最小化转对偶后目标方向要调整符号不要写反这里要特别提醒子问题求出来的目标值在最小化问题中对应上界 UB主问题的目标值对应下界 LB。迭代条件是 gap (UB-LB)/UB 小于某个阈值我通常取 0.5% 或 1%。很多初学者把上下界方向搞反导致明明收敛了却一直迭代这个必须检查清楚。4.5 CCG迭代主循环主循环不复杂就是把 4.3 和 4.4 反复串起来tol 0.005; % 收敛阈值 maxIter 30; UB inf; LB -inf; u_set ones(1, T); % 初始场景可以用预测场景 for iter 1:maxIter % 1. 把u_set加入主问题求解当前MP optimize(Con_MP, Obj_MP, options); x_val [value(x_w), value(x_pv), value(x_ess), value(x_de)]; LB value(Obj_MP); % 2. 固定x_val求解子问题SP [obj_sp, u_new] solve_SP(x_val); UB min(UB, obj_sp); % 上界取历史最小可行值 if (UB-LB)/UB tol break; end % 3. 把u_new作为新场景加入主问题 u_set [u_set, u_new]; end这里有个我踩过好几遍的坑上界是历史可行解中的最优点所以要取 min下界是当前近似问题的最优下界理论上随迭代单调不减。如果你看到 LB 比 UB 还大先别急着怀疑算法检查是不是主问题忘记加某条关键约束或者子问题对偶方向写反了。CCG 的收敛曲线应该是上界快速下降、下界缓慢上升最终贴合在一起。5. 跑通之后结果怎么读、参数怎么调5.1 收敛判据与迭代曲线我自己的经验是容量配置这类问题用 CCG 通常 5~15 轮就能收敛到 1% 以内。如果超过 30 轮还不收敛问题多半不在算法而在模型本身。常见的病根有三个子问题对偶少写了约束导致上界虚高不确定性集合参数 Γ 设得过大导致最坏场景重复跳跃大 M 系数设置不当导致子问题数值不稳定、场景来回反复。调参之前先画出 UB 和 LB 随迭代次数的曲线一眼就能分辨是哪种病UB 迟迟不下来通常是子问题建模问题LB 和 UB 各自在跳动一般是数值问题UB 快速下降但 LB 纹丝不动往往是主问题约束不足。5.2 不确定性水平与容量配置结果的对应关系我跑一个 24 时段的小算例时把 Γ 从 0 调到 6、12、24得到一组很有代表性的结果趋势以下是演示数据不是真实工程数值不确定预算 Γ风机(台)光伏(组)储能(等级)柴油机(台)年综合成本(万元)0确定性31221约 385631121约 412122932约 45824最保守2842约 527Γ0 时光伏装得多、储能配得少因为确定性模型觉得白天光照稳稳的Γ24每个时段都允许极端偏差时储能容量几乎翻倍、柴油机也加了一台同时光伏反而减少——因为光照不可靠时光伏的投资收益比下降了。这个趋势提醒我们鲁棒性不是白送的它表现为投资结构的变化而不仅仅是总成本上升。做项目汇报的时候我最怕别人只给一套数字所以现在习惯性地把不确定预算-总成本-各类电源容量放在同一张表里让决策者能看到每多一分鲁棒性要花多少钱换。5.3 求解时间与规模控制容量配置问题的规模很容易失控T 取 24 还好一旦做全年 8760 小时第二阶段变量数直接翻 365 倍求解时间立刻变得不可接受。工程上常用典型日聚类把全年压缩成 4~8 个典型日并分配权重学术上也有用两阶段情景缩减的做法。还有几个能明显提速的经验主问题里的整数变量很多时可以先用松弛解给 Cplex 提供热启动开启 Cplex/Gurobi 的 MIP emphasis 和流覆盖预处理选项别在最开始就追求 8760 全时段全场景先用 24 时段把模型跑通比什么都重要跑通后再延展到全年。5.4 鲁棒性代价怎么评估最后说说怎么向别人解释为什么贵了。把鲁棒方案的总成本减去确定性方案的总成本这部分差额叫鲁棒性代价本质是为应对不确定性支付的保险费。评估值不值得可以拿历史实测数据回代把两套方案分别放进真实的历史场景里跑模拟统计年成本、切负荷量、储能过放次数。如果确定性方案在极端年份切负荷严重鲁棒方案就值得如果历史数据里根本没有极端情况那确实可以不花这个钱。这个回代验证步骤看起来多花三五天但在项目验收时价值很大它能直接回答你这套鲁棒方案到底比拍脑袋方案强在哪。6. 实操中容易踩的坑与个人经验6.1 可行性割与最优性割必须分清楚CCG 迭代过程中子问题可能出现对某个容量配置 x 根本无可行调度的情况比如储能容量配得太小在最坏场景下无法同时满足负荷和 SOC 下限。这时候子问题返回的不应该是最优性割而是一条可行性割需要额外求解一个带松弛变量的辅助问题找到让约束条件被破坏到最小的场景把这条场景也加入主问题。很多公开代码不处理可行性问题遇到不可行直接报错或返回 inf结果迭代流程照样跑出来的配置却是错的。我的做法是在子问题里预设几个松弛变量切负荷、SOC 越限求解后先看它们是否为 0。如果最优解里任意松弛变量大于阈值就加可行性割而不是最优性割。6.2 整数变量和不确定变量的耦合处理如果第一阶段容量变量是整数台数第二阶段出力上限会写成x_w × P_W_max这类耦合形式。这本身没问题但有些同学会在子问题里把 x 当成连续变量处理偷偷放松整数性结果回代出来的最坏场景在真实整数容量下并不存在。我吃过这个亏后来规定子问题里 x 完全固定为当前整数解取值用value(x)直接代入而不是写成约束再让求解器去调。另一个相关坑是可行性割返回的新场景可能让主问题又生成一个完全不同的整数解导致场景集合来回震荡。这时可以在割里加上 x 的邻域限制或者在目标里加一个小的正则项让迭代曲线平滑下来。6.3 大M参数的取值是数值稳定性的命门对偶大M线性化是目前工程落地的主流方案但它对大M极其敏感。M 太小会砍掉可行域M 太大会让单纯形表的数值条件恶化在 Cplex 日志里通常表现为Numerical difficulties警告或者场景割来回震荡。我建议分三步处理对每个双线性项根据 λ 和 u 各自的取值范围分别算出局部大M——不要全局共用一个 1e6给 Cplex 开启数值强调参数numerical emphasis用小规模随机数据验证子问题目标值和用穷举台风场景算出的目标值做对比确认一致后再上完整模型。第三步很土但确实能救命。我之前有一次迭代 100 轮不收敛最后发现就是某个双线性项的 M 设小了砍掉了真正最坏的场景。6.4 我现在的固定调试顺序被这些坑教育过之后我形成了一套固定的开工顺序现在做类似项目都是这么推的把 Γ 设为 0整个模型退化成确定性规划先确保主问题单独能跑出合理结果手动固定一个容量方案单独调试子问题观察给定最坏场景时调度成本是否合理以及找出的场景是不是你直觉中的恶劣组合再跑完整 CCG重点看 UB/LB 曲线是否平滑收敛换不同 Γ 和不同预测曲线做参数敏感性测试。这套流程看起来慢实际比直接写完整代码然后反复 debug 快得多。如果你现在正在做类似的项目我建议至少留出两三天给前两步后面会省下好几周。两阶段鲁棒优化的容量配置方案不是跑完代码就结束的它需要你反复回到最坏场景是什么样子这个直觉问题上模型才真正可信。