潮流计算（直角坐标，极坐标以及快速解耦）matlab程序

%潮流计算主程序 function PowerFlow1(Y, Info, eps) if nargin == 2 eps = 1.0e-3; end disp('潮流计算-直角坐标 Copyright@2010 Southeast University'); Init = SetInit1(Info); DeltaPQ = UnBalance1(Init, Info, Y); tol = 1; time = 0; while tol disp('迭代次数'); time = time + 1 order = 0; n = length(Init(:, 1)); for i = 1:n DeltaUTheta(i, 1) = 0; DeltaUTheta(i, 2) = 0; if Init(i, 5) == 1 continue; elseif Init(i, 5) == 2 order = order + 1; PQ(order) = DeltaPQ(i, 1); order = order + 1; PQ(order) = DeltaPQ(i, 2); else order = order + 1; PQ(order) = DeltaPQ(i, 1); order = order + 1; PQ(order) = DeltaPQ(i, 2); end end disp('功率、电压不平衡量'); PQ' for i = 1:n if Init(i, 5) == 1 continue; elseif Init(i, 5) == 2 label_1 = translabel1(i); for j = 1:n if Init(j, 5) == 1 continue; elseif Init(j, 5) == 2 label_2 = translabel1(j); if label_1 == label_2 Jacobi(label_1, label_2) = PEii(i, Y, Init); Jacobi(label_1, label_2 + 1) = PFii(i, Y, Init); Jacobi(label_1 + 1, label_2 + 1) = UFii(i, Y, Init); Jacobi(label_1 + 1, label_2) = UEii(i, Y, Init); else Jacobi(label_1, label_2) = PEij(i, j, Y, Init); Jacobi(label_1, label_2 + 1) = PFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2 + 1) = UFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2) = UEij(i, j, Y, Init); end else label_2 = translabel1(j); Jacobi(label_1, label_2) = PEij(i, j, Y, Init); Jacobi(label_1, label_2 + 1) = PFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2 + 1) = UFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2) = UEij(i, j, Y, Init); end end else label_1 = translabel1(i); for j = 1:n if Init(j, 5) == 1 continue; elseif Init(j, 5) == 2 label_2 = translabel1(j); Jacobi(label_1, label_2) = PEij(i, j, Y, Init); Jacobi(label_1, label_2 + 1) = PFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2 + 1) = QFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2) = QEij(i, j, Y, Init); else label_2 = translabel1(j, Init); if label_1 == label_2 Jacobi(label_1, label_2) = PEii(i, Y, Init); Jacobi(label_1, label_2 + 1) = PFii(i, Y, Init); Jacobi(label_1 + 1, label_2 + 1) = QFii(i, Y, Init); Jacobi(label_1 + 1, label_2) = QEii(i, Y, Init); else Jacobi(label_1, label_2) = PEij(i, j, Y, Init); Jacobi(label_1, label_2 + 1) = PFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2 + 1) = QFij(i, j, Y, Init); Jacobi(label_1 + 1, label_2) = QEij(i, j, Y, Init); end end end end end disp('雅可比矩阵'); Jacobi UE = inv(-Jacobi) * (PQ'); disp('修正量'); UE sign = 0; tol = 0; n = length(UE); for i = 1:2:n if (UE(i) > eps)|(UE(i + 1) > eps) tol = 1; end label = restorelabel1(i); DeltaUE(label, 1) = UE(i); DeltaUE(label, 2) = UE(i + 1); end Init = Amend1(Init, DeltaUE); disp('修正电压实部和电压虚部'); Init(:,1:2) DeltaPQ = UnBalance1(Init, Info, Y); end