67 lines
2.4 KiB
Python
67 lines
2.4 KiB
Python
__author__ = 'dmy'
|
|
|
|
#解潮流方程
|
|
|
|
import numpy
|
|
import scipy.sparse as sparse
|
|
import scipy.io as sio
|
|
|
|
|
|
def pfSolution(pfEquFunc,balance,y,Pd,Qd,Pg,Qg,QGi,QGVm):
|
|
busNum=Pd.size
|
|
Vm=numpy.ones((busNum,1) )
|
|
Vm[QGi-1]=QGVm.reshape(Vm[QGi-1].shape)
|
|
Va=numpy.zeros((busNum,1) )
|
|
#先生成雅克比矩阵
|
|
yNonZeros=y.nonzero()
|
|
maxIte=30
|
|
Loop=0
|
|
convergence=True
|
|
while 1:
|
|
if Loop>maxIte:
|
|
break
|
|
Loop+=1
|
|
ta=Va[yNonZeros[0]]-Va[yNonZeros[1]]-numpy.angle(y[yNonZeros] ).transpose()
|
|
ta=sparse.coo_matrix( (ta.flatten(), yNonZeros),(busNum,busNum) )
|
|
YdotSin=numpy.abs(y).toarray()*numpy.sin(ta.toarray())
|
|
YdotCos=numpy.abs(y).toarray()*numpy.cos(ta.toarray())
|
|
diag_Volt_YdotSin=numpy.dot( numpy.diag(Vm.flatten()),YdotSin)
|
|
diag_Volt_YdotCos=numpy.dot( numpy.diag(Vm.flatten()),YdotCos)
|
|
YdotSinVolt=numpy.dot(YdotSin, Vm)
|
|
YdotCosVolt=numpy.dot(YdotCos, Vm)
|
|
diag_YdotSinVolt_=numpy.diag(YdotSinVolt.flatten())
|
|
diag_YdotCosVolt_=numpy.diag(YdotCosVolt.flatten())
|
|
dPdTyta=numpy.dot(diag_Volt_YdotSin,numpy.diag(Vm.flatten()))-numpy.dot(diag_YdotSinVolt_,numpy.diag(Vm.flatten()) )
|
|
dQdTyta=-numpy.dot(diag_Volt_YdotCos,numpy.diag(Vm.flatten()))+ numpy.dot(diag_YdotCosVolt_,numpy.diag(Vm.flatten()))
|
|
dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos
|
|
dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin
|
|
# 平衡节点相角不变
|
|
dPdTyta[balance-1,:]=0
|
|
dPdTyta[:,balance-1]=0
|
|
dPdTyta=dPdTyta+sparse.coo_matrix( ( [1,] ,([balance-1,],[balance-1,])) ,shape=(busNum,busNum) )
|
|
dQdTyta[:,balance-1]=0
|
|
dPdV[balance-1,:]=0
|
|
# PV节点电压不变
|
|
dQdV[QGi-1,:]=0
|
|
dQdV[:,QGi-1]=0
|
|
dQdV=dQdV+sparse.coo_matrix( ( [1 for x in range(QGi.size)] ,(QGi-1,QGi-1)) ,shape=(busNum,busNum) )
|
|
dQdTyta[QGi-1,:]=0
|
|
dPdV[:,QGi-1]=0
|
|
jacobi=numpy.concatenate( (numpy.concatenate( (dPdV,dPdTyta),axis=1), numpy.concatenate((dQdV,dQdTyta),axis=1) ),axis=0)
|
|
#计算不平衡量
|
|
(dP,dQ)=pfEquFunc(y,Vm,Va,Pd,Qd,Pg,Qg)
|
|
dP[balance-1]=0
|
|
dQ[QGi-1]=0
|
|
dPQ=numpy.concatenate((dP,dQ),axis=0 )
|
|
dX=numpy.linalg.solve(jacobi,-dPQ)
|
|
if numpy.max(numpy.abs(dX))<1e-5:
|
|
break
|
|
Vm+=dX[0:dX.size/2]
|
|
Va+=dX[dX.size/2:]
|
|
if Loop>maxIte:
|
|
convergence=False
|
|
return (Vm,Va,convergence)
|
|
|
|
|
|
|