节点1为平衡节点，节点4为PV节点，节点2、3为PQ节点,写出这个两机四节点系统的matlab潮流计算程序，要求分别用N-R法和PQ分解法实现

下面分别给出两种不同的方法，一种是基于牛顿-拉夫逊迭代法（N-R法），另一种是基于PQ分解法。

N-R法：

% 定义电路参数
V_b = 1; % 电压基准值
S_b = 1; % 功率基准值
Z_b = V_b^2/S_b; % 阻抗基准值

% 定义节点导纳矩阵
Y = [3-3i, -1+i, 0, -2+i;
     -1+i, 2-2i, -1+i, 0;
     0, -1+i, 2-2i, -1+i;
     -2+i, 0, -1+i, 3-3i];
 
% 定义节点类型
% 1: 平衡节点, 2: PQ节点, 3: PV节点
type = [1; 2; 2; 3];

% 定义初始节点电压值
V = [1; 1; 1; 1];

% 定义PV节点有功功率和电压值
Pv = 0.3;
Vp = 1.05;

% 定义迭代参数
max_iter = 100; % 最大迭代次数
tol = 1e-6; % 收敛容限

% 开始迭代
for iter = 1:max_iter
    % 构造节点注入功率矩阵
    P_inj = zeros(4, 1);
    Q_inj = zeros(4, 1);
    for i = 1:4
        for j = 1:4
            P_inj(i) = P_inj(i) + V(i)*V(j)*real(Y(i,j));
            Q_inj(i) = Q_inj(i) + V(i)*V(j)*imag(Y(i,j));
        end
    end
    
    % 构造雅可比矩阵
    J = zeros(4, 4);
    for i = 1:4
        for j = 1:4
            if i == j
                J(i,j) = -imag(Y(i,j))*V(i) + Q_inj(i);
            else
                J(i,j) = V(i)*V(j)*imag(Y(i,j))*sin(angle(V(i))-angle(V(j))) - V(i)*V(j)*real(Y(i,j))*cos(angle(V(i))-angle(V(j)));
            end
        end
    end
    
    % 求解节点电压增量
    delta_V = J \ (-[P_inj(2); Q_inj(2); Q_inj(3)]);
    
    % 更新节点电压值
    V(2) = V(2) + delta_V(1);
    V(3) = V(3) + delta_V(2);
    V(4) = Vp;
    
    % 判断是否收敛
    if max(abs(delta_V)) < tol
        break;
    end
end

% 计算各节点功率及其它参数
P = zeros(4, 1);
Q = zeros(4, 1);
for i = 1:4
    for j = 1:4
        P(i) = P(i) + V(i)*V(j)*real(Y(i,j))*cos(angle(V(i))-angle(V(j)));
        Q(i) = Q(i) + V(i)*V(j)*imag(Y(i,j))*sin(angle(V(i))-angle(V(j)));
    end
end
V_ang = angle(V);
V_mag = abs(V);

% 输出结果
fprintf('迭代次数：%d\n', iter);
fprintf('节点电压：\n');
disp(V);
fprintf('节点电压幅值：\n');
disp(V_mag);
fprintf('节点电压相角：\n');
disp(V_ang);
fprintf('节点有功功率：\n');
disp(P);
fprintf('节点无功功率：\n');
disp(Q);

PQ分解法：

% 定义电路参数
V_b = 1; % 电压基准值
S_b = 1; % 功率基准值
Z_b = V_b^2/S_b; % 阻抗基准值

% 定义节点导纳矩阵
Y = [3-3i, -1+i, 0, -2+i;
     -1+i, 2-2i, -1+i, 0;
     0, -1+i, 2-2i, -1+i;
     -2+i, 0, -1+i, 3-3i];
 
% 定义节点类型
% 1: 平衡节点, 2: PQ节点, 3: PV节点
type = [1; 2; 2; 3];

% 定义初始节点电压值
V = [1; 1; 1; 1];

% 定义PV节点有功功率和电压值
Pv = 0.3;
Vp = 1.05;

% 分解节点导纳矩阵
G = real(Y);
B = imag(Y);
B_pq = B(2:3, 2:3);
B_pv = B([2,4], [2,4]);
G_pv = G([2,4], [2,4]);

% 计算各节点注入功率
S = V .* conj(Y * V);
P_inj = real(S);
Q_inj = imag(S);

% 构造PQ节点注入功率矩阵
PQ_inj = [P_inj(2:3); Q_inj(2:3)];

% 求解PQ节点电压幅值
V_pq = (B_pq * B_pq' - G_pv * inv(B_pv) * G_pv') \ (B_pq * PQ_inj - G_pv * inv(B_pv) * [Pv-Vp*real(Y(2,4)); 0]);

% 更新节点电压值
V(2:3) = V_pq;
V(4) = Vp;

% 计算各节点功率及其它参数
S = V .* conj(Y * V);
P = real(S);
Q = imag(S);
V_ang = angle(V);
V_mag = abs(V);

% 输出结果
fprintf('节点电压：\n');
disp(V);
fprintf('节点电压幅值：\n');
disp(V_mag);
fprintf('节点电压相角：\n');
disp(V_ang);
fprintf('节点有功功率：\n');
disp(P);
fprintf('节点无功功率：\n');
disp(Q);