前两天帮师弟调一套配电网日前调度的MATLAB代码核心是用CPLEX求解二阶锥规划模型把风电Wind、并联电容器组CB、静止无功发生器SVG、有载调压变压器OLTC和储能系统ESS全部塞进24小时联合优化。问题本身不算复杂但来来回回折腾了快一天最后发现是数据维度没对齐导致的报错。这让我想起自己第一次碰这类模型时的情形——光是把每个设备该写到哪个约束、为什么用二阶锥而不是直接解非线性规划就花了好几周。所以我决定把这套题的建模思路、代码实现和踩过的坑整理出来既给师弟留个备忘也方便正被主动配电网调度、电压优化、储能运行这些问题折磨的同学参考。这篇文章不会只贴一堆代码完事而是把每个环节的“为什么”讲清楚。你看完以后应该能自己把模型搭起来而不是复制别人的代码然后祈祷它能跑。1. 先搞明白问题这套模型到底在优化什么1.1 需求拆解多设备协同的日前调度先说题目里的几个缩写。Wind是分布式风电通常作为有功源接入配电网某个节点CB是并联电容器组分组投切提供无功补偿SVG是静止无功发生器能连续调节无功响应速度快OLTC是有载调压变压器通过改变变比来调节电压ESS是储能系统既能充电也能放电具备时间上的转移能力。为什么要把这五类设备放在一起优化传统配电网里无功补偿和电压调节基本是各管各的电容器手动投切变压器靠调度员远程调节。但分布式电源大量接入以后配电网从被动单向网络变成了主动网络风电出力在一天里波动很大如果只用单一设备去应对很容易出现电压越限、无功倒送、设备动作过于频繁等问题。这套模型解决的正是这个问题在日前提前一天根据预测数据决定未来24小时每个时段各设备的运行状态。这个“24小时”很关键因为ESS的荷电状态会跨时段延续OLTC和CB的动作次数也受机械寿命限制不能只看单断面所以必须把时间耦合放进模型。1.2 为什么选二阶锥规划而不是非线性规划配电网潮流本质上是非线性的。经典DistFlow支路潮流方程里电压降和电流平方项会产生非凸约束直接建模成混合整数非线性规划MINLPCPLEX这类商业求解器帮不上忙自己写算法又很难保证收敛到全局最优。二阶锥规划的思路是引入电压幅值平方U和电流幅值平方L两个替代变量把原本的非凸等式约束松弛成一个凸锥约束。具体来说支路潮流中有一条约束是U_i × L_ij P_ij² Q_ij²这是一个旋转二次曲面非凸。二阶锥松弛把它放宽成U_i × L_ij ≥ P_ij² Q_ij²这个不等式可以用一个标准的旋转锥约束表达而旋转锥是凸的CPLEX可以直接处理。你可能会问把等式放宽成不等式结果还可靠吗这里有个业内都知道的结论对辐射状配电网在目标函数是网损最小、购电成本最小这类单调函数的情况下松弛通常是“紧”的意思是求解结果中这个不等式会自然取等号因此解和原问题实际上一致。打个比方如果某样东西越少越好你又给它设了一个下限那它就总会贴着下限走不会放任自己变大。这也是为什么这套框架在配电网优化里越来越流行——它有全局最优保证求解快还能塞进整数变量CB、OLTC论文里写出来也站得住脚。2. 五大设备逐一建模细节都在这里2.1 风电与负荷数据先把输入做对模型跑不跑得通一半以上取决于输入数据组织得对不对。这里说的输入主要是节点负荷和风电预测出力。我给你一个可以“抄作业”的构建方式。假设用IEEE 33节点配电网做算例先把24小时负荷组织成一个矩阵行是时段、列是节点每个元素表示该节点该时段的有功负荷。实际算例里通常还会给无功负荷一般按功率因数0.85到0.95折算即可。构造方式很简单从典型日负荷曲线峰平谷三段提取24个标幺值再乘上每个节点的峰值负荷。风电方面新手最容易犯的错是把风电预测值直接当作必须发出的固定功率。更合理的做法是把它当成“最大可发功率”引入可弃风变量约束写成0 ≤ P_wind(t) ≤ P_wind_forecast(t)这样模型可以在系统不需要那么多电时主动弃掉一部分风电同时在目标函数里给弃风加一个惩罚项。这个处理既贴近实际风电场运行又不会因为强制消纳导致电压越限或潮流不可行。如果你不想引入可弃风也可以直接把风电当成负的有功负荷但那样少了一个对照维度论文里没得写。2.2 CB与OLTC离散变量进入SOCP的两种方式CB的模型比较简单。电容器组按组投切每组容量Q_step固定投入组数是一个整数变量Q_cb(t) n_cb(t) × Q_step0 ≤ n_cb(t) ≤ n_cb_maxn_cb(t) ∈ Z这个约束进入模型后是线性的不会破坏SOCP的凸性只是让问题从SOCP变成混合整数SOCP即MISOCP。组数少的时候也可以把n_cb(t)展开成一组0-1变量这样更符合CPLEX处理二进制的习惯但组数多的时候直接用整数变量更省空间。OLTC的相对麻烦因为变比和电压是相乘关系直接乘就非线性了。工程上常见的处理方法是给变压器支路单独建模设变比tap(t)是离散档位取值范围比如0.9到1.1每档步长0.0125对应17档。写约束时可以写成U_sec(t) tap_ratio(t)² × U_pri(t)tap_ratio(t)是离散变量但它的平方依然是常数表里查出来的常数所以这个约束本质上是“常数×变量”的线性形式不会产生非凸项。还有一种做法是把OLTC等效成一个理想变压器串联一个短路阻抗在支路上增加一个虚拟节点两种方式数学上等价看你自己习惯哪种。2.3 SVG与ESS连续调节设备的约束边界SVG是这里面最简单的设备它就是节点上一个连续可调的无功源-Q_svg_max ≤ Q_svg(t) ≤ Q_svg_max你只需要把它放在节点无功平衡方程里不需要引入任何额外变量。有的模型还会给SVG加一个爬坡约束比如相邻时段无功变化率有限制实际工程里确实有这种需求但大多数论文算例里不写看场景需要。ESS的建模是重点。它有两个连续变量充电功率P_ch(t)和放电功率P_dis(t)再加一个状态变量SOC(t)。核心约束是时段间的能量递推SOC(t1) SOC(t) η_ch × P_ch(t) × Δt - P_dis(t) / η_dis × Δt这里η_ch和η_dis分别是充放电效率典型值在0.9到0.95之间。注意Δt的单位如果功率是kW、容量是kWhΔ t必须用小时比如1小时。很多人在这上面翻车计算出来SOC一会儿超上限一会儿为负全是单位换算问题。另外一个必须加的约束是“不能同时充放电”。从数学上看如果允许同时充放电模型会出现一个明显漏洞储能一边低价充电、一边按放电价卖电目标函数里形成虚假套利。解决办法是用两个二进制变量u_ch(t) u_dis(t) ≤ 10 ≤ P_ch(t) ≤ P_max × u_ch(t)0 ≤ P_dis(t) ≤ P_max × u_dis(t)这样充放电状态互斥物理上才说得通。2.4 目标函数不是只有网损一项很多初学者以为目标函数就是网损最小真算起来就会发现光最优化网损ESS基本不会动作OLTC和CB也不愿意频繁调节因为设备动作在目标里没有代价。真实可用的目标函数通常包含四部分向上级电网购电成本分时电价 × 根节点注入功率网损成本Σ r_ij × L_ij(t)这个量在SOCP变量里是线性的设备动作惩罚CB投切次数、OLTC挡位变化次数的绝对值之和弃风惩罚弃风量 × 惩罚系数写成表达式大致是min Σ_t [ price(t) × P_sub(t) c_loss × Σ_ij r_ij L_ij(t) c_tap × |tap(t1)-tap(t)| c_cb × |n_cb(t1)-n_cb(t)| c_wind_curtail × P_curtail(t) ]注意绝对值项不是线性函数需要引入辅助变量做线性化。比如用Δ_tap(t)表示相邻时段挡位差加两个不等式约束Δ_tap(t) ≥ tap(t1) - tap(t) Δ_tap(t) ≥ tap(t) - tap(t1)然后在目标函数里加c_tap × Δ_tap(t)。这个线性化技巧虽然基础但很重要写漏了模型就会变成非线性。权重系数怎么设我的经验是购电成本权重按实际电价网损成本给一个相对较小的值设备动作惩罚系数取一个相对值让模型不会因为频繁调设备而被惩罚也不会为了省一点点网损就疯狂动作。弃风惩罚系数一般设得比购电价高这样模型只在必要时才弃风。3. MATLABCPLEX实现过程从建模到求解3.1 技术选型为什么用YALMIP包一层直接调用CPLEX的MATLAB接口也能写但要把MISOCP问题转成CPLEX内部的稀疏矩阵结构写起来极其痛苦一个规模稍大的模型几百个变量手写矩阵容易把人写崩溃。YALMIP是一个MATLAB建模工具箱作用相当于“翻译官”你用高级语义把变量和约束声明出来它自动转换成求解器需要的标准格式再传给CPLEX。这种方式维护起来方便改模型也灵活论文里需要反复调约束的时候优势尤其明显。安装组合可以参考MATLAB R2020a及以上 YALMIP release 2021 CPLEX 12.10或20.1。装完之后在MATLAB里运行“yalmiptest”看到CPLEX被识别为可用求解器就说明环境没问题。在YALMIP中声明优化问题只需要三样东西sdpvar连续变量、binvar/intvar二进制/整数变量、约束列表F然后调用optimize(F, obj, ops)。3.2 核心代码骨架变量声明、约束组装与求解下面是这套模型的核心代码骨架我做了简化但结构是完整的。假设节点数N_bus支路数N_branch时段数T24。% 1. 变量声明 U sdpvar(N_bus, T, full); % 节点电压幅值平方 L sdpvar(N_branch, T, full); % 支路电流幅值平方 Pij sdpvar(N_branch, T, full); % 支路有功 Qij sdpvar(N_branch, T, full); % 支路无功 tap intvar(N_tap, T); % OLTC档位整数 n_cb intvar(N_cb, T); % CB投入组数整数 P_ch sdpvar(N_bus, T, full); % 储能充电功率 P_dis sdpvar(N_bus, T, full); % 储能放电功率 SOC sdpvar(N_bus, T, full); % 储能荷电状态 u_ch binvar(N_bus, T); % 充电状态 u_dis binvar(N_bus, T); % 放电状态 % 2. 目标函数 obj 0; for t 1:T obj obj price(t) * P_sub(t); % 购电成本 obj obj c_loss * sum(r_ij .* L(:,t)); % 网损 obj obj c_tap * sum(Delta_tap(:,t)); % 挡位动作惩罚 obj obj c_cb * sum(Delta_cb(:,t)); % 电容器投切惩罚 obj obj c_w * sum(P_curtail(:,t)); % 弃风惩罚 end % 3. 二阶锥约束 F []; for t 1:T for ij 1:N_branch % ||[2P; 2Q; U_i - L_ij]|| U_i L_ij F [F, cone([2*Pij(ij,t); 2*Qij(ij,t); ... U(from(ij),t) - L(ij,t)], ... U(from(ij),t) L(ij,t))]; end end % 4. 节点功率平衡、OLTC约束、ESS约束... % 这里省略按2.1-2.4节的公式展开 % 5. 求解 ops sdpsettings(solver, cplex, verbose, 2, ... cplex.mip.tolerances.mipgap, 0.001); result optimize(F, obj, ops); % 6. 输出 if result.problem 0 U_val value(U); Pij_val value(Pij); % 画图、导出表格... else disp(求解失败 result.info); end重点说con ()那一段。YALMIP里con e(x, y)表示约束‖x‖≤y而旋转锥约束U_i × L_ij ≥ P_ij² Q_ij²可以通过代数变换写成标准形式‖[2P_ij; 2Q_ij; U_i - L_ij]‖ ≤ U_i L_ij这段变换是很多人的知识盲区知道这个写法二阶锥约束就迎刃而解了。你可以自己验算一下两边平方左边展开是4(P²Q²)(U_i-L_ij)²右边是(U_iL_ij)²展开后约掉U_i²L_ij²就得到U_i × L_ij ≥ P_ij²Q_ij²。3.3 24小时数据组织与结果导出数据组织是实测中最花时间的部分。我的建议是写一个脚本文件单独处理数据不要什么都堆在主模型里。结构大致是Load_data.mat24×N_bus负荷矩阵Wind_data.mat24×1风电最大出力曲线标幺或实际值Price.mat24×1分时电价Network.mat节点编号、支路编号、电阻电抗、拓扑连接关系运行主程序之前先画一下负荷曲线和风电曲线肉眼确认趋势合理。这个习惯能省去很多排查故障的时间——我见过有人把负荷曲线乘错了系数电压约束怎么调都越限最后发现是数据放大100倍。求解完之后用value()提取变量整理成表格输出。可以在MATLAB里生成类似下面这个结果表时段根节点购电功率(kW)ESS放电(kW)CB投入组数OLTC档位最低节点电压(p.u.)132000430.982229800430.985..................导出表格以后再画三张图各时段设备出力堆叠图、24小时电压包络图、储能SOC曲线。这三张图基本是这类论文的标配。4. 常见问题与排查技巧实录4.1 求解器直接报infeasible怎么找原因我做仿真最怕见到“infeasible problem”因为它不告诉你是哪条约束出了问题。根据经验90%的情况出在以下三处一是变量维度不匹配。YALMIP里如果某个变量矩阵声明成N_bus×T约束里却用到了N_branch×T的维度模型会直接变成inf feasible而且报错信息很难定位。排查方法是把每个变量在用到之前的size()函数打出来核对一遍。二是约束写出来本身矛盾。比如OLTC变比范围太窄节点电压上下限又卡得很死两个约束叠加导致没有任何可行解。这个也很常见风大发时段电压被风电顶上去OLTC又压不下来就成这样了。三是储能SOC末时段约束太苛刻。如果强制要求SOC(24) SOC(1)但充放电效率和功率限制使得一天内“充回来”的量达不到初值模型就不可行。解决办法是把末时段约束放宽为SOC(24) ≥ SOC_lower或者给SOC一个较小的偏差范围。定位infeasible的方法我用的最多的是“注释法”把约束列表从头开始每跑一次加一组约束看到底哪组约束加入后问题变成不可行。虽然笨但有效。“注释法”虽然慢但配合YALMIP的diagnostics输出基本一小时之内能找到问题。4.2 二阶锥松弛不紧结果不可靠怎么办这是个学术上更敏感的问题。虽然辐射状配电网的SOCP松弛通常紧但不是所有场景都保证紧。判断方法很简单求解后计算所有支路的残差res_ij(t) U_i(t) × L_ij(t) - P_ij(t)² - Q_ij(t)²如果这个值都接近0比如小于1e-4说明松弛紧解是可信的。如果某些时段残差明显大于0说明模型出现了松弛间隙结果不等价于原问题。松弛不紧的常见原因目标函数里网损权重太小或者罚项权重设置导致没有动力把潮流推向边界。处理办法有几个把锥约束右侧的U_iL_ij稍微乘以一个略大于1的系数或者给松弛变量加一个小惩罚项又或者调整目标函数权重让网损在目标里占更大比例。要注意的是这个问题不能靠“把锥约束改回等式”来解决因为那会重新引入非凸性。正确做法是调整罚项后重算并在论文里附上残差统计表让审稿人看到你确实检查过松弛精度。4.3 求解太慢MISOCP整数变量太多了这套模型引入整数变量的地方有CB和OLTC如果每个时段都用整数变量总整数变量数是T × (N_cb N_oltc)24时段、CB接3个节点、OLTC有2台那就是24×5120个整数变量再加上ESS的充放电状态u_ch和u_dis又是48个二进制变量。CPLEX处理一百多个整数变量的MISOCP通常几十秒能出解但如果你加的CB组数特别多或者把设备数量扩大求解时间会明显上升。我的经验是三个优化方向设相对MIPGap比如0.001或0.0001不一定非要0。对工程调度来说gap在0.1%以内完全可接受。提供初始可行解。先用连续松弛版本所有整数变量放松跑一遍把结果里整数变量四舍五入后作为整数解的初始值传给CPLEX。减少整数变量数量。比如CB不按每个节点单独的组数变量而是用等值容量阶梯方式把多组投切等价成一个离散档位变量。sdpsettings里可以这样设置ops sdpsettings(solver,cplex,verbose,2, ... cplex.mip.tolerances.mipgap,0.001, ... cplex.mip.tolerances.integrality,1e-5);另外提醒一点CPLEX的MISOCP求解器中连续SOCP部分是内点法求解的warm start不一定总是有效。如果初始解给得不好反而可能拖慢速度所以初始解策略要自己对比测试。4.4 版本兼容与环境问题YALMIP和CPLEX的版本兼容是个老生常谈的问题但真的很影响体验。最常见的是YALMIP报“Unknown solver cplex”或者“No suitable solver for this problem”通常原因是CPLEX的MATLAB接口没有被正常加载。解决办法先确认你安装的CPLEX版本里有cplexinteractive和对应MATLAB接口文件夹然后在MATLAB里运行cplex.setupWindows下一般是setup_cplex.m。如果还不行检查CPLEX版本和MATLAB版本的对应关系。比如R2021b配CPLEX 20.1比较稳R2018a则建议配CPLEX 12.9。还有一个细节如果你的MATLAB是64位CPLEX也要装64位版本32位接口装上去能识别但会崩溃。这个坑我踩过一次报警方式极其诡异——某些模型能解某些模型一跑就闪退。5. 扩展思路与实操体会5.1 从确定性模型走向不确定性优化这套模型里所有的风电和负荷都用的是预测值属于确定性优化。如果要做更深入的研究可以考虑两个方向的扩展。一个是鲁棒优化。把风电预测误差建模成不确定集合比如盒式集合或椭球集合然后求解鲁棒问题。这里要注意鲁棒对偶转换后模型可能会变成双层结构处理起来比SOCP复杂很多但配电网这种规模较小的网络通常还能应付。另一个是随机规划。用蒙特卡洛生成若干风电场景每个场景对应一组约束目标函数变成所有场景的期望成本。这个模型规模会成倍增长但CPLEX也能处理中小规模的随机MISOCP。不过我的建议是先把确定性模型吃透再把简单的不确定性分析加进去。很多人一上来就想做鲁棒结果对偶推导错误最后论文算例根本解释不通。基础模型跑通以后不确定性扩展就只是工作量问题而不是难度问题。5.2 我在反复调试这套模型后的几点体会第一一定要先跑单时段、再跑24时段。先固定t1把模型简化成单断面优化确认潮流约束、各类设备约束没问题再扩展到全天。直接上24时段数据量一大报错以后根本分不清是潮流写错还是时序约束写错。第二量纲检查要养成习惯。电压标幺值、功率kW/MW、容量kWh这三者之间经常出现数量级不一致的情况。比如电压幅值平方U在标幺值下大约是1.0而功率可能是几千kW二阶锥约束里U_i×L_ij这个量纲要跟P²Q²匹配数值比例差距太大会让求解器精度崩溃。建议功率全部用标幺值或者全部用kW和kVar统一不要在同一个模型里混用。第三CB、OLTC这类设备的动作惩罚系数不能设太大也不能设太小。设太大设备一动不动电压调节能力闲置设太小设备每个时段都在动机械寿命直接报废。我一般会跑一组参数扫描看看动作频率随系数的变化曲线选一个合理的转折点。第四这套模型对论文写作帮助很大。MISOCP有全局最优解不像启发式算法一样需要反复调参数才能保证解的质量审稿人对这个模型通常比较认可。而且只要把“二阶锥松弛的紧性验证”和“与MINLP方法的对比”两件事做了文章的理论完整性一下就上去了。我印象最深的一次调试经历是某次把所有约束都加齐之后模型怎么都不可行我把所有约束逐个注释折腾了四个小时最后发现是储能SOC递推公式里效率η放错了位置——应该是充电时乘η、放电时除以η我写反了。这个错误用眼睛根本看不出来只有检查约束量纲或者看SOC曲线才会发现SOC一直在缓慢下降怎么充都充不满。所以如果你拿着代码怎么调都不对不妨先把储能那一段单独拎出来设一个最简单的场景给定初始SOC和固定充放电功率看看SOC曲线是否符合物理直觉。这一类“单元测试”式的验证方法比死磕整个大模型要高效得多。下次你再拿到类似的配电网调度题目按这个顺序走先理清设备和数据再写模型约束接着用小规模算例验证最后再全时段跑结果。每一步都确认无误最后的结果基本不会出大问题。