尧图网站建设 
简介面向配电网优化研究者的MATLAB程序包聚焦基于二阶锥规划的主动配电网最优潮流求解。资源以IEEE33节点系统为算例在YALMIPCPLEX求解环境下实现WindCBSVGOLTCESS多设备多时段24h协同运行优化对应《主动配电网多源协同运行优化研究》等参考文献。程序提供骨灰级注释与代码复现说明便于初学者理解二阶锥建模、潮流约束松弛及求解器调用全过程也可作为主动配电网调度、电压无功优化等方向的教学与科研参考。压缩包共2个文件含1个m主程序与1个txt说明文档m主程序包含从数据输入、模型建立到结果输出的完整流程txt说明文档提供参考文献与运行注意事项整体仅4KB结构精简上手门槛低。已有37人学习浏览适合具备一定MATLAB与优化基础、希望快速复现典型算例的读者。资源虽小但完整覆盖从问题建模到结果输出的关键代码尤其适合用于课程设计、论文复现或算法对比验证。1. 二阶锥规划求解主动配电网最优潮流一次把IEEE33算例跑通的实战笔记拿到一套IEEE33节点主动配电网最优潮流算例代码里同时接了风电、电容器组、SVG、有载调压变压器和储能24小时多时段联合优化求解器用的YALMIPCPLEX。以前遇到这种问题第一反应是内点法硬解非线性潮流但OLTC档位是整数、储能充放电是0/1状态混合整数非线性规划能把人折腾到怀疑人生。换成基于二阶锥规划的SOCP松弛之后非凸潮流方程变成可解凸锥CPLEX几分钟内给全局解这套思路在这份MATLAB程序里落得很扎实。本文把锥松弛怎么建、元件约束怎么写、求解器参数怎么设、哪几个位置最容易翻车讲透适合做配电网调度、分布式资源接入、最优潮流方向的研究生和一线工程师。2. 最优潮流问题怎么变成二阶锥DistFlow松弛与紧性条件2.1 DistFlow潮流方程的困境电压平方项打破凸性配电网是辐射状拓扑Baran和Wu在1989年提出的DistFlow方程至今仍是主动配电网最优潮流的主流建模基础。对每条支路ij潮流方程写成节点注入功率平衡P_j(t) ΣP_ij(t) − ΣP_jk(t) r_ij·l_ij(t)Q_j(t) ΣQ_ij(t) − ΣQ_jk(t) x_ij·l_ij(t)支路电压降落v_j(t) v_i(t) − 2(r_ij·P_ij(t) x_ij·Q_ij(t)) (r_ij² x_ij²)·l_ij(t)这里的v_i是节点电压幅值平方l_ij是支路电流幅值平方。表面上看这些方程是线性的真正的麻烦在于支路功率与电流之间还有一个定义式P_ij² Q_ij² v_i·l_ij。这是一个典型的二次非凸约束正是它让整个最优潮流问题掉进非凸优化的坑里。直接在这个非凸空间里做牛顿法初值稍微给偏一点就容易收敛到局部最优拓扑稍微复杂一些雅可比矩阵还可能出现奇异。主动配电网场景更麻烦OLTC的档位是离散整数CB电容器组的投切是整组逻辑ESS储能的充放电是0/1状态。这些离散变量一旦和上面的非凸二次约束混在一起整个模型就是混合整数非线性规划MINLP。实际工程里几乎不会指望用一个通用求解器直接啃MINLP因为分支定界的下界松弛不够紧计算量根本扛不住。2.2 二阶锥松弛的数学写法与紧性条件解决办法是把非凸等式约束松弛成凸不等式。引入u_i v_i、L_ij l_ij两个辅助变量之后把P_ij² Q_ij² v_i·l_ij改写成P_ij² Q_ij² ≤ v_i·l_ij这个不等式等价于一个标准二阶锥约束|| (2P_ij, 2Q_ij, v_i − l_ij) || ≤ v_i l_ij数学上用欧几里得范数表示就是norm([2P; 2Q; v−l]) ≤ vl。之所以强调锥形式而不是等价的二次约束形式是因为主流求解器对二阶锥有专门的内点算法和预处理收敛行为和数值稳定性都好得多。这一点在YALMIP建模时直接关系到能不能被CPLEX识别后面避坑章会细说。松弛是否可行取决于一个关键性质松弛之后的最优解是否仍然落在原问题的等号面上。对辐射状配电网只要目标函数在支路电流上单调递增比如目标包含网损项Σr_ij·l_ij而且电压幅值约束在合理范围锥松弛就会被顶到边界上取等号。这个性质叫exact conic relaxation文献里对辐射状网络有严格证明。高红均的《主动配电网最优潮流研究及其应用实例》里也专门讨论过这个前提条件做复现时先看目标函数里有没有网损项基本能判断这套算例会不会收敛到紧解。2.3 元件模型加入后的约束形态变化把风电、CB、SVG、OLTC、ESS加进来之后整个模型的形态会发生变化。风电在配电网里通常按可调无功处理P_wind由预测给出Q_wind在功率因数允许范围内连续调节SVG本质是一个连续无功源输出直接作为变量CB是离散无功源一组一组的投切逻辑用整数变量表达。OLTC是最容易写错的地方。变压器变比k_t变化会直接乘进电压约束如果不做处理就会出现变量乘积项。常见做法是引入二进制档位变量通过one-hot编码把变比的平方写成线性组合。ESS则需要同时追踪SOC状态递推关系和充放电互斥逻辑否则解出来的SOC曲线可能同时充电又放电。这些元件约束本质上是把原来的纯连续非凸问题变成混合整数二阶锥规划MISOCPCPLEX的主分支引擎恰好能处理二进制变量加上锥优化内点法内核正好能胜任这种规模的问题。3. IEEE33节点24小时算例落地YALMIPCPLEX建模实战3.1 系统配置与决策变量定义IEEE33节点系统是配电网优化领域最常用的标准算例基准电压12.66kV基准功率10MW33个节点、32条支路带有联络开关但正常运行方式是辐射状。本算例的分布式资源布置方式参考了配电网多源协同运行的常见配置风电接入节点18CB和SVG放在节点22OLTC装在根节点与首端之间ESS放在节点33。先定义时间尺度和变量矩阵。YALMIP里sdpvar声明的是连续决策变量binvar声明0/1决策变量整个模型按24个时段展开。T 24; % 时段数 N 33; B 32; % 节点数与支路数 u sdpvar(N, T, full); % 节点电压幅值平方 l sdpvar(B, T, full); % 支路电流幅值平方 Pij sdpvar(B, T, full); % 支路有功 Qij sdpvar(B, T, full); % 支路无功 Pgen sdpvar(N, T, full);% 节点注入有功 Qgen sdpvar(N, T, full);% 节点注入无功这里的u、l分别对应前面推导里的v_i和l_ij。sdpvar第三个参数full明确声明为完整矩阵避免YALMIP对矩阵形状产生歧义。Pgen和Qgen是广义节点注入功率根节点PCC的注入、风电出力、储能充放电最终都通过节点功率平衡方程汇入到这两个变量里。3.2 目标函数与潮流约束的YALMIP写法目标函数采用网损最小为主目标、附加惩罚项的加权形式。网损用r_ij·l_ij在支路上累加弃风惩罚和储能SOC偏差惩罚通过大系数压住。r_br 1e-3 * ones(B, 1); % 支路电阻标幺值 obj sum(sum(r_br .* l)); % 网损项 obj obj 1e3 * sum(sum(Pwind_curtail)); % 弃风惩罚 obj obj 100 * sum(sum((soc - 0.5).^2));% SOC偏移惩罚防止储能闲置网损项是锥松弛紧性的锚点不能省略。弃风惩罚系数取1e3是工程里常用的量级既能压制弃风量又不会反过来扭曲网损的目标主导地位。SOC偏移惩罚是为了让储能在电价平段不至于完全闲置100这个权重也是经验值如果目标里没这个项很多时候优化结果会把SOC钉在边界上。潮流约束是核心YALMIP里的cone()函数专门用来声明二阶锥Constraints []; for t 1:T % 支路电压降落方程 Constraints [Constraints, ... u(to, t) u(from, t) ... - 2*(r_br .* Pij(:,t) x_br .* Qij(:,t)) ... (r_br.^2 x_br.^2) .* l(:,t)]; % 二阶锥约束|| [2P; 2Q; u(i)-l] || u(i) l Constraints [Constraints, ... cone([2*Pij(:,t); 2*Qij(:,t); u(from,t) - l(:,t)], ... u(from,t) l(:,t))]; endcone(x, t)的语义是norm(x) ≤ t。把锥写成这种形式CPLEX会直接识别为二阶锥约束走锥优化内核如果误写成Pij.^2 Qij.^2 u.*l这种二次约束形式CPLEX会当成非凸二次约束拒绝求解。from和to是IEEE33支路表的首末端节点编号我这里直接用了向量索引真实场景里如果YALMIP对向量索引报错就拆成逐支路的循环写本质没有差别。3.3 多时段耦合约束OLTC、ESS、CB的写法OLTC是最容易翻车的地方。直接用变比平方参与运算会出现变量乘积必须做二进制展开。这里用9档位one-hot编码n_tap binvar(9, T, full); u_root sdpvar(T, 1); for t 1:T Constraints [Constraints, sum(n_tap(:,t)) 1]; % 每次只能选一档 tap_val 0.90 0.025 * (0:8); % 9档变比步长2.5% Constraints [Constraints, ... u_root(t) sum(tap_val.^2 .* n_tap(:,t))]; % 电压平方线性组合 endn_tap是9行24列的二进制矩阵每一列都是one-hot也就是说每个时段OLTC必须且只能落在一个档位。tap_val是9档变比向量范围0.90到1.10步长2.5%。由于n_tap是0/1变量tap_val.^2乘上n_tap之后自然得到该档位对应的电压平方整个过程是线性的。ESS建模分三块SOC状态递推、充放电功率上下限、充放电互斥逻辑。soc sdpvar(T1, 1); p_ch sdpvar(T, 1); % 充电功率 p_dch sdpvar(T, 1); % 放电功率 z_ch binvar(T, 1); % 充电状态 Constraints [Constraints, soc(1) 0.5]; % 初始SOC 50% for t 1:T % SOC递推 Constraints [Constraints, ... soc(t1) soc(t) 0.9*p_ch(t) - p_dch(t)/0.9]; % 充放电互斥用bigM Constraints [Constraints, 0 p_ch(t) 0.2*z_ch(t)]; Constraints [Constraints, 0 p_dch(t) 0.2*(1-z_ch(t))]; endSOC递推里充电效率0.9、放电效率0.9这是锂电池的常见假设0.2是单时段最大充放电功率单位是标幺值对应基准功率10MW下的2MW功率限值。z_ch是0/1变量把一个时段的充放电强行拆成互斥区间。bigM系数这里直接用功率上限0.2不要随手给个1e5的BigM过大的BigM会严重拖慢分支定界的收敛速度。CB投切建模思路和OLTC一致用整数变量表达组数即可不再单独列代码。把这些约束全部拼进Constraints集合最后一句ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 1e-4; optimize(Constraints, obj, ops);mipgap设到1e-4是精度和求解时间的折中点。设成0虽然更严谨但24时段MISOCP的分支搜索会成倍增加。4. 避坑记录SOCP最优潮流常见问题与排查办法4.1 求解结果不满足原始潮流方程锥松弛不紧现象优化目标数值异常小回代DistFlow方程时支路电压残差达到10⁻²级别潮流根本对不上。原因目标函数只写了弃风或储能成本没有网损项Σr·l锥约束缺少电流越大越不利的驱动力松弛面没有被顶到等号处。解决在目标函数里显式加入网损项并检查最优解处P_ij² Q_ij²与u_i·l_ij的相对残差。我一般用max(|P_ij² Q_ij² − u_i·l_ij| / (u_i·l_ij 1e-6))来判断小于1e-4就认为是紧的。4.2 CPLEX报错或求解卡死锥约束的建模方式现象YALMIP报Second order cone not supported或者CPLEX内部错误直接退出。原因把锥约束写成了P_ij² Q_ij² ≤ u_i·l_ij这种二次约束形式老版本CPLEX不认非凸二次约束而YALMIP又没帮你自动转成锥形式。解决统一用cone()函数显式建模并确认CPLEX版本在12.7以上。还有一个隐蔽坑如果YALMIP版本太老cone()可能被解析成非线性约束而不是锥约束可以在model export(Constraints, obj, ops)之后手动看model.K向量有没有锥的维度记录。4.3 储能SOC跳变与OLTC频繁动作现象SOC曲线锯齿状乱跳OLTC一天动作几十上百次明显不符合实际设备寿命要求。原因目标函数里没有状态变化惩罚。SOC在低电价时段充满、高电价时段放空是正常的但如果充放电损耗系数设得太小算法会利用SOC反复小幅波动的空隙蹭网损收益。解决在目标函数里加入SOC相对0.5的偏差惩罚同时给OLTC加动作次数约束用辅助变量和大M线性化Σ|k_t − k_{t−1}| ≤ N_maxN_max一般取5到10次。4.4 24小时变量维数大导致求解过慢现象24时段MISOCP一跑就是几十分钟分支树停不下来。原因二进制变量多one-hot编码的OLTC每组9个二进制变量在24个时段上就是216个0/1变量加上ESS充放电状态48个、CB投切变量上百个CPLEX默认mipgap0时分支压力很大。解决把sdpsettings里mipgap放宽到1e-4或1e-3工程上1e-3精度已经很好给BigM找紧上界而不是拍脑袋写大数CB投切改成整数变量而非独立二进制也能砍掉一整块分支维度。4.5 与牛顿法结果对不上平启动与基准值问题现象SOCP解的电压幅值在1.0pu附近波动很小但和MATPOWER牛顿法算出的潮流结果相差很大。原因根节点电压约束没加或者标幺值基准搞混了。IEEE33的基准电压是12.66kV基准功率10MW负荷和电源都要换算到标幺值任何一处基准不一致全网潮流都会偏。解决显式加u(1, t) 1约束根节点电压固定在1.0pu并统一检查负荷、风电出力、储能功率全部除以了10MW基准值。5. 进阶验证解完先别收工回代一次交流潮流SOCP和MINLP最大的区别在于SOCP解出来的是松弛问题的解必须验证它是否满足原始非凸方程。一个让人信服的验证办法是把最优解回代进DistFlow方程逐支路检查残差。% 取SOCP结果 u_val value(u); l_val value(l); Pij_val value(Pij); Qij_val value(Qij); % 逐支路回代 res_v zeros(B, T); for t 1:T for k 1:B res_v(k, t) u_val(from(k), t) - u_val(to(k), t) ... - 2*(r_br(k)*Pij_val(k, t) x_br(k)*Qij_val(k, t)) ... (r_br(k)^2 x_br(k)^2) * l_val(k, t); end end max_res max(abs(res_v(:))); fprintf(最大电压残差%.2e (pu)\n, max_res);当max_res小于1e-4时说明锥松弛是紧的每一个支路电压降落方程都被严格满足这套解可以直接当作真实可执行的调度方案使用。如果残差在10⁻³量级回去看目标函数是否包含网损项或者把网损权重从1.0逐步提到1.05再跑一次通常能把松弛压紧。另外一个更深的技巧是检查对偶变量。CPLEX求解锥优化时会输出每个锥约束的对偶乘子如果某个锥约束的对偶乘子接近0说明这个锥约束没有被激活松弛在那里是松的对应支路就是踢回非凸边界的候选对象。工程上我一般直接用回代残差定位因为它最直观。做这套算例复现时有一次我为了加快求解速度把网损项从目标函数里删了结果电压残差到了3×10⁻²曲线画出来一片平滑回代全崩。从那以后我每次跑完SOCP最优潮流都强制走一遍DistFlow回代校验残差不小于1e-4就坚决不收工。这套验证习惯救了我好几次希望帮到你。本文还有配套的精品资源点击获取 二阶锥规划求解主动配电网最优潮流：IEEE33算例实战与锥松弛要点
简介面向配电网优化研究者的MATLAB程序包聚焦基于二阶锥规划的主动配电网最优潮流求解。资源以IEEE33节点系统为算例在YALMIPCPLEX求解环境下实现WindCBSVGOLTCESS多设备多时段24h协同运行优化对应《主动配电网多源协同运行优化研究》等参考文献。程序提供骨灰级注释与代码复现说明便于初学者理解二阶锥建模、潮流约束松弛及求解器调用全过程也可作为主动配电网调度、电压无功优化等方向的教学与科研参考。压缩包共2个文件含1个m主程序与1个txt说明文档m主程序包含从数据输入、模型建立到结果输出的完整流程txt说明文档提供参考文献与运行注意事项整体仅4KB结构精简上手门槛低。已有37人学习浏览适合具备一定MATLAB与优化基础、希望快速复现典型算例的读者。资源虽小但完整覆盖从问题建模到结果输出的关键代码尤其适合用于课程设计、论文复现或算法对比验证。1. 二阶锥规划求解主动配电网最优潮流一次把IEEE33算例跑通的实战笔记拿到一套IEEE33节点主动配电网最优潮流算例代码里同时接了风电、电容器组、SVG、有载调压变压器和储能24小时多时段联合优化求解器用的YALMIPCPLEX。以前遇到这种问题第一反应是内点法硬解非线性潮流但OLTC档位是整数、储能充放电是0/1状态混合整数非线性规划能把人折腾到怀疑人生。换成基于二阶锥规划的SOCP松弛之后非凸潮流方程变成可解凸锥CPLEX几分钟内给全局解这套思路在这份MATLAB程序里落得很扎实。本文把锥松弛怎么建、元件约束怎么写、求解器参数怎么设、哪几个位置最容易翻车讲透适合做配电网调度、分布式资源接入、最优潮流方向的研究生和一线工程师。2. 最优潮流问题怎么变成二阶锥DistFlow松弛与紧性条件2.1 DistFlow潮流方程的困境电压平方项打破凸性配电网是辐射状拓扑Baran和Wu在1989年提出的DistFlow方程至今仍是主动配电网最优潮流的主流建模基础。对每条支路ij潮流方程写成节点注入功率平衡P_j(t) ΣP_ij(t) − ΣP_jk(t) r_ij·l_ij(t)Q_j(t) ΣQ_ij(t) − ΣQ_jk(t) x_ij·l_ij(t)支路电压降落v_j(t) v_i(t) − 2(r_ij·P_ij(t) x_ij·Q_ij(t)) (r_ij² x_ij²)·l_ij(t)这里的v_i是节点电压幅值平方l_ij是支路电流幅值平方。表面上看这些方程是线性的真正的麻烦在于支路功率与电流之间还有一个定义式P_ij² Q_ij² v_i·l_ij。这是一个典型的二次非凸约束正是它让整个最优潮流问题掉进非凸优化的坑里。直接在这个非凸空间里做牛顿法初值稍微给偏一点就容易收敛到局部最优拓扑稍微复杂一些雅可比矩阵还可能出现奇异。主动配电网场景更麻烦OLTC的档位是离散整数CB电容器组的投切是整组逻辑ESS储能的充放电是0/1状态。这些离散变量一旦和上面的非凸二次约束混在一起整个模型就是混合整数非线性规划MINLP。实际工程里几乎不会指望用一个通用求解器直接啃MINLP因为分支定界的下界松弛不够紧计算量根本扛不住。2.2 二阶锥松弛的数学写法与紧性条件解决办法是把非凸等式约束松弛成凸不等式。引入u_i v_i、L_ij l_ij两个辅助变量之后把P_ij² Q_ij² v_i·l_ij改写成P_ij² Q_ij² ≤ v_i·l_ij这个不等式等价于一个标准二阶锥约束|| (2P_ij, 2Q_ij, v_i − l_ij) || ≤ v_i l_ij数学上用欧几里得范数表示就是norm([2P; 2Q; v−l]) ≤ vl。之所以强调锥形式而不是等价的二次约束形式是因为主流求解器对二阶锥有专门的内点算法和预处理收敛行为和数值稳定性都好得多。这一点在YALMIP建模时直接关系到能不能被CPLEX识别后面避坑章会细说。松弛是否可行取决于一个关键性质松弛之后的最优解是否仍然落在原问题的等号面上。对辐射状配电网只要目标函数在支路电流上单调递增比如目标包含网损项Σr_ij·l_ij而且电压幅值约束在合理范围锥松弛就会被顶到边界上取等号。这个性质叫exact conic relaxation文献里对辐射状网络有严格证明。高红均的《主动配电网最优潮流研究及其应用实例》里也专门讨论过这个前提条件做复现时先看目标函数里有没有网损项基本能判断这套算例会不会收敛到紧解。2.3 元件模型加入后的约束形态变化把风电、CB、SVG、OLTC、ESS加进来之后整个模型的形态会发生变化。风电在配电网里通常按可调无功处理P_wind由预测给出Q_wind在功率因数允许范围内连续调节SVG本质是一个连续无功源输出直接作为变量CB是离散无功源一组一组的投切逻辑用整数变量表达。OLTC是最容易写错的地方。变压器变比k_t变化会直接乘进电压约束如果不做处理就会出现变量乘积项。常见做法是引入二进制档位变量通过one-hot编码把变比的平方写成线性组合。ESS则需要同时追踪SOC状态递推关系和充放电互斥逻辑否则解出来的SOC曲线可能同时充电又放电。这些元件约束本质上是把原来的纯连续非凸问题变成混合整数二阶锥规划MISOCPCPLEX的主分支引擎恰好能处理二进制变量加上锥优化内点法内核正好能胜任这种规模的问题。3. IEEE33节点24小时算例落地YALMIPCPLEX建模实战3.1 系统配置与决策变量定义IEEE33节点系统是配电网优化领域最常用的标准算例基准电压12.66kV基准功率10MW33个节点、32条支路带有联络开关但正常运行方式是辐射状。本算例的分布式资源布置方式参考了配电网多源协同运行的常见配置风电接入节点18CB和SVG放在节点22OLTC装在根节点与首端之间ESS放在节点33。先定义时间尺度和变量矩阵。YALMIP里sdpvar声明的是连续决策变量binvar声明0/1决策变量整个模型按24个时段展开。T 24; % 时段数 N 33; B 32; % 节点数与支路数 u sdpvar(N, T, full); % 节点电压幅值平方 l sdpvar(B, T, full); % 支路电流幅值平方 Pij sdpvar(B, T, full); % 支路有功 Qij sdpvar(B, T, full); % 支路无功 Pgen sdpvar(N, T, full);% 节点注入有功 Qgen sdpvar(N, T, full);% 节点注入无功这里的u、l分别对应前面推导里的v_i和l_ij。sdpvar第三个参数full明确声明为完整矩阵避免YALMIP对矩阵形状产生歧义。Pgen和Qgen是广义节点注入功率根节点PCC的注入、风电出力、储能充放电最终都通过节点功率平衡方程汇入到这两个变量里。3.2 目标函数与潮流约束的YALMIP写法目标函数采用网损最小为主目标、附加惩罚项的加权形式。网损用r_ij·l_ij在支路上累加弃风惩罚和储能SOC偏差惩罚通过大系数压住。r_br 1e-3 * ones(B, 1); % 支路电阻标幺值 obj sum(sum(r_br .* l)); % 网损项 obj obj 1e3 * sum(sum(Pwind_curtail)); % 弃风惩罚 obj obj 100 * sum(sum((soc - 0.5).^2));% SOC偏移惩罚防止储能闲置网损项是锥松弛紧性的锚点不能省略。弃风惩罚系数取1e3是工程里常用的量级既能压制弃风量又不会反过来扭曲网损的目标主导地位。SOC偏移惩罚是为了让储能在电价平段不至于完全闲置100这个权重也是经验值如果目标里没这个项很多时候优化结果会把SOC钉在边界上。潮流约束是核心YALMIP里的cone()函数专门用来声明二阶锥Constraints []; for t 1:T % 支路电压降落方程 Constraints [Constraints, ... u(to, t) u(from, t) ... - 2*(r_br .* Pij(:,t) x_br .* Qij(:,t)) ... (r_br.^2 x_br.^2) .* l(:,t)]; % 二阶锥约束|| [2P; 2Q; u(i)-l] || u(i) l Constraints [Constraints, ... cone([2*Pij(:,t); 2*Qij(:,t); u(from,t) - l(:,t)], ... u(from,t) l(:,t))]; endcone(x, t)的语义是norm(x) ≤ t。把锥写成这种形式CPLEX会直接识别为二阶锥约束走锥优化内核如果误写成Pij.^2 Qij.^2 u.*l这种二次约束形式CPLEX会当成非凸二次约束拒绝求解。from和to是IEEE33支路表的首末端节点编号我这里直接用了向量索引真实场景里如果YALMIP对向量索引报错就拆成逐支路的循环写本质没有差别。3.3 多时段耦合约束OLTC、ESS、CB的写法OLTC是最容易翻车的地方。直接用变比平方参与运算会出现变量乘积必须做二进制展开。这里用9档位one-hot编码n_tap binvar(9, T, full); u_root sdpvar(T, 1); for t 1:T Constraints [Constraints, sum(n_tap(:,t)) 1]; % 每次只能选一档 tap_val 0.90 0.025 * (0:8); % 9档变比步长2.5% Constraints [Constraints, ... u_root(t) sum(tap_val.^2 .* n_tap(:,t))]; % 电压平方线性组合 endn_tap是9行24列的二进制矩阵每一列都是one-hot也就是说每个时段OLTC必须且只能落在一个档位。tap_val是9档变比向量范围0.90到1.10步长2.5%。由于n_tap是0/1变量tap_val.^2乘上n_tap之后自然得到该档位对应的电压平方整个过程是线性的。ESS建模分三块SOC状态递推、充放电功率上下限、充放电互斥逻辑。soc sdpvar(T1, 1); p_ch sdpvar(T, 1); % 充电功率 p_dch sdpvar(T, 1); % 放电功率 z_ch binvar(T, 1); % 充电状态 Constraints [Constraints, soc(1) 0.5]; % 初始SOC 50% for t 1:T % SOC递推 Constraints [Constraints, ... soc(t1) soc(t) 0.9*p_ch(t) - p_dch(t)/0.9]; % 充放电互斥用bigM Constraints [Constraints, 0 p_ch(t) 0.2*z_ch(t)]; Constraints [Constraints, 0 p_dch(t) 0.2*(1-z_ch(t))]; endSOC递推里充电效率0.9、放电效率0.9这是锂电池的常见假设0.2是单时段最大充放电功率单位是标幺值对应基准功率10MW下的2MW功率限值。z_ch是0/1变量把一个时段的充放电强行拆成互斥区间。bigM系数这里直接用功率上限0.2不要随手给个1e5的BigM过大的BigM会严重拖慢分支定界的收敛速度。CB投切建模思路和OLTC一致用整数变量表达组数即可不再单独列代码。把这些约束全部拼进Constraints集合最后一句ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 1e-4; optimize(Constraints, obj, ops);mipgap设到1e-4是精度和求解时间的折中点。设成0虽然更严谨但24时段MISOCP的分支搜索会成倍增加。4. 避坑记录SOCP最优潮流常见问题与排查办法4.1 求解结果不满足原始潮流方程锥松弛不紧现象优化目标数值异常小回代DistFlow方程时支路电压残差达到10⁻²级别潮流根本对不上。原因目标函数只写了弃风或储能成本没有网损项Σr·l锥约束缺少电流越大越不利的驱动力松弛面没有被顶到等号处。解决在目标函数里显式加入网损项并检查最优解处P_ij² Q_ij²与u_i·l_ij的相对残差。我一般用max(|P_ij² Q_ij² − u_i·l_ij| / (u_i·l_ij 1e-6))来判断小于1e-4就认为是紧的。4.2 CPLEX报错或求解卡死锥约束的建模方式现象YALMIP报Second order cone not supported或者CPLEX内部错误直接退出。原因把锥约束写成了P_ij² Q_ij² ≤ u_i·l_ij这种二次约束形式老版本CPLEX不认非凸二次约束而YALMIP又没帮你自动转成锥形式。解决统一用cone()函数显式建模并确认CPLEX版本在12.7以上。还有一个隐蔽坑如果YALMIP版本太老cone()可能被解析成非线性约束而不是锥约束可以在model export(Constraints, obj, ops)之后手动看model.K向量有没有锥的维度记录。4.3 储能SOC跳变与OLTC频繁动作现象SOC曲线锯齿状乱跳OLTC一天动作几十上百次明显不符合实际设备寿命要求。原因目标函数里没有状态变化惩罚。SOC在低电价时段充满、高电价时段放空是正常的但如果充放电损耗系数设得太小算法会利用SOC反复小幅波动的空隙蹭网损收益。解决在目标函数里加入SOC相对0.5的偏差惩罚同时给OLTC加动作次数约束用辅助变量和大M线性化Σ|k_t − k_{t−1}| ≤ N_maxN_max一般取5到10次。4.4 24小时变量维数大导致求解过慢现象24时段MISOCP一跑就是几十分钟分支树停不下来。原因二进制变量多one-hot编码的OLTC每组9个二进制变量在24个时段上就是216个0/1变量加上ESS充放电状态48个、CB投切变量上百个CPLEX默认mipgap0时分支压力很大。解决把sdpsettings里mipgap放宽到1e-4或1e-3工程上1e-3精度已经很好给BigM找紧上界而不是拍脑袋写大数CB投切改成整数变量而非独立二进制也能砍掉一整块分支维度。4.5 与牛顿法结果对不上平启动与基准值问题现象SOCP解的电压幅值在1.0pu附近波动很小但和MATPOWER牛顿法算出的潮流结果相差很大。原因根节点电压约束没加或者标幺值基准搞混了。IEEE33的基准电压是12.66kV基准功率10MW负荷和电源都要换算到标幺值任何一处基准不一致全网潮流都会偏。解决显式加u(1, t) 1约束根节点电压固定在1.0pu并统一检查负荷、风电出力、储能功率全部除以了10MW基准值。5. 进阶验证解完先别收工回代一次交流潮流SOCP和MINLP最大的区别在于SOCP解出来的是松弛问题的解必须验证它是否满足原始非凸方程。一个让人信服的验证办法是把最优解回代进DistFlow方程逐支路检查残差。% 取SOCP结果 u_val value(u); l_val value(l); Pij_val value(Pij); Qij_val value(Qij); % 逐支路回代 res_v zeros(B, T); for t 1:T for k 1:B res_v(k, t) u_val(from(k), t) - u_val(to(k), t) ... - 2*(r_br(k)*Pij_val(k, t) x_br(k)*Qij_val(k, t)) ... (r_br(k)^2 x_br(k)^2) * l_val(k, t); end end max_res max(abs(res_v(:))); fprintf(最大电压残差%.2e (pu)\n, max_res);当max_res小于1e-4时说明锥松弛是紧的每一个支路电压降落方程都被严格满足这套解可以直接当作真实可执行的调度方案使用。如果残差在10⁻³量级回去看目标函数是否包含网损项或者把网损权重从1.0逐步提到1.05再跑一次通常能把松弛压紧。另外一个更深的技巧是检查对偶变量。CPLEX求解锥优化时会输出每个锥约束的对偶乘子如果某个锥约束的对偶乘子接近0说明这个锥约束没有被激活松弛在那里是松的对应支路就是踢回非凸边界的候选对象。工程上我一般直接用回代残差定位因为它最直观。做这套算例复现时有一次我为了加快求解速度把网损项从目标函数里删了结果电压残差到了3×10⁻²曲线画出来一片平滑回代全崩。从那以后我每次跑完SOCP最优潮流都强制走一遍DistFlow回代校验残差不小于1e-4就坚决不收工。这套验证习惯救了我好几次希望帮到你。本文还有配套的精品资源点击获取