电力系统三相短路计算的MATLAB代码.doc
MATLAB在电力系统三相短路故障分析中的应用_高正中
![MATLAB在电力系统三相短路故障分析中的应用_高正中](https://img.taocdn.com/s3/m/d9643ec85fbfc77da269b19f.png)
MATLAB 在电力系统三相短路故障分析中的应用
图 1 无穷大功率电源供电系统三相短路故障的仿真图
1.2 参数设置
数,按照仿真设置,当模型中所设故障点发生三相短路
(1)变压器 T 采用如图 1 所示的模型,根据给定的数
据,计算折算到 110kV 侧的参数如下:
变压器的电阻为:
RT=
ΔPSUN2 SN2
CLC n u m b e r: TP391.9;TM743
Do cu m e n t co d e : A
Article ID:1003- 0107(2013)10 - 0017- 03
0 引言 MATLAB 因为绘图功能强大和计算能力强,配以友
好的动态仿真环境,主要用于数值计算及可视化图形处 理,其优越的开发性、数据仿真分析高效的优点越来越 成为从事电力网络系统学习和研究的重要仿真工具。作 为一款优秀的综合性应用软件,利用其提供的 Simulink 集成环境,可以方便地对电力系统进行模型的搭建和仿 真[1]。Simulink 提供了充足的子模块库,我们可以根据相 应模型搭建的需要,从各个子模块库中选用合适的模 块。Simulink 中提供了各种基本模块,它们根据其主要 应用领域和实现功能进行了分类化管理,给用户查找使 用提供了便利。模块库的数量取决于用户安装,在电力 系统仿真中,标准 Simulink 模块库和电力系统模块库是 必不可少的。本文将通过三相短路实例具体讲解其模块 结构及应用。
合保护系统研究[J].煤炭科学技术,2010,38(12):80- 83. [4]马立国,公茂法,夏文华,等.基于 DSP 的矿井下电网漏
电保护配置方法的研究[J].煤矿机械,2012,33(5):206208. [5]孙勇.煤矿井下电网漏电保护系统设计[J].煤炭科学技 术,2012,40(5):81- 85.
matlab实验 电力系统短路分析
![matlab实验 电力系统短路分析](https://img.taocdn.com/s3/m/d2096e83aaea998fcd220e73.png)
页脚内容- 9 -实验二 短路电流计算程序的实现一、三相短路电流计算程序计算短路电流周期分量,如I ''(I ')时,实际上就是求解交流电路的稳态电流,其数学模型也就是网络的线性代数方程,一般选用节点电压方程。
方程的系数矩阵是对称的。
在短路电流计算中变化的量往往是方程的常数项,需要多次求解线性方程组。
1.等值网络图2-1给出了不计负荷情况下计算短路电流I ''的等值网络。
在图2-1(a )中G 代表发电机端电压节点,发电机等值电势和电抗分别为E '' 和dx '',D 表示负荷节点,f 点为直接短路点。
应用叠加原理如图2-1所示。
正常运行方式为空载运行,网络中各点电压均为1;在故障分量网络中。
只需作故障分量的计算。
由图2-1的故障分量网络可见,这个网络与潮流计算的网络的差别在于发电机节点上多接了对地电抗dx ''。
当然如果短路计算中可以忽略线路电阻和电纳,而且不计变压器的实际变比,则短路计算网络较潮流计算网络简化,而且网络本身是纯感性的。
1E ''x1E '' x 1E '' x 1-=1=图2-1 在不计负荷情况下计算短路电流I ″的等值电路2. 用节电阻抗矩阵计算短路电流如果已经形成了故障分量网络的节点阻抗矩阵,则矩阵中的对角元素就是网络从f 点看进去的等值阻页脚内容- 10 -抗,又称为f 点的自阻抗。
fi Z 为f 点与i 点的互阻抗,均用大写Z 表示。
由节点方程中的第f 个方程:n fn f ff f f I Z I Z I Z U ++++=11。
ff Z 为其它节电电流为零时,节点f 的电压和电流之比,即网络对f 点的等值阻抗。
根据故障分量网络,直接应用戴维南定理可求得直接短路电流(由故障点流出)为fff ffz Z U I +=0(2-1)式中,f z 为接地阻抗;0f U 为f 点短路前的电压。
复杂电网三相短路计算的MATLAB仿真
![复杂电网三相短路计算的MATLAB仿真](https://img.taocdn.com/s3/m/80a996e6910ef12d2af9e741.png)
请输入支路起始结点号(由1开始)s=4
请输入支路阻抗标么值,形如“0.3i”0.03i
阻抗矩阵:
0.0000 + 0.1500i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
Y=
则相应的导纳型节点方程为;
则相应的导纳型节点方程为:
则相应的导纳型节点方程为:
则相应的导纳型节点方程为:
则相应的导纳型节点方程为:
则相应的导纳型节点方程为:
则相应的导纳型节点方程为:
因子表的标准形式为:
对U求解的回代运算:
求解5节点短路时的短路故障电流和线路56,57,58的支路故障电流:
图1-2短路电流的计算过程
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0750i
请确定增加的支路种类1—接地树支;2—不接地树支;3—不接地链支;4—接地链支:type=3
请输入支路起始结点号(由1开始)s=2
请输入支路终止结点号(由1开始,比起始结点号大)e=3
请输入支路阻抗标么值,形如“0.3i”0.1i
第一阶段:分析原始资料、明确设计内容、制定详细计划;
学习阻抗矩阵及导纳矩阵,学习利用支路追加法形成阻抗矩阵的方法,学习三角分解法求解导纳矩阵的方法;
第二阶段分析目标电网,生成计算模型,确定设计方法;
第三阶段:根据前两阶段的成果,生成复杂电网的计算矩阵,完成矩阵的三角分解过程,并完成手工计算过程;
电力系统短路故障分析的MATLAB辅助程序设计短路计算程序
![电力系统短路故障分析的MATLAB辅助程序设计短路计算程序](https://img.taocdn.com/s3/m/671ca63530b765ce0508763231126edb6f1a7608.png)
电力系统短路故障分析的MATLAB辅助程序设计短路计算程序电力系统短路故障分析是电力系统设计和运行过程中非常重要的一环。
短路故障会导致电力系统各个部分的电压、电流和功率的突然变化,对设备的保护和稳定运行产生不利影响甚至引起事故。
因此,进行短路计算和故障分析非常必要。
MATLAB是一种功能强大的数值计算和数据可视化工具,对于电力系统短路计算和故障分析也可以发挥重要的作用。
下面将介绍如何使用MATLAB设计一个简单的电力系统短路计算程序。
首先,我们需要建立一个电力系统的模型。
电力系统可以用图模型表示,其中节点表示发电机、变压器、负荷等设备,边表示导线、变压器等电力连接。
我们可以使用MATLAB中的图模型工具箱创建电力系统模型,并且设置各个节点和边的属性,例如电压、电流、阻抗等。
然后,我们需要编写短路计算程序。
短路计算可以分为对称故障和不对称故障两种情况。
对称故障是指短路故障发生在电力系统的正常运行条件下,例如三相短路。
不对称故障是指短路故障发生在电力系统的不正常运行条件下,例如单相接地短路。
对于对称故障,我们可以使用节点电流法进行计算。
首先,应用基尔霍夫电流定律,根据电压和阻抗计算电流。
然后,根据节点电流方程和电流方程计算电流分布。
最后,根据电流分布计算短路电流和故障点的电压。
对于不对称故障,我们可以使用仿真方法进行计算。
首先,需要设置故障位置和故障类型,例如A相到地短路。
然后,根据故障位置和类型修改节点和边的参数,例如将故障位置的阻抗设置为零。
最后,使用数值方法求解电力系统的动态响应,得到短路电流和故障点的电压。
在MATLAB中,可以使用矩阵运算和数值求解函数实现短路计算。
例如,可以使用矩阵乘法和矩阵求逆函数计算节点电流和电流分布。
可以使用ODE求解器求解动态响应方程。
可以使用MATLAB的绘图函数绘制电力系统的电流分布和故障点的电压。
总结起来,电力系统短路故障分析的MATLAB辅助程序设计涉及建立电力系统模型、编写短路计算程序并使用MATLAB的数值计算和数据可视化工具进行计算和分析。
三相短路和单相接地短路MATLAB
![三相短路和单相接地短路MATLAB](https://img.taocdn.com/s3/m/11d17785f705cc175527097f.png)
电力系统三相短路和单相接地短路
实验目的:
1.理解掌握短路的类型及在短路故障后的影响。
2.运用 MATLAB 的电力系统工具箱对三相短路和单相接地短路进行建模,
并分别观察分析其电压、电流的波形,并得出结论。
实验内容:
在二相电力系统中,大多数故障都是由于短路故障引起的,在发生短路故障的情况下,电力系统从一种状态剧烈变化到另一种状态,产生复杂的暂态现象。
在三相系统中,可能发生的短路有:三相短路、两相短路、一相短路接地和单相接地短路。
建立如图 3.1 理想情况下小型电力系统的模型。
图 3.1
参数设置:
电源模块图 3.2:
图 3.2三相电压电流测量模块图 3.3
图 3.3输电线路图 3.4:
图 3.4输电线路模块图 3.5:
图 3.5
万用表模块图 3.6:
(电压)
(电流)
图 3.6万用表 1 图 3.7
图 3.7三相序分量分析图 3.8:
图 3.8算法模块图 3.9:
图 3.9
单项接地:
故障模块 3.10:
图 3.10
对于测量模块和三相模块相同。
matlab电力系统快速解耦法潮流计算及短路计算程序
![matlab电力系统快速解耦法潮流计算及短路计算程序](https://img.taocdn.com/s3/m/3380bf77f56527d3240c844769eae009581ba23e.png)
电力系统快速解耦法潮流分析及短路计算一.程序设计的基本思想:(1)由于电力系统潮流分析中要利用到矩阵运算,复数运算,故采用matlab编程。
采用文件输入,将系统的各个参数以文件的形式输入,便于程序的通用化。
(2)本程序共有两个输入文件,分别为线路参数的文件,和已知的节点状态文件(PQ)(3)为了使程序不仅仅局限于计算9节点网络,在形成节点导纳的函数Yn()中,利用循环,找出线路首节点中的最大编号数,自动确定节点导纳矩阵的维数。
故对于任意n节点网络,均可以计算出节点导纳矩阵(4)在(3)的前提下,为了使程序支持系统增加节点,增加负荷等造成的PQ参数改变,或者PQ表的加长。
对程序做了如下优化。
首先,程序执行的基础是PQ表中平衡节点在第一行,接下来是PV节点,最后是PQ节点,如果系统添加节点,或者删除节点,均在PQ表的末端操作,会造成PQ表的顺序不是平衡节点、PV节点、PQ节点的顺序。
故引入了seqencing()函数,其作用就是不论输入的PQ表是什么顺序,在程序读入后均按平衡-PV-PQ的顺序排列。
其次,顺序打乱的PQ表必须与支路参数表对应,故在Yn()函数中加入了两段循环体,使之对应(见相应函数体注释)(5)在满足了上述4个条件后,程序便可以通用化了。
当然,由于水平有限,且程序未能由大量数据测试,故缺陷在所难免,这里仅是做了通用化的尝试。
在本文最末附加了该程序通用化的实例。
二、潮流计算框图三.定义相应的函数1.形成节点导纳矩阵的函数Yn()function Y=Yn(x,y)%定义一名为Yn的函数,其功能是自动识别输入表中节点的个数,形成相应的节点导纳矩阵[fid,message]=fopen(x,'r') ; %从x文件中读入支路参数if fid==-1; %判断文件是否正确打开error(message);end;[HeadPoint,HeadNumber, EndPoint,EndNumber,R,X,B,k]=textread(x,'%s %d %s %d %f %f %f %f'); %将读入的参数处理为以列为向量的数组fclose(fid);%关闭文件L=length(HeadNumber); %确定输入表的行数[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);%调用seqencing函数,引入y文件中的PQ参数A=PointNumber;for i=1:L; %通过以下两循环体,实现PQ参数与支路参数的编号对应for j=1:L;if HeadNumber(i)==j;HeadNumber(i)=A(j);break;end;end;end;for i=1:L;for j=1:L;if EndNumber(i)==j;EndNumber(i)=A(j);break;end;end;end;Y=zeros(L,L); %根据txt文件中数据表的长度建立空的节点导纳矩阵for i=1:Lm=HeadNumber(i);n=EndNumber(i);if k(i)==0; %判断是否何种元件,为输电线元件if n~=0;Y(m,m)=Y(m,m)+1j*B(i)+1/(R(i)+1j*X(i));Y(n,n)=Y(n,n)+1j*B(i)+1/(R(i)+1j*X(i));Y(m,n)=Y(m,n)-1/(R(i)+1j*X(i));Y(n,m)=Y(n,m)-1/(R(i)+1j*X(i));elseY(m,m)=Y(m,m)+R(i)+1j*X(i);end;else %为变压器元件if n~=0;Y(m,m)=Y(m,m)+1/(R(i)+1j*X(i));Y(m,n)=Y(m,n)-1/(k(i)*(R(i)+1j*X(i)));Y(n,n)=Y(n,n)+1/(k(i)*k(i)*(R(i)+1j*X(i)));Y(n,m)=Y(n,m)-1/(k(i)*(R(i)+1j*X(i)));elseY(m,m)=Y(m,m)+R(i)+1j*X(i);end;end;end;maxm=HeadNumber(1);%通过下面两个循环体,确定输入表中节点编号的最大值,及为节点导纳矩阵的维数for i=1:L;if maxm<=HeadNumber(i);maxm=HeadNumber(i);end;end;maxn=EndNumber(1);for i=1:L;if maxn<=EndNumber(i);maxn=EndNumber(i);end;end;Y=Y(1:max(maxm,maxn),1:max(maxm,maxn));%形成导纳矩阵2.对不满足要求的PQ参数表进行排序的函数seqencing()function [Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y)%定义名为seqencing的函数,其功能是在系统添加节点,或输入的PQ参数的顺序不满足要求时,对PQ参数表进行重新排序,保证平衡节点放在第一行,接下来是PV节点,最后是PQ节点[fid,message]=fopen(y,'r'); %从y文件中读入PQ参数if fid==-1; %判断文件是否正确打开error(message);end;[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=textread(y,'%f %f %f %f %f %f');fclose(fid);L=length(PointNumber);%通过以下两个循环体,完成对PQ输入表的重新排序,其思想是,在PQ参数之前加入一列Pointstyle用于标识节点类型,平衡节点为0,PV节点为1,PQ节点为2,以Pointstyle列为基准进行排序for i=1:L;for j=1:L-i;if Pointstyle(j)>Pointstyle(j+1);t=Pointstyle(j+1);Pointstyle(j+1)=Pointstyle(j);Pointstyle(j)=t;t=PointNumber(j+1);PointNumber(j+1)=PointNumber(j);PointNumber(j)=t;t=Ps(j+1);Ps(j+1)=Ps(j);Ps(j)=t;t=Qs(j+1);Qs(j+1)=Qs(j);Qs(j)=t;t=Uk(j+1);Uk(j+1)=Uk(j);Uk(j)=t;end;end;end;3、形成解耦算法B’矩阵的函数 formB1()function B1=formB1(x,y)%定义名为B1的函数形成解耦算法中的B’矩阵,得到的B’矩阵用B1表示[fid,message]=fopen(x,'r') ; %从x文件中读入支路参数if fid==-1; %判断文件是否正确打开error(message);end;[HeadPoint,HeadNumber,EndPoint,EndNumber,R,X,B,k]=textread(x,'%s %d %s %d %f %f %f %f'); %将读入的参数处理为以列为向量的数组fclose(fid);L=length(HeadNumber);[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);%调用seqencing函数,引入y文件中的PQ参数A=PointNumber;%通过以下两循环体,实现PQ参数与支路参数的编号对应for i=1:L;for j=1:L;if HeadNumber(i)==j;HeadNumber(i)=A(j);break;end;end;end;for i=1:L;for j=1:L;if EndNumber(i)==j;EndNumber(i)=A(j);break;end;end;end;B1=zeros(L,L);for i=1:L %以行为单位,通过循环,用支路参数对B1进行修改,形成B’矩阵 m=HeadNumber(i);n=EndNumber(i);B1(m,m)=B1(m,m)-1/X(i);B1(n,n)=B1(n,n)-1/X(i);B1(m,n)=B1(m,n)+1/X(i);B1(n,m)=B1(n,m)+1/X(i);endmaxm=HeadNumber(1);for i=1:L;if maxm<=HeadNumber(i);maxm=HeadNumber(i);end;end;maxn=EndNumber(1);for i=1:L;if maxn<=EndNumber(i);maxn=EndNumber(i);end;end;B1=B1(2:max(maxm,maxn),2:max(maxm,maxn)); %形成B’矩阵4、形成解耦算法B’’矩阵的函数 formB11()function B11=formB11(x,y)%定义名为B11的函数形成解耦算法中B’'矩阵,用B11表示从x文件中读入支路参数确定Y,从y文件中读入PQ参数确定B11的维数,即除去平衡节点和pv节点,此处要求PQ参数录入时,将平衡节点和PQ节点放在前排,这一要求在Yn函数中通过seqencing函数已经满足Y=Yn(x,y);B=imag(Y);[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);i=1;j=1;while Pointstyle(i)<=1;i=i+1;j=j+1;end;B11=B(j:end,j:end); %形成B’’矩阵5、计算正常情况下系统节点电压的函数 powerflow()function [U0,O0]=powerflow(x,y)%定义名为powerflow的函数,利用快速解耦算法来计算正常情况下系统内各个节点的电压和相角[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);% 调用seqencing函数对PQ参数表进行排序Y=Yn(x,y); %形成节点导纳矩阵,Yn为n维B1=formB1(x,y); %形成解耦算法中的B矩阵,B1为n-1维B11=formB11(x,y); %形成解耦算法中的B'矩阵,B'为m维G=real(Y); %取Y的实部B=imag(Y); %取Y的虚部U0=Uk;O0=Ok;L=length(PointNumber);P=zeros(L,1);Q=zeros(L,1);dP=zeros(L,1);dQ=zeros(L,1);number=1;i=1;k=1;while Pointstyle(i)<=1; %通过k值确定系统中PQ节点的个数i=i+1;k=k+1;end;while number<100 %定义迭代次数上限为100次for i=2:L;sum1=0;for j=1:L;sum1=sum1+U0(j)*(G(i,j)*cos(O0(i)-O0(j))+B(i,j)*sin(O0(i)-O0(j))); %潮流方程,n-1维end;dP(i)=Ps(i)-U0(i)*sum1;endfor i=k:L;sum2=0;for j=1:L;sum2=sum2+U0(j)*(G(i,j)*sin(O0(i)-O0(j))-B(i,j)*cos(O0(i)-O0(j))); %潮流方程,m 维end;dQ(i)=Qs(i)-U0(i)*sum2;enddP1=dP(2:L)./U0(2:L);dQ1=dQ(k:L)./U0(k:L);a=max(norm(dP1,inf));b=max(norm(dQ1,inf));if max(a,b)<0.00001 %判断是否收敛break;disp(‘迭代’)disp(k);disp(‘次后收敛’);else %如不收敛,dO=-inv(B1)*dP1; %dO为n-1维dU=-inv(B11)*dQ1; %dU为m维zero1=zeros(k-1,1);zero2=[0];DU=[zero1;dU];DO=[zero2;dO];U0=U0+DU;O0=O0+DO;number=number+1;end;if number==100;disp('迭代100次后不收敛,迭代结束');end;end;四.对相应的系统进行潮流分析和短路计算定义完上述函数之后,可直接调用函数形成导纳矩阵,计算正常情况下的节点电压,进行短路计算计算短路电流,短路后各个节点电压以及支路潮流分布。
电力系统三相短路的编程计算
![电力系统三相短路的编程计算](https://img.taocdn.com/s3/m/1bc40426793e0912a21614791711cc7930b7787e.png)
电力系统三相短路的编程计算简介:电力系统三相短路是指电力系统发生故障时,三相电流短路连接,导致系统电流急剧增加,可能引发设备烧毁和电力系统瘫痪的危险。
因此,进行电力系统短路计算是非常重要的。
本文将介绍如何利用编程语言进行电力系统三相短路计算。
1.三相短路计算原理三相短路计算是指在给定电网拓扑、负荷数据以及故障点位置的情况下,通过计算得到短路电流的数值和相位。
计算过程主要包括以下几个步骤:1.确定故障位置:根据电网的拓扑结构确定故障点的位置。
2.分解复杂电网为简化电路:将复杂的电网拓扑结构分解为简化的电路模型。
3.计算电流分布:通过求解电流方程,得到各支路的电流分布。
4.计算短路电流:根据电流分布和故障点位置,计算短路电流大小和相位。
2.编程实现步骤2.1数据输入首先,需要从用户获取输入数据,包括电网拓扑结构、负荷数据和故障点位置。
可以通过编程语言提供的输入函数或者读取数据文件的方式获取数据。
2.2电网拓扑解析将输入的电网拓扑结构进行解析,构建电路模型。
可以采用图论算法进行解析,比如深度优先遍历或广度优先遍历,将电网拓扑结构转化为图的邻接矩阵或邻接表。
2.3电流分布计算利用电路分析的方法,根据电路模型和负荷数据,计算每个节点的电压值。
可以采用基尔霍夫电压法或者毕奥萨法进行计算。
2.4短路电流计算根据故障点位置和电流分布,计算各个支路的电流值和相位。
可以采用基尔霍夫电流法来计算短路电流。
2.5结果输出将计算得出的短路电流值和相位进行输出,可以利用编程语言提供的输出函数或者将结果保存到文件中。
3.编程计算实例以下是一个使用Python语言实现电力系统三相短路计算的示例代码:```pythonimport numpy as np#输入电网拓扑结构、负荷数据和故障点位置等数据topology = np.array([[0, 1, 1, 0],[1,0,1,1],[1,1,0,1],[0,1,1,0]])loads = np.array([100, 200, 150, 100])fault_location = 2#电压计算voltages = np.linalg.solve(topology, loads)#短路电流计算short_circuit_current = loads[fault_location] /np.sum(voltages)#输出结果print("短路电流大小:", short_circuit_current)```以上代码使用numpy库进行矩阵计算,首先输入电网拓扑结构、负荷数据和故障点位置等数据,然后利用numpy库提供的线性方程求解函数linalg.solve(计算电压值。
短路故障计算
![短路故障计算](https://img.taocdn.com/s3/m/342ef8334b35eefdc9d33310.png)
一、三相短路故障计算MATLAB程序n=4; %独立节点数nl=3; %支路数B1=[1 3 0.51i 0 1 0;2 3 0.59i 0 1 0;3 4 1.43i 0 1 0]; %线路参数形成的矩阵X=[1 0.2i;2 4i;3 0;4 0]; %对地阻抗形成的矩阵V0=[1;1;1;1]; %各节点的初电压标幺值形成的列矩阵D=[4 0]; %短路点阻抗组成的矩阵NF=1; %短路点数目B=[0;0;0;1]; %常数项矩阵(短路点对应为1,其余节点对应为0)Y=zeros(n); %初始化节点导纳矩阵for i=1:n %初始化节点是否非直接接地if X(i,2)~=0; %判断i节点是否非直接接地p=X(i,1);Y(p,p)=1./X(i,2); %算自导中接地阻抗的导纳endend%形成节点导纳矩阵for i=1:nlif B1(i,6)==0 %判断变压器的变比是否在低压侧(it=0说明i侧为低压侧;it=1说明i侧为高压侧)p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./B1(i,3)*B1(i,5); %节点互导纳Y(q,p)=Y(p,q); %节点互导纳Y(p,p)=Y(p,p)+1./B1(i,3)*B1(i,5)^2 +B1(i,4)./2; %节点导纳=自导纳+X修正量Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)./2;enddisp('导纳矩阵');disp(Y); %输出导纳矩阵A=Y;[n,m]=size(A);%解线性方程组,形成矩阵因子表Afor i=1:nA(i,i)=1./A(i,i);for j=i+1:nA(i,j)=A(i,j)*A(i,i); %矩阵规格化endfor k=i+1:nfor j=i+1:nA(k,j)=A(k,j)-A(k,i)*A(i,j);endendenddisp('矩阵A的因子表为');disp(A);B=input('请输入常数项矩阵(短路点对应为1,其余节点对应为0),:B='); %利用因子表对常数项进行的前代过程(修正矩阵B)for i=1:nB(i)=B(i)*A(i,i);for j=i+1:nB(j)=B(j)-A(j,i)*B(i);endenddisp('利用因子表对常数项进行会带的结果为: B=') ;disp(B)%利用因子表的上三角回带过程for i=n-1:-1:1for j=i+1:-1:2B(j-1)=B(j-1)-A(j-1,i+1)*B(i+1);endenddisp('在因子表的基础上求解线性方程组的解为:X=');disp(B); %解出各节点电压即得到Zik,此时B为列向量V0=input('请输入由各节点的初始电压标幺值形成的列矩阵:V0=');D=input('请输入由短路号、短路点阻抗组成的矩阵:D=');NF=input('请输入短路点的数目:NF=');Z=zeros(n);V=zeros(n);I=zeros(nl);for k=1:NF %求各短路点的电流标幺值disp('短路点与其他各节点的互阻抗');for i=1:nZ(i,D(k,1))=B(i,1);disp(Z(i,D(k,1)));endI(D(k,1),D(k,1))=V0(D(k,1),1)./(Z(D(k,1),D(k,1))+D(k,2));ft=num2str(D(k,1));ts1=('点短路时');ts2=('电流的标幺值If=');dn=strcat(ft,ts1,ts2);disp(dn);disp(I(D(k,1),D(k,1)));for i=1:nV(i,i)=V0(i,1)-I(D(k,1),D(k,1))*Z(i,D(k,1)); %求各节点的电压标幺值电压end%求非接地支路的短路电流标幺值for i=1:nlif B1(i,6)==0 %判断该支路是否为含有变压器支路且变比在低压侧k=B1(i,5);elsek=1./B1(i,5);endp=B1(i,1);q=B1(i,2);I(i,i)=(V(p,p)-V(q,q)./k)./B1(i,3); %各支路短路电流=支路残压/支路阻抗(标幺值)enddisp('各节点的电压标幺值U为(节点号从小到大排):');for i=1:ndisp(V(i,i));enddisp('各接地支路短路电流的标幺值I为(顺序同您输入B时一样');for i=1:nif X(i,2)~=0; %判断i节点是否为直接接地e=0;b=X(i,2);Ii0=(e-V(i,i))./b;disp(Ii0);endenddisp('各非接地支路短路电流的标幺值I为(顺序同您输入B时一样):');for i=1:nldisp(I(i,i));endend二、程序输出结果导纳矩阵0 - 6.9608i 0 0 + 1.9608i 00 0 - 1.9449i 0 + 1.6949i 00 + 1.9608i 0 + 1.6949i 0 - 4.3550i 0 + 0.6993i0 0 0 + 0.6993i 0 - 0.6993i矩阵A的因子表为0 + 0.1437i 0 -0.2817 00 0 + 0.5142i -0.8715 00 + 1.9608i 0 + 1.6949i 0 + 0.4300i -0.30070 0 0 + 0.6993i 0 + 2.0449i请输入常数项矩阵(短路点对应为1,其余节点对应为0),:B=[0;0;0;1]利用因子表对常数项进行会带的结果为: B=0 + 2.0449i在因子表的基础上求解线性方程组的解为:X=0 + 0.1732i0 + 0.5358i0 + 0.6149i0 + 2.0449i请输入由各节点的初始电压标幺值形成的列矩阵:V0=[1;1;1;1]请输入由短路号、短路点阻抗组成的矩阵:D=[4,0]请输入短路点的数目:NF=1短路点与其他各节点的互阻抗0 + 0.1732i0 + 0.5358i0 + 0.6149i0 + 2.0449i4点短路时电流的标幺值If=0 - 0.4890i各节点的电压标幺值U为(节点号从小到大排):0.91530.73800.6993各接地支路短路电流的标幺值I为(顺序同您输入B时一样0 + 4.5765i0 + 0.1845i各非接地支路短路电流的标幺值I为(顺序同您输入B时一样):0 - 0.4235i0 - 0.0655i0 - 0.4890i三、程序说明1、短路电流计算的计算机方法(MATLAB算法)短路电流计算实质上就是求解交流电路的稳态电流,其数学模型也就是网络的线性代数方程组。
电力系统短路故障分析的MATLAB辅助程序设计短路计算程序
![电力系统短路故障分析的MATLAB辅助程序设计短路计算程序](https://img.taocdn.com/s3/m/f35f47a484868762caaed556.png)
电力系统短路故障分析的MATLAB辅助程序设计电力系统短路故障可分为三相对称短路故障(three-phase balanced faults)和不对称短路故障(unbalanced faults )。
不对称短路故障又分为单相接地短路故障(single line-to-ground fault)、两相短路故障 (line-to-line fault)以及两相接地短路故障(double line-to-ground fault)。
根据故障分析结果可以对继电保护装置、自动装置进行整定计算,我们可以建立算法来形成节点阻抗矩阵,利用节点阻抗矩阵来计算短路故障情况下的节点电压和线路电流。
一、三相对称短路故障进行三相短路计算需要两个程序:zbuild /zbuildpi和symfault程序,zbuild、zbuildpi用来在MATLAB中形成节点阻抗矩阵,symfault用来计算三相对称故障。
Zbus=zbuild(zdata)这里的参数zdata是一个(e×4)阶矩阵,e是拓扑图的总支路数目。
第一列和第二列为元素两端的节点编号,第三列和第四列分别是线路的电阻、电抗的标幺值。
连接在0节点和发电机节点之间的发电机阻抗可能是次暂态电抗、暂态电抗或同步电抗,而且这个矩阵中还包含并联电抗器和负荷阻抗。
Zbus=zbuildpi(linedata,gendata,yload)这个函数与潮流计算程序是相容的,第一个参数linedat a与潮流计算程序中的文件是一致的。
第一列和第二列为节点编号;第三列到第五列分别是线路的电阻、电抗以及1/2线路电纳值,这三项都为在统一基准容量下的标幺值;最后一列是变压器分接头位置,对线路来说,必须输入1;线路无输入顺序。
发电机参数不包含在Linedata参数中,而是包含在第二个参数gendata中,gendata是一个g×4阶矩阵,g是发电机总数。
第一列和第二列为0节点、发电机节点编号,第三列和第四列为发电机的暂态电阻和暂态电抗。
电力系统三相短路的编程计算
![电力系统三相短路的编程计算](https://img.taocdn.com/s3/m/f2643a71ad51f01dc281f1ae.png)
电力系统分析课程专题报告学生姓名:班级:学号:指导教师:所在单位:提交日期:评分3电力系统三相短路的编程计算某某某(学院, )摘要:在电力系统中,三相短路故障造成的危害是最大的,发生的几率也最高,故短路计算对电力系统的稳定运行具有十分重要的意义。
作为电力系统三大计算之一,计算三相短路故障发生时的短路电流、各节点电压、各支路电流是短路计算的基本内容。
在电力系统短路电流的工程计算中,由于快速继电保护的应用,最重要的是计算短路电流基频交流分量的初始值,即次暂态电流I ''。
在给定电源电势时,实际上就是稳态交流电路的求解。
本文基于教材例3-2应用MATLAB 编程计算三相短路故障的电流电压情况,并于例3-3进行进一步的验证和完善。
关键词:三相短路;MATLAB中图分类号:TM 713 文献标识码:A0 引言在电力系统的四种短路类型中,三相短路是其中最严重的,其短路电流可达数万安以至十几万安,随之产生的热效应和电动力效应将使电气设备遭受严重破坏。
因此,计算短路电流的主要应用目的是电力系统设计中的电气设备选择,短路计算已成为电力系统运行分析、设计计算的重要环节。
实际电力系统短路电流交流分量初始值的计算,小型系统可以手算,而对于结构复杂的大型系统,短路电流计算量较大,用计算机进行辅助计算成为大势所趋。
1 解析法求三相短路电流1.1 参数说明(1) 为了元件参数标幺值计算方便,取基准容量B S 为A MV ⋅60,可设任意值,但必须唯一值参与计算。
(2) 取基准电压B U 为平均额定电压AV U ,基于例3-2中系统的额定电压等级有kv 10、kv 110,平均额定电压分别为kv 115、kv 5.10,平均额定电压与线路额定电压相差5%,为简化计算,故取平均额定电压。
(3) I ''为次暂态短路电流有效值,短路电流周期分量的初值等于时间0=t 时的有效值。
满足产生最大短路电流的三个条件下的最大次暂态短路电流作为计算依据。
发电机出口三相短路计算matlab
![发电机出口三相短路计算matlab](https://img.taocdn.com/s3/m/0465550fff4733687e21af45b307e87101f6f8ea.png)
发电机出口三相短路计算matlab发电机出口三相短路是指在电机出口端发生发电机绕组之间的相间短路故障。
这种故障可能会导致电机发生过电流和电磁力矩的异常增加,严重情况甚至会导致设备损坏。
因此,对发电机出口三相短路进行计算和分析非常重要。
在进行发电机出口三相短路计算时,可以借助Matlab工具进行相应的数值计算和仿真分析。
Matlab工具提供了一系列的电力系统仿真和计算函数,可以帮助我们进行发电机出口三相短路计算。
首先,我们可以使用Matlab中的复数运算功能来计算短路电流。
假设发电机出口处的电动势为E,短路电阻为Z,短路电流为I。
根据欧姆定律,我们可以得到I = E / Z。
在Matlab中,可以使用复数数学函数来计算电压和阻抗的复数表示,然后进行除法计算得到电流。
其次,我们需要考虑发电机出口三相短路的影响范围和电流传播路径。
可以使用Matlab中的电网等值方法来建立发电机的等值电路模型。
例如,可以使用Matlab中的电网等值工具箱来建立配电网的等值模型,并将发电机的等值电路与之相连。
然后,可以使用示波器模块来显示电流的传播路径和短路电流的分布情况。
此外,我们还可以借助Matlab的图形绘制功能来进行电流和电压波形的展示和分析。
通过绘制发电机出口端的电流和电压波形,可以清楚地观察到短路故障对电机运行的影响。
可以使用Matlab中的信号处理工具箱来进行频谱分析和谐波分析,以计算短路电流的谐波含量和幅值。
此外,我们还可以使用Matlab中的仿真工具来模拟发电机出口三相短路故障,并对故障后的电流和电压进行仿真分析。
例如,可以使用Matlab中的Simulink工具箱来建立电力系统的仿真模型,将发电机等值电路和负载模型相连,并对短路故障进行仿真。
可以通过仿真结果来判断发电机短路故障对系统稳定性和设备保护的影响。
综上所述,使用Matlab进行发电机出口三相短路计算非常方便。
通过Matlab提供的复数运算功能、电网等值方法、图形绘制功能和仿真工具,可以进行短路电流的计算、电流传播路径的分析、波形展示和仿真分析。
Matlab-Three-PhaseSource(三相电源)
![Matlab-Three-PhaseSource(三相电源)](https://img.taocdn.com/s3/m/481273c1250c844769eae009581b6bd97f19bcdb.png)
Matlab-Three-PhaseSource(三相电源)内部R-L 感性阻抗的三相电源的实现库电源模块库电源描述三相电源模块执⾏带有内部R L 阻抗的平衡的三相电压来源。
三电压来源在中与中性点连接成Y 或者其他⽅式。
你可以通过规定电源短路容量或者是直接规定电源内部的电阻和电感R 和L 或者是间接地设置X 与R ⽐率。
对话框和参数参数表线电压有效值内部线电压单位:VA 相初相⾓指定内部电压的相⾓为A 相相⾓,⽤度数表⽰。
三相电压的相序为正序。
这样,B 相和C 相的电压分别落后A 相120°和240°。
频率三相电源的频率单位:Hz内部连接⽅式三相电源的内部连接。
改变电源连接,电源模块图标就跟着改变。
选择其中⼀个连接⽅式:通过短路电流容量指定阻抗选择它表明你想通过短路电流容量和X/R 的⽐值来指定阻抗值。
基准电压下三相短路电流容量基准电压下三相电感短路功率单位:VA ,⽤来计算内部电感值。
只有选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
在给定基准相电压V base (单位:V )和电感三相短路功率Psc (单位:VA )时,内部电感L (单位:H )可由下式得到:212base sc V L P fπ=? 基准电压相基准电压,单位:V ,⽤来指定三相短路电流容量。
基准电压值通常就是名义电压源的⼤⼩。
只有选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
X/R ⽐值在标称电源频率或者X/R ⽐值内部电源阻抗品质因数下的X/R ⽐值。
只有选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
在给定内部电阻R (单位:Ω)和电感三相短路功率X (单位:Ω)时,X/R ⽐值可由下式得到:2(/)/X fL R X R X Rπ==只有不选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
电源电阻单位:Ω电感只有不选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
电感单位:H潮流计算表表中的这些参数是 Powergui模块潮流计算的⼯具,这些潮流计算参数只能给模型进⾏初始化,对块模型和仿真参数没有影响。
(完整版)短路计算matlab程序
![(完整版)短路计算matlab程序](https://img.taocdn.com/s3/m/2cd5d9e4f46527d3240ce0ce.png)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
(完整版)短路计算matlab程序
![(完整版)短路计算matlab程序](https://img.taocdn.com/s3/m/2cd5d9e4f46527d3240ce0ce.png)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
电力系统短路电流计算matlab程序
![电力系统短路电流计算matlab程序](https://img.taocdn.com/s3/m/81749eb1ed3a87c24028915f804d2b160b4e86ff.png)
电力系统短路电流计算matlab程序%电力系统极坐标下的牛顿-拉夫逊法潮流计算disp('电力系统极坐标下的牛顿-拉夫逊法潮流计算:');clearn=input('请输入结点数:n=');n1=input('请输入PV结点数:n1=');n2=input('请输入PQ结点数:n2=');isb=input('请输入平衡结点:isb=');pr=input('请输入精确度:pr=');K=input('请输入变比矩阵看:K=');C=input('请输入支路阻抗矩阵:C=');y=input('请输入支路导纳矩阵:y=');U=input('请输入结点电压矩阵:U=');S=input('请输入各结点的功率:S=');Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros( 1,n1+n2); for m=1:nfor R=1:nC(m,m)=C(m,m)+y(m,R);if K(m,R)~=0C(m,m)=C(m,m)+1/((K(m,R)*C(m,R))/(K(m,R)-1));C(R,R)=C(R,R)+1/((K(m,R)^2*C(m,R))/(1-K(m,R)));C(m,R)=C(m,R)*K(m,R);C(R,m)=C(m,R);endendendfor m=1:nfor R=1:nif m~=RZ(m)=Z(m)+1/C(m,R); endendendfor m=1:nfor R=1:nif m==RY(m,m)=C(m,m)+Z(m); elseY(m,R)=-1/C(m,R);endendenddisp('结点导纳矩阵:');disp(Y);disp('迭代中的雅克比矩阵:'); G=real(Y);B=imag(Y);O=angle(U);U1=abs(U);k=0;PR=1;P=real(S);Q=imag(S);while PR>prfor m=1:n2UD(m)=U1(m);endfor m=1:n1+n2for R=1:nPT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endPT1(m)=sum(PT);PP(m)=P(m)-PT1(m);PP1(k+1,m)=PP(m);endfor m=1:n2for R=1:nQT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R)));endQT1(m)=sum(QT);QQ(m)=Q(m)-QT1(m);QQ1(k+1,m)=QQ(m);endPR1=max(abs(PP));PR2=max(abs(QQ));PR=max(PR1,PR2);for m=1:n1+n2for R=1:n1+n2if m==RH(m,m)=U1(m)^2*B(m,m)+QT1(m);elseH(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R)));endendendfor m=1:n1+n2for R=1:n2if m==RN(m,m)=-U1(m)^2*G(m,m)-PT1(m);elseN(m,R)=-U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endendendfor m=1:n2for R=1:n1+n2if m==RJ(m,m)=U1(m)^2*G(m,m)-PT1(m);elseJ(m,R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endendendfor m=1:n2for R=1:n2if m==RL(m,m)=U1(m)^2*B(m,m)-QT1(m);elseL(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R)));endendendJJ=[H N;J L];PQ=[PP';QQ'];DA=-inv(JJ)*PQ;DA1=DA';for m=1:n1+n2OO(m)=DA1(m);endfor m=n:n1+n2+n2UU1(m-n1-n2)=DA1(m); endUD2=diag(UD);UU=UU1*UD2;for m=1:n1+n2O(m)=O(m)+OO(m); endfor m=1:n2U1(m)=U1(m)+UU(m); endfor m=1:n1+n2o(k+1,m)=180/pi*O(m); endfor m=1:n2u(k+1,m)=U1(m);endk=k+1;endfor m=1:nb(m)=U1(m)*cos(O(m)); c(m)=U1(m)*sin(O(m)); endfor R=1:nPH1(R)=U(isb)*conj(Y(isb,R))*conj(U(R));endPH=sum(PH1);for m=1:nfor R=1:nif m~=RC1(m,R)=1/C(m,R);elseC1(m,m)=C(m,m);endendendfor m=1:nfor R=1:nif (C(m,R)~=inf)&(m~=R)SS(m,R)=U1(m)^2*conj(C1(m,m))+U(m)*(conj(U(m))-conj(U(R)))*conj(C1(m,R));endendenddisp('迭代中的△P:');disp(PP1);disp('迭代中的△Q:');disp(QQ1);disp('迭代中相角:');disp(o);disp('迭代中电压的模:');disp(u);disp('平衡结点的功率:');disp(PH);disp('全部线路功率分布:');disp(SS);。
matlab三相断路器模块测量参数
![matlab三相断路器模块测量参数](https://img.taocdn.com/s3/m/996d03b1bdeb19e8b8f67c1cfad6195f312be8a3.png)
matlab三相断路器模块测量参数一、引言三相断路器是电力系统中常用的电气设备,用于保护电力系统免受过载、短路等故障的影响。
为了确保三相断路器的可靠性和安全性,需要对其进行各种参数的测量。
本文将介绍如何使用MATLAB三相断路器模块进行参数测量。
二、MATLAB三相断路器模块介绍MATLAB是一款常用的数学软件,具有强大的计算能力和良好的数据可视化功能。
在电力系统中,MATLAB可以通过三相断路器模块实现对三相断路器各种参数的测量。
三、使用MATLAB进行参数测量1. 安装MATLAB三相断路器模块首先需要安装MATLAB三相断路器模块,该模块可以从官方网站下载并安装。
2. 连接测试设备连接测试设备,包括电源、信号发生器和示波器等。
根据测试要求设置信号发生器输出信号频率和幅值,并将示波器连接到测试点上。
3. 打开MATLAB软件打开MATLAB软件,并加载已安装的三相断路器模块。
4. 设置测试参数根据测试要求设置测试参数,包括输入电压、负载电流等。
5. 进行参数测量在MATLAB中运行三相断路器模块,开始进行参数测量。
根据测试结果分析三相断路器的各项参数,包括开断时间、闭合时间、过流保护动作时间等。
四、MATLAB三相断路器模块应用案例以下是一个使用MATLAB三相断路器模块进行参数测量的应用案例。
1. 测试目标:测量三相断路器的开关时间和电流保护动作时间。
2. 测试设备:信号发生器、示波器、电源等。
3. 测试步骤:(1) 设置信号发生器输出频率和幅值。
(2) 将示波器连接到测试点上,并设置示波器触发方式。
(3) 在MATLAB中加载三相断路器模块,并设置测试参数,包括输入电压和负载电流等。
(4) 运行MATLAB程序,开始进行参数测量。
根据测试结果分析三相断路器的开关时间和电流保护动作时间是否符合要求。
4. 测试结果:经过多次测试,得出以下结论:(1) 三相断路器的开关时间稳定,平均开关时间为10ms左右。
电力系统短路故障的Matlab算法
![电力系统短路故障的Matlab算法](https://img.taocdn.com/s3/m/3e846eff0242a8956bece4b5.png)
2011-2012学年度下学期电力系统分析课程设计电力系统短路故障的计算机算法程序设计专业:电气工程及其自动化班级:学号:学生姓名:指导教师:2012年 2 月 28 日信息工程学院课程设计任务书目录课程设计说明 (4)1 任务提出与方案论证 (5)1.1 任务提出 (5)1.2 方案论证 (6)2 设计思路 (7)2.1 短路电流的计算方法 (7)2.2 利用节点阻抗矩阵计算短路电流 (8)3 详细设计 (10)3.1分析与计算 (10)3.2程序主框图及主要数据变量说明 (13)4 总结 (18)参考文献 (19)课程设计说明本文根据电力系统三相对称短路的特点,建立了合理的三相短路的数学模型,在此基础上,形成电力系统短路电流实用汁算方法即节点阻抗矩阵的支路追加法。
编制了对任意一个电力系统在任意点发生短路故障时三相短路电流及其分布的通用计算机程序,该办法适用予各种复杂结构的电力系统,从一个侧面展示了计算机应用于电力系统的广阔前景。
电力系统的短路故障是严重的,而又是发生几率最多的故障,一般说来,最严重的短路是三相短路。
当发生短路时,其短路电流可达数万安以至十几万安,它们所产生的热效应和电动力效应将使电气设备遭受严重破环。
为此,当发生短路时,继电保护装置必须迅速切除故障线路,以避免故障部分继续遭受危害,并使非故障部分从不正常运行情况下解脱出来,这要求电气设备必须有足够的机械强度和热稳定度,开关电气设备必须具备足够的开断能力,即必须经得起可能最大短路的侵扰而不致损坏。
因此,电力系统短路电流计算是电力系统运行分析,设计计算的重要环节,许多电业设计单位和个人倾注极大精力从事这一工作的研究。
在电力设计中选择电气设备必须计算短路电流, 短路电流的计算是电气专业设计不可缺少的环节,是电力设计中最重要的计算之一。
传统的设计方法中, 短路电流的计算是以手工计算形式进行的 ,先通过手工方法化简电力网络, 求出各电源点对短路点的转移阻抗, 从而求出计算用电抗XJS , 再查找运算曲线, 以求得短路电流的周期分量。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
电力系统三相短路计算a main.mclear tim%打开文件[dfile,pathname]=uigetfile('*.m' , 'Select Data File');if pathname == 0error( ' you must select a valid data file')elselfile =length(dfile);eval(dfile(1:lfile-2));end%定义输出文件output_file=fopen('output.dat', 'w' );%开始计时tic;%求解节点导纳矩阵,其中Ymatrix1 是考虑了变比,且支路未近似的导纳矩阵;Ymatrix2 是近似变比为 1,但是支路未近似计算的节点导纳矩阵;Ymatrix3 是近似变比为 1,采取近似支路参数 1的导纳矩阵; Ymatrix4 是近似变比为 1,采取近似支路参数2的导纳矩阵。
Y = Ymatrix2(bus,line);%对故障点进行导纳修正fixY = FixY(Y,bus,fault);%求注入电流Iinj = Inode(bus,calcSettings);U = fixY\Iinj;%得到故障支路与其他支路电流Bcurrent = Ibranch( line,U,fault,Y );%如果发生支路三相短路,那么对应该支路的电流修正为-999999-j999999Ib = ReviseBcurrent( fault,Bcurrent );%结束计时tim=toc;fprintf( '\n 程序运行结果 ' );fprintf( '\n 计算完成,共用时 %4.4fs, 相关结果已保存在 output.dat\n',tim);%输出结果fprintf_result(output_file, Ib);fprintf_result1(Ib);b FixY.mfunction fixY = FixY( Y,bus,fault )%对形成的导纳矩阵进行故障点的修正[nb,mb]=size(bus);[nf,mf]= size(fault);fixY = Y;%对发电机节点导纳修正for k=1:nbbusType=bus(k,7);if (busType==1)fixY(bus(k,1),bus(k,1)) = fixY(bus(k,1),bus(k,1)) + 1/1i/bus(k,8);endend%对节点短路和支路短路的导纳矩阵进行修正for k=1:nfnodeI=fault(k,1);nodeJ=fault(k,2);dis=fault(k,3);if (nodeI==0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue ;endif(nodeJ==0)fixY(nodeI,nodeI) = 999999+1i*999999;continue ;if(dis==0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 999999+1i*999999;continue ;endif(dis==1)&&(nodeI*nodeJ~=0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue ;endif(dis~=1)&&(dis~=0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = fixY(nodeI,nodeI) - fixY(nodeI,nodeJ)/dis;fixY(nodeJ,nodeJ) = fixY(nodeJ,nodeJ) - fixY(nodeI,nodeJ)/(1-dis); fixY(nodeI,nodeJ)=0;fixY(nodeJ,nodeI)=0;endendc fprintf_result.mfunction [ output_args ] = fprintf_result( output_file, Ib )%将得到的短路电流输入到输出文件中[n,m]=size(Ib);fprintf(output_file,' No. No. vector of I value of I\n' );for k=1:nI=Ib(k,1);J=Ib(k,2);I01=real(Ib(k,3));I02=imag(Ib(k,3));I1=Ib(k,4);if (I02>=0)fprintf( output_file,'%3d %3d %10.6f+j%10.6f %10.6f',I,J,I01,I02,I1);endif (I02<0)I02=abs(I02);fprintf( output_file,'%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf( output_file,'\n' );endendd fprintf_result1.mfunction [ output_args ] = fprintf_result1( Ib )%UNTITLED? ? ? ú′? ′|ê? è? oˉêy ? ? òa[n,m]=size(Ib);fprintf( ' No. No. vector of I value of I\n' );for k=1:nI=Ib(k,1);J=Ib(k,2);I01=real(Ib(k,3));I02=imag(Ib(k,3));I1=Ib(k,4);if (I02>=0)fprintf( '%3d %3d %10.6f+j%10.6f %10.6f' ,I,J,I01,I02,I1);endif (I02<0)I02=abs(I02);fprintf( '%3d %3d %10.6f-j%10.6f %10.6f' ,I,J,I01,I02,I1);endfprintf('\n' );endende Ibranch.mfunction Bcurrent = Ibranch( line,U,fault,Y )%计算短路电流%记录短路故障参数,如短路节点,如为支路短路,记录距离节点的距离%此段计算采用的支路参数未近似,如果计算近似的时候需要修改[nl,ml]=size(line);Bcurrent=zeros(nl+1,4);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);faultNode = 0;if (faultI==0)faultNode = faultJ;endif (faultJ==0)faultNode = faultI;endif (dis==1)&&(faultI*faultJ~=0) faultNode = faultJ;endif (dis==0)&&(faultI*faultJ~=0) faultNode = faultI;endif (faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Iij = 0;Iij1=0;end%计算非故障支路的短路电流for k=1:nli=line(k,1);j=line(k,2);Ui=U(i);if j~=0Uj=U(j);elseUj=0;endif line(k,2)==0Ym=line(k,5)+1i*line(k,6);Iij=Ui*Ym;Iij1=abs(Iij);endif line(k,2)~=0Zt=line(k,3)+1i*line(k,4);Yt=1/Zt;Ym=line(k,5)+1i*line(k,6);Iij=(Ui-Uj)*Yt+Ui*Ym;Iij1=abs(Iij);endBcurrent(k,1)=i;Bcurrent(k,2)=j;Bcurrent(k,3)=Iij;Bcurrent(k,4)=Iij1;end%如果为节点短路,修正短路点的电流大小if (faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Ifault = 0;branchCurrent=0;for k=1:nlI=line(k,1);J=line(k,2);if (I*J==0)continue ;endbranchCurrent = (U(I)-U(J))/(line(k,3)+1i*line(k,4));if(I==faultNode)Ifault = Ifault - branchCurrent ;else if (J==faultNode)Ifault = Ifault + branchCurrent ;endendendBcurrent(nl+1,3) = Ifault;Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3)); end%如果为支路短路,修正短路支路的短路电流大小if (dis~=0)&&(dis~=1)&&(faultI*faultJ~=0)Bcurrent(nl+1,1) = faultI;Bcurrent(nl+1,2) = faultJ;Bcurrent(nl+1,3) = U(faultI)*Y(faultI,faultJ)/dis + U(faultJ)*Y(faultI,faultJ)/(1-dis);Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));endendf Inode.mfunction Iinj = Inode( bus,calcSettings )%计算节点注入电流[nb,mb]=size(bus);Iinj = zeros(nb,1);for k=1:nbbusType=bus(k,7);if (calcSettings(1)==1)v = 1;elsev = bus(k,2);end%对发电机节点电流进行修正if (busType==1)Iinj(bus(k,1),1) = Iinj(bus(k,1),1) + v/1i/bus(k,8);endendendg ReviseBcurrent.mfunction Ib = ReviseBcurrent( fault,Bcurrent )%如果发生支路短路,对原来的计算电流进行修正,使该支路短路电流输出为 -999999-j999999clear faultI faultJ dis[nt,mt]=size(Bcurrent);Ib=zeros(nt,mt);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);for k=1:nt-1i=Bcurrent(k,1);j=Bcurrent(k,2);Ib(k,:)=Bcurrent(k,:);if(faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultI)&&(j==faultJ) Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endif(faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultJ)&&(j==faultI) Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endIb(nt,:)=Bcurrent(nt,:);endh Ymatrix1.mfunction Y = Ymatrix1( bus,line )%考虑变压器,并且支路参数不近似的节点导纳矩阵[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=1/Zt;Ym=line(k,5)+1i*line(k,6);K=line(k,7);if (K==0)&&(J~=0) Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif (K==0)&&(J==0) Y(I,I)=Y(I,I)+Ym;endif K>0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/(K*K); Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endif K<0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)+K*Yt;Y(J,I)=Y(I,J);endendendi Ymatrix2.mfunction Y = Ymatrix2( bus,line )%考虑变压器变比近似为1,支路参数不等效[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4); Yt=1/Zt;Ym=line(k,5)+1i*line(k,6);if J~=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif J==0Y(I,I)=Y(I,I)+Ym;endendendj Ymatrix3.mfunction Y = Ymatrix3( bus,line )%考虑变压器变比为 1,采用支路参数近似 1 [nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=imag(1/Zt);Ym=imag(line(k,5)+1i*line(k,6));if J~=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif J==0Y(I,I)=Y(I,I)+Ym;endendendk Ymatrix4.mfunction Y = Ymatrix4( bus,line )%变压器变比近似为 1,采用支路等效参数 2 [nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=1i*line(k,4);Yt=1/Zt;Ym=1i*line(k,6);if J~=0Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif J==0Y(I,I)=Y(I,I)+Ym;endend end。