求解物理系统运动方程ODE时遇TypeError问题求助

报错信息 #

Cell In[8], line 8, in f(X, t, m, g, alpha, I1, I3, d1, d2, r1, r2, mu)
      6 def f(X, t, m, g, alpha, I1, I3, d1, d2, r1, r2, mu):
      7     # Unpack the variables from X
----> 8     xCM, yCM, phi, omega3, N1, N2 = X
     11     # Calculate the k values for each wheel
     12     k1x = -mu * np.linalg.norm(np.arctan((yCM + r1*omega3*np.cos(phi) - d1*phi_dot*np.cos(phi)*np.sin(psi)) / (xCM - r1*omega3*np.sin(phi) + d1*phi_dot*np.sin(phi)*np.sin(psi)))) * np.sign(xCM - r1*omega3*np.sin(phi) + d1*phi_dot*np.sin(phi)*np.sin(psi))

TypeError: cannot unpack non-iterable float object

原始代码 #

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt


def f(X, t, m, g, alpha, I1, I3, d1, d2, r1, r2, mu):
    # Unpack the variables from X
    xCM, yCM, phi, omega3, N1, N2 = X
    
    
    # Calculate the k values for each wheel
    k1x = -mu * np.linalg.norm(np.arctan((yCM + r1*omega3*np.cos(phi) - d1*phi_dot*np.cos(phi)*np.sin(psi)) / (xCM - r1*omega3*np.sin(phi) + d1*phi_dot*np.sin(phi)*np.sin(psi)))) * np.sign(xCM - r1*omega3*np.sin(phi) + d1*phi_dot*np.sin(phi)*np.sin(psi))
    k1y = -mu * np.linalg.norm(np.arctan((yCM + r1*omega3*np.cos(phi) - d1*phi_dot*np.cos(phi)*np.sin(psi)) / (xCM - r1*omega3*np.sin(phi) + d1*phi_dot*np.sin(phi)*np.sin(psi)))) * np.sign(xCM + r1*omega3*np.sin(phi) - d1*phi_dot*np.sin(phi)*np.sin(psi))
    k2x = -mu * np.linalg.norm(np.arctan((yCM + r2*omega3*np.cos(phi) + d2*phi_dot*np.cos(phi)*np.sin(psi)) / (xCM - r2*omega3*np.sin(phi) - d2*phi_dot*np.sin(phi)*np.sin(psi)))) * np.sign(xCM - r2*omega3*np.sin(phi) - d2*phi_dot*np.sin(phi)*np.sin(psi))
    k2y = -mu * np.linalg.norm(np.arctan((yCM + r2*omega3*np.cos(phi) + d2*phi_dot*np.cos(phi)*np.sin(psi)) / (xCM - r2*omega3*np.sin(phi) - d2*phi_dot*np.sin(phi)*np.sin(psi)))) * np.sign(xCM + r2*omega3*np.sin(phi) + d2*phi_dot*np.sin(phi)*np.sin(psi))

    xCM_dot = phi_dot*(r1*k1y - r2*k2y)/m
    yCM_dot = (m*g*np.sin(alpha) + N1*k1x + N2*k2x - phi_dot**2*(r1*k1x + r2*k2x))/m
    phi_dot = omega3
    psi_dot = omega3*np.cos(phi)/np.sin(psi)
    omega3_dot = ((r1*k1y - r2*k2y)*np.cos(phi) - (r1*k1x + r2*k2x)*np.sin(phi) + (d1*k1y - d2*k2y)*np.cos(psi)*np.sin(phi) - (d1*k1x + d2*k2x)*np.sin(psi)*np.sin(phi) - I1*phi_dot*omega3*np.sin(phi) + I3*psi_dot*omega3*np.cos(phi))/I3

    # Return the derivatives as a numpy array
    return np.array([xCM_dot, yCM_dot, phi_dot, omega3_dot, N1, N2])

X0 = np.array([0, 0, 0, 0, 0, 0]) # initial conditions for [xCM, yCM, phi, psi, omega3]
t_span = [0, 10] # time interval to solve over
t_eval = np.linspace(t_span[0], t_span[1], 1000) # times to evaluate the solution at

m = 1 # mass of the unicycle
g = 9.81 # acceleration due to gravity
alpha = np.pi/4 # slope angle
I1 = 0.1 # moment of inertia about the x-axis
I3 = 0.2 # moment of inertia about the z-axis
d1 = 0.3 # distance between the center of mass and the axis of rotation for the left wheel
d2 = 0.3 # distance between the center of mass and the axis of rotation for the right wheel
r1 = 0.3 # radius of the left wheel
r2 = 0.3 # radius of the right wheel
mu = 0.2 # coefficient of friction
psi = 0.3

sol = solve_ivp(f, t_span, X0, t_eval=t_eval, args=(m, g, alpha, I1, I3, d1, d2, r1, r2, mu))

问题分析 #

1. 参数顺序错误：scipy.integrate.solve_ivp要求ODE函数的第一个参数是时间t，第二个是状态变量X，但原函数f将X放在首位，导致solve_ivp传入的单值时间t被当成多元素状态变量解包，触发类型错误。
2. 状态变量逻辑混乱：
   • 初始条件注释标注为5个变量，但实际定义了6个初始值，解包时遗漏了psi，同时将约束反力N1、N2作为状态变量（此类变量无需积分，应通过动力学方程求解）。
3. 未定义变量引用：计算k_i时提前使用了phi_dot，但该变量在解包后才定义，会触发未定义错误；psi作为全局变量传入，不符合函数封装规范。
4. 冗余计算：np.arctan返回标量，np.linalg.norm对单值无意义，属于冗余调用。

修正方案 #

1. 调整ODE函数参数顺序为f(t, X, ...)，匹配solve_ivp的要求。
2. 明确状态变量：将psi加入状态变量，移除N1、N2，通过静力学平衡方程求解约束反力。
3. 提前定义phi_dot，避免未定义错误。
4. 移除冗余的np.linalg.norm调用，简化k_i计算。
5. 将psi作为参数传入函数，避免全局依赖。

修正后代码 #

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt


def f(t, X, m, g, alpha, I1, I3, d1, d2, r1, r2, mu, psi):
    # 解包状态变量：xCM, yCM, phi, omega3, psi
    xCM, yCM, phi, omega3, psi = X
    
    # 提前定义导数变量，避免未定义错误
    phi_dot = omega3
    # 避免除以零的异常处理
    psi_dot = omega3 * np.cos(phi) / np.sin(psi) if np.abs(np.sin(psi)) > 1e-6 else 0

    # 计算k值，移除冗余的np.linalg.norm
    denominator1 = xCM - r1 * omega3 * np.sin(phi) + d1 * phi_dot * np.sin(phi) * np.sin(psi)
    numerator1 = yCM + r1 * omega3 * np.cos(phi) - d1 * phi_dot * np.cos(phi) * np.sin(psi)
    arctan_val1 = np.arctan(numerator1 / denominator1) if np.abs(denominator1) > 1e-6 else 0
    
    k1x = -mu * arctan_val1 * np.sign(xCM - r1 * omega3 * np.sin(phi) + d1 * phi_dot * np.sin(phi) * np.sin(psi))
    k1y = -mu * arctan_val1 * np.sign(xCM + r1 * omega3 * np.sin(phi) - d1 * phi_dot * np.sin(phi) * np.sin(psi))

    denominator2 = xCM - r2 * omega3 * np.sin(phi) - d2 * phi_dot * np.sin(phi) * np.sin(psi)
    numerator2 = yCM + r2 * omega3 * np.cos(phi) + d2 * phi_dot * np.cos(phi) * np.sin(psi)
    arctan_val2 = np.arctan(numerator2 / denominator2) if np.abs(denominator2) > 1e-6 else 0
    
    k2x = -mu * arctan_val2 * np.sign(xCM - r2 * omega3 * np.sin(phi) - d2 * phi_dot * np.sin(phi) * np.sin(psi))
    k2y = -mu * arctan_val2 * np.sign(xCM + r2 * omega3 * np.sin(phi) + d2 * phi_dot * np.sin(phi) * np.sin(psi))

    # 通过静力学平衡求解约束反力N1、N2（可根据实际物理模型调整）
    N1 = (m * g * np.cos(alpha) * d2) / (d1 + d2)
    N2 = (m * g * np.cos(alpha) * d1) / (d1 + d2)

    # 计算各状态变量的导数
    xCM_dot = phi_dot * (r1 * k1y - r2 * k2y) / m
    yCM_dot = (m * g * np.sin(alpha) + N1 * k1x + N2 * k2x - phi_dot**2 * (r1 * k1x + r2 * k2x)) / m
    omega3_dot = ((r1 * k1y - r2 * k2y) * np.cos(phi) - (r1 * k1x + r2 * k2x) * np.sin(phi) + 
                  (d1 * k1y - d2 * k2y) * np.cos(psi) * np.sin(phi) - 
                  (d1 * k1x + d2 * k2x) * np.sin(psi) * np.sin(phi) - 
                  I1 * phi_dot * omega3 * np.sin(phi) + I3 * psi_dot * omega3 * np.cos(phi)) / I3

    # 返回导数数组，顺序与状态变量完全匹配
    return np.array([xCM_dot, yCM_dot, phi_dot, omega3_dot, psi_dot])

# 初始条件：[xCM, yCM, phi, omega3, psi]
X0 = np.array([0, 0, 0, 0, 0.3])
t_span = [0, 10]
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# 物理参数
m = 1
g = 9.81
alpha = np.pi/4
I1 = 0.1
I3 = 0.2
d1 = 0.3
d2 = 0.3
r1 = 0.3
r2 = 0.3
mu = 0.2

# 求解ODE，将psi作为参数传入
sol = solve_ivp(f, t_span, X0, t_eval=t_eval, args=(m, g, alpha, I1, I3, d1, d2, r1, r2, mu, X0[4]))

# 验证求解结果并绘图
print(sol.success)
if sol.success:
    plt.figure(figsize=(12, 8))
    plt.subplot(2,2,1)
    plt.plot(sol.t, sol.y[0])
    plt.title('xCM vs Time')
    plt.xlabel('t')
    plt.ylabel('xCM')
    
    plt.subplot(2,2,2)
    plt.plot(sol.t, sol.y[1])
    plt.title('yCM vs Time')
    plt.xlabel('t')
    plt.ylabel('yCM')
    
    plt.subplot(2,2,3)
    plt.plot(sol.t, sol.y[2])
    plt.title('phi vs Time')
    plt.xlabel('t')
    plt.ylabel('phi')
    
    plt.subplot(2,2,4)
    plt.plot(sol.t, sol.y[3])
    plt.title('omega3 vs Time')
    plt.xlabel('t')
    plt.ylabel('omega3')
    
    plt.tight_layout()
    plt.show()

内容的提问来源于stack exchange，提问作者Gustavo Valente