前言
闲着没事编了一个潮流计算,极坐标,牛-拉法,就当练手了,供大家参考~
1. 主函数
# coding=UTF-8
import NetData
from YData import GetY
from InitateU import PolarU
from NewtonRaphson import PolarNR
from OutTxt import SingleTxt,StringTxt
# 获取支路参数
RootPath = 'F:\\FunctionLibrary\\runpfM\\' # 读取文件,得到节点和支路数据
FilePath_Node = RootPath+'NodeData.txt'
FilePath_Line = RootPath+'LineData.txt'
Out_Path = 'Cal-Process.txt'
S0 = '-'
for i in range(130):
S0 = S0+'-'
StringTxt(path=Out_Path,string=S0+'\n初始化参数:\n',fmt='w')
NodeData = NetData.GetNodeData(FilePath_Node,show=1)
LineData = NetData.GetLineData(FilePath_Line,show=1)
# 节点导纳矩阵
Y = GetY(NodeData,LineData,path=Out_Path,width=6)
# 初始化
U,Angle = PolarU(NodeData[:,8],NodeData[:,9],path=Out_Path,width=9)
# 开始计算
PQNode,PVNode,SlackNode,P_Real,Q_Real = NetData.GetNetData(NodeData,LineData,path=Out_Path)
Iter = 0
MaxIter = 10
Tol = 1e-4
while True:
Iter = Iter+1
StringTxt(path=Out_Path,string='\n\n'+S0)
SingleTxt(Iter,path=Out_Path,string='迭代次数:')
U,Angle,MaxError = PolarNR(U,Angle,Y,PQNode,PVNode,SlackNode,P_Real,Q_Real,Tol,path=Out_Path,width=9)
if Iter>MaxIter or MaxError<Tol:
break
# 结束循环,输出结果
if MaxError<Tol:
SingleTxt(Iter-1,path=Out_Path,string='迭代完成,更新次数为:')
SingleTxt(MaxError,path=Out_Path,string='最大误差为:')
else:
SingleTxt(MaxError,path=Out_Path,string='结果不收敛!')
2.获取网络参数
# coding=UTF-8
import numpy as np
#----------------------------NodeData-------------------------#
# NodeData: bus_i type Pd Qd Pg Qg Gs Bs Vm Va
# bus_i为节点编号
# Pd和Qd为节点的负荷功率
# Pg和Qg为节点的发电机输出功率
# Gs和Bs为节点对地电导和电纳
# Vm为电压幅值,Va为电压相角
def GetNumNode(FilePath_Node):
with open(FilePath_Node,'r') as file:
NumNode = 0
for line in file:
NumNode = NumNode+1 # 得到节点数目
print('节点数目:',NumNode)
return(NumNode)
def GetNodeData(FilePath_Node,*NumNode,**PlotOption):
if not NumNode: # 空的
NumNode = GetNumNode(FilePath_Node)
NodeData = np.zeros([NumNode,10]) # 初始化
i = 0
with open(FilePath_Node,'r') as file:
for line in file:
NodeData[i,:] = np.array([float(i) for i in line.split()])
i = i+1
if 'show' in PlotOption and PlotOption['show'] == 1:
print('NodeData=',NodeData)
return(NodeData)
#----------------------------LineData------------------------#
# LineData:fbus tbus r x b ratio
# fbus为起始节点
# tbus为终止节点
# r为支路电阻
# x为支路电抗
# b为节点对地导纳
# ratio为变压器变比
def GetNumLine(FilePath_Line):
with open(FilePath_Line,'r') as file:
NumLine = 0
for line in file:
NumLine = NumLine+1 # 得到节点数目
print('支路数目:',NumLine)
return(NumLine)
def GetLineData(FilePath_Line,*NumLine,**PlotOption):
if not NumLine: # 空的
NumLine = GetNumLine(FilePath_Line)
LineData = np.zeros([NumLine,6]) # 初始化
i = 0
with open(FilePath_Line,'r') as file:
for line in file:
LineData[i,:] = np.array([float(i) for i in line.split()])
i = i+1
if 'show' in PlotOption and PlotOption['show'] == 1:
print('LineData = ',LineData)
return(LineData)
#-输出网络参数--------#
def GetNetData(NodeData,LineData,**Option):
PQNode = NodeData[np.where(NodeData[:,1]==1),0] # PQ节点
PVNode = NodeData[np.where(NodeData[:,1]==2),0] # PV节点
SlackNode = NodeData[np.where(NodeData[:,1]==3),0] # 平衡节点
# SlackNode = SlackNode[0]
# print(SlackNode)
P_Real = -NodeData[:,2]+NodeData[:,4] # 节点输入有功功率
Q_Real = -NodeData[:,3]+NodeData[:,5] # 节点输入无功功率
from OutTxt import Real
flag = 0
if 'string' not in Option:
flag = 1
Option['string'] = 'PQ节点:\n'
Real(PQNode,**Option)
if flag:
Option['string'] = 'PV节点:\n'
Real(PVNode,**Option)
if flag:
Option['string'] = '平衡节点:\n'
Real(SlackNode,**Option)
if flag:
Option['string'] = '注入节点有功功率:\n'
Real(P_Real,**Option)
if flag:
Option['string'] = '注入节点无功功率:\n'
Real(Q_Real,**Option)
return(PQNode,PVNode,SlackNode,P_Real,Q_Real)
3.节点导纳矩阵
# coding=UTF-8
# 形成节点导纳矩阵
# NodeData: bus_i type Pd Qd Pg Qg Gs Bs Vm Va
# LineData: fbus tbus r x b_2 ratio
import numpy as np
import OutTxt
def GetY(NodeData,LineData,**Option):
NumNode = NodeData.shape[0]
NumLine = LineData.shape[0]
Y = np.zeros([NumNode,NumNode])+np.zeros([NumNode,NumNode])*1j
G0 = NodeData[:,6] # 节点对地电导
B0 = NodeData[:,7] # 节点对地电纳
for i in range(NumLine):
Node1 = int(LineData[i,0]-1)
Node2 = int(LineData[i,1]-1)
# print(Node1,Node2)
R = LineData[i,2]
X = LineData[i,3]
if LineData[i,5]==0: # 普通线路,无变压器
B_2 = LineData[i,4]
Y[Node1,Node1] = Y[Node1,Node1]+B_2*1j+1/(R+1j*X)
Y[Node2,Node2] = Y[Node2,Node2]+B_2*1j+1/(R+1j*X)
Y[Node1,Node2] = Y[Node1,Node2]-1/(R+1j*X)
Y[Node2,Node1] = Y[Node2,Node1]-1/(R+1j*X)
else: # 有变压器支路
K = LineData[i,5]
YT = 1/(R+1j*X)
Y[Node1,Node1] = Y[Node1,Node1]+(K-1)/K*YT+YT/K
Y[Node2,Node2] = Y[Node2,Node2]+(1-K)/K**2*YT+YT/K
Y[Node1,Node2] = Y[Node1,Node2]-1/K*YT
Y[Node2,Node1] = Y[Node2,Node1]-1/K*YT
for i in range(NumLine):
Node = int(NodeData[i,0]-1) # 第一列为节点编号
Y[Node,Node] = Y[Node,Node]+G0[i]+1j*B0[i]
#----------输出到txt文件--------#
if 'string' not in Option:
Option['string'] = 'Y矩阵:\n'
OutTxt.Complex(Y,**Option) # 元组
return(Y)
4.电压初始化
# coding=UTF-8
#---------------极坐标形式----------------#
from numpy import cos,sin,angle
import OutTxt
def PolarU(U_Amplitude,U_Angle,**Option):
U_Complex = U_Amplitude*cos(U_Angle)+U_Amplitude*sin(U_Angle)*1j
U = abs(U_Complex)
Angle = angle(U_Complex)
if 'string' not in Option:
Option['string'] = '\n电压初始化为:\n'
OutTxt.Complex(U_Complex,**Option)
return(U,Angle)
#---------------直角坐标形式--------------#
def RectanU(U_Amplitude,U_Angle,**Option):
U_e = U_Amplitude*cos(U_Angle)
U_f = U_Amplitude*cos(U_Angle)
if 'string' not in Option:
Option['string'] = '\n 电压实部初始化为:\n'
OutTxt.Real(U_e,**Option)
Option['string'] = '\n 电压虚部初始化为:\n'
OutTxt.Real(U_f,**Option)
return(U_e,U_f)
4.牛-拉法求解
# coding=UTF-8
# 形成节点导纳矩阵
import numpy as np
from OutTxt import Real
def PolarNR(U,Angle,Y,PQNode,PVNode,SlackNode,P_Real,Q_Real,Tol,**Option):
P_iter = 0 # 为形成雅可比矩阵
Q_iter = 0 # 为形成雅可比矩阵
PQNode = PQNode-1
PVNode = PVNode-1
SlackNode = SlackNode-1
NumNode = Y.shape[0] # 节点数目
NumPQ = max(PQNode.shape) # PQ节点数目
G = Y.real
B = Y.imag
P = np.zeros([NumNode,1])
Q = np.zeros([NumNode,1])
DeltaP = np.zeros([NumNode-1,1])
DeltaQ = np.zeros([NumPQ,1])
for i in range(NumNode): # 求解功率不平衡量
P[i] = U[i]*np.sum(U*(G[i,:]*np.cos(Angle[i]-Angle) + B[i,:]*np.sin(Angle[i]-Angle)))
Q[i] = U[i]*np.sum(U*(G[i,:]* np.sin(Angle[i]-Angle) - B[i,:]*np.cos(Angle[i]-Angle)))
if i!=SlackNode: # 不是平衡节点
DeltaP[P_iter] = P_Real[i]-P[i] # NumPQ+NumPV
if i in PQNode: # PQ节点
DeltaQ[Q_iter] = Q_Real[i]-Q[i] # NumPQ
Q_iter = Q_iter+1
P_iter = P_iter+1
DeltaPQ = np.vstack([DeltaP,DeltaQ]) # 功率不平衡量
Option['string'] = '功率不平衡量为:\n'
Real(DeltaPQ,**Option)
MaxError = np.max(np.abs(DeltaPQ))
print(MaxError)
if MaxError<Tol:
return(U,Angle,MaxError)
## H and N
HN_iter = -1 # 初始化雅可比矩阵
H = np.zeros([NumNode-1,NumNode-1])
N = np.zeros([NumNode-1,NumPQ])
# H and N
for i in range(NumNode):
if i!=SlackNode: # PQ或PV节点
H_iter_y = -1
N_iter_y = -1
HN_iter = HN_iter+1 # 记录H和N的行数
for j in range(NumNode):
if j != SlackNode:
H_iter_y = H_iter_y+1 # 记录H列数
if i != j: # 非平衡节点计算H矩阵
Angleij = Angle[i]-Angle[j]
H[HN_iter,H_iter_y] = -U[i]*U[j]*(G[i,j]*np.sin(Angleij)-B[i,j]*np.cos(Angleij))
else:
H[HN_iter,H_iter_y] = Q[i]+U[i]**2*B[i,i]
if j in PQNode:
N_iter_y = N_iter_y+1 # 记录N的列数
if i != j:
Angleij = Angle[i]-Angle[j]
N[HN_iter,N_iter_y] = -U[i]*U[j]*(G[i,j]*np.cos(Angleij)+B[i,j]*np.sin(Angleij))
else:
N[HN_iter,N_iter_y] = -P[i]-G[i,i]*U[i]**2
## J and L
JL_iter = -1 # 初始化雅可比矩阵
J = np.zeros([NumPQ,NumNode-1])
L = np.zeros([NumPQ,NumPQ])
for i in range(NumNode):
if i in PQNode: # PQ节点
JL_iter = JL_iter+1 # J和L的行数
J_iter_y = -1
L_iter_y = -1
for j in range(NumNode):
if j!=SlackNode: # 非平衡节点
J_iter_y = J_iter_y+1
if i!=j:
Angleij = Angle[i]-Angle[j]
J[JL_iter,J_iter_y] = U[i]*U[j]*(G[i,j]*np.cos(Angleij)+B[i,j]*np.sin(Angleij))
else:
J[JL_iter,J_iter_y] = -P[i]+G[i,i]*U[i]**2
if j in PQNode: # PQ节点
L_iter_y = L_iter_y+1
if i!=j:
Angleij = Angle[i]-Angle[j]
L[JL_iter,L_iter_y] = -U[i]*U[j]*(G[i,j]*np.sin(Angleij)-B[i,j]*np.cos(Angleij))
else:
L[JL_iter,L_iter_y] = -Q[i]+B[i,i]*U[i]**2
# 修正
Jaccobi = np.vstack([np.hstack([H,N]),np.hstack([J,L])])
Option['string'] = 'jacobi矩阵为:\n'
Real(Jaccobi,**Option)
Delta = np.linalg.solve(Jaccobi,DeltaPQ)
Option['string'] = '方程组求解结果:\n'
Real(Delta,**Option)
DeltaAngle = Delta[0:NumNode-1] # 注意下标
DeltaU_U = Delta[NumNode-1:]
DA_iter = -1
U_U_iter = -1
for i in range(NumNode):
if i!=SlackNode:
DA_iter = DA_iter+1
Angle[i] = Angle[i]-DeltaAngle[DA_iter]
if i in PQNode:
U_U_iter = U_U_iter+1
U[i] = U[i]-U[i]*DeltaU_U[U_U_iter]
Option['string'] = '更新之后的电压幅值为:\n'
Real(U,**Option)
Option['string'] = '相角为:\n'
Real(Angle,**Option)
return(U,Angle,MaxError)
5.输出文件
# coding=UTF-8
def ArgOption(**Option): # 解析可选参数
if 'string' in Option:
String = Option['string']
else:
String = '输出的矩阵:\n' # 默认值
if 'path' in Option:
FilePath = Option['path']
else:
FilePath = 'untitled.txt' # 默认为未命名的
if 'width' in Option:
N = Option['width']
else:
N = 8 # 默认值
if 'fmt' in Option:
Fmt = Option['fmt']
else:
Fmt = 'a' # 默认值
return(String,FilePath,N,Fmt)
#---------------------实数矩阵-----------------------#
def Real(Matrix,**Option):
string,FilePath,N,Fmt = ArgOption(**Option)
with open(FilePath,Fmt) as file:
file.write(string)
Ndim = Matrix.ndim
if Ndim == 1: # (Num)
Num = Matrix.shape[0]
for i in range(Num):
file.write(str(Dot10(Matrix[i])).ljust(N+1)[0:N-1])
file.write(' ')
file.write('\n')
elif Ndim == 2: # ([NumR,NumC])
NumR,NumC = Matrix.shape
for i in range(NumR):
for j in range(NumC):
file.write(str(Dot10(Matrix[i,j])).ljust(N+1)[0:N-1])
file.write(' ')
file.write('\n')
#-------------------复数矩阵-----------------------#
def Complex(Matrix,**Option):
string,FilePath,N,Fmt = ArgOption(**Option)
RealM = Matrix.real
ImagM = Matrix.imag
# print(FilePath)
with open(FilePath,Fmt) as file:
file.write(string)
Ndim = Matrix.ndim
if Ndim == 1: # (Num)
Num = Matrix.shape[0]
for i in range(Num):
file.write(str(Dot10(RealM[i])).ljust(N+1)[0:N-1])
file.write('+')
file.write('j'+str(Dot10(ImagM[i])).ljust(N+1)[0:N-1])
file.write(' ')
file.write('\n')
elif Ndim == 2: # ([NumR,NumC])
NumR,NumC = Matrix.shape
for i in range(NumR):
for j in range(NumC):
file.write(str(Dot10(RealM[i,j])).ljust(N+1)[0:N-1])
file.write('+')
file.write('j'+str(Dot10(ImagM[i,j])).ljust(N+1)[0:N-1])
file.write(' ')
file.write('\n')
#--单变量-----------------------#
def SingleTxt(Number,**Option):
string,FilePath,N,Fmt = ArgOption(**Option)
with open(FilePath,Fmt) as file:
file.write(string)
file.write(str(Dot10(Number)).ljust(N+1)[0:N-1])
file.write('\n')
#--字符串-----------------------#
def StringTxt(**Option):
string,FilePath,N,Fmt = ArgOption(**Option)
with open(FilePath,Fmt) as file:
file.write(string)
file.write('\n')
#--小数点后保留位数---
def Dot10(x):
if isinstance(x,int): # 整数的话,就不用变为浮点数
return(x)
else:
return('{:.10f}'.format(x))
6.屏幕输出
节点数目: 9
NodeData= [[1. 1. 1.25 0.5 0. 0. 0. 0. 1. 0. ]
[2. 1. 0. 0. 0. 0. 0. 0. 1. 0. ]
[3. 1. 0.9 0.3 0. 0. 0. 0. 1. 0. ]
[4. 1. 0. 0. 0. 0. 0. 0. 1. 0. ]
[5. 1. 1. 0.35 0. 0. 0. 0. 1. 0. ]
[6. 1. 0. 0. 0. 0. 0. 0. 1. 0. ]
[7. 2. 0. 0. 1.63 0. 0. 0. 1.025 0. ]
[8. 2. 0. 0. 0.85 0. 0. 0. 1.025 0. ]
[9. 3. 0. 0. 0. 0. 0. 0. 1.04 0. ]]
支路数目: 9
LineData = [[1. 2. 0.01 0.085 0.088 0. ]
[1. 6. 0.032 0.161 0.153 0. ]
[2. 3. 0.017 0.092 0.079 0. ]
[2. 9. 0. 0.058 0. 1. ]
[3. 4. 0.039 0.17 0.179 0. ]
[4. 5. 0.012 0.101 0.105 0. ]
[4. 8. 0. 0.059 0. 1. ]
[5. 6. 0.019 0.072 0.075 0. ]
[6. 7. 0. 0.063 0. 1. ]]
1.63
0.18620555300855945
0.002147363829061549
3.6215137443425437e-07
7.txt文件输出
-----------------------------------------------------------------------------------------------------------------------------------
初始化参数:
Y矩阵:
2.552+j-17.3 -1.36+j11.60 0.000+j0.000 0.000+j0.000 0.000+j0.000 -1.18+j5.975 0.000+j0.000 0.000+j0.000 0.000+j0.000
-1.36+j11.60 3.307+j-39.1 -1.94+j10.51 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j17.24
0.000+j0.000 -1.94+j10.51 3.224+j-15.8 -1.28+j5.588 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000
0.000+j0.000 0.000+j0.000 -1.28+j5.588 2.441+j-32.0 -1.15+j9.763 0.000+j0.000 0.000+j0.000 0.000+j16.94 0.000+j0.000
0.000+j0.000 0.000+j0.000 0.000+j0.000 -1.15+j9.763 4.586+j-22.5 -3.42+j12.98 0.000+j0.000 0.000+j0.000 0.000+j0.000
-1.18+j5.975 0.000+j0.000 0.000+j0.000 0.000+j0.000 -3.42+j12.98 4.614+j-34.6 0.000+j15.87 0.000+j0.000 0.000+j0.000
0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j15.87 0.000+j-15.8 0.000+j0.000 0.000+j0.000
0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j16.94 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j-16.9 0.000+j0.000
0.000+j0.000 0.000+j17.24 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j0.000 0.000+j-17.2
电压初始化为:
1.000000+j0.000000 1.000000+j0.000000 1.000000+j0.000000 1.000000+j0.000000 1.000000+j0.000000 1.000000+j0.000000 1.025000+j0.000000 1.025000+j0.000000 1.040000+j0.000000
PQ节点:
1.00000 2.00000 3.00000 4.00000 5.00000 6.00000
PV节点:
7.00000 8.00000
平衡节点:
9.00000
注入节点有功功率:
-1.2500 0.00000 -0.9000 0.00000 -1.0000 0.00000 1.63000 0.85000 0.00000
注入节点无功功率:
-0.5000 0.00000 -0.3000 0.00000 -0.3500 0.00000 0.00000 0.00000 0.00000
-----------------------------------------------------------------------------------------------------------------------------------
迭代次数:1
功率不平衡量为:
-1.25000
0.000000
-0.90000
0.000000
-1.00000
0.000000
1.630000
0.850000
-0.25900
0.856655
-0.04200
0.707728
-0.17000
0.624825
jacobi矩阵为:
-17.5792 11.60409 -0.00000 -0.00000 -0.00000 5.975134 -0.00000 -0.00000 -2.55279 1.365187 -0.00000 -0.00000 -0.00000 1.187604
11.60409 -40.0458 10.51068 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 1.365187 -3.30737 1.942191 -0.00000 -0.00000 -0.00000
-0.00000 10.51068 -16.0989 5.588244 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 1.942191 -3.22420 1.282009 -0.00000 -0.00000
-0.00000 -0.00000 5.588244 -32.7242 9.763170 -0.00000 -0.00000 17.37288 -0.00000 -0.00000 1.282009 -2.44198 1.159980 -0.00000
-0.00000 -0.00000 -0.00000 9.763170 -22.7478 12.98467 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 1.159980 -4.58649 3.426510
5.975134 -0.00000 -0.00000 -0.00000 12.98467 -35.2296 16.26984 -0.00000 1.187604 -0.00000 -0.00000 -0.00000 3.426510 -4.61411
-0.00000 -0.00000 -0.00000 -0.00000 -0.00000 16.26984 -16.2698 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000
-0.00000 -0.00000 -0.00000 17.37288 -0.00000 -0.00000 -0.00000 -17.3728 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000
2.552792 -1.36518 0.000000 0.000000 0.000000 -1.18760 0.000000 0.000000 -17.0972 11.60409 -0.00000 -0.00000 -0.00000 5.975134
-1.36518 3.307378 -1.94219 0.000000 0.000000 0.000000 0.000000 0.000000 11.60409 -38.3325 10.51068 -0.00000 -0.00000 -0.00000
0.000000 -1.94219 3.224200 -1.28200 0.000000 0.000000 0.000000 0.000000 -0.00000 10.51068 -15.5829 5.588244 -0.00000 -0.00000
0.000000 0.000000 -1.28200 2.441989 -1.15998 0.000000 0.000000 0.000000 -0.00000 -0.00000 5.588244 -31.3088 9.763170 -0.00000
0.000000 0.000000 0.000000 -1.15998 4.586491 -3.42651 0.000000 0.000000 -0.00000 -0.00000 -0.00000 9.763170 -22.3878 12.98467
-1.18760 0.000000 0.000000 0.000000 -3.42651 4.614114 0.000000 0.000000 5.975134 -0.00000 -0.00000 -0.00000 12.98467 -33.9799
方程组求解结果:
0.067047
0.037365
0.063033
-0.04188
-0.01900
-0.07300
-0.17318
-0.09080
-0.00902
-0.03343
-0.02187
-0.03861
-0.02232
-0.03884
更新之后的电压幅值为:
1.009020 1.033434 1.021873 1.038614 1.022326 1.038845 1.025000 1.025000 1.040000
相角为:
-0.06704 -0.03736 -0.06303 0.041882 0.019005 0.073003 0.173189 0.090808 0.000000
-----------------------------------------------------------------------------------------------------------------------------------
迭代次数:2
功率不平衡量为:
0.039962
-0.01068
0.042708
-0.02549
0.052979
-0.04120
-0.06048
-0.03247
-0.07164
-0.05002
-0.04221
-0.08189
-0.02740
-0.18620
jacobi矩阵为:
-18.0808 12.05267 0.000000 0.000000 0.000000 6.028138 0.000000 0.000000 -1.30909 1.782033 -0.00000 -0.00000 -0.00000 2.106986
12.13716 -41.8034 11.14866 0.000000 0.000000 0.000000 0.000000 0.000000 1.063829 -3.54292 1.765476 -0.00000 -0.00000 -0.00000
-0.00000 11.04338 -16.7992 5.755881 0.000000 0.000000 0.000000 0.000000 -0.00000 2.335229 -2.42408 1.974267 -0.00000 -0.00000
-0.00000 -0.00000 6.040862 -34.4550 10.39202 0.000000 0.000000 18.02213 -0.00000 -0.00000 0.732043 -2.65971 0.994216 -0.00000
-0.00000 -0.00000 -0.00000 10.33567 -23.9093 13.57371 0.000000 0.000000 -0.00000 -0.00000 -0.00000 1.468480 -3.74059 4.378073
6.375689 -0.00000 -0.00000 -0.00000 13.96653 -37.1593 16.81710 0.000000 0.358370 -0.00000 -0.00000 -0.00000 2.889494 -5.02075
-0.00000 -0.00000 -0.00000 -0.00000 -0.00000 16.81710 -16.8171 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -1.69048
-0.00000 -0.00000 -0.00000 18.02213 -0.00000 -0.00000 0.000000 -18.0221 -0.00000 -0.00000 -0.00000 -0.88247 -0.00000 -0.00000
3.889020 -1.78203 0.000000 0.000000 0.000000 -2.10698 0.000000 0.000000 -17.2241 12.05267 0.000000 0.000000 0.000000 6.028138
-1.06382 3.521544 -1.76547 0.000000 0.000000 0.000000 0.000000 0.000000 12.13716 -41.9034 11.14866 0.000000 0.000000 0.000000
0.000000 -2.33522 4.309497 -1.97426 0.000000 0.000000 0.000000 0.000000 -0.00000 11.04338 -16.2837 5.755881 0.000000 0.000000
0.000000 0.000000 -0.73204 2.608730 -0.99421 0.000000 0.000000 -0.88247 -0.00000 -0.00000 6.040862 -34.6188 10.39202 0.000000
0.000000 0.000000 0.000000 -1.46848 5.846554 -4.37807 0.000000 0.000000 -0.00000 -0.00000 -0.00000 10.33567 -23.2642 13.57371
-0.35837 0.000000 0.000000 0.000000 -2.88949 4.938352 -1.69048 0.000000 6.375689 -0.00000 -0.00000 -0.00000 13.96653 -37.5317
方程组求解结果:
0.003254
0.001876
0.002173
0.008650
0.007655
0.009333
0.011816
0.010095
0.012547
0.007307
0.009379
0.007281
0.010495
0.011073
更新之后的电压幅值为:
0.996360 1.025882 1.012288 1.031051 1.011596 1.027341 1.025000 1.025000 1.040000
相角为:
-0.07030 -0.03924 -0.06520 0.033231 0.011350 0.063670 0.161372 0.080713 0.000000
-----------------------------------------------------------------------------------------------------------------------------------
迭代次数:3
功率不平衡量为:
0.000539
0.000159
0.000522
-0.00067
0.000757
-0.00119
-0.00045
-0.00018
-0.00120
-0.00043
-0.00062
-0.00068
-0.00043
-0.00214
jacobi矩阵为:
-17.7110 11.81205 0.000000 0.000000 0.000000 5.898982 0.000000 0.000000 -1.28370 1.763095 -0.00000 -0.00000 -0.00000 2.021687
11.89872 -41.2436 10.96391 0.000000 0.000000 0.000000 0.000000 0.000000 1.026408 -3.48064 1.732883 -0.00000 -0.00000 -0.00000
-0.00000 10.85918 -16.5320 5.672834 0.000000 0.000000 0.000000 0.000000 -0.00000 2.299646 -2.40340 1.904805 -0.00000 -0.00000
-0.00000 -0.00000 5.935842 -34.0350 10.20710 0.000000 0.000000 17.89215 -0.00000 -0.00000 0.758361 -2.59667 0.986777 -0.00000
-0.00000 -0.00000 -0.00000 10.15415 -23.4438 13.28969 0.000000 0.000000 -0.00000 -0.00000 -0.00000 1.432382 -3.69272 4.261856
6.223731 -0.00000 -0.00000 -0.00000 13.66215 -36.5208 16.63497 0.000000 0.387795 -0.00000 -0.00000 -0.00000 2.850435 -4.87107
-0.00000 -0.00000 -0.00000 -0.00000 -0.00000 16.63497 -16.6349 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -1.63045
-0.00000 -0.00000 -0.00000 17.89215 -0.00000 -0.00000 0.000000 -17.8921 -0.00000 -0.00000 -0.00000 -0.85018 -0.00000 -0.00000
3.784782 -1.76309 0.000000 0.000000 0.000000 -2.02168 0.000000 0.000000 -16.7134 11.81205 0.000000 0.000000 0.000000 5.898982
-1.02640 3.480961 -1.73288 0.000000 0.000000 0.000000 0.000000 0.000000 11.89872 -41.2444 10.96391 0.000000 0.000000 0.000000
0.000000 -2.29964 4.204452 -1.90480 0.000000 0.000000 0.000000 0.000000 -0.00000 10.85918 -15.9332 5.672834 0.000000 0.000000
0.000000 0.000000 -0.75836 2.595327 -0.98677 0.000000 0.000000 -0.85018 -0.00000 -0.00000 5.935842 -34.0364 10.20710 0.000000
0.000000 0.000000 0.000000 -1.43238 5.694238 -4.26185 0.000000 0.000000 -0.00000 -0.00000 -0.00000 10.15415 -22.7447 13.28969
-0.38779 0.000000 0.000000 0.000000 -2.85043 4.868688 -1.63045 0.000000 6.223731 -0.00000 -0.00000 -0.00000 13.66215 -36.5251
方程组求解结果:
0.000061
0.000028
0.000047
0.000162
0.000151
0.000182
0.000196
0.000169
0.000176
0.000092
0.000121
0.000083
0.000133
0.000141
更新之后的电压幅值为:
0.996184 1.025787 1.012166 1.030965 1.011461 1.027196 1.025000 1.025000 1.040000
相角为:
-0.07036 -0.03926 -0.06525 0.033068 0.011198 0.063488 0.161176 0.080544 0.000000
-----------------------------------------------------------------------------------------------------------------------------------
迭代次数:4
功率不平衡量为:
0.000000
0.000000
0.000000
-0.00000
0.000000
-0.00000
-0.00000
-0.00000
-0.00000
-0.00000
-0.00000
-0.00000
-0.00000
-0.00000
迭代完成,更新次数为:3
最大误差为:0.00000
7.总结
- (1)python的function的参数传递挺好用的,可以学一学;
- (2)numpy很强大(让matlab情何以堪~);
- (3)感觉可以专门编一个类似于matpower的python包,期待~;
- (4)整体来说并不难,但是稍微有些细节需要琢磨一下;
注:也可以下载
代码