电力仿真的核心问题是解大规模线性方程组——电网的节点导纳矩阵Y 矩阵是稀疏矩阵求解节点电压需要解 Y×V I。elec-ops-simulation 把原本跑在 CPU 上的电网仿真计算牛顿-拉夫逊法、快速解耦法搬到昇腾 NPU 上利用 NPU 的高性能矩阵运算能力加速计算。电网仿真的基本数学模型电网仿真分三个层次电网计算层次 ├─ 1. 潮流计算Power Flow │ 求解节点电压幅值和相角 │ 数学本质解非线性方程组牛顿-拉夫逊迭代 │ ├─ 2. 暂态稳定计算Transient Stability │ 求解发电机转子运动方程微分方程组 │ 数学本质时域积分隐式梯形法 / 龙格-库塔 │ └─ 3. 短路电流计算Short Circuit 求解故障后各支路电流 数学本质解线性稀疏矩阵LU 分解 前代/回代三者的算量差异计算类型矩阵规模单次求解复杂度迭代次数潮流计算10K×10KO(N^1.5)稀疏 Cholesky5-20 次暂态稳定50K×50KO(N^2)稠密矩阵时步 0.01s × 10s 1000 步短路电流100K×100KO(N^1.3)稀疏 LU1 次直接求解NPU 的优势在暂态稳定计算——矩阵规模大、反复求解、适合 Cube 单元。牛顿-拉夫逊法在 NPU 上的实现潮流计算的核心迭代公式ΔP P_measured - P_calculated(V, θ) ΔQ Q_measured - Q_calculated(V, θ) 修正方程 [∂P/∂θ ∂P/∂V] [Δθ] [ΔP] [∂Q/∂θ ∂Q/∂V] [ΔV] [ΔQ] → J × Δx Δf → Δx J⁻¹ × Δf雅可比矩阵 J 是 2N×2N 的稀疏矩阵N 是节点数。每次迭代需要计算 J → 解线性方程组 → 更新 V 和 θ → 检查收敛。// elec-ops-simulation/kernels/nr_iteration.cpp__aicore__voidNewtonRaphsonIteration(GlobalTensorfloatV,// 节点电压幅值 [N]GlobalTensorfloattheta,// 节点电压相角 [N]GlobalTensorfloatP_meas,// 有功功率测量值 [N]GlobalTensorfloatQ_meas,// 无功功率测量值 [N]GlobalTensorfloatY_real,// 导纳矩阵实部稀疏 CSRGlobalTensorintY_col_idx,// 列索引GlobalTensorintY_row_ptr,// 行指针intN,// 节点数floattolerance// 收敛判据){intiter0;floatmax_mismatch1e9;while(max_mismatchtoleranceiter50){// 第一步计算 P_calc 和 Q_calc前代计算// P_calc[i] V[i] × Σ_j (V[j] × |Y_ij| × cos(θ_i - θ_j - φ_ij))// Q_calc[i] V[i] × Σ_j (V[j] × |Y_ij| × sin(θ_i - θ_j - φ_ij))//// 这是 O(N^2) 的稠密计算——NPU Cube 单元的优势区间LocalTensorfloatP_calc(N);LocalTensorfloatQ_calc(N);ComputePowerInjections(V,theta,Y_real,Y_col_idx,Y_row_ptr,P_calc,Q_calc,N);// 第二步计算失配量 ΔP, ΔQLocalTensorfloatdP(N),dQ(N);Sub(dP,P_meas,P_calc);// Vector 单元逐元素Sub(dQ,Q_meas,Q_calc);// 收敛判据max(|ΔP|, |ΔQ|) tolerancemax_mismatchMax(Abs(dP),Abs(dQ));// 第三步计算雅可比矩阵 J稀疏// J_ij ∂P_i/∂θ_j, ∂P_i/∂V_j, ∂Q_i/∂θ_j, ∂Q_i/∂V_j// J 是稀疏的导纳矩阵 Y 稀疏 → J 稀疏SparseMatrixJ(2*N,2*N);ComputeJacobian(V,theta,Y_real,Y_col_idx,Y_row_ptr,J,N);// 第四步解线性方程组 J × Δx Δf// 用稀疏 LU 分解更适合 NPU 的 SDMA 流水线LocalTensorfloatdx(2*N);SparseLUSolve(J,Concat(dP,dQ),dx,2*N);// 第五步更新 V 和 θ// V[i] ΔV[i], θ[i] Δθ[i]Add(V,V,Slice(dx,N,2*N));Add(theta,theta,Slice(dx,0,N));iter;}}牛顿-拉夫逊迭代的瓶颈在第三步计算雅可比矩阵和第四步解线性方程组。雅可比矩阵的计算是 O(N²) 稠密——NPU 的 Cube 单元在这里大幅领先 CPU。稀疏 LU 分解的 NPU 优化短路电流计算只需求解一次线性方程组不需要迭代但矩阵规模大100K×100K。直接 LU 分解的复杂度是 O(N³)——即使稀疏矩阵fill-in填充元也会让有效非零元膨胀到 O(N^1.5) 或更高。NPU 上的稀疏 LU 分解用三步流水线// elec-ops-simulation/kernels/sparse_lu.cpp__aicore__voidSparseLUFactorization(GlobalTensorfloatL_val,// 下三角CSRGlobalTensorintL_col_idx,GlobalTensorintL_row_ptr,GlobalTensorfloatU_val,// 上三角CSRGlobalTensorintU_col_idx,GlobalTensorintU_row_ptr,constSparseMatrixA,// 输入导纳矩阵intN){// 第一步符号分析在 CPU 上做只做一次// 确定 LU 分解中 fill-in 的位置非零结构SymbolicAnalysis(A,L_row_ptr,U_row_ptr);// 第二步数值分解在 NPU 上做每次故障计算都要跑// 左看右看每个主元消元步骤是稀疏的 GEMMfor(intk0;kN;k){// 选取主元带Pivot 的稀疏消元intpivotSelectPivot(A,k,N);// 加载第 k 列的非零元到 L1LocalTensorfloata_k;SparseGatherColumn(A,k,pivot,a_k);// 归一化L_ik A_ik / U_kk// 逐元素除法Vector 单元Div(a_k,a_k,A(k,k));// 秩 1 更新A_ij - L_ik × U_kj// 稀疏 GEMMCube 单元// 只有 a_k 的非零位置对应的行需要更新SparseRankOneUpdate(A,a_k,k,pivot);}// 第三步前向代入解 Ly b和后向代入解 Ux y// 稀疏三角求解复杂度 O(nnz(L) nnz(U))LocalTensorfloaty(N),x(N);SparseForwardSolve(L_val,L_col_idx,L_row_ptr,b,y,N);SparseBackwardSolve(U_val,U_col_idx,U_row_ptr,y,x,N);}稀疏 LU 分解的关键优化fill-in 控制重新排序节点编号减少非零元注入和稀疏 GEMMCube 单元处理非零块矩阵乘。暂态稳定的时域积分暂态稳定计算是电网仿真中最重的部分——需要以 0.01s 的步长模拟 10s 的故障过程每一步都要解一次微分方程组。发电机转子运动方程dδ/dt ω - ω₀ dω/dt (P_m - P_e - D(ω - ω₀)) / M其中 δ 是转子角ω 是转速P_m 是机械功率P_e 是电磁功率需要从潮流计算得来。// elec-ops-simulation/kernels/transient_sim.cpp__aicore__voidTransientSimulator(GlobalTensorfloatgen_delta,// 发电机转子角 [N_gen]GlobalTensorfloatgen_omega,// 发电机转速 [N_gen]GlobalTensorfloatbus_voltage,// 母线电压 [N_bus]floatdt,// 时步0.01sfloattotal_time,// 总模拟时间10sintN_gen,intN_bus){intnum_stepsint(total_time/dt);// 1000 步for(intt0;tnum_steps;t){// 第一步用当前电压计算电磁功率 P_e// P_e Σ_j (V_i × V_j × |Y_ij| × cos(δ_i - δ_j - φ_ij))// 每个发电机的 P_e 需要全电网电压 → 稠密矩阵乘LocalTensorfloatP_e(N_gen);ComputeElectricalPower(gen_delta,bus_voltage,Y_bus,P_e,N_gen,N_bus);// 第二步时域积分隐式梯形法// 梯形法x_{n1} x_n (dt/2) × (f(x_n) f(x_{n1}))// 需要解隐式方程 → 每个时步内迭代 3-5 次for(intinner0;inner5;inner){// 预测步显式欧拉LocalTensorfloatdelta_pred(N_gen),omega_pred(N_gen);EulerStep(gen_delta,gen_omega,P_e,delta_pred,omega_pred,dt,N_gen);// 校正步用预测值算 f(x_{n1})LocalTensorfloatP_e_pred(N_gen);ComputeElectricalPower(delta_pred,bus_voltage,Y_bus,P_e_pred,N_gen,N_bus);// 梯形法合并Axpy(gen_delta,0.5f*dt,Add(f(x_n),f(x_{n1})),gen_delta,N_gen);Axpy(gen_omega,0.5f*dt/M,Add((P_m-P_e)/M,(P_m-P_e_pred)/M),gen_omega,N_gen);}// 第三步更新母线电压解线性方程组// 每个时步都要解一次Y_bus × V IComputeBusVoltage(gen_delta,gen_omega,Y_bus,bus_voltage,N_bus);}}每个时步内的主要计算量P_e 计算O(N²) 稠密矩阵乘 线性方程组求解稀疏 LU。1000 个时步加在一起总计算量达到 10⁸ 量级——CPU 上需要数小时NPU 上压缩到分钟级。踩坑一导纳矩阵的重排序优化电网的导纳矩阵Y 矩阵的非零结构取决于电网的拓扑连接。糟糕的节点编号顺序会产生大量 fill-inLU 分解时原本为零的位置变成非零稀疏性被破坏。错误直接用原始节点编号做 LU 分解。// 原始节点编号按地理位置编号// 节点 1-100城区电网// 节点 101-500郊区电网// 节点 501-1000远郊电网// 三个区域之间只有 2-3 条联络线// → Y 矩阵是近乎分块对角的 极少的非零元在块外//// 但 LU 分解时这 2-3 条联络线会导致 fill-in 蔓延到整个矩阵// nnz 从 5000 膨胀到 500000100×正确用 AMDApproximate Minimum Degree算法重排序节点编号最小化 fill-in。// 重排序把联络线「压缩」到矩阵右下角// fill-in 被限制在很小的子矩阵里Permutation pAMDOrdering(Y_bus);Y_bus_permutedApplyPermutation(Y_bus,p);// LU 分解后的 nnz// 原始500000// 重排序后8000减少 98.4%踩坑二暂态稳定计算的数值发散隐式梯形法本身是数值稳定的A 稳定性但实现时如果 P_e 计算精度不够FP16 量化误差积累 1000 个时步后数值发散——发电机转速 ω 变成 NaN。根因FP16 的动态范围不够最大值 65504而 P_e 计算的中间结果V_i × V_j 的乘积可能超过这个范围如果电压基准值设得太大。正确做法暂态稳定计算用 FP32不用 FP16 量化。// 错误整个 pipeline 用 FP16// 电压基准值 1000 kV → V_i 1.5标幺值// V_i × V_j 2.25 → 正常// 但 P_e 计算中用到了 Y_ij 的幅值可能 100 Siemens// 2.25 × 100 225 → FP16 可以表示// 问题在累加Σ_j 有 500 项 → 部分和可能超过 65504 → 溢出// 正确// 牛顿-拉夫逊迭代FP32精度要求高// 暂态稳定时域积分FP32数值稳定性要求// 短路电流计算FP16 可以但要检查动态范围踩坑三NPU 上的稀疏格式选择稀疏矩阵有三种常用格式CSRCompressed Sparse Row、CSCCompressed Sparse Column、COOCoordinate。NPU 的 SDMA 单元对 CSR 格式的访问模式最友好按行连续存储但有些算子比如稀疏矩阵转置用 CSC 更快。踩坑所有算子都用 CSR 格式。// 稀疏矩阵乘法 C A × B// A 是 CSR行主 按行连续 → Cube 单元读取 A 的行快// B 需要按列访问 → CSR 不连续要转成 CSC// 每次矩阵乘都隐含一次 CSR→CSC 转换 → 额外 O(nnz(B)) 开销正确根据算子特性选格式。// 导纳矩阵 Y 同时存 CSR 和 CSC 两份// 占用 2× 存储空间但避免运行时转换开销Y_csrToCSR(Y);Y_cscToCSC(Y);// 离线转换只做一次// 矩阵乘A_csr × B_csc → 不需要运行时转换// 前代/回代解线性方程组只需要 A_csr按行访问// 矩阵转置只需要 A_csc按列访问elec-ops-simulation 把电网仿真从「离线计算」推进到「在线仿真」——故障发生后数秒内完成仿真给调度中心提供决策依据。NPU 的价值不只是训练大模型也在这些传统 HPC高性能计算场景里提供实用加速。