tensorflow 堆叠自编码器SDAE对凯斯西储大学轴承数据进行故障分类

之前的一些介绍：tensorflow自编码器+softmax对凯斯西储大学轴承数据进行故障分类
前面研究了一下自编码器，并顺利利用自编码器对轴承数据进行了故障分类，但我们在实际应用时，应用更多的还是堆叠自编码器，而网上关于堆叠自编码器的tensorflow代码很少，我前面几篇博客的代码中虽然也有两个隐含层，但其实也并不是真正意义上的堆叠，所以自己根据前期的论文资料和自己的理解，捣腾出了这一篇Blog，仅供参考。
关于堆叠自编码器的解释，网上也有很多大神讲的很好，这里我就放一张论文里的图片，图片把SDAE的结构解释地很清楚了。这里说一下，为了更好地提取特征，通常我们会在自编码器中加入噪声，形成降噪自编码器（DAE）模型，再堆叠成堆叠降噪自编码器（SDAE）模型，而我的代码里还没有加入噪声，所以算不上完整意义的SDAE，等后期我会把噪声也加上的，那时候再po一份代码上来。
参考论文：基于改进堆叠降噪自编码网络的轴承故障诊断研究_侯文擎

代码

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 19 14:33:56 2020

@author: Administrator
"""

import numpy as np
import scipy.io as scio
import tensorflow as tf
import os


"""
预处理有以下几个步骤：
1.以400个采样点为一个样本，取300个样本（取120000个采样点）
2.取200个作为训练集，100个作为测试集。
3.将六种类别数据堆叠成numpy数组
3.把训练集和测试集样本变形为(1620,20,20,1)，这样就可以喂入卷积神经网络里进行训练了。
4.打乱数据
5.将标签转化为one-hot
"""
#重置图，否则后面跑第二遍进行预测时，会提示找不到变量（因为我们第二遍运行的时候重复定义了）
tf.reset_default_graph()

def yuchuli(data,label):
    a = scio.loadmat(data)
    a=a['X']
    a = a[0:120000]
    a  = a.reshape(300,400)
    np.random.shuffle(a) #打乱函数顺序
    train = a[:200,:]
    test = a[200:,:]
    label_test = np.array([label for i in range(0,100)])
    #将标签填入后，label_test形状为（100，），维度并没有固定，这样我们后面无法天赋到y_true中，因此想要用reshape把它的形状固定下来
    label_test=label_test.reshape([100,1])
    label_train = np.array([label for i in range(0, 200)])
    label_train=label_train.reshape([200,1])
    return train,test ,label_train ,label_test

#叠加
def stackkk(a,b,c,d,e,f,g,h):
    aa = np.vstack((a, e)) #np.vstack():在竖直方向上堆叠
    bb = np.vstack((b, f)) 
    cc = np.vstack((c, g)) #源代码中这里是水平堆叠的，但是由于我们上面规定了label的形状为1列，为了统一，这里也改成垂直堆叠
    dd = np.vstack((d, h))
    return aa,bb,cc,dd

#在读入数据之前，需要对数据进行处理：1.每个mat文件中包含了3-4个列表，只保留想要的哪一个，其他的删掉 2.将保留下来的那一个的名字统一改为X，
#这样是方便预处理，否则就需要读取每个列表的名字，再进行数据读取，会麻烦一点
x_tra0,x_tes0,y_tra0,y_tes0 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Normal\Normal_0.mat',0)
x_tra1,x_tes1,y_tra1,y_tes1 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Inner\IR007_0.mat',1)
x_tra2,x_tes2,y_tra2,y_tes2 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Inner\IR014_0.mat',1)
x_tra3,x_tes3,y_tra3,y_tes3 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Inner\IR021_0.mat',1)
x_tra4,x_tes4,y_tra4,y_tes4 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Outer\OR007@6_0.mat',2)
x_tra5,x_tes5,y_tra5,y_tes5 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Outer\OR014@6_0.mat',2)
x_tra6,x_tes6,y_tra6,y_tes6 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Outer\OR021@6_0.mat',2)
x_tra7,x_tes7,y_tra7,y_tes7 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Roller\B007_0.mat',3)
x_tra8,x_tes8,y_tra8,y_tes8 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Roller\B014_0.mat',3)
x_tra9,x_tes9,y_tra9,y_tes9 = yuchuli(r'C:\Users\Administrator\Desktop\SDAE\6205-12K\Roller\B021_0.mat',3)

#对数据和标签进行分别的组合
tr1,te1,yr1,ye1=stackkk(x_tra0,x_tes0,y_tra0,y_tes0 ,x_tra1,x_tes1,y_tra1,y_tes1)
tr2,te2,yr2,ye2=stackkk(tr1,te1,yr1,ye1,x_tra2,x_tes2,y_tra2,y_tes2)
tr3,te3,yr3,ye3=stackkk(tr2,te2,yr2,ye2,x_tra3,x_tes3,y_tra3,y_tes3)
tr4,te4,yr4,ye4=stackkk(tr3,te3,yr3,ye3,x_tra4,x_tes4,y_tra4,y_tes4)
tr5,te5,yr5,ye5=stackkk(tr4,te4,yr4,ye4,x_tra5,x_tes5,y_tra5,y_tes5)
tr6,te6,yr6,ye6=stackkk(tr5,te5,yr5,ye5,x_tra6,x_tes6,y_tra6,y_tes6)
tr7,te7,yr7,ye7=stackkk(tr6,te6,yr6,ye6,x_tra7,x_tes7,y_tra7,y_tes7)
tr8,te8,yr8,ye8=stackkk(tr7,te7,yr7,ye7,x_tra8,x_tes8,y_tra8,y_tes8)
tr9,te9,yr9,ye9=stackkk(tr8,te8,yr8,ye8,x_tra9,x_tes9,y_tra9,y_tes9)
#print('测试标签:',ye9,'\n','训练标签：',yr9)

y_train = yr9
y_test = ye9
x_train=tr9
x_test=te9
#源代码中这里将tr9和te9都转换为了[2000,20,20,1]的矩阵，是为了方便传入CNN中
#但是我们的SDAE中只需要二维的矩阵，而这一点我们在预处理中就做好了，因此并不需要进行形状的改变
#转化为二维矩阵
#x_train = tr9.reshape(2000,20,20,1) 
#每一个类别数据取120000个数据，组成[300,400]的矩阵，即300个样本，其中前200行为训练样本，共十个类别，就是400*200*10=800000
#x_test = te9.reshape(1000,20,20,1)


state = np.random.get_state()
#获取随机生成器 np.random的状态 ,常与np.random.set_state() 搭配使用。使随机生成器random保持相同的状态（state）
np.random.shuffle(x_train) #打乱顺序
np.random.set_state(state)
np.random.shuffle(y_train)

def to_one_hot(labels,dimension=4):#4是指四个种类
    results = np.zeros((len(labels),dimension)) # len() 方法返回对象（字符、列表、元组等）长度或项目个数
    for i,label in enumerate(labels): #将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标，一般用在 for 循环当中
        results[i,label] = 1 
        #我们传入的y_train是一个200*1的矩阵，值在0~3之间，将这个值看做m，这个代码表示将results的第i行的第m+1（因为索引从0开始计数）个值变为1
        #比方说，y_true的第一行为1，则m=1，results的第一行就应该是第2个值为1，即第一行为[0,1,0,0]
    return results
one_hot_train_labels = to_one_hot(y_train)#将标签变为了[2000，4]的矩阵
one_hot_test_labels = to_one_hot(y_test)
#print(one_hot_train_labels)



#建立模型
x=tf.placeholder(tf.float32,[None,400])#占位符
y_true = tf.placeholder(tf.int32, [None,4])

def add_layer(input_data,input_num,output_num,activation_function=None):
    w=tf.Variable(initial_value=tf.random_normal(shape=[input_num,output_num]))
    b=tf.Variable(initial_value=tf.random_normal(shape=[1,output_num]))
    output=tf.add(tf.matmul(input_data,w),b)
    
    if activation_function:
        output=activation_function(output)
    return output,w,b

#单个的AE模型
def ae_nn(data,note_data,note):
    hidden_layer,hidden_w,hidden_b=add_layer(data,note_data,note,activation_function=tf.nn.sigmoid)
    output_layer,output_w,output_b=add_layer(hidden_layer,note,note_data,activation_function=tf.nn.sigmoid)
    return output_layer,hidden_layer,hidden_w,hidden_b

def loss_math(data,outputdata):
     loss_dae=tf.reduce_mean(tf.square(outputdata-data)) 
     optimizer_dae=tf.train.AdamOptimizer(learning_rate=0.01).minimize(loss_dae)
     return loss_dae,optimizer_dae

#堆叠自编码器结构，特点是权重都用的是之前单个DAE模型训练好的，隐含层的数据直接传递给下一个隐含层，没有输出层，最后的隐含层直接连接softmax分类器
def build_nn(data,note3,hidden_layer_w1,hidden_layer_w2,hidden_layer_w3,hidden_layer_b1,hidden_layer_b2,hidden_layer_b3):
    #编码
    hidden_layer1=tf.nn.sigmoid(tf.add(tf.matmul(data,hidden_layer_w1),hidden_layer_b1))
    hidden_layer2=tf.nn.sigmoid(tf.add(tf.matmul(hidden_layer1,hidden_layer_w2),hidden_layer_b2))
    hidden_layer3=tf.nn.sigmoid(tf.add(tf.matmul(hidden_layer2,hidden_layer_w3),hidden_layer_b3))
    #分类器
    softmax_layer,w,b=add_layer(hidden_layer3,note3,4,activation_function=tf.nn.softmax)
    return softmax_layer

def train_nn(data):
    #初始化种群
    pop=np.array([200,100,50])#三个隐含层的节点数
    #每个个体训练一次SDAE网络，得到输出结果，此时不是用的分类层的输出，而是直接用输出层的输出结果
    output_dae1,hidden_dae1,hidden_w1,hidden_b1=ae_nn(data,400,pop[0])#先假设note为100
    output_dae2,hidden_dae2,hidden_w2,hidden_b2=ae_nn(hidden_dae1,pop[0],pop[1])
    output_dae3,hidden_dae3,hidden_w3,hidden_b3=ae_nn(hidden_dae2,pop[1],pop[2])
    #利用输出结果计算损失值大小
    loss_dae1,optimizer_dae1=loss_math(data,output_dae1)
    loss_dae2,optimizer_dae2=loss_math(hidden_dae1,output_dae2)
    loss_dae3,optimizer_dae3=loss_math(hidden_dae2,output_dae3)

    y_predict=build_nn(data,pop[2],hidden_w1,hidden_w2,hidden_w3,hidden_b1,hidden_b2,hidden_b3)
    #softmax分类层的交叉熵函数
    cross=tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=y_true,logits=y_predict))
    
    optimizer=tf.train.AdamOptimizer(learning_rate=0.01).minimize(cross)
    #分类结果对比
    correct_prediction = tf.equal(tf.argmax(y_predict, 1), tf.argmax(y_true, 1))
    #计算预测准确率，它们都是Tensor
    accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
    saver=tf.train.Saver()
    
    with tf.Session() as sess:
        sess.run(tf.global_variables_initializer())
        
        #得到每个DAE模型的参数
        output_1,hidden_1,hidden_w1,hidden_b1,loss_1,optimizer_1=sess.run([output_dae1,hidden_dae1,hidden_w1,hidden_b1,loss_dae1,optimizer_dae1],feed_dict={x:x_train,y_true:one_hot_train_labels})
        output_2,hidden_2,hidden_w2,hidden_b2,loss_2,optimizer_2=sess.run([output_dae2,hidden_dae2,hidden_w2,hidden_b2,loss_dae2,optimizer_dae2],feed_dict={x:x_train,y_true:one_hot_train_labels})
        output_3,hidden_3,hidden_w3,hidden_b3,loss_3,optimizer_3=sess.run([output_dae3,hidden_dae3,hidden_w3,hidden_b3,loss_dae3,optimizer_dae3],feed_dict={x:x_train,y_true:one_hot_train_labels})
        loss_1,_=sess.run([loss_dae1,optimizer_dae1],feed_dict={x:x_train,y_true:one_hot_train_labels})
        loss_2,_=sess.run([loss_dae2,optimizer_dae2],feed_dict={x:x_train,y_true:one_hot_train_labels})
        loss_3,_=sess.run([loss_dae3,optimizer_dae3],feed_dict={x:x_train,y_true:one_hot_train_labels})
        
        if not os.path.exists("./SDAE_model_save/checkpoint"):
            for i in range(20):
                epoch_softmax_cost=0
                #x_data,y_data=mnist.train.next_batch(batch_size)
                softmax_cost,_=sess.run([cross,optimizer],feed_dict={x:x_train,y_true:one_hot_train_labels})
                epoch_softmax_cost+=softmax_cost
                # 正确的预测结果
                acc_train=sess.run(accuracy,feed_dict={x:x_train,y_true:one_hot_train_labels})
                # 在Session中运行Tensor可以得到Tensor的值
                # 这里是获取最终模型的正确率
                #print(sess.run(accuracy, feed_dict={x: mnist.test.images, y_true: mnist.test.labels}))
                print("Epoch",i,":",epoch_softmax_cost,acc_train)
     
            saver.save(sess,"./SDAE_model_save/model")
        else:
            saver.restore(sess,"./SDAE_model_save/model")
            predict=sess.run(y_predict,feed_dict={x:x_test,y_true:one_hot_test_labels})
            acc_test=sess.run(accuracy,feed_dict={x:x_test,y_true:one_hot_test_labels})
            #predict=np.argmax(predict,1)

            print(predict,acc_test)
        #print("测试集准确率：",)


if __name__=="__main__":
    train_nn(x)
        
  
  

写到有点high，预告一波：下一篇博客要么是更新加了噪声的堆叠自编码器，要么是用遗传算法优化SDAE ٩(๑´0`๑)۶