《两机五节点网络潮流计算方法牛拉法和pq法电力系统稳态分析课程设计.doc》由会员分享，可在线阅读，更多相关《两机五节点网络潮流计算方法牛拉法和pq法电力系统稳态分析课程设计.doc（25页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流两机五节点网络潮流计算方法牛拉法和pq法电力系统稳态分析课程设计.精品文档.电力系统稳态分析课程设计题 目 名 称 两机五节点网络潮流计算方法牛拉法和pq法 学 生 姓 名 学 号 系 、专 业 电气工程系电气工程及其自动化 指 导 教 师 2012年12月28日目 录摘要.第一章 原理简介.3 1.1对潮流分析的简介.3 1.1.1 潮流计算方法分析比较.31.2 MATLAB简介.4 1.2.1 矩阵的运算51.3牛顿拉夫逊法计算潮流分布.6第二章 程序及结果.102.1 设计资料及参数.10 2.1.1 牛顿拉夫逊法的程序框图132.2
2、 用Matlab设计程序.14 2.2.1 程序的编写142.2.2程序运行结果.192.23p_q法程序编写22 总结.32 参考文献.32 电力系统稳态分析课程设计 1.1对潮流分析的简介 潮流分析是研究电力系统的一种最基本和最重要的计算。最初，电力系统潮流计算是通过人工手算的，后来为了适应电力系统日益发展的需要，采用了交流计算台。随着电子数字计算机的出现，1956 年 Ward等人编制了实际可行的计算机潮流计算程序。这样，就为日趋复杂的大规模电力系统提供了极其有力的计算手段。经过几十年的时间，电力系统潮流计算已经发展得十分成熟。潮流计算是研究电力系统稳态运行情况的一种计算，是根据给定的运
3、行条件及系统接线情况确定整个电力系统各个部分的运行状态，如各母线的电压、各元件中流过的功率、系统的功率损耗等等。电力系统潮流计算是计算系统动态稳定和静态稳定的基础。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用电力系统潮流计算来定量的比较供电方案或运行方式的合理性、可靠性和经济性。1.1.1潮流计算方法分析比较 高斯 - 赛德尔潮流计算法原理简单，编程实现容易，特别是对于配网潮流有其独特优势。但是高斯 - 赛德尔潮流计算法在牛顿法以及各种解耦法出现以后似乎成了一种边缘性的方法。牛顿 - 拉夫逊法的优点是收敛速度快，若初值选择较好，算法将具有平方收敛特性，一般迭代 45 次便可以收
4、敛到一个非常精确的解，而且其迭代次数与所计算的网络规模基本无关。牛顿 - 拉夫逊法也具有良好的收敛可靠性 对于呈病态的系统，牛顿-拉夫逊法均能可靠地收敛。牛顿法的缺点是每次迭代的计算量和所需的内存量较大。这是因为雅可比阵元素的数目约为 2(n- 1)2(n- 1)个（直角坐标），且其数值在迭代过程中不断变化。不过，内存占用量及每次迭代所需的时间与程序设计技巧密切相关。牛顿-拉夫逊法的可靠收敛取决于一个良好的启动初值，如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的解点上。对于正常运行的系统，各节电电压一般均在额定值附近，偏移不会太大，并且各节电的相角差也不大，所以对各节电可以采用统
5、一的电压初值。P- Q分解法是为了改进牛顿-拉夫逊法在内存占用量及计算速度方面的不足，P- Q分解法根据电力系统实际运行状态的物理特点，对极坐标形式的牛顿 - 拉夫逊法修正方程式进行了合理的简化。它无论在内存占用量还是计算速度方面都比牛顿-拉夫逊法有较大的改进，是目前计算速度最快的潮流算法。1.2 Matlab 的简介 MATLAB是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率
6、地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。另外，MATLAB提供了一种特殊的工具：工具箱（TOOLBOXES）.这些工具箱主要包括：信号处理（SIGNAL PROCESSING）、控制系统（CONTROL SYSTEMS）、神经网络（NEURAL NETWORKS）、模糊逻辑(FUZZY LOGIC)、小波(WAVELETS)和模拟（SIMULATION）等等。不同领域、不同层次的用户通过相应工具的
7、学习和应用，可以方便地进行计算、分析及设计工作。MATLAB设计中，原始数据的填写格式是很关键的一个环节，它与程序使用的方便性和灵活性有着直接的关系。原始数据输入格式的设计，主要应从使用的角度出发，原则是简单明了，便于修改。 1.2.1 矩阵的运算常数与矩阵的运算即是同该矩阵的每一元素进行运算。但需注意进行数除时，常数通常只能做除数。基本函数运算中，矩阵的函数运算是矩阵运算中最实用的部分，常用的主要有以下几个：det(a) 求矩阵a的行列式eig(a) 求矩阵a的特征值inv(a)或a (-1) 求矩阵a的逆矩阵rank(a) 求矩阵a的秩trace(a) 求矩阵a的迹（对角线元素之和）我们在
8、进行工程计算时常常遇到矩阵对应元素之间的运算。这种运算不同于前面讲的数学运算，为有所区别，我们称之为数组运算。数组的加、减与矩阵的加、减运算完全相同。而乘除法运算有相当大的区别，数组的乘除法是指两同维数组对应元素之间的乘除法，它们的运算符为“.*”和“./”或“.”。前面讲过常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定，数组与常数之间的除法运算没有任何限制。另外，矩阵的数组运算中还有幂运算（运算符为 . ）、指数运算（exp）、对数运算（log）、和开方运算（sqrt）等。有了“对应元素”的规定，数组的运算实质上就是针对数组内部的每个元素进行的。矩阵的幂运算与数组的
9、幂运算有很大的区别。1.3牛顿拉夫逊法计算潮流分布一、牛顿拉夫逊法求解过程大致可以分为以下步骤：（1）形成节点导纳矩阵（2）将各节点电压设初值U（3）将节点初值代入相关求式，求出修正方程式的常数项向量（4）将节点电压初值代入求式，求出雅可比矩阵元素（5）求解修正方程，求修正向量（6）求取节点电压的新值（7）检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步 重新开始进行狭义次迭代，否则转入下一步（8）计算支路功率分布，PV节点无功功率和平衡节点柱入功率。二、直角坐标系计算 (1)牛顿-拉夫逊法潮流计算的公式。把牛顿法用于潮流计算，采用直角坐标形式。其中电压和支路导纳可表示为： （1-
10、1） (2)PQ节点的输出有功功率和无功功率是给定的，则第i节点的给定功率设为和（称为注入功率）。假定系统中的第1、2、m节点为PQ节点，对其中每一个节点的N-R法表达式: =（1、2、m） （1-2） (3)PV节点的有功功率和节点电压幅值是给定的。假定系统中的第m+1、m+2、n-1节点为PV节点，则对其中每一PV节点可以列写方程： （1-3） =（m+1、m+2、n-1）(4)形成雅可比矩阵。 当j=i时，对角元素为（1-4） 当时,矩阵非对角元素为： （1-5）三、极坐标计算 对于潮流计算中待于求出功率的节点功率方程组，在某个近似解附近用泰勒级数展开略去二阶及以上的高阶得到已矩阵表示的
11、修正方程：式中节点n为节点数，m为PV节点，雅阁比矩阵是（2n-2m-2）阶非奇异矩阵，雅可比矩阵各元素表示如下：牛顿拉夫逊极坐标潮流计算的修正方程的迭代方程为：2.1设计资料及参数：课程名称电力系统稳态分析课程设计设计题目两机五节点网络潮流计算牛拉法指导教师刘景霞时间1周一、教学要求电力系统稳态分析课程设计以设计和优化电力系统的潮流分析为重点，提高学生综合能力为目标，尽可能结合实际工程进行。设计内容的安排要充分考虑学校现有的设备，设计时间及工程实际需要，并使学生初步学会运用所学知识解决工程中的实际问题。二、设计资料及参数（一）设计原始资料1、待设计电气设备系统图2、电力系统网络各元件参数3、
12、电力系统电气元件的使用规范4、电力工程电气设计手册（二）设计参考资料1、电力系统稳态分析，陈珩，中国电力出版社，2007，第三版2、电力系统分析，韩祯祥，浙江大学出版社，2005，第三版3、电力系统分析课程实际设计与综合实验，祝书萍，中国电力出版社，2007，第一版三、设计要求及成果1.根据给定的参数或工程具体要求，收集和查阅资料；学习相关软件（软件自选）。2.在给定的电力网络上画出等值电路图。3.运用计算机进行潮流计算。4.编写设计说明书。基本要求：1.编写潮流计算程序；2.在计算机上调试通过（？）；3.运行程序并计算出正确结果（？）；4.写出课程设计报告（包括以下内容）（1份）（1）程序框
13、图；（2）源程序；（3）符号说明表；（4）算例及计算结果5.编写计算说明书（1份）。四、进度安排根据给定的参数或工程具体要求，收集和查阅资料（半天）学习软件（MATLAB或C语言等）（一天半）编程计算复杂系统潮流计算（三天）编写计算设计书（一天）五、评分标准课程设计成绩评定依据包括以下几点：1） 工作态度（占10%）；2） 基本技能的掌握程度（占20%）；3） 程序编写是否合理是否有运行结果(40%)；4） 课程设计说明书编写水平(占30%)。5） 分为优、良、中、合格、不合格五个等级。系统接线图其中节点1为平衡节点，节点2、3、4、5为PQ节点。2.1.1 牛顿拉夫逊法程序框图 电力系统稳态
14、分析课程设计 2.2.1 程序的编写%说明：为了使节点按照先PQ，再PV节点，最后平衡节点的次序编号，以便与公式对照，节点1与节点5对调。%节点阻抗矩阵Z=0,0.04+0.12i,0,0.08+0.24i,0;0.04+0.12i,0,0.06+0.18i,0.06+0.18i,0.02+0.06i;0,0.06+0.18i,0,0.01+0.03i,0.08+0.24i;0.08+0.24i,0.06+0.18i,0.01+0.03i,0,0;0,0.02+0.06i,0.08+0.24i,0,0;%求互导纳for m=1:5 for n=1:5 if Z(m,n)=0 y(m,n)=0;
15、 else y(m,n)=1/Z(m,n); end endend%求导纳for m=1:5 for n=1:5 if m=n Y(m,n)=-y(m,n); else Y(m,n)=sum(y(m,:); end end end%导纳矩阵G=real(Y);B=imag(Y);%计算各节点功率的不平衡量设U=E+jF ；Y=G+Bj；Y U=1,1,1,1,1.06; E=real(U); F=imag(U); %设S=P+Qj； S=-0.60-0.10i;0.20+0.20i;-0.45-0.15i;-0.40-0.05i;0 P=real(S);Q=imag(S); k=0;C=1;w
16、hile C 0.00001E(5)=1.06;F(5)=0;for m=1:4 for n=1:5%计算Pi,Qi,设Pi=Pt;Qi=Qt,按照书上的公式：Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n)+F(m)*(G(m,n)*F(n)+B(m,n)*E(n);Qt(n)=(F(m)*(G(m,n)*E(n)-B(m,n)*F(n)-E(m)*(G(m,n)*F(n)+B(m,n)*E(n);end%设P,Q的改变量为dP,dQdP(m)=P(m)-sum(Pt);dQ(m)=Q(m)-sum(Qt);end%计算Hii,Nii,Jii,Lii,由公式4-41b
17、 左侧公式实现,sum(Ai)，sum(Bi)用于实现公式中的sigerma从j到n的求和;for m=1:4 for n=1:5 Bi(n)=G(m,n)*F(n)+B(m,n)*E(n); Ai(n)=G(m,n)*E(n)-B(m,n)*F(n); end H(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m)+2*G(m,m)*F(m); N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)+2*G(m,m)*E(m); J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m); L(
18、m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m); end%设雅可比矩阵为JJ，以下语句用来实现雅可比矩阵中对角线上元素H N J L 的排列for m=1:4JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m);JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m);end%以下语句用于实现雅可比矩阵非对角线上元素的排列for m=1:4 for n=1:4 if m=n else H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m); N(m,n)=G(m,n)*E
19、(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m); L(m,n)=G(m,n)*F(m)-B(m,n)*E(m);JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n); end endendJJ%设由P,Q的改变量组成的81矩阵为PQ，由E,F的改变量组成的81矩阵为dUfor m=1:4 PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m);enddU=inv(JJ)*PQ;%求逆矩阵C=max(abs(dU);%abs绝对
20、值或复数绝对值for n=1:4 F(n)=F(n)+dU(2*n-1); E(n)=E(n)+dU(2*n);end for n=1:5 U(n)=E(n)+(F(n)*i; endk=k+1;k, dU=dU,PQ,Uend%计算S(5),也就是题目中的S1，即平衡节点功率for m=1:5 I(m)=Y(5,m)*U(m); end S(5)=U(5)*sum(conj(I)Y = 3.7500 -11.2500i -2.5000 + 7.5000i 0 -1.2500 + 3.7500i 0 -2.5000 + 7.5000i 10.8333 -32.5000i -1.6667 + 5
21、.0000i -1.6667 + 5.0000i -5.0000 +15.0000i 0 -1.6667 + 5.0000i 12.9167 -38.7500i -10.0000 +30.0000i -1.2500 + 3.7500i -1.2500 + 3.7500i -1.6667 + 5.0000i -10.0000 +30.0000i 12.9167 -38.7500i 0 0 -5.0000 +15.0000i -1.2500 + 3.7500i 0 6.2500 -18.7500iS = -0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.
22、1500i -0.4000 - 0.0500i 0 JJ = 11.2500 3.7500 -7.5000 -2.5000 0 0 -3.7500 -1.2500 -3.7500 11.2500 2.5000 -7.5000 0 0 1.2500 -3.7500 -7.5000 -2.5000 33.4000 10.5333 -5.0000 -1.6667 -5.0000 -1.6667 2.5000 -7.5000 -11.1333 31.6000 1.6667 -5.0000 1.6667 -5.0000 0 0 -5.0000 -1.6667 38.9750 12.8417 -30.00
23、00 -10.0000 0 0 1.6667 -5.0000 -12.9917 38.5250 10.0000 -30.0000 -3.7500 -1.2500 -5.0000 -1.6667 -30.0000 -10.0000 38.7500 12.9167 1.2500 -3.7500 1.6667 -5.0000 10.0000 -30.0000 -12.9167 38.7500k = 1dU = -0.1076 0.0093 -0.0473 0.0430 -0.0863 0.0154 -0.0922 0.0141PQ = -0.6000 -0.1000 0.5000 1.1000 -0
24、.3750 0.0750 -0.4000 -0.0500U = 1.0093 - 0.1076i 1.0430 - 0.0473i 1.0154 - 0.0863i 1.0141 - 0.0922i 1.0600 JJ = 11.0516 4.3956 -7.3011 -3.3304 0 0 -3.6505 -1.6652 -5.5956 10.8516 3.3304 -7.3011 0 0 1.6652 -3.6505 -7.7040 -2.9621 33.1594 13.0913 -5.1360 -1.9747 -5.1360 -1.9747 2.9621 -7.7040 -12.5804
25、 33.6083 1.9747 -5.1360 1.9747 -5.1360 0 0 -4.9331 -2.1238 38.3848 16.0299 -29.5988 -12.7427 0 0 2.1238 -4.9331 -16.8888 38.0788 12.7427 -29.5988 -3.6876 -1.6135 -4.9168 -2.1513 -29.5009 -12.9079 38.1553 16.2727 1.6135 -3.6876 2.1513 -4.9168 12.9079 -29.5009 -17.0727 38.0553k = 2dU = 0.0032 -0.0131
26、-0.0004 -0.0075 0.0017 -0.0101 0.0021 -0.0108PQ = 0.0164 -0.0636 -0.0770 -0.0220 -0.0008 -0.0317 0.0103 -0.0362U = 0.9963 - 0.1044i 1.0355 - 0.0477i 1.0053 - 0.0846i 1.0033 - 0.0902i 1.0600 JJ = 10.9774 4.3253 -7.2110 -3.2738 0 0 -3.6055 -1.6369 -5.4960 10.6556 3.2738 -7.2110 0 0 1.6369 -3.6055 -7.6
27、466 -2.9466 32.9334 12.9530 -5.0977 -1.9644 -5.0977 -1.9644 2.9466 -7.6466 -12.5843 33.3371 1.9644 -5.0977 1.9644 -5.0977 0 0 -4.8857 -2.0983 38.0494 15.8299 -29.3142 -12.5898 0 0 2.0983 -4.8857 -16.6938 37.6790 12.5898 -29.3142 -3.6499 -1.5922 -4.8665 -2.1230 -29.1988 -12.7380 37.7997 16.0621 1.592
28、2 -3.6499 2.1230 -4.8665 12.7380 -29.1988 -16.8443 37.6306k = 3dU = 1.0e-003 * 0.0241 -0.1731 -0.0025 -0.0892 0.0098 -0.1231 0.0128 -0.1327PQ = 1.0e-003 * -0.0034 -0.8423 -0.5236 -0.2008 -0.0815 -0.3198 0.0241 -0.3906U = 0.9961 - 0.1044i 1.0354 - 0.0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600 JJ =
29、 10.9764 4.3243 -7.2098 -3.2731 0 0 -3.6049 -1.6366 -5.4951 10.6529 3.2731 -7.2098 0 0 1.6366 -3.6049 -7.6459 -2.9464 32.9307 12.9517 -5.0973 -1.9643 -5.0973 -1.9643 2.9464 -7.6459 -12.5839 33.3340 1.9643 -5.0973 1.9643 -5.0973 0 0 -4.8851 -2.0980 38.0451 15.8278 -29.3106 -12.5883 0 0 2.0980 -4.8851
30、 -16.6919 37.6740 12.5883 -29.3106 -3.6494 -1.5920 -4.8658 -2.1227 -29.1950 -12.7363 37.7952 16.0599 1.5920 -3.6494 2.1227 -4.8658 12.7363 -29.1950 -16.8421 37.6252k = 4dU = 1.0e-007 * 0.0232 -0.2889 -0.0013 -0.1352 0.0054 -0.1876 0.0087 -0.2038PQ = 1.0e-006 * -0.0240 -0.1486 -0.0424 -0.0189 -0.0165
31、 -0.0415 -0.0042 -0.0525U = 0.9961 - 0.1044i 1.0354 - 0.0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600 S = -0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.1500i -0.4000 - 0.0500i 1.2982 + 0.2445i2.2.3pq法程序编写基于MATLAB 的PQ 分解法潮流计算程序MATLAB 已广泛应用于自动控制、数学运算、信号分析、计算机技术、图像信号处理、财务分析、航天工业和生物医学工程等领域。由于MATLAB语言功
32、能强大、人际界面友好、编程效率高、强大而智能化的作图功能，且具有编程语句简洁、灵活、表达和运算能力强等显著特点。程序清单%本程序的功能是用PQ分解法进行潮流计算n=input(请输入节点数：n=);nl=input(请输入支路数：nl=);isb=input(请输入平衡母线节点号：isb=);pr=input(请输入误差精度:pr=);B1=input(请输入由支路参数形成的矩阵:B1=);B2=input(请输入由节点参数形成的矩阵:B2=);X=input(请输入由节点号及其对地阻抗形成的矩阵:X=);na=input(请输入PQ节点号:na=);Y=zeros(n);YI= zeros(
33、n);e= zeros(1,n);f= zeros(1,n);V= zeros(1,n);O= zeros(1,n);for i=1:nif X(i,2)=0; p=x(i,1); Y(p,p)=1./X(i,2); endendfor i=1:nlif B1(i,6)= =0 p=B1(i,1);q=B1(i,2);else p= B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5);YI(p,q)= YI(p,q) -1./B1(i,3);Y(q,p) =Y(p,q);YI(q,p) = YI(p,q);Y(q,q)= Y(q,q)+
34、1./(B1(i,3)* B1(i,5)2)+B1(i,4)./2;YI(q,q)= YI(q,q) +1./B1(i,3);Y(p,p)= Y(p,p) +1./B1(i,3)+ B1(i,4)./2;YI(p,p)= YI(p,p) +1./B1(i,3);end %求导纳矩阵G=real(Y);B=imag(YI);BI=imag(Y);for i=1:n S(i)=B2(i,1)-B2(i,2); BI(i,i)= BI(i,i)+ B2(i,5);endP=real(S);Q=imag(S);for i=1:ne(i)=real(B2(i,3);f(i)=imag(B2(i,3);V
35、(i)=B2(i,4);endfor i=1:nif B2(i,6)= =2 V(i)=sqrt(e(i)2+f(i)2); O(i)=atan(f(i)./e(i);endendfor i=2:nif i= =n B(i,i)=1./B(i,i);else IC1=i+1; for j1=IC1:n B(i,j1)= B(i,j1)./B(i,i); end B(i,i)=1./B(i,i); for k=i+1:n for j1=i+1:n B(k,j1)=B(k,j1)-B(k,i)*B(i,j1); end endendendp=0;q=0;for i=1:nif B2(i,6)= =
36、2 p=p+1;k=0; for j1=1:n if B2(i1,6)= =2 k=k+1; A(p,k)=BI(i,j1); end endendendfor i=1:naif i= =na A(i,i)=1./A(i,i);else k=i+1; for j1=k:na A(i,j1)= A(i,j1)./A(i,i); end A(i,i)=1. /A(i,i); for k=i+1:na for j1=i+1:na A(k,j1)= A(k,j1)-A(k,i)*A(i,j1); endendendendICT2=1;ICT1=0;kp=1,kq=1;K=1;DET=0;ICT3=1;
37、while ICT2=0|ICT3=0ICT2=0;ICT3=0; for i=1:nif i=isb C(i)=0; for k=1:n C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k)+BI(i,k)*sin(O(i)-O(k); end DP1(i)=P(i)-V(i)*C(i); DP(i)=DP1(i)./V(i); DET=abs(DP1(i); if DET=pr ICT2=ICT2+1; endendendNp(K)=ICT2; if ICT2=0 for i=2:n DP(i)=B(i,i)*DP(i); if i =n IC1=i+1; for k=ic1:n DP(k)=DP(k)-B(k,i)*DP(i); end else for LZ=3:i L=i+3-LZ; IC4=L-1; for MZ=2:IC4 I=IC4+2-MZ; DP(I)=DP(I)-B(I,L)*DP(L); end end endendfor i=2:n O(i)=O(i)-DP(i);endkq=1;L=0;for i=1:n if B2(i,6)= =2C(i)=0;L=L+1;for k=1:n C(i)= C(i) +V(k)*(G(i,k)*sin(O(i)-O(k)+BI(i,k)*