我正在尝试解决逆弹道问题，其中必须在给定到达点
(az, el)
和初始弹射速度
posTarget
的情况下找到初始发射方向
v0
为此，我有几个类：
此类编译后续计算中所需的一些物理常量
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
class Constant:
_g: tuple[float, float, float] = (0.0, 0.0, -sp.constants.g)
_vSound: float = 340.0
_density: float = 1.22
此类包含后续计算中使用的一些实用程序
class Utils:
def __init__(self) -> None:
pass
@staticmethod
def make_vec(p1: tuple[float, float, float],
p2: tuple[float, float, float],
) -> tuple[float, float, float]:
return tuple(p1val - p2val for p1val, p2val in zip(p1, p2))
@classmethod
def vec_norm(cls, v: tuple[float, float, float]) -> float:
return np.linalg.norm(v)
@classmethod
def get_dist(cls,
p1: tuple[float, float, float],
p2: tuple[float, float, float],
) -> float:
return cls.vec_norm(cls.make_vec(p1 = p1, p2 = p2))
@classmethod
def speed_to_Ma(cls, v: tuple[float, float, float]) -> float:
return cls.vec_norm(v) / Constant._vSound
@staticmethod
def get_LOS_angles(p1: tuple[float, float, float],
p2: tuple[float, float, float],
) -> tuple[float, float]:
az: float = np.arctan2( (p2[1] - p1[1]), (p2[0] - p1[0]) )
el: float = np.arctan2((p2[2] - p1[2]), (p2[1] - p1[1]))
return az, el
此类代表传播的射弹
class Projectile:
def __init__(self,
m: float = None,
cal: float = None,
S: float = None,
CD0subsonic: float = None,
CD0max: float = None,
k: float = None,
) -> None:
self._m: float = m
self._cal: float = cal
self._CD0subsonic: float = CD0subsonic
self._CD0max: float = CD0max
self._k: float = k
if S is None:
self._S: float = .25 * np.pi * self._cal**2
else:
self._S: float = S
self._x: list[float] = []
self._y: list[float] = []
self._z: list[float] = []
self._pos: list[tuple[float, float, float]] = list(zip(self._x, self._y, self._z))
self._vx: list[float] = []
self._vy: list[float] = []
self._vz: list[float] = []
self._v: list[tuple[float, float, float]] = list(zip(self._vx, self._vy, self._vz))
@property
def m(self) -> float:
return self._m
@m.setter
def m(self, m: float) -> None:
self._m: float = m
@property
def cal(self) -> float:
return self._cal
@cal.setter
def cal(self, cal: float) -> None:
self._cal: float = cal
@property
def CD0subsonic(self) -> float:
return self._CD0subsonic
@CD0subsonic.setter
def CD0subsonic(self, CD0subsonic: float) -> None:
self._CD0subsonic: float = CD0subsonic
@property
def CD0max(self) -> float:
return self._CD0max
@CD0max.setter
def CD0max(self, CD0max: float) -> None:
self._CD0max: float = CD0max
@property
def k(self) -> float:
return self._k
@k.setter
def k(self, k: float) -> None:
self._k: float = k
@property
def S(self) -> float:
return self._S
@S.setter
def S(self, S: float) -> None:
self._S: float = S
@property
def x(self) -> list[float]:
return self._x
@x.setter
def x(self, x: float) -> None:
self._x.append(x)
def get_xidx(self, i: int) -> float:
return self.x[i]
def set_xidx(self, i: int, x: float) -> None:
self.x[i] = x
@property
def y(self) -> list[float]:
return self._y
@y.setter
def y(self, y: float) -> None:
self._y.append(y)
def get_yidx(self, i: int) -> float:
return self.y[i]
def set_yidx(self, i: int, y: float) -> None:
self.y[i] = y
@property
def z(self) -> list[float]:
return self._z
@z.setter
def z(self, z: float) -> None:
self._z.append(z)
def get_zidx(self, i: int) -> float:
return self.z[i]
def set_zidx(self, i: int, z: float) -> None:
self.z[i] = z
@property
def pos(self) -> list[tuple[float, float, float]]:
return list(zip(self._x, self._y, self._z))
@pos.setter
def pos(self, pos: tuple[float, float, float]) -> None:
self._pos.append(pos)
def get_posidx(self, i: int) -> tuple[float, float, float]:
return self.pos[i]
def set_posidx(self, i: int, pos: tuple[float, float, float]) -> None:
self.pos[i] = pos
@property
def vx(self) -> list[float]:
return self._vx
@vx.setter
def vx(self, vx: float) -> None:
self._vx.append(vx)
def get_vxidx(self, i: int) -> float:
return self.vx[i]
def set_vxidx(self, i: int, vx: float) -> None:
self.vx[i] = vx
@property
def vy(self) -> list[float]:
return self._vy
@vy.setter
def vy(self, vy: float) -> None:
self._vy.append(vy)
def get_vyidx(self, i: int) -> float:
return self.vy[i]
def set_vyidx(self, i: int, vy: float) -> None:
self.vy[i] = vy
@property
def vz(self) -> list[float]:
return self._vz
@vz.setter
def vz(self, vz: float) -> None:
self._vz.append(vz)
def get_vzidx(self, i: int) -> float:
return self.vz[i]
def set_vzidx(self, i: int, vz: float) -> None:
self.vz[i] = vz
@property
def v(self) -> list[tuple[float, float, float]]:
return list(zip(self._vx, self._vy, self._vz))
@v.setter
def v(self, v: tuple[float, float, float]) -> None:
self._v.append(v)
def get_vidx(self, i: int) -> tuple[float, float, float]:
return self.v[i]
def set_vidx(self, i: int, v: tuple[float, float, float]) -> None:
self.v[i] = v
和此类包含传播所需的方法和函数
class PMM:
def __init__(self,
projectile: Projectile = None,
tFinal: float = None,
dt: float = None,
nIterMax: int = None,
) -> None:
self._projectile: Projectile = projectile
self._tFinal: float = tFinal
self._dt: float = dt
self._nIterMax: int = nIterMax
@property
def projectile(self) -> Projectile:
return self._projectile
@projectile.setter
def projectile(self, projectile: Projectile) -> None:
self._projectile: Projectile = projectile
@property
def tFinal(self) -> float:
return self._tFinal
@tFinal.setter
def tFinal(self, tFinal: float) -> None:
self._tFinal: float = tFinal
@property
def dt(self) -> float:
return self._dt
@dt.setter
def dt(self, dt: float) -> None:
self._dt: float = dt
@property
def nIterMax(self) -> int:
return self._nIterMax
@nIterMax.setter
def nIterMax(self, nIterMax: int) -> None:
self._nIterMax: int = nIterMax
def gravity_force(self) -> tuple[float, float, float]:
return tuple(self.projectile.m * gval for gval in Constant._g)
def drag_force(self,
v: tuple[float, float, float]
) -> tuple[float, float, float]:
return tuple(
-0.5
* Constant._density
* self.projectile.S
* self.CD0(v = v)
* Utils.vec_norm(v)
* vval for vval in v
)
def CD0(self, v: tuple[float, float, float]) -> float:
vMa: float = Utils.speed_to_Ma(v)
if vMa < 1.0:
return self.projectile.CD0subsonic
elif vMa == 1.0:
return self.projectile.CD0max
else:
return (self.projectile.CD0max / (vMa ** self.projectile.k))
def init_proj_pos(self, pos: tuple[float, float, float]) -> None:
self.projectile.x = pos[0]
self.projectile.y = pos[1]
self.projectile.z = pos[2]
def init_proj_vel(self, v0: float, az: float, el: float) -> None:
self.projectile.vx = v0 * np.cos(az) * np.cos(el)
self.projectile.vy = v0 * np.sin(az) * np.cos(el)
self.projectile.vz = v0 * np.sin(el)
def run_bare(self,
target: Target = None,
pos: tuple[float, float, float] = None,
tolDist: float = None,
v0: float = None,
az: float = None,
el: float = None,
) -> float:
posTarget: tuple[float, float, float] = target.get_posidx(i = 0)
dt: float = self.dt
t: float = dt
self.projectile.x.clear()
self.projectile.y.clear()
self.projectile.z.clear()
self.init_proj_pos(pos = pos)
self.projectile.v.clear()
self.init_proj_vel(v0 = v0, az = az, el = el)
gForce: tuple[float, float, float] = self.gravity_force()
i: int = 1
nIterMax: int = self.nIterMax
while (t <= self.tFinal and i < nIterMax):
vOld: tuple[float, float, float] = self.projectile.get_vidx(i = i - 1)
posOld: tuple[float, float, float] = self.projectile.get_posidx(i = i - 1)
cd0Force: tuple[float, float, float] = self.drag_force(vOld)
v = tuple(
vval + dt * (cd0 + gval) / self.projectile.m
for vval, cd0, gval in zip(vOld, cd0Force, gForce)
)
self.projectile.vx = v[0]
self.projectile.vy = v[1]
self.projectile.vz = v[2]
pos = tuple(x + dt * vval for x, vval in zip(posOld, v))
self.projectile.x = pos[0]
self.projectile.y = pos[1]
self.projectile.z = pos[2]
distOldY: float = np.abs(posOld[1] - posTarget[1])
distNewY: float = np.abs(pos[1] - posTarget[1])
deltaDistY: float = distNewY - distOldY
if distNewY < tolDist:
z: float = self.projectile.get_zidx(i = -1)
self.projectile.x.clear()
self.projectile.y.clear()
self.projectile.z.clear()
self.projectile.v.clear()
return z - posTarget[2]
else:
if deltaDistY < 0.0:
i += 1
t += dt
continue
else:
del self.projectile.pos[i]
del self.projectile.v[i]
dt *= .1
continue
i += 1
t += dt
self.projectile.x.clear()
self.projectile.y.clear()
self.projectile.z.clear()
self.projectile.v.clear()
return np.inf
def run(self,
pos: tuple[float, float, float],
v0: float,
az: float,
el: float,
) -> None:
self.projectile.x.clear()
self.projectile.y.clear()
self.projectile.z.clear()
self.projectile.v.clear()
t: float = self.dt
self.init_proj_pos(pos = pos)
self.init_proj_vel(v0 = v0, az = az, el = el)
print(f"el (run): {np.degrees(el)}")
gForce: tuple[float, float, float] = self.gravity_force()
i: int = 1
while (t <= self.tFinal):
vOld: tuple[float, float, float] = self.projectile.get_vidx(i = i - 1)
posOld: tuple[float, float, float] = self.projectile.get_posidx(i = i - 1)
cd0Force: tuple[float, float, float] = self.drag_force(vOld)
v = tuple(
vval + self.dt * (cd0 + gval) / self.projectile.m
for vval, cd0, gval in zip(vOld, cd0Force, gForce)
)
self.projectile.vx = v[0]
self.projectile.vy = v[1]
self.projectile.vz = v[2]
pos = tuple(x + self.dt * vval for x, vval in zip(posOld, v))
self.projectile.x = pos[0]
self.projectile.y = pos[1]
self.projectile.z = pos[2]
t += self.dt
i += 1
def get_init_orientation(self,
target: Target,
pos: tuple[float, float, float],
tolAng: float,
tolDist: float,
v0: float,
) -> tuple[float, float]:
pos_proj: tuple[float, float, float] = self.projectile.get_posidx(i = 0)
pos_targ: tuple[float, float, float] = target.get_posidx(i = 0)
az, el = Utils.get_LOS_angles(p1 = pos_proj, p2 = pos_targ)
el_min, el_max = 0.3 * el, 1.5 * el
el = .5 * (el_min + el_max)
nIter: int = 0
distZ = self.run_bare(target = target,
pos = pos,
tolDist = tolDist,
v0 = v0,
az = az,
el = el,
)
while (abs(distZ) > tolDist) and (abs(el - el_min) > tolAng):
if distZ < 0.0:
el_min = el
el = .5 * (el_min + el_max)
else:
el_max = el
el = .5 * (el_min + el_max)
distZ = self.run_bare(target = target,
pos = pos,
tolDist = tolDist,
v0 = v0,
az = az,
el = el,
)
nIter += 1
return (az, el)
然后通过
if __name__ == "__main__":
t: Target = Target()
tInit: tuple[float, float, float] = (50.0, 50.0, 50.0)
t.x = tInit[0]
t.y = tInit[1]
t.z = tInit[2]
p: Projectile = Projectile(m = 4.04e-3,
cal = 5.56e-3,
CD0subsonic = .15,
CD0max = .85,
k = 1.0,
)
pInit: tuple[float, float, float] = (0.0, 0.0, 0.0)
p.x = pInit[0]
p.y = pInit[1]
p.z = pInit[2]
az_LOS, el_LOS = Utils.get_LOS_angles(p1 = p.get_posidx(i = 0),
p2 = t.get_posidx(i = 0),
)
dist_LOS: float = Utils.get_dist(p1 = p.get_posidx(i = 0), p2 = t.get_posidx(i = 0))
dt: float = 1e-3
v0: float = 750.0
tFinal: float = 1.5 * (dist_LOS / v0)
nIterMax: int = int(tFinal / dt)
pmm: PMM = PMM(projectile = p,
tFinal = tFinal,
dt = dt,
nIterMax = nIterMax,
)
tolDist: float = 1.0e-3
tolAngle: float = 1.0e-3
az, el = pmm.get_init_orientation(target = t,
pos = pInit,
tolAng = tolAngle,
tolDist = tolDist,
v0 = v0,
)
pmm.run(v0 = v0, az = az, el = el, pos = pInit)
fig = plt.figure()
ax = plt.axes(projection = '3d')
ax.scatter(pmm.projectile.x, pmm.projectile.y, pmm.projectile.z)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.scatter(t.x, t.y, t.z, label = "Target")
plt.legend()
plt.show()
函数执行测试
run_bare
计算给定
(az, el)
的传播，并且如果射弹在该点的一定容差范围内到达（沿
y
轴），它返回目标和射弹之间相应的高度差。该值用于
get_init_orientation
来寻找好的
el
经检查，方法
el
中传递给
run_bare
的
get_init_orientation
值似乎正在变化，但初始速度，即
self.projectile.get_vidx(i = 0)
始终相同，即使方法
init_proj_vel
接收到新的速度
el
在
get_init_elevation
中的每个循环处。因此，代码无法找到好的
el
为什么
el
更新了，但相应的初始弹射速度却没有更新？
代码没有更新
self.projectile.v
数组，因为在
run_bare
函数开始时调用了
self.projectile.v.clear()
，这会清除掉数组中的所有元素，包括你在
init_proj_vel
中新添加的速度。
为了解决这个问题，可以将
self.projectile.v.clear()
和其他类似的清空数组的操作移到
run_bare
函数的外部，例如在
get_init_orientation
函数的开头。这样，每次调用
run_bare
时，它都会使用上一次调用
init_proj_vel
设置的最新速度。
下面是修改后的
get_init_orientation
函数：
def get_init_orientation(self,
target: Target,
pos: tuple[float, float, float],
tolAng: float,
tolDist: float,
v0: float,
) -> tuple[float, float]:
self.projectile.x.clear()
self.projectile.y.clear()
self.projectile.z.clear()
self.projectile.v.clear() # 清空数组的操作移到此处
pos_proj: tuple[float, float, float] = self.projectile.get_posidx(i = 0)
pos_targ: tuple[float, float, float] = target.get_posidx(i = 0)
az, el = Utils.get_LOS_angles(p1 = pos_proj, p2 = pos_targ)
# ... 其他代码 ...
此外，可能还需要检查代码的其他部分，确保没有其他地方意外清空了
self.projectile.v
数组，导致初始速度没有被正确更新。