33节点配电网潮流分区matlab程序

以下是一个简单的33节点配电网潮流分区的Matlab程序，仅供参考：

% 节点数据
busdata = [1 0 0 0.95 0 0 0;
           2 3 1 1.02 0 0 0;
           3 2 1 1.01 0 0 0;
           4 3 1 1.03 0 0 0;
           5 2 1 0.98 0 0 0;
           6 5 1 0.97 0 0 0;
           7 1 1 1.05 0 0 0;
           8 6 1 0.99 0 0 0;
           9 7 1 1.07 0 0 0;
           10 8 1 1.06 0 0 0;
           11 10 1 1.04 0 0 0;
           12 9 1 1.01 0 0 0;
           13 12 1 1.015 0 0 0;
           14 13 1 1.03 0 0 0;
           15 14 1 1.01 0 0 0;
           16 15 1 1.02 0 0 0;
           17 16 1 1.01 0 0 0;
           18 9 1 1.01 0 0 0;
           19 18 1 1.02 0 0 0;
           20 19 1 1.01 0 0 0;
           21 20 1 1.02 0 0 0;
           22 21 1 1.01 0 0 0;
           23 22 1 1.02 0 0 0;
           24 23 1 1.01 0 0 0;
           25 24 1 1.02 0 0 0;
           26 25 1 1.01 0 0 0;
           27 26 1 1.02 0 0 0;
           28 27 1 1.01 0 0 0;
           29 28 1 1.02 0 0 0;
           30 29 1 1.01 0 0 0;
           31 30 1 1.02 0 0 0;
           32 31 1 1.01 0 0 0;
           33 32 1 1.02 0 0 0];

% 支路数据
linedata = [1 2 0.02 0.06 0.03 0;
            2 3 0.045 0.135 0.02 0;
            3 4 0.02 0.06 0.03 0;
            2 5 0.024 0.072 0.01 0;
            5 6 0.011 0.033 0.01 0;
            6 4 0.015 0.045 0.02 0;
            1 7 0.018 0.036 0.01 0;
            7 8 0.02 0.04 0.02 0;
            8 9 0.022 0.044 0.02 0;
            9 10 0.02 0.04 0.02 0;
            10 11 0.02 0.06 0.03 0;
            9 12 0.02 0.06 0.03 0;
            12 13 0.022 0.066 0.01 0;
            13 14 0.02 0.06 0.03 0;
            14 15 0.02 0.06 0.03 0;
            15 16 0.02 0.06 0.03 0;
            16 17 0.02 0.06 0.03 0;
            12 18 0.02 0.06 0.03 0;
            18 19 0.02 0.06 0.03 0;
            19 20 0.02 0.06 0.03 0;
            20 21 0.02 0.06 0.03 0;
            21 22 0.02 0.06 0.03 0;
            22 23 0.02 0.06 0.03 0;
            23 24 0.02 0.06 0.03 0;
            24 25 0.02 0.06 0.03 0;
            25 26 0.02 0.06 0.03 0;
            26 27 0.02 0.06 0.03 0;
            27 28 0.02 0.06 0.03 0;
            28 29 0.02 0.06 0.03 0;
            29 30 0.02 0.06 0.03 0;
            30 31 0.02 0.06 0.03 0;
            31 32 0.02 0.06 0.03 0;
            32 33 0.02 0.06 0.03 0];

% 构造节点导纳矩阵
nbus = size(busdata,1);
Ybus = zeros(nbus);
for k=1:size(linedata,1)
    f = linedata(k,1);
    t = linedata(k,2);
    R = linedata(k,3);
    X = linedata(k,4);
    B = linedata(k,5);
    Y = 1/(R+j*X);
    Ybus(f,f) = Ybus(f,f) + Y + j*B/2;
    Ybus(t,t) = Ybus(t,t) + Y + j*B/2;
    Ybus(f,t) = Ybus(f,t) - Y;
    Ybus(t,f) = Ybus(t,f) - Y;
end

% 潮流计算
Pbusinj = busdata(:,6);
Qbusinj = busdata(:,7);
Vbus = busdata(:,4);
theta = zeros(nbus,1);
iter = 0;
converged = 0;
while ~converged && iter<20
    iter = iter + 1;
    Pbus = Vbus.*conj(Ybus*Vbus);
    Qbus = Vbus.*conj(Ybus*Vbus);
    mis = [Pbus-Pbusinj; Qbus-Qbusinj];
    J = zeros(2*nbus,2*nbus);
    for k=1:nbus
        for m=1:nbus
            if k==m
                dPdVkm = real(conj(Vbus(k))*(Ybus(k,:)*Vbus) + conj(Ybus(k,k))*Vbus(k));
                dPdthetakm = -imag(conj(Vbus(k))*(Ybus(k,:)*Vbus) + conj(Ybus(k,k))*Vbus(k));
                dQdVkm = -imag(conj(Vbus(k))*(Ybus(k,:)*Vbus) + conj(Ybus(k,k))*Vbus(k));
                dQdthetakm = -real(conj(Vbus(k))*(Ybus(k,:)*Vbus) + conj(Ybus(k,k))*Vbus(k));
            else
                dPdVkm = real(conj(Vbus(k))*Ybus(k,m)*Vbus(m));
                dPdthetakm = -imag(conj(Vbus(k))*Ybus(k,m)*Vbus(m));
                dQdVkm = -imag(conj(Vbus(k))*Ybus(k,m)*Vbus(m));
                dQdthetakm = -real(conj(Vbus(k))*Ybus(k,m)*Vbus(m));
            end
            J(k,nbus+m) = dPdVkm;
            J(nbus+k,nbus+m) = dQdVkm;
            J(k,m) = dPdthetakm;
            J(nbus+k,m) = dQdthetakm;
        end
    end
    dx = -J\mis;
    dtheta = dx(1:nbus);
    dV = dx(nbus+1:2*nbus);
    theta = theta + dtheta;
    Vbus = Vbus + dV;
    converged = max(abs(mis)) < 1e-7;
end

% 进行分区，定义分区1和分区2的节点
partition1 = [1 2 3 4 5 6 7 8 9 10 11 12 18];
partition2 = [13 14 15 16 17 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33];

% 输出分区1的潮流结果
fprintf('节点     电压幅值      电压相角      有功功率      无功功率\n');
for i=1:length(partition1)
    fprintf('%4d %12.4f %12.4f %12.4f %12.4f\n', partition1(i), abs(Vbus(partition1(i))), ...
            angle(Vbus(partition1(i)))*180/pi, real(Pbus(partition1(i))), imag(Qbus(partition1(i))));
end

% 输出分区2的潮流结果
fprintf('\n节点     电压幅值      电压相角      有功功率      无功功率\n');
for i=1:length(partition2)
    fprintf('%4d %12.4f %12.4f %12.4f %12.4f\n', partition2(i), abs(Vbus(partition2(i))), ...
            angle(Vbus(partition2(i)))*180/pi, real(Pbus(partition2(i))), imag(Qbus(partition2(i))));
end

程序中先定义了节点数据和支路数据，然后根据支路数据构造了节点导纳矩阵。接着进行潮流计算，直到收敛为止。最后按照给定的分区将节点分为两个部分，并输出每个节点的潮流结果。

注意，这里的潮流计算只是一个简单的示例，实际情况下可能需要更复杂的算法和数据处理。