使用深度学习，对低氧舱滞后时间变化的估计（一）

最近希望通过ansys仿真结合深度学习对低氧舱气体扩散过程中存在的滞后时间进行学习，通过学习模型，通过输入控制系统的控制参数，获得系统的滞后时间，用于未来控制低氧舱决策系统中。
此为该项目的第一部分。~~~~~

一、整体框架

1、实现目标

首先，低氧舱气体输入过程中，由于气体扩散的原因，会产生一定的滞后，并且这个滞后时间很长，对控制系统的精确控制造成一定的困扰。如果可以精确的估计滞后时间，将滞后时间带入到控制系统当中，那么控制精度会极大的提高。

实现目标： 通过学习控制系统中的部分输入参数，和低氧舱的状态参数，和输出的滞后时间（或者扩散系数）来估计低氧舱气体扩散过程中存在的滞后时间。

2、基本思路和框架

低氧舱系统通过控制不同比例的气体输入，不同的气体注入速度。以达到控制舱内氧气浓度的目的。而影响滞后时间的因素有很多。
根据富勒公式：

式中，D为二元气体的扩散系数，MA，MB分别为扩散的两种气体摩尔质量，T为温度，P为压强，va，vb为扩散体积。
各参数单位如下：

扩散体积通过查每一种原子的扩散体积，根据分子式，同过求和可以获得分子的扩散体积。以上式中，所有出现的量均为影响滞后时间的相关量。因此，我推测即使在多元复杂的气体扩散中，这些量依然和滞后时间（扩散系数）存在相关的关系。查阅资料发现多元气体的扩散系数较为复杂，切为非线性，因此没办法计算出具体的公式，所以可以尝试通过机器学习获得这个复杂的模型。

①、首先，通过ansys等仿真软件建立低氧舱的模型。

通过这个模型，获得气体均匀扩散到整个舱内的时间。这个扩散时间就为滞后时间（滞后系数相关量）。设定不同的边界条件（注气孔和出气孔的进出速度，注入气体的混合比例（混合气体的平均分子质量），舱内原有气体的比例（气体的平均分子质量），运行温度，运行压力）。获得对应参数下的滞后时间，这些映射关系作为样本用于机器学习。

②、深度学习的基本思路
深度学习中，输入的数据主要为上述二元气体扩散中的6个指标。输出为扩散系数。在输入前和输入后分别使用其对应的转化公式进行转换（如：输入后扩散系数需要通过低氧舱的体积和公式转换成扩散时间（滞后时间））。

1）、设定一个全连接神经网络用于提高输入数据的维度。输入数据为6个参数，提高数据规模至36。（原来的每一个参数，扩充至由6个参数来表示）
注：扩充后的参数和所有的输入相关，每一个神经网络中每一层的每一个神经元还相同，因此经过全连接层输出的数据为统一的量纲。
2）、将数据进行归一化，softmax转化成概率的形式
3）、将全连接神经网络的数据转成二维，类似张量图6×6规模
4）、将新的学习数据的每一行使用归一化处理。
5）、设定一个卷积神经网络。用于挖掘数据关系。网络将6×6的张量图转换成5×5。
6）、设定一个池化层。用于对特征进行压缩，池化采用平均的办法，池化后张量图转换成4×4。
7）、将处理过的数据转换成一维通过一个全连接神经网络最后获得一个结果。

二、模拟滞后时间数据

1、模拟代码

模拟代码用于生成符合富勒公式的二元气体各项指标和扩散系数。用于测试深度学习的模型。

import pandas as pd
import math
import random
import csv
def read_to_write(data,name):
    data = pd.DataFrame(data)
    # index = [1]则只写入一行
    print(data)
    renovation = data
    print(renovation)
    renovation.to_csv(name,index = False ,sep= ',',header=None)
    return renovation
# 函数用来读取csv文件和写入csv文件
# name是字符串，分别是读取的csv文件名和写入的csv文件名
# data是写入的数据，格式为字典{}，字典内部每一列数据用逗号隔开，每一列数据表示方法为：‘索引’，数组（列数据）

def D_calculate(T,P,Ma,Mb,Va,Vb):
    T = int(T)
    P = int(P)
    Ma = int(Ma)
    Mb = int(Mb)
    va = float(Va)
    vb = float(Vb)
    D1 = 0.0101*T**1.25
    D2 = (1/Ma) + (1/Mb)
    D3 = math.sqrt(D2)
    D4 = D1 * D3
    D5 = P * (va**(1/3) + vb**(1/3))
    D = D4/D5
    print(D)
    return D
# 函数为二元气体扩散系数的计算公式
# D = D_calculate(T,P,Ma,Mb,va,vb)

Diffusion_coefficient = []
T = []
P = []
Ma = []
Mb = []
Va = []
Vb = []
D = []

for i in range(640):
    t = random.randint(273,300)
    # 随机生成一个范围内的整数
    T.append(t)
    p = random.randint(1010000,1020000)
    P.append(p)
    ma = random.randint(28,32)
    Ma.append(ma)
    mb = random.randint(28,32)
    Mb.append(mb)
    va = random.uniform(9.64,12.44)
    # 随机生成范围内的随机数
    Va.append(va)
    vb = random.uniform(9.64,12.44)
    Vb.append(vb)
    d = D_calculate(t,p,ma,mb,va,vb)
    D.append(d)
# 循环通过伪随机数不断地生成符合要求的温度压强等模拟的随机数，然后通过上面的D_calculate计算出D的值
# 将计算的值和伪随机数全部append到一个列表中

name = 'simulation'
data = {'temperature':T,
        'presure':P,
        'air in cabin':Ma,
        'air input':Mb,
        'Molecular diffusion volume 1':Va,
        'Molecular diffusion volume 2':Vb,
        'diffusion coeffcient':D}
read_to_write(data,name)

2、模拟结果

生成一个csv文件在程序所在文件夹下。


这里是每一列所对应的label。

三、深度学习部分

1、全连接神经网络部分

1、处理数据，数据分割

with open(r'simulation') as f:
    reader = csv.reader(f)
    for row in reader:
             data1.append(row[:])

# 模拟的数据在文件里读取后以字符串形式在数组中，因此需要将字符串转换成数字的形式。
for i in range(len(data1)):
    Data = data1[i]
    Data_out = []
    for j in range(len(Data)):
        if j < 4:
            Data_in = int(Data[j])
        else:
            Data_in = float(Data[j])
        Data_out.append(Data_in)
    data.append(Data_out)
data = np.array(data)
"""print(data)"""
# m为样本数量，n为样本长度
m,n = np.shape(data)
print(m,n)
for i in range(m):
    for j in range(n):
       data[i][j] = data[i][j].astype('float64')#是从第三列开始的
# 字段类型转换
data = data.astype('float64')
# 训练样本得输入数据
set1 = data[:,:-1]
# 训练样本得输出数据
set2 = data[:,-1]
"""print(set1 ,"输入和输出数据分隔符",set2 )"""
# 数据分段化处理，切片
# 输出的dataX是每组xback个，常氧数据每个（每行）序列有m个项，dataX中每个是由m个数据构成的向量
def create_interval_dataset(dataset1, dataset2, xback):
    dataX, dataY = [], []
    for i in range(0, len(dataset1) - xback, 32):
        # dataset1 全部的训练数据，每32个为一组
        print(len(dataset1) - xback)
        print(i)
        dataX.append(dataset1[i:i + xback])
        dataY.append(np.reshape(dataset2[i:i + xback], [32, -1]))
        # 将数据转变成32行得维度向量形式，再加到dataY中

    return np.asarray(dataX), np.asarray(dataY)

# dataX，dataY是每32个为一组构成一个向量，全部数组共640组分割成32组构成一个向量
dataX,dataY = create_interval_dataset(set1, set2, 32)
"""print(dataX ,"分割后数据组输入和输出分隔符",dataY)"""
# 至此数据处理和分割完毕
n_set1, n_set2 = np.asarray(set1), np.asarray(set2)
X_train = []
Y_train = []
for i in range(len(dataX)):
    X_tr, Y_tr = dataX[i],dataY[i]
    # print("11", X_tr, Y_tr)
    X_tr = np.reshape(X_tr, [-1,6])
    X_train.append(X_tr)
    Y_tr = np.reshape(Y_tr, [-1,1])
    Y_train.append(Y_tr)
    # print("train:",X_tr , Y_tr )

2、全连接神经网络及其学习结果

step = 5000
# 学习率
learning_rate = 0.00001
# 每次迭代输入网络的数据大小
batch = 32
x = tf1.placeholder(tf.float32,shape = [None,6])
y = tf1.placeholder(tf.float32,shape = [None,1])

w1 = tf1.Variable(tf1.truncated_normal([6,36],stddev = 0.1))
# 训练样本的输入是6个数，而经过全连接层输出的层值为36个，因此是【6，36】
# 如果不设置stddev参数，训练精度最多只能到20%，stddev用于设置正态分布被截断的标准差
w_out = tf1.Variable(tf1.random_normal([36,1]))
b1 = tf1.Variable(tf1.truncated_normal([36]) )
# 为该层节点偏置，值和层输出值的数量一样
b_out = tf1.Variable(tf1.random_normal([1]))
# 定义全连接神经网络，relu为激活函数，使其非线性化

def Fully_neural_network(X):
    layer_1 = tf.nn.relu(tf.add(tf.matmul(X, w1), b1))
    layer_out = tf.nn.relu(tf.add(tf.matmul(layer_1, w_out),  b_out))
    """print("中间层的数据情况:",layer_1 ,"隐藏层到输出的情况：",layer_out )"""
    return layer_out

# 全连接层的输出
net1_out = Fully_neural_network(x)
# Softmax简单的说就是把一个N*1的向量归一化为（0，1）之间的值，采用指数运算
pre = tf.nn.softmax(net1_out)
# tf.nn.softmax_cross_entropy_with_logits_v2（）。计算 softmax(logits) 和 labels 之间的交叉熵
loss = tf1.reduce_mean(tf1.nn.softmax_cross_entropy_with_logits_v2(logits=net1_out,labels = y ))
# Adam算法的优化器，自适应矩估计
# apply_gradients将梯度应用于变量。
# compute_gradients计算 var_list 中变量的 loss 的梯度。
# get_slot返回由 Optimizer 为 var 创建的名为name的slot。
# get_slot_names返回由 Optimizer 创建的 slot 名称的列表。
# variables编码Optimizer当前状态的变量列表。
optimizer = tf1.train.AdamOptimizer(learning_rate=learning_rate)
# minimize通过更新var_list添加操作以最大限度地最小化loss。
train_op = optimizer.minimize(loss)
# equal判断两个元素是否相等，判断方法为元素逐个判断，相等为True，不相等为False
# 得到的值为一个只有1，0的数组
correct_pre = tf1.equal(tf.argmax(pre, 1), tf.argmax(y, 1))
# tf.cast()函数的作用是执行 tensorflow 中张量数据类型转换
# 比如读入的图片如果是int8类型的，一般在要在训练前把图像的数据格式转换为float32。
# 在tensor的某一维度上求值的函数reduce_
# reduce_mean为求平均值
accuracy = tf1.reduce_mean(tf.cast(correct_pre, tf.float64))

# 初始化函数
init = tf1.global_variables_initializer()
print(len(X_train))
with tf1.Session() as sess:
    sess.run(init)
    for i in range(1,step+1):
        m = random.randint(0,18)
        n = random.randint(0,18)
        j = random.randint(0,30)
        test_x = [X_train[m][j]]
        test_y = [Y_train[m][j]]
        batch_x = X_train[n]
        batch_y = Y_train[n]
        sess.run(train_op, feed_dict={x: batch_x, y: batch_y})
        if i % 100 == 0 or i == 1:
            print(batch_x, batch_y)
            l, acc = sess.run([loss, accuracy], feed_dict={x: batch_x,
                                                           y: batch_y})
            print("损失：" , l, acc)
            print("Step " + str(i) + ", Minibatch Loss= " + "{:.4f}".format(l) + ", Training Accuracy= " + "{:.3f}".format(
                acc))
            print("Optimization Finished!")
            # Calculate accuracy for MNIST test images
            print("测试数据：",test_x ,test_y )
            print("Testing Accuracy:", \
                  sess.run(accuracy, feed_dict={x: test_x ,
                                                y: test_y}))

3、该部分的全部代码

import os
import tensorflow as tf
import tensorflow.compat.v1 as tf1
import csv
import numpy as np
import random
import matplotlib.pyplot as plt
tf1.disable_eager_execution()

data1 = []
data = []
# 打开模拟生成的伪随机数文件
# 提取出每一行分别append至data中形成数组，每一项包含输入参数和输出参数
with open(r'simulation') as f:
    reader = csv.reader(f)
    for row in reader:
             data1.append(row[:])

# 模拟的数据在文件里读取后以字符串形式在数组中，因此需要将字符串转换成数字的形式。
for i in range(len(data1)):
    Data = data1[i]
    Data_out = []
    for j in range(len(Data)):
        if j < 4:
            Data_in = int(Data[j])
        else:
            Data_in = float(Data[j])
        Data_out.append(Data_in)
    data.append(Data_out)
data = np.array(data)
"""print(data)"""
# m为样本数量，n为样本长度
m,n = np.shape(data)
print(m,n)
for i in range(m):
    for j in range(n):
       data[i][j] = data[i][j].astype('float64')#是从第三列开始的
# 字段类型转换
data = data.astype('float64')
# 训练样本得输入数据
set1 = data[:,:-1]
# 训练样本得输出数据
set2 = data[:,-1]
"""print(set1 ,"输入和输出数据分隔符",set2 )"""
# 数据分段化处理，切片
# 输出的dataX是每组xback个，常氧数据每个（每行）序列有m个项，dataX中每个是由m个数据构成的向量
def create_interval_dataset(dataset1, dataset2, xback):
    dataX, dataY = [], []
    for i in range(0, len(dataset1) - xback, 32):
        # dataset1 全部的训练数据，每32个为一组
        print(len(dataset1) - xback)
        print(i)
        dataX.append(dataset1[i:i + xback])
        dataY.append(np.reshape(dataset2[i:i + xback], [32, -1]))
        # 将数据转变成32行得维度向量形式，再加到dataY中

    return np.asarray(dataX), np.asarray(dataY)

# dataX，dataY是每32个为一组构成一个向量，全部数组共640组分割成32组构成一个向量
dataX,dataY = create_interval_dataset(set1, set2, 32)
"""print(dataX ,"分割后数据组输入和输出分隔符",dataY)"""
# 至此数据处理和分割完毕
n_set1, n_set2 = np.asarray(set1), np.asarray(set2)
X_train = []
Y_train = []
for i in range(len(dataX)):
    X_tr, Y_tr = dataX[i],dataY[i]
    # print("11", X_tr, Y_tr)
    X_tr = np.reshape(X_tr, [-1,6])
    X_train.append(X_tr)
    Y_tr = np.reshape(Y_tr, [-1,1])
    Y_train.append(Y_tr)
    # print("train:",X_tr , Y_tr )

# 定义训练参数
# 网络迭代次数
step = 5000
# 学习率
learning_rate = 0.00001
# 每次迭代输入网络的数据大小
batch = 32
x = tf1.placeholder(tf.float32,shape = [None,6])
y = tf1.placeholder(tf.float32,shape = [None,1])

w1 = tf1.Variable(tf1.truncated_normal([6,36],stddev = 0.1))
# 训练样本的输入是6个数，而经过全连接层输出的层值为36个，因此是【6，36】
# 如果不设置stddev参数，训练精度最多只能到20%，stddev用于设置正态分布被截断的标准差
w_out = tf1.Variable(tf1.random_normal([36,1]))
b1 = tf1.Variable(tf1.truncated_normal([36]) )
# 为该层节点偏置，值和层输出值的数量一样
b_out = tf1.Variable(tf1.random_normal([1]))
# 定义全连接神经网络，relu为激活函数，使其非线性化

def Fully_neural_network(X):
    layer_1 = tf.nn.relu(tf.add(tf.matmul(X, w1), b1))
    layer_out = tf.nn.relu(tf.add(tf.matmul(layer_1, w_out),  b_out))
    """print("中间层的数据情况:",layer_1 ,"隐藏层到输出的情况：",layer_out )"""
    return layer_out

# 全连接层的输出
net1_out = Fully_neural_network(x)
# Softmax简单的说就是把一个N*1的向量归一化为（0，1）之间的值，采用指数运算
pre = tf.nn.softmax(net1_out)
# tf.nn.softmax_cross_entropy_with_logits_v2（）。计算 softmax(logits) 和 labels 之间的交叉熵
loss = tf1.reduce_mean(tf1.nn.softmax_cross_entropy_with_logits_v2(logits=net1_out,labels = y ))
# Adam算法的优化器，自适应矩估计
# apply_gradients将梯度应用于变量。
# compute_gradients计算 var_list 中变量的 loss 的梯度。
# get_slot返回由 Optimizer 为 var 创建的名为name的slot。
# get_slot_names返回由 Optimizer 创建的 slot 名称的列表。
# variables编码Optimizer当前状态的变量列表。
optimizer = tf1.train.AdamOptimizer(learning_rate=learning_rate)
# minimize通过更新var_list添加操作以最大限度地最小化loss。
train_op = optimizer.minimize(loss)
# equal判断两个元素是否相等，判断方法为元素逐个判断，相等为True，不相等为False
# 得到的值为一个只有1，0的数组
correct_pre = tf1.equal(tf.argmax(pre, 1), tf.argmax(y, 1))
# tf.cast()函数的作用是执行 tensorflow 中张量数据类型转换
# 比如读入的图片如果是int8类型的，一般在要在训练前把图像的数据格式转换为float32。
# 在tensor的某一维度上求值的函数reduce_
# reduce_mean为求平均值
accuracy = tf1.reduce_mean(tf.cast(correct_pre, tf.float64))

# 初始化函数
init = tf1.global_variables_initializer()
print(len(X_train))
with tf1.Session() as sess:
    sess.run(init)
    for i in range(1,step+1):
        m = random.randint(0,18)
        n = random.randint(0,18)
        j = random.randint(0,30)
        test_x = [X_train[m][j]]
        test_y = [Y_train[m][j]]
        batch_x = X_train[n]
        batch_y = Y_train[n]
        sess.run(train_op, feed_dict={x: batch_x, y: batch_y})
        if i % 100 == 0 or i == 1:
            print(batch_x, batch_y)
            l, acc = sess.run([loss, accuracy], feed_dict={x: batch_x,
                                                           y: batch_y})
            print("损失：" , l, acc)
            print("Step " + str(i) + ", Minibatch Loss= " + "{:.4f}".format(l) + ", Training Accuracy= " + "{:.3f}".format(
                acc))
            print("Optimization Finished!")
            # Calculate accuracy for MNIST test images
            print("测试数据：",test_x ,test_y )
            print("Testing Accuracy:", \
                  sess.run(accuracy, feed_dict={x: test_x ,
                                                y: test_y}))