设为首页收藏本站|繁體中文 快速切换版块

电力研学论坛

 找回密码
 立即加入

QQ登录

只需一步,快速开始

搜索
查看: 1787|回复: 5
收起左侧

[分享] 牛拉法潮流程序

[复制链接]
最佳答案
0 

该用户从未签到

发表于 2010-7-15 09:31:48 | 显示全部楼层 |阅读模式
快速获得阅读权限、下载权限、提升会员组攻略【新手必看】2016-03-03更新

马上加入,结交更多好友,共享更多资料,让你轻松玩转电力研学社区!

您需要 登录 才可以下载或查看,没有帐号?立即加入

x
希望对大家有所帮助哦
  1. clc
  2. disp('此程序为牛拉法解潮流')
  3. nPQ=input('请输入PQ节点的个数:');
  4. nPV=input('请输入PV节点的个数:');
  5. n=nPQ+nPV+1;
  6. Ps=[0;-0.5;0.2];
  7. Qs=[0;-0.3];
  8. Us=[1.0+j*0;1.0+j*0;1.05+j*0;1.05+j*0];
  9. % nl nr R X Bl Br
  10. zdata=[1 2 0 0.1880 -0.6815 0.6040;
  11. 1 3 0.1302 0.2479 0.0129 0.0129;
  12. 1 4 0.1736 0.3306 0.0172 0.0172;
  13. 3 4 0.2603 0.4959 0.0259 0.0259];
  14. % nx B
  15. xdata=[2 0.05];
  16. dPQU=1;
  17. %计算导纳矩阵
  18. nl=zdata(:,1);
  19. nr=zdata(:,2);
  20. R=zdata(:,3);
  21. X=zdata(:,4);
  22. Bl=zdata(:,5);
  23. Br=zdata(:,6);
  24. nx=xdata(:,1);
  25. Bx=xdata(:,2);
  26. nbr=length(zdata(:,1));
  27. nbrx=length(xdata(:,1));
  28. Z=R+j*X;
  29. y=ones(nbr,1)./Z;
  30. Y=zeros(n,n);
  31. %计算非对角元素
  32. for ii=1:nbr
  33. Y(nl(ii),nr(ii))= Y(nl(ii),nr(ii))-y(ii);
  34. Y(nr(ii),nl(ii))= Y(nl(ii),nr(ii));
  35. end
  36. %计算对角元素
  37. for ii=1:n
  38. for jj=1:nbr
  39. if nl(jj)==ii|nr(jj)==ii
  40. Y(ii,ii)=Y(ii,ii)+y(jj);
  41. end
  42. end
  43. end
  44. for ii=1:nbr
  45. Y(nl(ii),nl(ii))= Y(nl(ii),nl(ii))+j*Bl(ii);
  46. Y(nr(ii),nr(ii))= Y(nr(ii),nr(ii))+j*Br(ii);
  47. end
  48. for ii=1:nbrx
  49. Y(nx(ii),nx(ii))=Y(nx(ii),nx(ii))+j*Bx(ii);
  50. end
  51. %分离G、B
  52. G=real(Y);
  53. B=imag(Y);
  54. disp('导纳矩阵:');
  55. Y
  56. e=real(Us);
  57. f=imag(Us);
  58. k=0;
  59. while dPQU>0.00001
  60. %求dP
  61. dP=zeros(n-1,1);
  62. for ii=1:n-1
  63. t=0;
  64. for jj=1:n
  65. t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
  66. end
  67. dP(ii)=Ps(ii)-real(t*(e(ii)+j*f(ii)));
  68. end
  69. %求dQ
  70. dQ=zeros(nPQ,1);
  71. for ii=1:nPQ
  72. t=0;
  73. for jj=1:n
  74. t=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));
  75. end
  76. dQ(ii)=Qs(ii)-imag(t*(e(ii)+j*f(ii)));
  77. end
  78. %求dU^2
  79. dU2=zeros(nPV,1);
  80. ii=1:nPV;
  81. dU2(ii)=abs(Us(ii+nPQ))^2-abs(e(ii+nPQ)+j*f(ii+nPQ))^2;
  82. dS=[dP;dQ;dU2];
  83. dPQU=max(abs(dS));
  84. if(dPQU>0.00001)
  85. k=k+1
  86. %形成雅克比行列式
  87. Jacob=zeros(2*(n-1),2*(n-1));
  88. %P部分
  89. for ii=1:n-1
  90. mid1=0;
  91. mid2=0;
  92. for jj=1:n
  93. if ii~=jj&&jj<n
  94. Jacob(ii,2*jj-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));
  95. Jacob(ii,2*jj)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
  96. end
  97. mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
  98. mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
  99. end
  100. Jacob(ii,2*ii-1)=-mid2-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);
  101. Jacob(ii,2*ii)=-mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
  102. end
  103. %Q部分
  104. for ii=1:nPQ
  105. mid1=0;
  106. mid2=0;
  107. for jj=1:n
  108. if ii~=jj&&jj<n
  109. Jacob(ii+n-1,2*jj-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);
  110. Jacob(ii+n-1,2*jj)=G(ii,jj)*e(ii)+B(ii,jj)*f(ii);
  111. end
  112. mid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj);
  113. mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj);
  114. end
  115. Jacob(ii+n-1,2*ii-1)=mid1+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);
  116. Jacob(ii+n-1,2*ii)=-mid2+G(ii,ii)*e(ii)+B(ii,ii)*f(ii);
  117. end
  118. %U2部分
  119. for ii=nPQ+1:n-1
  120. Jacob(ii+n-1,2*ii-1)=-2*e(ii);
  121. Jacob(ii+n-1,2*ii)=-2*f(ii);
  122. end
  123. dU=-inv(Jacob)*dS;
  124. de=zeros(n-1,1);
  125. df=zeros(n-1,1);
  126. ii=1:n-1;
  127. de(ii)=dU(2*ii-1);
  128. df(ii)=dU(2*ii);
  129. e(ii)=e(ii)+de(ii);
  130. f(ii)=f(ii)+df(ii);
  131. end
  132. end%迭代结束
  133. U=e+j*f;
  134. %计算PV节点的Q
  135. P=zeros(n,1);
  136. Q=zeros(n,1);
  137. for ii=1:nPV
  138. t=0;
  139. for jj=1:n
  140. t=t+conj(Y(ii+nPQ,jj))*(e(jj)-j*f(jj));
  141. end
  142. Q(ii+nPQ)=imag(t*(e(ii+nPQ)+j*f(ii+nPQ)));
  143. end
  144. %计算平衡节点
  145. t=0;
  146. for jj=1:n
  147. t=t+conj(Y(n,jj))*(e(jj)-j*f(jj));
  148. end
  149. P(n)=real(t*(e(n)+j*f(n)));
  150. Q(n)=imag(t*(e(n)+j*f(n)));
  151. ii=1:n-1;
  152. P(ii)=Ps(ii);
  153. ii=1:nPQ;
  154. Q(ii)=Qs(ii);
  155. %计算线路潮流
  156. Sij=zeros(nbr,1);
  157. Sji=zeros(nbr,1);
  158. dSij=zeros(nbr,1);
  159. for ii=1:nbr
  160. Sij(ii)=U(nl(ii))*(conj(U(nl(ii)))*(-j*Bl(ii))+(conj(U((nl(ii))))-conj(U((nr(ii)))))*conj(y(ii)));
  161. Sji(ii)=U(nr(ii))*(conj(U(nr(ii)))*(-j*Br(ii))+(conj(U((nr(ii))))-conj(U((nl(ii)))))*conj(y(ii)));
  162. dSij(ii)=Sij(ii)+Sji(ii);
  163. end
  164. nn=[1:n]';
  165. disp(' n e f P Q');
  166. Display1=[nn e f P Q]
  167. disp(' nl nr Sij Sji dSij');
  168. Display2=[nl nr Sij Sji dSij]

复制代码
楼主热帖
真诚赞赏,手留余香
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
回复 呼我

使用道具 举报

最佳答案
0 

该用户从未签到

发表于 2010-12-11 19:35:42 | 显示全部楼层
我正好要用,虽然不知道这个能不能行,先感谢楼主分享啦
真诚赞赏,手留余香
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
最佳答案
0 

该用户从未签到

发表于 2010-12-24 11:19:59 | 显示全部楼层
不行吧,不能保证雅可比矩阵是非奇异矩阵,不能直接求逆吧
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
最佳答案
0 
  • TA的每日心情
    奋斗
    2017-8-18 11:17
  • 签到天数: 52 天

    连续签到: 1 天

    [LV.5]常住居民I

    发表于 2016-3-29 16:46:35 | 显示全部楼层
    正在学习中
    真诚赞赏,手留余香
    [发帖际遇]: 花尹儿特爱帮助研友,研友一致同意奖励他 学分1 点,帅呆了. 幸运榜 / 衰神榜
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
    最佳答案
    0 
  • TA的每日心情
    奋斗
    2018-12-28 08:24
  • 签到天数: 11 天

    连续签到: 4 天

    [LV.3]偶尔看看II

    发表于 2018-12-27 08:50:46 | 显示全部楼层
    原创,永远支持~~~~~~~~~~~~
    真诚赞赏,手留余香
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
    您需要登录后才可以回帖 登录 | 立即加入

    本版积分规则

    招聘斑竹

    小黑屋|手机APP(beta)|Archiver|电力研学网 ( 赣公网安备 36040302000210号 赣ICP备12000811号-1) |网站地图

    GMT+8, 2019-4-22 09:25

    Powered by Discuz! X3.4 Licensed

    © 2001-2017 Comsenz Inc.

    快速回复 返回顶部 返回列表