__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)