WeldingSpotPerformance/OptTimeCalculator.py

from numpy import sqrt, arcsin, arccos, cos, sin

from scipy.optimize import minimize
import numpy

from Params import Params



class OptTimeCalculator():
    
    def __init__(self, params : Params):
        self.Ts = {}
        self.params = params
        self.Tmax = 1
        
    def tGrow(self, F):
        return arcsin(F/(self.params.Ftg)) * sqrt(self.params.m/self.params.k)

    def T(self, s, s2 = 0, debug = False):
        if hasattr(s, "__len__"):
            s = s[0]
        v0q = min(sqrt(2 * self.params.a * s), self.params.v1limit)
        v0 = min(v0q, sqrt(1/(self.params.k*self.params.m))* self.params.Ftg)
        t1 = v0 / self.params.a
        
        #test second
        t21 = sqrt(s2/self.params.a2)
        t21 = min(self.params.v2limit/self.params.a2, t21)
        t22 = max(0, (s2 - (self.params.a2 * t21 * t21)) / self.params.v2limit)
        T2 = t22 + 2 * t21
        
        t2t = max(0, (s - (self.params.a * t1 * t1 /2)) / v0q)
        #if debug:
            #print("!", t1, t2t, "|", t21, t22, "|", t1+t2t, 2*t21+t22, t2t + t1 < T2)
        if t2t + t1 < T2:
            if debug:
                print("bad time")
            return 2 * self.Tmax
        #print("t2", t2t) 
        s3 = s+((self.params.Fd - self.params.fl)/2 + self.params.fl)/self.params.k
        t3 = sqrt(s3 / self.params.a)
        v3max = min(t3 * self.params.a, self.params.v1limit)
        t3 = v3max / self.params.a
        t3q = max(0, (s3 - t3*t3*self.params.a)/v3max)
        v0 = v0 * self.params.k
        l = sqrt(self.params.awork**2 * self.params.m**2 + self.params.m/self.params.k * v0**2)
        #print("v0", v0, v0/k)
        Fmeet = (self.params.Ftg**2 - self.params.m/self.params.k*v0**2)/(2 * self.params.awork * self.params.m)
        print(Fmeet)
        #if debug:
            #print(Fmeet)
        Fq = self.params.Fq
        if Fmeet > Fq:
            if debug:
                print("cant reach", Fmeet)
            return self.Tmax
        #print(Ftg - Fmeet)
        if Fmeet - self.params.awork*self.params.m > l:
            if debug:
                print("wut")
            return 3*self.Tmax
        beta = self.params.awork * self.params.m
        t2 = sqrt(self.params.m/self.params.k) * (arcsin((Fmeet - beta)/(l)) + arccos(sqrt(self.params.m/self.params.k) * v0/l))
        if debug:
            qtq = sqrt(self.params.k/self.params.m)
            Fqq = -self.params.awork*self.params.m*cos(qtq*t2) + self.params.awork*self.params.m + v0 / qtq * sin(qtq*t2)
            print("angles", (Fmeet / (2*beta) -beta)/(l), (sqrt(self.params.m/self.params.k) * v0/l), Fqq, Fmeet)
        tend = self.tGrow(Fq) - self.tGrow(Fmeet)
        vp = 1/sqrt(self.params.k * self.params.m) * sqrt(self.params.Ftg**2 - self.params.Fq**2)
        ap = Fq / self.params.m
        tprop = 2*vp / ap
        #print(s, t1, t2, Fmeet, tend)
        T = t1 + t2t + t2 + 2*t3 + t3q + tprop + tend
        
        self.Ts = {"t1":t1, "t2":t2t, "tm":t2, "tend":tend, "t4":tprop, "to1":t3, "to2": t3q, "T":T,
                   "tclose":t1+t2t, "tgrow":t2+tend+tprop, "topen": 2*t3+t3q}
        if debug:
            print("phases", t1, t2t, "|", t2, tend, tprop, "|", 2*t3+t3q, "|", T)
            print("tclose: ",t1+t2t, "tgrow:", t2+tend+tprop, "topen:", 2*t3+t3q, "all:", T)
        return T
    
    def calcSecond(self, Tstart, s2):
        T1 = self.Ts["tclose"]
        t1 = T1 / 2 - sqrt(T1**2 - 4 * s2 / self.params.a2) / 2
        t2 = sqrt(T1 * T1 - 4 * s2 / self.params.a2)
        
        T2 = self.Ts["topen"] - Tstart
        to1 = T2 / 2 - sqrt(T2**2 - 4 * s2 / self.params.a2) / 2
        to2 = sqrt(T2 * T2 - 4 * s2 / self.params.a2)
        
        Ts = {"t1":t1, "t2":t2, "to1":to1, "to2": to2}
        self.T2 = Ts
        return self.T2

    def calcOpt(self):
        s0 = [0.04]
        bnds = [(0.00001, None)]
        res = minimize(self.T, s0, method='Nelder-Mead', tol=1e-7, bounds=bnds)
        nonbound = res.x
        bound = res.x
        if nonbound < self.params.smin1:
            #print("IS SMIN")
            bound = max(nonbound, self.params.smin1)
        #T(q, debug = True)
        return nonbound, bound

    def showPlot(self):

        N = 50000
        ss = numpy.linspace(0.0001, 0.030, N)
        Ts = numpy.array([self.T(si) for si in ss])

        import matplotlib.pyplot as plt
        from matplotlib import use

        q, qq  = self.calcOpt()
        q = 1000*q
        qq = 1000 * qq
        use("GTK4Agg", force = True)
        plt.plot(ss*1000, Ts*1000, "red", label = "T(S)")
        plt.plot(numpy.ones(N) * qq, Ts*1000, label = "factOpt")
        plt.plot(numpy.ones(N) * q, Ts*1000, "black", label = "trueOpt")
        plt.title("T(s)")
        plt.xlabel("s, мм")
        plt.ylabel("T, мс")
        plt.legend()
        plt.savefig("model.png")
        plt.show()
        
    def getTimes(self, s):
        self.T(s, debug=False)
        return self.Ts
   

if __name__ == "__main__":
    opt = OptTimeCalculator() 
    nonbound, bound = opt.calcOpt()
    
    print("nonbound", nonbound*1000, "мм")
    print(opt.getTimes(nonbound))
    
    print("bound", bound*1000, "мм")
    print(opt.getTimes(bound))