电力系统潮流计算及Matlab编程实现

大家好，欢迎来到IT知识分享网。

2.1.3 节点电压用极坐标表示时的牛顿-拉夫逊法潮流计算

2.1.4 潮流计算程序框图

2.2 PQ分解法

3. MATLAB实例计算

1. 潮流计算：

复杂电力系统分析计算的一般方法是对整个电力系统建立数学模型，并通过计算机编程求出个节点的电压及电力系统中的功率分布。

2. 潮流计算常用算法：

2.1 牛顿-拉夫逊算法

2.1.1 牛顿-拉夫逊法的基本原理

牛顿–拉夫逊法(简称牛顿法)在数学上是求解非线性代数方程式的有效方
法。其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求
解的过程。即通常所称的逐次线性化过程。

对于非线性代数方程组：

$f(x^{(0)})+\dot{f}(x^{(0)})\Delta x^{(0)}=0$

上式称为牛顿法的修正方程式。由此可以求得第一次迭代的修正量：

$\Delta x^{(0)}=-[\dot{f}(x^{(0)})]^{-1}f(x^{(0)})$

$\dot{f}(x^{(k)})\Delta x^{(k)}=-f(x^{(k)})$

$x^{(k+1)=x^{(k)}}+\Delta x^{(k)}$

2.1.2 潮流计算的修正方程

n个独立节点的电力系统，其功率方程为：

$P_{i}-jQ_{i}=U_{i}\sum_{j=1}^{n}Y_{ij}U_{j}(i=1,2,3...n)$

将节点电压用极坐标表示，导纳矩阵用复数的形式表示并展开时，可以得到：

令

近似得到n+m-1个修正方程式：

2.1.3 节点电压用极坐标表示时的牛顿-拉夫逊法潮流计算

修正方程如下：

2.1.4 潮流计算程序框图

2.2 PQ分解法

略详情可见《电力系统分析》课本中。

3. MATLAB实例计算

系统框图：

原始数据：

（1）电力线路数据

I侧母线	J侧母线	电阻	电抗	电纳的1/2
1	2	0.0192	0.0575	0.0264
1	3	0.0452	0.1852	0.0204
2	4	0.057	0.1737	0.0184
3	4	0.0132	0.0379	0.0042
2	5	0.0472	0.1983	0.0209
2	6	0.0581	0.1763	0.0187
4	6	0.0119	0.0414	0.0045
5	7	0.046	0.116	0.0102
6	7	0.0267	0.082	0.0085
6	8	0.012	0.042	0.0045
9	11	0	0.208	0
9	10	0	0.11	0
12	13	0	0.14	0
12	14	0.1231	0.2559	0
12	15	0.0662	0.1304	0
12	16	0.0945	0.1987	0
14	15	0.221	0.1997	0
16	17	0.0824	0.1932	0
15	18	0.107	0.2185	0
18	19	0.0639	0.1292	0
19	20	0.034	0.068	0
10	20	0.0936	0.209	0
10	17	0.0324	0.0845	0
10	21	0.0348	0.0749	0
10	22	0.0727	0.1499	0
21	22	0.0116	0.0236	0
15	23	0.1	0.202	0
22	24	0.115	0.179	0
23	24	0.132	0.27	0
24	25	0.1885	0.3292	0
25	26	0.2554	0.38	0
25	27	0.1093	0.2087	0
27	29	0.2198	0.4153	0
27	30	0.3202	0.6027	0
29	30	0.2399	0.4533	0
8	28	0.0636	0.2	0.0214
6	28	0.0169	0.0599	0.0065

表2-1 电力线路数据

（2）变压器数据

I侧母线	J侧母线	电抗	变比
6	9	0.208	1.0155
6	10	0.556	0.9629
4	12	0.256	1.0129

表2-2 变压器数据

（3）发电机数据

母线名	母线类型	有功功率	无功功率	电压幅值	电压相角
2	PV	0.5756	–	1.0338	–
5	PV	0.2456	–	1.0058	–
8	PV	0.35	–	1.023	–
11	PV	0.1793	–	1.0913	–
13	PV	0.1691	–	1.0883	–
1	STACK	–	–	1.05	0

表2-3 发电机数据

（4）负荷数据

母线名	母线类型	有功功率	无功功率	电压幅值	电压相角
3	PQ	-0.024	-0.012	–	–
4	PQ	-0.076	-0.016	–	–
6	PQ	0	0	–	–
7	PQ	-0.228	-0.109	–	–
9	PQ	0	0	–	–
10	PQ	-0.058	-0.02	–	–
12	PQ	-0.112	-0.075	–	–
14	PQ	-0.062	-0.016	–	–
15	PQ	-0.082	-0.025	–	–
16	PQ	-0.035	-0.018	–	–
17	PQ	-0.09	-0.058	–	–
18	PQ	-0.032	-0.009	–	–
19	PQ	-0.095	-0.034	–	–
20	PQ	-0.022	-0.007	–	–
21	PQ	-0.175	-0.112	–	–
22	PQ	0	0	–	–
23	PQ	-0.032	-0.016	–	–
24	PQ	-0.087	-0.067	–	–
25	PQ	0	0	–	–
26	PQ	-0.035	-0.023	–	–
27	PQ	0	0	–	–
28	PQ	0	0	–	–
29	PQ	-0.024	-0.009	–	–
30	PQ	-0.106	-0.019	–	–
2	PV	-0.217	–	1.0338	–
5	PV	-0.942	–	1.0058	–
8	PV	-0.3	–	1.023	–

表2-4 负荷数据

代码实现：

chaoliu_freedom.m %复杂潮流计算 %编辑时间：2022/11/17 11：30 %编辑人：葛仪菲 %说明： %节点类型 标号 %PQ节点 1 %PV节点 2 %slack节点 3 %编号任意 不强制节点编号顺序 clear; clc %总结点数 IEEE30 n=30; %计数位 用于节点重编 Gen_PV = 1; Load_PQ = 0; Load_PV = 1; Line_PQ = 0; Line_PV = 0; %Change-----记录前编号 new_Change-----记录新编号 （按照PQ节点---PV节点---平衡节点） Change = zeros(1,n); new_Change = zeros(1,n); %导入发电机 负荷 电力线路 变压器 补偿电容原始数据 （编号任意） load 'Gen_free.mat'; load 'Load_free.mat'; load 'Line_data_free.mat'; load 'Trans_data_free.mat'; load 'Cap_data_free.mat' %%导入发电机 负荷 电力线路 变压器 补偿电容原始数据 （编号按照PQ节点---PV节点---平衡节点顺序） % load 'Gen.mat'; % load 'Load.mat'; % load 'Line_data.mat'; % load 'Trans_data.mat'; % load 'Cap_data.mat' %重新编号 PQ节点 0--m PV节点 m+1----n-1 平衡节点 n for num=1:size(Gen,1) if(Gen(num,2)==2) Change(num) = Gen(num,1); Gen(num,1)=n-size(Gen,1)+Gen_PV; Gen_PV = Gen_PV+1; new_Change(num) = Gen(num,1); elseif(Gen(num,2)==3) Change(num) = Gen(num,1); Gen(num,1)=n; new_Change(num) = Gen(num,1); end end for num=1:size(Load,1) if(Load(num,2)==1) Change(num+size(Gen,1)) = Load(num,1); Load(num,1)=1+Load_PQ; Load_PQ=Load_PQ+1; new_Change(num+size(Gen,1)) = Load(num,1); end end %编号列表对比 yes = [Change;new_Change]; %负荷更新编号 for num = 1:size(Load,1) if(Load(num,2)==2) for j = 1:length(Change) if(Load(num,1)==Change(j)) Load(num,1) = new_Change(j); break; end end end end %电力线路更新编号 for num=1:size(Line_data,1) for j = 1:length(Change) if(Line_data(num,1)==Change(j)) Line_data(num,1) = new_Change(j); break; end end for k = 1:length(Change) if(Line_data(num,2)==Change(k)) Line_data(num,2) = new_Change(k); break; end end end %变压器更新编号 for num=1:size(Trans_data,1) for j = 1:length(Change) if(Trans_data(num,1)==Change(j)) Trans_data(num,1) = new_Change(j); end if(Trans_data(num,2)==Change(j)) Trans_data(num,2) = new_Change(j); end end end %补偿电容更新编号 for num=1:size(Cap_data,1) for j = 1:length(Change) if(Cap_data(num,1)==Change(j)) Cap_data(num,1) = new_Change(j); end if(Cap_data(num,2)==Change(j)) Cap_data(num,2) = new_Change(j); end end end %电力线路参数 : I侧母线 J侧母线 阻抗 1/2接地导纳 Line = [Line_data(:,1:2) Line_data(:,3)+i*Line_data(:,4) Line_data(:,5)*i]; % %变压器参数： I侧母线 J侧母线 电抗 变比 Trans = [Trans_data(:,1:2) Trans_data(:,3)*i Trans_data(:,4)]; % %补偿电容参数： 母线 电纳 Cap = [Cap_data(:,1),Cap_data(:,2)*i]; %母线计数 Bus=ones(1,n); %1-牛顿拉夫逊法， 2-PQ分解法 mode=2; Tmax=10; %最大迭代次数 limit=1.0e-12; %要求精度 %变压器π型等效阻抗参数 Zt=zeros(size(Trans,1),3); Zt(:,1)=Trans(:,3)./Trans(:,4); Zt(:,2)=Trans(:,3)./(1-Trans(:,4)); Zt(:,3)=Trans(:,3)./(Trans(:,4).^2-Trans(:,4)); Trans_pi=[Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3)]; %提取同时为发电机和负荷的编号 num2 = 1; simple = zeros(1,n); for num=1:size(Gen,1) for num1=1:size(Load,1) if(Gen(num,1)==Load(num1,1)) simple(num2) = Gen(num,1); num2 = num2+1; end end end %计算PQ,PV节点数 n=numel(Bus); %总节点数 m=n-1; %PQ节点数 for num=1:size(Gen,1)%数组行数 if Gen(num,2)==2 %除去PV节点就是PQ节点 m=m-1; end end flag = 0; for num=1:size(Load,1) if Load(num,2)==2 for num1=1:length(simple) if(Load(num,1)==simple(num1)) flag = 0; break; else flag = 1; end end if(flag==1) m=m-1; end end end %PQ节点包含浮游节点，其PQ=0 %提取P,Q,U向量 P,Q为原始数据，Pi,Qi为计算结果 P=zeros(1,n); Q=zeros(1,n); U=ones(1,n); cita=zeros(1,n); %该处为角度制 for num=1:size(Gen,1) if Gen(num,2)==1 %PQ节点 P(Gen(num,1))=Gen(num,3); Q(Gen(num,1))=Gen(num,4); end if Gen(num,2)==2 %PV节点 P(Gen(num,1))=Gen(num,3); U(Gen(num,1))=Gen(num,4); end if Gen(num,2)==3 %slack节点 U(Gen(num,1))=Gen(num,3); cita(Gen(num,1))=Gen(num,4); end end for num=1:size(Load,1) if Load(num,2)==1 %PQ节点 P(Load(num,1))=P(Load(num,1))+Load(num,3); %考虑到节点可能同时存在发电机和负荷 进行累加 Q(Load(num,1))=Q(Load(num,1))+Load(num,4); end if Load(num,2)==2 %PV节点 P(Load(num,1))=P(Load(num,1))+Load(num,3); U(Load(num,1))=Load(num,4); end if Load(num,2)==3 %slack节点 U(Load(num,1))=Load(num,3); cita(Load(num,1))=Load(num,4); end end disp('初始条件：') disp('各节点有功：') disp(P); disp('各节点无功：') disp(Q); disp('各节点电压幅值：') disp(U); cita=(deg2rad(cita)); %角度转换成弧度 disp('各节点电压相角(度)：') disp(rad2deg(cita)); %节点导纳矩阵的计算 Y=zeros(n); %新建节点导纳矩阵 y=zeros(n); %网络中的真实导纳 %计算y(i,j) for num=1:size(Line,1) ii=Line(num,1); jj=Line(num,2); y(ii,jj)=1/Line(num,3); y(jj,ii)=y(ii,jj); end for num=1:size(Trans_pi,1) ii=Trans_pi(num,1); jj=Trans_pi(num,2); y(ii,jj)=1/Trans_pi(num,3); y(jj,ii)=y(ii,jj); end %计算y(i,i) for num=1:size(Line,1) ii=Line(num,1); jj=Line(num,2); y(ii,ii)=y(ii,ii)+Line(num,4); y(jj,jj)=y(jj,jj)+Line(num,4); end for num=1:size(Trans_pi,1) ii=Trans_pi(num,1); jj=Trans_pi(num,2); y(ii,ii)=y(ii,ii)+Trans_pi(num,4); y(jj,jj)=y(jj,jj)+Trans_pi(num,5); end %算上补偿电容 for num=1:size(Cap,1) ii=Cap(num,1); y(ii,ii)=y(ii,ii)+Cap(num,2); end %由y计算Y ysum=sum(y,1); for num=1:n for j=1:n if num==j Y(num,j)=ysum(num); %自导纳 else Y(num,j)=-y(num,j); %互导纳 end end end disp('节点导纳矩阵：'); disp(Y); G=real(Y); %电导矩阵 B=imag(Y); %电纳矩阵 %计算▲P ▲Q Pi Qi [dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita ); disp('有功不平衡量：'); disp(dP); disp('无功不平衡量：'); disp(dQ); for num=1:Tmax fprintf('第%d次迭代\n',num); %雅可比矩阵的计算 if(mode==1) J=Jacobi( n,m,U,cita,B,G,Pi,Qi ); disp('雅可比矩阵'); disp(J); end %求解节点电压修正量 if(mode==1) [dU,dcita]=Correct( n,m,U,dP,dQ,J ); else [dU,dcita]=PQ_LJ( n,m,dP,dQ,U,B ); end disp('电压、相角修正量：'); disp(dU); disp(rad2deg(dcita)); %修正节点电压 U=U+dU; cita=cita+dcita; disp('节点电压幅值：'); disp(U); disp('节点电压相角：'); disp(rad2deg(cita)); %计算功率不平衡量 [dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita ); disp('有功不平衡量：'); disp(dP); disp('无功不平衡量：'); disp(dQ); %收敛判断 if (max(abs(dP))<limit && max(abs(dQ))<limit ) break; end%if end%for %迭代结束，判断收敛 if (max(abs(dP))<limit && max(abs(dQ))<limit ) disp('计算收敛'); else disp('计算不收敛或未达到要求精度'); end %打印功率 fprintf('迭代总次数：%d\n', num); disp('节点电压幅值：'); disp(U); %电压阈值判断 for num=1:length(U) if(U(num)>1.05) disp('超过了节点电压的阈值'); elseif((U(num)<0.95)) disp('低于了节点电压的阈值'); end end subplot(2,2,1); plot([1:length(U)],U,'.b-','LineWidth',0.8); xlabel('节点号') ylabel('节点电压'); hold on; disp('节点电压相角：'); disp(rad2deg(cita)); deg_cita = rad2deg(cita); subplot(2,2,2); plot([1:length(deg_cita)],deg_cita,'*m-','LineWidth',0.8); deg_cita_range = max(deg_cita) - min(deg_cita); disp('相角的极差为：'); disp(deg_cita_range); xlabel('节点号') ylabel('节点电压相角'); disp('有功计算结果：'); disp(Pi); Pi_max = max(Pi); Pi_min = min(Pi); disp('有功功率最大值：'); disp(Pi_max); disp('有功功率最小值：'); disp(Pi_min); subplot(2,2,3); plot([1:length(Pi)],Pi,'+k-','LineWidth',0.8); xlabel('节点号') ylabel('节点有功功率'); disp('无功计算结果：'); disp(Qi); Qi_max = max(Qi); Qi_min = min(Qi); disp('无功功率最大值：'); disp(Qi_max); disp('无功功率最小值：'); disp(Qi_min); subplot(2,2,4); plot([1:length(Qi)],Qi,'.g-','LineWidth',0.8); xlabel('节点号') ylabel('节点无功功率'); %结果汇总 all = [P' Q' U' deg_cita']; 【2】Unbalanced.m function [ dP,dQ,Pi,Qi ] = Unbalanced( n,m,P,Q,U,G,B,cita ) %计算ΔPi有功的不平衡量 for i=1:n for j=1:n Pn(j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j))); end Pi(i)=sum(Pn); end dP=P(1:n-1)-Pi(1:n-1); %dP有n-1个 %计算ΔQi无功的不平衡量 for i=1:n for j=1:n Qn(j)=U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j))); end Qi(i)=sum(Qn); end dQ=Q(1:m)-Qi(1:m); %dQ有m个 end%func 【3】Jacobi.m function [ J ] = Jacobi( n,m,U,cita,B,G,Pi,Qi ) %雅可比矩阵的计算 %分块 H N K L %i!=j时 for i=1:n-1 for j=1:n-1 H(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j))); end end for i=1:n-1 for j=1:m N(i,j)=-U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j))); end end for i=1:m for j=1:n-1 K(i,j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j))); end end for i=1:m for j=1:m L(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j))); end end %i==j时 for i=1:n-1 H(i,i)=U(i).^2*B(i,i)+Qi(i); end for i=1:m N(i,i)=-U(i).^2*G(i,i)-Pi(i); end for i=1:m K(i,i)=U(i).^2*G(i,i)-Pi(i); end for i=1:m L(i,i)=U(i).^2*B(i,i)-Qi(i); end %合成雅可比矩阵 J=[H N;K L]; end 【4】Correct.m function [ dU,dcita ] = Correct( n,m,U,dP,dQ,J ) %求解节点电压修正量 for i=1:m Ud2(i,i)=U(i); end dPQ=[dP dQ]'; dUcita=(-inv(J)*dPQ)'; dcita=dUcita(1:n-1); dcita=[dcita 0]; dU=(Ud2*dUcita(n:n+m-1)')'; dU=[dU zeros(1,n-m)]; end 【5】PQ.m function [ dU,dcita ] = PQ_LJ( n,m,dP,dQ,U,B ) dP_U=dP./U(1:n-1); dQ_U=dQ./U(1:m); dUdcita=(-inv(B(1:n-1,1:n-1))*dP_U')'; dcita=dUdcita./U(1:n-1); dU=(-inv(B(1:m,1:m))*dQ_U')'; dU=[dU zeros(1,n-m)]; dcita=[dcita 0];%补零 end

潮流计算结果：

*负荷数据*		牛拉法				PQ法
母线名	母线类型	有功功率	无功功率	电压幅值	电压相角	有功功率	无功功率	电压幅值	电压相角
1	1	-0.024	-0.012	1.033449	-4.	-0.024	-0.012	1.033449	-4.7117
2	1	-0.076	-0.016	1.028907	-5.	-0.076	-0.016	1.028908	-5.
3	1	0	0	1.023545	-6.	0	0	1.023543	-6.
4	1	-0.228	-0.109	1.008546	-8.047579	-0.228	-0.109	1.008545	-8.047497
5	1	0	0	1.045219	-8.17787	0	0	1.04521	-8.
6	1	-0.058	-0.02	1.04232	-10.09924	-0.058	-0.02	1.042302	-10.09909
7	1	-0.112	-0.075	1.047429	-9.	-0.112	-0.075	1.047437	-9.
8	1	-0.062	-0.016	1.033616	-10.16719	-0.062	-0.016	1.033986	-10.15907
9	1	-0.082	-0.025	1.030007	-10.2991	-0.082	-0.025	1.029825	-10.303
10	1	-0.035	-0.018	1.037992	-9.	-0.035	-0.018	1.037985	-9.
11	1	-0.09	-0.058	1.03566	-10.24732	-0.09	-0.058	1.035645	-10.24711
12	1	-0.032	-0.009	1.022134	-10.93164	-0.032	-0.009	1.022007	-10.93398
13	1	-0.095	-0.034	1.020624	-11.11331	-0.095	-0.034	1.020531	-11.1148
14	1	-0.022	-0.007	1.025258	-10.91823	-0.022	-0.007	1.025188	-10.91923
15	1	-0.175	-0.112	1.030191	-10.57232	-0.175	-0.112	1.030165	-10.57228
16	1	0	0	1.030818	-10.56717	0	0	1.030792	-10.56714
17	1	-0.032	-0.016	1.022148	-10.77762	-0.032	-0.016	1.021975	-10.78122
18	1	-0.087	-0.067	1.020145	-11.06602	-0.087	-0.067	1.019992	-11.0693
19	1	0	0	1.023359	-11.0592	0	0	1.023281	-11.06047
20	1	-0.035	-0.023	1.005755	-11.47258	-0.035	-0.023	1.005342	-11.47838
21	1	0	0	1.033973	-10.79035	0	0	1.033944	-10.79036
22	1	0	0	1.018877	-6.	0	0	1.018871	-6.
23	1	-0.024	-0.009	1.014362	-11.99441	-0.024	-0.009	1.014295	-11.99606
24	1	-0.106	-0.019	1.003018	-12.85807	-0.106	-0.019	1.002936	-12.85988
25	2	0.3586	0	1.0338	-2.	0.3586	0	1.0338	-2.73305
26	2	-0.6964	0	1.0058	-8.	-0.6964	0	1.0058	-8.
27	2	0.05	0	1.023	-6.	0.05	0	1.023	-6.47136

*发电机数据*		牛拉法				PQ法
母线名	母线类型	有功功率	无功功率	电压幅值	电压相角	有功功率	无功功率	电压幅值	电压相角
25	2	0.3586	0	1.0338	-2.	0.3586	0	1.0338	-2.73305
26	2	-0.6964	0	1.0058	-8.	-0.6964	0	1.0058	-8.
27	2	0.05	0	1.023	-6.	0.05	0	1.023	-6.47136
28	2	0.1793	0	1.0913	-6.	0.1793	0	1.0913	-6.
29	2	0.1691	0	1.0883	-8.056039	0.1691	0	1.0883	-8.055297
30	3	0	0	1.05	0	0	0	1.05	0

参考：

潮流计算的matlab程序实现方法_北窗北川的博客-CSDN博客_csdnmatlab潮流计算

免责声明：本站所有文章内容,图片，视频等均是来源于用户投稿和互联网及文摘转载整编而成，不代表本站观点，不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益，请在线联系站长,一经查实,本站将立刻删除。本文来自网络,若有侵权，请联系删除，如若转载，请注明出处：https://haidsoft.com/139420.html

电力系统潮流计算及Matlab编程实现

1. 潮流计算：

2. 潮流计算常用算法：

2.1 牛顿-拉夫逊算法

2.1.1 牛顿-拉夫逊法的基本原理

2.1.2 潮流计算的修正方程

2.1.3 节点电压用极坐标表示时的牛顿-拉夫逊法潮流计算

2.1.4 潮流计算程序框图

2.2 PQ分解法

3. MATLAB实例计算

相关推荐

发表回复