背景
這裏將展示如何使用 Python 建立一個簡單但實用的數碼孿生。
鋰離子電池將成為我們的實物資產。
該數碼孿生將使我們能夠分析和預測電池行為，並且可以整合到任何虛擬資產管理工作流程中。
名詞1：虛擬系統和數碼孿生
數碼孿生（Digital Twin）是工業 4.0 的關鍵組成部份。
其基本原理是在虛擬世界中復制物理資產以對其動態進行建模。
數碼孿生由三部份組成：
- 物理實體 ：指實際存在的物理器材或系統。
- 數碼對映 ：指物理實體在數碼世界中的虛擬模型。
- 數據連線 ：指物理實體與數碼對映之間的即時數據傳輸和反饋。
數碼孿生透過傳感器和物聯網技術，將物理實體的即時數據傳輸到數碼對映中，並透過大數據分析、人工智能和仿真技術，對器材的執行狀態進行監測、分析和預測。
這個「twin」預計會以其「physical twin」相同的方式響應輸入變量。
該虛擬物件必須整合一個模型才能做到這一點。
數碼孿生最重要的特征是可以在數碼環境中模擬「物理」行為的模型。
當我們說「物理」時，我們指的是任何現實世界的實體（可以是鋰離子電池、水泵、人、城市或貓）。
任何可以建模的東西都可以虛擬化。
透過模型從物理世界到虛擬世界
名詞2：鋰離子電池數碼模型
可充電鋰離子電池是一種尖端電池技術，其電化學依賴於鋰離子。
除了便攜式技術器材之外，這些電池也是電動汽車和配電網絡儲能等套用的重要資產。
這些電池最關鍵的方面是它們的老化成本。
經過反復的充放電迴圈後，電池的電芯會退化，導致充電容量下降。
這種現象一直是開發更持久電池的關鍵研究領域。它的模型同樣也是一個有爭議的話題。
更多細節在下面程式碼中再詳細介紹。
完整程式碼
1、匯入包
import pandas as pdimport numpy as npfrom sklearn.model_selection import train_test_splitimport kerasimport tensorflow as tffrom keras.models importSequentialfrom keras.layers import LSTM,Densefrom keras.optimizers import SGDfrom tensorflow.keras.layers importDropoutfrom matplotlib import pyplot as pltimport plotly.express as pxfrom plotly.subplots import make_subplotsimport plotly.graph_objects as goimport plotly.figure_factory as ff
2、載入數據
df = pd.read_csv('discharge.csv')print( df.head())# 篩選出了電池編號為 'B0005' 的數據df = df[df['Battery']=='B0005']# 進一步篩選了溫度測量值大於 36 的數據df = df[df['Temperature_measured']>36]# 篩選後的 df 數據框按照 id_cycle 列進行分組，並對每個分組內的列取最大值dfb = df.groupby(['id_cycle']).max()# 添加了一個新列 Cumulated_T 到 dfb 數據框中，該列是 Time 列的累積和（累計時間）。# 這意味著 Cumulated_T 的每一行值是從第一個周期到當前周期結束的時間總和dfb['Cumulated_T']= dfb['Time'].cumsum()
數據集介紹：
- Voltage_measured ：測量的電壓值（單位：伏特）。
- Current_measured ：測量的電流值（單位：安培）。負值表示放電狀態。
- Temperature_measured ：測量的溫度值（單位：攝氏度）。
- Current_charge ：充電電流值（單位：安培）。
- Voltage_charge ：充電電壓值（單位：伏特）。
- Time ：測量時間點（單位：秒）。
- Capacity ：電池的容量（單位：安時）。
- id_cycle ：放電周期的編號。
- type ：操作類別（放電或充電）。
- ambient_temperature ：環境溫度（單位：攝氏度）。
- time ：數據記錄的年份。
輸出：
建立散點矩陣圖
展示 DataFrame 中各個變量之間的關系，透過散點矩陣圖可以一次性地看到多個變量之間的散點關系和分布情況，有助於數據探索和分析
import plotly.express as pxfig = px.scatter_matrix(dfb.drop(columns=['Time','type','ambient_temperature','time','Battery']),)fig.update_traces(marker=dict(size=2, color='crimson', symbol='square')),fig.update_traces(diagonal_visible=False)fig.update_layout( title='電池數據集', width=900, height=1200,)fig.update_layout({'plot_bgcolor':'#f2f8fd','paper_bgcolor':'white',},template='plotly_white', font=dict(size=7))fig.show()
建立散點矩陣圖
展示 DataFrame 中各個變量之間的關系，透過散點矩陣圖可以一次性地看到多個變量之間的散點關系和分布情況，有助於數據探索和分析
fig = go.Figure()fig.add_trace(go.Scatter(x=dfb['Cumulated_T']/3600, y=dfb['Capacity'], mode='lines', name='Capacity', marker_size=3, line=dict(color='crimson', width=3)))fig.update_layout( title="電池放電容量", xaxis_title="工作時間 [小時]", yaxis_title=f"電池容量（安時）")fig.update_layout({'plot_bgcolor':'#f2f8fd','paper_bgcolor':'white',},template='plotly_white')
3、定義物理模型
根據參考文獻 [1]，物理模型的基本方程式如下：
L =1−(1− L ′) e − fd
其中， L 表示電池的壽命， L ′ 表示初始電池壽命。 fd 是單位時間和迴圈的線性降解率。它可以描述為：
fd = fd ( t , δ , σ , Tc )
其中， t 表示充電時間， δ 表示迴圈放電深度， σ 表示迴圈平均充電狀態， Tc 表示電池溫度。電池容量的方程式可以寫為：
C = C 0 efd
經驗上發現， fd 可以近似表示為：
fd = tkTci
其中， k =0.13 ， i 表示迴圈數， t 表示每個迴圈的充電時間。
# 用來計算電池的壽命（L）和電池容量的變化from math import e# 計算電池壽命L = (dfb['Capacity']-dfb['Capacity'].iloc[0:1].values[0])/-dfb['Capacity'].iloc[0:1].values[0]K = 0.13# 來計算電池的壽命 L_1 = 1-e**(-K*dfb.index*dfb['Temperature_measured']/(dfb['Time']))dfb['C. Capacity'] = -(L_1*dfb['Capacity'].iloc[0:1].values[0]) + dfb['Capacity'].iloc[0:1].values[0]
比較了一個物理模型（Physical model）預測的電池容量變化和 NASA 數據集中記錄的實際電池容量
fig = go.Figure()fig.add_trace(go.Scatter(x=dfb.index, y=dfb['C. Capacity'], mode='lines', name='物理模型', line=dict(color='navy', width=2.5,)))fig.add_trace(go.Scatter(x=dfb.index, y=dfb['Capacity'], mode='markers', marker=dict( size=4, color='grey', symbol='cross'), name='NASA 數據集', line_color='navy'))fig.update_layout( title="物理模型比較", xaxis_title="迴圈數", yaxis_title=", 容量 [安時]")fig.update_layout(legend=dict( yanchor="top", y=0.9, xanchor="left", x=0.8))fig.update_layout({'plot_bgcolor':'#f2f8fd','paper_bgcolor':'white',},template='plotly_white')
4、將實驗數據與物理模型進行比較
評估物理模型預測的電池容量 (dfb['C. Capacity']) 與實際觀測數據 (dfb['Capacity']) 之間的平均絕對誤差。MAE 越小表示物理模型預測的精度越高。
# 平均絕對誤差M = pd.DataFrame()S = pd.DataFrame()def MAE(M,S): return np.sum(S-M)/len(S)print(f'平均絕對誤差 =', round(MAE(dfb['Capacity'], dfb['C. Capacity']), 3))
輸出：
平均絕對誤差 = 0.004
5、混合數碼孿生模型
# 輸入（X_in）和輸出（X_out）的定義，並將它們分割為訓練集和測試集。# 準備數據，以便將其用於訓練和測試混合數碼孿生模型或其他預測模型，然後評估模型的預測能力X_in = dfb['C. Capacity'] # input: the simulation time seriesX_out = dfb['Capacity'] - dfb['C. Capacity'] # output: difference between measurement and simulationX_in_train, X_in_test, X_out_train, X_out_test = train_test_split(X_in, X_out, test_size=0.33)
# 在Keras中，Dense函數構建一個全連線的神經網絡層，自動初始化權重和偏差。# 第一個隱藏層model = Sequential()model.add(Dense(64, activation='relu'))model.add(Dense(64, activation='relu'))model.add(Dense(1))
開始訓練
# 使用 Keras 訓練神經網絡模型，具體涉及模型的編譯和訓練過程。epochs =100loss ="mse"model.compile(optimizer='adam', loss=loss, metrics=['mae'],#Mean Absolute Error)history = model.fit(X_in_train, X_out_train, shuffle=True, epochs=epochs, batch_size=20, validation_data=(X_in_test, X_out_test), verbose=1)
建立了一個圖表，展示了神經網絡模型訓練過程中訓練集和驗證集的平均絕對誤差（MAE）隨著 epochs 的變化情況
fig = go.Figure()fig.add_trace(go.Scatter(x=np.arange(0, epochs,1), y=history.history['mae'], mode='lines', name=f'Training MAE', marker_size=3, line_color='orange'))fig.add_trace(go.Scatter(x=np.arange(0, epochs,1), y=history.history['val_mae'], mode='lines', name=f'Validation MAE', line_color='grey'))fig.update_layout( title="網絡訓練", xaxis_title="訓練輪數", yaxis_title=f"平均絕對誤差")fig.update_layout({'plot_bgcolor':'#f2f8fd','paper_bgcolor':'white',},template='plotly_white')
6、編譯混合數碼孿生模型
fig = go.Figure()fig.add_trace(go.Scatter(x=X_in_train, y=X_out_train, mode='markers', name=f'Modelled Capacity', marker=dict( size=4, color='grey', symbol='cross'), line_color='crimson'))fig.add_trace(go.Scatter(x = X_in_train, y=model.predict(X_in_train).reshape(-1), mode='lines', name=f'Trained Capacity', line=dict(color='navy', width=3)))fig.update_layout( title="網絡訓練", xaxis_title="模擬容量", yaxis_title="Δ (模擬容量 - 測量容量.)")fig.update_layout(legend=dict( yanchor="top", y=0.95, xanchor="left", x=0.45))fig.update_layout({'plot_bgcolor':'#f2f8fd',#or azure'paper_bgcolor':'white',},template='plotly_white')
X_twin = X_in + model.predict(X_in).reshape(-1)fig = go.Figure()fig.add_trace(go.Scatter(x=dfb.index, y=X_twin, mode='lines', name=f'Hybrid digial twin', line=dict(color='firebrick', width=3)))fig.add_trace(go.Scatter(x=dfb.index, y=dfb['C. Capacity'], mode='lines', name=f'Modelled capacity', line=dict(color='navy', width=3, dash='dash')))fig.add_trace(go.Scatter(x=dfb.index, y=dfb['Capacity'], mode='markers', marker=dict( size=4, color='grey', symbol='cross'), name=f'Observed capacity', line_color='navy'))fig.update_layout( title="混合數碼孿生模型與其他模型的比較", xaxis_title="迴圈數", yaxis_title="容量（安時）")fig.update_layout(legend=dict( yanchor="top", y=0.95, xanchor="left", x=0.77))fig.update_layout({'plot_bgcolor':'#f2f8fd','paper_bgcolor':'white',},template='plotly_white')
7、使用混合數碼孿生模型進行預測
import numpy as npfrom math import e# 定義迴圈數和溫度cycles = np.arange(168,500,1)temperature = dfb['Temperature_measured'].iloc[167]time = dfb['Time'].iloc[167]# 計算 L_eK =0.13L_e =1- np.exp(-K * cycles * temperature / time)# 計算模擬的初始容量X_in_e =-(L_e * dfb['Capacity'].iloc[0:1].values[0])+ dfb['Capacity'].iloc[0:1].values[0]# 轉換 X_in_e 的形狀和類別以符合模型的輸入要求X_in_e = X_in_e.reshape(-1,1).astype(np.float32)# 進行預測predicted_values = model.predict(X_in_e).reshape(-1)# 計算混合數碼孿生模型的預測容量C_twin_e = X_in_e + predicted_values# 打印結果，偵錯用print(f"Cycles: {cycles}")print(f"Temperature: {temperature}")print(f"Time: {time}")print(f"L_e shape: {L_e.shape}, X_in_e shape: {X_in_e.shape}, predicted_values shape: {predicted_values.shape}")print(f"C_twin_e: {C_twin_e}")
cycles = np.arange(168,500,1)temperature = dfb['Temperature_measured'].iloc[167]time = dfb['Time'].iloc[167]# 計算 L_eK =0.13L_e =1- np.exp(-K * cycles * temperature / time)# 計算模擬的初始容量X_in_e =-(L_e * dfb['Capacity'].iloc[0:1].values[0])+ dfb['Capacity'].iloc[0:1].values[0]# 轉換 X_in_e 的形狀和類別以符合模型的輸入要求X_in_e = X_in_e.reshape(-1,1).astype(np.float32)# 進行預測predicted_values = model.predict(X_in_e).reshape(-1)# 計算混合數碼孿生模型的預測容量C_twin_e = X_in_e + predicted_valuesprint(C_twin_e)
X_twin = X_in + model.predict(X_in).reshape(-1)fig = go.Figure()fig.add_trace(go.Scatter(x=cycles, y=X_in_e, mode='lines', name=f'C modelled (predicted)', line=dict(color='navy', width=3, dash='dash')))fig.add_trace(go.Scatter(x=cycles, y=C_twin_e, mode='lines', name=f'C Digital twin (predicted)', line=dict(color='crimson', width=3, dash='dash')))fig.add_trace(go.Scatter(x=dfb.index, y=X_twin, mode='lines', name=f'C Digital twin', line=dict(color='crimson', width=2)))fig.add_trace(go.Scatter(x=dfb.index, y=dfb['C. Capacity'], mode='lines', name=f'C modelled', line=dict(color='navy', width=2)))fig.update_layout( title="電池容量預測", xaxis_title="迴圈數", yaxis_title="電池容量 [安時]")fig.update_layout(legend=dict( yanchor="top", y=0.95, xanchor="left", x=0.72))fig.update_layout({'plot_bgcolor':'#f2f8fd','paper_bgcolor':'white',},template='plotly_white')
參考：
https://towardsdatascience.com/how-to-build-a-digital-twin-b31058fd5d3e
