You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
321 lines
13 KiB
321 lines
13 KiB
"""PCA的重构算法"""
|
|
|
|
|
|
import math
|
|
import numpy as np
|
|
|
|
|
|
def min_pos(X):
|
|
X[X <= 0] = np.max(X)
|
|
m = np.min(X)
|
|
re = np.where(X == m)
|
|
min_i = re[0]
|
|
min_j = re[1]
|
|
if m < 0:
|
|
m = 0
|
|
return m, min_i, min_j
|
|
|
|
|
|
def Lars(X, Y, D, DIAG, t, limit_line, LockVariable):
|
|
n, m = X.shape
|
|
beta = np.zeros((1, m))
|
|
beta_temp = np.zeros((m, 1))
|
|
A = []
|
|
N_Co_added = 0
|
|
i = 0
|
|
mse = []
|
|
limit_m = len(LockVariable)
|
|
for k in range(m - limit_m - 1):
|
|
i += 1
|
|
ero = np.array(Y - beta[-1, :].T)
|
|
# c=np.dot(P,DIAG,P.T,ero)
|
|
c = np.dot(D, DIAG).dot(D.T).dot(ero) # 计算相关性
|
|
c[LockVariable] = 0
|
|
C = np.max(np.abs(c))
|
|
mse.append(np.dot(ero.T, D).dot(DIAG).dot(D.T).dot(ero))
|
|
if mse[k] < limit_line:
|
|
break
|
|
elif k == 0:
|
|
addTndex = np.where(abs(c) == C)[-1][0]
|
|
A.append(addTndex) # 活动集
|
|
# 更新正在添加的相应协方差索引的数量
|
|
N_Co_added = N_Co_added + 1
|
|
A_c = list(set(range(0, m)).difference(set(A))) # 非活动集
|
|
s_A = np.diag(np.sign(c[A]))
|
|
num_Zero_Coeff = len(A_c)
|
|
## 计算 X_A， A_A ， u_A ，the inner product vecto
|
|
X_A = np.dot(X[:, A], s_A).reshape(n, -1)
|
|
G_A = np.dot(X_A.T, X_A)
|
|
One_A = np.ones((len(A), 1))
|
|
s = One_A.copy()
|
|
if G_A.shape == ():
|
|
inv_GA = 1 / G_A
|
|
else:
|
|
inv_GA = np.linalg.pinv(G_A)
|
|
# G_a_inv_red_cols = np.sum(inv_GA, 1)
|
|
A_A = 1 / np.sqrt(np.dot(s.T, inv_GA).dot(s))
|
|
w_A = (A_A * inv_GA).dot(s) # w_A: (less then 90%)构成等角的单位向量
|
|
u_A = np.dot(X_A, w_A) # .reshape(n)
|
|
a = X.T.dot(u_A) # inner product vector
|
|
gamma_Test = np.zeros((num_Zero_Coeff, 2))
|
|
# gamma=[]
|
|
# if k == m - 1 - limit_m:
|
|
if N_Co_added == m - 1:
|
|
gamma = C / A_A
|
|
else:
|
|
for j in range(num_Zero_Coeff):
|
|
j_p = A_c[j]
|
|
first_term = (C - c[j_p]) / (A_A - a[j_p])
|
|
second_term = (C + c[j_p]) / (A_A + a[j_p])
|
|
gamma_Test[j, :] = np.array([first_term, second_term]).reshape(1, -1)
|
|
gamma, min_i, min_j = min_pos(gamma_Test)
|
|
# gamma.append(m_s)
|
|
try:
|
|
addTndex = A_c[np.min(min_i)]
|
|
except:
|
|
print(123)
|
|
beta_temp = np.zeros((m, 1))
|
|
beta_temp[A] = beta[k, A].reshape(-1, 1) + np.dot(s_A, gamma * w_A)
|
|
beta = np.vstack((beta, beta_temp.transpose())) # 更新的系数即故障f
|
|
return beta, mse
|
|
|
|
|
|
def recon_fault_diagnosis_r(x, m, limit_line, beta_value, model, spe_recon, m_spe, rbc=None):
|
|
|
|
"""基于重构的故障诊断,重构故障>3
|
|
"""
|
|
k_p = int(model["K"])
|
|
count = 0
|
|
dimension = x.shape[0]
|
|
if rbc is None:
|
|
fault_para = sorted(beta_value, key=lambda item: item[1], reverse=True)
|
|
else:
|
|
fault_para = sorted(rbc, key=lambda item: item[1], reverse=True)
|
|
fault_para_index_temp = [item[0] for item in fault_para]
|
|
fault_para_index = [fault_para_index_temp[: i] for i in range(1, len(beta_value) + 1)][::-1]
|
|
fault_matrix, fault_index_recon = [], []
|
|
for item in fault_para_index:
|
|
ke_si_matrix = np.zeros(shape=(dimension, dimension))
|
|
k_matrix = np.zeros(shape=(dimension, 1))
|
|
for i in range(dimension):
|
|
ke_si = np.zeros(shape=(dimension, 1))
|
|
ke_si_i = np.zeros(shape=(dimension, 1)) # 每一行的fai对应的ke_si
|
|
if i in item:
|
|
ke_si_i[i] = 1 # 计算每一行的ke_si
|
|
fai = m @ ke_si_i @ np.linalg.pinv(ke_si_i.T @ m @ ke_si_i) # 计算每一行的fai
|
|
for j in range(dimension):
|
|
ke_si_j = ke_si.copy() # 每一行每一列的ke_si
|
|
if i == j:
|
|
ke_si_matrix[i, j] = 1 # ke_si矩阵对角线元素全为1
|
|
else:
|
|
if j in item:
|
|
ke_si_j[j] = 1
|
|
ke_si_matrix[i, j] = ke_si_j.T @ fai # 计算矩阵其他元素
|
|
k_matrix[i] = x.T @ fai # 计算k 矩阵
|
|
if np.linalg.det(ke_si_matrix) > 1e-15 and len(item) <= ke_si_matrix.shape[0] - k_p:
|
|
f_matrix = np.linalg.pinv(ke_si_matrix) @ k_matrix # 计算f 矩阵
|
|
else:
|
|
continue
|
|
f_matrix_1 = np.zeros(shape=(dimension, 1))
|
|
for k in item:
|
|
f_matrix_1[k] = f_matrix[k]
|
|
if ((x - f_matrix_1.T) @ m @ (x - f_matrix_1.T).T)[0][0] < limit_line:
|
|
count += 1
|
|
fault_matrix = f_matrix_1
|
|
fault_index_recon = item
|
|
else:
|
|
return fault_matrix, np.array(fault_index_recon)
|
|
return fault_matrix, np.array(fault_index_recon) # 没有找到故障
|
|
|
|
|
|
def recon_fault_diagnosis_r_l(x, m, fault_index: "list"):
|
|
|
|
"""基于重构的故障诊断,重构故障<=3,或者是直接用于fai后的spe重构
|
|
"""
|
|
dimension = x.shape[0]
|
|
item = fault_index.copy()
|
|
ke_si_matrix = np.zeros(shape=(dimension, dimension))
|
|
k_matrix = np.zeros(shape=(dimension, 1))
|
|
for i in range(dimension):
|
|
ke_si = np.zeros(shape=(dimension, 1))
|
|
ke_si_i = np.zeros(shape=(dimension, 1)) # 每一行的fai对应的ke si
|
|
if i in item:
|
|
ke_si_i[i] = 1 # 计算每一行的ke_si
|
|
fai = m @ ke_si_i @ np.linalg.pinv(ke_si_i.T @ m @ ke_si_i) # 计算每一行的fai
|
|
for j in range(dimension):
|
|
ke_si_j = ke_si.copy() # 每一行每一列的ke si
|
|
if i == j:
|
|
ke_si_matrix[i, j] = 1 # ke_si矩阵对角线元素全为1
|
|
else:
|
|
if j in item:
|
|
ke_si_j[j] = 1
|
|
ke_si_matrix[i, j] = ke_si_j.T @ fai # 计算矩阵其他元素
|
|
k_matrix[i] = x.T @ fai # 计算k 矩阵
|
|
f_matrix = np.linalg.pinv(ke_si_matrix) @ k_matrix # 计算f 矩阵
|
|
f_matrix_1 = np.zeros(shape=(dimension, 1))
|
|
for k in fault_index:
|
|
f_matrix_1[k] = f_matrix[k]
|
|
fault_matrix = f_matrix_1
|
|
return fault_matrix, np.array(fault_index) # 没有找到故障
|
|
|
|
|
|
def recon_fault_diagnosis_r_c(x, m, limit_line, beta_value, model, spe_recon, m_spe, lock_variable,
|
|
selectVec, rbc=None):
|
|
|
|
"""基于重构的故障诊断,重构故障>3
|
|
"""
|
|
k_p = int(model["K"])
|
|
count = 0
|
|
dimension = x.shape[0]
|
|
if rbc is None:
|
|
fault_para = sorted(beta_value, key=lambda item: item[1], reverse=True)
|
|
else:
|
|
fault_para = sorted(rbc, key=lambda item: item[1], reverse=True)
|
|
fault_para_index_temp = [item[0] for item in fault_para]
|
|
fault_para_index = [fault_para_index_temp[: i] for i in range(1, len(beta_value) + 1)][::-1]
|
|
fault_matrix, fault_index_recon = [], []
|
|
for item in fault_para_index:
|
|
ke_si_matrix = np.zeros(shape=(dimension, dimension))
|
|
k_matrix = np.zeros(shape=(dimension, 1))
|
|
for i in range(dimension):
|
|
ke_si = np.zeros(shape=(dimension, 1))
|
|
ke_si_i = np.zeros(shape=(dimension, 1)) # 每一行的fai对应的ke_si
|
|
if i in item:
|
|
ke_si_i[i] = 1 # 计算每一行的ke_si
|
|
fai = m @ ke_si_i @ np.linalg.pinv(ke_si_i.T @ m @ ke_si_i) # 计算每一行的fai
|
|
for j in range(dimension):
|
|
ke_si_j = ke_si.copy() # 每一行每一列的ke_si
|
|
if i == j:
|
|
ke_si_matrix[i, j] = 1 # ke_si矩阵对角线元素全为1
|
|
else:
|
|
if j in item:
|
|
ke_si_j[j] = 1
|
|
ke_si_matrix[i, j] = ke_si_j.T @ fai # 计算矩阵其他元素
|
|
k_matrix[i] = x.T @ fai # 计算k 矩阵
|
|
if np.linalg.det(ke_si_matrix) > 1e-15 and len(item) <= ke_si_matrix.shape[0] - k_p:
|
|
f_matrix = np.linalg.pinv(ke_si_matrix) @ k_matrix # 计算f 矩阵
|
|
else:
|
|
# 如果矩阵不可逆则使用贡献图法进行重构
|
|
res = contribution_plot(x, model, lock_variable, selectVec)
|
|
return res, "plot" # 返回故障幅值和故障参数
|
|
f_matrix_1 = np.zeros(shape=(dimension, 1))
|
|
for k in item:
|
|
f_matrix_1[k] = f_matrix[k]
|
|
if ((x - f_matrix_1.T) @ m @ (x - f_matrix_1.T).T)[0][0] < limit_line:
|
|
count += 1
|
|
fault_matrix = f_matrix_1
|
|
fault_index_recon = item
|
|
else:
|
|
return fault_matrix, np.array(fault_index_recon)
|
|
return fault_matrix, np.array(fault_index_recon) # 没有找到故障
|
|
|
|
|
|
def contribution_plot(x, model, lock_variable, select_vec):
|
|
|
|
"""贡献图法直接计算重构值"""
|
|
|
|
final_data = np.dot(x, select_vec).dot(select_vec.T)
|
|
final_data[lock_variable] = x[lock_variable]
|
|
recon_data = np.add(np.multiply(final_data, model["Train_X_std"]), model["Train_X_mean"]) # 重构值
|
|
return recon_data
|
|
|
|
|
|
# 前后项序列
|
|
|
|
def sequence_d(test_data, limit, m, f_m, model, ke_si_matrix):
|
|
"""
|
|
采用前后项序列进行诊断
|
|
:param ke_si_matrix: 重构矩阵
|
|
:param model: 模型信息
|
|
:param test_data: 故障数据
|
|
:param limit: 限值
|
|
:param m: 指标矩阵
|
|
:param f_m: 故障方向矩阵
|
|
:return:
|
|
"""
|
|
f_ds_matrix_a = np.zeros(shape=test_data.shape) # 诊断结果的故障方向矩阵
|
|
st = time.time()
|
|
f_d = np.zeros(shape=test_data.shape)
|
|
for i in range(test_data.shape[0]):
|
|
data = test_data[i]
|
|
if data @ m @ data.T > limit:
|
|
# 超限
|
|
f_index, f_a = sequence(data, limit, m, ke_si_matrix) # 故障方向以及故障幅值
|
|
f_d[i] = f_a
|
|
f_ds_matrix_a[i, f_index] = 1
|
|
et = time.time()
|
|
fd1 = (f_m == 1) & (f_ds_matrix_a == 1) # 诊出率
|
|
fa1 = (f_m == 0) & (f_ds_matrix_a == 1) # 误诊率
|
|
fdr_variable1 = fd1[fd1].shape[0] / f_m[f_m == 1].shape[0]
|
|
far_variable1 = fa1[fa1].shape[0] / f_m[f_m == 0].shape[0]
|
|
recon_data = np.add(np.multiply(test_data - f_d, model["Train_X_std"]), model["Train_X_mean"])
|
|
return {"fdr_variable": fdr_variable1, "far_variable": far_variable1, "time": round(et - st, 2), "recon_data": recon_data.tolist()}
|
|
|
|
|
|
def sequence(data, limit, m, ke_si_matrix):
|
|
"""
|
|
前后项序列
|
|
:param ke_si_matrix: ke_si矩阵
|
|
:param data: 样本
|
|
:param limit: 限值
|
|
:param m: 指标矩阵
|
|
:return:
|
|
"""
|
|
t_c = [] # 故障的目标集合
|
|
c_c = list(range(data.shape[0])) # 候选集合
|
|
flag = False # 记录前项序列是否诊断出
|
|
# 初始化矩阵
|
|
dimension = data.shape[0]
|
|
k_matrix = np.zeros(shape=(dimension, 1))
|
|
f_a = np.zeros(shape=(1, dimension)) # 故障幅值
|
|
for i in range(dimension):
|
|
ke_si = np.zeros(shape=(dimension, 1))
|
|
ke_si[i, 0] = 1
|
|
fai = m @ ke_si @ np.linalg.inv(ke_si.T @ m @ ke_si) # 计算每一行的fai
|
|
k_matrix[i] = data.T @ fai # 计算k 矩阵
|
|
# 前项序列
|
|
while True:
|
|
if len(t_c) == 0:
|
|
index_c = list(product(c_c)) # 第一次计算
|
|
else:
|
|
index_c = generate_index(t_c, c_c) # 生成故障参数的排列
|
|
c = {} # 每个index的限值
|
|
for index in index_c:
|
|
f_r = ke_si_matrix[list(index)][:, list(index)]
|
|
f_matrix = k_matrix[index, 0] @ np.linalg.inv(f_r) # 计算f 矩阵
|
|
f = np.zeros(shape=(1, dimension))
|
|
f[0, index] = f_matrix
|
|
c[index] = ((data - f) @ m @ (data - f).T)[0][0] # 计算每个index的重构限值
|
|
l = min(c.values()) # 计算限值的最小值
|
|
l_index = [k for k in c.keys() if c[k] == l][0] # 计算最小值的k即为此次的index
|
|
t_c = list(l_index) # 更新目标集合
|
|
c_c = list(set(c_c) - set(t_c)) # 更新候选集合
|
|
if l < limit or len(c_c) == 0:
|
|
flag = True
|
|
break
|
|
if flag and len(t_c) > 1:
|
|
# 前项序列诊断出来
|
|
while True:
|
|
index_c = list(combinations(t_c, len(t_c) - 1)) # 生成故障参数的集合
|
|
c = {} # 每个index的限值
|
|
for index in index_c:
|
|
f_r = ke_si_matrix[list(index)][:, list(index)]
|
|
f_matrix = k_matrix[index, 0] @ np.linalg.inv(f_r) # 计算f 矩阵
|
|
f = np.zeros(shape=(1, dimension))
|
|
f[0, index] = f_matrix
|
|
c[index] = ((data - f) @ m @ (data - f).T)[0][0] # 计算每个index的重构限值
|
|
l = min(c.values()) # 计算限值的最小值
|
|
l_index = [k for k in c.keys() if c[k] == l][0] # 计算最小值的k即为此次的index
|
|
if l < limit:
|
|
# 比限值小
|
|
t_c = l_index
|
|
else:
|
|
# 比限值大
|
|
break
|
|
# 计算故障幅值
|
|
f_r = ke_si_matrix[list(t_c)][:, list(t_c)]
|
|
f_matrix = k_matrix[t_c, 0] @ np.linalg.inv(f_r) # 计算f 矩阵
|
|
f = np.zeros(shape=(1, dimension))
|
|
f[0, t_c] = f_matrix
|
|
return t_c, f
|
|
|
|
|