• 基于变分自动编码器VAE的电池剩余使用寿命RUL估计

    加载模块

    1. import math
    2. import itertools
    3. import numpy as np
    4. import pandas as pd
    5. import seaborn as sns
    6. import tensorflow as tf
    7. from keras import layers
    8. from sklearn.svm import SVR
    9. from tensorflow import keras
    10. from keras import backend as K
    11. import matplotlib.pyplot as plt
    12. from keras.regularizers import l2
    13. from keras.regularizers import l1_l2
    14. from keras.models import load_model
    15. from keras.callbacks import ReduceLROnPlateau
    16. from sklearn.model_selection import KFold, GridSearchCV
    17. from sklearn.preprocessing import StandardScaler
    18. from sklearn.pipeline import make_pipeline
    19. from sklearn.model_selection import GroupShuffleSplit
    20. from sklearn.metrics import mean_squared_error, r2_score
    21. from sklearn.model_selection import train_test_split, cross_val_score, cross_validate
    22. from keras.callbacks import Callback, EarlyStopping, ModelCheckpoint, TensorBoard, LambdaCallback
    23. from keras.layers import Input, Dense, Lambda, LSTM, RepeatVector, Bidirectional, Masking, Dropout, BatchNormalization

    读取数据

    1. # Data reading
    2. data_dir = './original_data/'
    3.  
    4. # Specify the file path and column names
    5. file_path = os.path.join(data_dir, 'battery_RUL.txt')
    6.  
    7. # Specify the column names
    8. column_names = [
    9.     'unit_nr',
    10.     'time_cycles',
    11.     's_discharge_t',
    12.     's_decrement_3.6-3.4V',
    13.     's_max_voltage_discharge',
    14.     's_min_voltage_charge',
    15.     'Time_at_4.15V_s',
    16.     's_time_constant_current',
    17.     's_charging_time',
    18.     'RUL'
    19. ]
    20.  
    21. # Read the text file into a DataFrame
    22. df = pd.read_csv(file_path, sep='\s+', header=None, skiprows=0, names=column_names)
    23.  
    24. # Calculate the unique number of batteries
    25. unique_batteries = df['unit_nr'].nunique()
    26. print("\n Unique number of batteries:", unique_batteries)
    27.  
    28. # Check for missing values
    29. if df.isnull().any().any():
    30.     print("\n There are missing values in the DataFrame.")
    31.  
    32. # Print the shape of the DataFrame
    33. print("\n DataFrame shape:", df.shape)
    34.  
    35. # print df info
    36. print(df.info())
    37.  
    38. # Print the first few rows of the DataFrame
    39. print("\n First few rows of the DataFrame:")
    40. df.head()

    Unique number of batteries: 14

    DataFrame shape: (15064, 10)

    RangeIndex: 15064 entries, 0 to 15063
    Data columns (total 10 columns):
    # Column Non-Null Count Dtype
    --- ------ -------------- -----
    0 unit_nr 15064 non-null float64
    1 time_cycles 15064 non-null float64
    2 s_discharge_t 15064 non-null float64
    3 s_decrement_3.6-3.4V 15064 non-null float64
    4 s_max_voltage_discharge 15064 non-null float64
    5 s_min_voltage_charge 15064 non-null float64
    6 Time_at_4.15V_s 15064 non-null float64
    7 s_time_constant_current 15064 non-null float64
    8 s_charging_time 15064 non-null float64
    9 RUL 15064 non-null int64
    dtypes: float64(9), int64(1)
    memory usage: 1.1 MB
    None

    First few rows of the DataFrame:

    图片

    数据前处理

    1. # Filter rows where 'time_cycles' is equal to 1
    2. df_at_cycle_1 = df[df['time_cycles'] == 1]
    3.  
    4. # Create a bar chart with all battery units on the X-axis
    5. plt.figure(figsize=(12, 6))
    6. plt.bar(df_at_cycle_1['unit_nr'], df_at_cycle_1['RUL'])
    7. plt.xlabel('Battery Number')
    8. plt.ylabel('RUL at Cycle 1')
    9. plt.title('RUL Values at the Beginning of Cycle 1 for Each Battery Unit')
    10. plt.xticks(df_at_cycle_1['unit_nr'])  # Set X-axis ticks explicitly
    11. plt.tight_layout()
    12. plt.show()

    图片

    1. plt.figure(figsize = (8,8))
    2. sns.heatmap(df.corr(),annot=True, cbar=False, cmap='Blues', fmt='.1f')

    图片

    Correlation between RUL and:

    Max. Voltage Dischar. (V) is 0.8

    Min. Voltage Charg. (V) is -0.8

    Time at 4.15V (s) is 0.2

    Cycle index is -1.0

    Discharge Time (s), Decrement 3.6-3.4V (s), Time constant current (s) and Charging time (s) are 0.

    And correlation between Time at 4.15V and these four features are 0.8, 0.5,0.6 and 0.7.

    1. df1=df.drop(['s_discharge_t', 's_decrement_3.6-3.4V', 's_time_constant_current','s_charging_time'], axis=1)
    2. plt.figure(figsize = (4,4))
    3. sns.heatmap(df1.corr(),annot=True, cbar=False, cmap='Blues', fmt='.1f')

    图片

    df1.head()

    图片

    1. def exponential_smoothing(df, sensors, n_samples, alpha=0.4):
    2.     df = df.copy()
    3.     # first, take the exponential weighted mean
    4.     df[sensors] = df.groupby('unit_nr', group_keys=True)[sensors].apply(lambda x: x.ewm(alpha=alpha).mean()).reset_index(level=0, drop=True)
    5.  
    6.     # second, drop first n_samples of each unit_nr to reduce filter delay
    7.     def create_mask(data, samples):
    8.         result = np.ones_like(data)
    9.         result[0:samples] = 0
    10.         return result
    11.  
    12.     mask = df.groupby('unit_nr')['unit_nr'].transform(create_mask, samples=n_samples).astype(bool)
    13.     df = df[mask]
    14.  
    15.     return df
    16. #-----------------------------------------------------------------------------------------------------------------------
    17. def data_standardization(df, sensors):
    18.     df = df.copy()
    19.     # Apply StandardScaler to the sensor data
    20.     scaler = StandardScaler()
    21.     df[sensors] = scaler.fit_transform(df[sensors])
    22.  
    23.     return df
    24.  
    25.     # MMS_X = MinMaxScaler()
    26.     # mms_y = MinMaxScaler()
    27. #-----------------------------------------------------------------------------------------------------------------------
    28. def gen_train_data(df, sequence_length, columns):
    29.     data = df[columns].values
    30.     num_elements = data.shape[0]
    31.  
    32.     # -1 and +1 because of Python indexing
    33.     for start, stop in zip(range(0, num_elements-(sequence_length-1)), range(sequence_length, num_elements+1)):
    34.         yield data[start:stop, :]
    35. #-----------------------------------------------------------------------------------------------------------------------
    36. def gen_data_wrapper(df, sequence_length, columns, unit_nrs=np.array([])):
    37.     if unit_nrs.size <= 0:
    38.         unit_nrs = df['unit_nr'].unique()
    39.  
    40.     data_gen = (list(gen_train_data(df[df['unit_nr']==unit_nr], sequence_length, columns))
    41.                for unit_nr in unit_nrs)
    42.     data_array = np.concatenate(list(data_gen)).astype(np.float32)
    43.     return data_array
    44. #-----------------------------------------------------------------------------------------------------------------------
    45. def gen_labels(df, sequence_length, label):
    46.     data_matrix = df[label].values
    47.     num_elements = data_matrix.shape[0]
    48.  
    49.     # -1 because I want to predict the rul of that last row in the sequence, not the next row
    50.     return data_matrix[sequence_length-1:num_elements, :]
    51. #-----------------------------------------------------------------------------------------------------------------------
    52. def gen_label_wrapper(df, sequence_length, label, unit_nrs=np.array([])):
    53.     if unit_nrs.size <= 0:
    54.         unit_nrs = df['unit_nr'].unique()
    55.  
    56.     label_gen = [gen_labels(df[df['unit_nr']==unit_nr], sequence_length, label)
    57.                 for unit_nr in unit_nrs]
    58.     label_array = np.concatenate(label_gen).astype(np.float32)
    59.     return label_array
    60. #---------Original code------------------------------------------------------------------------------------------------------
    61. def gen_test_data(df, sequence_length, columns, mask_value):
    62.     if df.shape[0] < sequence_length:
    63.         data_matrix = np.full(shape=(sequence_length, len(columns)), fill_value=mask_value) # pad
    64.         idx = data_matrix.shape[0] - df.shape[0]
    65.         data_matrix[idx:,:] = df[columns].values  # fill with available data
    66.     else:
    67.         data_matrix = df[columns].values
    68.  
    69.     # specifically yield the last possible sequence
    70.     stop = data_matrix.shape[0]
    71.     start = stop - sequence_length
    72.     for i in list(range(1)):
    73.         yield data_matrix[start:stop, :]
    74. #---------------------------------------------------------------------------------------------------------------------------
    75. def change_test_index(df, initial_unit_number, start_index, end_index, min_rows, max_rows):
    76.     df.reset_index(drop=True, inplace=True)
    77.     y_test = []  # Initialize an empty list to store y_test values
    78.     y_test.append(df.loc[0, 'RUL'])
    79.  
    80.     while end_index < len(df):
    81.         # Calculate the number of rows to be assigned to the current unit_nr
    82.         num_rows = min_rows + (end_index % (max_rows - min_rows + 1))
    83.         end_index = end_index + num_rows
    84.  
    85.         # Update the unit_nr for the current block of rows
    86.         df.loc[start_index:end_index, 'unit_nr'] = initial_unit_number
    87.  
    88.         # Update the "time_cycles" column starting from the first row with "time_cycles" == 1
    89.         time_cycle_to_change = end_index + 1
    90.         df.loc[time_cycle_to_change, 'time_cycles'] = 1
    91.  
    92.         # Append the values to y_test (replace 'column_name' with the actual column name you want to append)
    93.         y_test.append(df.loc[time_cycle_to_change, 'RUL'])
    94.  
    95.         # Update the starting and ending index for the next block of rows
    96.         start_index = end_index + 1
    97.         initial_unit_number += 1
    98.  
    99.     # Drop rows with NaN values at the end of DataFrame 'df'
    100.     df.dropna(axis=0, how='any', inplace=True)
    101.  
    102.     # Drop any NaN values at the end of list 'y_test'
    103.     while len(y_test) > 0 and pd.isnull(y_test[-1]):
    104.         y_test.pop()
    105.  
    106.     return df, y_test

    获取数据

    1. def get_data(df, sensors, sequence_length, alpha):
    2.     # List of battery units
    3.     battery_units = range(1, 15)
    4.     sensor_names = ['s_{}'.format(i+1) for i in range(0,10)]
    5.     # Define the number of batteries for training and testing
    6.     train_batteries = 12
    7.     test_batteries = 2
    8.  
    9.     # Extract the batteries for training and testing
    10.     train_units = battery_units[:train_batteries]
    11.     test_units = battery_units[train_batteries:train_batteries + test_batteries]
    12.  
    13.     # Create the training and testing datasets
    14.     train = df[df['unit_nr'].isin(train_units)].copy()  # Use copy to avoid the SettingWithCopyWarning
    15.     test = df[df['unit_nr'].isin(test_units)].copy()  # Use copy to avoid the SettingWithCopyWarning
    16.  
    17.     X_test_pre, y_test = change_test_index(test, 13, 0, 0, 230, 240)
    18.  
    19.     # y_test = pd.Series(y_test)  # Convert y_test list to a pandas Series
    20.  
    21.     # remove unused sensors
    22.     drop_sensors = [element for element in sensor_names if element not in sensors]
    23.  
    24.     # Apply standardization to the training and testing data using data_standardization function
    25.     standard_train = data_standardization(train, sensors)
    26.     standard_test = data_standardization(test, sensors)
    27.  
    28.     # Exponential smoothing of training and testing data
    29.     X_train_pre = exponential_smoothing(standard_train, sensors, 0, alpha)
    30.     X_test_pre = exponential_smoothing(standard_test, sensors, 0, alpha)
    31.  
    32.     # Train-validation split
    33.     gss = GroupShuffleSplit(n_splits=1, train_size=0.85, random_state=42)
    34.     # Generate the train/val for each sample by iterating over the train and val units
    35.     for train_unit, val_unit in gss.split(X_train_pre['unit_nr'].unique(), groups=X_train_pre['unit_nr'].unique()):
    36.         train_unit = X_train_pre['unit_nr'].unique()[train_unit]  # gss returns indexes and index starts at 1
    37.         val_unit = X_train_pre['unit_nr'].unique()[val_unit]
    38.  
    39.         x_train = gen_data_wrapper(X_train_pre, sequence_length, sensors, train_unit)
    40.         y_train = gen_label_wrapper(X_train_pre, sequence_length, ['RUL'], train_unit)
    41.  
    42.         x_val = gen_data_wrapper(X_train_pre, sequence_length, sensors, val_unit)
    43.         y_val = gen_label_wrapper(X_train_pre, sequence_length, ['RUL'], val_unit)
    44.  
    45.     # create sequences for test
    46.     test_gen = (list(gen_test_data(X_test_pre[X_test_pre['unit_nr']==unit_nr], sequence_length, sensors, -99.))
    47. 			  for unit_nr in X_test_pre['unit_nr'].unique())
    48.  
    49.     x_test = np.concatenate(list(test_gen)).astype(np.float32)
    50.  
    51.     return x_train, y_train, x_val, y_val, x_test, y_test

    训练回调

    1. # --------------------------------------- TRAINING CALLBACKS  ---------------------------------------
    2. class save_latent_space_viz(Callback):
    3. 	def __init__(self, model, data, target):
    4. 		self.model = model
    5. 		self.data = data
    6. 		self.target = target
    7.  
    8. 	def on_train_begin(self, logs={}):
    9. 		self.best_val_loss = 100000
    10.  
    11. 	def on_epoch_end(self, epoch, logs=None):
    12. 		encoder = self.model.layers[0]
    13. 		if logs.get('val_loss') < self.best_val_loss:
    14. 			self.best_val_loss = logs.get('val_loss')
    15. 			viz_latent_space(encoder, self.data, self.target, epoch, True, False)
    16.  
    17. def get_callbacks(model, data, target):
    18. 	model_callbacks = [
    19. 		EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=30),
    20.         ModelCheckpoint(filepath='./checkpoints/checkpoint',monitor='val_loss', mode='min', verbose=1, save_best_only=True, save_weights_only=True),
    21. 		TensorBoard(log_dir='./logs'),
    22.         save_latent_space_viz(model, data, target)
    23. 	]
    24. 	return model_callbacks
    25.  
    26. def viz_latent_space(encoder, data, targets=[], epoch='Final', save=False, show=True):
    27.     z, _, _  = encoder.predict(data)
    28.     plt.figure(figsize=(3, 3))  # Smaller figsize value to reduce the plot size
    29.     if len(targets) > 0:
    30.         plt.scatter(z[:, 1], z[:, 0],  c=targets)
    31.     else:
    32.         plt.scatter(z[:, 1], z[:, 0])
    33.     plt.xlabel('z - dim 1')
    34.     plt.ylabel('z - dim 2')
    35.     plt.colorbar()
    36.     if show:
    37.         plt.show()
    38.     if save:
    39.         plt.savefig('./images/latent_space_epoch' + str(epoch) + '.png')
    40.     return z

    最佳学习率

    1. # ----------------------------------------- FIND OPTIMAL LR  ----------------------------------------
    2. class LRFinder:
    3.     """
    4.     Cyclical LR, code tailored from:
    5.     https://towardsdatascience.com/estimating-optimal-learning-rate-for-a-deep-neural-network-ce32f2556ce0
    6.     """
    7.  
    8.     def __init__(self, model):
    9.         self.model = model
    10.         self.losses = []
    11.         self.lrs = []
    12.         self.best_loss = 1e9
    13.  
    14.     def on_batch_end(self, batch, logs):
    15.         # Log the learning rate
    16.         lr = K.get_value(self.model.optimizer.lr)
    17.         self.lrs.append(lr)
    18.  
    19.         # Log the loss
    20.         loss = logs['loss']
    21.         self.losses.append(loss)
    22.  
    23.         # Check whether the loss got too large or NaN
    24.         if batch > 5 and (math.isnan(loss) or loss > self.best_loss * 4):
    25.             self.model.stop_training = True
    26.             return
    27.  
    28.         if loss < self.best_loss:
    29.             self.best_loss = loss
    30.  
    31.         # Increase the learning rate for the next batch
    32.         lr *= self.lr_mult
    33.         K.set_value(self.model.optimizer.lr, lr)
    34.  
    35.     def find(self, x_train, y_train, start_lr, end_lr, batch_size=64, epochs=1, **kw_fit):
    36.         # If x_train contains data for multiple inputs, use length of the first input.
    37.         # Assumption: the first element in the list is single input; NOT a list of inputs.
    38.         N = x_train[0].shape[0] if isinstance(x_train, list) else x_train.shape[0]
    39.  
    40.         # Compute number of batches and LR multiplier
    41.         num_batches = epochs * N / batch_size
    42.         self.lr_mult = (float(end_lr) / float(start_lr)) ** (float(1) / float(num_batches))
    43.         # Save weights into a file
    44.         initial_weights = self.model.get_weights()
    45.  
    46.         # Remember the original learning rate
    47.         original_lr = K.get_value(self.model.optimizer.lr)
    48.  
    49.         # Set the initial learning rate
    50.         K.set_value(self.model.optimizer.lr, start_lr)
    51.  
    52.         callback = LambdaCallback(on_batch_end=lambda batch, logs: self.on_batch_end(batch, logs))
    53.  
    54.         self.model.fit(x_train, y_train,
    55.                        batch_size=batch_size, epochs=epochs,
    56.                        callbacks=[callback],
    57.                        **kw_fit)
    58.  
    59.         # Restore the weights to the state before model fitting
    60.         self.model.set_weights(initial_weights)
    61.  
    62.         # Restore the original learning rate
    63.         K.set_value(self.model.optimizer.lr, original_lr)
    64.  
    65.     def find_generator(self, generator, start_lr, end_lr, epochs=1, steps_per_epoch=None, **kw_fit):
    66.         if steps_per_epoch is None:
    67.             try:
    68.                 steps_per_epoch = len(generator)
    69.             except (ValueError, NotImplementedError) as e:
    70.                 raise e('`steps_per_epoch=None` is only valid for a'
    71.                         ' generator based on the '
    72.                         '`keras.utils.Sequence`'
    73.                         ' class. Please specify `steps_per_epoch` '
    74.                         'or use the `keras.utils.Sequence` class.')
    75.         self.lr_mult = (float(end_lr) / float(start_lr)) ** (float(1) / float(epochs * steps_per_epoch))
    76.  
    77.         # Save weights into a file
    78.         initial_weights = self.model.get_weights()
    79.  
    80.         # Remember the original learning rate
    81.         original_lr = K.get_value(self.model.optimizer.lr)
    82.  
    83.         # Set the initial learning rate
    84.         K.set_value(self.model.optimizer.lr, start_lr)
    85.  
    86.         callback = LambdaCallback(on_batch_end=lambda batch,
    87.                                                       logs: self.on_batch_end(batch, logs))
    88.  
    89.         self.model.fit_generator(generator=generator,
    90.                                  epochs=epochs,
    91.                                  steps_per_epoch=steps_per_epoch,
    92.                                  callbacks=[callback],
    93.                                  **kw_fit)
    94.  
    95.         # Restore the weights to the state before model fitting
    96.         self.model.set_weights(initial_weights)
    97.  
    98.         # Restore the original learning rate
    99.         K.set_value(self.model.optimizer.lr, original_lr)
    100.  
    101.     def plot_loss(self, n_skip_beginning=10, n_skip_end=5, x_scale='log'):
    102.         """
    103.         Plots the loss.
    104.         Parameters:
    105.             n_skip_beginning - number of batches to skip on the left.
    106.             n_skip_end - number of batches to skip on the right.
    107.         """
    108.         plt.ylabel("loss")
    109.         plt.xlabel("learning rate (log scale)")
    110.         plt.plot(self.lrs[n_skip_beginning:-n_skip_end], self.losses[n_skip_beginning:-n_skip_end])
    111.         plt.xscale(x_scale)
    112.         plt.show()
    113.  
    114.     def plot_loss_change(self, sma=1, n_skip_beginning=10, n_skip_end=5, y_lim=(-0.01, 0.01)):
    115.         """
    116.         Plots rate of change of the loss function.
    117.         Parameters:
    118.             sma - number of batches for simple moving average to smooth out the curve.
    119.             n_skip_beginning - number of batches to skip on the left.
    120.             n_skip_end - number of batches to skip on the right.
    121.             y_lim - limits for the y axis.
    122.         """
    123.         derivatives = self.get_derivatives(sma)[n_skip_beginning:-n_skip_end]
    124.         lrs = self.lrs[n_skip_beginning:-n_skip_end]
    125.         plt.ylabel("rate of loss change")
    126.         plt.xlabel("learning rate (log scale)")
    127.         plt.plot(lrs, derivatives)
    128.         plt.xscale('log')
    129.         plt.ylim(y_lim)
    130.         plt.show()
    131.  
    132.     def get_derivatives(self, sma):
    133.         assert sma >= 1
    134.         derivatives = [0] * sma
    135.         for i in range(sma, len(self.lrs)):
    136.             derivatives.append((self.losses[i] - self.losses[i - sma]) / sma)
    137.         return derivatives
    138.  
    139.     def get_best_lr(self, sma, n_skip_beginning=10, n_skip_end=5):
    140.         derivatives = self.get_derivatives(sma)
    141.         best_der_idx = np.argmin(derivatives[n_skip_beginning:-n_skip_end])
    142.         return self.lrs[n_skip_beginning:-n_skip_end][best_der_idx]

    结果

    1. # --------------------------------------------- RESULTS  --------------------------------------------
    2. def get_model(path):
    3.     saved_VRAE_model = load_model(path, compile=False)
    4.  
    5.     # return encoder, regressor
    6.     return saved_VRAE_model.layers[1], saved_VRAE_model.layers[2]
    7.  
    8. def evaluate(y_true, y_hat, label='test'):
    9.     mse = mean_squared_error(y_true, y_hat)
    10.     rmse = np.sqrt(mse)
    11.     variance = r2_score(y_true, y_hat)
    12.     print('{} set RMSE:{}, R2:{}'.format(label, rmse, variance))
    13.     return rmse, variance
    14.  
    15. def score(y_true, y_hat):
    16.   res = 0
    17.   for true, hat in zip(y_true, y_hat):
    18.     subs = hat - true
    19.     if subs < 0:
    20.       res = res + np.exp(-subs/10)[0]-1
    21.     else:
    22.       res = res + np.exp(subs/13)[0]-1
    23.   print("score: ", res)
    24.  
    25. def results(path, x_train, y_train, x_test, y_test):
    26.     # Get model
    27.     encoder, regressor = get_model(path)
    28.     # Latent space
    29.     train_mu = viz_latent_space(encoder, x_train, y_train)
    30.     test_mu = viz_latent_space(encoder, x_test, y_test)
    31.     # Evaluate
    32.     y_hat_train = regressor.predict(train_mu)
    33.     y_hat_test = regressor.predict(test_mu)
    34.  
    35.     evaluate(y_train, y_hat_train, 'train')
    36.     evaluate(y_test, y_hat_test, 'test')
    37.     score(y_test, y_hat_test)

    模型结构

    1. class Sampling(keras.layers.Layer):
    2.   """Uses (z_mean, sigma) to sample z, the vector encoding an engine trajetory."""
    3.   def call(self, inputs):
    4.     mu, sigma = inputs
    5.     batch = tf.shape(mu)[0]
    6.     dim = tf.shape(mu)[1]
    7.     epsilon = tf.keras.backend.random_normal(shape=(batch, dim))
    8.     return mu + tf.exp(0.5 * sigma) * epsilon

    1. class RVE(keras.Model):
    2.     def __init__(self, encoder, regressor, decoder=None, **kwargs):
    3.         super(RVE, self).__init__(**kwargs)
    4.         self.encoder = encoder
    5.         self.regressor = regressor
    6.         self.total_loss_tracker = keras.metrics.Mean(name="total_loss")
    7.         self.kl_loss_tracker = keras.metrics.Mean(name="kl_loss")
    8.         self.reg_loss_tracker = keras.metrics.Mean(name="reg_loss")
    9.         self.decoder = decoder
    10.         if self.decoder!=None:
    11.           self.reconstruction_loss_tracker = keras.metrics.Mean(name="reconstruction_loss")
    12.  
    13.  
    14.     @property
    15.     def metrics(self):
    16.       if self.decoder!=None:
    17.         return [
    18.             self.total_loss_tracker,
    19.             self.kl_loss_tracker,
    20.             self.reg_loss_tracker,
    21.             self.reconstruction_loss_tracker
    22.         ]
    23.       else:
    24.         return [
    25.             self.total_loss_tracker,
    26.             self.kl_loss_tracker,
    27.             self.reg_loss_tracker,
    28.         ]
    29.  
    30.     def train_step(self, data):
    31.         x, target_x = data
    32.         with tf.GradientTape() as tape:
    33.             # kl loss
    34.             mu, sigma, z = self.encoder(x)
    35.             kl_loss = -0.5 * (1 + sigma - tf.square(mu) - tf.exp(sigma))
    36.             kl_loss = tf.reduce_mean(tf.reduce_sum(kl_loss, axis=1))
    37.             # Regressor
    38.             reg_prediction = self.regressor(z)
    39.             reg_loss = tf.reduce_mean(
    40.                 keras.losses.mse(target_x, reg_prediction)
    41.             )
    42.             # Reconstruction
    43.             if self.decoder!=None:
    44.               reconstruction = self.decoder(z)
    45.               reconstruction_loss = tf.reduce_mean(
    46.                   keras.losses.mse(x, reconstruction)
    47.               )
    48.               total_loss = kl_loss + reg_loss + reconstruction_loss
    49.               self.reconstruction_loss_tracker.update_state(reconstruction_loss)
    50.             else:
    51.               total_loss = kl_loss + reg_loss
    52.         grads = tape.gradient(total_loss, self.trainable_weights)
    53.         self.optimizer.apply_gradients(zip(grads, self.trainable_weights))
    54.         self.total_loss_tracker.update_state(total_loss)
    55.         self.kl_loss_tracker.update_state(kl_loss)
    56.         self.reg_loss_tracker.update_state(reg_loss)
    57.         return {
    58.             "loss": self.total_loss_tracker.result(),
    59.             "kl_loss": self.kl_loss_tracker.result(),
    60.             "reg_loss": self.reg_loss_tracker.result(),
    61.         }
    62.  
    63.  
    64.     def test_step(self, data):
    65.         x, target_x = data
    66.  
    67.         # kl loss
    68.         mu, sigma, z = self.encoder(x)
    69.         kl_loss = -0.5 * (1 + sigma - tf.square(mu) - tf.exp(sigma))
    70.         kl_loss = tf.reduce_mean(tf.reduce_sum(kl_loss, axis=1))
    71.         # Regressor
    72.         reg_prediction = self.regressor(z)
    73.         reg_loss = tf.reduce_mean(
    74.             keras.losses.mse(target_x, reg_prediction)
    75.         )
    76.         # Reconstruction
    77.         if self.decoder!=None:
    78.           reconstruction = self.decoder(z)
    79.           reconstruction_loss = tf.reduce_mean(
    80.               keras.losses.mse(x, reconstruction)
    81.           )
    82.  
    83.           total_loss = kl_loss + reg_loss + reconstruction_loss
    84.         else:
    85.           total_loss = kl_loss + reg_loss
    86.  
    87.         return {
    88.             "loss": total_loss,
    89.             "kl_loss": kl_loss,
    90.             "reg_loss": reg_loss,
    91.         }

    1. # Set hyperparameters
    2. sequence_length = 200
    3. alpha = 0.2
    4.  
    5. # Load and preprocess data
    6. df = df  # Load dataset
    7. sensors = ['s_max_voltage_discharge', 's_min_voltage_charge', "Time_at_4.15V_s"] # Define the sensors
    8.  
    9. # Call get_data_with_kfold to get the necessary data
    10. x_train, y_train, x_val, y_val, x_test, y_test = get_data(df1, sensors, sequence_length, alpha)
    11.  
    12. # from scipy.signal import savgol_filter
    13. # # Apply Savitzky-Golay filter
    14. # x_val_smoothed = savgol_filter(x_val, window_length=4, polyorder=2, axis=0)
    15.  
    16. # Setup the network parameters:
    17. timesteps = x_train.shape[1]
    18. input_dim = x_train.shape[2]
    19. intermediate_dim = 32
    20. batch_size = 256
    21. latent_dim = 2
    22. masking_value = -99 # used to mask values in sequences with less than 250 cycles until 250 is reached
    23. kernel_regularizer=l1_l2(l1=0.001, l2=0.001)
    24. dropout_rate = 0.1
    25. # --------------------------------- Encoder --------------------------------------
    26. inputs = Input(shape=(timesteps, input_dim,), name='encoder_input')
    27. mask = Masking(mask_value=masking_value)(inputs)
    28. h = Bidirectional(LSTM(intermediate_dim))(mask) # LSTM encoding
    29.  
    30. mu = Dense(latent_dim, kernel_regularizer = kernel_regularizer)(h) # VAE Z layer
    31. mu = Dropout(dropout_rate)(mu)
    32.  
    33. sigma = Dense(latent_dim, kernel_regularizer = kernel_regularizer)(h)
    34. sigma = Dropout(dropout_rate)(sigma)
    35.  
    36. z = Sampling()([mu, sigma])
    37.  
    38. # Instantiate the encoder model:
    39. encoder = keras.Model(inputs, [mu, sigma, z], name='encoder')
    40. # ------------------------------- Regressor --------------------------------------
    41. reg_latent_inputs = Input(shape=(latent_dim,), name='z_sampling_reg')
    42.  
    43. reg_intermediate = Dense(16, activation='tanh', kernel_regularizer = kernel_regularizer)(reg_latent_inputs)
    44. reg_intermediate = BatchNormalization()(reg_intermediate)
    45. reg_intermediate = Dropout(dropout_rate)(reg_intermediate)
    46.  
    47. reg_outputs = Dense(1, name='reg_output', kernel_regularizer = kernel_regularizer)(reg_intermediate)
    48. reg_outputs = Dropout(dropout_rate)(reg_outputs)
    49.  
    50. # Instantiate the classifier model:
    51. regressor = keras.Model(reg_latent_inputs, reg_outputs, name='regressor')
    52.  
    53. print("Shape of x_train:", x_train.shape)
    54. print("Shape of y_train:", y_train.shape)
    55. print("Shape of x_val:", x_val.shape)
    56. print("Shape of y_val:", y_val.shape)
    57. print("Shape of x_test:", x_test.shape)
    58. print("Shape of y_test:", len(y_test))

    Shape of x_train: (8795, 200, 3)
    Shape of y_train: (8795, 1)
    Shape of x_val: (1758, 200, 3)
    Shape of y_val: (1758, 1)
    Shape of x_test: (10, 200, 3)
    Shape of y_test: 10

    1. rve = RVE(encoder, regressor)
    2. lr_finder = LRFinder(rve)
    3.  
    4. rve.compile(optimizer=keras.optimizers.Adam(learning_rate=0.0000001))
    5.  
    6. # with learning rate growing exponentially from 0.0000001 to 1
    7. lr_finder.find(x_train, y_train, start_lr=0.0000001, end_lr=1, batch_size=batch_size, epochs=10)
    8.  
    9. # Plot the loss
    10. lr_finder.plot_loss()

    图片

    图片

    开始训练

    1. # Instantiate the RVE model
    2. rve = RVE(encoder, regressor)
    3.  
    4. #Compile the RVE model with the Adam optimizer
    5. rve.compile(optimizer=keras.optimizers.Adam(0.01))
    6.  
    7. # Define the early stopping callback
    8. early_stopping = EarlyStopping(monitor='loss', min_delta=3, patience=5,  verbose=1, mode='min', restore_best_weights=True)
    9.  
    10. # Call get_data_with_kfold to get the necessary data
    11. x_train, y_train, x_val, y_val, x_test, y_test = get_data(df1, sensors, sequence_length, alpha)
    12.  
    13. # Fit the RVE model with the callbacks
    14. rve.fit(x_train, y_train, epochs=500, batch_size=batch_size, validation_data=(x_val, y_val), callbacks=[early_stopping])

    RUL估计

    1. train_mu = viz_latent_space(rve.encoder, np.concatenate((x_train, x_val)), np.concatenate((y_train, y_val)))
    2. test_mu = viz_latent_space(rve.encoder, x_test, y_test)
    3.  
    4. # Evaluate
    5. y_hat_train = rve.regressor.predict(train_mu)
    6. y_hat_test = rve.regressor.predict(test_mu)
    7.  
    8. evaluate(np.concatenate((y_train, y_val)), y_hat_train, 'train')
    9. evaluate(y_test, y_hat_test, 'test')

    图片

    图片

    330/330 [==============================] - 0s 1ms/step 1/1 [==============================] - 0s 20ms/step train set RMSE:26.346721649169922, R2:0.9899616368349312 test set RMSE:246.51723899515346, R2:0.5192274671464132


    (246.51723899515346, 0.5192274671464132)

    工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。

  • 相关阅读:
    常用SQL语法总结
    嵌入式Linux入门—Linux多线程编程、互斥量、信号量、条件变量
    iOS15.4来袭：新增“男妈妈”表情及口罩面容解锁、AirTags反跟踪等新功能
    五、《图解HTTP》报文首部和HTTP缓存
    【Buildroot】工具包使用
    Android面试题 - 01
    js实现整数和小数分开并添加不同的样式
    「网页开发｜前端开发｜Vue」09 Vue状态管理Vuex：让页面根据用户登录状态渲染不同内容
    数据结构之八大排序——简单选择排序
    关于DDD和COLA的一些总结和思考
  • 原文地址：https://blog.csdn.net/weixin_39402231/article/details/139598235