flowcompute.bas

PQ分解法

Attribute VB_Name = "FlowCompute"

Option Explicit
Public xzjdmc As String
Public xzjdlx As String
Public SQL_name() As String
Public SQL_jdwz() As String
Public SQL_node_jdwz As String
Public I_SQL  As Integer
Public DieDai_CiShu_Hope As Integer
Public DieDai_JingDu As Double '迭代精度
Public Number_Pingheng_Node As Integer
Public ZFNumber_Fdjys As Integer
'*********移植过来的******************
Public Number_ShortNode As Integer '短路节点编号
Public Name_ShortNode As String '网络中节点名称
Public Number_Fuzl As Integer '父子路的编号
Public ZDn_djys_s() As Double '正序导纳对角元素实部
Public ZDn_djys_x() As Double '正序导纳对角元素虚部
Public ZDn_fdjys_s() As Double '正序导纳非对角元素实部
Public ZDn_fdjys_x() As Double '正序导纳非对角元素虚部
Public ZDn_fdjys_hh() As Integer '正序导纳非对角元素行号
Public ZDn_fdjys_lh() As Integer '正序导纳非对角元素列号
Public FDn_djys_s() As Double '负序
Public FDn_djys_x() As Double
Public FDn_fdjys_s() As Double
Public FDn_fdjys_x() As Double
Public FDn_fdjys_hh() As Integer '
Public FDn_fdjys_lh() As Integer '
Public LDn_djys_s() As Double '零序
Public LDn_djys_x() As Double
Public LDn_fdjys_s() As Double
Public LDn_fdjys_x() As Double
Public LDn_fdjys_hh() As Integer '
Public LDn_fdjys_lh() As Integer '
Public ZZk_ys_s() As Double  '正序阻抗
Public ZZk_ys_x() As Double
Public fzk_ys_s() As Double '负序阻抗
Public fzk_ys_x() As Double
Public LZk_ys_s() As Double '零序阻抗
Public LZk_ys_x() As Double
Public Sb_jzgl As Double '基准功率
Public LNumber_Fdjys As Integer 'LING序导纳非对角元素的个数
Public FlagCl As String
Public E() As Double '节点电压实部
Public F() As Double '节点电压虚部
Public ES() As String '节点类型
Public EF() As Double   'ef的记录值
Public Ps() As Double  'p初值
Public Qs() As Double 'qv初值
Public vv() As Double
Public qq() As Double
Public subpi() As Double
Public subqi() As Double
Public kp As Integer
Public kq As Integer
Public subpiflag As Boolean
Public subqiflag As Boolean
Public pqqd As Double
Public pqjd As Double

'*********移植过来的******************
Public Sub Flow_Compute()
Dim rstjdxx As ADODB.Recordset
Dim rstjdxx_flow As ADODB.Recordset
Dim rstzlxx As ADODB.Recordset
Dim rstzlxx_flow As ADODB.Recordset
Dim rstsbtz_fh As ADODB.Recordset
Dim rstsbtz_fDJ As ADODB.Recordset
Dim No_Xunhuan As Integer
Dim A_Xunhuan As Integer
Dim B_Xunhuan As Integer
Dim C_Xunhuan As Integer
Dim D_Xunhuan As Integer
Dim A_Daihuan As Double
Dim B_Daihuan As Double
Dim C_Daihuan As Double
Dim D_Daihuan As Double
Dim E_Daihuan As Double
Dim F_Daihuan As Double
Dim Chongdian_s As Double
Dim Zl_I As Double
Dim Zno_fdj As Integer
Dim Flag_Record_Fsjys As String
Dim Two_SameRecord  As Integer
Dim I_Record_ZFDNZk As Integer
Dim DieDai_Flag As Boolean '已经达到 精度
Dim DieDai_CiShu As Integer
Dim Str_Msgbox As Integer
Dim Flag_Msgboxretry  As Boolean
Dim Jpp() As Double '雅可比矩阵
Dim Ww() As Double '不平衡量
Dim Sjdy As Double '实际电压
Dim Dybyz As Double '电压标幺值
Dim Dyjiaodu As Double '电压角度
Dim P_Load As Double '统计负荷的有功
Dim Q_Load As Double '统计负荷的无功
Dim i As Integer
For i = 1 To Node_Number
    If ES(i) = "0" Then
        E(i) = 1
        F(i) = 0
    End If
    If ES(i) = "1" Then
        F(i) = 0
        Qs(i) = E(i) * E(i)
    End If
    If ES(i) = "2" Then
        F(i) = 0
    End If
Next i
DieDai_CiShu_Hope = 10 '初始值用5  （暂时）
DieDai_JingDu = 0.00001 '??????????????????????????????????????
ExitSub_flag = False
DieDai_CiShu = 0
DieDai_Flag = False
Do Until DieDai_CiShu > DieDai_CiShu_Hope
    DieDai_CiShu = DieDai_CiShu + 1 '迭代的次数加一
    ReDim Ww(1 To 2 * (Node_Number - 1)) As Double '计算不平衡量
    For No_Xunhuan = 1 To Node_Number - 1
        If ES(No_Xunhuan) = "0" Then 'PQ节点
            For A_Xunhuan = 1 To Node_Number
                Ww(2 * No_Xunhuan - 1) = Ww(2 * No_Xunhuan - 1) - (-E(No_Xunhuan) * (gg(No_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - bb(No_Xunhuan, A_Xunhuan) * F(A_Xunhuan)) - F(No_Xunhuan) * (gg(No_Xunhuan, A_Xunhuan) * F(A_Xunhuan) + bb(No_Xunhuan, A_Xunhuan) * E(A_Xunhuan)))
                Ww(2 * No_Xunhuan) = Ww(2 * No_Xunhuan) + (F(No_Xunhuan) * (gg(No_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - bb(No_Xunhuan, A_Xunhuan) * F(A_Xunhuan)) - E(No_Xunhuan) * (gg(No_Xunhuan, A_Xunhuan) * F(A_Xunhuan) + bb(No_Xunhuan, A_Xunhuan) * E(A_Xunhuan)))
            Next A_Xunhuan
            Ww(2 * No_Xunhuan - 1) = Ps(No_Xunhuan) - Ww(2 * No_Xunhuan - 1)
            Ww(2 * No_Xunhuan) = Qs(No_Xunhuan) - Ww(2 * No_Xunhuan)
        End If
        If ES(No_Xunhuan) = "1" Then 'PV节点
            For A_Xunhuan = 1 To Node_Number
                    Ww(2 * No_Xunhuan - 1) = Ww(2 * No_Xunhuan - 1) - (-E(No_Xunhuan) * (gg(No_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - bb(No_Xunhuan, A_Xunhuan) * F(A_Xunhuan)) - F(No_Xunhuan) * (gg(No_Xunhuan, A_Xunhuan) * F(A_Xunhuan) + bb(No_Xunhuan, A_Xunhuan) * E(A_Xunhuan)))
            Next A_Xunhuan
            Ww(2 * No_Xunhuan - 1) = Ps(No_Xunhuan) - Ww(2 * No_Xunhuan - 1)
            Ww(2 * No_Xunhuan) = Qs(No_Xunhuan) - E(No_Xunhuan) * E(No_Xunhuan) - F(No_Xunhuan) * F(No_Xunhuan)
        End If
    Next No_Xunhuan
    DieDai_Flag = True
    ZDWC = 0
    For A_Xunhuan = 1 To 2 * (Node_Number - 1) '检查迭代的精度
        If Abs(Ww(A_Xunhuan)) > DieDai_JingDu Then '
            DieDai_Flag = False
            If Abs(Ww(A_Xunhuan)) > Abs(ZDWC) Then '最大误差找出来
                ZDWC = Ww(A_Xunhuan)
            End If
        End If
    Next A_Xunhuan
    If DieDai_Flag = True Then
        Exit Do
    End If
    ReDim Jpp(1 To 2 * (Node_Number - 1), 1 To 2 * (Node_Number - 1)) As Double '计算雅可比矩阵(结果是相反数)
    For A_Xunhuan = 1 To Node_Number - 1
        For B_Xunhuan = 1 To Node_Number - 1
            If A_Xunhuan = B_Xunhuan Then '对角元素
                Select Case ES(A_Xunhuan)
                    Case "0" 'pq节点
                        Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = -(-gg(A_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - bb(A_Xunhuan, A_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = -(bb(A_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - gg(A_Xunhuan, A_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) = -(bb(A_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - gg(A_Xunhuan, A_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) = -(gg(A_Xunhuan, A_Xunhuan) * E(A_Xunhuan) + bb(A_Xunhuan, A_Xunhuan) * F(A_Xunhuan))
                        For C_Xunhuan = 1 To Node_Number
                            If A_Xunhuan = C_Xunhuan Then
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) + gg(C_Xunhuan, C_Xunhuan) * E(C_Xunhuan) - bb(C_Xunhuan, C_Xunhuan) * F(C_Xunhuan)
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) + gg(C_Xunhuan, C_Xunhuan) * F(C_Xunhuan) + bb(C_Xunhuan, C_Xunhuan) * E(C_Xunhuan)
                                Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) = Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) - gg(C_Xunhuan, C_Xunhuan) * F(C_Xunhuan) - bb(C_Xunhuan, C_Xunhuan) * E(C_Xunhuan)
                                Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) = Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) + gg(C_Xunhuan, C_Xunhuan) * E(C_Xunhuan) - bb(C_Xunhuan, C_Xunhuan) * F(C_Xunhuan)
                            End If
                            If A_Xunhuan <> C_Xunhuan Then
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) + gg(A_Xunhuan, C_Xunhuan) * E(C_Xunhuan) - bb(A_Xunhuan, C_Xunhuan) * F(C_Xunhuan)
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) + gg(A_Xunhuan, C_Xunhuan) * F(C_Xunhuan) + bb(A_Xunhuan, C_Xunhuan) * E(C_Xunhuan)
                                Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) = Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) - gg(A_Xunhuan, C_Xunhuan) * F(C_Xunhuan) - bb(A_Xunhuan, C_Xunhuan) * E(C_Xunhuan)
                                Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) = Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) + gg(A_Xunhuan, C_Xunhuan) * E(C_Xunhuan) - bb(A_Xunhuan, C_Xunhuan) * F(C_Xunhuan)
                            End If
                        Next C_Xunhuan
                    Case "1" 'pv节点
                        Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = -(-gg(A_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - bb(A_Xunhuan, A_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = -(bb(A_Xunhuan, A_Xunhuan) * E(A_Xunhuan) - gg(A_Xunhuan, A_Xunhuan) * F(A_Xunhuan))
                        For C_Xunhuan = 1 To Node_Number
                            If A_Xunhuan = C_Xunhuan Then
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) + gg(C_Xunhuan, C_Xunhuan) * E(C_Xunhuan) - bb(C_Xunhuan, C_Xunhuan) * F(C_Xunhuan)
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) + gg(C_Xunhuan, C_Xunhuan) * F(C_Xunhuan) + bb(C_Xunhuan, C_Xunhuan) * E(C_Xunhuan)
                            End If
                            If A_Xunhuan <> C_Xunhuan Then
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) + gg(A_Xunhuan, C_Xunhuan) * E(C_Xunhuan) - bb(A_Xunhuan, C_Xunhuan) * F(C_Xunhuan)
                                Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) + gg(A_Xunhuan, C_Xunhuan) * F(C_Xunhuan) + bb(A_Xunhuan, C_Xunhuan) * E(C_Xunhuan)
                            End If
                        Next C_Xunhuan
                        Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) = 2 * E(A_Xunhuan)
                        Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) = 2 * F(A_Xunhuan)
                End Select
            End If
            
            If A_Xunhuan <> B_Xunhuan Then '非对角元素
                Select Case ES(A_Xunhuan)
                    Case "0" 'pq节点
                        Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = -(-gg(A_Xunhuan, B_Xunhuan) * E(A_Xunhuan) - bb(A_Xunhuan, B_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = -(bb(A_Xunhuan, B_Xunhuan) * E(A_Xunhuan) - gg(A_Xunhuan, B_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) = -(bb(A_Xunhuan, B_Xunhuan) * E(A_Xunhuan) - gg(A_Xunhuan, B_Xunhuan) * F(A_Xunhuan))
                        Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) = -(gg(A_Xunhuan, B_Xunhuan) * E(A_Xunhuan) + bb(A_Xunhuan, B_Xunhuan) * F(A_Xunhuan))
                    Case "1" 'pv节点
                            Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan - 1) = -(-gg(A_Xunhuan, B_Xunhuan) * E(A_Xunhuan) - bb(A_Xunhuan, B_Xunhuan) * F(A_Xunhuan))
                            Jpp(2 * A_Xunhuan - 1, 2 * B_Xunhuan) = -(bb(A_Xunhuan, B_Xunhuan) * E(A_Xunhuan) - gg(A_Xunhuan, B_Xunhuan) * F(A_Xunhuan))
                            Jpp(2 * A_Xunhuan, 2 * B_Xunhuan - 1) = 0
                            Jpp(2 * A_Xunhuan, 2 * B_Xunhuan) = 0
                End Select
            End If
        Next B_Xunhuan
    Next A_Xunhuan
    Call Line_F(2 * (Node_Number - 1), Jpp(), Ww(), EF())
    If ExitSub_flag = True Then
        Unload FrmSysProc: Exit Sub
    End If
    For No_Xunhuan = 1 To Node_Number - 1
        E(No_Xunhuan) = E(No_Xunhuan) + EF(2 * No_Xunhuan - 1)
        F(No_Xunhuan) = F(No_Xunhuan) + EF(2 * No_Xunhuan)
    Next No_Xunhuan
'    If DieDai_CiShu > 1 Then
'        For i = 1 To Node_Number - 1
''            Dim var1 As Single
''            var1 = Sqr(E(i) ^ 2 + F(i) ^ 2)
'            If (Sqr(E(i) ^ 2 + F(i) ^ 2) < 0.9 Or Sqr(E(i) ^ 2 + F(i) ^ 2) > 1.1) Then
'                E(i) = 1
'                F(i) = 0
'            End If
'        Next i
'    End If
Loop
ReDim unode(1 To Node_Number) As Double
ReDim jdnode(1 To Node_Number) As Double
For i = 1 To Node_Number
    unode(i) = Sqr(E(i) ^ 2 + F(i) ^ 2)
    jdnode(i) = Atn(F(i) / E(i)) * 180 / pi
Next i
If DieDai_Flag = True Then
    Ps(Node_Number) = 0
    Qs(Node_Number) = 0
    For No_Xunhuan = 1 To Node_Number '平衡节点的有功功率和无功功率
        If No_Xunhuan <> Node_Number Then
            Ps(Node_Number) = Ps(Node_Number) + gg(Node_Number, No_Xunhuan) * E(No_Xunhuan) - bb(Node_Number, No_Xunhuan) * F(No_Xunhuan)
            Qs(Node_Number) = Qs(Node_Number) - gg(Node_Number, No_Xunhuan) * F(No_Xunhuan) - bb(Node_Number, No_Xunhuan) * E(No_Xunhuan)
        End If
        If No_Xunhuan = Node_Number Then
            Ps(Node_Number) = Ps(Node_Number) + gg(Node_Number, Node_Number) * E(No_Xunhuan) - bb(Node_Number, Node_Number) * F(No_Xunhuan)
            Qs(Node_Number) = Qs(Node_Number) - gg(Node_Number, Node_Number) * F(No_Xunhuan) - bb(Node_Number, Node_Number) * E(No_Xunhuan)
        End If
    Next No_Xunhuan
    Ps(Node_Number) = E(Node_Number) * Ps(Node_Number) ' * Sb_jzgl
    Qs(Node_Number) = E(Node_Number) * Qs(Node_Number) '* Sb_jzgl
    For i = 1 To Node_Number 'pv节点的有功功率和无功功率
        If ES(i) = "1" Then
            Qs(i) = 0
            For No_Xunhuan = 1 To Node_Number
                If No_Xunhuan <> i Then
                    Qs(i) = Qs(i) - E(i) * (gg(i, No_Xunhuan) * F(No_Xunhuan) + bb(i, No_Xunhuan) * E(No_Xunhuan)) + F(i) * (gg(i, No_Xunhuan) * E(No_Xunhuan) - bb(i, No_Xunhuan) * F(No_Xunhuan))
                End If
                If No_Xunhuan = i Then
                    Qs(i) = Qs(i) - E(i) * (gg(i, i) * F(No_Xunhuan) + bb(i, i) * E(No_Xunhuan)) + F(i) * (gg(i, No_Xunhuan) * E(No_Xunhuan) - bb(i, No_Xunhuan) * F(No_Xunhuan))
                End If
            Next No_Xunhuan
            Qs(i) = Qs(i) '* Sb_jzgl
        End If
    Next i
End If


End Sub
Public Function Line_F(num As Integer, a() As Double, b() As Double, ByRef X() As Double)
    'On Error Resume Next
    Dim i, j, k, n, p, m, q  As Integer
    Dim Daihuan As Double
    ReDim X(1 To num) As Double
    For k = 1 To num '消去运算,消去的结果是改变了上三角的元素，对角和非对角都没有改变
        For n = k + 1 To num '列主消元法
            If Abs(a(n, k)) > Abs(a(k, k)) Then
              For p = k To num
                 Daihuan = a(k, p): a(k, p) = a(n, p): a(n, p) = Daihuan
              Next p
              Daihuan = b(k): b(k) = b(n): b(n) = Daihuan
           End If
       Next n
       If a(k, k) = 0 Then
          MsgBox "出现异常，请检查！", vbOKOnly + vbInformation, "提示信息"
'          Screen.MousePointer = vbArrow: Unload FrmSysProc: ExitSub_flag = True
          Exit Function
       End If
       For i = k + 1 To num
          For m = k + 1 To num
             a(i, m) = a(i, m) - a(k, m) * a(i, k) / a(k, k)
          Next m
          b(i) = b(i) - b(k) * a(i, k) / a(k, k)
       Next i
    Next k
    X(num) = b(num) / a(num, num)
    For i = 1 To num - 1
       Daihuan = 0
       For m = 1 To i
          Daihuan = Daihuan + a(num - i, num - m + 1) * X(num - m + 1)
       Next m
       X(num - i) = (b(num - i) - Daihuan) / a(num - i, num - i)
    Next i
End Function
Public Sub Conndatabase()
    Dim rstjdxx_flow_code As ADODB.Recordset
    ReDim E(1 To Node_Number) As Double
    ReDim F(1 To Node_Number) As Double
    ReDim Ps(1 To Node_Number) As Double
    ReDim Qs(1 To Node_Number) As Double
    ReDim QsPv(1 To Node_Number) As Double
    ReDim ES(1 To Node_Number) As String
    Set rstjdxx_flow_code = New ADODB.Recordset '赋予初值
    rstjdxx_flow_code.Open "select * from jdxx_flow_code  where  code ='" & Operate_Code & "'", cnn
    rstjdxx_flow_code.MoveFirst