电力系统稳态实验报告

电力系统稳态潮流计算上机实验报告

一、问题

如下图所示的电力系统网络，分别用牛顿拉夫逊法、PQ解耦法、高斯赛德尔法、保留非线性法计算该电力系统的潮流。

发电机的参数如下，

*表示任意值 负荷参数如下，

如上图所示的电力系统，可以看出，节点1、2、3是PQ节点，节点4是PV节点，而将节点5作为平衡节点。根据问题所需，采用牛顿拉夫逊法、PQ解耦法、高斯赛德尔法、保留非线性法，通过对每次修正量的收敛判据的判断，得出整个电力系统的潮流，并分析这四种方法的收敛速度等等。 算法分析

1.牛顿拉夫逊法

节点5为平衡节点，不参加整个的迭代过程，节点1、2、3为PQ节点，节点4为PV节点，计算修正方程中各量，进而得到修正量，判断修正量是否收敛，如果不收敛，迭代继续，如果收敛，算出PQ节点的电压幅值以及电压相角，得出PV节点的无功量以及电压相角，得出平衡节点的输出功率。

潮流方程的直角坐标形式，

Pi ei Gijej Bijfj fi Gijfj Bijej

j i

j i

Qi fi Gijej Bijfj ei Gijfj Bijej

j i

j i

直角坐标形式的修正方程式，

P

Qn m 1

m2 U

n 1

H

M RN

e

L

f S

n 1n 1

修正方程式中的各量值的计算，

Pi pis [ei Gijej Bijfj fi Gijfj Bijej ]

j i

j i

Qi Qis [fi Gijej Bijfj ei Gijfj Bijej ]

j i

j i

222

Ui Uis (ei fi2)

Jacobi矩阵的元素计算，

Mij

Bijei Gijfi (j i) Qi

Gijfj Bijej Biiei Giifi(j i) ej

j i

Lij

Gijei Bijfi (j i) Qi

Gijej Bijfj Giiei Biifi(j i) fj

j i

Rij

U ej

2

i

0

2ei

2

(j i)(j i)

Sij

Ui fj

0

2fi

(j i)(j i)

牛顿拉夫逊法潮流计算的流程图如下，

2.PQ解耦法

如同牛顿拉夫逊法，快速解耦法的前提是，输电线路的阻抗要比电阻大得多，并且输电线路两端的电压相角相差不大，此时可利用PQ快速解耦法，来计算整个电力系统网络的潮流。

快速解耦法的迭代方程组，

P＝－H θ Q＝－L( U/U)

快速解耦法潮流计算的流程图如下，

图1-2 快速解耦法的程序原理框图

3.高斯赛德尔法

高斯赛德尔法原理较前两种方法简单，程序设计十分容易，占内存小，是所有的潮流计算方法中迭代计算量最小的。

高斯赛德尔法的迭代格式为，

Ui

.(k 1)

ss1 Pi jQi Y

(k)Yii

Ui

.(k 1) i 1

i1U1 YijUj

j 2

.s

j i 1

Y

n

ij

U

.(k)j

(1 17)

高斯赛德尔法的收敛判据如下，

maxi

.(k 1)

i

U

.(k)

i

在高斯赛德尔法中，不应对PV节点的幅值进行修正，只对其电压相角进行一定的修正。

高斯赛德尔法潮流计算的流程图为，

4.保留非线性法

保留非线性法主要是在牛顿拉夫逊法的基础上，通过泰勒展开，保留到二阶项，由于三阶导数值等于零，所以泰勒展开式是准确的，无截断误差，与牛顿拉夫逊法不同的是，保留非线性法只需算一次jacobi矩阵，每次迭代得到的修正量都是在初始值上的修正量，因此，保留非线性法的计算量小于牛顿拉夫逊法的计算量，大大节约计算机的内存空间，提高计算机的计算速度。

保留非线性的迭代格式为，

x(k 1) (J(x(0))) 1[y(x(0)) ys+y( x(k))]

x(k 1) x(0) x(k 1)

式中，k表示迭代次数；J为按x＝x(0)估计而得。

收敛判据为，

max xi(k 1) xi(k)

i

也可采用相继二次迭代的二阶项之差作为收敛判据(更合理)，相应的收敛判据如下，

maxyi( x(k 1)) yi( x(k))

i

保留非线性法的流程图如下，

三、MATLAB仿真结果

1.牛顿拉夫逊法

迭代次数k=6；

各节点的电压值、有功功率以及无功功率见下表。

2.PQ解耦法

迭代次数k=13；

迭代次数k=137;

4.保留非线性法 迭代次数k=11。

各节点的电压值、有功功率以及无功功率见下表。 从以上的MATLAB仿真结果可以看出，牛顿拉夫逊法只需迭代6次，迭代次数最少，保留非线性迭代11次，大约是牛拉法的两倍，PQ解耦法迭代次数13次，收敛速度相比于保留非线性稍慢，而高斯赛德尔法迭代次数达到137次，高斯赛德尔法算法简单，占用内存小，但是牺牲迭代次数。

从以上的仿真结果可以得出，牛拉法收敛速度快，算法具有平方收敛特性，是所有算法中收敛最快的，具有良好的收敛可靠性，并且牛顿法所需的内存量及每次迭代的时间均较高斯赛德尔法多。

在PQ解耦法中，用解两个阶数几乎减半的方程组(一个n-1及一个n-m-1)代替牛顿法的结一个2n-m-2阶方程组，显著地减少了内存需求量及计算量，系数矩阵B’及B’’是两个常数阵，为此只需在迭代循环前一次形成并进行三角分解组成因子表，在迭代过程中反复应用，大大缩短了每次迭代所需时间。快速解耦法达到收敛所需的迭代次数比牛顿法多，快速解耦法的程序设计较牛顿法简单，但从牛顿法到快速解耦法的演化时在元件的R<<X以及线路两端相角差比较小等假设基础上进行的，当系统不符合这些假设时，迭代就会出现问题。

高斯赛德尔法中，原理简单，程序设计十分容易，线性非线性方程组均适用，并且导纳矩阵是一个对称且高度稀疏的矩阵，因此占用内存非常节省，每次迭代的计算量也小，是各种潮流算法中最小的。但是收敛速度很慢，迭代次数将随所计算网络节点数的增加而直线上升，从上文的仿真结果就能看出，收敛速度是四种方法中最慢的。

保留非线性法中的雅可比矩阵，只需一次形成，并由三角分解构成因子表，而牛顿法中，每次重新形成因子表，保留非线性与牛拉法最大的区别在于 x(k)的含义，在保留非线性中， x(k)是相对于始终不变的初始估计值x(0)的修正量，而在牛拉法中， x(k)是相对于上一次迭代所得到的迭代点x(k)的修正量，但是保留非线性法达到收敛所需迭代次数多，收敛特性为直线但总计算速度较快。保留非线性法在收敛性方面，属于“等斜率法”的范畴，和牛顿法的平方收敛特性相比，达到收敛的迭代次数较牛顿法多，较快速解耦法，收敛的可靠性更好，计算速度可以接近快速解耦法。 五、证明

= +

证明：

因为 + = +

=

=

( )2= ( + ) = +

所以 + =Re + = = 附录： 1. 牛拉法

2. clear; 3. clc;

4. yb=zeros(5,5);

5. yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i; 6. yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i; 7. yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i; 8. yb(4,4)=(-66.66667i)/2; 9. yb(5,5)=(-33.33333i)/2; 10. yb=yb+conj(yb'); 11. k=0;

12. eps1=10^-4; 13. jeps=1; 14. G=real(yb); 15. B=imag(yb);

16. e=[1; 1 ;1 ;1.05 ;1.05]; 17. f=zeros(5,1);

18. pis=[-1.6; -2 ;-3.7; 5]; 19. qis=[-0.8; -1; -1.3]; 20. deta_p=zeros(4,1); 21. deta_q=zeros(3,1); 22. deta=zeros(8,1); 23. deta_ef=zeros(8,1); 24. deta_e=zeros(4,1); 25. deta_f=zeros(4,1); 26. U=zeros(5,1); 27. while (jeps>eps1) 28. p=zeros(4,1); 29. q=zeros(3,1); 30. for i=1:4 31. for j=1:5

32. p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j))+f(i)*(G(i,j)*f(j)+B(i,j)*e(j));

2

2

34. deta_p(i)=pis(i)-p(i); 35. end 36. for i=1:3 37. for j=1:5

38. q(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j))-e(i)*(G(i,j)*f(j)+B(i,j)*e(j)); 39. end

40. deta_q(i)=qis(i)-q(i); 41. end

42. deta_UU=1.05*1.05-(e(4)*e(4)+f(4)*f(4)); 43. jacobi=zeros(8,8); 44. for i=1:4 45. for j=1:4

46. jacobi(2*i-1,2*j-1)=-G(i,j)*e(i)-B(i,j)*f(i); 47. jacobi(2*i-1,2*j)=B(i,j)*e(i)-G(i,j)*f(i); 48. end 49. end 50. for i=1:3 51. for j=1:4

52. jacobi(2*i,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i); 53. jacobi(2*i,2*j)=G(i,j)*e(i)+B(i,j)*f(i); 54. end 55. end

56. for i=1:2:7

57. jacobi(i,i)=0; 58. jacobi(i,i+1)=0; 59. for j=1:5

60. jacobi(i,i)=jacobi(i,i)-(G((i+1)/2,j)*e(j)-B((i+1)/2,j)*f(j)); 61. jacobi(i,i+1)=jacobi(i,i+1)-(G((i+1)/2,j)*f(j)+B((i+1)/2,j)*e(j)); 62. end 63.

jacobi(i,i)=jacobi(i,i)-G((i+1)/2,(i+1)/2)*e((i+1)/2)-B((i+1)/2,(i+1)/2)*f((i+1)/2); 64.

jacobi(i,i+1)=jacobi(i,i+1)+B((i+1)/2,(i+1)/2)*e((i+1)/2)-G((i+1)/2,(i+1)/2)*f((i+1)/2); 65. end

66. for i=2:2:6

67. jacobi(i,i-1)=0; 68. jacobi(i,i)=0; 69. for j=1:5

70. jacobi(i,i-1)=jacobi(i,i-1)-(G(i/2,j)*f(j)-B(i/2,j)*e(j)); 71. jacobi(i,i)=jacobi(i,i)-(G(i/2,j)*e(j)-B(i/2,j)*f(j)); 72. end

73. jacobi(i,i-1)=jacobi(i,i-1)+B(i/2,i/2)*e(i/2)-G(i/2,i/2)*f(i/2); 74. jacobi(i,i)=jacobi(i,i)+G(i/2,i/2)*e(i/2)+B(i/2,i/2)*f(i/2);

76. jacobi(8,7)=-2*e(4); 77. jacobi(8,8)=-2*f(4); 78. for i=1:2:7

79. deta(i)=deta_p((i+1)/2); 80. end

81. for i=2:2:6

82. deta(i)=deta_q(i/2); 83. end

84. deta(8)=deta_UU;

85. deta_ef=-inv(jacobi)*deta; 86. for i=1:4

87. deta_e(i)=deta_ef(2*i-1); 88. deta_f(i)=deta_ef(2*i); 89. end 90. for i=1:4

91. e(i)=e(i)+deta_e(i); 92. f(i)=f(i)+deta_f(i); 93. end

94. jeps=max(max(abs(deta_ef)),max(abs(deta_p))); 95. jeps=max(jeps,max(abs(deta_q))); 96. k=k+1; 97. end 98. for i=1:5

99. U(i)=e(i)+sqrt(-1)*f(i); 100. end

101. disp('µçѹµÄ·ùÖµ£º'); 102. disp(abs(U)); 103. disp('µçѹ£º'); 104. disp(U);

105. theta=zeros(5,1); 106. for i=1:5

107. theta(i)=angle(U(i)); 108. end

109. disp(theta);

110. y0=[-4i;-8i;-4i;0;0]; 111. s=zeros(5,5); 112. for i=1:5 113. for j=1:5

114. s(i,j)=U(i)*(conj(U(i))*conj(y0(i))+(conj(U(i))-conj(U(j)))*conj(-yb(i,j))); 115. end 116. end 117. for i=1:5 118. s(i,i)=0;

120. deta_s=zeros(5,5); 121. p=zeros(5,1); 122. q=zeros(5,1); 123. for i=1:5 124. for j=1:5

125. p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j))+f(i)*(G(i,j)*f(j)+B(i,j)*e(j)); 126. end 127. end 128. for i=1:5 129. for j=1:5

130. q(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j))-e(i)*(G(i,j)*f(j)+B(i,j)*e(j)); 131. end 132. end 133. for i=1:5 134. for j=1:5

135. deta_s(i,j)=s(i,j)+s(j,i); 136. end 137. end 138. disp(s); 139. disp(deta_s); 140. disp(p);

disp(q);

2.PQ解耦法

clear; clc;

yb=zeros(5,5);

yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i; yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i; yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i; yb(4,4)=(-66.66667i)/2; yb(5,5)=(-33.33333i)/2; yb=yb+conj(yb'); k=1; kp=0; kq=0; eps1=10^-5; G=real(yb); B=imag(yb); B1=B(1:4,1:4); B2=B(1:3,1:3);

U=[1; 1 ;1 ;1.05 ;1.05]; theta=zeros(5,1); deta_theta=zeros(4,1);

deta_U=zeros(3,1); pis=[-1.6; -2 ;-3.7; 5]; qis=[-0.8; -1; -1.3]; deta_p=zeros(4,1); deta_q=zeros(3,1); while (kp==0||kq==0) p=zeros(4,1); for i=1:4 for j=1:5

p(i)=p(i)+U(i)*U(j)*(G(i,j)*cos(theta(i)-theta(j))+B(i,j)*sin(theta(i)-theta(j))); end

deta_p(i)=pis(i)-p(i); end for i=1:4

deta_p(i)=deta_p(i)/U(i); end

deta_theta=-inv(B1)*deta_p; for i=1:4

deta_theta(i)=deta_theta(i)/U(i); theta(i)=theta(i)+deta_theta(i); end

eps_theta=max(max(abs(deta_theta)),max(abs(deta_p))); if eps_theta<eps1 kp=1; else

kp=0; end

q=zeros(3,1); for i=1:3 for j=1:5

q(i)=q(i)+U(i)*U(j)*(G(i,j)*sin(theta(i)-theta(j))-B(i,j)*cos(theta(i)-theta(j))); end

deta_q(i)=qis(i)-q(i); end for i=1:3

deta_q(i)=deta_q(i)/U(i); end

deta_U=-inv(B2)*deta_q; for i=1:3

U(i)=U(i)+deta_U(i); end

eps_U=max(max(abs(deta_U)),max(abs(deta_q))); if eps_U<eps1 kq=1;

else

kq=1; end

if max(eps_theta,eps_U)>eps1 k=k+1; end end for i=1:5

U(i)=U(i)*(cos(theta(i))+sqrt(-1)*sin(theta(i))); end disp(U); disp(theta); p=zeros(5,1); q=zeros(5,1); for i=1:5 for j=1:5

p(i)=p(i)+abs(U(i))*abs(U(j))*(G(i,j)*cos(theta(i)-theta(j))+B(i,j)*sin(theta(i)-theta(j))); end end disp(p); for i=1:5 for j=1:5

q(i)=q(i)+abs(U(i))*abs(U(j))*(G(i,j)*sin(theta(i)-theta(j))-B(i,j)*cos(theta(i)-theta(j))); end end disp(q);

y0=[-4i;-8i;-4i;0;0]; s=zeros(5,5); for i=1:5 for j=1:5

s(i,j)=U(i)*(conj(U(i))*conj(y0(i))+(conj(U(i))-conj(U(j)))*conj(-yb(i,j))); end end for i=1:5 s(i,i)=0; end

deta_s=zeros(5,5); for i=1:5 for j=1:5

deta_s(i,j)=s(i,j)+s(j,i); end end

disp(s); disp(deta_s);

3.高斯赛德尔法

clear; clc;

yb=zeros(5,5);

yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i; yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i; yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i; yb(4,4)=(-66.66667i)/2; yb(5,5)=(-33.33333i)/2; yb=yb+conj(yb'); G=real(yb); B=imag(yb);

s=[-1.6-0.8i;-2-i;-3.7-1.3i;5]; U0=[1;1;1;1.05;1.05]; U1=zeros(4,1); deta_U=zeros(4,1); eps1=10^-4; eps_U=1; theta=0; k=0;

while (eps_U>eps1) q=0; for j=1:5

q=q+U0(4)*conj(yb(4,j))*conj(U0(j)); end

q=imag(q);

s(4)=5+sqrt(-1)*q; p=0; for j=1:5

p=p+yb(4,j)*U0(j); end

p=p-yb(4,4)*U0(4);

U1(4)=(conj(s(4))/conj(U0(4))-p)/yb(4,4); theta=angle(U1(4));

U1(4)=1.05*(cos(theta)+sqrt(-1)*sin(theta)); for i=1:3 a=0; b=0; for j=1:i-1

a=a+yb(i,j)*U1(j); end

for j=i+1:4

b=b+yb(i,j)*U0(j); end

U1(i)=(conj(s(i))/conj(U0(i))-yb(i,5)*U0(5)-a-b)/yb(i,i); end for i=1:4

deta_U=U1(i)-U0(i); U0(i)=U1(i); end

eps_U=max(abs(deta_U)); k=k+1; end

disp(abs(U0)); disp(U0);

theta=zeros(5,1); for i=1:5

theta(i)=angle(U0(i)); end

disp(theta);

y0=[-4i;-8i;-4i;0;0]; s=zeros(5,5); for i=1:5 for j=1:5

s(i,j)=U0(i)*(conj(U0(i))*conj(y0(i))+(conj(U0(i))-conj(U0(j)))*conj(-yb(i,j))); end end for i=1:5 s(i,i)=0; end

deta_s=zeros(5,5); for i=1:5 for j=1:5

deta_s(i,j)=s(i,j)+s(j,i); end end disp(s); disp(deta_s);

4.保留非线性法

clear; clc;

yb=zeros(5,5);

yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i; yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i; yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i; yb(4,4)=(-66.66667i)/2;

yb(5,5)=(-33.33333i)/2; yb=yb+conj(yb'); k=0;

eps1=10^-4; unlineeps=1; G=real(yb); B=imag(yb); U=zeros(5,1);

e=[1; 1 ;1 ;1.05 ;1.05]; f=zeros(5,1);

pis=[-1.6; -2 ;-3.7; 5]; qis=[-0.8; -1; -1.3]; deta_p=zeros(4,1); deta_q=zeros(3,1); deta=zeros(8,1); deta_ef=zeros(8,1); deta_e0=zeros(5,1); deta_e1=zeros(5,1); deta_f0=zeros(5,1); deta_f1=zeros(5,1); p=zeros(4,1); q=zeros(3,1); for i=1:4 for j=1:5

p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j))+f(i)*(G(i,j)*f(j)+B(i,j)*e(j)); end

deta_p(i)=pis(i)-p(i); end for i=1:3 for j=1:5

q(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j))-e(i)*(G(i,j)*f(j)+B(i,j)*e(j)); end

deta_q(i)=qis(i)-q(i); end

deta_UU=1.05*1.05-(e(4)*e(4)+f(4)*f(4)); jacobi=zeros(8,8); for i=1:4 for j=1:4

jacobi(2*i-1,2*j-1)=-G(i,j)*e(i)-B(i,j)*f(i); jacobi(2*i-1,2*j)=B(i,j)*e(i)-G(i,j)*f(i); end end for i=1:3 for j=1:4

jacobi(2*i,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i); jacobi(2*i,2*j)=G(i,j)*e(i)+B(i,j)*f(i); end end

for i=1:2:7 jacobi(i,i)=0; jacobi(i,i+1)=0; for j=1:5

jacobi(i,i)=jacobi(i,i)-(G((i+1)/2,j)*e(j)-B((i+1)/2,j)*f(j)); jacobi(i,i+1)=jacobi(i,i+1)-(G((i+1)/2,j)*f(j)+B((i+1)/2,j)*e(j)); end

jacobi(i,i)=jacobi(i,i)-G((i+1)/2,(i+1)/2)*e((i+1)/2)-B((i+1)/2,(i+1)/2)*f((i+1)/2); jacobi(i,i+1)=jacobi(i,i+1)+B((i+1)/2,(i+1)/2)*e((i+1)/2)-G((i+1)/2,(i+1)/2)*f((i+1)/2); end

for i=2:2:6

jacobi(i,i-1)=0; jacobi(i,i)=0; for j=1:5

jacobi(i,i-1)=jacobi(i,i-1)-(G(i/2,j)*f(j)-B(i/2,j)*e(j)); jacobi(i,i)=jacobi(i,i)-(G(i/2,j)*e(j)-B(i/2,j)*f(j)); end

jacobi(i,i-1)=jacobi(i,i-1)+B(i/2,i/2)*e(i/2)-G(i/2,i/2)*f(i/2); jacobi(i,i)=jacobi(i,i)+G(i/2,i/2)*e(i/2)+B(i/2,i/2)*f(i/2); end

jacobi(8,7)=-2*e(4); jacobi(8,8)=-2*f(4); while (unlineeps>eps1) deta_pp=zeros(4,1); for i=1:4 for j=1:5

deta_pp(i)=deta_pp(i)+deta_e0(i)*(G(i,j)*deta_e0(j)-B(i,j)*deta_f0(j))+deta_f0(i)*(G(i,j)*deta_f0(j)+B(i,j)*deta_e0(j));

end end

deta_qq=zeros(3,1); for i=1:3 for j=1:5

deta_qq(i)=deta_qq(i)+deta_f0(i)*(G(i,j)*deta_e0(j)-B(i,j)*deta_f0(j))-deta_e0(i)*(G(i,j)*deta_f0(j)+B(i,j)*deta_e0(j));

end end

deta_UUU=-(deta_e0(4)*deta_e0(4)+deta_f0(4)*deta_f0(4));

deta(i)=deta_p((i+1)/2)-deta_pp((i+1)/2); end

for i=2:2:6

deta(i)=deta_q(i/2)-deta_qq(i/2); end

deta(8)=deta_UU+deta_UUU; deta_ef=-inv(jacobi)*deta; for i=1:4

deta_e1(i)=deta_ef(2*i-1); deta_f1(i)=deta_ef(2*i); end

deta_ppp=zeros(4,1); for i=1:4 for j=1:5

deta_ppp(i)=deta_ppp(i)+deta_e1(i)*(G(i,j)*deta_e1(j)-B(i,j)*deta_f1(j))+deta_f1(i)*(G(i,j)*deta_f1(j)+B(i,j)*deta_e1(j));

end end

deta_qqq=zeros(3,1); for i=1:3 for j=1:5

deta_qqq(i)=deta_qqq(i)+deta_f1(i)*(G(i,j)*deta_e1(j)-B(i,j)*deta_f1(j))-deta_e1(i)*(G(i,j)*deta_f1(j)+B(i,j)*deta_e1(j));

end end

unlineeps=max(max(abs(deta_e1-deta_e0)),max(abs(deta_f1-deta_f0))); unlineeps=max(unlineeps,max(abs(deta_qq-deta_qqq))); unlineeps=max(unlineeps,max(abs(deta_pp-deta_ppp))); for i=1:4

deta_e0(i)=deta_e1(i); deta_f0(i)=deta_f1(i); end k=k+1; end for i=1:4

e(i)=e(i)+deta_e0(i); f(i)=f(i)+deta_f0(i); end for i=1:5

U(i)=e(i)+sqrt(-1)*f(i); end

disp(U);

theta=zeros(5,1); for i=1:5

theta(i)=angle(U(i)); end

disp(theta); p=zeros(5,1); q=zeros(5,1); for i=1:5 for j=1:5

p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j))+f(i)*(G(i,j)*f(j)+B(i,j)*e(j)); end end for i=1:5 for j=1:5

q(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j))-e(i)*(G(i,j)*f(j)+B(i,j)*e(j)); end end disp(p); disp(q);

y0=[-4i;-8i;-4i;0;0]; s=zeros(5,5); for i=1:5 for j=1:5

s(i,j)=U(i)*(conj(U(i))*conj(y0(i))+(conj(U(i))-conj(U(j)))*conj(-yb(i,j))); end end for i=1:5 s(i,i)=0; end

deta_s=zeros(5,5); for i=1:5 for j=1:5

deta_s(i,j)=s(i,j)+s(j,i); end end disp(s); disp(deta_s);

本文来源：https://www.bwwdw.com/article/nia1.html