用matlab写一个基于加权最小二乘法的电压状态估计并生成图的代码

以下是一个基于加权最小二乘法的电压状态估计的 MATLAB 代码示例，其中使用 IEEE 14 节点系统模型进行演示。

% 读取 IEEE 14 节点系统的数据
busdata = [
     1   1.0600   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     2   1.0450   0.0000  21.7000  12.7000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     3   1.0100   0.0000  94.2000  19.0000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     4   1.0000   0.0000  47.8000 -3.9000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     5   1.0000   0.0000   7.6000   1.6000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     6   1.0000   0.0000  11.2000   7.5000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     7   1.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     8   1.0000   0.0000  71.4000  17.4000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
     9   1.0000   0.0000  17.4000   1.8000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
    10   1.0000   0.0000  20.0000   5.8000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
    11   1.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   1.0000
    12   1.0000   0.0000   0.0000   0.0000   0.0000  12.0000   0.0000  10.6000   0.0000   1.0000
    13   1.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   8.2000   0.0000   1.0000
    14   1.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000  22.8000   0.0000   1.0000
];
linedata = [
     1     2   0.0192   0.0575   0.0264
     1     5   0.0540   0.2230   0.0492
     2     3   0.0460   0.1160   0.0204
     2     4   0.0580   0.1760   0.0320
     2     5   0.0560   0.1490   0.0264
     3     4   0.0670   0.1710   0.0296
     4     5   0.0132   0.0420   0.0072
     4     7   0.0000   0.2090   0.0000
     4     9   0.0000   0.5560   0.0000
     5     6   0.0000   0.2520   0.0000
     6    11   0.0940   0.1980   0.0338
     6    12   0.0000   0.2550   0.0000
     6    13   0.0000   0.1400   0.0000
     7     8   0.0000   0.1760   0.0000
     7     9   0.0000   0.1100   0.0000
     9    10   0.0310   0.0840   0.0140
    9     14   0.0000   0.1200   0.0000
    10    11   0.0530   0.1830   0.0264
    12    13   0.0000   0.0840   0.0000
    13    14   0.0000   0.0580   0.0000
];

% 构建节点导纳矩阵 Ybus
nl = size(linedata, 1);
nb = size(busdata, 1);
Ybus = zeros(nb, nb);
for k=1:nl
    f = linedata(k, 1);
    t = linedata(k, 2);
    r = linedata(k, 3);
    x = linedata(k, 4);
    b = linedata(k, 5);
    z = r + 1j*x;
    y = 1/z;
    Ybus(f, f) = Ybus(f, f) + y + 1j*b/2;
    Ybus(t, t) = Ybus(t, t) + y + 1j*b/2;
    Ybus(f, t) = Ybus(f, t) - y;
    Ybus(t, f) = Ybus(t, f) - y;
end

% 定义测量数据
% 第一列为节点编号，第二列为电压幅值，第三列为相角，第四列为有功功率，第五列为无功功率
% 测量数据包括节点 2、3、6、8、10、11、12、13、14 共 9 个节点
z = [
    2   1.045   0.000   0.000   0.000
    3   1.010   0.000   0.000  94.200
    6   1.000   0.000   0.000  11.200
    8   1.000   0.000   0.000  71.400
   10   1.000   0.000   0.000  20.000
   11   1.000   0.000   0.000   0.000
   12   1.000  10.600  12.000   0.000
   13   1.000   8.200   0.000   0.000
   14   1.000  22.800   0.000   0.000
];

% 计算节点注入功率
P = real(z(:, 4)) / 100;
Q = imag(z(:, 4)) / 100;

% 初始化电压幅值和相角
V = busdata(:, 2);
theta = zeros(nb, 1);

% 构建 W 矩阵
sigma = 1; % 调整权重的参数
W = zeros(2*nb, 2*nb);
for k=1:size(z, 1)
    i = z(k, 1);
    Pm = P(k);
    Qm = Q(k);
    Vm = z(k, 2);
    theta_m = z(k, 3);
    G = real(Ybus(i, i));
    B = imag(Ybus(i, i));
    W(2*i-1, 2*i-1) = (sigma*Vm*G)^(-2);
    W(2*i, 2*i) = (sigma*Vm*B)^(-2);
    W(2*i-1, 2*i) = -sigma*Vm*Qm/(Vm^2*G^2 + Vm^2*B^2);
    W(2*i, 2*i-1) = -sigma*Vm*Pm/(Vm^2*G^2 + Vm^2*B^2);
end

% 迭代计算电压状态估计
maxiter = 50; % 最大迭代次数
tol = 1e-5; % 收敛容差
for iter=1:maxiter
    % 计算注入功率和注入电流
    Pinj = zeros(nb, 1);
    Qinj = zeros(nb, 1);
    I = zeros(nb, 1);
    for k=1:nb
        for m=1:nb
            Vkm = V(k) - V(m);
            theta_km = theta(k) - theta(m);
            Ykm = Ybus(k, m);
            Pinj(k) = Pinj(k) + V(k)*V(m)*(real(Ykm)*cos(theta_km) + imag(Ykm)*sin(theta_km));
            Qinj(k) = Qinj(k) + V(k)*V(m)*(real(Ykm)*sin(theta_km) - imag(Ykm)*cos(theta_km));
        end
        I(k) = Pinj(k) - 1j*Qinj(k);
    end

    % 构建雅可比矩阵
    J = zeros(2*nb, 2*nb);
    for k=1:nb
        for m=1:nb
            if k == m
                J(2*k-1, 2*m-1) = -Qinj(k) - V(k)^2*imag(Ybus(k, k));
                J(2*k-1, 2*m) = Pinj(k) - V(k)^2*real(Ybus(k, k));
                J(2*k, 2*m-1) = Pinj(k) + V(k)^2*real(Ybus(k, k));
                J(2*k, 2*m) = Qinj(k) - V(k)^2*imag(Ybus(k, k));
            else
                Vkm = V(k) - V(m);
                theta_km = theta(k) - theta(m);
                Ykm = Ybus(k, m);
                J(2*k-1, 2*m-1) = -V(k)*V(m)*(real(Ykm)*sin(theta_km) - imag(Ykm)*cos(theta_km));
                J(2*k-1, 2*m) = V(m)*(real(Ykm)*cos(theta_km) + imag(Ykm)*sin(theta_km));
                J(2*k, 2*m-1) = V(k)*V(m)*(real(Ykm)*cos(theta_km) + imag(Ykm)*sin(theta_km));
                J(2*k, 2*m) = V(k)*(real(Ykm)*sin(theta_km) - imag(Ykm)*cos(theta_km));
            end
        end
    end

    % 计算电压状态估计增量
    dX = inv(J'*W*J)*J'*W*[P; Q] - [V.*cos(theta); V.*sin(theta)];

    % 更新电压幅值和相角
    V = V + dX(1:nb);
    theta = theta + dX(nb+1:end);

    % 判断是否收敛
    if max(abs(dX)) < tol
        break;
    end
end

% 画出节点电压幅值和相角的估计结果
figure;
subplot(2, 1, 1);
bar(1:nb, V);
title('Estimated voltage magnitude');
xlabel('Bus number');
ylabel('Voltage magnitude (p.u.)');
subplot(2, 1, 2);
bar(1:nb, theta*180/pi);
title('Estimated voltage angle');
xlabel('Bus number');
ylabel('Voltage angle (degree)');

该代码将输出一个包含电压幅值和相角的图表，显示节点电压状态估计的结果。