做微网运行调度的人应该都有过这种纠结手里的分布式光伏、储能、柔性负荷越来越多但用传统集中式优化算出来的调度结果拿到实际现场根本执行不下去。原因很简单——微网里的各个主体都有自己的利益诉求一个优化模型里硬塞进去的目标函数根本代表不了那些真实决策者。后来我开始接触合作型Stackelberg博弈配合Matlab/Cplex做微网运行策略优化算是把这层利益关系在数学模型里理顺了。这篇内容就是把整个思路和工程实现细节拆开来讲适合正在做微网双层优化、主从博弈或者打算用Cplex求解混合整数规划的同学参考。1. 为什么微网运行问题会演变成主从博弈1.1 微网运行的本质不是单目标优化微网可以看成一个小型电力系统有分布式电源光伏、风电、燃气轮机、储能、各类负荷还能通过公共连接点和配电网交换功率。传统做法是把微网看成一个整体以运行成本最低或新能源消纳最大为目标用最优潮流或混合整数线性规划求解。这个思路在微网规模小、产权单一的时候没问题比如园区微网整体由一家公司建设运营大家的目标天然一致。但当分布式光伏走进千家万户电动汽车和柔性负荷参与市场互动后微网内部的利益主体开始分化。微网运营商管理公共母线、储能和联络线追求运营收益最大化用户或者产消者自己也有光伏、储能同时还要用电他们追求的是自身成本最小和用能舒适度。这时候你要是再强行用一个目标函数刻画所有主体的行为要么忽略用户响应导致调度结果失真要么忽视运营商收益导致策略不具备可持续性。这就是为什么微网运行问题会从“单目标优化”演变成“多主体博弈”的根本原因。1.2 Stackelberg博弈在微网里的层级关系Stackelberg博弈在微网运行中的典型结构是这样的上层领导者微网运营商。先制定内部电价购电/售电价格下发到下层用户。下层跟随者多个用户/产消者。在给定电价下做自己的优化决策比如调整光伏上网功率、储能充放电、柔性负荷的运行时段。下层把最优响应返回上层上层再根据这个响应更新电价策略反复迭代直到达到均衡。这个结构的内在逻辑是“领导者有先动优势”运营商先出牌、定电价用户根据电价做出理性响应。运营商虽然是先出牌但也不能乱定价因为你定的价格太高或太低用户会用“少买电”或者“多卖电”的实际行动来回应你最终双方会在某个价格—电量组合上达成均衡。这种一主一从、按顺序决策的建模方式和微网“运营商定价—用户响应用能”的真实互动几乎一一对应。1.3 为什么还要加“合作型”三个字这里很多刚接触的朋友会误解以为“合作型”指的是上下层联合起来优化成一个大模型。实际不是。合作型Stackelberg博弈里的“合作”体现在收益分配机制上各主体参与博弈之后系统会产生一个总收益增量这个增量相比大家各自不互动、各用各的基准场景多出来的部分需要在成员之间按某种公平规则分掉。我最初做这个模型时也想过既然Stackelberg均衡本身就能提高效率为什么还要多此一举做收益分配后来在算例里对比才明白单纯的主从博弈解出来的均衡点确实能达到整体较优但收益分配如果不合理某些主体的收益可能还不如不合作。比如储能被运营商频繁调用充放电损耗全由储能承担但利润大头被运营商拿走储能所有者就没有动力继续参与合作。引入合作剩余和Shapley值分配之后可以保证每个参与主体分到的收益都不低于单独运行的基准值策略才真正具备可落地性。2. 从物理约束到博弈模型把运行策略“翻译”成数学规划2.1 决策变量怎么分层做博弈模型的第一步是把“谁决策什么”搞清楚变量设计直接决定后面求解的复杂度。我的做法是严格分两层上层变量运营商内部购电电价、售电电价联络线功率运营商侧储能的充放电功率燃气轮机出力。下层变量用户i光伏实际上网功率、可调负荷功率、用户侧储能充放电功率、与运营商的交互功率。上下层变量的耦合点就是电价。模型只通过价格信号传递信息既符合市场化交易机制也在数学上让双层模型有一个清晰的耦合结构。在实际代码里我会把这组耦合变量单独列出来方便后续构建互补约束和线性化表达式。这里有一个容易踩的坑运营商侧储能和用户侧储能要分开建模。很多同学偷懒把整个微网的所有储能放在一个状态变量里结果Shapley分配时根本算不清“这个储能的收益到底该归谁”因为储能产权主体不一样。我现在的习惯是每个储能节点单独建立SOC递推模型变量名里就能看出归属。2.2 上层目标函数和约束上层运营商的目标是收益最大化我常用的写法是收益 售电收入向用户卖电、向主网卖电 补贴收入 燃气轮机电费收入 - 从主网购电成本 - 储能运行维护成本 - 燃气轮机燃料成本 - 向用户支付的购电成本用户光伏卖给运营商 - 合作剩余分配支出其中燃气轮机的燃料成本一般写成二次函数或分段线性函数储能运维成本简化成单位功率成本乘以充放电功率绝对值。注意功率绝对值项在Cplex里需要处理要么引入辅助变量p_pos、p_neg要么改写成不等式约束不能让绝对值直接出现在目标函数里因为线性规划处理器不认这种非线性项。约束条件主要包含微网功率平衡约束电源出力 储能放电 主网购电 负荷 储能充电 主网售电联络线功率上下限燃气轮机出力上下限和爬坡约束储能SOC递推关系以及充放电状态互斥约束SOC递推是很基础但容易出错的部分。SOC(k1) SOC(k) - P_bat(k)*Δt/E_cap充电为正时这里就是减号取决于你定义的正方向一定要统一好充放电正方向我见过不止一次因为正方向定义不一致导致计算结果出现“一边充电一边放电”的矛盾场景。2.3 下层用户目标函数下层用户的优化目标我一般写成效用最大化形式max 用电效用 U_i(x) - 购电成本这里的用电效用函数可以用分段线性函数近似也可以用二次函数。二次函数的优点是能反映边际效用递减缺点是KKT条件会引入二次项线性化略微麻烦分段线性函数的好处是直接可用线性规划处理我倾向先用分段线性跑通以后想提升精度再换二次项。用户的约束包括光伏出力上限可上网功率不超过当前光伏预测出力可调负荷的可转移区间比如洗衣机、热水器可以在电价低谷时段运行用户侧储能SOC递推和容量约束与运营商的交互功率上下限用户侧变量和运营商侧变量之间通过电价相乘形成目标函数里的双线性项这是整个模型里必须重点处理的部分后面第3章会说。2.4 合作剩余与Shapley值的设计Shapley值的计算公式是这样的φ_i Σ_{S⊆N{i}} [ |S|! (n-|S|-1)! / n! ] * [ v(S∪{i}) - v(S) ]其中v(S)是联盟S的特征函数值也就是“联盟S里这些主体合作时系统的最大总收益”。n是参与主体总数。每个主体分到的收益φ_i等于他在所有可能联盟里的边际贡献的加权平均。这个分配方法具有公平性和可验证性是合作博弈里最常被引用的分配机制。在微网场景里计算Shapley值我对n3或n4的小规模算例直接枚举所有联盟每一个联盟对应一个优化问题重新求解一次得到v(S)然后把数值代入公式。这个过程别手工去做写一个自动化脚本自动生成联盟组合并调用求解器否则联盟数量稍微多一点就会在重复劳动里出错。原则上主体数量超过5个以后枚举法计算量大增那时候可以考虑改用班茨哈夫权力指数、核仁或者近似算法来分摊但对于学术验证场景Shapley值仍然是解释力最强的选择。3. MatlabCplex求解实现双层模型怎么落地3.1 环境配置里的实际顺序Matlab配Cplex网上版本差异很大这里只讲我验证过的稳定路径。先确认Matlab版本和Cplex版本支持关系一般来说Cplex 12.9/12.10配合Matlab R2020以后版本问题不大。安装完Cplex后在Matlab里执行addpath(genpath(Cplex安装目录/cplex/matlab))添加路径然后运行cplex.getVersion验证接口是否启动。注意有时候提示找不到cplex类大概率是路径没加进去或者Java环境变量问题别急着重装Cplex先把路径设置成永久路径保存。分布在大陆的同学们如果遇到许可证问题请走学术许可证通道申请这是正规且长期稳定的做法我不建议去搜索任何非官方许可证方案那些东西容易把Cplex安装搞坏得不偿失。3.2 双层模型转单层KKT条件线性化这是全文最核心的工程步骤。下层问题是一个以电价为参数的优化问题如果下层是凸优化线性规划或二次规划就可以用KKT条件把下层问题“等价替换”为一组等式和不等式约束然后嵌入上层问题把双层模型变成一个单层混合整数优化问题。具体来说对于一个标准线性规划形式的下层问题min cxs.t. Ax ≤ b, x ≥ 0对应的KKT条件包括原问题可行Ax ≤ b, x ≥ 0对偶问题可行Ay ≤ c, y ≥ 0y为对偶变量互补松弛条件y(b - Ax)0或者写成逐项形式梯度条件一般隐式满足其中互补松弛条件是非线性的。工程上最常用的大M法把逐项互补条件线性化。对于一对互补变量u和v即满足u ≥ 0, v ≥ 0, u*v 0引入二进制变量z写成0 ≤ u ≤ M * z 0 ≤ v ≤ M * (1 - z)这样就得到一个混合整数规划MIP直接用Cplex求解。M的取值有讲究后面第5章专门说。需要强调KKT条件转换的前提是下层问题凸且满足约束规格比如Slater条件。我手头的微网模型里下层都是线性约束这个条件能满足。如果你的下层模型里带了非线性等式约束比如潮流方程那就不能直接套KKT需要先把潮流方程线性化或者改用其他方法。3.3 KKT转化的关键代码结构用Matlab调用Cplex求解MIP我喜欢把矩阵构建和求解分开这样查错方便% 变量编号约定x[上层连续变量; 下层连续变量; 对偶变量; 二进制辅助变量] nVar nUpper nLower nDual nBin; f zeros(nVar, 1); % 填目标函数系数注意从Cplex返回的结果索引要对齐 % 不等式约束 Aineq * x bineq Aineq sparse(...); bineq ...; % 等式约束 Aeq * x beq Aeq sparse(...); beq ...; % 变量边界 lb -inf*ones(nVar,1); ub inf*ones(nVar,1); % 变量类型连续为C二进制为B ctype char(ones(1,nVar)*C); ctype(binIdx) B; % 调用cplexmilp求解 options cplexoptimset(cplex); options.Display on; options.Timing 1; options.MaxTime 3600; [x, fval, exitflag] cplexmilp(f, Aineq, bineq, Aeq, beq, lb, ub, ctype, [], options);这段代码框架看起来简单但真正写起来比较费劲的是构造互补松弛的大M约束。我建议把互补松弛条件的索引单独写成一段函数输入是下层约束编号输出是对应的大M约束矩阵块这样比在主矩阵里手工拼要少错很多。调试的时候经常遇到“运行结果非整数但变量类型标了B”的问题这一般是变量类型数组没对齐注意检查ctype变量的长度是否和总变量数一致。3.4 迭代博弈法和单层KKT法的选择有些同学更喜欢用迭代法上层先给定电价下层求解得到用电量上层再根据用电量调整电价反复迭代直到电价不再变化。这个方法的优点是思路直接、代码量小适合上课演示但实际用起来有两个麻烦第一迭代未必收敛电价可能震荡尤其在储能参与时容易出现周期振荡第二即使收敛通常只能保证收敛到局部均衡你很难向审稿人或者客户证明你求出来的点是全局均衡。而KKT单层化之后用Cplex直接求解MIP理论上能得到全局最优解代价是模型规模显著增大尤其是二进制辅助变量数量等于互补松弛条件的个数微网规模一大变量数量轻轻松松上万。我的取舍原则是微网内用户数不超过8个、变量规模在可控范围时首选单层KKT方法如果规模实在太大再考虑能不能用迭代或者Benders分解。实际做学术仿真场景单层KKT是主流选择。3.5 Cplex求解参数的实用设置Cplex默认参数在大规模MIP上表现一般我总结出一套适合微网博弈模型的参数组合MIP容差默认1e-4即可不需要刻意收紧如果发现结果出现微小违反约束可以把options.mip.tolerances.mipgap设为1e-6再跑一次验证。时间限制建议设置MaxTime为15-30分钟不要一个模型无限跑下去审稿流程等不起。线程数个人电脑建议设4或6设满反而因为内存带宽问题让求解变慢。数值强调打开solutiontarget和numericalemphasis微网模型里SOC递推系数常常跨数量级小时级别和兆瓦级别混在一起数值问题很常见。还有个很实用的技巧把求解结果导出成诊断文件options.saverm 1Cplex会生成模型文件配合Cplex自带的交互界面可以放大查看哪个约束在作怪这个在调试复杂模型时几乎是救命稻草。4. 算例设计从均衡解到运行曲线4.1 典型算例结构我常用的验证算例这样设置一个单母线微网包含1台风电机组、1组光伏阵列、1台燃气轮机、1套运营商侧储能、3个产消者用户每户自带光伏其中1户还带储能微网通过联络线接入上级配电网。调度周期24小时单位时段1小时。分时电价由博弈内生决定不对用户强制定价。参数可以按下面这个框架设置具体数值看你的研究场景调整参数数值说明峰时段09:00-12:00, 18:00-22:00负荷高峰、电价可能上浮平时段07:00-09:00, 12:00-18:00过渡时段谷时段22:00-次日07:00负荷低谷、电价可能下探储能容量1000 kWh运营商/ 300 kWh用户初始SOC设为0.2储能功率上限200 kW运营商/ 50 kW用户充放电效率95%燃气轮机150-400 kW爬坡率60 kW/h燃料成本二次函数联络线-300 kW 到 300 kW负值为购电正值为售电光伏出力用典型日曲线按晴天场景确定用户可调负荷每户约30 kW可转移转移时段跨度为4小时4.2 求解流程的完整执行顺序我自己跑算例会按照固定顺序执行每一步都验证结果无误再进入下一步建上层目标函数和约束先临时忽略下层直接求解一个“纯单层经济调度”问题确认功率平衡和储能SOC递推没问题。加入下层模型但不加博弈下层作为参数输入比如电价固定验证用户响应逻辑。引入KKT条件替换下层并加入大M线性化生成完整的单层MIP模型。调用Cplex求解检查exitflag和目标函数值。根据均衡解计算合作剩余枚举联盟求解Shapley值。把所有结果整理成24小时运行曲线对比合作与非合作场景。这套流程实际操作下来每一步都能快速定位问题来源。我见过太多人直接跳到第3步结果模型报错之后根本分不清是上层约束问题还是KKT代码问题只能一遍一遍盲调。4.3 结果展示的三种典型图表结果通常展示三类内容内部分时电价曲线、储能充放电运行曲线、各主体收益对比表。内部分时电价曲线是最能体现博弈价值的结果。观察这个曲线你会发现运营商在光伏出力大的午间会把收购电价压低但不敢压太低因为电价太低会导致用户光伏不卖给运营商在晚高峰时段运营商倾向于提高售电价格但同样有约束因为电价过高会导致用户更多使用自建储能来供电运营商卖电量反而下降。这就是主从博弈均衡的直观体现很有解释力。储能曲线要重点看“电价低充电、电价高放电”的峰谷套利规律。如果发现储能和电价的互动不符合这个规律先查SOC递推公式和充电效率参数是否录入正确。收益对比表是最能打动非技术读者的部分主体单独运行收益合作Stackelberg收益合作剩余分配后收益提升比例微网运营商5200610056007.7%用户16806407205.9%用户27407107906.8%用户352056061017.3%注意这个表里有一行有意思用户1、用户2在合作Stackelberg场景下的“原始收益”反而低于单独运行这正是我前面提到的“单纯主从博弈收益分配不公”的体现。经过Shapley值重新分配后所有主体都实现了收益提升这个结果才是“合作型”的验证点。4.4 均衡解唯一性问题的处理很多人在分析结果时会忽略一个事MIP可能会有多个最优解也就是同一个目标函数值对应多组不同的电价和功率曲线。如果直接拿Cplex第一次求出来的解进入Shapley计算得到的分配结果可能不稳定。我建议在目标函数里加一个极小的二次正则项比如电价变化平方项的0.001倍让Cplex在多个最优解里有倾向性地选择更平滑的一组。加了正则项之后结果复现性会好很多图和数据也更好看。这个方法不影响最优值太多属于工程处理上值得做的事。5. 模型调试的完整链路与工程经验5.1 第一批跑出来的结果为什么我不信我自己第一次跑通完整模型时看到目标函数值比预期高了不少以为博弈模型带来巨大收益提升差点直接写进文章里。后来做敏感性分析才发现是因为我在Shapley值计算时忘了把储能维护成本纳入特征函数导致合作剩余被系统性高估。这件事给我留下的习惯是第一批计算出来的结果先不要急着分析第一件事是对着单独运行场景、非合作博弈场景、合作博弈场景三组结果做一遍合理性检查。比如储能SOC曲线必须落在[0.1,0.9]区间功率平衡约束在任意时段残差小于1e-6各主体收益加总应该等于微网总收益这个不等式核对如果对不上一定哪里编码错了。5.2 大M取值一个直接影响求解质量的细节大M取值这个细节直接决定互补松弛条件的线性化成不成立。M如果取得太小小于变量本身可能达到的最大值会错误地把某些可行解空间裁掉M如果取太大比如动不动就1e8Cplex求解过程的数值稳定性会大幅下降出现各种莫名其妙的不可行或者震荡。我的经验是先求出每个变量的合理边界。比如电价变量可以从主网电价和成本参数估算一个理论上下限设为[0.1, 2.0]元/kWh储能功率上限直接由P_max决定。然后M取这个边界上限的10倍左右比如电价上限是2M取20。对于SOC类变量上限就是1M取5就够了。不要对所有约束用同一个M分组设置数值稳定性会好很多。5.3 互补松弛条件线性化之后模型无解的排查思路如果你把KKT和大M都写对之后Cplex却报没有可行解最常见的两个原因第一大M取小了。检查方法把M放大10倍重新求解如果可行解出现说明M确实取小了。第二互补松弛条件里遗漏了对偶变量的非负约束。很多KKT表述里会把对偶变量非负性当作隐含条件但在代码里必须显式写进约束否则Cplex不受限制对偶变量可能取负值导致整个模型实际被错误约束撑爆。我自己的排查顺序是先去掉互补松弛条件但保留KKT其他约束求解看是否可行如果去掉后可行问题一定在互补松弛相关约束再逐个检查下半部分的松弛变量和大M约束是否正确通常能在30分钟内定位到问题。5.4 迭代解法在什么情况下仍然值得用虽然我推荐KKT单层化但迭代解法在某些场景仍然有意义。比如你只是想快速验证用户响应机制是否合理或者需要和分布式算法做对比迭代法的实现成本低很多。我常用的一种稳健迭代叫“不动点迭代 超松弛”把上层电价更新公式写成ρ_new (1-α)ρ_old αρ_response其中α取0.3左右可以有效避免电价震荡。收敛判据我不会只看电价变化绝对值还要看用户总用电量的变化。有时候电价已经收敛但用户侧光伏上网功率还在缓慢漂移这通常是因为下层模型里用了带滞后的目标函数导致对偶解不稳定。5.5 结果可解释性的三个验证标尺最后不管用什么方法求解我都会用三把尺子检查结果的可解释性第一功率平衡尺24个时段里每个时段的发电用电必须严格平衡误差超过1e-4就要追溯。 第二收益合理性尺各主体收益必须为正且合作分配的收益必须高于单独运行。出现负收益说明有约束或者分配方式出了问题。 第三价格响应尺电价升高时用户购电量应该下降或持平电价降低时用户购电应该上升或持平这个“需求随价格递减”的基本经济规律如果被违背大概率是Shapley或KKT代码里有逻辑错误。这三把尺子跑完都通过我才会认为这组结果是可用的、可以放进论文或项目报告里的。我一个人从头搭这套模型的时候最大的感受是博弈模型本身不复杂复杂的全在细节里。KKT条件、大M法、Shapley值这三块单独拿出来每一个都是经典知识点但真正把它们组合成一套能稳定运行、结果还能有解释力的微网运行策略需要反复调试的功夫远超预期。我建议你先做一个三用户单体微网的小算例把这三块分别调通再逐步扩展规模。最后分享一个小技巧Cplex装完以后每次打开Matlab先运行一句cplex.getVersion确认接口正常再开始写模型代码这个小习惯能帮你省掉很多“模型没问题但接口死活不响应”的诡异晚上。