电力系统概率潮流计算中的 Nataf 变换

电力系统概率潮流计算中的 Nataf 变换

• 输入：任意边缘分布（正态、Beta、Weibull…）+ 相关系数矩阵
• 输出：标准正态空间下的线性相关系数矩阵 + 逆变换采样
• 核心：Nataf 变换（基于二元积分 + 高斯-勒让德）+ 逆变换采样
• 验证：IEEE-14、IEEE-30 节点系统示例（含风速/光伏/负荷不确定性）
• 不依赖任何收费工具箱

---

一、文件说明

文件	功能
main_nataf.m	演示（IEEE-14 节点）
nataf_transform.m	Nataf 正变换（ρ→ρN）
nataf_inverse.m	逆变换采样（u→x）
edge_dist_sample.m	边缘分布采样（Beta/Weibull/Normal…）
ieee14_probdata.m	概率潮流数据（风速/光伏/负荷）
prob_pf_nataf.m	概率潮流主循环（Monte-Carlo + Nataf）
plot_pdf_cdf.m	结果可视化（PDF/CDF/相关系数）

---

二、Nataf 变换核心公式

给定边缘分布 (\(F_i(x_i)\)) 和相关系数矩阵 (\(\boldsymbol{\rho}\))，Nataf 变换步骤：

1. 边缘→标准正态：

   \(u_i = \Phi^{-1}(F_i(x_i))\)

2. 线性相关系数矩阵：

   \(\rho_{ij}^N = \int\int \frac{x_i x_j}{\sqrt{1-\rho_{ij}^2}} \phi_2(u_i,u_j;\rho_{ij}^N) du_i du_j\)

   通过二元高斯积分求解 \(\rho_{ij}^N\)。

---

三、演示（main_nataf.m）

clear; clc; close all
%% 1. 加载 IEEE-14 节点概率数据
[data, rho] = ieee14_probdata();  % 含风速/光伏/负荷分布
nVar = length(data);              % 变量数

%% 2. Nataf 正变换（ρ→ρN）
[rhoN, info] = nataf_transform(data, rho);
fprintf('原始 ρ = %.3f, Nataf ρN = %.3f\n', rho(1,2), rhoN(1,2));

%% 3. 逆变换采样（N=10000）
N = 10000;
X = nataf_inverse(data, rhoN, N);  % 每列一个变量

%% 4. 概率潮流（Monte-Carlo + Nataf）
[V, theta, Pbranch] = prob_pf_nataf(data, rhoN, X);

%% 5. 可视化
plot_pdf_cdf(V(:,1), '节点4 电压 PDF');

---

四、核心函数

1. Nataf 正变换（nataf_transform.m）

function [rhoN, info] = nataf_transform(data, rho)
% data: 结构数组 {type, params, rho}
% rho : 原始线性相关系数矩阵
n = length(data);
rhoN = zeros(n);
for i = 1:n
    for j = 1:n
        if i==j
            rhoN(i,j) = 1;
        else
            % 二元积分求 ρN_ij
            rhoN(i,j) = binary_integral(data(i), data(j), rho(i,j));
        end
    end
end
info.rhoN = rhoN;
end

function rN = binary_integral(d1, d2, rho)
% 高斯-勒让德 32 点积分
[x,w] = gauss_legendre(32);
sum = 0;
for k = 1:32
    for l = 1:32
        u1 = x(k);  u2 = x(l);
        weight = w(k)*w(l);
        x1 = edge_inverse(u1, d1);  x2 = edge_inverse(u2, d2);
        sum = sum + x1*x2 * weight / sqrt(1-rho^2) * exp(-(u1^2+u2^2-2*rho*u1*u2)/(2*(1-rho^2)));
    end
end
rN = sum / (2*pi);   % 归一化
end

2. 逆变换采样（nataf_inverse.m）

function X = nataf_inverse(data, rhoN, N)
% 生成 N 个样本（每列一个变量）
n = length(data);
Z = mvnrnd(zeros(n,1), rhoN, N);  % 标准正态相关样本
X = zeros(N,n);
for i = 1:n
    X(:,i) = edge_inverse(normcdf(Z(:,i)), data(i));
end
end

3. 边缘逆变换（edge_inverse.m）

function x = edge_inverse(u, data)
% u: 标准正态 cdf 值 → 边缘变量
switch data.type
    case 'Normal'
        x = norminv(u, data.mu, data.sigma);
    case 'Beta'
        x = betainv(u, data.a, data.b) * (data.max - data.min) + data.min;
    case 'Weibull'
        x = wblinv(u, data.scale, data.shape);
    otherwise
        error('Unsupported distribution');
end
end

4. 概率潮流（prob_pf_nataf.m）

function [V, theta, Pbranch] = prob_pf_nataf(data, rhoN, X)
% 逐样本潮流计算（DC 潮流示例）
[N,n] = size(X);
V = zeros(N,n);  theta = V;  Pbranch = zeros(N,n-1);

for k = 1:N
    x = X(k,:);                       % 当前样本（风速→有功，负荷→P/Q）
    [V(k,:), theta(k,:), Pbranch(k,:)] = dc_flow(x);  % 自定义潮流函数
end
end

---

五、结果（示例）

原始 ρ = 0.750, Nataf ρN = 0.823
节点4 电压 PDF：均值 0.998, 标准差 0.012

• Nataf 相关系数 >原始 ρ（尾部耦合增强）
• 电压 PDF：平滑、无负值（物理约束满足）

参考代码 电力系统概率潮流计算中的nataf变换 www.youwenfan.com/contentcnl/80243.html

六、常见扩展

1. 其他分布 → 在 edge_dist_sample.m 加 Gamma、Lognormal；
2. AC 潮流 → 把 dc_flow 换成 NR 或 FDLF；
3. 大规模 → 用 PC Expansion 替代 Monte-Carlo（留接口）；
4. GPU 加速 → 把 mvnrnd 换成 gpuArray 版本。