loadflow.m

电力系统潮流计算

% YBUS Matrix

YS01 = 1/(0.02 + 0.06i);
YS02 = 1/(0.04 + 0.08i);
YS03 = 1/(0.01 + 0.03i);
YS13 = 1/(0.02 + 0.06i);
YS12 = 1/(0.03 + 0.09i);

YC01 = 0.03i;
YC02 = 0.04i;
YC03 = 0.05i;
YC13 = 0.06i;
YC12 = 0.07i;

YBUS = zeros(4, 4);
YBUS = [YS01 + YS02 + YS03 + YC01 + YC02 + YC03, -YS01, -YS02, -YS03; -YS01, YS01 + YS12 + YS13 + YC01 + YC12 + YC13, -YS12, -YS13 ; -YS02, -YS12, YS02 + YS12 + YC02 + YC12, 0; -YS03, -YS13, 0, YS03 + YS13 + YC03 + YC13];
%YBUS
% 函荐沥狼
syms v0;
syms v1;
syms v2;
syms v3;

syms th_0;
syms th_1;
syms th_2;
syms th_3;

%傈仿规沥侥
P1 = (v1 * v1 * (real(YBUS(2,2)) * cos(th_1 - th_1) + imag(YBUS(2,2)) * sin(th_1 - th_1))) + (v1 * v0 * (real(YBUS(2,1)) * cos(th_1 - th_0) + imag(YBUS(2,1)) * sin(th_1 - th_0))) + (v1 * v2 * (real(YBUS(2,3)) * cos(th_1 - th_2) + imag(YBUS(2,3)) * sin(th_1 - th_2))) + (v1 * v3 * (real(YBUS(2,4)) * cos(th_1 - th_3) + imag(YBUS(2,4)) * sin(th_1 - th_3)));
P2 = (v2 * v2 * (real(YBUS(3,3)) * cos(th_2 - th_2) + imag(YBUS(3,3)) * sin(th_2 - th_2))) + (v2 * v0 * (real(YBUS(3,1)) * cos(th_2 - th_0) + imag(YBUS(3,1)) * sin(th_2 - th_0))) + (v2 * v1 * (real(YBUS(3,2)) * cos(th_2 - th_1) + imag(YBUS(3,2)) * sin(th_2 - th_1)));
P3 = (v3 * v3 * (real(YBUS(4,4)) * cos(th_3 - th_3) + imag(YBUS(4,4)) * sin(th_3 - th_3))) + (v3 * v0 * (real(YBUS(4,1)) * cos(th_3 - th_0) + imag(YBUS(4,1)) * sin(th_3 - th_0))) + (v3 * v1 * (real(YBUS(4,2)) * cos(th_3 - th_1) + imag(YBUS(4,2)) * sin(th_3 - th_1)));
Q2 = (v2 * v2 * (real(YBUS(3,3)) * sin(th_2 - th_2) + imag(YBUS(3,3)) * cos(th_2 - th_2))) + (v2 * v0 * (real(YBUS(3,1)) * sin(th_2 - th_0) + imag(YBUS(3,1)) * cos(th_2 - th_0))) + (v2 * v1 * (real(YBUS(3,2)) * sin(th_2 - th_1) + imag(YBUS(3,2)) * cos(th_2 - th_1)));
Q3 = (v3 * v3 * (real(YBUS(4,4)) * sin(th_3 - th_3) + imag(YBUS(4,4)) * cos(th_3 - th_3))) + (v3 * v0 * (real(YBUS(4,1)) * sin(th_3 - th_0) + imag(YBUS(4,1)) * cos(th_3 - th_0))) + (v3 * v1 * (real(YBUS(4,2)) * sin(th_3 - th_1) + imag(YBUS(4,2)) * cos(th_3 - th_1)));

%imag(YS01)

% Jacobian matrix
J = zeros(5,5);

%h = zeros(5,1);
%h = [diff(Q3, th_1); diff(Q3, th_2); diff(Q3, th_3); diff(Q3, v2); diff(Q3, v3)]

J = [diff(P1, th_1), diff(P1, th_2), diff(P1, th_3), diff(P1, v2), diff(P1, v3); diff(P2, th_1), diff(P2, th_2), diff(P2, th_3), diff(P2, v2), diff(P2, v3); diff(P3, th_1), diff(P3, th_2), diff(P3, th_3), diff(P3, v2), diff(P3, v3); diff(Q2, th_1), diff(Q2, th_2), diff(Q2, th_3), diff(Q2, v2), diff(Q2, v3); diff(Q3, th_1), diff(Q3, th_2), diff(Q3, th_3), diff(Q3, v2), diff(Q3, v3);];

%扁瘤樊
v0 = 1;
v1 = 0.98;
th_0 = 0 / 180 *  3.14;
P1t = -0.2;
P2t = 0.2;
P3t = 0.3;
Q2t = 0.1;
Q3t = 0.4;


%檬扁蔼

th_1 = 3 / 180 *  3.14;
th_2 = 6 / 180 *  3.14;
th_3 = 9 / 180 *  3.14;
v2 = 0.9;
v3 = 1.1;

dx = zeros(5, 1);

%subs(v0 * v1 * (real(YS01) * cos(th_1 - th_0) + imag(YS01) * sin(th_1 - th_0)))
%subs(P1)

for i=1: 3
    dp = [P1t - subs(P1); P2t - subs(P2); P3t - subs(P3); Q2t - subs(Q2); Q3t - subs(Q3)];
    Jt = subs(J);
    Ji = Jt^-1;
    dx = Ji * dp;
    
    dp
    subs(Q3);
    Ji
    dx
    
    th_1 = th_1 + dx(1,1)
    th_2 = th_2 + dx(2,1)
    th_3 = th_3 + dx(3,1)
    v2 = v2 + dx(4,1)
    v3 = v3 + dx(5,1)
end



%dp = zeros(5, 1);
%dp(1,1) = P1t - subs(P1);
%dp(2,1) = P2t - subs(P2);
%dp(3,1) = P3t - subs(P3);
%dp(4,1) = Q2t - subs(Q2);
%dp(5,1) = Q3t - subs(Q3);
%Jt = subs(J^-1)
%dx = -1 * Jt * dp;


%subs(J)