Merge branch 'dev' of http://git.fast-fun.cn:92/lmstack/data_analyze_platform into dev

LIB/MIDDLE/Anomaly Detection/V1_0_0/anomalyPCA.py

@@ -0,0 +1,119 @@
+#热失控预警：PCA异常指数
+
+import pandas as pd
+import numpy as np
+from scipy.signal import savgol_filter
+from sklearn.preprocessing import RobustScaler
+from sklearn.decomposition import PCA
+
+#筛选特征
+def makedataset1(df_data):
+    df_data=df_data.drop(['Unnamed: 0','总电流[A]','GSM信号','外电压','SOH[%]','开关状态','充电状态','故障等级','故障代码','绝缘电阻','上锁状态','加热状态','单体均衡状态','总输出状态'],axis=1,errors='ignore')
+    df_data=df_data.drop(["单体温度"+str(i) for i in range(1,5)],axis=1,errors='ignore')
+    df_data=df_data.drop(["其他温度"+str(i) for i in range(1,7)],axis=1,errors='ignore')
+    listV=[s for s in list(df_data) if '单体电压' in s]
+    for i in range(1,len(listV)+1):
+        df_data=df_data[(df_data['单体电压'+str(i)]>2200) & (df_data['单体电压'+str(i)]<4800)]
+    df_data=df_data[df_data['SOC[%]']>20]
+    df_data['时间']=[df_data.loc[i,'时间戳'][0:15] for i in df_data.index]
+    data_set=df_data.groupby('时间').mean()
+    for k in data_set.columns:
+        data_set[k]=savgol_filter(data_set[k],3,2)
+    return data_set
+
+#新建统计特征
+def makedataset2(df_data):
+    data_set=makedataset1(df_data)
+    listV=[s for s in list(df_data) if '单体电压' in s]
+    data_set["最低单体电压"]=data_set[["单体电压"+str(i) for i in range(1,len(listV)+1)]].min(axis=1)
+    data_set["最高单体电压"]=data_set[["单体电压"+str(i) for i in range(1,len(listV)+1)]].max(axis=1)
+    data_set["平均单体电压"]=data_set[["单体电压"+str(i) for i in range(1,len(listV)+1)]].mean(axis=1)
+    data_set["最大单体压差"]=[data_set.loc[k,"最高单体电压"]-data_set.loc[k,"最低单体电压"] for k in data_set.index]
+    data_set["低压差"]=[data_set.loc[k,"平均单体电压"]-data_set.loc[k,"最低单体电压"] for k in data_set.index]
+    data_set=data_set.drop(["单体电压"+str(i) for i in range(1,len(listV)+1)],axis=1)
+    return data_set
+
+#标准化
+def process(data_set):
+    features=data_set.columns
+    sX=RobustScaler(copy=True)
+    data_set2=data_set.copy()
+    data_set2.loc[:,features]=sX.fit_transform(data_set2[features])
+    return data_set2
+
+#异常指数函数
+def anomalyScores(originalDF,reducedDF):
+    loss=np.sum((np.array(originalDF)-np.array(reducedDF))**2,axis=1)
+    loss=pd.Series(data=loss,index=originalDF.index)
+    loss=(loss-np.min(loss))/(np.max(loss)-np.min(loss))
+    return loss
+
+#建立PCA模型
+def anomalyPCA(x_train_pro):
+    n_components=4
+    whiten=True
+    random_state=2
+    pca=PCA(n_components=n_components,whiten=whiten,random_state=random_state)
+    pca.fit(x_train_pro)
+    return pca
+
+#判断PCA异常指数
+def transform(df_data_pro,model,df_data):
+    #降维
+    X_train=model.transform(df_data_pro)
+    X_train=pd.DataFrame(data=X_train,index=df_data_pro.index)
+    #还原
+    X_train_inverse=model.inverse_transform(X_train)
+    X_train_inverse=pd.DataFrame(data=X_train_inverse,index=df_data_pro.index)
+    #异常指数
+    anomalyScoresModel=anomalyScores(df_data_pro,X_train_inverse)
+    anomalyScoresModel=savgol_filter(anomalyScoresModel,15,3)
+    df_data2=df_data.copy()
+    df_data2['anomalyScores_'+str(model)]=anomalyScoresModel
+    return df_data2
+
+#判断离群
+def detect_outliers(data,pred,threshold=3):
+    anomaly=data['anomalyScores_PCA(n_components=4, random_state=2, whiten=True)']
+    anomalypred=pred['anomalyScores_PCA(n_components=4, random_state=2, whiten=True)']
+    mean_d=np.mean(anomaly.values)
+    std_d=np.std(anomaly.values)
+    max_score=np.max(anomaly.values)
+    outliers2=pd.DataFrame()
+    for k in anomalypred.index:
+        z_score= (anomalypred[k]-mean_d)/std_d
+        if (np.abs(z_score) >threshold) & (anomalypred[k]>max_score):
+            outliers2=outliers2.append(pred[anomalypred.values==anomalypred[k]])
+    return outliers2
+
+#训练模型
+def train_model(data_train):
+    x_train1=makedataset1(data_train) 
+    x_train2=makedataset2(data_train)  
+    x_train_pro1=process(x_train1) 
+    x_train_pro2=process(x_train2) 
+    pca1=anomalyPCA(x_train_pro1) 
+    pca2=anomalyPCA(x_train_pro2) 
+    res1=transform(x_train_pro1,pca1,x_train1)
+    res2=transform(x_train_pro2,pca2,x_train2)
+    return pca1,pca2,res1,res2
+
+#预测
+def prediction(data_test,pca1,pca2):
+    x_test1=makedataset1(data_test) 
+    x_test2=makedataset2(data_test) 
+    x_test_pro1=process(x_test1) 
+    x_test_pro2=process(x_test2) 
+    pred1=transform(x_test_pro1,pca1,x_test1)
+    pred2=transform(x_test_pro2,pca2,x_test2)
+    return pred1,pred2
+
+#判定异常
+def check_anomaly(outliers1,outliers2):
+    if (len(outliers1)>0) & (len(outliers2)>0):
+        outliers=pd.merge(outliers1,outliers2,on='时间')
+        outliers=outliers[outliers['SOC[%]_x']>45]
+        outliers=outliers.drop(['总电压[V]_y','单体压差_y','SOC[%]_y'],axis=1)
+        return outliers
+    
+

LIB/MIDDLE/Anomaly Detection/V1_0_0/main_anomalyPCA.py

@@ -0,0 +1,42 @@
+#热失控预警：PCA异常指数
+#训练模型
+
+from LIB.BACKEND import DBManager
+dbManager = DBManager.DBManager()
+from LIB.MIDDLE.CellStateEstimation.Common import log
+import pandas as pd
+from anomalyPCA import *
+import joblib
+import datetime
+
+dataSOH = pd.read_excel('sn-20210903.xlsx',sheet_name='sn-20210903')
+fileNames = dataSOH['sn']
+fileNames = list(fileNames)
+l = len(fileNames)
+
+now_time=datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')   #type: str
+now_time=datetime.datetime.strptime(now_time,'%Y-%m-%d %H:%M:%S')     #type: datetime
+start_time=now_time-datetime.timedelta(days=365)
+end_time=str(now_time)
+start_time=str(start_time)
+
+mylog=log.Mylog('log.txt','error')
+mylog.logcfg()
+
+for k in range(l): 
+    try: 
+        sn = fileNames[k]
+        df_data = dbManager.get_data(sn=sn, start_time=start_time, end_time=end_time, data_groups=['bms'])
+        data_train = df_data['bms']
+        
+        if len(data_train)>0:
+            pca1,pca2,res1,res2=train_model(data_train) 
+            joblib.dump(pca1,'pca1_'+sn+'.m')  
+            joblib.dump(pca2,'pca2_'+sn+'.m')  
+            res1.to_csv('res1_'+sn+'.csv',encoding='gbk')
+            res2.to_csv('res2_'+sn+'.csv',encoding='gbk')
+    
+    except Exception as e:
+        print(repr(e))
+        mylog.logopt(sn,e)
+        pass

LIB/MIDDLE/CellStateEstimation/BatSafetyWarning/V1_0_1/CBMSBatInterShort.py

@@ -120,10 +120,10 @@ class BatInterShort():
     def _celldeltsoc_get(self,cellvolt_list,dict_baltime,capacity): 
         cellsoc=[]
         celldeltsoc=[]
-        for j in range(self.param.CellVoltNums):    #获取每个电芯电压对应的SOC值
+        for j in range(self.param.CellVoltNums):   #获取每个电芯电压对应的SOC值
             cellvolt=cellvolt_list[j]
             ocv_soc=np.interp(cellvolt,self.param.LookTab_OCV,self.param.LookTab_SOC)
-            if j in dict_baltime.keys():
+            if j+1 in dict_baltime.keys():
                 ocv_soc=ocv_soc+dict_baltime[j+1]*self.param.BalCurrent/(capacity*3600)   #补偿均衡电流
             else:
                 pass
@@ -337,7 +337,7 @@ class BatInterShort():
 
         #容量初始化
         if self.df_soh.empty:
-            batsoh=self.df_bms.loc[0,'SOH[%]']
+            batsoh=(self.df_bms.loc[0,'SOH[%]'])
             capacity=self.param.Capacity*batsoh/100
         else:
             batsoh=self.df_soh.loc[len(self.df_soh)-1,'soh']
@@ -540,7 +540,7 @@ class BatInterShort():
 
         #容量初始化
         if self.df_soh.empty:
-            batsoh=self.df_bms.loc[0,'SOH[%]']
+            batsoh=(self.df_bms.loc[0,'SOH[%]'])
             capacity=self.param.Capacity*batsoh/100
         else:
             batsoh=self.df_soh.loc[len(self.df_soh)-1,'soh']
@@ -639,7 +639,7 @@ class BatInterShort():
                         pass                
                 elif standingtime>3600*12:
                     standingtime=0
-                    cellvolt_now=np.array(self._avgvolt_get(i))
+                    cellvolt_now=self._avgvolt_get(i)
                     if not cellvolt_now.empty:
                         cellvolt_min=min(cellvolt_now)
                         cellvolt_max=max(cellvolt_now)

LIB/MIDDLE/CellStateEstimation/BatSafetyWarning/V1_0_1/CBMSSafetyWarning.py

@@ -3,6 +3,7 @@ import numpy as np
 import datetime
 import time
 from matplotlib import pyplot as plt
+import pymannkendall as mk
 from LIB.MIDDLE.CellStateEstimation.Common.V1_0_1 import BatParam
 
 class SafetyWarning:
@@ -31,6 +32,7 @@ class SafetyWarning:
         df_res=pd.DataFrame(columns=['start_time', 'end_time', 'product_id', 'code', 'level', 'info','advice'])
 
         time_now=datetime.datetime.now()
+        time_now=time_now.strftime('%Y-%m-%d %H:%M:%S')
         time_sp='0000-00-00 00:00:00'
 
         #参数初始化.......................................
@@ -41,6 +43,13 @@ class SafetyWarning:
         R2_list=[]
         voltsigmafault_list=[]
         uniformfault_list=[]
+        mk_trend_list=[]
+        mk_p_list=[]
+        mk_z_list=[]
+        mk_Tau_list=[]
+        mk_slope_list=[]
+        mk_s_list=[]
+        mk_svar_list=[]
 
         if not self.df_short.empty:
             short_current=self.df_short['short_current']
@@ -80,8 +89,7 @@ class SafetyWarning:
             cellvolt_rank=self.df_uniform['cellvolt_rank']
             cellvolt_rank=cellvolt_rank.str.replace("[", '')
             cellvolt_rank=cellvolt_rank.str.replace("]", '')
-
-        # plt.figure()
+    
         for i in range(self.param.CellVoltNums):
             #漏电流故障判断...........................................................................
             if not self.df_short.empty:
@@ -103,11 +111,11 @@ class SafetyWarning:
             #电压变化率及电压离群度.................................................................................
             if not self.OutLineVol_Rate.empty and VoltChange.iloc[-1]['time']*36000>18*3600 and len(VoltChange)>5:
 
+                VoltChange[volt_column[i]]=VoltChange[volt_column[i]].map(lambda x:eval(x))
+                y=VoltChange[volt_column[i]]
                 volt3sigma=np.array(Volt_3Sigma[volt_column[i]].map(lambda x:eval(x)))
                 volt3sigma_sum=np.sum(volt3sigma<-3)
                 #电压变化率
-                VoltChange[volt_column[i]]=VoltChange[volt_column[i]].map(lambda x:eval(x))
-                y=VoltChange[volt_column[i]]
                 a1,b1=np.polyfit(VoltChange['time'].tolist(),y.tolist(),1)
                 y1=a1*VoltChange['time']+b1
                 y_mean=y.mean()
@@ -115,22 +123,55 @@ class SafetyWarning:
                 R2_list.append(R2)
                 
                 volt_rate.append(a1)
-                # plt.plot(xtime1,y1,label=a1)
+                # plt.plot(xtime1,y1,label='单体电压'+str(i+1))
+                # plt.xlabel('时间', fontsize=25)
+                # plt.ylabel('SOC差', fontsize=25)
+                # plt.xticks(fontsize=20)
+                # plt.yticks(fontsize=20)
                 # plt.scatter( xtime1,VoltChange[volt_column[i]],marker='o')
-                # plt.legend(loc='best')
-                # plt.title('单体电压'+str(i+1))
+                # plt.legend(bbox_to_anchor=(1, 0), loc=3, fontsize=16)
+                # plt.title(self.celltype)
                 # plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
                 # plt.rcParams['axes.unicode_minus']=False #用来正常显示负号  
-                # plt.show()
+                # # plt.show()
                 
                 if volt3sigma_sum>len(volt3sigma)/2:
                     voltsigmafault=1
                 else:
                     voltsigmafault=0
+                voltsigmafault_list.append(voltsigmafault)
+                
+                #mana-kendell趋势检验
+                mk_res=mk.regional_test(np.array(y)) 
+                mk_trend_list.append(mk_res.trend)
+                mk_p_list.append(mk_res.p)
+                mk_z_list.append(mk_res.z)
+                mk_Tau_list.append(mk_res.Tau)
+                mk_slope_list.append(mk_res.slope)
+                mk_s_list.append(mk_res.s)
+                mk_svar_list.append(mk_res.var_s)
+                """
+                trend：趋势；
+                h：有无趋势；
+                p：趋势的显著水平，越小趋势越明显；
+                z：检验统计量，正代表随时间增大趋势，负代表随时间减小趋势；
+                Tau：反映两个序列的相关性，接近1的值表示强烈的正相关，接近-1的值表示强烈的负相关；
+                s：Mann-Kendal的分数,如果S是一个正数，那么后一部分的观测值相比之前的观测值会趋向于变大；如果S是一个负数，那么后一部分的观测值相比之前的观测值会趋向于变小
+                slope：趋势斜率
+                """
+                # print('单体电压{}：\n'.format(i+1), mk_res)
+
             else:
                 volt_rate.append(0)
                 R2_list.append(0)
-            voltsigmafault_list.append(voltsigmafault)
+                voltsigmafault_list.append(0)
+                mk_trend_list.append(0)
+                mk_p_list.append(0)
+                mk_z_list.append(0)
+                mk_Tau_list.append(0)
+                mk_slope_list.append(0)
+                mk_s_list.append(0)
+                mk_svar_list.append(0)
 
             #电芯SOC排名判断.............................................................................
             if not self.df_uniform.empty:
@@ -151,13 +192,29 @@ class SafetyWarning:
                 faultinfo='电芯{}发生热失控安全预警'.format(i+1)
                 faultadvice='1联系用户远离电池，立刻召回电池'
                 df_res.loc[0]=[time_now, time_sp, self.sn, faultcode, faultlv, faultinfo, faultadvice]
+                break
             else:
-                pass       
+                pass  
         
-        #电池电压变化率热失控预警...............................................................................
+        # plt.show()   
+        #电池电压变化率离群度计算...............................................................................
         volt_rate_std=np.std(volt_rate)
         volt_rate_mean=np.mean(volt_rate)
         volt_rate_3sigma=(np.array(volt_rate)-volt_rate_mean)/volt_rate_std
+        
+        #mk离群度计算
+        mk_slope_std=np.std(mk_slope_list)
+        mk_slope_mean=np.mean(mk_slope_list)
+        mk_slope_3sigma=(np.array(mk_slope_list)-mk_slope_mean)/mk_slope_std
+
+        # mk_s_std=np.std(mk_s_list)
+        # mk_s_mean=np.mean(mk_s_list)
+        # mk_s_3sigma=(np.array(mk_s_list)-mk_s_mean)/mk_s_std 
+
+        mk_z_std=np.std(mk_z_list)
+        mk_z_mean=np.mean(mk_z_list)
+        mk_z_3sigma=(np.array(mk_z_list)-mk_z_mean)/mk_z_std  
+
         if not self.df_soh.empty and self.celltype<50:
             cellsoh=eval(self.df_soh.loc[0,'cellsoh'])
             cellsoh_std=np.std(cellsoh)
@@ -165,8 +222,35 @@ class SafetyWarning:
             cellsoh_3sigma=((np.array(cellsoh)-cellsoh_mean)/cellsoh_std)
         else:
             cellsoh_3sigma=[0]*self.param.CellVoltNums
-        for i in range(len(volt_rate)):
-            if volt_rate[i]<self.param.TrwVoltRate and volt_rate_3sigma[i]<-3 and abs(cellsoh_3sigma[i])<2.5 and  voltsigmafault_list[i]==1:
+        
+        #电压/SOC变化率
+        # for i in range(len(volt_rate)):
+        #     if volt_rate[i]<self.param.TrwVoltRate and volt_rate_3sigma[i]<-3 and abs(cellsoh_3sigma[i])<2.5 and  voltsigmafault_list[i]==1:
+        #         faultcode=110
+        #         faultlv=4
+        #         faultinfo='电芯{}发生热失控安全预警'.format(i+1)
+        #         faultadvice='2联系用户远离电池，立刻召回电池'
+        #         df_res.loc[0]=[time_now, time_sp, self.sn, faultcode, faultlv, faultinfo, faultadvice]
+        #     else:
+        #         pass
+
+        #mana-kendall趋势检测
+        for i in range(len(mk_p_list)):
+            #适用动态工况判断
+            if mk_trend_list[i]=='decreasing' and mk_p_list[i]<self.param.mk_p and mk_z_list[i]<self.param.mk_z and mk_Tau_list[i]<self.param.mk_Tau and mk_slope_3sigma[i]<-3 and mk_slope_list[i]<self.param.mk_slope and abs(cellsoh_3sigma[i])<2.5 and volt_rate_3sigma[i]<-3:
+                faultcode=110
+                faultlv=4
+                faultinfo='电芯{}发生热失控安全预警'.format(i+1)
+                faultadvice='2联系用户远离电池，立刻召回电池'
+                df_res.loc[0]=[time_now, time_sp, self.sn, faultcode, faultlv, faultinfo, faultadvice]
+            #适用静态工况判断
+            elif self.celltype<=50 and mk_trend_list[i]=='decreasing' and mk_p_list[i]<self.param.mk_p and mk_z_list[i]<-6 and (mk_Tau_list[i]<-0.9 or mk_z_3sigma[i]<-3) and mk_slope_3sigma[i]<-3.5 and mk_slope_list[i]<-0.025 and volt_rate_3sigma[i]<-3:
+                faultcode=110
+                faultlv=4
+                faultinfo='电芯{}发生热失控安全预警'.format(i+1)
+                faultadvice='2联系用户远离电池，立刻召回电池'
+                df_res.loc[0]=[time_now, time_sp, self.sn, faultcode, faultlv, faultinfo, faultadvice]
+            elif self.celltype>50 and mk_trend_list[i]=='decreasing' and mk_p_list[i]<self.param.mk_p and mk_z_list[i]<-6 and (mk_Tau_list[i]<-0.9 or mk_z_3sigma[i]<-3) and mk_slope_3sigma[i]<-3.5 and mk_slope_list[i]<-0.25 and volt_rate_3sigma[i]<-3:
                 faultcode=110
                 faultlv=4
                 faultinfo='电芯{}发生热失控安全预警'.format(i+1)
@@ -174,7 +258,5 @@ class SafetyWarning:
                 df_res.loc[0]=[time_now, time_sp, self.sn, faultcode, faultlv, faultinfo, faultadvice]
             else:
                 pass
-        
-        # plt.show()
 
         return df_res

LIB/MIDDLE/CellStateEstimation/BatSafetyWarning/V1_0_1/VoltStray.py

@@ -79,8 +79,8 @@ def main(sn,df_bms,celltype):
     df_bms['time']=pd.to_datetime(df_bms['时间戳'], format='%Y-%m-%d %H:%M:%S')
     volt_column = ['单体电压'+str(i) for i in range(1,param.CellVoltNums+1)]
     columns=['time']+volt_column
-    df_bms=df_bms[df_bms['SOC[%]']>10]
-    # df_standing=df_bms[(df_bms['PackCrnt']>-0.1) & (df_bms['PackCrnt']<0.1)]
+    df_bms=df_bms[(df_bms['SOC[%]']>10)]
+    df_bms=df_bms[(df_bms['PackCrnt']<1)]
     # df_chrg=df_bms[(df_bms['PackCrnt']<-1)]
 
     #电压/SOC变化率计算
@@ -99,7 +99,6 @@ def main(sn,df_bms,celltype):
 
         VolChng = cal_volt_change(df_soc,volt_column)
     else:
-        df_bms=df_bms[df_bms['SOC[%]']<95]
         # df_bms=df_bms[(df_bms['PackCrnt']>-0.1) & (df_bms['PackCrnt']<0.1)]
         df_ori = df_bms[columns]
         df = df_ori.drop_duplicates(subset=['time']) # 删除时间相同的数据
@@ -117,7 +116,7 @@ def main(sn,df_bms,celltype):
 
     OutLineVol=DataFrame(columns=['time','sn','VolOl_Uni','VolChng_Uni'])
 
-    #静置电压变化率和离群度计算
+    #电压变化率和离群度计算
     if len(VolChng)>5 and len(VolSigma)>5:
         VolChng['time'] = time_list1
         VolChng= VolChng.set_index('time')

LIB/MIDDLE/CellStateEstimation/Common/V1_0_1/BatParam.py

@@ -29,6 +29,12 @@ class BatParam:
         self.SohLow=70
         self.SohDiff=15
 
+        #mana-kendall趋势检验参数
+        self.mk_p=0.001
+        self.mk_z=-5
+        self.mk_Tau=-0.7
+        self.mk_slope=-0.1
+
         self.OcvWeight_Temp=[-30,-20,-10,0,10,20,30,40,50]
         self.OcvWeight_StandingTime=[0,500,600,1200,1800,3600,7200,10800]
         self.OcvWeight            =[[0,0,  0,  0,    0,   0.1,0.3, 1],
@@ -91,7 +97,7 @@ class BatParam:
             self.CellTempDiffLv1=10
             self.CellTempDiffLv2=15
 
-            self.TrwVoltRate=-1
+            # self.TrwVoltRate=-1
             
             self.DifVolGap = 3
             self.CellOVlmt=5
@@ -136,7 +142,7 @@ class BatParam:
             self.CellTempDiffLv1=10
             self.CellTempDiffLv2=15
 
-            self.TrwVoltRate=-1
+            # self.TrwVoltRate=-1
             
             self.DifVolGap = 3
             self.CellOVlmt=5
@@ -181,7 +187,7 @@ class BatParam:
             self.CellTempDiffLv1=10
             self.CellTempDiffLv2=15  
 
-            self.TrwVoltRate=-1
+            # self.TrwVoltRate=-1
             
             self.DifVolGap = 3
             self.CellOVlmt=5
@@ -201,8 +207,8 @@ class BatParam:
             self.PackCrntDec=-1
             self.BalCurrent=0.015
 
-            self.LookTab_SOC = [-5,   0,	    5,	    10,	    15,	    20,	    25,	    30,	    35,	    40,	    45,	    50,	    55,	    60,	    65,	    70,	    75,	    80,	    85,	    90,	    95,	    100,   105,     110]
-            self.LookTab_OCV = [3.0,  3.152, 	3.397, 	3.438, 	3.481, 	3.523, 	3.560, 	3.586, 	3.604, 	3.620, 	3.638, 	3.661, 	3.693, 	3.748, 	3.803, 	3.853, 	3.903, 	3.953, 	4.006, 	4.063, 	4.121, 	4.183, 4.253,  4.50]
+            self.LookTab_SOC = [0,	    5,	    10,	    15,	    20,	    25,	    30,	    35,	    40,	    45,	    50,	    55,	    60,	    65,	    70,	    75,	    80,	    85,	    90,	    95,	    100,   105,     110]
+            self.LookTab_OCV = [3.152, 	3.397, 	3.438, 	3.481, 	3.523, 	3.560, 	3.586, 	3.604, 	3.620, 	3.638, 	3.661, 	3.693, 	3.748, 	3.803, 	3.853, 	3.903, 	3.953, 	4.006, 	4.063, 	4.121, 	4.183, 4.253,  4.50]
 
 
             self.CellOvLv1=4.3
@@ -229,7 +235,7 @@ class BatParam:
             self.CellTempDiffLv1=10
             self.CellTempDiffLv2=15   
 
-            self.TrwVoltRate=-1    
+            # self.TrwVoltRate=-1    
             
             self.DifVolGap = 3
             self.CellOVlmt=5
@@ -285,7 +291,8 @@ class BatParam:
             self.CellTempDiffLv1=10
             self.CellTempDiffLv2=15   
 
-            self.TrwVoltRate=-8
+            # self.TrwVoltRate=-8
+            self.mk_slope=-0.3
             
             
             self.DifVolGap = 3
@@ -339,7 +346,8 @@ class BatParam:
             self.CellTempDiffLv1=28
             self.CellTempDiffLv2=32   
 
-            self.TrwVoltRate=-8  
+            # self.TrwVoltRate=-8 
+            self.mk_slope=-0.3
             
                                    
             self.DifVolGap = 3

LIB/MIDDLE/odo/DailyMileageEstimation/V1_0_4/calcul_mileage.py

@@ -0,0 +1,321 @@
+import pandas as pd
+import numpy as np
+import os
+import math
+import datetime
+import time
+
+#将时间戳由 "%Y-%m-%d %H:%M:%S" 切换为 sec
+def timeconvert(df_in,column_name):  
+    df=df_in.copy()
+    df.index=range(len(df))
+    time=df[column_name]
+    timeInSeries=[]
+    time2=datetime.datetime.strptime(str(time[0]),"%Y-%m-%d %H:%M:%S")
+    for k in range(len(time)):
+        time1=datetime.datetime.strptime(str(time[k]),"%Y-%m-%d %H:%M:%S")    
+        t=(time1-time2)
+        timeInSeries.append(t.days*86400+t.seconds)
+    df.loc[:,'相对时间[s]']=pd.DataFrame(timeInSeries,columns=['相对时间[s]'])
+    return df
+
+#计算累积能量
+def cal_accumKwh(df_in):
+    df_in1=df_in.copy()
+    df1=df_in1[['总电流[A]','总电压[V]','相对时间[s]']]
+    df1=df1.dropna(axis=0,how='any')
+    I=df1['总电流[A]'].values
+    V=df1['总电压[V]'].values
+    t=df1['相对时间[s]'].values
+    accumAh=[0.0]
+    for k in range(1,len(I)):
+        accumAh_temp=(t[k]-t[k-1])*((I[k]+I[k-1])/(2*3600))*(V[k]+V[k-1])/2/1000
+        accumAh.append(accumAh[-1]+accumAh_temp)
+    df1.loc[:,'累积能量[Kwh]']=accumAh
+    df_out=pd.merge(df_in1,df1[['累积能量[Kwh]']],how='left',left_index=True, right_index=True)
+    return(df_out)
+
+#将时间格式化为整数
+def str_data_to_num(str_data):
+    # 格式时间成毫秒
+    strptime = time.strptime(str_data,"%Y-%m-%d %H:%M:%S")
+    # print("strptime",strptime)
+    mktime = int(time.mktime(strptime)*1000)
+    # print("mktime",mktime)
+    return mktime
+
+def df_date_To_int(df_in):
+    df=df_in.copy()
+    for k in range(len(df)):
+        df.loc[k,'绝对时间[ms]']=str_data_to_num(df['时间戳'][k])
+    return df
+
+#根据经纬度获取两点之间的距离
+def cal_dis_meters(radius,latitude1, longitude1,latitude2, longitude2):  
+    radLat1 = (math.pi/180)*latitude1  
+    radLat2 = (math.pi/180)*latitude2  
+    radLng1 = (math.pi/180)*longitude1  
+    radLng2= (math.pi/180)*longitude2      
+    d=2*math.asin(math.sqrt(math.pow(math.sin((radLat1-radLat2)/2.0),2)+math.cos(radLat1)*math.cos(radLat2)*math.pow(math.sin((radLng1-radLng2)/2.0),2)))*radius
+    return d
+
+#第一部分距离计算
+def cal_distance1(part1,avg_cost):
+    part1.index=[i for i in range(len(part1))]
+    part1=part1.copy()
+    part1.loc[0,'△距离']=0
+    AccumEnergy1=part1['累积能量[Kwh]'].values
+    for k1 in range(1,len(part1)):
+        if (AccumEnergy1[k1]-AccumEnergy1[k1-1])>0:
+            part1.loc[k1,'△距离']=(AccumEnergy1[k1]-AccumEnergy1[k1-1])*avg_cost
+            part1.loc[k1,'方法']=4
+        else:
+            part1.loc[k1,'△距离']=0
+            part1.loc[k1,'方法']=5
+    return part1
+
+#第二部分计算
+def cal_distance2(part2,avg_cost):
+    part2_gps=part2[part2['纬度']>0]
+    part2_gps.index=[i for i in range(len(part2_gps))]
+    times_list=part2_gps['相对时间[s]'].values
+    lat=part2_gps['纬度'].values
+    lng=part2_gps['经度'].values
+    AccumEnergy2=part2_gps['累积能量[Kwh]'].values
+    part2_gps=part2_gps.copy()
+    part2_gps.loc[0,'△距离']=0
+    for k2 in range(1,len(part2_gps)):
+        delta_energy=AccumEnergy2[k2]-AccumEnergy2[k2-1]
+        delta_span=cal_dis_meters(6378.137,lat[k2], lng[k2],lat[k2-1], lng[k2-1])
+        v_spd=3600*delta_span/(times_list[k2]-times_list[k2-1]) 
+        if times_list[k2]-times_list[k2-1]<60:
+            if v_spd>70 :
+                part2_gps.loc[k2,'△距离']=delta_energy*avg_cost
+                part2_gps.loc[k2,'方法']=1_1
+            else:
+                part2_gps.loc[k2,'△距离']=delta_span
+                part2_gps.loc[k2,'方法']=1_2
+        else:            
+            if delta_energy<=0:
+                part2_gps.loc[k2,'△距离']=0
+                part2_gps.loc[k2,'方法']=2
+            else:
+                if v_spd>70:
+                    part2_gps.loc[k2,'△距离']=delta_energy*avg_cost
+                    part2_gps.loc[k2,'方法']=3_1
+                else:
+                    part2_gps.loc[k2,'△距离']=max(delta_span,delta_energy*avg_cost)
+                    part2_gps.loc[k2,'方法']=3_2
+    return part2_gps
+
+#第三部分距离计算
+def cal_distance3(part3,avg_cost):
+    part3.index=[i for i in range(len(part3))]
+    part3=part3.copy()
+    part3.loc[0,'△距离']=0
+    AccumEnergy3=part3['累积能量[Kwh]'].values
+    for k3 in range(1,len(part3)):
+        if (AccumEnergy3[k3]-AccumEnergy3[k3-1])>0:
+            part3.loc[k3,'△距离']=(AccumEnergy3[k3]-AccumEnergy3[k3-1])*avg_cost
+            part3.loc[k3,'方法']=4
+        else:
+            part3.loc[k3,'△距离']=0
+            part3.loc[k3,'方法']=5
+    return part3
+
+def real_odo(df_in,avg_cost):
+    # df_handle=df_in.copy()
+    df=timeconvert(df_in,"时间戳")#计算相对时间
+    df=cal_accumKwh(df)#计算累积能量
+    positive_lat_index=df[df['纬度']>0].index
+    if len(positive_lat_index)>2:
+        first_index=positive_lat_index[0]
+        end_index=positive_lat_index[-1]
+        #将数据分割为有第一个GPS数据之前（part1）、
+        #有最后一个GPS数据之后(part3)以及
+        #有第一个GPS数据和最后一个GPS数据中间的数据(part3)
+        part1=df[0:first_index+1]
+        part2=df[first_index:end_index+1]
+        part3=df[end_index:-1]
+        part1=cal_distance1(part1,avg_cost)
+        part2_gps=cal_distance2(part2,avg_cost)
+        part3=cal_distance3(part3,avg_cost) 
+        df_out=pd.concat([part1,part2_gps[1:],part3[1:]])
+    else:
+        part1=df
+        part2=df[0:0]
+        part3=df[0:0]
+        part1=cal_distance1(part1,avg_cost)
+        df_out=part1.copy()
+    
+    #计算累积里程
+    df_out.loc[0,'累积里程[km]']=0
+    df_out.index=[i for i in range(len(df_out))]
+
+    for k in range(1,len(df_out)):
+        df_out.loc[k,'累积里程[km]']=df_out.loc[k-1,'累积里程[km]']+df_out.loc[k,'△距离']
+    return(df_out)
+
+def calcul_mileage(sn,data_bms,data_gps):
+
+    #合并两张表格
+    df_bms=data_bms.copy()
+    df_gps=data_gps.copy()
+    #删除纬度小于10的点
+    for k in range(1,len(df_gps)):
+        if df_gps.loc[k,'纬度']<10:
+            df_gps=df_gps.drop(k)
+    df_bms.set_index(["时间戳"], inplace=True)
+    df_gps.set_index(["时间戳"], inplace=True)
+    merge_df1=pd.merge(df_bms, df_gps,how='outer', left_index=True, right_index=True)
+    merge_df1.loc[:,'时间戳']=merge_df1.index
+    merge_df1.index=[i for i in range(len(merge_df1))]
+
+    #参数定义
+    cost_min=34 #最小能耗km/kwh
+    cost_max=45 #最高能耗km/kwh
+
+    df_input=merge_df1[['时间戳','总电流[A]', '总电压[V]','SOC[%]','充电状态','纬度', '经度']]
+    #电流电压数据修复
+    df_input=timeconvert(df_input,"时间戳")#计算相对时间
+    related_times=df_input['相对时间[s]'].values
+    lat_list=df_input['纬度'].values
+    lng_list=df_input['经度'].values
+    #Vlt_list=df_input['总电压[V]'].values
+    for k in range(1,len(df_input)):
+        if math.isnan(df_input.loc[k,'总电流[A]']):
+            if related_times[k]-related_times[k-1]<5:
+                df_input.loc[k,'总电流[A]']=df_input.loc[k-1,'总电流[A]']
+                df_input.loc[k,'总电压[V]']=df_input.loc[k-1,'总电压[V]']
+                df_input.loc[k,'修复']=1
+            elif lat_list[k-1]>0:
+                delta_odo=cal_dis_meters(6378.137,lat_list[k], lng_list[k],lat_list[k-1], lng_list[k-1])
+                df_input.loc[k,'总电流[A]']=1000*3600*(delta_odo/cost_min)/70/(related_times[k]-related_times[k-1])
+                df_input.loc[k,'总电压[V]']=70
+                df_input.loc[k,'修复']=2
+            elif related_times[k]-related_times[k-1]<30:
+                df_input.loc[k,'总电流[A]']=df_input.loc[k-1,'总电流[A]']
+                df_input.loc[k,'总电压[V]']=df_input.loc[k-1,'总电压[V]']
+                df_input.loc[k,'修复']=3
+            else:
+                df_input=df_input.drop(k)
+
+    df_input.index=[i for i in range(len(df_input))]
+    df_input['日期']=[df_input.loc[i,'时间戳'][0:10] for i in range(len(df_input))]     
+    date_list=np.unique(df_input['日期'].values)
+
+    #list_result=[]
+    input_res=pd.DataFrame()
+    res = pd.DataFrame()
+    df_input_copy=df_input.copy()
+    for date_str in date_list:
+        df_input_copy=df_input[df_input['日期']==date_str]
+        data_len=len(df_input_copy)
+        time_list=df_input_copy['时间戳'].values
+        #初始化参数
+        start_time=datetime.datetime.strptime(time_list[0], "%Y-%m-%d %H:%M:%S")
+        last_time=time_list[0]
+        #reset_time=start_time
+        start_index=0
+        end_index=1
+        cal_count=0
+        #df_deltaData_last=df_input_copy[0:0]
+        #初始化输出
+        df_result=pd.DataFrame()
+        
+        df_result.loc[0,'时间']=start_time
+        df_result.loc[0,'时间']=str(df_result.loc[0,'时间'])
+        df_result.loc[0,'时间'] =df_result.loc[0,'时间'].rpartition(':')[0]
+        
+        df_result.loc[0,'累积里程[km]']=0
+        df_result.loc[0,'能耗[km/kwh]']=cost_min
+        #process_store=[]
+        avg_cost=34
+        cross_odo=0
+        for k in range(1,data_len):
+            
+            target_Time=start_time+datetime.timedelta(minutes=5)
+            timenow=datetime.datetime.strptime(time_list[k], "%Y-%m-%d %H:%M:%S")
+            if target_Time>=timenow:
+                delta_time=(target_Time-timenow).seconds
+                cal_flag=0
+                if delta_time<60:
+                    end_index=k
+                    cal_flag=1
+            else:        
+                end_index=k
+                cal_flag=1
+            
+            if cal_flag==1:
+                cal_count=cal_count+1
+                
+                df_deltaData=df_input_copy[max(start_index-100,0):end_index+1]
+                avg_cost=max(cost_min,df_result.loc[cal_count-1,'能耗[km/kwh]'])
+                #avg_cost=cost_min
+                df_span=real_odo(df_deltaData,avg_cost)
+                st_str=df_span['时间戳'].values[0]
+                st_time=datetime.datetime.strptime(st_str, "%Y-%m-%d %H:%M:%S")
+                last_time_s=datetime.datetime.strptime(last_time, "%Y-%m-%d %H:%M:%S")
+                if last_time_s>=st_time:
+                    last_time=last_time
+                    cross_odo=cross_odo
+                else:
+                    last_time=start_time
+                    cross_odo=0
+                delta_odo=df_span['累积里程[km]'].values[-1]-df_span[df_span['时间戳']==str(last_time)]['累积里程[km]'].values[0]
+                delta_energy=df_span['累积能量[Kwh]'].values[-1]-df_span[df_span['时间戳']==str(last_time)]['累积能量[Kwh]'].values[0]
+                df_result.loc[cal_count,'时间']=timenow
+                df_result.loc[cal_count,'时间']=str(df_result.loc[cal_count,'时间'])
+                df_result.loc[cal_count,'时间'] = df_result.loc[cal_count,'时间'].rpartition(':')[0]
+                
+                df_result.loc[cal_count,'累积里程[km]']=df_result.loc[cal_count-1,'累积里程[km]']+delta_odo-cross_odo
+                if delta_energy>0:
+                    df_result.loc[cal_count,'能耗[km/kwh]']=min(max(delta_odo/delta_energy,cost_min),cost_max)
+                else:
+                    df_result.loc[cal_count,'能耗[km/kwh]']=cost_min
+                
+                # df_span[df_span['纬度']>0][-1:]
+                if len(df_span[df_span['纬度']>0])>0:
+                    last_gps=df_span[df_span['纬度']>0][-1:]
+                    idx=last_gps.index.values[0]
+                    last_time=df_span.loc[idx,'时间戳']        
+                    cross_odo=df_span['累积里程[km]'].values[-1]-df_span.loc[idx,'累积里程[km]']
+                else:
+                    last_time=df_span['时间戳'].values[-1]        
+                    cross_odo=0
+                
+                start_index=end_index
+                start_time=timenow
+                cal_flag=0
+    
+        res=res.append(df_result)
+        input_res=input_res.append(df_input_copy)
+
+    l= len(input_res['时间戳'])
+    for k in range(l):
+        input_res.loc[k,'时间戳'] =input_res.loc[k,'时间戳'].rpartition(':')[0]
+
+    res.index=[i for i in range(len(res))]
+    res['日期']=[res.loc[i,'时间'][0:10] for i in range(len(res))]     
+
+    merge_res=pd.DataFrame()
+    for date_str in date_list:
+        input_res_copy=input_res[input_res['日期']==date_str]
+        input_res_copy=input_res_copy.drop_duplicates('时间戳')
+        res_copy=res[res['日期']==date_str]
+        df_merge=pd.DataFrame()
+        df_merge=pd.merge(input_res_copy,res_copy,left_on='时间戳',right_on='时间',how='right')
+        merge_res=merge_res.append(df_merge)
+
+    merge_res['sn']=sn
+    merge_res=merge_res[['时间戳','sn','总电流[A]','总电压[V]','SOC[%]','充电状态','纬度','经度','能耗[km/kwh]','累积里程[km]']]
+    merge_res=merge_res.rename(columns={'累积里程[km]':'每日累积里程[km]'})
+
+    return merge_res
+
+
+
+
+
+
+

LIB/MIDDLE/odo/DailyMileageEstimation/V1_0_4/main_daily_mileage.py

@@ -0,0 +1,42 @@
+from LIB.BACKEND import DBManager
+import calcul_mileage
+from LIB.MIDDLE.CellStateEstimation.Common import log
+import datetime
+import pandas as pd
+
+dbManager = DBManager.DBManager()
+
+now_time=datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')   #type: str
+now_time=datetime.datetime.strptime(now_time,'%Y-%m-%d %H:%M:%S')     #type: datetime
+start_time=now_time-datetime.timedelta(days=1)
+end_time=str(now_time)
+start_time=str(start_time)
+
+dataSOH = pd.read_excel('sn-20210903.xlsx',sheet_name='sn-20210903')
+fileNames = dataSOH['sn']
+fileNames = list(fileNames)
+l = len(fileNames)
+
+#log信息配置
+mylog=log.Mylog('log.txt','error')
+mylog.logcfg()
+
+
+for k in range(l):
+    try:
+        sn = fileNames[k]
+        df_data = dbManager.get_data(sn=sn, start_time=start_time, end_time=end_time, data_groups=['bms','gps'])
+        data_bms = df_data['bms']
+        data_gps = df_data['gps']
+        if sn[:4] in ['MGMC','UDO2']:
+            data_bms['总电流[A]'] = -data_bms['总电流[A]']
+
+        #...............每日累积里程............................................................................
+        if len(data_bms['时间戳'])>0 and len(data_gps['时间戳'])>0:
+            df_res = calcul_mileage.calcul_mileage(sn,data_bms,data_gps)
+            df_res.to_csv('Mileage_'+sn+'.csv')
+        
+    except Exception as e:
+        print(repr(e))
+        mylog.logopt(sn,e)
+        pass

LIB/MIDDLE/odo/DailyMileageEstimation/V1_0_4/main_daily_mileage_sn.py

@@ -0,0 +1,26 @@
+from LIB.BACKEND import DBManager
+import calcul_mileage
+from LIB.MIDDLE.CellStateEstimation.Common import log
+import datetime
+import pandas as pd
+
+dbManager = DBManager.DBManager()
+
+#now_time=datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')   #type: str
+now_time='2021-12-27 00:00:00'
+now_time=datetime.datetime.strptime(now_time,'%Y-%m-%d %H:%M:%S')     #type: datetime
+start_time=now_time-datetime.timedelta(days=1)
+end_time=str(now_time)
+start_time=str(start_time)
+
+sn='MGMCLN750N215N180'
+df_data = dbManager.get_data(sn=sn, start_time=start_time, end_time=end_time, data_groups=['bms','gps'])
+data_bms = df_data['bms']
+data_gps = df_data['gps']
+if sn[:4] in ['MGMC','UDO2']:
+    data_bms['总电流[A]'] = -data_bms['总电流[A]']
+
+        #...............每日累积里程............................................................................
+if len(data_bms['时间戳'])>0 and len(data_gps['时间戳'])>0:
+    df_res = calcul_mileage.calcul_mileage(sn,data_bms,data_gps)
+    df_res.to_csv('Mileage_'+sn+'.csv')

LIB/MIDDLE/odo/DailyMileageEstimation/V1_0_4/main_mileage.py

@@ -0,0 +1,39 @@
+
+from LIB.BACKEND import DBManager
+import calcul_mileage
+from LIB.MIDDLE.CellStateEstimation.Common import log
+import datetime
+import pandas as pd
+
+dbManager = DBManager.DBManager()
+
+dataSOH = pd.read_excel('sn-20210903.xlsx',sheet_name='sn-20210903')
+fileNames = dataSOH['sn']
+fileNames = list(fileNames)
+l = len(fileNames)
+
+#log信息配置
+mylog=log.Mylog('log.txt','error')
+mylog.logcfg()
+
+month='202110'
+
+for k in range(l):
+    try:
+        sn = fileNames[k]
+
+        df_data = dbManager.get_data(sn=sn, start_time='2021-12-01 00:00:00', end_time='2022-01-01 00:00:00', data_groups=['bms','gps'])
+        data_bms = df_data['bms']
+        data_gps = df_data['gps']
+        if sn[:4] in ['MGMC','UDO2']:
+            data_bms['总电流[A]'] = -data_bms['总电流[A]']
+
+        #...............每日累积里程............................................................................
+        if len(data_bms['时间戳'])>0 and len(data_gps['时间戳'])>0:
+            df_res = calcul_mileage.calcul_mileage(sn,data_bms,data_gps)
+            df_res.to_csv('Mileage_'+sn+'_'+month+'.csv')
+        
+    except Exception as e:
+        print(repr(e))
+        mylog.logopt(sn,e)
+        pass