Conformance Relations and Hyperproperties for Doping Detection in Time and Space — Case Study Details

The Jupyter notebook file and raw data is available here.

For our analysis, we use the Pandas data analysis and manipulation tool to represent and manipulate data. For visualisation, we show data graphically using the matplotlib library.

%matplotlib inline
import pandas as pd
pd.set_option('display.max_columns', None)
import matplotlib.pyplot as plt
import seaborn as sns 
sns.set(style="darkgrid")
import numpy as np

The results of our experiments are stored in individual CSV files. The raw data is loaded in variables with a leading underscore. For the main analysis part, we remove the data that was recordered before and after each test. A test starts at time point 0 and ends with time point 1179 (except for DoubleNEDC, which is twice as long). The trimmed data is stored in the appropriate variable without the leading underscore.

# Load Data
_nedc1 = pd.read_table('nedc1.csv')
_nedc2 = pd.read_table('nedc2.csv')
_perm1 = pd.read_table('perm1.csv')
_perm2 = pd.read_table('perm2.csv')
_sine1 = pd.read_table('sine1.csv')
_sine2 = pd.read_table('sine2.csv')
_double1 = pd.read_table('double1.csv')
_double2 = pd.read_table('double2.csv')

# Remove data points before and after each test cycle
def trimData(df, cycleLen=1180):
    firstRowIndex = df['Time'].gt(-0.5).idxmax()
    return df[firstRowIndex:firstRowIndex+cycleLen].reset_index(drop=True)

nedc1 = trimData(_nedc1)
nedc2 = trimData(_nedc2)
perm1 = trimData(_perm1)
perm2 = trimData(_perm2)
double1 = trimData(_double1, cycleLen=2*1180)
double2 = trimData(_double2, cycleLen=2*1180)
sine1 = trimData(_sine1)
sine2 = trimData(_sine2)

# Constants
VELOCITY_KEY = 'Velocity ECU'
NOX_KEY = 'mdot_NOx DC'

Experiments are represented by tables. Retimings are applied on two experiments, resulting in a merged single table. The following functions provide auxiliary tools to merge and separate (merged) tables.

• All experiments use the same column labels. When merging tables (e.g., when applying retimings) we use a renaming functions to rename the NEDC columns by adding ' Ref' to each column name. This allows to identify after the merge to which of the original experiments data belongs to. When splitting a merged table, the renaming is reversed. Functions rename and renameBack allow easy generation of renaming functions (and their inverse), that rename by appending a string. renameRef and renameBackRef are the renaming functions that add ' Ref' to the names, and which we use in the following.
• joinExperiments creates the merged table of two single tables (without applying a retiming).
• sortFrame is a convenient way to permanently reorder the order of table rows through sorting by a particular column.
• splitExperiments splits merged tables into two single tables. The function assumes that the order of the columns is as it would be for tables merged using joinExperiments.

# Merging and splitting experiments

def rename(additionalString):
    return (lambda x: str(x) + str(additionalString))

def renameBack(additionalString):
    def f(s):
        if s[-len(additionalString):] == additionalString:
            return s[:-len(additionalString)]
        else:
            raise Exception('String was not renamed by adding "additionalString" previously. Cannot rename back!')
    return f

renameRef = rename(' Ref')
renameBackRef = renameBack(' Ref')


def joinExperiments(df1, df2, f_rename=renameRef):
    df1 = df1.rename(f_rename, axis='columns')
    joined = pd.concat([df1, df2], axis=1)
    return joined

def sortFrame(df, key, dropOldIndex=False):
    return df.sort_values(key).reset_index(drop=dropOldIndex)

def splitExperiments(joinedDf):
    if len(joinedDf.columns) % 2 != 0:
        raise Exception('Joined DFs must have even number of columns!')
    i = int(len(joinedDf.columns)/2)
    renamedColumns = joinedDf.columns[:i]
    originalColumns = joinedDf.columns[i:]
    def assertCorrectColumnNames(o, r):
        if r.find(o) == 0 and len(r) > len(o):
            return True
        else:
            raise Exception("Column names don't match!")
    list(map(assertCorrectColumnNames, originalColumns, renamedColumns))
    df1 = joinedDf.drop(columns=originalColumns)
    df1.columns = originalColumns
    df2 = joinedDf.drop(columns=renamedColumns)
    return (df1, df2)

Like all experiments, we executed NEDC twice. For simplicity, we compare all non-NEDC experiments only to one reference NEDC, which is the average of NEDC-1 and NEDC-2. Unfortunately, due to technical problems no OBD data (velocity, engine speed) is available for NEDC-1, so we take only data from NEDC-2. The averaged NEDC will be available as nedc in the following cells.

# Take avergae of NEDC1 and NEDC2 data
def averageCycles(df1, df2):
    if (df1.columns != df2.columns).any():
        raise Exception('Can only average cycles with identical columns!')
    cols = df1.columns
    name = rename(' (Cycle 2)')
    joined = joinExperiments(df1, df2, name)
    avg = joined.apply(lambda row: pd.Series(map(lambda col: row[col] if col == 'Time' else (row[col] + row[name(col)])/2.0, cols), index=cols, name=row.name), axis=1)
    return avg

nedc = averageCycles(nedc1, nedc2)
# There is no OBD data for NEDC1, hence we keep the data from NEDC2; averaging results in incorrect values
for col in ['Velocity ECU', 'Engine Speed']:
    nedc[col] = nedc2[col]

For the evaluation of our experiments, we accumulate NOx values over 1180 seconds. The accumulated value is given in milligrams per kilometre. Function integrateData computes the sum of the values in a column, w.r.t. a given step function. It can compute the sum of all NOx values (which are given in grams) and can compute the total distance that the test vehicle travelled during the test execution. For the computation of the total distance, the speed values must be converted from km/h to m/s, and the result from m in km. Function accumulatedNOxPerTotalDistance is doing these conversions and computations.

# Integrating data
def integrateData(df, key, step=lambda prev, this: prev + this):
    df = df.copy()
    previous = 0.0
    for i in range(0, len(df)):
        previous = step(previous, df.loc[i, key])
        df.loc[i, key] = previous
    return df

def accumulatedNOxPerTotalDistance(df):
    distDf = integrateData(df, VELOCITY_KEY, lambda prev, v: prev + v/3.6)  # convert speed to m/s
    totalDistance = (distDf.loc[len(distDf)-1, VELOCITY_KEY])/1000  # convert distance to km
    return integrateData(df, NOX_KEY, lambda prev, v: prev + v*1000/totalDistance)

For later comparisons, we need accumulated NOx values and prepare the appropriate values here.

nedcAccumulated = accumulatedNOxPerTotalDistance(nedc)
perm1Accumulated = accumulatedNOxPerTotalDistance(perm1)
perm2Accumulated = accumulatedNOxPerTotalDistance(perm2)
double1Accumulated = accumulatedNOxPerTotalDistance(double1)
double2Accumulated = accumulatedNOxPerTotalDistance(double2)
sine1Accumulated = accumulatedNOxPerTotalDistance(sine1)
sine2Accumulated = accumulatedNOxPerTotalDistance(sine2)

Another utility function we will need is getMaximumValueError, which, for a retiming given as a merged table, computes the row-wise maximum absolute difference of values in column key and the appropriate renamed column. It returns the largest of these differences.

# Find maximum value error
def getMaximumValueError(retimingDf, key, f_rename=renameRef):
    return max(retimingDf.apply(lambda row: abs(row[key] - row[f_rename(key)]), axis=1))

In order to show test results graphically, we use function plotCycles to compare velocity and NOx values of two experiments. The blue and orange lines show speed and NOx values, respectively, for the first table. Green and red are used for the second. Time ticks, NOx ticks and time values can be modified if necessary.

NOX_TICKS_DEFAULT = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1]
NOX_TICKS_ACCUMULATED = [0.0, 80.0, 182.0, 316.0, 483.0, 632.0]

TIME_TICKS_DEFAULT = [0, 200,400,600,800,1000,1179]
TIME_TICKS_DOUBLE = [0, 200,400,600,800,1000,1180, 2359]

def plotCycles(df1, df2, title, timeTicks=TIME_TICKS_DEFAULT, noxTicks=NOX_TICKS_DEFAULT, timeValues=None):
    if isinstance(timeValues, type(None)):
        timeValues = df1['Time']
    speed1 = df1[VELOCITY_KEY]
    speed2 = df2[VELOCITY_KEY]
    nox1 = df1[NOX_KEY]
    nox2 = df2[NOX_KEY]
    
    fig, ax1 = plt.subplots(nrows=1)
    
    ax1.set_title(title)
    ax1.set_xlabel("Time")
    ax1.set_ylabel("Speed")
    ax1.plot(timeValues, speed1, color='tab:blue')
    ax1.plot(timeValues, speed2, color='tab:green')
    ax1.margins(x=0, y=0)
    ax1.get_yaxis().tick_left()
    ax1.tick_params(grid_linestyle='', labelleft=True, labelright=False, left=True, right=False)
    ax1.set_xticks(timeTicks)

    ay1 = ax1.twinx()
    ay1.set_ylabel("NOx")
    ay1.plot(timeValues, nox1, color='tab:orange')
    ay1.plot(timeValues, nox2, color='tab:red')
    ay1.set_yticks(noxTicks)
    ay1.tick_params(grid_linestyle='', labelleft=False, labelright=True, left=False, right=True)
    ay1.get_yaxis().tick_right()
    ay1.set_xticks(timeTicks)
    ay1.margins(x=0)

    fig.set_size_inches(16, 5)
    fig.set_dpi(144.0)

For example, the following command plots a comparison of NEDC and PermNEDC-1 (with NOx given in g/s).

plotCycles(nedc, perm1, "Perm1 and averaged NEDC")

Retiming Experiments

In the following, we will investigate each of the three experiment types mentioned in the paper: NEDC permutations, lengthened test cycles and hybrid conformance.

All functions below assume that the data is sampled by 1Hz, which is asserted with the following function.

# For correctness of the following algorithms, it is necessary that the data is sampled by 1Hz
def assertEquiDistantDiscreteTime(df):
    time = df['Time']
    for i in range(1, len(time)):
        if (abs(time[i] - time[i-1] - 1.0)) > 10E-9:
            raise Exception('Data not sampled as expected!')

assertEquiDistantDiscreteTime(nedc1)
assertEquiDistantDiscreteTime(nedc2)
assertEquiDistantDiscreteTime(perm1)
assertEquiDistantDiscreteTime(perm2)
assertEquiDistantDiscreteTime(double1)
assertEquiDistantDiscreteTime(double2)
assertEquiDistantDiscreteTime(sine1)
assertEquiDistantDiscreteTime(sine2)

NEDC Permutations

The PermNEDC cycle is generated by splitting the UDC subcycle into seven pieces. The length of each piece is given by euc0 to euc6 below. The length of EUDC is given in extra. We encode the NEDC and PermNEDC in lists, that store sequences of index and length of UDC pieces. The original EUC is stored in ops (length of the pieces) and [0,1,2,3,4,5,6] (the index of the pieces). The four UDC repetitions are encoded by nedcOps and nedcOpsi. Permutation NEDC is defined by the permuted ops1234 and opsi1234 lists given below.

The functions nedc2perm and perm2nedc use this encoding to map a time point t of NEDC to the appropriate t' of PermNEDC (i.e., t has been shifted to the new position t'). Similarly, perm2nedc computes the inverse of nedc2perm.

# r_p
euc0 = 6+5
euc1 = 0+4+8+2+3
euc2 = 16+5
euc3 = 0+5+2+5+24+8+3
euc4 = 16+5
euc5 = 0+5+2+9+2+8+12+8+13+2+7+3
euc6 = 7
euc = euc0 + euc1 + euc2 + euc3 + euc4 + euc5 + euc6
extra = 20+5+2+9+2+8+2+13+50+4+4+69+13+50+35+30+20+10+16+8+10+20
nedcLen = 4*euc + extra

ops = [euc0, euc1, euc2, euc3, euc4, euc5, euc6]

# NEDC ops
nedcOpsi = list(range(28))
nedcOps = ops + ops + ops + ops

# PermNEDC ops
ops1 = [euc0, euc3, euc2, euc5, euc4, euc1, euc6]
opsi1 = [0, 3, 2, 5, 4, 1, 6]
ops2 = [euc0, euc5, euc2, euc1, euc4, euc3, euc6]
opsi2 = [0+7, 5+7, 2+7, 1+7, 4+7, 3+7, 6+7]
ops3 = [euc0, euc5, euc2, euc3, euc4, euc1, euc6]
opsi3 = [0+14, 5+14, 2+14, 3+14, 4+14, 1+14, 6+14]
ops4 = [euc0, euc3, euc2, euc1, euc4, euc5, euc6]
opsi4 = [0+21, 3+21, 2+21, 1+21, 4+21, 5+21, 6+21]
ops1234 = ops1 + ops2 + ops3 + ops4
opsi1234 = opsi1 + opsi2 + opsi3 + opsi4

import functools
import operator

foldl = lambda func, acc, xs: functools.reduce(func, xs, acc)

totalOps1234 = foldl(lambda a, e: a + [a[-1]+e], [0], ops1234)
totalNedcOps = foldl(lambda a, e: a + [a[-1]+e], [0], nedcOps)

def getOpsLength(idx):
    return ops[idx % 7]

def findOpsIndexForTime(t, totalOps):
    for i in range(1, len(totalOps)):
        if t < totalOps[i]:
            return (i-1, t-totalOps[i-1])
    raise Exception("Time beyond ops")

def findOpsIndexForEUCIndex(idx, opsi):
    for i in range(len(opsi)):
        if opsi[i] == idx:
            return i
    raise Exception("Ops index not found")
    
def nedc2perm(t):
    if t >= 4*euc:
        # Extra urban -> identity
        return t
    
    (nedcOpsIndex, offset) = findOpsIndexForTime(t, totalNedcOps)
    nedcEucIndex = nedcOpsi[nedcOpsIndex] # in case of NEDC trivial, because identity
    permOpsIndex = findOpsIndexForEUCIndex(nedcEucIndex, opsi1234)
    return totalOps1234[permOpsIndex] + offset
    
def perm2nedc(t):
    if t >= 4*euc:
        # Extra urban -> identity
        return t
    
    (permOpsIndex, offset) = findOpsIndexForTime(t, totalOps1234)
    permEucIndex = opsi1234[permOpsIndex]
    nedcOpsIndex = findOpsIndexForEUCIndex(permEucIndex, nedcOpsi)
    return totalNedcOps[nedcOpsIndex] + offset

# Verify correctness
for i in range(780):
    if nedc2perm(perm2nedc(i)) != i:
        raise Exception("Oh no! Something is wrong at index " + str(i))

The pair of retiming functions (nedc2perm, perm2nedc) represents the retiming $\mathsf{Ret}_p$. The functions doPermRetiming1 and doPermRetiming2 below, apply the retiming functions. Retiming functions are applied to two experiment tables and result in a single, merged, table. To show the result visually, splitExperiments is used to separate the original and retimed version of the original tables.

For example, doPermRetiming1(nedc, perm1) keeps nedc untouched, but reorders perm1, such that it looks like the NEDC. The result is shown in the first plot. doPermRetiming2(nedc, perm1) does the opposite. It leaves perm1 unchanged but reorders nedc to match the time structure of perm1. The result is shown in third and fourth plot.

# Applying retiming function `nedc2perm` to transform PermNEDC results back into NEDC order
def doPermRetiming1(nedc, perm, nedcRenaming=renameRef):
    nedc = nedc.rename(nedcRenaming, axis='columns')
    retiming = nedc.apply(lambda row: pd.concat([row, perm.loc[nedc2perm(row.name)]]), axis=1)
    return retiming

# Applying retiming function `perm2nedc` to transform NEDC results into PermNEDC order
def doPermRetiming2(nedc, perm, nedcRenaming=renameRef):
    nedc = nedc.rename(nedcRenaming, axis='columns')
    retiming = perm.apply(lambda row: pd.concat([nedc.loc[perm2nedc(row.name)], row]), axis=1)
    return retiming

perm1Retiming1 = doPermRetiming1(nedc, perm1)
recoveredPerm1 = splitExperiments(perm1Retiming1)[1]
recoveredPerm1Accumulated = accumulatedNOxPerTotalDistance(recoveredPerm1)
plotCycles(nedcAccumulated, recoveredPerm1Accumulated, "Perm1 and averaged NEDC (nedc2perm)", noxTicks=NOX_TICKS_ACCUMULATED)

perm2Retiming1 = doPermRetiming1(nedc, perm2)
recoveredPerm2 = splitExperiments(perm2Retiming1)[1]
recoveredPerm2Accumulated = accumulatedNOxPerTotalDistance(recoveredPerm2)
plotCycles(nedcAccumulated, recoveredPerm2Accumulated, "Perm2 and averaged NEDC (nedc2perm)", noxTicks=NOX_TICKS_ACCUMULATED)

nedcPerm1Retiming = doPermRetiming2(nedc, perm1)  # Doesn't matter if we pick perm1 or perm2 as both use the same retiming function
permutedNEDC = splitExperiments(nedcPerm1Retiming)[0]
permutedNEDCAccumulated = accumulatedNOxPerTotalDistance(permutedNEDC)
plotCycles(permutedNEDCAccumulated, perm1Accumulated, "Perm1 and averaged NEDC (perm2nedc)", noxTicks=NOX_TICKS_ACCUMULATED, timeValues=perm1['Time'])

nedcPerm2Retiming = doPermRetiming2(nedc, perm2)  # Provide for future applications
plotCycles(permutedNEDCAccumulated, perm2Accumulated, "Perm2 and averaged NEDC (perm2nedc)", noxTicks=NOX_TICKS_ACCUMULATED, timeValues=perm2['Time'])