<jupyter_start><jupyter_text>EngLyser Research <NAME> | http://phsheth.github.io<jupyter_code>import numpy as np
from scipy.fftpack import fft, ifft
from scipy.signal import butter, lfilter
import scipy as sp
import matplotlib.pyplot as plt
from pylab import *

from time import perf_counter
import pandas as pd

from sklearn.cross_validation import train_test_split

from sklearn.ensemble import RandomForestClassifier


%matplotlib qt
def thetagen(cycles,thetastepsize, forwardoffset):
    thetasingle = np.array(np.arange(0,720,thetastepsize))
    thetarel = np.tile(thetasingle,cycles)
    #thetacum = np.cumsum(thetarel)
    thetacum = np.array(np.arange(0,(720*cycles),thetastepsize))
    #thetacum = thetacum
    return thetasingle, thetarel, thetacum
def gasforcegen(gasforceraw, cycles, thetastepsize,forwardoffset,misfiring):
    thetasingle_deg, thetarel_deg, thetacum_deg = thetagen(cycles,thetastepsize,forwardoffset)
    gasforcerawdata = pd.read_csv(gasforceraw,sep=',',header=None)
    funcakimainterp_gf = sp.interpolate.Akima1DInterpolator(gasforcerawdata[0],gasforcerawdata[1])
    gasforce_theta_deg = funcakimainterp_gf(thetasingle_deg)
    gasforce_theta_deg_cycles = np.tile(gasforce_theta_deg,cycles)
    gasforce_theta_deg_cycles = gasforce_theta_deg_cycles*misfiring
    return gasforce_theta_deg_cycles, thetarel_deg,thetacum_deg
def torquegen(l1_m, l2_m,
              mrod, Izz, lg, mp, engine_speed_rpm, rp_m, 
              gasforcefile, cycles, thetastepsize,forwardoffset,misfiring):
    gasforce_theta_deg_cycles, thetarel_deg,thetacum_deg = gasforcegen(gasforcefile, cycles, thetastepsize,forwardoffset,misfiring)
    #gasforce_theta_deg_cycles = 35*gasforce_theta_deg_cycles/max(gasforce_theta_deg_cycles)
    thetacum_rad = thetacum_deg*np.pi/180
    theta_d = engine_speed_rpm*2*np.pi/60
    theta_dd = 0
    beta = np.arcsin(l1_m*np.sin(thetacum_rad)/l2_m)
    con1 = (1-(np.square(l1_m*np.sin(thetacum_rad)/l2_m)))
    beta_d = (theta_d*l1_m*np.cos(thetacum_rad))/(l2_m*np.sqrt(con1))
    
    a_beta_dd = l1_m*((theta_dd*np.cos(thetacum_rad))-((np.square(theta_d))*np.sin(thetacum_rad)))
    b_beta_dd = l2_m*np.sqrt(con1)
    c_beta_dd = np.square(theta_d)*np.power(l1_m,3)*np.square(np.cos(thetacum_rad))*np.sin(thetacum_rad)
    d_beta_dd = np.power(l2_m,3)*con1
    beta_dd= (a_beta_dd/b_beta_dd)+(c_beta_dd/d_beta_dd)

    a_vpis = (l1_m*theta_d*np.sin(thetacum_rad))
    b_vpis = (np.square(l1_m)*theta_d*np.sin(2*thetacum_rad))
    c_vpis = 2*l2_m*np.sqrt(con1)
    vpis = -a_vpis-(b_vpis/c_vpis)
    
    a_arx = l1_m*theta_dd*np.sin(thetacum_rad)
    b_arx = l1_m*(np.square(theta_d))*np.cos(thetacum_rad)
    c_arx = theta_dd*lg*np.square(l1_m)*np.sin(2*thetacum_rad)
    d_arx = 2*np.square(l2_m)*np.sqrt(con1)
    e_arx = np.square(l2_m)*(2*con1)
    f_arx = np.square(theta_d*l1_m)*lg
    g_arx = 2*np.cos(2*thetacum_rad)*np.sqrt(con1)
    h_arx = np.square(l1_m*np.sin(2*thetacum_rad))
    i_arx = 2*np.square(l2_m)*np.sqrt(con1)
    arx = -a_arx-b_arx-(c_arx/d_arx)-(1/e_arx)*f_arx*(g_arx+(h_arx/i_arx))
    
    a_ary = l1_m*theta_dd*np.cos(thetacum_rad)
    b_ary = l1_m*np.square(theta_d)*np.sin(thetacum_rad)
    c_ary = lg*l1_m*theta_dd*np.cos(thetacum_rad)/l2_m
    d_ary = lg*l1_m*theta_d*np.sin(thetacum_rad)/l2_m
    ary = a_ary - b_ary - c_ary+d_ary
    
    a_apx = l1_m*theta_dd*np.sin(thetacum_rad)
    b_apx = l1_m*(np.square(theta_d))*np.cos(thetacum_rad)
    c_apx = theta_dd*np.square(l1_m)*np.sin(2*thetacum_rad)
    d_apx = 2*l2_m*np.sqrt(con1)
    e_apx = l2_m*(2*con1)
    f_apx = np.square(theta_d*l1_m)
    g_apx = 2*np.cos(2*thetacum_rad)*np.sqrt(con1)
    h_apx = np.square(l1_m*np.sin(2*thetacum_rad))
    i_apx = 2*np.square(l2_m)*np.sqrt(con1)
    apx = -a_apx-b_apx-(c_apx/d_apx)-(1/e_apx)*f_apx*(g_apx+(h_apx/i_apx))
    
    fpx = (mp*apx)+((np.pi)*np.square(rp_m)*gasforce_theta_deg_cycles*0.1*1000000/4)
    
    fax = mrod*arx+fpx
    fay = (1/l2_m)*(((Izz*beta_dd-(fax*lg-fpx*(l2_m-lg)*np.sin(beta)))/np.cos(beta))+mrod*ary*(l2_m-lg))
    
    fx = fax*np.cos(thetacum_rad)+fay*np.sin(thetacum_rad)
    fy = fay*np.cos(thetacum_rad)-fax*np.sin(thetacum_rad)
    #forwardoffset = forwardoffset*np.pi/180
    #thetacum_rad = thetacum_rad
    rpm2radpersec = 6*sp.pi/180
    omega1_radpsec = engine_speed_rpm*rpm2radpersec
    time = thetacum_rad/omega1_radpsec
    torque = -fy*l1_m*1000
    forwardoffset_rad = forwardoffset*np.pi/180
    time_offset = forwardoffset_rad/omega1_radpsec
    offsetindex = (np.where(time == time_offset))
    torque = np.roll(torque,offsetindex[0][0])
    return torque, time, thetacum_deg
def torquegen_nomcyl():
    torque_cyl1, time_cyl1, theta_cyl1 = torquegen(float(60/1000), float(140/1000), 
                                                   float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080), 
                                                   'gaspressures.dat', 1, 1,0,1)
    torque_cyl5, time_cyl5, theta_cyl5 = torquegen(float(60/1000), float(140/1000), 
                                                   float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080), 
                                                   'gaspressures.dat', 1, 1,90,1)
    torque_cyl4, time_cyl4, theta_cyl4 = torquegen(float(60/1000), float(140/1000), 
                                                   float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080), 
                                                   'gaspressures.dat', 1, 1,180,1)
    torque_cyl8, time_cyl8, theta_cyl8 = torquegen(float(60/1000), float(140/1000), 
                                                   float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080), 
                                                   'gaspressures.dat', 1, 1,270,1)
    torque_cyl6, time_cyl6, theta_cyl6 = torquegen(float(60/1000), float(140/1000), 
                                                   float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080), 
                                                   'gaspressures.dat', 1, 1,360,1)
    torque_cyl3, time_cyl3, theta_cyl3 = torquegen(float(60/1000), float(140/1000), 
                                                   float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080), 
                                                   'gaspressures.dat', 1, 1,360,1)   
fy_1, time_1, theta_1 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,0,1)
fy_5, time_5, theta_5 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,90,1)
fy_4, time_4, theta_4 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,180,1)
fy_8, time_8, theta_8 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,270,1)
fy_6, time_6, theta_6 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,360,1)
fy_3, time_3, theta_3 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,450,1)
fy_7, time_7, theta_7 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,540,1)
fy_2, time_2, theta_2 = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,630,1)
fy_1 = fy_1 + np.random.normal(0,((((max(fy_1))-(min(fy_1)))/2))*0.01,len(fy_1))
fy_5 = fy_5 + np.random.normal(0,((((max(fy_5))-(min(fy_5)))/2))*0.01,len(fy_5))
fy_4 = fy_4 + np.random.normal(0,((((max(fy_4))-(min(fy_4)))/2))*0.01,len(fy_4))
fy_8 = fy_8 + np.random.normal(0,((((max(fy_8))-(min(fy_8)))/2))*0.01,len(fy_8))
fy_6 = fy_6 + np.random.normal(0,((((max(fy_6))-(min(fy_6)))/2))*0.01,len(fy_6))
fy_3 = fy_3 + np.random.normal(0,((((max(fy_3))-(min(fy_3)))/2))*0.01,len(fy_3))
fy_7 = fy_7 + np.random.normal(0,((((max(fy_7))-(min(fy_7)))/2))*0.01,len(fy_7))
fy_2 = fy_2 + np.random.normal(0,((((max(fy_2))-(min(fy_2)))/2))*0.01,len(fy_2))
plt.rcParams['figure.figsize'] = [10, 10]
plt.subplot(2,1,1)

plt.plot(time_1,fy_1,label='cyl_1')
plt.plot(time_1,fy_5,"--",label='cyl_5')
plt.plot(time_1,fy_4,"--",label='cyl_4')
plt.plot(time_1,fy_8,"--",label='cyl_8')
plt.plot(time_1,fy_6,"--",label='cyl_6')
plt.plot(time_1,fy_3,"--",label='cyl_3')
plt.plot(time_1,fy_7,"--",label='cyl_7')
plt.plot(time_1,fy_2,"--",label='cyl_2')
plt.legend()
plt.xlabel("time [sec]")
plt.ylabel("torque [N-mm]")
plt.xlim(0,0.8)
plt.title("V8 Engine | Firing Order: 1-5-4-8-6-3-7-2")

plt.subplot(2,1,2)
plt.plot(theta_1,fy_1,label='cyl_1')
plt.plot(theta_1,fy_5,"--",label='cyl_5')
plt.plot(theta_1,fy_4,"--",label='cyl_4')
plt.plot(theta_1,fy_8,"--",label='cyl_8')
plt.plot(theta_1,fy_6,"--",label='cyl_6')
plt.plot(theta_1,fy_3,"--",label='cyl_3')
plt.plot(theta_1,fy_7,"--",label='cyl_7')
plt.plot(theta_1,fy_2,"--",label='cyl_2')
plt.legend(loc = 'upper left')
plt.xlim(0,720)
plt.xlabel("crank angle (theta) [degree]")
plt.ylabel("torque [N-mm]")


plt.show()
torque_cylnom = fy_1+fy_5+fy_4+fy_8+fy_6+fy_3+fy_7+fy_2
fy_1_mf, time_1_mf, theta_1_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,0,0.1)
fy_5_mf, time_5_mf, theta_5_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,90,0.1)
fy_4_mf, time_4_mf, theta_4_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,180,0.1)
fy_8_mf, time_8_mf, theta_8_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,270,0.1)
fy_6_mf, time_6_mf, theta_6_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,360,0.1)
fy_3_mf, time_3_mf, theta_3_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,450,0.1)
fy_7_mf, time_7_mf, theta_7_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,540,0.1)
fy_2_mf, time_2_mf, theta_2_mf = torquegen(float(60/1000), float(140/1000),
                 float(0.6833038489), float( 1e-9), float(40/1000), float(0.3282204182), float(2000), float(0.080),
                 'gaspressures.dat', 20, 0.1,630,0.1)
fy_1_mf = fy_1_mf + np.random.normal(0,((((max(fy_1_mf))-(min(fy_1_mf)))/2))*0.01,len(fy_1_mf))
fy_5_mf = fy_5_mf + np.random.normal(0,((((max(fy_5_mf))-(min(fy_5_mf)))/2))*0.01,len(fy_5_mf))
fy_4_mf = fy_4_mf + np.random.normal(0,((((max(fy_4_mf))-(min(fy_4_mf)))/2))*0.01,len(fy_4_mf))
fy_8_mf = fy_8_mf + np.random.normal(0,((((max(fy_8_mf))-(min(fy_8_mf)))/2))*0.01,len(fy_8_mf))
fy_6_mf = fy_6_mf + np.random.normal(0,((((max(fy_6_mf))-(min(fy_6_mf)))/2))*0.01,len(fy_6_mf))
fy_3_mf = fy_3_mf + np.random.normal(0,((((max(fy_3_mf))-(min(fy_3_mf)))/2))*0.01,len(fy_3_mf))
fy_7_mf = fy_7_mf + np.random.normal(0,((((max(fy_7_mf))-(min(fy_7_mf)))/2))*0.01,len(fy_7_mf))
fy_2_mf = fy_2_mf + np.random.normal(0,((((max(fy_2_mf))-(min(fy_2_mf)))/2))*0.01,len(fy_2_mf))
torque_cyl1mf = fy_1_mf+fy_5   +fy_4   +fy_8   +fy_6   +fy_3   +fy_7   +fy_2
torque_cyl5mf = fy_1   +fy_5_mf+fy_4   +fy_8   +fy_6   +fy_3   +fy_7   +fy_2
torque_cyl4mf = fy_1   +fy_5   +fy_4_mf+fy_8   +fy_6   +fy_3   +fy_7   +fy_2
torque_cyl8mf = fy_1   +fy_5   +fy_4   +fy_8_mf+fy_6   +fy_3   +fy_7   +fy_2
torque_cyl6mf = fy_1   +fy_5   +fy_4   +fy_8   +fy_6_mf+fy_3   +fy_7   +fy_2
torque_cyl3mf = fy_1   +fy_5   +fy_4   +fy_8   +fy_6   +fy_3_mf+fy_7   +fy_2
torque_cyl7mf = fy_1   +fy_5   +fy_4   +fy_8   +fy_6   +fy_3   +fy_7_mf+fy_2
torque_cyl2mf = fy_1   +fy_5   +fy_4   +fy_8   +fy_6   +fy_3   +fy_7   +fy_2_mf

torque_cyl1234mf = fy_1_mf+fy_5   +fy_4_mf+fy_8   +fy_6   +fy_3_mf+fy_7   +fy_2_mf
torque_cyl5678mf = fy_1   +fy_5_mf+fy_4   +fy_8_mf+fy_6_mf+fy_3   +fy_7_mf+fy_2
plt.rcParams['figure.figsize'] = [25, 10]

plt.plot(time_1,torque_cylnom, label ="Nominal")
plt.plot(time_1,torque_cyl1mf, label ="cyl1mf")
plt.plot(time_1,torque_cyl5mf, label ="cyl5mf")
plt.plot(time_1,torque_cyl4mf, label ="cyl4mf")
plt.plot(time_1,torque_cyl8mf, label ="cyl8mf")
plt.plot(time_1,torque_cyl6mf, label ="cyl6mf")
plt.plot(time_1,torque_cyl3mf, label ="cyl3mf")
plt.plot(time_1,torque_cyl7mf, label ="cyl7mf")
plt.plot(time_1,torque_cyl2mf, label ="cyl2mf")
plt.plot(time_1,torque_cyl1234mf, label ="cyl1234mf")
plt.plot(time_1,torque_cyl5678mf, label ="cyl5678mf")
plt.xlim(0,0.1)
plt.xlabel("time [sec]")
plt.ylabel("torque [N-mm]")
plt.title("Nominal Torque")
plt.legend()
plt.show()
x = np.block([[vibr_cylnom],
              [vibr_cyl1mf],
              [vibr_cyl5mf],
              [vibr_cyl4mf],
              [vibr_cyl8mf],
              [vibr_cyl6mf],
              [vibr_cyl3mf],
              [vibr_cyl7mf],
              [vibr_cyl2mf],
              [vibr_cyl1234mf],
              [vibr_cyl5678mf]])
x
y = np.array([[0,0,0,0,0,0,0,0],
            [1,0,0,0,0,0,0,0],
            [0,1,0,0,0,0,0,0],
            [0,0,1,0,0,0,0,0],
            [0,0,0,1,0,0,0,0],
            [0,0,0,0,1,0,0,0],
            [0,0,0,0,0,1,0,0],
            [0,0,0,0,0,0,1,0],
            [0,0,0,0,0,0,0,1],
            [1,0,1,0,0,1,0,1],
            [0,1,0,1,1,0,1,0]], dtype = float)
y
dataframe_np = np.block([x,y])
dataframe_np
v8_df = pd.DataFrame(dataframe_np)
v8_df
x,y
x_train, x_test, y_train, y_test = train_test_split(x,y, test_size = 0.1, random_state=0)
#x_train = np.reshape(x_train, (-1,1))
#x_test = np.reshape(x_test, (-1,1))
#y_train = np.reshape(y_train, (-1,1))
#y_test = np.reshape(y_test, (-1,1))
np.shape(x_train), np.shape(y_train)
np.shape(x_test), np.shape(y_test )
from sklearn.svm import SVC
classifier = RandomForestClassifier(n_jobs=2, random_state = 0 )
classifier.fit(x_train, y_train)
y_pred = classifier.predict(x_test)
y_pred,y_test<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup