用Python编写潮流程序（牛拉直角坐标法）

我的2021年的开篇作~~，无非是把自己之前用到的写进博客里面。~~
今天写的是牛顿-拉夫逊法的直角坐标方式求解潮流程序，这类代码在m文件里面已经出现很多，但是Python中比较少，所以今天我先占个位置。

首先理解一下牛顿-拉夫逊法，它其实是一种以逼近的思想去求解最终解，最简单的公式是：

在求解潮流计算时，比较好用的是直角坐标系的方法，一般网络可以求解的话，15次内就可以迭代出最终解，其方程是这样子的（这里偷一下懒，直接照搬文库里的书本内容）：



公式就以上几个，用起来的话其实还是容易的，不要把逻辑混乱就可以了，那么接下来讲一下代码部分。

num_Y = size(Y,1)
#Y为网络的节点矩阵
    n=size(Y,1)-1
    if n == 5:
        PVlist = [4]
    else:
        PVlist = [29,30,31,32,33,34,35,36,37]
    #有一个是平衡节点，故减去1,n=5,num_Y=6，这里用PVlist是保存平衡节点的序号，在下面会用别的公式计算
    J=mat(zeros((2*n,2*n)))
    #J矩阵为雅克比矩阵，大小为2*n
    for i in range(n):
        #i为行数
        for j in range(n):
            #j为列数
            if i in PVlist:
                if i==j:
                    J[2*i,j*2] = sum([-(Y[i,k].real*U[k].real-Y[i,k].imag*U[k].imag) for k in range(num_Y)]) - Y[i,i].real*U[i].real-Y[i,i].imag*U[i].imag
                    J[i*2,j*2+1] = sum([-(Y[i,k].real*U[k].imag + Y[i,k].imag*U[k].real) for k in range(num_Y)]) + Y[i,i].imag*U[i].real-Y[i,i].real*U[i].imag
                    J[i*2+1,j*2]=-2*U[i].real
                    J[i*2+1,j*2+1] = -2*U[i].imag
                else:
                    J[i * 2, j * 2] = -(Y[i, j].real * U[i].real + Y[i, j].imag * U[i].imag)
                    J[i * 2, j * 2 + 1] = Y[i, j].imag * U[i].real - Y[i, j].real * U[i].imag
                    J[i*2+1,j*2] = 0
                    J[i*2+1,j*2+1] = 0
            elif i == j:
                J[i*2,j*2] = sum([-(Y[i,k].real*U[k].real-Y[i,k].imag*U[k].imag) for k in range(num_Y)]) - Y[i,i].real*U[i].real-Y[i,i].imag*U[i].imag
                J[i*2,j*2+1] = sum([-(Y[i,k].real*U[k].imag + Y[i,k].imag*U[k].real) for k in range(num_Y)]) + Y[i,i].imag*U[i].real-Y[i,i].real*U[i].imag
                J[i*2+1,j*2] = sum([Y[i,k].real*U[k].imag+Y[i,k].imag*U[k].real for k in range(num_Y)]) + Y[i,i].imag*U[i].real - Y[i,i].real*U[i].imag
                J[i*2+1,j*2+1] = sum([-(Y[i,k].real *U[k].real - Y[i,k].imag*U[k].imag) for k in range(num_Y)]) + Y[i,i].real*U[i].real+Y[i,i].imag*U[i].imag
            else:
                J[i*2,j*2] = -(Y[i,j].real*U[i].real+Y[i,j].imag*U[i].imag)
                J[i*2+1,j*2+1] = -J[i*2,j*2]
                J[i*2,j*2+1] = Y[i,j].imag*U[i].real - Y[i,j].real*U[i].imag
                J[i*2+1,j*2] = J[i*2,j*2+1]
    J_inv = J.I
    # print(J_inv)
    delta_U = -J_inv*delta_w
    #delta_U为节点电压修正量
# print(delta_U)

在求得最后的电压修正量后，得到迭代的电压，又把P,Q向量求出来，如下：

   n=size(mat(P))
    num_Y = size(Y, 1)
    U=fix_U
    if n == 5:
        PVlist = [4]
    else:
        PVlist = [29,30,31,32,33,34,35,36,37]
    delta_w = mat(ones((n*2,1)))
    for i in range(n):
        if i in PVlist:
            delta_w[i*2,0] = P[i] - U[i].real*sum([Y[i,k].real*U[k].real -Y[i,k].imag*U[k].imag for k in range(num_Y)])-U[i].imag*sum([Y[i,k].real*U[k].imag + Y[i,k].imag*U[k].real for k in range(num_Y)])
            delta_w[i*2+1,0] = statU[i]**2 -  abs((U[i])**2)
            #奇数是PU节点的幅值
        else:
            delta_w[i*2,0] = P[i] - U[i].real*sum([Y[i,k].real*U[k].real -Y[i,k].imag*U[k].imag for k in range(num_Y)])-U[i].imag*sum([Y[i,k].real*U[k].imag + Y[i,k].imag*U[k].real for k in range(num_Y)])
            delta_w[i*2+1,0] = Q[i] - U[i].imag*sum([Y[i,k].real*U[k].real - Y[i,k].imag*U[k].imag for k in range(num_Y)]) + U[i].real * sum([Y[i,k].real*U[k].imag+Y[i,k].imag*U[k].real for k in range(num_Y)])
        #delta_w为8*1向量```

在新的一年，希望大家牛气冲天。