Commit fa2497ae authored by Chetan Sharma's avatar Chetan Sharma
Browse files

Saving progress

parent 61bc7546
%% Cell type:markdown id: tags:
# Data Verification
%% Cell type:code id: tags:
``` python
import numpy as np
import time
import shelve
import logging
import queue
from collect_data import DataPoint, RawDataPacket
# from collect_data import DataPoint, RawDataPacket
from matplotlib import pyplot as plt
from scipy.signal import medfilt
TEST_NAMES = ['6']
results = list()
with shelve.open('results') as db:
for test_name in TEST_NAMES:
results += db[test_name]
data_points = list(map(lambda x: x.decompose_packet(), results))
print(data_points[0].W)
x = list(map(lambda x: x.f_r, data_points))
y = list(map(lambda x: x.T, data_points))
plt.figure()
plt.plot(x, y)
```
%% Output
0.003175
[<matplotlib.lines.Line2D at 0x7f9bf7983e50>]
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
~/anaconda3/envs/ampd/lib/python3.8/shelve.py in __getitem__(self, key)
110 try:
--> 111 value = self.cache[key]
112 except KeyError:
KeyError: '6'
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
<ipython-input-2-aae045e25021> in <module>
13 with shelve.open('results') as db:
14 for test_name in TEST_NAMES:
---> 15 results += db[test_name]
16
17 data_points = list(map(lambda x: x.decompose_packet(), results))
~/anaconda3/envs/ampd/lib/python3.8/shelve.py in __getitem__(self, key)
112 except KeyError:
113 f = BytesIO(self.dict[key.encode(self.keyencoding)])
--> 114 value = Unpickler(f).load()
115 if self.writeback:
116 self.cache[key] = value
AttributeError: Can't get attribute 'RawDataPacket' on <module '__main__'>
%% Cell type:code id: tags:
``` python
%matplotlib notebook
new_data = list(map(lambda x: RawDataPacket(x.D, x.W, x.R, x.N, x.f_r, x.w, x.measurements), results))
data_points = list(map(lambda x: x.decompose_packet_sorta(), new_data))
plt.figure()
last_time = 0
for tTs, tFys, Ts, Fys in data_points:
times = [t + last_time for t in tTs]
stuff = Ts
filter_size = 3
fc = int((filter_size + 1) / 2)
plt.plot(times[fc: -fc], medfilt(stuff, filter_size)[fc: -fc])
last_time = times[-1]
```
%% Output
%% Cell type:markdown id: tags:
## Plan
We have three datasets. Two variables are tweaked across dataset. We want to perform some kind of regression to figure out our cutting force coefficients.
Regression is linear if we only vary feed, but becomes nonlinear once we incorporate WOC.
First, maybe just try using a linear regressor for both?
%% Cell type:code id: tags:
``` python
from scipy.optimize import least_squares
def single_residual(d, K_tc, K_te, p):
return T_exp_lin(d.D, d.N, d.R, d.W, d.f_r, d.w, K_tc, K_te, p) - d.T
def gen_train_pairs():
X = [[d.D, d.N, d.R, d.W, d.f_r, d.w] for d in data_points]
y = [d.T for d in data_points]
return (X, y)
residual_array = np.vectorize(single_residual)
def objective(x):
K_tc = x[0]
K_te = x[1]
p = x[2]
return residual_array(data_points, K_tc, K_te, p)
class NonLinRegressor:
def fit(self, X, Y):
self.coef = least_squares(lambda coef: [self.evaluate(x, coef) - y for x, y in zip(X, Y)], self.coef).x
return self
def predict(self, X):
return [self.evaluate(x, self.coef) for x in X]
def get_params(self, deep=True):
return {}
def set_params(self, **params):
...
def score(self, X, Y):
u = sum([(y_true - y_pred) ** 2 for y_true, y_pred in zip(Y, self.predict(X))])
y_true_mean = np.mean(Y)
v = sum([(y_true - y_true_mean) ** 2 for y_true in Y])
return (1 - u/v)
class TRegressor(NonLinRegressor):
def __init__(self):
self.coef = [0,0,1]
def evaluate(self, x, coef):
return T_exp_lin(*x, coef[0], coef[1], coef[2])
class FRegressor(NonLinRegressor):
def __init__(self):
self.coef = [0,0,0,0, 0, 0]
def evaluate(self, x, coef):
return F_exp_lin(*x, *coef)[0][1]
estimator = TRegressor()
regressor = RANSACRegressor(base_estimator = estimator, min_samples = 5)
X, y = gen_train_pairs()
# y = [d.Fy * 9.81 ** 2 for d in data_points]
regressor.fit(X, y)
# K_tc, K_te, K_rc, K_re, p, q = regressor.estimator_.coef
K_tc, K_te, p = regressor.estimator_.coef
```
%% Cell type:code id: tags:
``` python
# print(K_tc, K_te, K_rc, K_re, p, q)
y_est = list(map(lambda d : T_exp_lin(d.D, d.N, d.R, d.W, d.f_r, d.w, K_tc, K_te, p), data_points))
# y_est = list(map(lambda d : F_exp_lin(d.D, d.N, d.R, d.W, d.f_r, d.w, K_tc, K_te, K_rc, K_re, p, q), data_points))
# y_est = np.array(y_est)
# y_est = np.squeeze(y_est)[:, 1]
# print(y_est)
plt.figure()
plt.plot(x, y)
plt.plot(x, y_est)
```
%% Output
[<matplotlib.lines.Line2D at 0x7fc984255c10>]
%% Cell type:code id: tags:
``` python
from sklearn.linear_model import RANSACRegressor, LinearRegression
from models import T_x_vector, T_lin, T_exp_lin, calc_h_a, calc_f_t, F_exp_lin
def dataPoint_to_X(d):
return T_x_vector(d.D, d.N, d.R, d.W, d.f_r, d.w)
def single_residual(d, K_tc, K_te, p):
return T_exp_lin(d.D, d.N, d.R, d.W, d.f_r, d.w, K_tc, K_te, p) - d.T
residual_array = np.vectorize(single_residual)
def objective(x):
K_tc = x[0]
K_te = x[1]
p = x[2]
return residual_array(data_points, K_tc, K_te, p)
X = list(map(dataPoint_to_X, data_points))
estimator = LinearRegression(fit_intercept=False)
regressor = RANSACRegressor(base_estimator = estimator)
regressor.fit(X, y)
K_tc, K_te = regressor.estimator_.coef_
print("K_tc = ", K_tc, ", K_te = ", K_te)
x_h_a = np.linspace(0.1e-3, 4e-3, 100)
y_h_a = [calc_h_a(6.35e-3 / 2 , x, calc_f_t(3, 0.0016, 500)) for x in x_h_a]
```
%% Output
K_tc = 858494934.9591976 , K_te = -696.3150933941946
%% Cell type:markdown id: tags:
Now we check if we actually fit anything useful
%% Cell type:code id: tags:
``` python
y_est = list(map(lambda d : T_lin(d.D, d.N, d.R, d.W, d.f_r, d.w, K_tc, K_te), data_points))
plt.figure()
plt.plot(x, y)
plt.plot(x, y_est)
plt.figure()
plt.plot(x_h_a, y_h_a)
plt.plot([0.00079375, 0.0015875, 0.003175], [0,0,0], 'd')
```
%% Output
[<matplotlib.lines.Line2D at 0x7f22f8fe0640>]
......
from ml import LinearModel
from fake_cut import Fake_Cut
from objects import EndMill, Conditions
PARAMS = [858494934.9591976, -696.3150933941946, 858494934.9591976, -696.3150933941946]
# ERROR = [0.01, 0.01]
# NOISE = [0.1, 0.1]
ERROR = [0, 0]
NOISE = [0, 0]
CONDITIONS = Conditions(1e-3, 3.175e-3, 5e-3, 5000, EndMill(3, 3.175e-3, 3.175e-3, 10e-3, 10e-3))
model = LinearModel()
cut = Fake_Cut(PARAMS, ERROR, NOISE)
data = cut.cut(CONDITIONS)
model.ingest_datum(data)
print("{0:.20f}".format(model.params[1]))
\ No newline at end of file
......@@ -62,10 +62,11 @@ class Cut:
self._logger.info("Layer prepared for clearing")
def cut(self, W, f_r, w):
def cut(self, conditions):
"""
Performs a stroke of facing. Returns a data blob.
"""
_, W, f_r, w, _ = conditions.unpack()
X_START = x_cut - self.endmill.r_c + W
if X_START > self.X_END:
raise MachineCrash("Cutting too far in X direction: X = " + str(X_START))
......
from models import T_lin, F_lin
from objects import Conditions, Data
import numpy as np
class Fake_Cut:
def __init__(self, params, error, noise):
"""
Args:
params: list of format [K_tc, K_te, K_rc, K_re]
error: list of standard deviations of format [o_T, o_Fy]. Simulates the sensor being "wrong" sometimes.
error: list of standard deviations of format [o_T, o_Fy]. Simulates the sensor being noisy.
"""
self.params = params
self.error = error
self.noise = noise
def cut(self, conditions : Conditions):
# use prediction as output
T = T_lin(conditions, *self.params[:2])
_, Fy = F_lin(conditions, *self.params)
# add sensor error
T_error = np.random.normal(T, self.error[0])
Fy_error = np.random.normal(Fy, self.error[1])
# add sensor noise
T_noisy = np.random.normal(T_error, self.noise[0], (100))
Fy_noisy = np.random.normal(Fy_error, self.noise[1], (100))
# generate fake times
t = np.linspace(0, 1, 100)
# return fake reading
return Data(*conditions.unpack(), np.array([t, T_noisy]), np.stack([t, Fy_noisy]))
\ No newline at end of file
......@@ -5,7 +5,7 @@ from sklearn import linear_model
from matplotlib import pyplot as plt
from models import T_lin, F_lin, T_x_vector, Fy_x_vector
from objects import Data, Reading, EndMill, Prediction
from objects import Data, Conditions, EndMill, Prediction
class Model(abc.ABC):
"""
......@@ -40,19 +40,21 @@ class LinearModel(Model):
self.training_T_y = list()
self.training_Fy_x = list()
self.training_Fy_y = list()
self.regressor_T = linear_model.Lasso(fit_intercept = False, warm_start=True)
self.regressor_Fy = linear_model.Lasso(fit_intercept = False, warm_start=True)
self.regressor_T = linear_model.LinearRegression(fit_intercept = False)
self.regressor_Fy = linear_model.LinearRegression(fit_intercept = False)
self.params = np.array([0,0,0,0])
def ingest_datum(self, datum):
# decompose
D, W, f_r, w, endmill, Ts, Fys = datum.D, datum.W, datum.f_r, datum.w, datum.endmill, datum.Ts, datum.Fys
R, N = endmill.R, endmill, N
T, Fy = np.median(Ts[:, 1]), np.median(Fys[:, 1])
_, _, _, _, _, Ts, Fys = datum.unpack()
T, Fy = np.median(Ts[1, :]), np.median(Fys[1, :])
# get linear coefficients
T_x = T_x_vector(D, N, R, W, f_r, w)
Fy_x = Fy_x_vector(D, N, R, W, f_r, w)
T_x = T_x_vector(datum.conditions())
Fy_x = Fy_x_vector(datum.conditions())
print("T", T)
print("coef times real K_tc", T_x[1] * 858494934.9591976)
# add to training set
self.training_T_x.append(T_x)
......@@ -68,24 +70,25 @@ class LinearModel(Model):
# calculate best fit from data
self.regressor_T.fit(self.training_T_x, self.training_T_y)
K_tc, K_te = self.regressor.coef_
K_tc, K_te = self.regressor_T.coef_
print(K_tc, K_te)
self.params[0], self.params[1] = K_tc, K_te
# transform Fy into a smaller linear problem and fit
intercepts = training_Fy_x @ np.array([K_tc, K_te, 0, 0])[np.newaxis].T
training_Fy_y_no_intercepts = training_Fy_y - intercepts
self.regressor_Fy.fit(training_Fy_x[:, 2:], training_Fy_y_no_intercepts)
K_rc, K_re = self.regressor_Fy.coef_
self.params[0], self.params[1] = K_rc, K_re
def predict_one(self, conditions):
# decompose
D, W, f_r, w, endmill = conditions.D, conditions.W, conditions.f_r, conditions.w, conditions.endmill
R, N = endmill.R, endmill, N
# evaluate
T = T_lin(D, N, R, W, f_r, w, self.params[0], self.params[1])
F = F_lin(D, N, R, W, f_r, w, *self.params)
T = T_lin(conditions, *self.params[:2])
F = F_lin(conditions, *self.params)
# repack and return
return Prediction(D, W, f_r, w, endmill, T, F)
return Prediction(conditions, T, F)
import numpy as np
def calc_h_a(R, W, f_t):
from objects import Conditions, Data, Prediction, MachineChar
def calc_f_t(conditions : Conditions):
D, W, f_r, w, endmill = conditions.unpack()
N = endmill.N
return (2 * np.pi * f_r) / (N * w)
def calc_h_a(conditions : Conditions):
D, W, f_r, w, endmill = conditions.unpack()
R = endmill.r_c
f_t = calc_f_t(conditions)
phi_st = np.pi - np.arccos(1 - W / R)
phi_ex = np.pi
return -f_t * (np.cos(phi_ex) - np.cos(phi_st)) / (phi_ex - phi_st)
def calc_f_t(N, f_r, w):
return (2 * np.pi * f_r) / (N * w)
# exponential models, account for strange decrease in cutting forces as chip thickness increases.
def T_exp_lin(D, N, R, W, f_r, w, K_TC, K_te, p):
f_t = calc_f_t(N, f_r, w)
h_a = calc_h_a(R, W, f_t)
def T_exp_lin(conditions : Conditions, K_TC, K_te, p):
D, W, f_r, w, endmill = conditions.unpack()
N, R, _, _, _ = endmill.unpack()
f_t = calc_f_t(conditions)
h_a = calc_h_a(conditions)
K_tc = K_TC * h_a ** -p
return ((D * N) * (K_te * R * np.arccos(1 - W/R) + K_tc * W * f_t)) / (2 * np.pi)
def F_exp_lin(D, N, R, W, f_r, w, K_TC, K_te, K_RC, K_re, p, q):
f_t = calc_f_t(N, f_r, w)
h_a = calc_h_a(R, W, f_t)
def F_exp_lin(conditions : Conditions, K_TC, K_te, K_RC, K_re, p, q):
D, W, f_r, w, endmill = conditions.unpack()
N, R, _, _, _ = endmill.unpack()
f_t = calc_f_t(conditions)
h_a = calc_h_a(conditions)
K_tc = K_TC * h_a ** -p
K_rc = K_RC * h_a ** -q
Fx = - ((D*N*(K_tc*W**2*f_t + K_rc*W**(3/2)*f_t*(2*R - W)**(1/2)))/4 - (D*N*R*(2*K_tc*W*f_t - 2*K_re*W + 2*K_te*W**(1/2)*(2*R - W)**(1/2) + K_rc*W**(1/2)*f_t*(2*R - W)**(1/2)))/4)/(R**2*np.pi) - (D*K_rc*N*f_t*np.arccos((R - W)/R))/(4*np.pi)
......@@ -26,36 +38,42 @@ def F_exp_lin(D, N, R, W, f_r, w, K_TC, K_te, K_RC, K_re, p, q):
# linear models, easier to fit and more simple
def T_lin(D, N, R, W, f_r, w, K_tc, K_te):
f_t = calc_f_t(N, f_r, w)
def T_lin(conditions, K_tc, K_te):
D, W, f_r, w, endmill = conditions.unpack()
N, R, _, _, _ = endmill.unpack()
f_t = calc_f_t(conditions)
return ((D * N) * (K_te * R * np.arccos(1 - W/R) + K_tc * W * f_t)) / (2 * np.pi)
def F_lin(D, N, R, W, f_r, w, K_tc, K_te, K_rc, K_re):
f_t = calc_f_t(N, f_r, w)
def F_lin(conditions, K_tc, K_te, K_rc, K_re):
D, W, f_r, w, endmill = conditions.unpack()
N, R, _, _, _ = endmill.unpack()
f_t = calc_f_t(conditions)
Fx = - ((D*N*(K_tc*W**2*f_t + K_rc*W**(3/2)*f_t*(2*R - W)**(1/2)))/4 - (D*N*R*(2*K_tc*W*f_t - 2*K_re*W + 2*K_te*W**(1/2)*(2*R - W)**(1/2) + K_rc*W**(1/2)*f_t*(2*R - W)**(1/2)))/4)/(R**2*np.pi) - (D*K_rc*N*f_t*np.arccos((R - W)/R))/(4*np.pi)
Fy = (D*K_tc*N*f_t*np.arccos((R - W)/R))/(4*np.pi) - ((D*N*(K_rc*W**2*f_t - K_tc*W**(3/2)*f_t*(2*R - W)**(1/2)))/4 - (D*N*R*(2*K_te*W + 2*K_rc*W*f_t + 2*K_re*W**(1/2)*(2*R - W)**(1/2) - K_tc*W**(1/2)*f_t*(2*R - W)**(1/2)))/4)/(R**2*np.pi)
return np.array([Fx, Fy])
# coefficient calculators for linear regression
def T_x_vector(D, N, R, W, f_r, w):
def T_x_vector(conditions : Conditions):
"""
Outputs vector that corresponds to coefficients in order:
K_tc, K_te
"""
f_t = calc_f_t(N, f_r, w)
K_tc_C = (D * N * W * f_t) / (2 * np.pi)
K_te_C = (D * N * R * np.arccos(1 - (W / R))) / (2 * np.pi)
return [K_tc_C, K_te_C]
D, W, f_r, w, endmill = conditions.unpack()
N, R, _, _, _ = endmill.unpack()
f_t = calc_f_t(conditions)
return[(D * N * W * f_t) / (2 * np.pi), (D * N * R * np.arccos(1 - (W / R))) / (2 * np.pi)]
def Fy_x_vector(D, N, R, W, f_r, w):
f_t = calc_f_t(N, f_r, w)
def Fy_x_vector(conditions : Conditions):
D, W, f_r, w, endmill = conditions.unpack()
N, R, _, _, _ = endmill.unpack()
f_t = calc_f_t(conditions)
return [(D*N*f_t*(W**(3/2)*(2*R - W)**(1/2) + R**2*np.arccos((R - W)/R) - R*W**(1/2)*(2*R - W)**(1/2)))/(4*R**2*np.pi), (D*N*W)/(2*R*np.pi), (D*N*W*f_t*(2*R - W))/(4*R**2*np.pi), (D*N*W**(1/2)*(2*R - W)**(1/2))/(2*R*np.pi)]
def deflection_load(D_a, F, endmill, E_e = 650e9):
def deflection_load(D_a, prediction : Prediction, E_e = 650e9):
"""
Calculates ratio of tool deflection to allowed deflection.
Uses FEA model described in doi:10.1016/j.ijmachtools.2005.09.009.
......@@ -67,11 +85,10 @@ def deflection_load(D_a, F, endmill, E_e = 650e9):
Returns:
Ratio of tool deflection to allowed deflection
"""
r_c = endmill.r_c
r_s = endmill.r_s
l_c = endmill.l_c
l_s = endmill.l_s
N = endmill.N
D, W, f_r, w, endmill, T, F = prediction.unpack()
F = np.norm(F)
N, r_c, r_s, l_c, l_s = endmill.unpack()
# prepare variables, this must be done since the FEA model has arbitrary constants defined for certain units
d_1 = r_c * 2 * 1e3 # convert to diameter in mm
d_2 = r_s * 2 * 1e3
......@@ -98,7 +115,7 @@ def deflection_load(D_a, F, endmill, E_e = 650e9):
D_m = C * (F / E) * ((l_1 ** 3 / d_1 ** 4) + ((l_2 ** 3 - l_1 ** 3) / d_2 ** 4)) ** G * 1e-3
return D_m / D_a
def failure_prob_milling(F, T, endmill, D, roughing = False, m = 4, o = 1500e6, a_c = 0.8):
def failure_prob_milling(prediction : Prediction, roughing = False, m = 4, o = 1500e6, a_c = 0.8):
"""
Calculates failure probability according to Weibull distribution. Method adapted from https://doi.org/10.1016/S0007-8506(07)62072-1
Args:
......@@ -113,11 +130,10 @@ def failure_prob_milling(F, T, endmill, D, roughing = False, m = 4, o = 1500e6,
Returns:
Probability of failure
"""
# establish base variables
r_c = endmill.r_c
r_s = endmill.r_s
l_c = endmill.l_c
l_s = endmill.l_s
D, W, f_r, w, endmill, T, F = prediction.unpack()
F = np.norm(F)
N, r_c, r_s, l_c, l_s = endmill.unpack() # establish base variables
r_ceq = r_c * a_c # equivalent cutter radius
I_ceq = np.pi / 4 * r_ceq ** 4 # equivalent cutter 2nd inertia
I_s = np.pi / 4 * r_s ** 4 # shank 2nd inertia
......@@ -137,7 +153,7 @@ def failure_prob_milling(F, T, endmill, D, roughing = False, m = 4, o = 1500e6,
return failure_prob
def motor_torque_load(T_m, w, K_T, R_w, V_max, I_max):
def motor_torque_load(prediction : Prediction, machinechar : MachineChar):
"""
Calculates ratio of current motor torque to max motor torque.
Args:
......@@ -150,12 +166,16 @@ def motor_torque_load(T_m, w, K_T, R_w, V_max, I_max):
Returns:
Ratio of current motor torque to max achievable torque
"""
_, _, _, w, _, T, _ = prediction.unpack()
r_e, K_T, R_w, V_max, I_max, T_nom, _ = machinechar.unpack()
w_m = w * r_e
T_m = (T + T_nom) / r_e
# store K_V for convenience
K_V = 1 / K_T
# max torque is either determined by max current that can be supplied or max current achievable given winding resistance and whatnot
return T_m / min(K_T * I_max, (V_max - K_V * w) / R_w)
return T_m / min(K_T * I_max, (V_max - K_V * w_m) / R_w)
def motor_speed_load(T_m, w, K_T, R_w, V_max):
def motor_speed_load(prediction : Prediction, machinechar : MachineChar):
"""
Calculates ratio of
Args:
......@@ -167,11 +187,16 @@ def motor_speed_load(T_m, w, K_T, R_w, V_max):
Returns:
Ratio of current motor speed to max achievable speed
"""
_, _, _, w, _, T, _ = prediction.unpack()
r_e, K_T, R_w, V_max, I_max, T_nom, _ = machinechar.unpack()
w_m = w * r_e
T_m = (T + T_nom) / r_e
# store K_V for convenience
K_V = 1 / K_T
# max speed is affected by winding resistive voltage drop along with backemf
return w / (K_V * (V_max - T_m * R_w / K_T))
return w_m / (K_V * (V_max - T_m * R_w / K_T))
def optimality(D, W, f_r):
def optimality(conditions : Conditions):
D, W, f_r, _, _ = conditions.unpack()
# returns MMR in units of m^3 / s
return D * W * f_r
......@@ -2,20 +2,47 @@ import numpy as np
class EndMill:
def __init__(self, N, r_c, r_s, l_c, l_s):
self.R = R
self.N = N
self.r_c = r_c
self.r_s = r_s
self.l_c = l_c
self.l_s = l_s
def unpack(self):
return [self.N, self.r_c, self.r_s, self.l_c, self.l_s]
def __str__(self):
return str(self.unpack())
class MachineChar:
def __init__(self, r_e, K_T, R_w, V_max, I_max, T_nom, f_r_max):
self.r_e = r_e
self.K_T = K_T
self.R_w = R_w
self.V_max = V_max
self.I_max = I_max
self.T_nom = T_nom
self.f_r_max = f_r_max
def unpack(self):
return [self.r_e, self.K_T, self.R_w, self.V_max, self.I_max, self.T_nom]
def __str__(self):
return str(self.unpack())
class Conditions:
def __init__(self, D : float, W : float, f_r : float, w : float, endmill : EndMi