MATLAB牛顿拉夫逊法算潮流分析

1、MATLA井顿拉夫逊法算潮流分析The function of this program is to use Newton, the rafson methodThe % B1 matrix: 1, the first branch of the branch; No. 2, terminal;Branch resistance; 4, the circuit (or transformer admittance);% 5, the change of the branch; The branch head is at the K side of 1, 1 side 0;The % 7, li

2、ne/transformer id (0/1) transformer parameter is calculated to the end of the branch when the end of the branch is on the K side, 0 to the first end% B2 matrix: 1, power of the node generator; The load power of the node;% 3, PQ node voltage initial value, or a given value of the balance node and PV

3、node voltageThe connection of the capacitor (inductance) of the shunt capacitance (inductance) is received by the node% 5, node classification label: 1 is the balance node (should be number 1); 2 for the PQ node; 3 for PV nodes;The clear;Isb = 1; % input ( please enter the balance bus node number: i

4、sb = );Pr = 1 e - 5; % input ( please input the error accuracy: pr = );N = 10; % input ( please enter the number of nodes: n = );Nl = 10; % input ( please enter a number of directions: nl = );B1 = 1, 2, 4, 4, 1, 4, 1, 1, 0, 0, 1, 1, 1, 1, 4, 1, 1, 1, 0, 1, 0,1, 04, 5.1 + 19.2 I 2.1e-4 I 1 0 0;3, 4.2

5、5 + 16i 1.75e-4 I 1, 0, 0;5, 4.25 + 16i 7e minus 4i 1 0;3, 5.1 + 19.2 I 2.1e-4 I 1 0 0;6, 4, 5.95 + 22.4 I 9.8e-4 I 1, 0, 0;7, 1.78 + 53.89 I 0, 38.5/231, 0, 1;8, 1.49 + 48.02 I 0, 11/231, 0, 1;9, 9, 9, 9, 4, 9, 1, 9, 1, 0, 11, 0, 1.5 10 2.46 + 70.17 I 0 38.5/231 0 1 % input ( please enter the matri

6、x formed by the branch parameter: B1 = );B2 is equal to 0, 0, 25.5, 0, 1;220 2 0 0 0;220 2 0 0 0;220 2 0 0 0;220 2 0 0 0;120, 231 3;0, 61.11 + 37.87 I 35, 0, 2;0, 47.53 + 29.46 I 10, 0, 2;0, 54.32 + 33.66 I 10, 0, 2;0 40.74 + 25.25 I 35 0 2) % input ( please enter the matrix of the parameter paramet

7、ers of each node: B2 = );% % n = 4; % input ( please enter the number of nodes: n = );% nl = 4; % input ( please enter a number of directions: nl= );% B1 is equal to 1, 2, 4 + 16i 0, 0, 0;% 1, 3, 4 + 16i 0, 0, 0;% 2, 3, 2 + 8i 0, 0, 0;% 2, 41.49 + 48.02 I 0, 11/110, 1 % input ( please enter the matr

8、ix formed by the branch parameter: B1 = );% B2 is equal to 0, 0, 115, 0, 1;% 0, 0, 110, 0, 2;% 0, 20 + 4i 110, 0, 2;% 0, 10 + 6i 100 2 % input ( please input the matrix of the parameters of each node: B2 = );% Y = zeros (n); E = zeros (1, n); F = zeros (1, n); V = zeros (1, n);Sida = zeros (1, n); S

9、1 = zeros (nl);The admittance matrix For I = 1: nl % from 1 to n1If B1 (I, 7) = 1 % If B1 (I, 6) = = 0 % left (the first end) is on one sideP = B1 (I, 1); Q = B1 (I, 2);The else % left node (the first end) is on the K sideP = B1 (I, 2); Q = B1 (I, 1);The endY (p, q) = Y (p, q) - 1. / (B1 (I, 3) * B1

10、 (I, 5). % the diagonal elementY of q, p is equal to Y of p, q. % the diagonal elementY (q, q) = Y (q, q) + 1)/(B1 (I, 3) * B1 (I, 5) A 2); % diagonal element K sideY (p) = Y (p) = Y (p) + 1. / B1 (I, 3) + B1 (I, 4); % diagonal element 1 + excitation admittanceElse P = B1 (I, 1); Q = B1 (I, 2);Y (p,

11、 q) = Y (p, q) - 1. / B1 (I, 3); % the diagonal elementY of q, p is equal to Y of p, q. % the diagonal elementY (q) = Y (q) = Y (q) + 1. / B1 (I, 3) + B1 (I, 4). / 2.0; % diagonally equal to half of the lineY (p) = Y (p) = Y (p) + 1. / B1 (I, 3) + B1 (I, 4). / 2.0; % diagonally, half of the power of

12、 the lineThe endThe endDisp ( the admittance matrix Y = ); Disp (Y);% given the initial node voltage andgiven each node power injection G = real (Y); B = imag (Y); The % decomposes the real and imaginary parts of the admittance matrixFor I = 1: n % given the real and imaginary parts of the initial v

13、oltage of each nodeE (I) = real (B2 (I, 3);F (I) = imag (B2 (I, 3);V (I) = abs (B2 (I, 3); % PV, balance node, and PQ node voltage modulusThe endFor I = 1: n % given each node injection powerS (I) = B2 (I, 1) -b2 (I, 2); The % I node is injected with the power SG - SLB (I, I) = B (I, I) + B2 (I, 4);

14、 The % I node has no amount of work compensation.The end% = = = = = = = = = = = = = = = = = = with Newton iteration - Ralph monson method is used to solve the nonlinear algebraic equation (power equation) = = = = = = = = = = = = = = = = = = = = = = =P = real (S); Q = imag (S); The % decomposition of

15、 the active and reactive power of each nodeICT1 = 0; IT2 = 1; N0 = 2 * n; N1 = N0 + 1; A = 0; % iterations,ICT1, a; The number of nodes that are not satisfied with the convergenceis IT2While IT2 = 0 % N0 = 2 * n jacobian; N1 + 1 = N0 extensionIT2 = 0; A = a + 1;JZ = Jacobi matrix (, num2str (a), ) c

16、ancellation operations;JZ1 = Jacobi matrix (, num2str (a), ) the second iteration ; JZ0= power equation number (, num2str (a), ) the time difference: % to calculate thepower of each node and power deviation and voltage deviation of PV nodesFor I = 1: n % n nodes 2n rows (two equations P and Q or U f

17、or each node)P = 2 * I - 1; M = p + 1; C (I) = 0; D (I) = 0;For j1 = 1: n % I row n columns (n nodes interadmittance and node voltage multiplied by the current)C (I) = C (I) + G (I) * e (j1) -B (I, j1) * f (j1); %( ej- Bij Gij * * fj) TOC o 1-5 h z D (I) = D (I) + G (I, j1) * f (j1) + B (I, j1) * e

18、(j1); %2(fj + Bij Gij * * ej)The endThe calculated value of the work and the reactive power P , Q P1 = C (I) * e (I) + f (I) * D (I); Power P % node calculation ei2(ej - Bij Gij * * fj) + fi2 (fj + Bij Gij * * ej)Q1 = C (I) * f (I) - e (I) * D (I); Power Q % node calculation fi2(ej - Bij Gij * * fj)

19、 -ei 2 (fj + Bij Gij * * ej)V2 = e (I) A 2 + f (I) A 2; % voltage die squaredThe difference between the power difference and the PV node is equal to = =If I = isb % non-equilibrium node (PQ or PV node)If B2 (I, 5) = 3 % non-pv nodes (only PQ nodes)J (m, N1) = P (I) -p1; The % PQ node has the power d

20、ifference J (m,N1) to extend column delta PJ (p, N1) = Q (I) -q1; The % PQ node has no work power difference J(p, N1) expands column delta QElse % PV nodes = = = = = =J (m, N1) = P (I) -p1; The % PV node has the power difference J (m,N1) to extend column delta PJ (p, N1) = V (I) A 2 - V2. The % PV n

21、ode voltage module (p, N1) extends column delta UThe endEnd % (if I = isb) non-equilibrium nodes (PQ or PV nodes)End % (for I = 1: n) n nodes 2n rows (two equations P and Q or U for each node)For m = 1: N0JJN1 (m) = J (m, N1);The endDisp (JZ0); Disp (JJN1);% deviation value judgmentpower and voltage

22、 deviation value whether meet the requirements of PV nodes For k = 3: N0 % removes all nodes other than the balanced node 1 and 2DET = abs (k, N1);If DET = pr % PQ nodes power deviation and the voltage deviation of the PV node are metIT2 = IT2 + 1; The number of nodes that does not meet the requirem

23、ent is plus 1The endThe endICT2 (a) = IT2; The number of nodes that do not meet therequirement; A is the number of iterationsICT1 = ICT1 + 1; % the number of iterationsIf ICT2 (a) = = 0; The number of nodes currently not satisfied iszeroBreak % exits the iteration operationThe end% aboveto calculate

24、 the various nodes of power and power deviation and voltagedeviation of PV nodes The Jacobi matrix formed the correct equation for the Jacobi matrixFor I = 2: n % n nodes 2n rows (two equations P and Q or U for each node)If I = isb % non-equilibrium node (PQ or PV node)If B2 (I, 5) = 3%, the element

25、 of the Jacobi matrix for the PQnode is equal to = = =C (I) = 0; D (I) = 0;For j1 = 1: n % I row n columns (n nodes interadmittance and nodevoltage multiplied by the current)C (I) = C (I) + G (I) * e (j1)-B (I, j1) * f (j1); %2 (ej- Bij Gij * * fj)D (I) = D (I) + G (I, j1) * f (j1) + B (I, j1) * e (

26、j1); %2(fj + Bij Gij * * ej)The endFor j1 = 2: n % I row n columns (2n Jacobi matrix elements dP/DE anddP/df or dQ/DE and dQ/df)If j1 = isb&j 1 = I % unbalanced node &non-diagonal elementX1 = -g (I, j1) * e (I) - B (I, j1) * f (I); % X1 = dP/DE = -dq/df = -x4X2 = B (I, j1) * e (I) - G (I, j1) * f (I

27、); X2 = dP/df = dQ/DE = X3The X3 = X2; % X2 = dp/df X3 = dQ/DEX4 = - X1; % X1 = dP/DE X4 = dQ/dfP = 2 * I - 1; Q = 2 * j1-1;J (p, q) = X3; M = p + 1; So its going to be dQ/DE J (p, N) = dQ.J (m, q) = X1; Q = q + 1; The dP/DE J (m, N) = dP.J (p, q) = X4; J (m, q) = X2; % X4 = dQ/df X2 = dp/dfElseif j

28、1 = = = i&j 1 = isb % non-equilibrium node & diagonal elementX1 = -c (I) -g (I) * e (I) - B (I) * f (I); % dP/DEX2 = -d (I) + B (I) * e (I) -g (I) * f (I); % dP/dfX3 = D (I) + B (I) * e (I) -g (I) * f (I); % dQ/DEX4 = -c (I) + G (I) * e (I) + B (I) * f (I); % dQ/dfP = 2 * I - 1; Q = 2 * j1-1; J (p,

29、q) = X3; % expansion column deltaQ J of p,N) = DQ;M = p + 1;J (m, q) = X1; Q = q + 1; J (p, q) = X4; % expansion column delta P J (m, N) = DP;J (m, q) = X2;The endElse % if B2 (I, 5) = 3 otherwise (that is the PV node)The element of the Jacobi matrix is equal to = = = = =For j1 = 1: nIf j1 = isb&j 1

30、 = I % unbalanced node &non-diagonal elementX1 = -g (I, j1) * e (I) - B (I, j1) * f (I); % dP/DEX2 = B (I, j1) * e (I) - G (I, j1) * f (I); % dP/dfX5 = 0; X6 = 0;P = 2 * I - 1; Q = 2 * j1-1; J (p, q) = x 5; % PV node voltage error.M = p + 1;J (m, q) = X1; Q = q + 1; J (p, q) = X6; The % PV node has

31、the work error J (m, N) = DP;J (m, q) = X2;Elseif j1 = = = i&j 1 = isb % non-equilibrium node & diagonal elementX1 = -c (I) -g (I) * e (I) - B (I) * f (I); % dP/DEX2 = -d (I) + B (I) * e (I) -g (I) * f (I); % dP/dfX5 = -2 * e (I);X6 is equal to minus 2 times f of I.P = 2 * I - 1; Q = 2 * j1-1; J (p,

32、 q) = x 5; % PV node voltage error.M = p + 1;J (m, q) = X1; Q = q + 1; J (p, q) = X6; The % PV node has the work error J (m, N) = DP;J (m, q) = X2;The endEnd % (if B2 (I, 5) = 3 else)End % (if I = isb)End % (for I = 1: n) n nodes 2n rows (two equations P and Q or U for each node)JZ0 = formed the fir

33、st (, num2str (a), ) the subjacobi matrix:;% disp (JZ0); Disp (J);% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = above to form a complete Jacobi matrix = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =The elimination of the modified equation formed by the Jacobi matrix

34、 is solved by the gaussian elimination method (by the column elimination and the return of the matrix)For k = 3: N0 % N0 = 2 * n (from the third row, the first, secondline is the equilibrium node)For k1 = k + 1: N1 % from the Jacobi element of k + 1 to the extension delta P, delta Q or delta UJ (k,

35、k1) = J (k, k1). / J (k, k); The non-diagonal elements of Krow K columns are normalized by the K line K column diagonal elementThe endJ (k, k) = 1; The % diagonally normalized K lines K column diagonal elements assignment 1in columns eliminationoperation =For k2 = k + 1: N0 % from k + 1 to 2 * n las

36、t rowFor k3 = k + 1: N1 % from k2 + 1 to the extension column cancels outthe triangle elements after the k + 1 rowJ (k2, k3) = J (k2, k3) -j (k2, k) * J (k, k3); % eliminationoperationEnd % USES the current row K3 column element minus the current row K column element times the K row K3 column elemen

37、tJ (k2, k) = 0; The k column element of the current row is 0The endThe end% JZ = Jacobi matrix (, num2str (a), ) cancellation operations; JZ1 = Jacobi matrix (, num2str (a), ) the second iteration ;% disp (JZ); Disp (J);% = = = = = = = = = = = = = = = = = = = = on column generation algorithm =For k

38、= N0: - 1:3For k1 = k - 1: -1:3J (k1, N1) = J (k1, N1) - J (k1, k) * J (k, N1);J (k1, k) = 0;The endThe endFor m = 1: N0JJN1 (m) = J (m, N1);The endDisp (JZ1); Disp (JJN1); % disp (J);% modify the node voltage For k = 3:2: N0-1L is equal to k plus 1.E (L) = e (L) -j (k, N1); % modification node volt

39、age real partK1 = k + 1;F (L) = f (L) -j (k1, N1); % modifies the voltage virtual part of the nodeU (L) = SQRT (e (L) A 2 + f (L) A 2);The endDisp ( each node voltage module ); Disp (U);% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = end of an iteration = = = = = = = = = = = = = = = =

40、= = = = = = = = = = = = = =The end*end of the iterative calculation of the output process * * * * * * * * *Disp ( iteration: );Disp (ICT1-1);Disp ( not the number of precision requirements: );Disp (ICT2);For k = 1: nV (k) = SQRT (e (k) A 2 + f (k) A 2); Percent calculates the magnitude of each nodes

41、 voltageSida (k) = atan (f (k). / e (k). The percent calculates the Angleof each nodes voltageE (k) = E (k) + f (k) * j; % will represent each node voltage in the pluralThe end% = = = = = = = = = = = = = = = calculated output = = = = = = = = =Disp ( the voltage complex value of the nodes is (the nod

42、e number is from small to large) : );Disp (E); The actual voltage value of each node is shown in the pluralDisp ( );Disp ( the value of the voltage modulus of each node is the size of the node number of nodes from small to large: );Disp (V); The % shows the magnitude of V for each nodeDisp ( );Disp

43、( the voltage phase of each node is sida (the node number is from small to large) : );Disp (sida); The % display the voltage phase Angle of each nodeFor p = 1: nC (p) = 0;For q = 1: nC (p) = C (p) + conj (p, q). % calculate the conjugate value of the injection current of each nodeThe endS (p) = E (p

44、) * C (p); % calculates the conjugate value of the powerS = voltage X injection current of each nodeThe endDisp ( power S of each node (node number from small to large) : );Disp (S); The % shows the injection power of each nodeDisp ( );Disp ( the first power Si on each branch of the branch is in the

45、 same order as you enter B1: );For I = 1: nlP = B1 (I, 1); Q = B1 (I, 2);If B1 (I, 7) = = 0Si (p, q) = E (p) * conj (E (p) * B1 (I, 4). / 2 + (E (p) -e(q). / B1 (I, 3).Siz (I) = Si (p, q);The elseIf B1 (I, 6) = = 0Si (p, q) = E (p) * (E (p) * B1 (I, 4).(p) * B1 (I, 5) -e (q) * (1. / (B1 (I, 3) * B1

46、(I, 5)Siz (I) = Si (p, q);The elseSi (p, q) = E (p) * conj (E (p) - E (q) * B1 (I, 5) * (1)/(B1 (I,3) * B1 (I, 5)入 2);Siz (I) = Si (p, q);The endThe endS ( num2str (p), , num2str (q),) = , num2str (Si (p, q)Disp (ZF);Disp ( );The endDisp ( the end power Sj of each branch is the same as when you ente

47、r B1) : );For I = 1: nlP = B1 (I, 1); Q = B1 (I, 2);If B1 (I, 7) = = 0(q, p) = E (q) * conj (E (q) * B1 (I, 4). / 2 + (E (q) -e (p). / B1(I, 3).Sjy (I) = Sj (q, p);The elseIf B1 (I, 6) = = 0Sj (q) (q, p) = E * conj (E (q) - E (p) * B1 (I, 5) * (1)/(B1 (I,3) * B1 (I, 5)入 2);Sjy (I) = Sj (q, p);The el

48、se(q, p) = E (q) * (E (q) * B1 (I, 4).+ (E (q) * B1 (I, 5) -e (p) * (1. / (B1 (I, 3) * B1 (I, 5)Sjy (I) = Sj (q, p);The endThe endS ( S ( num2str ), , num2str (p),) = , Disp (ZF);Disp ();The endDisp (the power loss DS for each branch of the branch is the same as when you enter B1 in the order);For I

49、 = 1: nlP = B1 (I, 1); Q = B1 (I, 2);DS (I) = Si (p) + Sj (q, p);The ZF = DS (, num2str (p), , , num2str (q), ) = , num2str(DS (I).Disp (ZF);Disp ( );The endZWS = 0; JDDY = ; JDP = ; JDQ = ; JDDYJD = ;For I = 1: n % total net loss for all nodes injected with power algebra andZWS = ZWS + S (I);JDDYJD = strcat (JDDYJD, num2str (I), (, num2str ( I), )JDDY = strcat (JDDY, num2str (I), (, num2str (V (I),);JDP = strcat (JDP, num2str (I), (, num2str),);(JDQ, num2str (I)