完成了潮流计算

Signed-off-by: dmy@lab <dmy@lab.lab>

YMatrix.py

@@ -11,15 +11,16 @@ def FormYMatrix(lineBlock,transBlock,busNum):
     x=lineBlock[:,4]
     gb=1/(r+1j*x)
     y=sparse.coo_matrix((-gb,(lineBlock[:,1]-1,lineBlock[:,2]-1)),shape=(busNum,busNum))
-    y=sparse.coo_matrix((-gb,(lineBlock[:,2]-1,lineBlock[:,1]-1)),shape=(busNum,busNum))
+
+    y+=sparse.coo_matrix((-gb,(lineBlock[:,2]-1,lineBlock[:,1]-1)),shape=(busNum,busNum))
     y-=sparse.coo_matrix((y.sum(1).A1,(numpy.arange(int(busNum)),numpy.arange(int(busNum)))),shape=(busNum,busNum))
     #加接地支路
     b2=lineBlock[:,5]
     y+=sparse.coo_matrix((1j*b2,(lineBlock[:,1]-1,lineBlock[:,1]-1)),shape=(busNum,busNum))
     y+=sparse.coo_matrix((1j*b2,(lineBlock[:,2]-1,lineBlock[:,2]-1)),shape=(busNum,busNum))
     #加变压器
-    transX=transBlock[:,3]
-    transR=transBlock[:,4]
+    transR=transBlock[:,3]
+    transX=transBlock[:,4]
     transGB=1/(transR+1j*transX)
     k=transBlock[:,5]
     y+=sparse.coo_matrix((-transGB/k,(transBlock[:,1]-1,transBlock[:,2]-1)),shape=(busNum,busNum))

pfSolver.py

@@ -0,0 +1,66 @@
+__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)
+
+
+

powerflowEquation.py

@@ -0,0 +1,10 @@
+__author__ = 'dmy'
+
+import numpy
+
+def getPowerflowEquation(y,Vm,Va,Pd,Qd,Pg,Qg):
+    cmpV=Vm*numpy.exp(1j*Va)#复数电压
+    S=cmpV*(numpy.dot(y.toarray(),cmpV).conj())
+    P=numpy.real(S)+Pd-Pg
+    Q=numpy.imag(S)+Qd-Qg
+    return (P,Q)

pso.py

@@ -1,9 +1,27 @@
 __author__ = 'dmy'
 
+import scipy.sparse as sparse
+import numpy
+
+from powerflowEquation import *
+from pfSolver import *
 from reader import *
 from YMatrix import *
 
 if __name__=='__main__':
     ieeeData=IEEEData('IEEE4.dat')
+    busNum=ieeeData.busNum
     yMatrix=FormYMatrix(ieeeData.lineBlock,ieeeData.transBlock,ieeeData.busNum)
+    pqBusi=ieeeData.PQbus[:,0]
+    Pd=sparse.coo_matrix((ieeeData.PQbus[:,3],(pqBusi-1,numpy.zeros((busNum,1),dtype=float).flatten() )),shape=(busNum,1))
+    Qd=sparse.coo_matrix((ieeeData.PQbus[:,4],(pqBusi-1,numpy.zeros((busNum,1),dtype=float).flatten() )),shape=(busNum,1))
+    Pg=sparse.coo_matrix((ieeeData.PQbus[:,1],(pqBusi-1,numpy.zeros((busNum,1),dtype=float).flatten() )),shape=(busNum,1))
+    Qg=sparse.coo_matrix((ieeeData.PQbus[:,2],(pqBusi-1,numpy.zeros((busNum,1),dtype=float).flatten() )),shape=(busNum,1))
+    Vm=numpy.ones((busNum,1),dtype=float )
+    Va=numpy.zeros((busNum,1),dtype=float )
+    getPowerflowEquation(yMatrix,Vm,Va,Pd,Qd,Pg,Qg)
+    balance=ieeeData.balance
+    QGi=numpy.array([int(x) for x in ieeeData.QG[:,0] ],dtype=int)
+    QGVm=ieeeData.QG[:,1]
+    pfSolution(getPowerflowEquation,balance,yMatrix,Pd,Qd,Pg,Qg,QGi,QGVm)
     print('Finished.')

reader.py

@@ -34,6 +34,7 @@ class IEEEData:
         self._toPG(array)
         self._toQG(array)
         self.busNum=array[0,1]
+        self.balance=array[2,1]
 
     def _toTrans(self,array):
         zeros=numpy.where(array[:,0]==0)[0]