当前位置：首页 > news >正文

基于无迹变换的电网概率潮流分析 MATLAB 实现

news 2026/5/5 1:39:44

基于无迹变换（Unscented Transformation, UT）的电网概率潮流分析 MATLAB 实现

一、整体思路（工程级）

随机输入变量（负荷、风电、光伏） ↓ 无迹变换（UT）生成 Sigma 点 ↓ 逐点执行确定性潮流计算（Newton‑Raphson） ↓ 加权统计输出变量（电压、功率） ↓ 概率指标（均值、方差、PDF、越限概率）

二、电网概率潮流数学模型

1、输入随机变量

x=[PL1,QL1,PW1,PPV1,… ]T \mathbf{x} = [P_{L1}, Q_{L1}, P_{W1}, P_{PV1}, \dots]^Tx=[PL1,QL1,PW1,PPV1,…]T

负荷：正态分布
风电：Weibull →Nataf 变换 → 高斯空间
光伏：Beta →Nataf 变换 → 高斯空间

2、输出变量

z=[V1,V2,…,Pline,Qline]T \mathbf{z} = [V_1, V_2, \dots, P_{line}, Q_{line}]^Tz=[V1,V2,…,Pline,Qline]T

三、MATLAB 主程序

主脚本：`ut_probabilistic_power_flow.m`

%% ===============================% 基于无迹变换的电网概率潮流分析% 适用：IEEE 标准算例（Matpower）% ===============================clear;clc;close all;%% 1. 载入电网数据（IEEE 14节点为例）mpc=loadcase('case14');baseMVA=mpc.baseMVA;%% 2. 定义随机变量（负荷 + 风电）% 格式：[均值, 标准差]% 假设：节点 9、14 为负荷随机波动% 节点 6 为风电机组random_vars=[1.0,0.05;% P9 负荷1.0,0.05;% Q9 负荷1.0,0.05;% P14 负荷1.0,0.05;% Q14 负荷0.8,0.2;% 风电 P6 (均值0.8pu，波动20%)];n=size(random_vars,1);% 随机变量维数kappa=1;% UT 缩放因子%% 3. 无迹变换 Sigma 点生成mean_x=random_vars(:,1);std_x=random_vars(:,2);cov_x=diag(std_x.^2);L=chol(cov_x,'lower');lambda=n+kappa;N_sigma=2*n+1;Wm=zeros(1,N_sigma);% 均值权重Wc=zeros(1,N_sigma);% 协方差权重sigma_pts=zeros(n,N_sigma);% 第0个 Sigma 点sigma_pts(:,1)=mean_x;Wm(1)=kappa/lambda;Wc(1)=Wm(1);% 其余 Sigma 点fori=1:nsigma_pts(:,i+1)=mean_x+sqrt(lambda)*L(:,i);sigma_pts(:,i+1+n)=mean_x-sqrt(lambda)*L(:,i);Wm(i+1)=0.5/lambda;Wm(i+1+n)=0.5/lambda;Wc(i+1)=0.5/lambda;Wc(i+1+n)=0.5/lambda;end%% 4. 确定性潮流计算nbus=size(mpc.bus,1);V_out=zeros(nbus,N_sigma);fprintf('开始 UT 概率潮流计算（共 %d 个场景）...\n',N_sigma);fork=1:N_sigma mpc_run=mpc;% ===== 注入随机功率 =====mpc_run.bus(9,3)=sigma_pts(1,k)*baseMVA;% P9mpc_run.bus(9,4)=sigma_pts(2,k)*baseMVA;% Q9mpc_run.bus(14,3)=sigma_pts(3,k)*baseMVA;% P14mpc_run.bus(14,4)=sigma_pts(4,k)*baseMVA;% Q14mpc_run.gen(3,2)=sigma_pts(5,k)*baseMVA;% 风电 P6% ===== 执行潮流 =====results=runpf(mpc_run);ifresults.successV_out(:,k)=results.bus(:,8);% 电压幅值elseV_out(:,k)=NaN;endend%% 5. 统计输出变量V_mean=zeros(nbus,1);V_cov=zeros(nbus,nbus);fori=1:nbusV_mean(i)=sum(Wm'.*V_out(i,:)');tmp=V_out(i,:)'-V_mean(i);V_cov(i,i)=sum(Wc.*tmp.^2);end%% 6. 结果展示fprintf('\n========== 概率潮流结果 ==========\n');fori=1:nbusfprintf('Bus %2d: 电压均值 = %.4f pu, 标准差 = %.4f pu\n',...i,V_mean(i),sqrt(V_cov(i,i)));end%% 7. 可视化figure;subplot(1,2,1)histogram(V_out(9,:),30,'Normalization','pdf');title('节点9电压概率密度');xlabel('电压(pu)');ylabel('概率密度');subplot(1,2,2)[f,x]=ksdensity(V_out(9,:));plot(x,f,'LineWidth',1.5);title('节点9电压PDF（核密度估计）');grid on;

四、非高斯风光出力处理

Weibull 风电 → 高斯空间（Nataf）

functionx_gauss=weibull_to_gauss(P_mean,P_std)% 将 Weibull 分布映射到标准正态空间c=(P_mean/P_std)^2;k=(0.9874/c)^0.25;lambda=P_mean/gamma(1+1/k);% 抽样 WeibullP_sample=wblrnd(lambda,k,[1000,1]);% Nataf 变换x_gauss=norminv(empirical_cdf(P_sample,P_sample));end