I have a problem to optimize the motion of a spacecraft. It starts relative to point O, moving in an undisturbed near-Earth circular orbit with a radius of 6871 km. The gravitational parameter of the Earth is equal to 3.9860044∙10^14 m^3/s^2. The problem is to fly in dt = 86400 seconds from a point in phase space r0 = (10 km, 100 km, -5 km)^T, v0 = (1 m/s, -10 m/s, 3 m/s)^T to the origin, that is, to the point rf = 0, vf = 0. The initial mass of the spacecraft is equal to 1000 kg, the specific impulse of the spacecraft propulsion system is 220 seconds (2157.463 m/s), and the thrust is varied in range from 0.362 to 100 N. I am trying to minimize the propellant mass index to carry out the meeting between two spacecraft to solve the problem of orbital transfer and spacecraft meeting. Consider the HCW equations as a mathematical model.

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

# Создание модели GEKKO
m = GEKKO()

# Дискретизация времени
nt = 201
m.time = np.linspace(0, 86400, nt)

# Константы
mu = 3.9860044e14  # гравитационный параметр Земли, м^3/с^2
R = 6871e3  # радиус орбиты, м
n = np.sqrt(mu / R**3)  # угловая скорость орбитального движения, рад/с
Isp = 2157.463  # удельный импульс, с
g0 = 9.81  # ускорение свободного падения, м/с^2

# Начальные условия
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  # начальная масса, кг

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', lb=10)
thrust = m.MV(value=50, lb=0.362, ub=100, name='thrust')
thrust.STATUS = 1
thrust.DMAX = 20
theta = [m.MV(value=1.0, lb=-np.pi, ub=np.pi, name=f'theta{i}') for i in range(3)]
for i in range(3):
    theta[i].STATUS = 1
    theta[i].DMAX = 0.1

# Уравнения движения Хилла
m.Equation(r[0].dt() == v[0])
m.Equation(v[0].dt() == 2 * n * v[1] + 3 * n**2 * r[0])
m.Equation(r[1].dt() == v[1])
m.Equation(v[1].dt() == -2 * n * v[0])
m.Equation(r[2].dt() == v[2])
m.Equation(v[2].dt() == -n**2 * r[2])

# Массовый расход по формуле Циолковского
delta_v = thrust / mass
m.Equation(mass.dt() == -delta_v / (Isp * g0))

# Конечные условия
#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])

# Final conditions
p = np.zeros(nt); p[-1]=1
final = m.Param(p)
for i in range(3):
    m.fix_final(r[i], val=rf[i])
    m.fix_final(v[i], val=vf[i])
    m.Minimize(1e-6*(r[i]-rf[i])**2)
    m.Minimize(1e-6*(v[i]-vf[i])**2)

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

# Целевая функция: минимизация конечной массы
#m.Obj(mass)

# Решение задачи оптимизации
#m.options.IMODE = 6
#m.solve(disp=True)

# Solve the optimization problem
m.options.IMODE = 6
m.options.NODES = 2
m.options.RTOL = 1e-6
m.options.OTOL = 1e-6
m.options.SOLVER = 3
m.options.MAX_ITER = 50

# initialize
m.options.COLDSTART = 1
m.solve(disp=True)

# optimize
m.options.TIME_SHIFT = 0
m.options.COLDSTART = 0
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)

# Построение траектории в 3D
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# Построение траектории
#ax.plot(r_data[:, 0], r_data[:, 1], r_data[:, 2], label='Trajectory', color='r')

# Plot the initial and final positions
ax.scatter(r0[0], r0[1], r0[2], color='r', label='Initial Position')
ax.scatter(r[0].value, r[1].value, r[2].value, color='b', label='Trajectory')
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.figure(1,figsize=(6,4))
plt.subplot(3,1,1)
plt.plot(m.time,thrust.value,'r-',label='thrust')
plt.plot(m.time,theta[0].value,'b:',label='theta[0]')
plt.plot(m.time,theta[1].value,'g--',label='theta[1]')
plt.plot(m.time,theta[2].value,'k-',label='theta[2]')
plt.grid(); plt.legend(); plt.ylabel('MVs')
plt.subplot(3,1,2)
plt.plot(m.time,r[0].value,'b:',label='r[0]')
plt.plot(m.time,r[1].value,'g--',label='r[1]')
plt.plot(m.time,r[2].value,'k-',label='r[2]')
plt.grid(); plt.legend(); plt.ylabel('Position')
plt.subplot(3,1,3)
plt.plot(m.time,v[0].value,'b:',label='v[0]')
plt.plot(m.time,v[1].value,'g--',label='v[1]')
plt.plot(m.time,v[2].value,'k-',label='v[2]')
plt.grid(); plt.legend(); plt.ylabel('Velocity')
plt.xlabel('Time (sec)')
plt.tight_layout()

plt.show()

# 2D Построение
fig, ax = plt.subplots(1, 3, figsize=(18, 6))
ax[0].plot(r_data[:, 0], r_data[:, 1], label='Trajectory')
ax[0].scatter(r0[0], r0[1], color='g', label='Initial Position')
ax[0].scatter(rf[0], rf[1], color='m', label='Final Position')
ax[0].set_xlabel('X (km)')
ax[0].set_ylabel('Y (km)')
ax[0].set_title('Trajectory in X-Y plane')
ax[0].legend()

ax[1].plot(r_data[:, 0], r_data[:, 2], label='Trajectory')
ax[1].scatter(r0[0], r0[2], color='g', label='Initial Position')
ax[1].scatter(rf[0], rf[2], color='m', label='Final Position')
ax[1].set_xlabel('X (km)')
ax[1].set_ylabel('Z (km)')
ax[1].set_title('Trajectory in X-Z plane')
ax[1].legend()

ax[2].plot(r_data[:, 1], r_data[:, 2], label='Trajectory')
ax[2].scatter(r0[1], r0[2], color='g', label='Initial Position')
ax[2].scatter(rf[1], rf[2], color='m', label='Final Position')
ax[2].set_xlabel('Y (km)')
ax[2].set_ylabel('Z (km)')
ax[2].set_title('Trajectory in Y-Z plane')
ax[2].legend()

plt.show()

I got results, but the spaceship does not reach the end point, i.e. the solution is not optimal. I need the spaceship to leave the starting point to the rendezvous point (end point) in the given flight time.Can you help me with this?