优化 c 值以进行曲线拟合

Question

我有下面的代码，可以求解梁的无量纲固有频率，然后对其进行量纲化。我需要使用第一个频率模式通过优化 c1、c2、c3 和 c4 值来将曲线拟合到实验数据，但我不确定从哪里开始，并且到目前为止我尝试过的所有操作都产生了错误和不正确的值.

pdata = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
fdata = np.array([31786.90305, 33344.70228, 34007.47212, 34388.51794, 34640.69891, 34822.00471, 34959.70061, 35068.44216, 35156.86945, 35230.43423])

L = 0.0067  # Protruding length in meters
rw = 0.005  # Shank radius in meters
Rw = 0.01848 / 2  # Nut radius in meters

E = 200e9  # Pa (Young's Modulus)
I = 4.908738521234052e-10  # m^4 (Moment of Inertia)
rho = 7850  # kg/m^3 (Density)
A = 7.853981633974483e-05  # m^2 (Cross-sectional Area)

def create_Amat(eqns, vars):
    Amat = sp.zeros(len(eqns), len(vars))
    for i, eqn in enumerate(eqns):
        for j, var in enumerate(vars):
            Amat[i, j] = sp.diff(eqn, var)
    return Amat

def CBmodel_T(kt, kr, p):
    # Defining symbolic variables that will be used
    a1, a2, a3, a4, x, beta, omega = sp.symbols('a1 a2 a3 a4 x beta omega')

    # The Genereal Solution for free vibrations in a beam
    GS = a1 * sp.cos(beta * x) + a2 * sp.sin(beta * x) + a3 * sp.cosh(beta * x) + a4 * sp.sinh(beta * x)

    # Defining differentials of the General Solution equation
    # Left Hand Side (x=0)
    Y = GS
    Y1 = sp.diff(Y, x)
    Y2 = sp.diff(Y1, x)
    Y3 = sp.diff(Y2, x)

    # Right Hand Side (x=1)
    D = GS
    D1 = sp.diff(D, x)
    D2 = sp.diff(D1, x)
    D3 = sp.diff(D2, x)

    # Substituting values of x into the differentials
    #Left Hand Side (x=0)
    Y_0 = Y.subs(x, 0)
    Y_1 = Y1.subs(x, 0)
    Y_2 = Y2.subs(x, 0)
    Y_3 = Y3.subs(x, 0)

    # Right Hand Side (x=1)
    D_0 = D.subs(x, 1)
    D_1 = D1.subs(x, 1)
    D_2 = D2.subs(x, 1)
    D_3 = D3.subs(x, 1)

    # Defining Boundary Conditions for the system
    BC1 = Y_3 + kt * Y_0 - p * Y_1
    BC2 = Y_2 - kr * Y_1
    BC3 = D_2
    BC4 = D_3

    # Forming a matrix that contains only the coefficients
    # from the G.S. equation to find the non-trivial solution
    BCs = [BC1, BC2, BC3, BC4]
    var = [a1, a2, a3, a4]

    Amat = create_Amat(BCs, var)

    # Defining the symbolic variable for beta
    beta_sym = sp.sqrt(abs(omega))

    # Substitute symbolic values for beta1, kt, and kr into Amat
    Amat_numeric = Amat.subs({beta: beta_sym, kt: kt, kr: kr})

    # Find the determinant of Amat
    Amat_numeric = sp.simplify(sp.det(Amat_numeric))

    # Convert Amat_numeric to a function handle
    f_omega = lambdify(omega, Amat_numeric, 'numpy')

    # An empty array to store the roots
    n = []

    # Iterate through ranges of initial values for omega
    for j in range(0, 350, 1):
        bracket = [j, j + 1]
        try:
            result = root_scalar(f_omega, bracket=bracket)
            if result.converged and np.sign(f_omega(bracket[0])) != np.sign(f_omega(bracket[1])):
                n.append(result.root)
        except ValueError:
            # Ignore ValueError (f(a) and f(b) must have different signs)
            pass

    # Sorts and removes nonreal numbers
    n = np.array(n)
    n = n[n == np.real(n)]
    n = np.round(n[n == np.real(n)], decimals=10)
    n = np.unique(n)

    if kt == 0 or kr == 0:
        n = n[n > -1]
    
    else:
        n = n[n > 0]
    
    # Extract the first 5 unique values
    w1, w2, w3, w4, w5 = n[:5]

    # Return the non-dimensional frequencies
    return w1, w2, w3, w4, w5

def dimensionalise(w1, w2, w3, w4, w5, E, I, L, rho, A):

    # Equation used to dimensionalise each nondimensional natural frequency
    f1 = w1 * (1/(2*np.pi)) * np.sqrt((E * I) / (rho * A * L**4))
    f2 = w2 * (1/(2*np.pi)) * np.sqrt((E * I) / (rho * A * L**4))
    f3 = w3 * (1/(2*np.pi)) * np.sqrt((E * I) / (rho * A * L**4))
    f4 = w4 * (1/(2*np.pi)) * np.sqrt((E * I) / (rho * A * L**4))
    f5 = w5 * (1/(2*np.pi)) * np.sqrt((E * I) / (rho * A * L**4))

    # Return the dimensionalised frequencies
    return f1, f2, f3, f4, f5

def stiffnessfunction(p,L,rw,Rw,c1,c2,c3,c4):
    #Pre-allocate arrays for stiffness values
    kt = np.zeros_like(p)
    kr = np.zeros_like(p)

    #Normal/tangential stiffness
    kn = (c1*np.tanh(c2*p))/(1+(c3*np.exp(-c4*p)))

    #Constant parameters
    ct = ((np.pi/2)*(1-0.3))/(2-0.3)
    cr = (1/4)*((Rw**2+rw**2)/(L**2))

    #Translational and rotational stiffness parameters
    kt = ct*kn
    kr = cr*kn

    return kt, kr

Answer 1

这是我尝试使用的，但没有成功 将 numpy 导入为 np 从 scipy.optimize 导入 Differential_evolution

# Constants for kt and kr calculations
ct = ((np.pi/2) * (1 - 0.3)) / (2 - 0.3)
cr = (1/4) * ((0.00924**2 + 0.005**2) / (0.0067**2))

# Define the determinant function using numpy operations
def det_function(c, beta, p):
    c1, c2, c3, c4 = c
    kn = (c1 * np.tanh(c2 * p)) / (1 + c3 * np.exp(-c4 * p))
    kt = ct * kn
    kr = cr * kn
    return abs(-(2 * beta**6 * np.cos(beta) * np.cosh(beta) - 2 * beta**6) * kt * kr +
           (2 * beta**7 * np.sin(beta) * np.cosh(beta) - 2 * beta**7 * np.cos(beta) * np.sinh(beta)) * kt +
           (2 * beta**9 * np.sin(beta) * np.cosh(beta) + 2 * beta**9 * np.cos(beta) * np.sinh(beta)) * kr +
           2 * beta**10 * (np.cos(beta) * np.cosh(beta) - 1) + 
           2 * beta**8 * p * np.sin(beta) * np.sinh(beta))

# Function to calculate the sum of determinants
def sum_of_determinants(c):
    total_det = 0
    for beta_val in beta_values:
        for p_val in pdata:
            total_det += det_function(c, beta_val, p_val)
    return total_det

# Bounds for the parameters
parameter_bounds = [(0.01, 25000), (0.01, 250), (0.01, 250), (0.01, 250)]

# Perform optimization
result = differential_evolution(sum_of_determinants, parameter_bounds, strategy='best1bin', maxiter=1000, popsize=15, tol=0.01, mutation=(0.5, 1), recombination=0.7)

# Print the result
print("Optimal parameters:", result.x)
print("Minimum sum of absolute determinants:", result.fun)