第一篇：第5章高等数学计算的MATLAB实现讲稿
第5章 高等数学计算的MATLAB实现
高等数学是大学数学学习的基本内容。利用MATLAB的符号工具箱，可以解决极限、导数、微分、积分、级数和微分方程等方面的问题。
5.1 函数和极限
5.1.1 函数
使用符号表达式，可以进行复合函数运算和反函数运算，下面分别予以介绍。
⒈ 复合函数运算
在MATLAB中，符号表达式的复合函数运算主要是通过函数compose来实现的。compose函数的调用格式如下：
●compose(f,g)：返回复合函数f(g(y))。在这里ff(x),gg(y)。其中，x是findsym定义的f函数的符号变量，y是findsym定义的g函数的符号变量。
●compose(f,g,z)：返回自变量为z的复合函数f(g(z))。在这里，f=f(x),g=g(y)，x，y分别是findsym定义的f函数和g函数的符号变量。
●compose(f,g,x,z)：返回复合函数f(g(z))，并且使x成为f函数的独立变量。如：f=cos(x/t),compose(f,g,x,z)返回cos(g(z)/t)，而compose(f,g,t,z)返回cos(x/g(z))。
●compose(f,g,x,y,z)：返回复合函数f(g(z))，并且使x与y分别成为函数f与g的独立变量。
例5－1 将fx和xtany复合到一个函数中，指定x和y为它们的独立变量，自变量为z。
程序：
syms x y t z;g=tan(y);f=x^t;compose(f,g,x,y,z)%求复合函数 运行结果： ans = tan(z)^t ⒉ 符号表达式的反函数运算 在MATLAB中，符号表达式的反函数运算主要是通过函数finverse来实现的。finverse函数的调用格式如下。
5.2 导数
5.2.1 求函数的导数
在MATLAB中，微分和求导都可以由函数diff实现。diff函数可同时处理数值和符号两种情况下的求导和微分。该函数的调用格式如下所示。
●diff(F)：对findsym函数返回独立变量求微分，F为符号表达式。●diff(F,'a')：对a变量求微分，F为符号表达式。
●diff(F,n)：对findsym函数返回的独立变量求n次微分，F为符号表达式。●diff(F,'a',n)或diff(F,n,'a')：对变量a求n次微分，F为符号表达式。
(x1)5例5－4 求f(x)的二阶导数。
x1程序： syms x f=(x-1)^5/(x+1);df=diff(f,1);%求导数 d2f=diff(f,2);df=simplify(df)%化简 d2f=simplify(d2f)运行结果： df = 2*(x-1)^4*(2*x+3)/(x+1)^2 d2f = 4*(x-1)^3*(3*x^2+9*x+8)/(x+1)^3 5.2.2 求隐函数的导数
例5－5 求隐函数
1F(x,y)xysiny
2所确定的导数dy。dx程序：
%求隐函数的导数
f=sym('x-y+1/2*sin(y)');fx=diff(f,'x');fy=diff(f,'y');
5.3.1 渐近线
求函数图形的水平渐近线，需要求x趋于无穷时f的极限，即
limit(f,inf)limit(f,-inf)求f的垂直渐近线，使分母等于0，用下面的命令进行求解。roots=solve(denom)%返回方程xx30的解。综合程序： syms x num=3*x^2+6*x-1;denom=x^2+x-3;f=num/denom;a=limit(f,inf);b=double(a);roots=solve(denom);ezplot(f)%符号函数作图命令。
hold on%在原有的图形上面叠加图形。plot([-2*pi 2*pi],[b b], 'g')%绘水平渐近线
plot(double(roots(1))*[1 1],[-5 10], 'r')%绘垂直渐近线 plot(double(roots(2))*[1 1],[-5 10], 'r')%绘垂直渐近线 title('水平渐近线和垂直渐近线')hold off%取消图形叠加 运行结果见图5－1。5.3.2 极值
23x26x1从图5－1可以看出，函数f(x)至少有2个极值点，求解
x2x3程序：
syms x num=3*x^2+6*x-1;denom=x^2+x-3;f=num/denom;a=limit(f,inf);b=double(a);roots=solve(denom);ezplot(f)%符号函数作图命令。
hold on%在原有的图形上面叠加图形。plot([-2*pi 2*pi],[b b], 'g')%绘水平渐近线
图5－2 表示函数的渐近线和极值
5.3.3 拐点
求函数的拐点，需要先求函数的2阶导数，后面的处理方法与求极值方法相似。
5.4 不定积分和定积分
MATLAB中，用符号工具箱的int函数求函数的不定积分和定积分。int函数的调用格式如下所示。
●int(F)：对findsym函数返回独立变量求不定积分，F为符号表达式。
●int(F,v)：对v变量求不定积分，F为符号表达式。
●int(F,a,b)：对findsym函数返回独立变量求从a到b的定积分，F为符号表达式。
●int(F,v,a,b)：对v变量求从a到b的定积分，F为符号表达式。
5.4.4 定积分的应用
例5－11 计算由两条抛物线y2x，yx2所围成的图形的面积。程序：
%求曲线的交点
[x1,y1]=solve('y^2=x','y=x^2');x1=double(x1);y1=double(y1);n=numel(x1);%下面寻找实数解 m=1;x0=[];y0=[];for k=1:n
if isreal(x1(k))&&isreal(x1(k))
x0(m)=x1(k);y0(m)=y1(k);
m=m+1;
end end x0=sort(x0);y0=sort(y0);%排序 %下面计算定积分 syms x f=sqrt(x)-x^2;A=int(f,x,x0(1),x0(2))运行结果： A = 1/3 x2y2例5－12 计算由椭圆221所围成的图形绕x轴旋转而成的旋转ab体的体积。
V程序：
syms a b x f=pi*b*b*(a*a-x*x)/a/a;V=int(f,x,-a,a)运行结果： V = 4/3*pi*b^2*a 例5－13 计算由曲线yaab2a2(a2x2)dx
23/2x上相应于x从a到b的一段弧的长度。3-9
0 d =
-3 f =
-1 g =
-3 例5－16 已知三点M(1,1,1)、A(2,2,1)和B(2,1,2)，求AMB。程序：
M=[1 1 1];A=[2 2 1];B=[2 1 2];ma=A-M;mb=B-M;c=dot(ma,mb)/sqrt(dot(ma,ma))/sqrt(dot(mb,mb));amb=acos(c)运行结果： amb =
1.0472 例5－17 已知三角形ABC的顶点是A(1,2,3)、B(3,4,5)和C(2,4,7)求三角形ABC的面积。
程序：
A=[1 2 3];B=[3 4 5];C=[2 4 7];ab=B-A;ac=C-A;S=sqrt(dot(cross(ab,ac),cross(ab,ac)))/2 运行结果： S =
3.7417 5.5.2 曲面及其方程
利用MATLAB提供的绘图函数，可以绘制给定函数的曲面。相关内容在前面已介绍过，请参见4.2.4小节。
5.6 多元函数的极限和求导
对于函数有多个变量的情况，求极限和导数时需要指定函数对哪个变量进行求取。在MATLAB中仍然使用limit和diff函数求多元函数的极限和导数。5.6.1 求多元函数的极限
例5－18 求极限limsin(xy)sin(x)。
y0y-11●symsum(s)：求符号表达式s相对于符号变量k的和，k由findsym函数确定，取值从0到k-1。
●symsum(s,v)：求符号表达式s相对于符号变量v的和，v从0到v-1。●symsum(s,a,b)和symsum(s,v,a,b)：指定符号表达式s从v=a累加到v=b。
1例5－21 求级数k、2和xk(x1)。
i0k1kk0n1程序：
syms x k n s1=symsum(n)s2=symsum(1/k^2,1,inf)s3=symsum(x^k,k,0,inf)运行结果： s1 = 1/2*n^2-1/2*n s2 = 1/6*pi^2 s3 =-1/(x-1)5.7.2 泰勒级数展开
用taylor函数进行泰勒级数展开。该函数的调用格式如下：
●taylor(f,n,v)：返回f的n-1阶马克劳林多项式近似。f为表示函数的符号表达式，v指定表达式中的独立变量。v可以是字符串或符号变量。
●taylor(f,n,v,a)：返回f关于a的n-1阶泰勒级数近似。变量a可以是数值、符号或表示数值值或未知值的字符串。n，v和a的顺序没有先后之分。taylor函数根据变量的位置和类型确定它们的用途。还可以忽略n,v,a等变量中的任何一个。如果不确定v，taylor函数用findsym函数确定函数的独立变量。n的默认值为6。
f(n)(a)泰勒级数：f(x)(xa)n。
n!n01例5－22 求函数f(x)的泰勒级数展开，取前9项。
54cosx程序： syms x f=1/(5+4*cos(x));t=taylor(f,9)运行结果：
3d2y2cos2xydx例5－24 求解微分方程yx01
dyx00dx程序：
y=dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0','x');%求解微分方程 y=simplify(y)%化简y的形式。运行结果： y = 4/3*cos(x)-2/3*cos(x)^2+1/3
习题五
1、完成实验指导书中的实验五的上半部分。
第二篇：计算讲稿
如何做好有效的计算教学
浅谈如何提高小学生数学计算能力
福州路小学 张成勋
《新课程标准》关于小学数学计算指出“使学生能够正确地进行整数、小数、分数的四则计算,对于其中一些基本的计算,要达到一定的熟练程度,并逐步做到计算方法合理、灵活”。
小学计算教学贯穿于小学数学教学的全过程, 计算教学直接关系着学生对数学基础知识与基本技能的掌握，关系着学生观察，记忆，注意，思维等能力的发展，关系着学生的学习习惯，情感，意志等非智力因素的培养。“计算”在教学中所占的比重相当大，无论是应用题、统计知识，还是几何题、简易方程，都离不开计算。计算的准确率和速度如何，将直接影响学生学习的质量，因此，小学阶段的计算教学就显得异常重要。很多数学教师都有同感:每次批阅学生作业或试卷,总会看到不少学生在计算中出现错写,多写,少写数字或符号,打错小数点等现象。下面就提高计算能力的实效性，我谈一点个人的体会。
一、激发学生的计算兴趣,提高学生的计算能力
1.适时列举中外数学家的典型事例。教学中教师可以适时列举中外数学家的典型事例,来激发学生对数学的兴趣和爱好,提高计算效果。比如:我国著名数学家陈景润为了攻克歌德巴赫猜想,草稿纸演算了几麻袋；“会移动的黑板”中上街散步还在马车的后背箱上专心致志演算的安培,这些生动典型的事例都可以激发学生对计算的兴趣。
2.有效的创设数学情境。计算比较枯燥,因此,教学时要根据小学生的心理特点,将童话、游戏、比赛等融入课堂教学中,同时要注意题目的灵活性、注意练习形式的多样性,从而激发小学生的计算兴趣,提高计算能力。
3.成立数学兴趣小组。成立数学兴趣小组不仅可以丰富小学生的课外生活,而且可以调动学生学习的积极性、主动性和创造性。可以定期或不定期举办数学讲座、速算、巧算比赛,从而使学生达到算得准、算得巧的目的,增强计算情趣。
二、使学生理解和牢固掌握有关基础知识
学生的计算离不开数学概念、运算定律、运算性质、运算法则和计算公式等内容。对学生不易理解的某些计算法则，往往成为教学的难点。在教学中教师不能急于求成，应帮助学生以掌握基础知识为突破口，分散、突破难点。例如教学同分母分数加减法时，首先要让学生领会分母相同即分数单位不同，而分数单位相同，就能直接相加减，懂得了这个道理，例如教学异分母分数加减法时，首先要让学生领会分母不同即分数单位不同，而分数单位不同，就不能直接相加减，懂得了这个道理，再引导学生运用通分的知识，化异分母分数为同分母分数，于是问题就转化为已学过的同分母分数相加减了。
乘法分配律的理解和掌握对我们计算四则混合运算很重要。怎样 让学生从实质上理解乘法分配律？
在乘法分配律的教学中，如果只求形式不求实质理解，一方面从认识的角度看是不严谨的，另一方面很容易造成学生不求甚解、囫囵吞枣的不良认知习惯。如果满足于从形式上掌握乘法分配律，对于学生的后续发展也极为不利。因此，在教学时结合生活中的购物问题。先出示了这样一道例题：一张桌子65元，一张椅子35元，学校买了15套课桌椅，一共要花多少元？学生用了两种解答方法65×15 +35×15，（65+35）×15，借助对同一实际问题的不同解决方法让学生体会乘法分配律的合理性，体会分配律中“分配”的含义
三、理解算理,优化算法，提高计算技巧
新课标要求在教学中要关注学生的学习过程，那么我们在计算教学中也应该关注过程教学。在计算过程中，理解算理很重要，他是计算的前提，而算法优化则是计算的关键。学生计算错误的原因常常是算理在学习的过程中没有理解到位。在计算教学中根据知识体系之间的联系可以在迁移中帮助学生理解算理。例如教学除数是小数的除法，学生已经学习了除数是整数的除法，积累了以下的两点认识：计算时就按整数除法的方法算出结果；商的小数点和被除数的小数点对齐。这些认识是学生学习除数是小数除法的基础，在实际教学中教师可以先复习除数是整数的除法，如：“38.4÷24”，在学生明确商的小数点是如何确定后，把复习题改成“3.84÷2.4”，在学生尝试计算中着重引导学生分析怎样把除数是小数转化成除数是整数的除法，在学生初步理解算理的基础上进行“试一试”的教学：
0.12÷3= 0.12÷0.03=
学生在两组题目的练习比较中发现：先运用商不变规律把除数是小数的转化成除数是整数的除法，再按除数是整数的除法的方法来计算。如果教师直接通过例题的教学就让学生尝试计算，学生将缺少再次理解算理的机会。所以“试一试”的教学为学生提供了自主迁移的机会，对学生更深刻地理解算理是十分必要的。在教学整数加法时要让学生理解为什么要“相同数位要对齐，相加满十向前一位进一”的算理。在教学小数加法时要让学生理解为什么小数点要对齐。
学生在理解算理的基础上，应该引导学生自主探索和合作交流出算法的多样化，教师还应该引领学生寻找最优的算法。以这种算法为主进行训练，从而提高学生的计算能力。
例如：在两位数乘整十数探索“24×10”的口算方法时，有的学生联系情境图，先算9箱有多少瓶：24×9=216，再加1箱的24瓶：216+24=240；先算5箱有多少瓶：24×
5=120再算10箱有多少瓶：120×2=240；把每箱中的24瓶分成20瓶和4瓶，先算10个20瓶是200瓶，再算4个10瓶是40瓶，再用200+40=240；还有利用24×1=24迁移出24×10=240。在发散的基础上引导学生着重理解最后一种算法“24乘1个十得24个十就是240”，在比较中引领学生进行算法的优化，在练习中重点运用这种算法，从而提高学生的计算能力。
为了计算简便，解题中要训练学生合理运用运算定律，灵活解题，增强计算技巧。如计算3.4×0.125 +4.6×0.125，学生一眼就能看出运用乘法分配律可以得出（3.4+4.6）×0.125。教学时，教师不应就此满足，可进一步深化，充分挖掘学生的潜能，如依次出示：1.25× 0.34 + 4.6 × 0.125 3.4 ÷8+ 4.6 × 0.125 这样，学生也就不会一遇到稍有变化的题目就不会解，同时学生的思维也得到了训练。教学中要减少学生计算的错误，提高计算的正确率，应根据学生的实际情况，因材施教，因人施教，采取相应的对策，才能提高学生计算的能力。
四、加强练习和基本技能训练（精讲精练）
教师设计练习时最好分层进行，形式多样。特别是练习的内容要注意有针对性，有层次，有梯度，练习的形式要多样，学生在进行计算练习时才不会觉得枯燥，才会觉得有兴趣。在设计练习题时要注意围绕重点与难点来设计一些有针对性的练习，尽量让学生能够练习有所收获。但要注意练习的数量要有个度，不能只要量不讲质，搞题海战术，就会适得其反。部分学生本身缺乏勤奋学习的精神，再加上计算本身又枯燥乏味，缺乏情节，学生遇到题量较多时，易产生抵触情绪，不愿计算，严重的可影响学生对学习数学的兴趣，教学中，首先要注意对题量的控制，其次计算形式要多样，除了计算题，可适当增加一些判断、选择题，趣味题，如：（）+（）－（）+（）=0。这样既减轻了学生的学业负担，又增加了学生的学习兴趣。
1、加强口算基本训练,提高学生的计算能力
任何一道四则混合运算题都是有若干道口算题综合而成的,口算的正确、迅速与否
直接关系到计算能力的提高。设计口算练习时,要有针对性,由易到难,逐步提高,包括一些简便运算题,经常进行口算练习,做到每堂课练习5分钟，有利于提高学生的计算能力。
2、加强简算训练，提高计算效率。简便计算是小学计算教学的重要组成部分，它要求学生充分运用学过的运算定律、性质、公式，合理改变运算的数据及运算顺序，使计算尽可能简便、快捷，提高计算效率。因此，在教学中，必须加强简算训练，逐步增强简算意识，提高简算能力。
计算中，学生容易套用、滥用一些性质、定律，要让学生进行一些对比练习，自己诊断错误，反思计算出错的症结点，防止再次出现同样的错误。如：300－175+25，300－175－25；125×8÷125×8，（125×8）÷（125×8）；„„让学生辨析，什么情况下运用性质、定律可以简便，明白为什么有些可以用，有些题不能用。
3、加强估算教学,提高学生的计算能力
估算可以培养学生的“数感”,可以引导学生深入理解“运算”,可以帮助学生检查计算的结果正确与否,运用估算的方法可以对计算的结果做预先定位,快速地确定计算结果的取值范围,通过计算前的估算和计算后的检查,可以避免由于粗心大意造成的错误。如:个位是2和7的两个数,计算结果的个位如果相加就肯定是9,相乘就一定是4。
4、重视错题分析,提高学生的计算能力
针对计算中存在的问题和一些常见而又典型的错例,教师要与学生一起分析、交流,通过集体“会诊”,达到既“治病”又“预防”的目的。至于那些形近又易错的题目,教师应该组织学生经常进行对比练习,培养学生比较、鉴别的能力。
五、培养学生良好的学习习惯,提高学生的计算能力
学生的计算错误,表面上看是“粗心”造成的,可是这些“粗心”其实就是没有养成良好的学习习惯。所以在平时的计算训练时,要让学生养成一对、二审、三算、四查的好习惯。
1.对:就是认真核对。题目都抄错了,结果又怎么能正确呢?所以,要求学生在抄题和每步计算时,都应当及时与原题或上一步算式进行核对,以免抄错数或运算符号。
2.审:就是认真审题。引导学生在做计算题时,不应拿起笔来就下手算,必须先审题,弄清这道题应该先算什么,再算什么,有没有简便的计算方法,然后才能动笔算。
3.算:就是认真书写、计算。作业、练习的书写都要工整,不能潦草,格式一定要规范,对题目中的数字、小数点、运算符号的书写尤其要符合规范,数字间有适当的间隔,草稿上的竖式也要数位对齐、条理清楚,计算时精力集中,不急不抢。另外,计算必须先求准,再求快。
4.查:就是仔细检查。计算完,首先要检查计算方法是不是合理;其次,检查数字、符号会不会抄错,小数点会不会错写或漏写;再次,对计算中途得到的每一个得数和最后的结果都要进行检查和验算。
总之，在计算教学中也要关注过程，让学生在理解算理的基础上掌握方法，掌握计算法则，优化算法。这样学生经历了过程，明白了道理，对学生的计算技能和创新能力的培养有好处，有利于学生的后续发展。
第三篇：MATLAB计算24点游戏代码
clear,closeall clc a=5;b=7;c=10;
d=4;%这里输入需要计算的四个数字a,b,c,d f=[a b c d];tic;g=perms(f);[m,n]=size(g);h='+-*/';fori=1:24 for k1=1:4 for k2=1:4 for k3=1:4
str11=[num2str(g(i,1)),h(k1),num2str(g(i,2)),h(k2),num2str(g(i,3)),h(k3),num2str(g(i,4))];
str22=['(',num2str(g(i,1)),h(k1),num2str(g(i,2)),')',h(k2),num2str(g(i,3)),h(k3),num2str(g(i,4))];
str33=['(',num2str(g(i,1)),h(k1),num2str(g(i,2)),h(k2),num2str(g(i,3)),')',h(k3),num2str(g(i,4))];
str44=['(',num2str(g(i,1)),h(k1),num2str(g(i,2)),')',h(k2),'(',num2str(g(i,3)),h(k3),num2str(g(i,4)),')',];A=str2num(str11);B=str2num(str22);C=str2num(str33);D=str2num(str44);if A==24||B==24||C==24||D==24 break else end end
if A==24||B==24||C==24||D==24 break else end end
if A==24||B==24||C==24||D==24 break else end end
if A==24||B==24||C==24||D==24 break else end end
if A==24
answer=str11;elseif B==24
answer=str22;elseif C==24
answer=str33;elseif D==24
answer=str44;else
answer='无解';end
disp(['计算方法',num2str(answer)])time=toc;
disp(['计算耗时',num2str(time),'s'])
第四篇：matlab计算AHP层次分析法
用matlab解决层次分析法AHP
1、求矩阵最大特征值及特征向量 用matlab求：
输入：A=[1 1/2 2 1/4;2 1 1 1/3;1/2 1 1 1/3;4 3 3 1]
[x,y]=eig(A)得出：特征向量x=[0.2688 0.3334 0.2373 0.8720]
最大特征值λmax=4.1964
2、一致性检验
CI=（λmax-n）/(n-1)=(4.1964-4)/(4-1)=0.0655 CR=CI/RI=0.0655/0.9=0.0727
(注：维数为4时，RI=0.9)CR=0.0727<0.1,矩阵一致性通过检验
3、对最大特征值进行归一化处理，即可得到各指标权重（归一化：分项/分项之和）W=[0.157 0.195 0.139 0.510]
第五篇：用matlab电力系统潮流计算
题目：潮流计算与matlab
教学单位 电气信息学院
姓 名 学 号
年 级
专 业 电气工程及其自动化
指导教师
职 称 副教授
摘 要
电力系统稳态分析包括潮流计算和静态安全分析。本文主要运用的事潮流计算，潮流计算是电力网络设计与运行中最基本的运算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中的各元件的电力损耗，进而求得电能损耗。本位就是运用潮流计算具体分析，并有MATLAB仿真。
关键词： 电力系统 潮流计算 MATLAB
Abstract Electric power system steady flow calculation and analysis of the static safety analysis.This paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power.The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation.Key words: Power system;Flow calculation;MATLAB simulation
目 录 任务提出与方案论证....................................................................................................................................2 2 总体设计........................................................................................................................................................3 2.1潮流计算等值电路.............................................................................................................................3 2.2建立电力系统模型.............................................................................................................................3 2.3模型的调试与运行.............................................................................................................................3 3 详细设计........................................................................................................................................................4 3.1 计算前提............................................................................................................................................4 3.2手工计算.............................................................................................................................................7 4设计图及源程序...........................................................................................................................................11 4.1MATLAB仿真.......................................................................................................................................11 4.2潮流计算源程序...............................................................................................................................11 5 总结.............................................................................................................................................................31 参考文献..........................................................................................................................................................32 任务提出与方案论证
潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况下确定系统稳态运行状态的一种基本方法,是电力系统规划和运营中不可缺少的一个重要组成部分。可以说,它是电力系统分析中最基本、最重要的计算,是系统安全、经济分析和实时控制与调度的基础。常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。是电力系统研究人员长期研究的一个课题。它既是对电力系统规划设计和运行方式的合理性、可靠性及经济性进行定量分析的依据，又是电力系统静态和暂态稳定计算的基础。
潮流计算经历了一个由手工到应用数字电子计算机的发展过程，现在的潮流算法都以计算机的应用为前提用计算机进行潮流计算主要步骤在于编制计算机程序，这是一项非常复杂的工作。对系统进行潮流分析，本文利用 MATLAB中的SimpowerSystems工具箱设计电力系统，在simulink 环境下，不仅可以仿真系统的动态过程，还可以对系统进行稳态潮流分析。
总体设计
SimpowerSystems使用Simulink环境，可以将该系统中的发电机、变压器，线路等模型联结起来，形成电力系统仿真模拟图。在加人测量模块，并对各元件的参数进行设置后，用measurement和sink中的仪器可以观察各元件的电压、电流、功率的大小。
2.1潮流计算等值电路
10MWYN,d11463MWVA15MWGGGG120MW10kVp015.7kWps73kWI%0.50Vs%10.5YN,d1116MWVA463MW“xd0.134x20.161x0.060cosN0.8510kVp011kWps50kWI%0.550Vs%10.5YN,d11210MWVA35kV32km25MW110kV80km25MW110kV70km110kVYN,d11220MWVA20MWGGG415MW”xd0.136x20.16x0.0730cosN0.8p018.6kWps89kWI%0.530MW0Vs%10.510kVp015.7kWps73kWI%0.5GG0V%10.5sYN,d11216MWVA63MWVAp044kWps121kW10kVI%0.350Vs%10.5GG35MWYN,Y,d11210MVA10kVGGG312MW150MW“xd0.128x20.154x0.0540cosN0.85225MW”xd0.128x20.157x0.05910cosN0.8x0.136x20.161x0.0750cosN0.8“d
2.2建立电力系统模型
在Simulink中按照电力系统原型选择元件进行建模。所建立的模型和建立的方法在详细设计中详述。
在电力系统模型的建立工程中主要涉及到的是：元器件的选择及其参数的设置；发电机选型；变压器选择；线路的选择；负荷模型的选择；母线选择。
2.3模型的调试与运行
建立系统模型，并设置好参数以后，就可以在Simulink环境下进行仿真运行。运行的具体结果和分析也在详细设计中详述。
30km35kV0km31p044kWps121kWI%0.350Vs%10.5p013.2kWps63kWI0%0.55V%10.5s(12)80MWVs(23)%6.55MWVs(13)%17.53 详细设计
3.1 计算前提
首先是发电机的参数计算，先对5个发电厂简化为5台发电机来计算。发电机G1：
P141560MWQ160tan(arccos0.8)45MVar发电机G2：
P2463252MWQ2252tan(arccos0.85)156MVar发电机G3：
P331236MWQ336tan(arccos0.8)27MVar发电机G4：
P415050MWQ450tan(arccos0.85)31MVar发电机G5：
P522550MWQ550tan(arccos0.8)37.5MVar
其次是变电站的参数计算，我们还是对7个变电站简化为7台变压器来计算。变压器T1：
2psVN7311023RT1101033.450232SN(1610)2Vs%VN10.51102XT1101079.406SN16103S01p0j变压器T2：（双并联）
I0%SN(0.0157j0.0800)MVA 100RT2XT221psVN18911023101031.3462322SN2(2010)21Vs%VN110.51102101031.7625 32SN22010S022(p0j变压器T3：（四并联）
I0%SN)(0.0372j0.2000)MVA100 1psVN2112111023RT3101030.0922324SN4(6310)XT31Vs%VN2110.5110210105.042 4SN463103I0%SN)(0.1760j0.8820)MVA100S034(p0j变压器T4：（双并联）
1RT11.725021 XT4XT139.70302S042S01(0.0314j0.1600)MVART4变压器T5：
RT54RT30.3680XT54XT320.168S051S03(0.0440j0.2205)MVA41633523100.386 322(1010)
变压器T6：（两个三绕组变压器并联）
RT61RT62RT631[Vs(12)%Vs(13)%Vs(23)%]10.7521Vs2%[Vs(12)%Vs(23)%Vs(13)%]0.25
21Vs3%[Vs(13)%Vs(23)%Vs(12)%]6.752Vs1%21Vs1%VNXT61106.5842SNXT62XT6321Vs2%VN100.1532SN21Vs3%VN104.134 2SNS062(P06j变压器T7：（双并联）
I0%10)(0.0264j0.1100)MVA 100 RT7XT721psVN1503523101030.3062322SN2(1010)21Vs%VN110.535210106.4312SN210103
S072(p0jI0%SN)(0.0220j0.1100)MVA100再次是传输线参数计算，5条传输线的具体计算如下。
根据教材查得r00.21/km x00.4/km b02.810S/km 6线路L1：
线路L2：
线路L3：（双回路）
线路L4：
线路L5：（双回路）RL1r0l10.21408.4XL1x0l10.44016B6L1b0l12.810401.12104S Q1L12BL1V2N0.6776MVarRL2r0l20.2113027.3XL2x0l20.413052B6L2b0l22.8101303.64104S Q1L22BL2V2N2.2022MVarR12rl1L30320.21707.35X11L32x0l320.47014 BL32b40l322.8106703.9210SQ1L32B2L3VN2.3716MVarRL4r0l40.216012.6XL4x0l40.46024BL4b0l42.8106601.68104S Q12L42BL4VN1.0164MVar
11RL5r0l50.21202.12211XL5x0l50.4204 22BL52b0l522.8106201.12104S1QL5BL3VN20.0686MVar23.2手工计算
FLR1：
P2102ST12(RT1jXT1)(3.450j74.406)(0.0285j0.6562)MVA2VN110Sa10MWST1S01jQL1(10.0442j0.1142)MVAP2Q210.044220.11422SL1(RL1jXL1)(8.4j16)(0.070j0.1334)MVAVN21102ST2P2Q2402452(RT2jXT2)(1.346j31.7625)(0.4032j9.5156)MVAVN21102FLR2SbSG120ST260j45200.4032j9.5156(39.5968j35.4844)MVAScSbSa25jQL1SL1(4.4826j35.9144)MVA：
ST3P2Q225221562(RT3jXT3)(0.092j5.042)(0.6679j36.6024)MVA22VN110PQ4.493134.1048(RjX)(27.3j52)(2.67j5.0854)MVAL2L222VN1102222Sc'(4.4931j34.1048)MVASL2FLRSdSG2Sc'120ST3S03jQL2SL2(132.9792j149.229)MVA3：
ST4P2Q262272(RT4jXT4)(1.725j39.703)(0.1091j2.5101)MVAVN21102P2Q2133.59552149.99562(RL3jXL3)(7.35j14)(24.51j46.682)MVA22VN110'Sd(133.5955j149.9956)MVASL3'SeSG3Sd3025ST4S04jQL3SL3(89.945j130.0151)MVAFLR4： ST5P2Q2502312(RT5jXT5)(0.368j20.168)(0.1052j5.7687)MVA22VN110P2Q292.74872133.99372(RL4jXL4)(12.6j24)(27.654j52.674)MVA22VN110Se'(92.7481j133.9937)MVASL4SfSG4Se'80ST5S05jQL4SL4(34.9449j107.3469)MVAFLR5: 152ST72(0.306j6.431)(0.0562j1.1812)MVA35Sh15ST7S07jQL5(15.0782j0.3422)MVA15.078220.34222SL5(2.1j4)(0.3899j0.743)MVA352SiShSL5S06jQL55(20.4945j1.1266)MVA 15237.52ST63(0.386j4.34)(0.514j5.7793)MVA35220.650520.54512ST62(0.386j0.153)(0.1345j0.0533)MVA23526.336298.73692ST61(0.386j6.584)(3.2905j56.1256)MVA352SgSfST61SG5ST62ST63Si35(25.5114j194.12)MVA计算每一个FLR的功率分布和电压分布计算如下： FLR1：
VT2PRQX401.3464531.762512.8970kVVN115 Vb115VT2102.1030kVPRQX10.04428.40.144216VL10.8489kVVb102.1030VaVbVL1101.2541kVFLR2：
功率分布：
SL2ZZZ*T3*L2*SdT3(VbVN)ZZL2****VN(0.092j5.042)(132.9792j149.229)1418.6727.392j57.042T3(4.8812j13.8097)MVAST3ZZZ*L2*L2*SdT3(VbVN)ZZL2VN(27.3j52)(132.9792j149.229)1418.6727.392j57.042T3(108.687j122.62)MVA 电压分布：
Sc1SL2SL2(4.8812j13.8097)(2.67j5.0854)(7.5512j8.7243)MVA7.551227.38.7243522.424kV102.1030VdVbVL2102.103(2.424)104.527kVVL2功率分布：
FLR3：
SL3ZZZ*T4*L3*SeT4(VG3Vd)ZZL3****VN(1.725j39.703)(89.945j130.0151)1037.9279.075j53.73T4(59.444j16.846)MVAST4ZZZ*L3*L3*SeT4(VbVN)ZZL3VN(7.35j14)(89.945j130.0151)1037.9279.075j53.73T4(31.811j60.1256)MVA 电压分布：
Se1SL3SL3(59.444j19.846)(24.51j46.682)(83.954j26.836)MVA83.9547.3526.836149.404kV105.5643VeVdVL396.16kVVL3功率分布：
FLR4：
SL4ZZZ*T5*L4*SfT5(VG3Vd)ZZL4****VN=(0.368j20.168)(34.9449j107.3469)1037.92712.968j44.168T5(20.843j19.689)MVAST4ZZZ*L4*L4*SfT5(VG3Vd)ZZL4VN=(12.6j24)(34.9449j107.3469)1037.92712.968j44.168T5(1.398j44.389)MVA 电压分布：
Se1SL3SL3(59.444j16.846)(24.51j46.682)(83.954j63.528)MVA83.95412.663.5282424.464kV105.5643VeVdVL381.10kVVL4FLR5: 这里我们先将f点和发电机G5当做电源，经过ZT61和ZT63构成两端供电网络以g点作为运算负荷进行计算。ST6ST4(0.386j4.134)(20.2656j70.9293)(22.093837)35(3.900j25.1175)MVA0.772j10.718(0.386j6.584)(20.2656j70.9293)(22.093837)35(16.5061j91.7905)MVA0.772j10.718电压分布：
ST631ST63ST63(16.6421j97.5698)MVA16.64210.38697.56984.13410.9186kV37Vg37VT6326.0814VVT6320.26560.38670.9293(0.153)0.1162kV
26.0814ViVgVT6226.1976VT6220.49452.11.126641.815kV26.1976VhViVL524.3826VL5
4设计图及源程序
4.1MATLAB仿真
相关的原始数据输入格式如下：
1、B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写。
2、对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗置于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为含变压器支路的变压器的变比，第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器。
3、B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
4、X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地参数。
4.2潮流计算源程序
%本程序的功能是用牛顿——拉夫逊法进行11节点潮流计算 clear;n=11;%input('请输入节点数：n=');nl=11;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;
0.05620+0.13289i
0.05382i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2];% B1矩阵：
1、支路首端号；
2、末端号；
3、支路阻抗；
4、支路对地电纳 %
5、支路的变比；
6、支路首端处于K侧为1，1侧为0 % B2矩阵：
1、该节点发电机功率；
2、该节点负荷功率；
3、节点电压初始值 %
4、PV节点电压V的给定值；
5、节点所接的无功补偿设备的容量 %
6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %
3为PV节点；
%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路数
if B1(i,6)==0
%左节点处于1侧
p=B1(i,1);q=B1(i,2);
else
%左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元
Y(q,p)=Y(p,q);
%非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
end %求导纳矩阵
disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出导纳阵的实部和虚部
for i=1:n
%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
end for i=1:n
%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量
end %===================== P=real(S);Q=imag(S);
%分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0
% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2;
%电压模平方
%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========
if B2(i,6)~=3
%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========
DP=P(i)-P1;
% PV节点有功误差
DV=V(i)^2-V2;
% PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
end
end
end
end
end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================
for k=3:N0
% N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U
for k2=k1:N1
% 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U
J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1;
% 对角元规格化K行K列对角元素赋1
%==================== 回代运算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行从第三行开始到当前行的前一行k4行消去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end
%用当前行K2列元素减去当前行k列元素乘以第k行K14 列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代运算
==================================
for k3=k1:N0
% 从k+1行到2*n最后一行
for k2=k1:N1
% 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
else
%是第三行k=3
%====================== 第三行k=3的前代运算 ========================
for k3=k1:N0
%从第四行到2n行（最后一行）
for k2=k1:N1
%从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end
%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0;%当前行第3列元素已消为0
end
end
end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
end
%------修改节点电压-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2(a)=IT2;
%不满足要求的节点数
ICT1=ICT1+1;
%迭代次数 end %用高斯消去法解”w=-J*V“
disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%计算各节点电压的角度
E(k)=e(k)+f(k)*j;
%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);
%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);
%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);
%显示各节点的电压相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);
%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);
%显示各节点的注入功率
disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
%本程序的功能是用牛顿——拉夫逊法进行10节点潮流计算 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 clear;n=10;%input('请输入节点数：n=');nl=10;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2];% B1矩阵：
1、支路首端号；
2、末端号；
3、支路阻抗；
4、支路对地电纳 %
5、支路的变比；
6、支路首端处于K侧为1，1侧为0 % B2矩阵：
1、该节点发电机功率；
2、该节点负荷功率；
3、节点电压初始值 %
4、PV节点电压V的给定值；
5、节点所接的无功补偿设备的容量 %
6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %
3为PV节点；
%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路数
if B1(i,6)==0
%左节点处于1侧
p=B1(i,1);q=B1(i,2);
else
%左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元
Y(q,p)=Y(p,q);
%非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
end %求导纳矩阵
disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出导纳阵的实部和虚部
for i=1:n
%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
end for i=1:n
%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量
end %===================== P=real(S);Q=imag(S);
%分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0
% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2;
%电压模平方
%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========
if B2(i,6)~=3
%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素
===========
DP=P(i)-P1;
% PV节点有功误差
DV=V(i)^2-V2;
% PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
end
end
end
end
end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================
for k=3:N0
% N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U
for k2=k1:N1
% 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U
J(k,k2)=J(k,k2)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1;
% 对角元规格化K行K列对角元素赋1
%==================== 回代运算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行从第三行开始到当前行的前一行k4行消
去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代运算
==================================
for k3=k1:N0
% 从k+1行到2*n最后一行
for k2=k1:N1
% 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
else
%是第三行k=3
%====================== 第三行k=3的前代运算 ========================
for k3=k1:N0
%从第四行到2n行（最后一行）
for k2=k1:N1
%从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end
%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0;%当前行第3列元素已消为0
end
end
end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
end
%------修改节点电压-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2(a)=IT2;
%不满足要求的节点数
ICT1=ICT1+1;
%迭代次数 end %用高斯消去法解”w=-J*V“ disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%计算各节点电压的角度
E(k)=e(k)+f(k)*j;
%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);
%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);
%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);
%显示各节点的电压相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);
%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);
%显示各节点的注入功率
disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
%本程序的功能是用牛顿——拉夫逊法进行12节点潮流计算 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 clear;n=12;%input('请输入节点数：n=');nl=12;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;
0.05620+0.13289i
0.05382i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0
0
0
2];% B1矩阵：
1、支路首端号；
2、末端号；
3、支路阻抗；
4、支路对地电纳 %
5、支路的变比；
6、支路首端处于K侧为1，1侧为0 % B2矩阵：
1、该节点发电机功率；
2、该节点负荷功率；
3、节点电压初始值 %
4、PV节点电压V的给定值；
5、节点所接的无功补偿设备的容量 %
6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %
3为PV节点；
%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路数
if B1(i,6)==0
%左节点处于1侧
p=B1(i,1);q=B1(i,2);
else
%左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元
Y(q,p)=Y(p,q);
%非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
end %求导纳矩阵
disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出导纳阵的实部和虚部
for i=1:n
%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
end for i=1:n
%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量
end %===================== P=real(S);Q=imag(S);
%分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0
% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2;
%电压模平方
%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========
if B2(i,6)~=3
%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========
DP=P(i)-P1;
% PV节点有功误差
DV=V(i)^2-V2;
% PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
end
end
end
end
end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================
for k=3:N0
% N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U
for k2=k1:N1
% 从k+1列的Jacobi元素到扩展列的△P、△Q
或 △U
J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1;
% 对角元规格化K行K列对角元素赋1
%==================== 回代运算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行从第三行开始到当前行的前一行k4行消去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代运算
==================================
for k3=k1:N0
% 从k+1行到2*n最后一行
for k2=k1:N1
% 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
else
%是第三行k=3
%====================== 第三行k=3的前代运算 ========================
for k3=k1:N0
%从第四行到2n行（最后一行）
for k2=k1:N1
%从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end
%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0;%当前行第3列元素已消为0
end
end
end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
end
%------修改节点电压-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2(a)=IT2;
%不满足要求的节点数
ICT1=ICT1+1;
%迭代次数 end %用高斯消去法解”w=-J*V" disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%计算各节点电压的角度
E(k)=e(k)+f(k)*j;
%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);
%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);
%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);
%显示各节点的电压相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);
%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);
%显示各节点的注入功率
disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
如果源程序的运行结果需要作图可用下面的程序
figure(1);subplot(1,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(1,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('节点号');ylabel('节点注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('节点号');ylabel('节点注入无功');grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路号');ylabel('支路首端注入有功');grid on;subplot(2,2,4);bar(Q1);xlabel('支路号');ylabel('支路首端注入无功');grid on;
总结
通过本次课程设计让我有复习了一次潮流计算的相关知识，跟家清晰了什么事潮流计算以及潮流计算的在电力系统的重要性。电力系统的稳定运行状况即是正常运行状况，是指电力系统在稳定运行条件下电压、功率的分布，也称为潮流分布。电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。同时潮流分布也是电力系统规划设计的一项基础工作。
整个计算过程的模型建立并不是十分复杂，但计算过程十分繁琐、计算量相当的大，而且由于枝节太多很容易算错。不过在计算潮流计算的过程中却对以往学过的电力系统分析的相关知识进行了一次较为深入的复习。而且整个计算对计算量的要求很大，锻炼了我们的计算能力。而且对细节的把握也得到了锻炼，做题的精细程度得到了提高。
参考文献
[1]何仰赞, 温增银《电力系统分析》（第三版）［M].华中科技大学，2002 [2]http://baike.baidu.com/view/627420.[3]王守相，刘玉田 电力系统潮流计算研究现状--《山东电力技术》1996年05期
[4]何仰赞，温增银.电力系统分析（上册）第三版［M].湖北：华中科技大学出版社，2002 [5] 刘同娟.MATLAB在电路分析中的应用.电气电子教学学报.2002 [6] 西安交通大学等.电力系统计算[M].北京：水利电力出版社，1993.12 [7] 李光琦.电力系统暂态分析[M].北京： 水利电力出版社，2002.5 [8]何仰赞，温增银.电力系统分析（下册）第三版［M].湖北：华中科技大学出版社，2002 [9]韦化，李滨，杭乃善，等．大规模水一火电力系统最优潮流的现代内点算法实现[J]．中国电机工程学报，2003．23(6)：13一l8．
[10]Chen Luo—nan，Suzuki Hideki，Katou Kazuo．Mean fieldtheory for optimal power flow[J]．IEEE Transactions OilPower Systems，1997，12(4)：1481·1486
点击咨询 