1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
| #-*- coding: utf-8 -*-
# ----------------------------------------------------
# Description: Find Optimal Feature Set algorithm
# Created by: Bendong Tan
# Created time: Saturday, August 04, 2018
# Last Modified: Monday, September 11, 2018
# Wuhan University
# ----------------------------------------------------
##############输入输出说明###############
# 所有矩阵均采用0-1编码
# 节点矩阵 bus_mat
# 线路矩阵 branch_mat
# pmu位置矩阵 pmu_loc
import numpy as np
class PMU_Place(object):
def __init__(self, bus_mat=None, branch_mat=None, load_mat=None, gen_mat=None, pmu_loc=np.array([8])):
if bus_mat == None or branch_mat == None or pmu_loc == None:
# 默认新英格兰39节点系统，输入请按照下述格式构造
# 节点矩阵,节点编号
self.bus_mat = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39])
# 线路矩阵，线路首端和末端编号
self.branch_mat = np.array([[1, 2], [1, 39], [2, 3], [2, 25],
[2, 30], [3, 4], [3, 18], [4, 5],
[4, 14], [5, 6], [5, 8], [6, 7],
[6, 11], [6, 31], [7, 8], [8, 9],
[9, 39], [10, 11], [10, 13], [10, 32],
[12, 11], [12, 13], [13, 14], [14, 15],
[15, 16], [16, 17], [16, 19], [16, 21],
[16, 24], [17, 18], [17, 27], [19, 20],
[19, 33], [20, 34], [21, 22], [22, 23],
[22, 35], [23, 24], [23, 36], [25, 26],
[25, 37], [26, 27], [26, 28], [26, 29],
[28, 29], [29, 38]])
# 负荷矩阵，负荷所在节点编号
self.load_mat = np.array([3, 4, 7, 8, 12, 15, 16, 18, 20, 21, 23, 24, 25, 26, 27, 28, 29, 31, 39])
# 发电机矩阵，发电机所在节点编号
self.gen_mat = np.array([30, 31, 32, 33, 34, 35, 36, 37, 38, 39])
# pmu位置矩阵，PMU所在节点编号
self.pmu_loc = pmu_loc
else:
self.bus_mat = bus_mat
self.branch_mat = branch_mat
self.pmu_loc = pmu_loc
# 构造邻接矩阵、pmu位置数组、零注入节点数组
def Construct_matrix(self):
# 初始化
N = len(self.bus_mat)
link_array = np.zeros((N, N))
pmu_array = np.zeros((N, 1))
ZIB_array = np.ones((N, 1))
# 构造邻接矩阵
branch = self.branch_mat
for i in range(len(branch)):
link_array[branch[i, 0] - 1, branch[i, 1] - 1] = 1
link_array[branch[i, 1] - 1, branch[i, 0] - 1] = 1
# 构造pmu位置数组
pmu_array[self.pmu_loc - 1] = 1 # 减1是因为numpy从0数起
# 构造零注入节点数组,除去发电机以及负荷节点之外的节点
ZIB_array[self.load_mat - 1] = 0
ZIB_array[self.gen_mat - 1] = 0
return link_array, pmu_array, ZIB_array
def Observation(self, link_array, pmu_array, ZIB_array):
# 初始化
N = len(self.bus_mat)
V = np.zeros((N, 1)) # 电压矩阵
I = np.zeros((N, N)) # 电流矩阵，电流已知那么线路功率也是已知的
G = np.zeros((len(self.gen_mat), 1)) # 发电机矩阵
L = np.zeros((len(self.load_mat), 1)) # 负荷矩阵
# 中间变量
A = np.zeros((N, N))
B = np.zeros((N, 1))
# 按照规则进行可观性描述
# 先确定安装有PMU的点
# 安装PMU的节点其节点电压和相连支路电流、相连节点电压可观
for i in range(N):
# if pmu_array[i] == 1:
for j in range(N):
if j != i:
V[i] = pmu_array[i] or (pmu_array[j] and link_array[i, j])
I[i, j] = (pmu_array[i] or pmu_array[j]) and link_array[i, j]
I[j, i] = (pmu_array[i] or pmu_array[j]) and link_array[i, j]
for i in range(100): # 进行迭代直至收敛
# 支路两端电压可观，则该支路电流可观
for i in range(N):
for j in range(N):
if j != i:
if V[i] == 1 and V[j] == 1 and I[i, j] == 0 and link_array[i, j] == 1:
I[i, j] = V[i] and V[j]
I[j, i] = V[i] and V[j]
# 支路电流和另一端电压可知，则该支路另一端电压已知
for i in range(N):
for j in range(N):
if j != i and V[i] == 0:
V[i] = V[i] or (I[i, j] and V[j])
# 零注入节点未配置PMU，且仅有一条相连支路电流未知，该支路电流可观
for i in range(N):
if ZIB_array[i] == 1:
for j in range(N):
if i != j:
for k in range(N):
A[k] = link_array[i, k] and I[i, k]
if A.sum(axis=0)[0] == link_array[i, :].sum(axis=0) - 1 and I[i, j] == 0:
I[i, j] = ZIB_array[i] and link_array[i, j]
I[j, i] = ZIB_array[i] and link_array[i, j]
# 零注入节点未配置PMU，且所有相邻节点电压已知，则该节点电压已知
for i in range(N):
if ZIB_array[i] == 1 and V[i] == 0:
for j in range(N):
if i != j:
for k in range(N):
B[k] = link_array[i, k] and V[k]
if B.sum(axis=0)[0] == link_array[i, :].sum(axis=0) and V[i] == 0:
V[i] = ZIB_array[i]
return V, I, G, L
# 物理量索引转换
def Index(self,V, I,link_array):
N = len(self.bus_mat)
G = [] # 发电机矩阵
L = [] # 负荷矩阵
V_Bus = np.nonzero(V)
I_Line = np.nonzero(I)
# 由于Python从0计数，节点电压矩阵元素值+1
for i in range(len(V_Bus[0])):
V_Bus[0][i] += 1
V_Bus = V_Bus[0].tolist()
for i in range(len(I_Line[0])):
I_Line[0][i] += 1
for i in range(len(I_Line[0])):
I_Line[1][i] += 1
I_Line = list(I_Line)
I_Line_1 = I_Line[0].tolist()
I_Line_2 = I_Line[1].tolist()
I_Line_Result = []
for i in range(len(I_Line_1)):
I_Line_Result.append((I_Line_1[i], I_Line_2[i]))
# 由于支路电流会计算两次，这里删去支路两端节点相同的一半元素，把不存在于原线路矩阵中的编号去除
for i in range(len(I_Line_1)):
if [I_Line_1[i],I_Line_2[i]] not in self.branch_mat.tolist():
I_Line_Result.remove((I_Line_1[i], I_Line_2[i]))
for i in range(len(self.gen_mat)):
if V[self.gen_mat[i]-1]!=0 and sum(I[self.gen_mat[i]-1,:])==sum(link_array[self.gen_mat[i]-1,:]):
G.append(self.gen_mat[i])
for i in range(len(self.load_mat)):
if V[self.load_mat[i]-1]!=0 and sum(I[self.load_mat[i]-1,:])==sum(link_array[self.load_mat[i]-1,:]):
L.append(self.load_mat[i])
return V_Bus,I_Line_Result,np.array(G),np.array(L)
if __name__ == '__main__':
PMU = PMU_Place(bus_mat=None, branch_mat=None, load_mat=None, gen_mat=None, pmu_loc=np.array([3,8,10,16,20,23,25,29]))
link_array, pmu_array, ZIB_array = PMU.Construct_matrix()
V, I, G, L = PMU.Observation(link_array, pmu_array, ZIB_array)
V_Bus, I_Line_Result, G, L=PMU.Index(V, I,link_array)
print('电压可观节点总数为：', sum(sum(V)))
print('电压可观节点为:', V_Bus)
print('发电机可观节点为:', G)
print('负荷可观节点为:', L)
print('电流可观支路总数为：', sum(sum(I)) / 2) # 由于支路电流会计算两次，这里删去支路两端节点相同的一半元素
print('电流可观支路两端节点编号为：', I_Line_Result)
|
