coal_materials/CBA_vad_fcad.ipynb

In [1]:

import os
os.environ['CUDA_DEVICE_ORDER'] = 'PCB_BUS_ID'
os.environ['CUDA_VISIBLE_DEVICES'] = '0, 1'

In [2]:

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
#新增加的两行
from pylab import mpl
# 设置显示中文字体
mpl.rcParams["font.sans-serif"] = ["SimHei"]

mpl.rcParams["axes.unicode_minus"] = False

In [3]:

data_0102 = pd.read_excel('./data/20240102/20240102.xlsx', header=[0,1,2])
data_0102

Out[3]:

	Unnamed: 0_level_0	氢	碳	氮	氧	弹筒发热量	挥发分	固定炭
	化验编号	Had	Cad	Nad	Oad	Qb,ad	Vad	Fcad
	Unnamed: 0_level_2	(%)	(%)	(%)	(%)	MJ/kg	(%)	(%)
0	2720110529	3.93	70.18	0.81	25.079	27.820	32.06	55.68
1	2720096883	3.78	68.93	0.77	26.512	27.404	29.96	54.71
2	2720109084	3.48	69.60	0.76	26.148	27.578	29.31	55.99
3	2720084708	3.47	66.71	0.76	29.055	26.338	28.58	53.87
4	2720062721	3.87	68.78	0.80	26.542	27.280	29.97	54.78
...	...	...	...	...	...	...	...	...
223	2720030490	4.12	68.85	0.97	26.055	27.864	32.94	51.89
224	2720028633	3.97	67.04	0.94	28.043	27.368	31.88	51.38
225	2720028634	4.12	68.42	0.96	26.493	27.886	33.16	52.00
226	2720017683	3.88	67.42	0.94	27.760	26.616	31.65	50.56
227	2720017678	3.81	66.74	0.92	28.530	26.688	31.02	50.82

228 rows × 8 columns

In [4]:

cols = [''.join([y for y in x if 'Unnamed' not in y]) for x in data_0102.columns]
cols

Out[4]:

['化验编号',
 '氢Had(%)',
 '碳Cad(%)',
 '氮Nad(%)',
 '氧Oad(%)',
 '弹筒发热量Qb,adMJ/kg',
 '挥发分Vad(%)',
 '固定炭Fcad(%)']

In [5]:

data_0102.columns = cols

In [6]:

data_0102

Out[6]:

	化验编号	氢Had(%)	碳Cad(%)	氮Nad(%)	氧Oad(%)	弹筒发热量Qb,adMJ/kg	挥发分Vad(%)	固定炭Fcad(%)
0	2720110529	3.93	70.18	0.81	25.079	27.820	32.06	55.68
1	2720096883	3.78	68.93	0.77	26.512	27.404	29.96	54.71
2	2720109084	3.48	69.60	0.76	26.148	27.578	29.31	55.99
3	2720084708	3.47	66.71	0.76	29.055	26.338	28.58	53.87
4	2720062721	3.87	68.78	0.80	26.542	27.280	29.97	54.78
...	...	...	...	...	...	...	...	...
223	2720030490	4.12	68.85	0.97	26.055	27.864	32.94	51.89
224	2720028633	3.97	67.04	0.94	28.043	27.368	31.88	51.38
225	2720028634	4.12	68.42	0.96	26.493	27.886	33.16	52.00
226	2720017683	3.88	67.42	0.94	27.760	26.616	31.65	50.56
227	2720017678	3.81	66.74	0.92	28.530	26.688	31.02	50.82

228 rows × 8 columns

In [7]:

# out_cols = ['挥发分Vad(%)']
# drop_cols = ['化验编号', '固定炭Fcad(%)']
out_cols = ['固定炭Fcad(%)']
drop_cols = ['挥发分Vad(%)', '化验编号']

In [8]:

out_cols

Out[8]:

['固定炭Fcad(%)']

In [9]:

data = data_0102.copy()

In [10]:

train_data = data.dropna(subset=out_cols).fillna(0)

In [11]:

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import tensorflow.keras.backend as K

2024-01-08 18:09:21.597754: I tensorflow/stream_executor/platform/default/dso_loader.cc:53] Successfully opened dynamic library libcudart.so.11.0

In [12]:

class TransformerBlock(layers.Layer):
    def __init__(self, embed_dim, num_heads, ff_dim, name, rate=0.1):
        super().__init__()
        self.att = layers.MultiHeadAttention(num_heads=num_heads, key_dim=embed_dim, name=name)
        self.ffn = keras.Sequential(
            [layers.Dense(ff_dim, activation="relu"), layers.Dense(embed_dim),]
        )
        self.layernorm1 = layers.LayerNormalization(epsilon=1e-6)
        self.layernorm2 = layers.LayerNormalization(epsilon=1e-6)
        self.dropout1 = layers.Dropout(rate)
        self.dropout2 = layers.Dropout(rate)

    def call(self, inputs, training):
        attn_output = self.att(inputs, inputs)
        attn_output = self.dropout1(attn_output, training=training)
        out1 = self.layernorm1(inputs + attn_output)
        ffn_output = self.ffn(out1)
        ffn_output = self.dropout2(ffn_output, training=training)
        return self.layernorm2(out1 + ffn_output)

In [13]:

from tensorflow.keras import Model

In [14]:

from tensorflow.keras.initializers import Constant

In [15]:

# Custom loss layer
class CustomMultiLossLayer(layers.Layer):
    def __init__(self, nb_outputs=2, **kwargs):
        self.nb_outputs = nb_outputs
        self.is_placeholder = True
        super(CustomMultiLossLayer, self).__init__(**kwargs)
        
    def build(self, input_shape=None):
        # initialise log_vars
        self.log_vars = []
        for i in range(self.nb_outputs):
            self.log_vars += [self.add_weight(name='log_var' + str(i), shape=(1,),
                                              initializer=tf.initializers.he_normal(), trainable=True)]
        super(CustomMultiLossLayer, self).build(input_shape)

    def multi_loss(self, ys_true, ys_pred):
        assert len(ys_true) == self.nb_outputs and len(ys_pred) == self.nb_outputs
        loss = 0
        for y_true, y_pred, log_var in zip(ys_true, ys_pred, self.log_vars):
            mse = (y_true - y_pred) ** 2.
            pre = K.exp(-log_var[0])
            loss += tf.abs(tf.reduce_logsumexp(pre * mse + log_var[0], axis=-1))
        return K.mean(loss)

    def call(self, inputs):
        ys_true = inputs[:self.nb_outputs]
        ys_pred = inputs[self.nb_outputs:]
        loss = self.multi_loss(ys_true, ys_pred)
        self.add_loss(loss, inputs=inputs)
        # We won't actually use the output.
        return K.concatenate(inputs, -1)

In [16]:

num_heads, ff_dim = 3, 16

In [17]:

def get_prediction_model():
    inputs = layers.Input(shape=(1,len(feature_cols)), name='input')
    x = layers.Conv1D(filters=64, kernel_size=1, activation='relu')(inputs)
    # x = layers.Dropout(rate=0.1)(x)
    lstm_out = layers.Bidirectional(layers.LSTM(units=64, return_sequences=True))(x)
    lstm_out = layers.Dense(128, activation='relu')(lstm_out)
    # transformer_block = TransformerBlock(128, num_heads, ff_dim, name='first_attn')
    # out = transformer_block(lstm_out)
    # out = layers.GlobalAveragePooling1D()(out)
    out = layers.Dropout(0.1)(lstm_out)
    out = layers.Dense(64, activation='relu')(out)
    bet = layers.Dense(1, activation='sigmoid', name='vad')(out)
    model = Model(inputs=inputs, outputs=bet)
    return model

In [18]:

from tensorflow.python.keras.utils.vis_utils import plot_model

In [19]:

maxs = train_data.max()
mins = train_data.min()
for col in train_data.columns:
    if maxs[col] - mins[col] == 0:
        continue
    train_data[col] = (train_data[col] - mins[col]) / (maxs[col] - mins[col])

In [20]:

train_data.head(2)

Out[20]:

	化验编号	氢Had(%)	碳Cad(%)	氮Nad(%)	氧Oad(%)	弹筒发热量Qb,adMJ/kg	挥发分Vad(%)	固定炭Fcad(%)
0	0.996547	0.773973	0.835414	0.456522	0.171463	0.811249	0.847737	0.828147
1	0.851118	0.671233	0.799943	0.369565	0.210254	0.782038	0.674897	0.794606

In [21]:

# feature_cols = [x for x in train_data.columns if x not in out_cols and '第二次' not in x]
feature_cols = [x for x in train_data.columns if x not in out_cols and x not in drop_cols]
use_cols = feature_cols + out_cols

In [22]:

use_data = train_data.copy()
for col in use_cols:
    use_data[col] = use_data[col].astype('float32')

In [23]:

from sklearn.model_selection import KFold, train_test_split
kf = KFold(n_splits=6, shuffle=True, random_state=42)

In [24]:

from tensorflow.keras import optimizers
from tensorflow.python.keras.utils.vis_utils import plot_model
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score, mean_absolute_percentage_error

In [25]:

from keras.callbacks import ReduceLROnPlateau
reduce_lr = ReduceLROnPlateau(monitor='val_loss', patience=10, mode='auto')

In [26]:

model = get_prediction_model()
model.summary()

2024-01-08 18:09:35.590279: I tensorflow/stream_executor/platform/default/dso_loader.cc:53] Successfully opened dynamic library libcuda.so.1
2024-01-08 18:09:35.656503: E tensorflow/stream_executor/cuda/cuda_driver.cc:328] failed call to cuInit: CUDA_ERROR_INVALID_DEVICE: invalid device ordinal
2024-01-08 18:09:35.656548: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:169] retrieving CUDA diagnostic information for host: zhaojh-yv621
2024-01-08 18:09:35.656557: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:176] hostname: zhaojh-yv621
2024-01-08 18:09:35.656758: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:200] libcuda reported version is: 520.61.5
2024-01-08 18:09:35.656795: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:204] kernel reported version is: 520.61.5
2024-01-08 18:09:35.656802: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:310] kernel version seems to match DSO: 520.61.5
2024-01-08 18:09:35.657280: I tensorflow/core/platform/cpu_feature_guard.cc:142] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 AVX512F FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.

Model: "model"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input (InputLayer)           [(None, 1, 5)]            0         
_________________________________________________________________
conv1d (Conv1D)              (None, 1, 64)             384       
_________________________________________________________________
bidirectional (Bidirectional (None, 1, 128)            66048     
_________________________________________________________________
dense (Dense)                (None, 1, 128)            16512     
_________________________________________________________________
dropout (Dropout)            (None, 1, 128)            0         
_________________________________________________________________
dense_1 (Dense)              (None, 1, 64)             8256      
_________________________________________________________________
vad (Dense)                  (None, 1, 1)              65        
=================================================================
Total params: 91,265
Trainable params: 91,265
Non-trainable params: 0
_________________________________________________________________

In [27]:

def print_eva(y_true, y_pred, tp):
    MSE = mean_squared_error(y_true, y_pred)
    RMSE = np.sqrt(MSE)
    MAE = mean_absolute_error(y_true, y_pred)
    MAPE = mean_absolute_percentage_error(y_true, y_pred)
    R_2 = r2_score(y_true, y_pred)
    print(f"COL: {tp}, MSE: {format(MSE, '.2E')}", end=',')
    print(f'RMSE: {round(RMSE, 3)}', end=',')
    print(f'MAPE: {round(MAPE * 100, 3)} %', end=',')
    print(f'MAE: {round(MAE, 3)}', end=',')
    print(f'R_2: {round(R_2, 3)}')
    return [MSE, RMSE, MAE, MAPE, R_2]

In [28]:

from keras.losses import mean_squared_error

In [30]:

vad_eva_list = list()
fcad_eva_list = list()
train_data = use_data[use_cols].copy()
for (train_index, test_index) in kf.split(train_data):
    train = train_data.loc[train_index]
    valid = train_data.loc[test_index]
    X = np.expand_dims(train[feature_cols].values, axis=1)
    Y = [x for x in train[out_cols].values.T]
    X_valid = np.expand_dims(valid[feature_cols].values, axis=1)
    Y_valid = [x for x in valid[out_cols].values.T]
    prediction_model = get_prediction_model()
    prediction_model.compile(optimizer='adam', loss=mean_squared_error)
    hist = prediction_model.fit(X, Y[0], epochs=120, batch_size=8, verbose=0, 
                               validation_data=(X_valid, Y_valid[0]),
                               callbacks=[reduce_lr]
                               )
    rst = prediction_model.predict(X_valid).squeeze(axis=1)
    pred_rst = pd.DataFrame.from_records(np.asarray(rst), columns=out_cols)
    real_rst = valid[out_cols].copy()
    for col in out_cols:
        pred_rst[col] = pred_rst[col] * (maxs[col] - mins[col]) + mins[col]
        real_rst[col] = real_rst[col] * (maxs[col] - mins[col]) + mins[col]
    y_pred_vad = pred_rst[out_cols].values.reshape(-1,)
    # y_pred_fcad = pred_rst['固定炭Fcad(%)'].values.reshape(-1,)
    y_true_vad = real_rst[out_cols].values.reshape(-1,)
    # y_true_fcad = real_rst['固定炭Fcad(%)'].values.reshape(-1,)
    vad_eva = print_eva(y_true_vad, y_pred_vad, tp='挥发分Vad')
    # fcad_eva = print_eva(y_true_fcad, y_pred_fcad, tp='固定炭Fcad')
    vad_eva_list.append(vad_eva)
    # fcad_eva_list.append(fcad_eva)
    del prediction_model

2024-01-08 18:03:50.956250: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:176] None of the MLIR Optimization Passes are enabled (registered 2)
2024-01-08 18:03:50.974801: I tensorflow/core/platform/profile_utils/cpu_utils.cc:114] CPU Frequency: 2200000000 Hz

COL: 挥发分Vad, MSE: 5.84E-01,RMSE: 0.764,MAPE: 2.111 %,MAE: 0.633,R_2: 0.874
COL: 挥发分Vad, MSE: 1.06E+00,RMSE: 1.028,MAPE: 2.941 %,MAE: 0.869,R_2: 0.872
COL: 挥发分Vad, MSE: 6.70E-01,RMSE: 0.819,MAPE: 2.217 %,MAE: 0.658,R_2: 0.829
COL: 挥发分Vad, MSE: 5.96E-01,RMSE: 0.772,MAPE: 2.07 %,MAE: 0.607,R_2: 0.896
WARNING:tensorflow:5 out of the last 9 calls to <function Model.make_predict_function.<locals>.predict_function at 0x7f6e8d6f8940> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has experimental_relax_shapes=True option that relaxes argument shapes that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.
COL: 挥发分Vad, MSE: 8.56E-01,RMSE: 0.925,MAPE: 2.335 %,MAE: 0.717,R_2: 0.805
WARNING:tensorflow:6 out of the last 11 calls to <function Model.make_predict_function.<locals>.predict_function at 0x7f6e8f6e4160> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has experimental_relax_shapes=True option that relaxes argument shapes that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.
COL: 挥发分Vad, MSE: 7.24E-01,RMSE: 0.851,MAPE: 2.435 %,MAE: 0.713,R_2: 0.851

In [ ]:

out_cols = ['固定炭Fcad(%)']
fcad_eva_list = list()
train_data = use_data[use_cols].copy()
for (train_index, test_index) in kf.split(train_data):
    train = train_data.loc[train_index]
    valid = train_data.loc[test_index]
    X = np.expand_dims(train[feature_cols].values, axis=1)
    Y = [x for x in train[out_cols].values.T]
    X_valid = np.expand_dims(valid[feature_cols].values, axis=1)
    Y_valid = [x for x in valid[out_cols].values.T]
    prediction_model = get_prediction_model()
    prediction_model.compile(optimizer='adam', loss=mean_squared_error)
    hist = prediction_model.fit(X, Y[0], epochs=120, batch_size=8, verbose=0, 
                               validation_data=(X_valid, Y_valid[0]),
                               callbacks=[reduce_lr]
                               )
    rst = prediction_model.predict(X_valid).squeeze(axis=1)
    pred_rst = pd.DataFrame.from_records(np.asarray(rst), columns=out_cols)
    real_rst = valid[out_cols].copy()
    for col in out_cols:
        pred_rst[col] = pred_rst[col] * (maxs[col] - mins[col]) + mins[col]
        real_rst[col] = real_rst[col] * (maxs[col] - mins[col]) + mins[col]
    y_pred = pred_rst[out_cols].values.reshape(-1,)
    y_true = real_rst[out_cols].values.reshape(-1,)
    fcad_eva = print_eva(y_true, y_pred, tp='固定炭Fcad')
    fcad_eva_list.append(fcad_eva)
    del prediction_model

2024-01-08 18:09:42.506363: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:176] None of the MLIR Optimization Passes are enabled (registered 2)
2024-01-08 18:09:42.522809: I tensorflow/core/platform/profile_utils/cpu_utils.cc:114] CPU Frequency: 2200000000 Hz

COL: 固定炭Fcad, MSE: 7.16E-01,RMSE: 0.846,MAPE: 1.337 %,MAE: 0.715,R_2: 0.972
COL: 固定炭Fcad, MSE: 1.04E+00,RMSE: 1.018,MAPE: 1.65 %,MAE: 0.847,R_2: 0.98
COL: 固定炭Fcad, MSE: 8.89E-01,RMSE: 0.943,MAPE: 1.294 %,MAE: 0.724,R_2: 0.936
COL: 固定炭Fcad, MSE: 5.17E-01,RMSE: 0.719,MAPE: 1.066 %,MAE: 0.545,R_2: 0.985

In [32]:

vad_df = pd.DataFrame.from_records(vad_eva_list, columns=['MSE', 'RMSE', 'MAE', 'MAPE', 'R_2'])
vad_df.sort_values(by='R_2').mean()

Out[32]:

MSE     0.747816
RMSE    0.859839
MAE     0.699474
MAPE    0.023513
R_2     0.854338
dtype: float64

In [ ]:

fcad_df = pd.DataFrame.from_records(fcad_eva_list, columns=['MSE', 'RMSE', 'MAE', 'MAPE', 'R_2'])
fcad_df.sort_values(by='R_2').mean()

In [ ]:

train, valid = train_test_split(use_data[use_cols], test_size=0.3, random_state=42, shuffle=True)
valid, test = train_test_split(valid, test_size=0.3, random_state=42, shuffle=True)

In [31]:

prediction_model = get_prediction_model()
trainable_model = get_trainable_model(prediction_model)

In [34]:

X = np.expand_dims(train[feature_cols].values, axis=1)
Y = [x for x in train[out_cols].values.T]
Y_valid = [x for x in valid[out_cols].values.T]

In [ ]:

trainable_model.compile(optimizer='adam', loss=None)
hist = trainable_model.fit([X, Y[0], Y[1]], epochs=120, batch_size=8, verbose=1, 
                           validation_data=[np.expand_dims(valid[feature_cols].values, axis=1), Y_valid[0], Y_valid[1]],
                           callbacks=[reduce_lr]
                           )

In [41]:

rst = prediction_model.predict(np.expand_dims(test[feature_cols], axis=1))
rst

Out[41]:

[array([[0.73740077],
        [0.89292204],
        [0.7599046 ],
        [0.67802393],
        [0.6815233 ],
        [0.88627005],
        [0.6121343 ],
        [0.7072234 ],
        [0.8561135 ],
        [0.52762157],
        [0.8325021 ],
        [0.50241977],
        [0.8242289 ],
        [0.68957335],
        [0.6980361 ],
        [0.82116604],
        [0.8566438 ],
        [0.53687835],
        [0.56832707],
        [0.78476715],
        [0.85638577]], dtype=float32),
 array([[0.68600863],
        [0.78454906],
        [0.8179163 ],
        [0.94351083],
        [0.86383885],
        [0.69705516],
        [0.6913491 ],
        [0.80277354],
        [0.93557894],
        [0.82278305],
        [0.82674253],
        [0.93518937],
        [0.8094449 ],
        [0.9206344 ],
        [0.7747319 ],
        [0.9137207 ],
        [0.9491073 ],
        [0.93225   ],
        [0.6185102 ],
        [0.8867341 ],
        [0.82890105]], dtype=float32)]

In [42]:

[np.exp(K.get_value(log_var[0]))**0.5 for log_var in trainable_model.layers[-1].log_vars]

Out[42]:

[0.9991559102070927, 0.9998196796918477]

In [46]:

real_rst.columns

Out[46]:

Index(['挥发分Vad(%)', '固定炭Fcad(%)'], dtype='object')

In [ ]:

for col in out_cols:
    pred_rst[col] = pred_rst[col] * (maxs[col] - mins[col]) + mins[col]
    real_rst[col] = real_rst[col] * (maxs[col] - mins[col]) + mins[col]

In [47]:

y_pred_vad = pred_rst['挥发分Vad(%)'].values.reshape(-1,)
y_pred_fcad = pred_rst['固定炭Fcad(%)'].values.reshape(-1,)
y_true_vad = real_rst['挥发分Vad(%)'].values.reshape(-1,)
y_true_fcad = real_rst['固定炭Fcad(%)'].values.reshape(-1,)

In [56]:

pm25_eva = print_eva(y_true_vad, y_pred_vad, tp='挥发分Vad')
pm10_eva = print_eva(y_true_fcad, y_pred_fcad, tp='固定炭Fcad')

COL: 挥发分Vad, MSE: 3.35E-01,RMSE: 0.579,MAPE: 1.639 %,MAE: 0.504,R_2: 0.87
COL: 固定炭Fcad, MSE: 1.11E+00,RMSE: 1.055,MAPE: 1.497 %,MAE: 0.814,R_2: 0.876

In [ ]:

In [ ]: