import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import pandas as pd
import glob
import time
from qore_sdk.client import WebQoreClient
from tqdm import tqdm_notebook as tqdm
from ipywidgets import IntProgress

df_35 = pd.read_csv('CS2_35.csv',index_col=0).dropna()
df_36 = pd.read_csv('CS2_36.csv',index_col=0).dropna()
df_37 = pd.read_csv('CS2_37.csv',index_col=0).dropna()
df_38 = pd.read_csv('CS2_38.csv',index_col=0).dropna()

fig = plt.figure(figsize=(9,8),dpi=150)
names = ['capacity','resistance','CCCT','CVCT']
titles = ['Discharge Capacity (mAh)', 'Internal Resistance (Ohm)',
          'Constant Current Charging Time (s)','Constant Voltage Charging Time (s)']
plt.subplots_adjust(hspace=0.25)
for i in range(4):
    plt.subplot(2,2,i+1)
    plt.plot(df_35[names[i]],'o',ms=2,label='#35')
    plt.plot(df_36[names[i]],'o',ms=2,label='#36')
    plt.plot(df_37[names[i]],'o',ms=2,label='#37')
    plt.plot(df_38[names[i]],'o',ms=2,label='#38')
    plt.title(titles[i],fontsize=14)
    plt.legend(loc='upper right')
    if(i == 3):
        plt.ylim(1000,5000)
    plt.xlim(-20,1100)

X_train,y_train = ConvertData([df_35,df_37,df_38],50)
X_test,y_test = ConvertData([df_36],50)

print(X_train.shape,X_test.shape,y_train.shape,y_test.shape)

(2588, 50) (873, 50) (2588,) (873,)

idx = np.arange(0,X_train.shape[0],1)
idx = np.random.permutation(idx)
idx_lim = idx[:500]


X_train = X_train[idx_lim]
y_train = y_train[idx_lim]

X_train = X_train.reshape([X_train.shape[0], X_train.shape[1], 1])
X_test = X_test.reshape([X_test.shape[0], X_test.shape[1], 1])
print(X_train.shape,X_test.shape,y_train.shape,y_test.shape)

(500, 50, 1) (873, 50, 1) (500,) (873,)

%%time
client = WebQoreClient(username, password, endpoint=endpoint)
time_ = client.regression_train(X_train, y_train)
print('Time:', time_['train_time'], 's')

Time: 1.8583974838256836 s
Wall time: 2.7 s

res = client.regression_predict(X_train)


fig=plt.figure(figsize=(12, 4),dpi=150)
plt.plot(res['Y'],alpha=0.7,label='Prediction')
plt.plot(y_train,alpha=0.7,label='True')
plt.legend(loc='upper right',fnotallow=12)

res = client.regression_predict(X_test[300:800])

x_range = np.linspace(301,800,500)
fig=plt.figure(figsize=(12, 4),dpi=150)
plt.plot(x_range,y_test[300:800],label='true',c='gray')
plt.plot(x_range,res['Y'],label='Qore',c='red',alpha=0.8)
plt.xlabel('Number of Cycle',fnotallow=13)
plt.ylabel('DIscharge Capacity (Ah)',fnotallow=13)
plt.title('Prediction of Discharge Capacity of Test Data (CS2-36)',fnotallow=14)
plt.legend(loc='upper right',fnotallow=12)

from sklearn.metrics import mean_squared_error
mean_squared_error(y_test[300:800], res['Y'])

initial = X_test[500]
results = []
for i in tqdm(range(50)):
    if(i == 0):
        res = client.regression_predict([initial])['Y']
        results.append(res[0])
        time.sleep(1)
    else:
        initial = np.vstack((initial[1:],np.array(res)))
        res = client.regression_predict([initial])['Y']
        results.append(res[0])
        time.sleep(1)

HBox(children=(IntProgress(value=0, max=50), HTML(value='')))

fig=plt.figure(figsize=(12,4),dpi=150)
plt.plot(np.linspace(501,550,50),results,'o-',ms=4,lw=1,label='predict')
plt.plot(np.linspace(401,550,150),y_test[400:550],'o-',lw=1,ms=4,label='true')
plt.legend(loc='upper right',fnotallow=12)
plt.xlabel('Number of Cycle',fnotallow=13)
plt.ylabel('Discharge Capacity (Ah)',fnotallow=13)

https://www.zhihu.com/consult/people/792359672131756032?isMe=1

final_pred = []
count = 0
for i in range(100):
    j = max_cycles[i] 
    temp = pred[count+j-1]
    count=count+j
    final_pred.append(int(temp))


print(final_pred)




fig = plt.figure(figsize=(18,10))
plt.plot(final_pred,color='red', label='prediction')
plt.plot(y_test,color='blue', label='y_test')
plt.legend(loc='upper left')
plt.grid()
plt.show()


print("mean_squared_error >> ", mean_squared_error(y_test,final_pred))
print("root_mean_squared_error >> ", math.sqrt(mean_squared_error(y_test,final_pred)))
print("mean_absolute_error >>",mean_absolute_error(y_test,final_pred))

fig = plt.figure(figsize=(18,10))
plt.plot(logistic_pred,color='red', label='prediction')
plt.plot(y_test,color='blue', label='y_test')
plt.legend(loc='upper left')
plt.grid()
plt.show()

fig = plt.figure(figsize=(18,10))
plt.plot(lstm_pred,color='red', label='prediction')
plt.plot(y_test,color='blue', label='y_test')
plt.legend(loc='upper left')
plt.grid()
plt.show()

完整代碼可通過(guò)知乎學(xué)術(shù)咨詢獲得：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import random
import warnings
warnings.filterwarnings('ignore')


from sklearn.metrics import mean_squared_error, r2_score,mean_absolute_error
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler,MinMaxScaler


from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor,ExtraTreesRegressor
from sklearn.neighbors import KNeighborsRegressor
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import random
import time
import warnings
warnings.filterwarnings('ignore')


from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler,MinMaxScaler
#from google.colab import drive
#drive.mount('/content/drive')


# model
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, LSTM, Conv1D
from tensorflow.keras.layers import BatchNormalization, Dropout
from tensorflow.keras.layers import TimeDistributed, Flatten
from tensorflow.keras.layers.experimental import preprocessing
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.callbacks import ReduceLROnPlateau,EarlyStopping

https://www.zhihu.com/consult/people/792359672131756032?isMe=1

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt 
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
import tensorflow as tf
from tensorflow.keras.layers import Dense

tensorflow=2.8.0
keras=2.8.0
sklearn=1.0.2

path = ''


# define column names 
col_names=["unit_nb","time_cycle"]+["set_1","set_2","set_3"] + [f's_{i}' for i in range(1,22)]


# read data

df_train = train_data = pd.read_csv(path+"train_FD001.txt", index_col=False, 
                                    sep= "\s+", header = None,nam
es=col_names )

df_test和y_test同理導(dǎo)入，看一下訓(xùn)練樣本

df_train.head()

df_train[col_names[1:]].describe().T

sns.set_style("darkgrid")
plt.figure(figsize=(16,10))
k = 1
for col in col_names[2:] :
    plt.subplot(6,4,k)
    sns.histplot(df_train[col],color='Green')
    k+=1
plt.tight_layout()
plt.show()

def remaining_useful_life(df):
    # Get the total number of cycles for each unit
    grouped_by_unit = df.groupby(by="unit_nb")
    max_cycle = grouped_by_unit["time_cycle"].max()
    
    # Merge the max cycle back into the original frame
    result_frame = df.merge(max_cycle.to_frame(name='max_cycle'), left_notallow='unit_nb', right_index=True)
    
    # Calculate remaining useful life for each row
    remaining_useful_life = result_frame["max_cycle"] - result_frame["time_cycle"]
    result_frame["RUL"] = remaining_useful_life
    
    # drop max_cycle as it's no longer needed
    result_frame = result_frame.drop("max_cycle", axis=1)
    return result_frame
df_train = remaining_useful_life(df_train)
df_train.head()

plt.figure(figsize=(10,5))
sns.histplot(max_ruls.RUL, color='r')
plt.xlabel('RUL')
plt.ylabel('Frequency')
plt.axvline(x=max_ruls.RUL.mean(), ls='--',color='k',label=f'mean={max_ruls.RUL.mean()}')
plt.axvline(x=max_ruls.RUL.median(),color='b',label=f'median={max_ruls.RUL.median()}')
plt.legend()
plt.show()

plt.figure(figsize=(20, 8))

heatmap = sns.heatmap(cor_matrix, vmin=-1, vmax=1, annot=True)
heatmap.set_title('Correlation Heatmap', fnotallow={'fontsize':12}, pad=10);

col = df_train.describe().columns
#we drop colummns with standard deviation is less than 0.0001
sensors_to_drop = list(col[df_train.describe().loc['std']<0.001]) + ['s_14']
print(sensors_to_drop)
# 
df_train.drop(sensors_to_drop,axis=1,inplace=True)
df_test.drop(sensors_to_drop,axis=1,inplace=True)


sns.set_style("darkgrid")
fig, axs = plt.subplots(4,4, figsize=(25, 18), facecolor='w', edgecolor='k')
fig.subplots_adjust(hspace = .22, wspace=.2)
i=0
axs = axs.ravel()
index = list(df_train.unit_nb.unique())
for sensor in df_train.columns[1:-1]:
    for idx in index[1:-1:15]:
        axs[i].plot('RUL', sensor,data=df_train[df_train.unit_nb==idx])
        
    axs[i].set_xlim(350,0)
    axs[i].set(xticks=np.arange(0, 350, 25))
    axs[i].set_ylabel(sensor)
    axs[i].set_xlabel('Remaining Use Life')
    i=i+1

X_train = df_train[df_train.columns[3:-1]]
y_train = df_train.RUL
X_test = df_test.groupby('unit_nb').last().reset_index()[df_train.columns[3:-1]]


y_train = y_train.clip(upper=155)


# create evalute function for train and test data
def evaluate(y_true, y_hat):
    RMSE = np.sqrt(mean_squared_error(y_true, y_hat))
    R2_score = r2_score(y_true, y_hat)
    return [RMSE,R2_score];


#Make Dataframe which will contain results
Results = pd.DataFrame(columns=['RMSE-Train','R2-Train','RMSE-Test','R

2-Test','time-train (s)'])

import time
Sc = StandardScaler()
X_train1 = Sc.fit_transform(X_train)
X_test1 = Sc.transform(X_test)
# create and fit model
start = time.time()
lm = LinearRegression()
lm.fit(X_train1, y_train)
end_fit = time.time()- start
# predict and evaluate
y_pred_train = lm.predict(X_train1)
y_pred_test = lm.predict(X_test1)
Results.loc['LR']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results


def plot_prediction(y_test,y_pred_test,score):
    plt.style.use("ggplot")  
    fig, ax = plt.subplots(1, 2, figsize=(17, 4), gridspec_kw={'width_ratios': [1.2, 3]})
    fig.subplots_adjust(wspace=.12)
    ax[0].plot([min(y_test),max(y_test)],
         [min(y_test),max(y_test)],lw=3,c='r')
    ax[0].scatter(y_test,y_pred_test,lw=3,c='g')
    ax[0].annotate(text=('RMSE: ' + "{:.2f}".format(score[0]) +'\n' +
                  'R2: ' + "{:.2%}".format(score[1])), xy=(0,140), size='large');
    ax[0].set_title('Actual vs predicted RUL')
    ax[0].set_xlabel('Actual')
    ax[0].set_ylabel('Predicted');
    
    ax[1].plot(range(0,100),y_test,lw=2,c='r',label = 'actual')
    ax[1].plot(range(0,100),y_pred_test,lw=1,ls='--', c='b',label = 'prediction')
    ax[1].legend()


    ax[1].set_title('Actual vs predicted RUL')
    ax[1].set_xlabel('Engine num')
    ax[1].set_ylabel('RUL');


plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test))

# create and fit model
start = time.time()
svr = SVR(kernel="rbf", gamma=0.25, epsilnotallow=0.05)
svr.fit(X_train1, y_train)
end_fit = time.time()-start 
# predict and evaluate
y_pred_train = svr.predict(X_train1)
y_pred_test = svr.predict(X_test1)
Results.loc['SVM']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results


plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test))

start=time.time()
dtr = DecisionTreeRegressor(random_state=42, max_features='sqrt', max_depth=10, min_samples_split=10)
dtr.fit(X_train1, y_train)
end_fit =time.time()-start
# predict and evaluate
y_pred_train = dtr.predict(X_train1)
y_pred_test = dtr.predict(X_test1)
Results.loc['DTree']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results


plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test

))

from sklearn.neighbors import KNeighborsRegressor
# Evaluating on Train Data Set
start = time.time()
Kneigh = KNeighborsRegressor(n_neighbors=7)
Kneigh.fit(X_train1, y_train)
end_fit =time.time()-start
# predict and evaluate
y_pred_train = Kneigh.predict(X_train1)
y_pred_test = Kneigh.predict(X_test1)
Results.loc['KNeigh']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results


plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test))

start = time.time()
rf = RandomForestRegressor(n_jobs=-1, n_estimators=130,max_features='sqrt', min_samples_split= 2, max_depth=10, random_state=42)
rf.fit(X_train1, y_train)
y_hat_train1 = rf.predict(X_train1)
end_fit = time.time()-start
# predict and evaluate
y_pred_train = rf.predict(X_train1)
y_pred_test = rf.predict(X_test1)
Results.loc['RF']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results

plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test))

from sklearn.ensemble import GradientBoostingRegressor
# Evaluating on Train Data Set
start = time.time()
xgb_r = GradientBoostingRegressor(n_estimators=45, max_depth=10, min_samples_leaf=7, max_features='sqrt', random_state=42,learning_rate=0.11)
xgb_r.fit(X_train1, y_train)
end_fit =time.time()-start
# predict and evaluate
y_pred_train = xgb_r.predict(X_train1)
y_pred_test = xgb_r.predict(X_test1)
Results.loc['XGboost']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results


plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test))

from sklearn.ensemble import VotingRegressor
start=time.time()
Vot_R = VotingRegressor([("rf", rf), ("xgb", xgb_r)],weights=[1.5,1],n_jobs=-1)
Vot_R.fit(X_train1, y_train)
end_fit =time.time()-start
# predict and evaluate
y_pred_train = Vot_R.predict(X_train1)
y_pred_test = Vot_R.predict(X_test1)
Results.loc['VotingR']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results


plot_prediction(y_test.RUL,y_pred_test,evaluate(y_test, y_pred_test))

star=time.time()
model = tf.keras.models.Sequential()
model.add(Dense(32, activatinotallow='relu'))
model.add(Dense(64, activatinotallow='relu'))
model.add(Dense(128, activatinotallow='relu'))
model.add(Dense(128, activatinotallow='relu'))
model.add(Dense(1, activatinotallow='linear'))
model.compile(loss= 'msle', optimizer='adam', metrics=['msle'])
history = model.fit(x=X_train1,y=y_train, epochs = 40, batch_size = 64)
end_fit = time.time()-star


# predict and evaluate
y_pred_train = model.predict(X_train1)
y_pred_test = model.predict(X_test1)
Results.loc['ANN']=evaluate(y_train, y_pred_train)+evaluate(y_test, y_pred_test)+[end_fit]
Results