Задача оптимизации движения космического корабля ⇐ Python

[Anonymous]

У меня задача оптимизировать движение космического корабля. Он стартует относительно точки О, двигаясь по невозмущенной околоземной круговой орбите радиусом 6871 км. Гравитационный параметр Земли равен 3,9860044∙10^14 м^3/с^2. Задача состоит в том, чтобы пролететь за dt = 86400 секунд из точки фазового пространства r0 = (10 км, 100 км, -5 км)^T, v0 = (1 м/с, -10 м/с, 3 м/с). )^Т в начало координат, то есть в точку rf = 0, vf = 0. Начальная масса КА равна 1000 кг, удельный импульс двигательной установки КА 220 секунд (2157,463 м/с). , а тяга варьируется в диапазоне от 0,362 до 100 Н. Пытаюсь минимизировать индекс массы топлива для проведения сближения двух космических кораблей, решения задачи перевода на орбиту и сближения и стыковки космических кораблей.

Код: Выделить всё

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create GEKKO model
m = GEKKO()

# Time discretization
nt = 101
m.time = np.linspace(0, 1, nt)

# Constants
mu = 3.9860044e14
R = 6871e3
Isp = 2157.463
g0 = 9.81

# Initial conditions
r0 = np.array([10e3, 100e3, -5e3])
v0 = np.array([1, -10, 3])
rf = np.array([0, 0, 0])
vf = np.array([0, 0, 0])
m0 = 1000

# Variables
r = [m.Var(value=r0[i], name=f'r{i}') for i in range(3)]
v = [m.Var(value=v0[i], name=f'v{i}') for i in range(3)]
mass = m.Var(value=m0, name='mass')
thrust = m.MV(value=50, lb=0.362, ub=100, name='thrust')
thrust.STATUS = 1
theta = [m.Var(value=0, lb=-np.pi, ub=np.pi, name=f'theta{i}') for i in range(3)]

# Equations of motion
for i in range(3):
m.Equation(r[i].dt() == v[i])
m.Equation(v[i].dt() == (thrust * m.cos(theta[i])) / mass - (mu / R**3) * r[i])

# Mass flow rate
m_flow = thrust / (Isp * g0)
m.Equation(mass.dt() == -m_flow)

# Final conditions
for i in range(3):
m.fix(r[i], pos=nt-1, val=rf[i])
m.fix(v[i], pos=nt-1, val=vf[i])

# Objective: Minimize the final mass
m.Obj(mass)

# Solve the optimization problem
m.options.IMODE = 6
m.solve(disp=True)

# Display results
print(f'Final mass: {mass.value[-1]} kg')
for i in range(3):
print(f'Final position r{i}: {r[i].value[-1]} m')
print(f'Final velocity v{i}: {v[i].value[-1]} m/s')

# Extract trajectory data
r_data = np.zeros((nt, 3))
for i in range(3):
r_data[:, i] = np.array(r[i].value)

# Plot the trajectory in 3D
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# Plot the trajectory
ax.plot(r_data[:, 0], r_data[:, 1], r_data[:, 2], label='Trajectory')

# Plot the Earth as a sphere
u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
x = R * np.cos(u) * np.sin(v)
y = R * np.sin(u) * np.sin(v)
z = R * np.cos(v)
ax.plot_surface(x, y, z, color='w', alpha=0.3)

# Plot the initial and final positions
ax.scatter(r0[0], r0[1], r0[2], color='r', label='Initial Position')
ax.scatter(rf[0], rf[1], rf[2], color='g', label='Final Position')

# Labels and legend
ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Z (km)')
ax.set_title('Spacecraft Trajectory')
ax.legend()

plt.show()

Результаты получены, но на графике я вижу только начальную и конечную точки, но не вижу траектории, описываемой космическим кораблем, выходящим из начальной точки в точку встречи (конец точку).
Не могли бы вы мне с этим помочь?


Подробнее здесь: https://stackoverflow.com/questions/786 ... spacecraft