潮流计算(Python语言描述)

单文件

import numpy as np
# 本程序采用牛顿——拉夫逊进行潮流计算
n = 5                        # input("请输入节点数,n=")
nl = 5                       # input("请输入支路数,nl=")
isb = 1                      # input("请输入平衡节点的节点号:isb=")
pr = 0.00001                  # input("请输入误差精度,pr=")
B1 = [[1, 2, 0.03j, 0, 1.05, 0],
      [2, 3, 0.08+0.3j, 0.5j, 1, 0],
      [2, 4, 0.1+0.35j, 0, 1, 0],
      [3, 4, 0.04+0.25j, 0.5j, 1, 0],
      [3, 5, 0.015j, 0, 1.05, 1]]
B2 = [[0, 0, 1.05, 1.05, 0, 1],
      [0, 3.7+1.3j, 1, 0, 0, 2],
      [0, 2+1j, 1, 0, 0, 2],
      [0, 1.6+0.8j, 1, 0, 0, 2],
      [5, 0, 1.05, 1.05, 0, 3]]
X = [[1, 0],
     [2, 0],
     [3, 0],
     [4, 0],
     [5, 0]]
Y = [[0 for i in range(n)] for i in range(n)]
e = np.zeros((1, n))
f = np.zeros((1, n))
V = np.zeros((1, n))
S1 = np.zeros((nl, nl))
S = [0 for i in range(n)]
for i in range(n):
    if X[i][1] != 0:
        Y[X[i][0], X[i][0]] = 1/X[i][1]

for k in range(nl):
    if B1[k][5] == 0:
        k1 = B1[k][0]-1
        k2 = B1[k][1]-1
    else:
        k1 = B1[k][1]-1
        k2 = B1[k][0]-1
    Y[k1][k2] = Y[k1][k2] - 1 / (B1[k][2] * B1[k][4])
    Y[k2][k1] = Y[k1][k2]
    Y[k2][k2] = Y[k2][k2] + 1 / (B1[k][2] * B1[k][4] ** 2) + B1[k][3] / 2
    Y[k1][k1] = Y[k1][k1] + 1 / B1[k][2] + B1[k][3] / 2
G = np.zeros((n, n))
B = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        G[i, j] = Y[i][j].real
        B[i, j] = Y[i][j].imag
for i in range(n):
    e[0, i] = B2[i][2].real
    f[0, i] = B2[i][2].imag
    V[0, i] = B2[i][3]
for i in range(n):
    S[i] = B2[i][0] - B2[i][1]
    B[i, i] = B[i, i] + B2[i][4]
print(S)
S = np.mat(S)
P = S.real
Q = S.imag
print(P)
print(Q)
ICT1 = 0
IT2 = 1
N0 = 2*n
N = N0+1
A = 0
C = [0 for i in range(n)]
D = [0 for i in range(n)]
while IT2 != 0:
    IT2 = 0;
    a = a+1
    for i in range(n):
        if i != isb-1:
            C[i] = 0
            D[i] = 0
            for j1 in range(n):
                C[i] = C[i] + G[i, j1] * e[0, j1] - B[i, j1] * f[0, j1]
                D[i] = D[i] + G[i, j1] * f[0, j1] + B[i, j1] * e[0, j1]
            P1 = C[i] * e[0, i] + f[0, i] * D[i]
            Q1 = -D[i] * e[0, i] + f[0, i] * C[i]
            V2 = e[0, i] * e[0, i] + f[0, i] * f[0, i]
            if B2[i, 5] != 3:
                DP = P[0, i] - P1
                DQ = Q[0, i] - Q1
                for j1 in range(1, n):
                    if j1 != isb-1 and j != i:
                        X1 = -G[i, j1] * e[0, i] - B[i, j1] * f[0, i]
                        X2 = -G[i, j1] * f[0, i] + B[i, j1] * e[0, i]
                        X3 = X2
                        X4 = -X1
                        p = 2*i-1
                        q = 2*j1-1
                        J = 

主函数

from input_data import input_net_args
from NL_iteration import NL_Iteration

# Tips :
# 序号PQ节点在前,PV节点在后
# 输入参数为标幺值

# 1.数据输入
Line_arg = [
#   导线首端    导线末端    串联电阻    串联电抗    并联电导    并联电纳
    [1  , 2  , 0.00037  , 0.00404 , 0 , 0.12264],
    [2  , 3  , 0.00023  , 0.00309 , 0 , 0.16046],
    [2  , 4  , 0.000585 , 0.07445 , 0 , 1.61600],
    [4  , 5  , 0.00006  , 0.00170 , 0 , 0.53379],
    [5  , 3  , 0.00052  , 0.00702 , 0 , 1.46965],
    [5  , 10 , 0.00030  , 0.00459 , 0 , 1.00998],
    [5  , 9  , 0.00068  , 0.00830 , 0 , 1.63571],
    [10 , 14 , 0.00005  , 0.00045 , 0 , 0.09629],
    [10 , 9  , 0.00013  , 0.00148 , 0 , 0.30356],
    [9  , 11 , 0.00243  , 0.00942 , 0 , 0.03080],
    [9  , 13 , 0.00034  , 0.00455 , 0 , 0.40547],
    [9  , 8  , 0.00070  , 0.00885 , 0 , 1.82756],
    [8  , 12 , 0.00025  , 0.00046 , 0 , 0.97515],
    [8  , 7  , 0.00046  , 0.00057 , 0 , 1.04222],
    [7  , 15 , 0.00049  , 0.00511 , 0 , 0.24952],
    [7  , 6  , 0.00046  , 0.00637 , 0 , 1.25707],
    [6  , 16 , 0.00015  , 0.00187 , 0 , 0.09495],
    [6  , 3  , 0.00113  , 0.01521 , 0 , 0.71247],
    # 变压器
    [17 , 14 , 0.00023 , 0.019500 , 0 , 0],
    [18 , 4  , 0.00012 , 0.012199 , 0 , 0],
    [19 , 15 , 0.00060 , 0.044110 , 0 , 0],
    [20 , 16 , 0.00022 , 0.019280 , 0 , 0]
]

Node_args = [
#       节点序号    类型    参数
    [1 , "pq"    , {
    
    "p":3.98600, "q":0.1640} ],
    [2 , "pq"    , {
    
    "P":3.46400, "Q":0.0000} ],
    [3 , "pq"    , {
    
    "P":8.31100, "Q":-0.686} ],
    [4 , "pq"    , {
    
    "P":0.00000, "Q":0.0000} ],
    [5 , "pq"    , {
    
    "P":1.17400, "Q":3.8850} ],
    [6 , "pq"    , {
    
    "p":15.0900, "q":-3.781} ],
    [7 , "pq"    , {
    
    "P":3.18600, "Q":0.7900} ],
    [8 , "pq"    , {
    
    "P":5.72100, "Q":3.6850} ],
    [9 , "pq"    , {
    
    "P":14.9710, "Q":-1.466} ],
    [10, "pq"    , {
    
    "P":9.07100, "Q":0.7590} ],
    [11, "pq"    , {
    
    "P":6.29200, "Q":-3.761} ],
    [12, "pq"    , {
    
    "P":5.47700, "Q":-2.305} ],
    [13, "pq"    , {
    
    "P":12.2850, "Q":1.2770} ],
    [14, "pq"    , {
    
    "P":0, "Q":0} ],
    [15, "pq"    , {
    
    "P":0, "Q":0} ],
    [16, "pq"    , {
    
    "P":0, "Q":0} ],
    [17, "pv"    , {
    
    "P":11.3330, "v":0.9931} ],
    [18, "pv"    , {
    
    "P":55.7030, "v":0.9884} ],
    [19, "pv"    , {
    
    "P":5.98600, "v":0.9815} ],
    [20, "slack" , {
    
    "V":1.0, "Theta":-0.7828} ]
]

Init_val = [
    #       节点序号    参数
    [1, {
    
    "e":1.02, "f":0}],
    [2, {
    
    "e":1, "f":0}],
    [3, {
    
    "e":1, "f":0}],
    [4, {
    
    "e":0.605533, "f":-0.781194}],
    [5, {
    
    "e":0.556615, "f":-0.813259}],
    [6, {
    
    "e":0.63727226, "f":-0.6927696}],
    [7, {
    
    "e":0.561588141, "f":-0.786478}],
    [8, {
    
    "e":0.4674631, "f":-0.8591612}],
    [9, {
    
    "e":0.428363, "f":-0.8940802}],
    [10, {
    
    "e":0.463429, "f":-0.87664366}],
    [11, {
    
    "e":0.4069631, "f":-0.9118022}],
    [12, {
    
    "e":0.4503565, "f":-0.8830939}],
    [13, {
    
    "e":0.3723258, "f":-0.9114882}],
    [14, {
    
    "e":0.4080243, "f":-0.9054081}],
    [15, {
    
    "e":0.5949248, "f":-0.7806451}],
    [16, {
    
    "e":1.02, "f":0}],
    [17, {
    
    "e":0.3723258, "f":-0.9114882}],
    [18, {
    
    "e":0.4080243, "f":-0.9054081}],
    [19, {
    
    "e":0.5949248, "f":-0.7806451}],
    [20, {
    
    "e":1.02, "f":0}]
]

# 数据处理
args = input_net_args(Line_arg, Node_args, Init_val)
args.gen_node_admittance_matrix()
args.gen_node_infos()
args.gen_init_values()
# 迭代

nl = NL_Iteration(args)
nl.start_iteration()

迭代求解过程

import numpy as np
import copy

class NL_Iteration(object):
    def __init__(self, infos):
        self.infos = infos
        self.init_value = self.infos.Init_val
        self._delta_P_PQ = []
        self._delta_Q_PQ = []
        self._delta_P_PV = []
        self._delta_U_PV = []
        self.delta_left = []
        self.J = []
        

    def __calc_delta_val(self):
        self._delta_P_PQ = []
        self._delta_Q_PQ = []
        self._delta_P_PV = []
        self._delta_U_PV = []
        for i, node in enumerate(self.infos.Node_infos):
            if node["node_type"] == "SLACK":
                continue
            e_i = 0
            f_i = 0
            P_i = 0
            Q_i = 0
            for item in self.init_value:
                if item[0] == (i+1):
                    e_i = item[1]["e"]
                    f_i = item[1]["f"]
                    break            
            
            if node["node_type"] == "PQ":
                P_i = node["P"]
                Q_i = node["Q"]

                temp = 0
                for j, init_val in enumerate(self.init_value):
                    temp += e_i * (self.infos.G[i][j] * init_val[1]["e"] - self.infos.B[i][j] * init_val[1]["f"]) + f_i * (self.infos.G[i][j] * init_val[1]["f"] + self.infos.B[i][j] * init_val[1]["e"])
                self._delta_P_PQ.append(P_i - temp)

                temp = 0
                for j, init_val in enumerate(self.init_value):
                    temp += f_i * (self.infos.G[i][j] * init_val[1]["e"] - self.infos.B[i][j] * init_val[1]["f"]) - e_i * (self.infos.G[i][j] * init_val[1]["f"] + self.infos.B[i][j] * init_val[1]["e"])
                self._delta_Q_PQ.append(Q_i - temp)


            if node["node_type"] == "PV":
                P_i = node["P"]

                temp = 0
                for j, init_val in enumerate(self.init_value):
                    temp += e_i * (self.infos.G[i][j] * init_val[1]["e"] - self.infos.B[i][j] * init_val[1]["f"]) + f_i * (self.infos.G[i][j] * init_val[1]["f"] + self.infos.B[i][j] * init_val[1]["e"])
                self._delta_P_PV.append(P_i - temp)

                temp = node["V"]**2 - (self.init_value[i][1]["e"]**2 + self.init_value[i][1]["f"]**2)
                self._delta_U_PV.append(temp)

        # 
        self.delta_left = []
        for i in range(0, len(self._delta_P_PQ)):
            self.delta_left.append(self._delta_P_PQ[i])
            self.delta_left.append(self._delta_Q_PQ[i])
        for i in range(0, len(self._delta_P_PV)):
            self.delta_left.append(self._delta_P_PV[i])
            self.delta_left.append(self._delta_U_PV[i])


        

    def __gen_J_mat(self):
        self.J = []
        for i, node_i in enumerate(self.infos.Node_infos):
            if node_i["node_type"] == "SLACK":
                continue
            e_i, f_i = 0, 0
            for item in self.init_value:
                if item[0] == (i+1):
                    e_i = item[1]["e"]
                    f_i = item[1]["f"]
                    break
            t1, t2 = [], []
            for j, node_j in enumerate(self.infos.Node_infos):
                if node_j["node_type"] == "SLACK":
                    continue

                H_ii, N_ii, J_ii, L_ii, R_ii, S_ii = 0, 0, 0, 0, 0, 0
                H_ij, N_ij, J_ij, L_ij, R_ij, S_ij = 0, 0, 0, 0, 0, 0
                
                if i == j :
                    for k, item in enumerate(self.init_value):
                        H_ii += self.infos.G[i][k] * item[1]["f"] + self.infos.B[i][k] * item[1]["e"]
                        N_ii += self.infos.G[i][k] * item[1]["e"] - self.infos.B[i][k] * item[1]["f"]
                    J_ii +=  N_ii - self.infos.B[i][i] * self.init_value[i][1]["f"] - self.infos.G[i][i] * self.init_value[i][1]["e"]
                    L_ii += -H_ii + self.infos.G[i][i] * self.init_value[i][1]["f"] - self.infos.B[i][i] * self.init_value[i][1]["e"]
                    H_ii += -self.infos.B[i][i] * self.init_value[i][1]["e"] + self.infos.G[i][i] * self.init_value[i][1]["f"]
                    N_ii +=  self.infos.G[i][i] * self.init_value[i][1]["e"] + self.infos.B[i][i] * self.init_value[i][1]["f"]
                    
                    t1.extend([H_ii, N_ii])
                    
                    if node_i["node_type"] == "PV":
                        R_ii = 2 * self.init_value[i][1]["f"]
                        S_ii = 2 * self.init_value[i][1]["e"]
                        t2.extend([R_ii, S_ii])
                    else:
                        t2.extend([J_ii, L_ii])
                    
                else:
                    H_ij = -self.infos.B[i][j] * e_i + self.infos.G[i][j] * f_i
                    N_ij =  self.infos.G[i][j] * e_i + self.infos.B[i][j] * f_i
                    J_ij = -N_ij
                    L_ij =  H_ij
                    R_ij = 0
                    S_ij = 0

                    t1.extend([H_ij, N_ij])

                    if node_i["node_type"] == "PV":
                        t2.extend([R_ij, S_ij])
                    else:
                        t2.extend([J_ij, L_ij])
            
            self.J.append(t1)
            self.J.append(t2)

    def __correction(self, correction_matrix):
        correction_list = correction_matrix.T.getA()[0]
        precision = 10**-3
        for item in correction_list:
            if abs(item) > precision :
                break
            return True


        slack_node_num = 0
        for i, node in enumerate(self.infos.Node_infos):
            if node["node_type"] == "SLACK":
                slack_node_num = i + 1
                break
        temp = copy.deepcopy(self.init_value)
        flag_slack = False
        for i, item in enumerate(temp):
            if item[0] == slack_node_num:
                flag_slack = True
                continue
            if flag_slack :
                item[1]["f"] += correction_list[2* (i - 1)]
                item[1]["e"] += correction_list[2* (i - 1) + 1]
            else:
                item[1]["f"] += correction_list[2* i]
                item[1]["e"] += correction_list[2* i + 1]
        self.init_value = copy.deepcopy(temp)

        return False


    def start_iteration(self):

        stop_flag = False
        iteration_time = 0

        while not stop_flag:
            iteration_time += 1
            self.__calc_delta_val()
            self.__gen_J_mat()
            J = np.matrix(self.J).I
            dt = np.transpose(np.matrix(self.delta_left))
            correction_value = J * dt
            stop_flag = self.__correction(correction_value)

            if iteration_time > 100000 :
                print("不收敛")
                break

        print("迭代完毕:次数", iteration_time, end="\n")
        for item in self.init_value:
            print("节点: ", item[0], "   电压: ", item[1]["e"], " + j", item[1]["f"], end="\n")
        
        

输入数据的生成

class input_net_args(object):

    def __init__(self, arg_line, arg_node, arg_init_val):
        self.Line = arg_line
        self.Node = arg_node
        self._Init_val = arg_init_val
        self.Init_val = []
        self.Node_infos = []
        self._line_admittance = []
        self.G = []
        self.B = []
        self.Order = 0
        self._line_args_conv()

    # 复数倒数    
    def _complx_reciprocal(self, real, imag):
        mod = real*real + imag*imag
        return real/mod, imag/mod
    # 阻抗转导纳
    def impedance2admittance(self, R, X):
        return self._complx_reciprocal(R, X)
    # 导纳转阻抗
    def admittance2impedance(self, G, B):
        return self._complx_reciprocal(G, B)

    def _line_args_conv(self):
        for item in self.Line:
            temp = item
            temp[2], temp[3] = self.impedance2admittance(item[2], item[3])
            self._line_admittance.append(temp)


    # 生成导纳矩阵
    def gen_node_admittance_matrix(self):
        temp = []
        for item in self.Line:
            temp.append(item[0])
            temp.append(item[1])

        # 导纳矩阵阶数
        self.Order = max(temp)

        for i in range(0, self.Order):
            temp_1 = []
            temp_2 = []
            for j in range(0, self.Order):
                if i == j:
                    val_real = 0
                    val_imag = 0
                    for item in self.Line:
                        if ( item[0] == (j + 1) or item[1] == (j + 1) ):
                            val_real += item[2] + item[4]
                            val_imag += item[3] + item[5]
                    temp_1.append(val_real)  
                    temp_2.append(-val_imag)                  
                else:
                    val_real = 0
                    val_imag = 0
                    for item in self.Line:
                        if ( item[0] == (i + 1) and item[1] == (j + 1) ) or ( item[1] == (i + 1) and item[0] == (j + 1) ):
                            val_real += -item[2]
                            val_imag += -item[3]
                    temp_1.append(val_real)
                    temp_2.append(-val_imag)                  
            self.G.append(temp_1)
            self.B.append(temp_2)
        return self.G, self.B
            
    def _gen_node_info(self, node_num, node_type, node_info):
        if (not isinstance(node_num, int)) or (node_num <= 0):
            return "argument node_num is not a vaild value.", False
        if (not isinstance(node_type, str)) or (not node_type.isalpha()) or (not ((node_type.upper() == "PQ") or (node_type.upper() == "PV") or (node_type.upper() == "SLACK"))):
            return "argument node_type is not a vaild value.", False

        node_value = {
    
    "node_num":node_num, "node_type":node_type.upper()}
        
        if node_type.upper() == "PQ":
            for k in node_info:
                if k.upper() == "P":
                    node_value.update({
    
    "P":node_info[k]})
                elif k.upper() == "Q":
                    node_value.update({
    
    "Q":node_info[k]})
          

        if node_type.upper() == "PV":
            for k in node_info:
                if k.upper() == "P":
                    node_value.update({
    
    "P":node_info[k]})
                elif k.upper() == "V":
                    node_value.update({
    
    "V":node_info[k]})

        if node_type.upper() == "SLACK":
            for k in node_info:
                if k.upper() == "V":
                    node_value.update({
    
    "V":node_info[k]})
                elif k.upper() == "THETA":
                    node_value.update({
    
    "THETA":node_info[k]})
        return node_value

    def gen_node_infos(self):
        for i in range(0, len(self.Node)):
            for node_info in self.Node:
                if node_info[0] == (i+1):
                    self.Node_infos.append(self._gen_node_info(node_info[0], node_info[1], node_info[2]))
                    break
    def gen_init_values(self):
        for i in range(0, len(self._Init_val)):
            for val in self._Init_val:
                if val[0] == (i+1):
                    self.Init_val.append(val)
                    break

        



if __name__ == "__main__":
    
    # 1.数据输入
    Line_arg = [
    #   导线首端    导线末端    串联电阻    串联电抗    并联电导    并联电纳
        [1 , 2 , 0.02 , 0.06 , 0 , 0],
        [1 , 3 , 0.08 , 0.24 , 0 , 0],
        [2 , 3 , 0.06 , 0.18 , 0 , 0],
        [2 , 4 , 0.06 , 0.18 , 0 , 0],
        [2 , 5 , 0.04 , 0.12 , 0 , 0],
        [3 , 4 , 0.01 , 0.03 , 0 , 0],
        [4 , 5 , 0.08 , 0.24 , 0 , 0]
    ]

    Node_args = [
#       节点序号    类型    参数
        [1, "slack" , {
    
    "V":0.98, "Theta":0}],
        [2, "pq"    , {
    
    "P":0.22, "Q":0.13} ],
        [4, "pv"    , {
    
    "P":0.22, "V":0.99} ],
        [3, "pq"    , {
    
    "P":0.12, "Q":0.13} ],
        [5, "pq"    , {
    
    "P":-0.22, "Q":0.13} ]
    ]

    Init_val = [
        #       节点序号    参数
        [3, {
    
    "e":1, "f":0}],
        [1, {
    
    "e":1, "f":0}],
        [4, {
    
    "e":1, "f":0}],
        [2, {
    
    "e":1, "f":0}],
        [5, {
    
    "e":1, "f":0}]
    ]

    # 数据处理
    test_args = input_net_args(Line_arg, Node_args, Init_val)
    G, B = test_args.gen_node_admittance_matrix()
    # print(G, "\n", B, "\n\n")    
    test_args.gen_node_infos()
    # print(test_args.Node_infos)
    test_args.gen_init_values()
    # print(test_args.Init_val)

求解结果

在这里插入图片描述

猜你喜欢

转载自blog.csdn.net/weixin_40653652/article/details/109908983