查找火箭和飞碟的拦截 - TypeError

如何解决查找火箭和飞碟的拦截 - TypeError

对于任务，我需要计算发射火箭的角度β。 这枚火箭应该在空中击中飞碟，并且火箭在前 0.5 秒内刚刚加速。 在练习的第一部分中，我已经计算了飞碟的轨迹。在给定参数旁边还有 x 方向的风。 我试图通过最小化两个轨迹距离的包装函数来计算角度β。 我总是得到 TypeError: 循环的 ufunc 不支持没有可调用 sqrt 方法的 int 类型的参数 0 我不知道我的错误是什么，我希望有人可以帮助我。 实验看起来像这样： Picture - rocket - skeet

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

v_wind = 3.5 #m / s in positive x direction
rho_l = 1.2 # kg/m3

# skeet params
cw_skeet_top = 1.11
cw_skeet_side = 0.5
A_skeet_top = 0.0086
A_skeet_side = 0.00432
m_skeet = 0.105 #kg

v0_skeet = 40
x0_skeet = 60
y0_skeet = 0

# rocket params
cw_rocket_side = 1
cw_rocket_front = 0.5
A_rocket_side = 0.08
A_rocket_front = 0.0312
m_rocket = 2

v0_rocket = 0
x0_rocket = 20
y0_rocket = 0
F_accel_rocket = 120 #in Newton,divide by m_rocket to find accel in m/s^2
t_accel_end = 0.5

def drag(v_rel,rho_l,cw,A):
    Fw = rho_l * cw * A * v_rel**2/2.0
    if (v_rel < 0):
        Fw = -Fw
    return Fw

def skeet(t,y,cw_skeet_side,cw_skeet_top,A_skeet_side,A_skeet_top,m_skeet):
    # Zur Übersichtlichkeit:
    # y[0] - x-position
    # y[1] - x-veLocity
    # y[2] - y-position
    # y[3] - y-veLocity

    d2xdt2 = -drag(y[1]-v_wind,A_skeet_side)/m_skeet   
    dxdt = -y[1]    
    d2ydt2 = -9.81 - drag(y[3],A_skeet_top)/m_skeet
    dydt = y[3]
    
    return dxdt,d2xdt2,dydt,d2ydt2

def hit_ground(t,*args):
    return y[2]
hit_ground.terminal = True
hit_ground.direction = -1

final_x = 0 #m
DGLargs=(rho_l,m_skeet)

def wrapperAlpha(alpha,final_x,DGLargs):
    tspan = [0,100]
    maxStep=tspan[1]/1000  
    start_vec = [x0_skeet,v0_skeet*np.cos(alpha/180*np.pi),y0_skeet,v0_skeet*np.sin(alpha/180*np.pi)]
    sol = solve_ivp(skeet,tspan,start_vec,args=(DGLargs),max_step=maxStep,events=(hit_ground))
    return np.abs(sol.y[0,-1] - final_x)

starting_guess = 10 #deg
res = fsolve(wrapperAlpha,x0 = starting_guess,args = (final_x,DGLargs))

angle = res
print(angle)

tspan = [0,100]
maxStep=tspan[1]/1000 

start_vec = [x0_skeet,v0_skeet*np.cos(angle/180*np.pi),v0_skeet*np.sin(angle/180*np.pi)]
sol = solve_ivp(skeet,events=(hit_ground))
x_skeet = sol.y[0,:]
y_skeet = sol.y[2,:]
plt.plot(x_skeet,y_skeet)

#important exercise:
RocketArgs = (rho_l,cw_rocket_side,cw_rocket_front,A_rocket_side,A_rocket_front,m_rocket,F_accel_rocket,t_accel_end)

def rocket(t,beta,t_accel_end):
    if t_accel_end > t:
        d2xdt2 = (-drag(y[1]-v_wind,A_rocket_side)+(F_accel_rocket*np.cos(beta/180*np.pi)))/m_rocket  
        dxdt = y[1]    
        d2ydt2 = (-9.81 - drag(y[3],A_rocket_front)+(F_accel_rocket*np.sin(beta/180*np.pi)))/m_rocket
        dydt = y[3]
    else:
        d2xdt2 = (-drag(y[1]-v_wind,A_rocket_side))/m_rocket  
        dxdt = y[1]    
        d2ydt2 = (-9.81 - drag(y[3],A_rocket_front))/m_rocket
        dydt = y[3]
        
    
    return dxdt,d2ydt2

def wrapperRocket(beta,x_skeet,y_skeet,RocketArgs):
    tspan = [0,100]
    maxStep = 0.01
    start_vec = [x0_rocket,v0_rocket*np.cos(beta/180*np.pi),y0_rocket,v0_rocket*np.sin(beta/180*np.pi)]
    sol2 = solve_ivp(rocket,args=(beta,t_accel_end),dense_output=True)
    for i in range(len(x_skeet)):
        distance = np.sqrt((sol2.y[0,i]-x_skeet)**2+(sol2.y[2,i]-y_skeet)**2)
    return distance 

starting_guess = 20 #deg
res = minimize(wrapperRocket,args = (x_skeet,RocketArgs))

angle = res
print(angle)