module1.bas

PQ分解法

Attribute VB_Name = "Module1"
Public Type data
    voltage As Double
    angle As Double
    p As Double
    q As Double
    nodetype As String
End Type



Public Sub mypq()
Dim i, j, jj, ii As Integer
Dim PQdata() As data
ReDim PQdata(1 To 5) As data
'---------------------原始数据输入
With PQdata(1)
    .nodetype = "2"
    .voltage = 1.06
    .angle = 0
End With
With PQdata(2)
    .nodetype = "0"
    .p = 0.2
    .q = 0.2
    .voltage = 1
    .angle = 0
End With
With PQdata(3)
    .nodetype = "0"
    .p = -0.45
    .q = -0.15
    .voltage = 1
    .angle = 0
End With
With PQdata(4)
    .nodetype = "0"
    .p = -0.4
    .q = -0.05
    .voltage = 1
    .angle = 0
End With
With PQdata(5)
    .nodetype = "0"
    .p = -0.6
    .q = -0.1
    .voltage = 1
    .angle = 0
End With
'****************
Dim G() As Double
Dim b() As Double
ReDim G(1 To 5, 1 To 5) As Double
ReDim b(1 To 5, 1 To 5) As Double
'----------------输入导纳数据
    G(1, 1) = 6.25
    G(1, 2) = -5
    G(1, 3) = -1.25
    G(2, 1) = -5
    G(2, 2) = 10.834
    G(2, 3) = -1.667
    G(2, 4) = -1.667
    G(2, 5) = -2.5
    G(3, 1) = -1.25
    G(3, 2) = -1.667
    G(3, 3) = 12.917
    G(3, 4) = -10
    G(4, 2) = -1.667
    G(4, 3) = -10
    G(4, 4) = 12.917
    G(4, 5) = -1.25
    G(5, 2) = -2.5
    G(5, 4) = -1.25
    G(5, 5) = 3.75
    b(1, 1) = -18.75
    b(1, 2) = 15
    b(1, 3) = 3.75
    b(2, 1) = 15
    b(2, 2) = -32.5
    b(2, 3) = 5
    b(2, 4) = 5
    b(2, 5) = 7.5
    b(3, 1) = 3.75
    b(3, 2) = 5
    b(3, 3) = -38.75
    b(3, 4) = 30
    b(4, 2) = 5
    b(4, 3) = 30
    b(4, 4) = -38.75
    b(4, 5) = 3.75
    b(5, 2) = 7.5
    b(5, 4) = 3.75
    b(5, 5) = -11.25
'***************
Dim B1() As Double   '----B'，不含平衡节点对应的行和列
Dim B2() As Double  '-------B''，不含平衡节点和PV节点对应的行和列
Dim X(1 To 4) As Double
Dim xmax As Double
ReDim B1(1 To 4, 1 To 4)
ReDim B2(1 To 4, 1 To 4)
Dim iii, jjj As Integer
iii = 1: jjj = 1
ii = 1: jj = 1
For i = 1 To 5
    If PQdata(i).nodetype <> "2" Then
        For j = 1 To 5
            If PQdata(j).nodetype <> "2" Then
                B1(iii, jjj) = b(i, j)
                jjj = jjj + 1
            Else
            End If
        Next j
        iii = iii + 1
        jjj = 1
    Else
    End If

    If PQdata(i).nodetype = "0" Then  '----------------形成B'',只含PQ节点
        For j = 1 To 5
            If PQdata(j).nodetype = "0" Then
                B2(ii, jj) = b(i, j)
                jj = jj + 1
            Else
            End If
        Next j
        ii = ii + 1
        jj = 1
    Else
    End If
Next i

Dim B1F() As Double
Dim b2f() As Double
ReDim B1F(1 To 4, 1 To 4) As Double
ReDim b2f(1 To 4, 1 To 4) As Double
Call IM(4, B1, B1F)
Call IM(4, B2, b2f)
'For i = 1 To 4
'    For j = 1 To 4
'        B1F(i, j) = Format(B1F(i, j), "0.000000")
'        b2f(i, j) = Format(b2f(i, j), "0.000000")
'    Next j
'Next i


'+++++++++++++++++++++++++++++++begin iteration---------
Dim k As Integer
k = 0
Dim kp, kq As Integer
Dim detP() As Double
Dim detangle() As Double
Dim detQU() As Double
Dim detU() As Double
Dim Udetangle() As Double
Dim detPU() As Double

  
Dim detQ() As Double
kp = 1: kq = 1
11:
ReDim detP(1 To 5) As Double
    For i = 1 To 5
        With PQdata(i)
            If .nodetype <> "2" Then
                For j = 1 To 5
                    detP(i) = detP(i) + PQdata(j).voltage * (G(i, j) * Cos(.angle - PQdata(j).angle) + b(i, j) * Sin(.angle - PQdata(j).angle))
                Next j
                detP(i) = detP(i) * .voltage
                detP(i) = .p - detP(i)
'                detP(i) = Format(detP(i), "0.000000")
            Else
            End If
        End With
    Next i

ReDim detangle(1 To 4, 0 To 1) As Double
ReDim Udetangle(1 To 4, 0 To 1) As Double
ReDim detPU(1 To 4, 0 To 1) As Double
ii = 1
For i = 1 To 5
    With PQdata(i)
        If .nodetype <> "2" Then
            detPU(ii, 1) = detP(i) / .voltage
            X(ii) = detPU(ii, 1)
            ii = ii + 1
        End If
    End With
Next i
Call find_max(4, X, xmax)
If xmax <= 0.00001 Then
    kp = 0
    If kq = 0 Then
        GoTo 33
    Else
        GoTo 44
    End If
Else
End If
        
Call MatrixMul(4, 4, 1, B1F, detPU, Udetangle)
For i = 1 To 4
'    Udetangle(i, 1) = Format(Udetangle(i, 1), "0.000000")
    With PQdata(i + 1)
        detangle(i, 1) = Udetangle(i, 1) / (-.voltage)
        .angle = .angle + detangle(i, 1)
    End With
Next i
k = k + 1
kq = 1
44:
ReDim detQ(1 To 5) As Double
For i = 1 To 5
    With PQdata(i)
        If .nodetype = "0" Then
            For j = 1 To 5
                detQ(i) = detQ(i) + PQdata(j).voltage * (G(i, j) * Sin(.angle - PQdata(j).angle) - b(i, j) * Cos(.angle - PQdata(j).angle))
            Next j
        End If
        detQ(i) = .voltage * detQ(i)
        detQ(i) = .q - detQ(i)
'        detQ(i) = Format(detQ(i), "0.000000")
    End With
Next i

ReDim detQU(1 To 4, 0 To 1) As Double
ReDim detU(1 To 4, 0 To 1) As Double
ii = 1
For i = 1 To 5
    With PQdata(i)
        If .nodetype = "0" Then
            detQU(ii, 1) = detQ(i) / .voltage
            X(ii) = detQU(ii, 1)
            ii = ii + 1
        End If
    End With
Next i
Call find_max(4, X, xmax)
If xmax <= 0.00001 Then
    kq = 0
    If kp = 0 Then
        GoTo 33
    Else
        GoTo 11
    End If
Else
End If

Call MatrixMul(4, 4, 1, b2f, detQU, detU)
For i = 1 To 4
'    detU(i, 1) = Format(detU(i, 1), "0.000000")
    With PQdata(i + 1)
        .voltage = .voltage + detU(i, 1) / (-1)
    End With
Next i
GoTo 11
33:
For j = 1 To 5
With PQdata(1)
.p = .p + PQdata(j).voltage * (G(1, j) * Cos(.angle - PQdata(j).angle) + b(1, j) * Sin(.angle - PQdata(j).angle))
.q = .q + PQdata(j).voltage * (G(1, j) * Sin(.angle - PQdata(j).angle) - b(1, j) * Cos(.angle - PQdata(j).angle))
End With
Next j
With PQdata(1)
.p = .voltage * .p
.q = .voltage * .q
End With

End Sub
Private Function find_max(n As Integer, X() As Double, ByRef x_Max As Double)
Dim i As Integer
x_Max = X(1)
For i = 1 To n
    If Abs(x_Max) < Abs(X(i)) Then
        x_Max = X(i)
    Else
    End If
Next i
x_Max = Abs(x_Max)
End Function