2025-AKO/v5_class.py

# Csar Fdez, UdL, 2025
# Changes from v1:   Normalization 
# IN v1, each failure type has its own normalization pars (mean and stdevs)
# In v2, mean and stdev is the same for all data
# v3.py trains the models looping in TIME_STEPS (4,8,12,16,20,24,....) finding the optimal Threshold factor

#  Derived from v3_class, derived from v3.py with code from v1_multifailure.py
#  This code don't train for multiple time steps !!

#  partial and total blocked condenser merged in one class.
#  Construction of train and test sets changed. Now is done by days

import pandas as pd
import matplotlib.pyplot as plt
import datetime
import numpy as np
import keras
import os.path
from keras import layers
from optparse import OptionParser
import copy
import pickle


parser = OptionParser()
parser.add_option("-t", "--train", dest="train", help="Trains the models (false)", default=False, action="store_true")
parser.add_option("-n", "--timesteps", dest="timesteps", help="TIME STEPS ", default=12)
#parser.add_option("-f", "--thresholdfactor", dest="TF", help="Threshold Factor ", default=1.4)
# threshold makes no sense when classifying, becaues we apply many models and decide class for the less MSE

(options, args) = parser.parse_args()


# data files arrays. Index:
# 0.  No failure
# 1.  Blocked evaporator
# 2.   Full Blocked condenser
# 3.   Partial Blocked condenser
# 4   Fan condenser not working
# 5.  Open door


NumberOfFailures=3  # So far, we have only data for the first 4 types of failures
datafiles=[[],[]]   # 0 for train,  1 for test
for i in range(NumberOfFailures+1):
    datafiles[0].append([])
    datafiles[1].append([])

# Next set of ddata corresponds to Freezer, SP=-26
datafiles[0][0]=['2024-08-07_5_','2024-08-08_5_','2025-01-25_5_','2025-01-26_5_'] 
datafiles[1][0]=['2025-01-27_5_','2025-01-28_5_'] 

datafiles[0][1]=['2024-12-11_5_', '2024-12-12_5_','2024-12-13_5_'] 
datafiles[1][1]=['2024-12-14_5_','2024-12-15_5_'] 

datafiles[0][2]=['2024-12-18_5_','2024-12-21_5_','2024-12-22_5_','2024-12-23_5_','2024-12-24_5_'] 
datafiles[1][2]=['2024-12-19_5_','2024-12-25_5_','2024-12-26_5_'] 


datafiles[0][3]=['2024-12-28_5_','2024-12-29_5_','2024-12-30_5_'] 
datafiles[1][3]=['2024-12-31_5_','2025-01-01_5_'] 

#r1s5 supply air flow temperature
#r1s1 inlet evaporator temperature
#r1s4 condenser outlet

# VAriables r1s4 and pa1 apiii  may not exists in cloud controlers


features=['r1 s1','r1 s4','r1 s5','pa1 apiii']
features=['r1 s1','r1 s4','r1 s5']
features=['r1 s5']
# Feature combination suggested by AKO
#features=['r1 s1','r1 s4','r1 s5','pa1 apiii']
#features=['r1 s1','r1 s4','r1 s5']
#features=['r1 s1','r1 s5','pa1 apiii']
#features=['r1 s5','pa1 apiii']
features=['r1 s1','r1 s5']
#features=['r1 s5']



featureNames={}
featureNames['r1 s1']='$T_{evap}$'
featureNames['r1 s4']='$T_{cond}$'
featureNames['r1 s5']='$T_{air}$'
featureNames['pa1 apiii']='$P_{elec}$'

unitNames={}
unitNames['r1 s1']='$(^{o}C)$'
unitNames['r1 s4']='$(^{o}C)$'
unitNames['r1 s5']='$(^{o}C)$'
unitNames['pa1 apiii']='$(W)$'


#features=['r1 s1','r1 s2','r1 s3','r1 s4','r1 s5','r1 s6','r1 s7','r1 s8','r1 s9','r1 s10','r2 s1','r2 s2','r2 s3','r2 s4','r2 s5','r2 s6','r2 s7','r2 s8','r2 s9','pa1 apiii','tc s1','tc s2']

#features=['r2 s2', 'tc s1','r1 s10','r1 s6','r2 s8']

NumFeatures=len(features)

df_list=[[],[]]
for i in range(NumberOfFailures+1):
    df_list[0].append([])
    df_list[1].append([])

for i in range(NumberOfFailures+1):
    dftemp=[]
    for f in datafiles[0][i]:
        print("                 ", f)
        df1 = pd.read_csv('./data/'+f+'.csv')
        dftemp.append(df1)
    df_list[0][i]=pd.concat(dftemp)

for i in range(NumberOfFailures+1):
    dftemp=[]
    for f in datafiles[1][i]:
        print("                 ", f)
        #df1 = pd.read_csv('./data/'+f+'.csv', parse_dates=['datetime'], dayfirst=True, index_col='datetime')
        df1 = pd.read_csv('./data/'+f+'.csv')
        dftemp.append(df1)
    df_list[1][i]=pd.concat(dftemp)



# subsampled to 5'  =  30 * 10"
# We consider smaples every 5' because in production, we will only have data at this frequency
subsamplingrate=30

dataframe=[[],[]]
for i in range(NumberOfFailures+1):
    dataframe[0].append([])
    dataframe[1].append([])

for i in range(NumberOfFailures+1):
    datalength=df_list[0][i].shape[0]
    dataframe[0][i]=df_list[0][i].iloc[range(0,datalength,subsamplingrate)][features]
    dataframe[0][i].reset_index(inplace=True,drop=True)
    dataframe[0][i].dropna(inplace=True)
for i in range(NumberOfFailures+1):
    datalength=df_list[1][i].shape[0]
    dataframe[1][i]=df_list[1][i].iloc[range(0,datalength,subsamplingrate)][features]
    dataframe[1][i].reset_index(inplace=True,drop=True)
    dataframe[1][i].dropna(inplace=True)


# Train data is first 2/3 of data
# Test data is: last 1/3 of data 
dataTrain=[]
dataTest=[]
for i in range(NumberOfFailures+1):
    dataTrain.append(dataframe[0][i].values)
    dataTest.append(dataframe[0][i])

# Calculate means and stdev
a=dataTrain[0]
for i in range(1,NumberOfFailures+1):
    a=np.vstack((a,dataTrain[i]))

means=a.mean(axis=0) 
stdevs=a.std(axis=0)
def normalize2(train,test):
    return( (train-means)/stdevs, (test-means)/stdevs )

dataTrainNorm=[]
dataTestNorm=[]
for i in range(NumberOfFailures+1):
    dataTrainNorm.append([])
    dataTestNorm.append([])

for i in range(NumberOfFailures+1):
    (dataTrainNorm[i],dataTestNorm[i])=normalize2(dataTrain[i],dataTest[i])

def plotData():    
    fig, axes = plt.subplots(
        nrows=NumberOfFailures+1, ncols=2, figsize=(15, 20), dpi=80, facecolor="w", edgecolor="k",sharex=True
    )
    for i in range(NumberOfFailures+1):
        axes[i][0].plot(np.concatenate((dataTrainNorm[i][:,0],dataTestNorm[i][:,0])),label="Fail "+str(i)+",  feature 0")
        axes[i][1].plot(np.concatenate((dataTrainNorm[i][:,1],dataTestNorm[i][:,1])),label="Fail "+str(i)+",  feature 1")
    #axes[1].legend()
    #axes[0].set_ylabel(features[0])
    #axes[1].set_ylabel(features[1])
    plt.show()

#plotData()
#exit(0)


NumFilters=64
KernelSize=7
DropOut=0.2
ThresholdFactor=1.4
def create_sequences(values, time_steps):
    output = []
    for i in range(len(values) - time_steps + 1):
        output.append(values[i : (i + time_steps)])
    return np.stack(output)



def listToString(l):
    r=''
    for i in l:
        r+=str(i)
    return(r.replace(' ',''))


model=[]
modelckpt_callback =[]
es_callback =[]
path_checkpoint=[]

timesteps=int(options.timesteps)
x_train=[]
for i in range(NumberOfFailures+1):
    x_train.append(create_sequences(dataTrainNorm[i],timesteps))
    model.append([])
    model[i] = keras.Sequential(
        [
            layers.Input(shape=(x_train[i].shape[1], x_train[i].shape[2])),
            layers.Conv1D(
                filters=NumFilters,
                kernel_size=KernelSize,
                padding="same",
                strides=2,
                activation="relu",
            ),
            layers.Dropout(rate=DropOut),
            layers.Conv1D(
                filters=int(NumFilters/2),
                kernel_size=KernelSize,
                padding="same",
                strides=2,
                activation="relu",
            ),
            layers.Conv1DTranspose(
                filters=int(NumFilters/2),
                kernel_size=KernelSize,
                padding="same",
                strides=2,
                activation="relu",
            ),
            layers.Dropout(rate=DropOut),
            layers.Conv1DTranspose(
                filters=NumFilters,
                kernel_size=KernelSize,
                padding="same",
                strides=2,
                activation="relu",
            ),
            layers.Conv1DTranspose(filters=x_train[i].shape[2], kernel_size=KernelSize, padding="same"),
        ]
    )
    model[i].compile(optimizer=keras.optimizers.Adam(learning_rate=0.001), loss="mse")
    model[i].summary()
    path_checkpoint.append("model_class_v5_"+str(i)+"_"+str(timesteps)+listToString(features)+"_checkpoint.weights.h5")
    es_callback.append(keras.callbacks.EarlyStopping(monitor="val_loss", min_delta=0, patience=15))
    modelckpt_callback.append(keras.callbacks.ModelCheckpoint( monitor="val_loss", filepath=path_checkpoint[i], verbose=1, save_weights_only=True, save_best_only=True,))


if options.train:
    history=[]    
    for i in range(NumberOfFailures+1):
        history.append(model[i].fit( x_train[i], x_train[i], epochs=400, batch_size=128, validation_split=0.3, callbacks=[  es_callback[i], modelckpt_callback[i]      ],))

        x_train_pred=model[i].predict(x_train[i])
else:
    for i in range(NumberOfFailures+1):
        model[i].load_weights(path_checkpoint[i])



# Let's plot some features

colorline=['black','violet','lightcoral','cyan','lime','grey']
colordot=['grey','darkviolet','red','blue','green','black']

#featuresToPlot=['r1 s1','r1 s2','r1 s3','pa1 apiii']
featuresToPlot=features

indexesToPlot=[]
for i in featuresToPlot:
    indexesToPlot.append(features.index(i))



#   2nd scenario. Go over anomalies and classify it by less error
datalist=[dataTestNorm[0],dataTestNorm[1],dataTestNorm[2],dataTestNorm[3]]
x_test=create_sequences(datalist[0],int(options.timesteps))
for i in range(1,len(datalist)):
    x_test=np.vstack((x_test,create_sequences(datalist[i],int(options.timesteps))))

# Define ranges for plotting in different colors
testRanges=[]
r=0
for i in range(len(datalist)):
    testRanges.append([r,r+datalist[i].shape[0]-int(options.timesteps)])
    r+=datalist[i].shape[0]-int(options.timesteps)

testClasses=[0,1,2,3]

if not len(testClasses)==len(testRanges):
    print("ERROR:  testClasses and testRanges must have same length")
    exit(0)

x_test_predict=[]
for m in range(NumberOfFailures+1):
    x_test_predict.append(model[m].predict(x_test))

x_test_predict=np.array((x_test_predict))
test_mae_loss =[]
for m in range(NumberOfFailures+1):
    test_mae_loss.append(np.mean(np.abs(x_test_predict[m,:,:,:] - x_test), axis=1))

test_mae_loss=np.array((test_mae_loss))
test_mae_loss_average=np.mean(test_mae_loss,axis=2)  # average over features
classes=np.argmin(test_mae_loss_average,axis=0)

x=[]
y=[]
for j in range(NumberOfFailures+1):
    x.append([])
    y.append([])
for j in range(NumberOfFailures+1):
    for k in range(testRanges[j][0],testRanges[j][1]):
        if not  classes[k]==testClasses[j]:
            x[classes[k]].append(k)
            y[classes[k]].append(x_test[k,0,indexesToPlot[0]]*stdevs[0]+means[0])


def plotData4():
    NumFeaturesToPlot=len(indexesToPlot)
    plt.rcParams.update({'font.size': 16})
    fig, axes = plt.subplots(
        nrows=NumFeaturesToPlot, ncols=1, figsize=(15, 10), dpi=80, facecolor="w", edgecolor="k",sharex=True
    )
    for i in range(NumFeaturesToPlot):
        init=0
        end=testRanges[0][1]
        for j in range(NumberOfFailures+1):
            if NumFeaturesToPlot==1:
                axes.plot(range(init,end),x_test[testRanges[j][0]:testRanges[j][1],0,indexesToPlot[i]]*stdevs[i]+means[i],label="Class "+str(j), color=colorline[j],linewidth=1)
            else:
                axes[i].plot(range(init,end),x_test[testRanges[j][0]:testRanges[j][1],0,indexesToPlot[i]]*stdevs[i]+means[i],label="Class "+str(j), color=colorline[j],linewidth=1)
            if j<NumberOfFailures:
                init=end
                end+=(testRanges[j+1][1]-testRanges[j+1][0])

            #if i==0:
            #    axes[0].plot(x[j],y[j] ,color=colordot[j],marker='.',markersize=10,linewidth=0,label="Fail detect  type "+str(j) )



        s=''
        s+=featureNames[features[indexesToPlot[i]]]
        s+=' '+unitNames[features[indexesToPlot[i]]]
        if NumFeaturesToPlot==1:
            axes.set_ylabel(s)
            axes.grid()
        else:
            axes[i].set_ylabel(s)
            axes[i].grid()

    for j in range(NumberOfFailures+1):
        if NumFeaturesToPlot==1:
            axes.plot(x[j],y[j] ,color=colordot[j],marker='.',markersize=10,linewidth=0,label="Fail detect  type "+str(j) )
        else:
            axes[0].plot(x[j],y[j] ,color=colordot[j],marker='.',markersize=10,linewidth=0,label="Fail detect  type "+str(j) )
            
    if NumFeaturesToPlot==1:
        axes.legend(ncol=4,loc=(0.1,0.98))
    else:
        axes[0].legend(ncol=4,loc=(0.1,0.98))

        
    #axes[NumFeaturesToPlot-1].set_xlabel("Sample number")
    plt.show()

def whichClass(k,ranges):
    for i in range(NumberOfFailures+1):
        if k in range(ranges[i][0],ranges[i][1]):
            return(i)
    print("Error:  Class not exists")
    exit(0)        

##   It remains to implemenent anomaly metrics for each failure type
def anomalyMetric(classes,testranges,testclasses):  
    # FP, TP: false/true positive
    # TN, FN: true/false negative
    # Sensitivity (recall): probab failure detection if data is fail: TP/(TP+FN)
    # Precision: Rate of positive results:  TP/(TP+FP)  
    # F1-score: predictive performance measure: 2*Precision*Sensitity/(Precision+Sensitity)
    TP=np.zeros(NumberOfFailures+1)
    FP=np.zeros(NumberOfFailures+1)
    FN=np.zeros(NumberOfFailures+1)
    Sensitivity=np.zeros(NumberOfFailures+1)
    Precision=np.zeros(NumberOfFailures+1)
    for i in range(len(testranges)):
        for k in range(testranges[i][0],testranges[i][1]):
            if classes[k]==testclasses[i]:
                TP[i]+=1
            else:
                FP[i]+=1
    for k in range(testranges[NumberOfFailures][1]):
        for i in range(len(testranges)):
            classK=whichClass(k,testranges)
            if not classK==testClasses[i]:
                if not classes[k]==classK:
                    FN[classes[k]]+=1

    for i in range(NumberOfFailures+1):
        Sensitivity[i]=TP[i]/(TP[i]+FN[i])
        Precision[i]=TP[i]/(TP[i]+FP[i])
    S=Sensitivity.mean()
    P=Precision.mean()
    F1=2*S*P/(S+P)
    print("Sensitivity: ",Sensitivity) 
    print("S: ",S) 
    print("Precision: ",Precision) 
    print("P: ",P) 
    print("F1-Score: ",F1)

anomalyMetric(classes,testRanges,testClasses)
plotData4()
exit(0)