# Looking for Navigators for the Famous Vega Classic

Do you have what it takes to compete in the famous Vega Classic?

Looking for skilled navigators. Apply below.

You have a craft that is capable of 10m/s$$^2$$ acceleration in any direction, and we're looking for the fasted time from Vega Starship 1 to Vega Starship 2 around the beacon.

The rules of the race:

1. You will start at SS1 (co-ordinates (0,0)), at 0 relative velocity.
2. You can accelerate in any direction (assume you can change direction instantaneously), but the total magnitude must be less than 10m/s$$^2$$.
3. You will go around the beacon (co-ordinates (1e6,0)) and avoid the red shaded area.
4. You will end at SS2 (co-ordinates (1e6,1e8)), at 0 relative velocity.
5. We re-iterate: You must finish at 0 relative velocity! Especially after the smear debacle of last year.

$$[[(t_0,t_1),(a_{x,0},a_{y,0})]$$,
$$[(t_1,t_2),(a_{x,1},a_{y,1})]$$,
...
$$[(t_{n-1},t_n),(a_{x,n-1},a_{y,n-1})]]$$

where $$t_0 = 0$$, $$t_n$$ is your total time, and $$t_k$$ is the time of your $$k$$th adjustment. And the $$k$$th acceleration component is such that $$a_{x,k}^2+a_{y,k}^2\leq a^2$$. Your acceleration can also be provided as an angle in radians, where 0 is to the right, and $$\pi/2$$ is straight up. Such a co-ordinate would look like this: $$[(t_k,t_{k+1}),\theta_k]$$. You can mix and match.

Note, we're assuming non-relativistic velocities here.

Here's python code for evaluating a solution in the given format. Final position should be within 1m of the target, and the velocity should be "close to zero" meaning substantially less than 1.

from math import cos, sin
def race(commands):

L1 = 1e6   #m - horizontal distance
L2 = 1e8   #m - vertical distance
a = 10     #m/s^2

x,y,vx,vy = 0,0,0,0
for cmd in commands:
if type(cmd[1]) is tuple or type(cmd[1]) is list:
#cmd[1] is acceleration components
assert abs(cmd[1][0]**2 + cmd[1][1]**2 - a**2)<0.00001
ax,ay = cmd[1]
else:
#cmd[1] is acceleration angle
ax,ay = a*cos(cmd[1]), a*sin(cmd[1])
burnTime = cmd[0][1]-cmd[0][0]
dt = burnTime / int(burnTime+1)
for i in range(int(burnTime+1)):
x,y,vx,vy = x + vx*dt + 0.5*ax*dt**2,y + vy*dt + 0.5*ay*dt**2,vx + ax*dt,vy + ay*dt

if x<L1 and y>0:
print("Warning:",x,y,vx,vy)

print("Time:",cmd[0][1])
print("Final position & velocity:",x,y,vx,by)


For reference, my current record has a current output of:

Final position & velocity: 1000000.0000009941 100000000.16834144 2.41069386675008e-10 -3.068626813984565e-06

• Does "go around the beacon" have any meaning beyond "avoid the shaded area" (or, more precisely, "avoid the shaded area and go to its southeast rather than its northwest", but the latter is obviously preferable anyway if you want to minimize your time)? That is, are we supposed to pass close to the beacon in any sense? I'm pretty sure the answer implied by the question as currently written, and the adjudication code, is no, but it seems worth checking. – Gareth McCaughan Mar 23 at 23:53
• Why not starting directly to the right? Is this just slower? – Moti Mar 25 at 3:06
• Definitely mathematics rather than physics. Solving it is not a routine exercise, so it ticks at least one of the boxes for being a puzzle in the sense we use around here. Yes, the fastest thing will certainly involve continuously varying thrust, and I think an ideal answer would (if there is one) give an explicit formula for the thrust at given time. I spent a little while trying to find, if not an actual answer, at least constraints that in principle determine the actual answer, but I am fairly sure I did something wrong because [... continues] – Gareth McCaughan Mar 30 at 19:11
• ... I ended up with what seem to be enough constraints that (in a slight generalization of the given problem) there don't seem to be enough remaining degrees of freedom to allow there always to be a solution to the problem :-). – Gareth McCaughan Mar 30 at 19:12
• @GarethMcCaughan that's why good Vega Classic navigators are in high demand – Dr Xorile Apr 1 at 2:58

## Solution 1 (6674.361 seconds): Discrete optimization

Break the trajectory into $$N$$ timesteps of length $$\Delta t$$. The time after the $$i$$-th timestep will be $$t_i=i\Delta t$$, with the starting time $$t_0=0$$ and the final time $$t_N$$. Given an acceleration of $$\vec a_i$$ during the $$i$$-th timestep, the velocity at time $$t_i$$ will be $$\vec v_i=\sum_{k=1}^{i}\vec a_k\Delta t$$ (assuming $$\vec v_0=\vec 0$$). The displacement during each timestep will be the average velocity times $$\Delta t$$: $$\vec x_i-\vec x_{i-1}=(\vec v_i+\vec v_{i-1})\Delta t/2$$. This means the total displacement will be:

\begin{align} \vec x_i &= \sum_{j=1}^i\left[\left(\sum_{k=1}^{j}\vec a_k\Delta t\right)+\left(\sum_{k=1}^{j-1}\vec a_k\Delta t\right)\right]\Delta t/2 \\ &= \sum_{k=1}^i[(2i-2k+1)\vec a_k]\Delta t^2/2 \end{align}

(Again assuming that $$\vec x_0=\vec 0$$.)

I assert without proof that an optimal trajectory will pass through the corner. This means that we can separate the trajectory into two parts or legs. Break the first leg into $$M=\lfloor N/2\rfloor$$ equal timesteps $$\Delta t_1$$, and the second leg into $$N-M$$ equal timesteps $$\Delta t_2$$. Since the initial and final velocities are zero, we can write the velocity at $$t_M$$ as either $$\vec v_M = \sum_{k=1}^{M}\vec a_k \Delta t_1 = -\sum_{k=M+1}^{N}\vec a_k \Delta t_2$$. The displacement of each leg is known and can be calculated as above. Note that in both the velocity and displacement equations, $$\Delta t_i$$ can be factored out. Finally, we restrict $$|\vec a_i|\le a_\textrm{max}$$. This gives us the following system of equations:

Minimize: $$t_N=M\Delta t_1+(N-M)\Delta t_2$$

Subject to: $$|\vec a_i| \le a_\textrm{max},\ i\in[1,N]\\\\$$ \begin{align} \left(\frac{\Delta \vec v_1}{\Delta t_1}\right) &= \sum_{k=1}^{M}\vec a_k \\ \left(\frac{\Delta \vec v_2}{\Delta t_2}\right) &= \sum_{k=M+1}^{N}\vec a_k \\ \left(\frac{\Delta \vec x_1}{\Delta t_1^2}\right) &= \sum_{k=1}^{M}\left(M-k+\frac{1}{2}\right)\vec a_k \\ \left(\frac{\Delta \vec x_2}{\Delta t_2^2}\right) &= \sum_{k=M+1}^{N}\left(M-k+\frac{1}{2}\right)\vec a_k \\ \end{align} $$\left(\frac{\Delta \vec v_1}{\Delta t_1}\right)\Delta t_1 = \left(\frac{\Delta \vec v_2}{\Delta t_2}\right)\Delta t_2$$ \begin{align} \left(\frac{\Delta \vec x_1}{\Delta t_1^2}\right)\Delta t_1^2 &= (10^6,0) \\ \left(\frac{\Delta \vec x_2}{\Delta t_2^2}\right)\Delta t_2^2 &= (0,10^8) \\ \end{align}

(The factor of $$\vec a_k$$ in the sum for $$\Delta\vec x_2$$ comes from the fact that we are summing in reverse from $$N$$ to $$M+1$$, so the factor is $$(N-M)-j+1/2$$ where $$j=N-k+1$$. This is equal to ($$N-M-N+k-1+1/2=-M+k-1/2$$), but since we are calculating backwards in time we need to negate the result, coincidentally giving us $$M-k+1/2$$. This is only possible since we know $$\vec v_N=\vec 0$$.)

The reason for factoring out $$\Delta t_i$$ is that this system has only six quadratic constraints, plus six linear constraints and $$N$$ convex cone constraints. This simplification lets a solver like SCIP handle the problem efficiently. Here is an example specification of the problem and a solution found by SCIP. And here is the solution in the format the validator expects; the validator output is:

Time: 6674.36
Final position & velocity: 1000048.6673564073 99999992.4373888 0.0033702956625456526 0.009344372255398525


Note my own calculations give a final position of (999 999.988, 99 999 999.984), well within the 1 meter tolerance, and a final velocity of (9.8×10-7, -7.9×10-7), which is 'substantially' less than 1. The position at $$t_M$$ is (999 999.982, -9.1×10-12), so I am technically cutting the corner, but only by 18 mm.

## Solution 2 (6674.291 seconds): Calculus of Variations

Following the formulation in D.G. Stechert, "On the Use of the Calculus of Variations in Trajectory Optimization Problems," 1963 (pdf), we can set up our equations of motion in terms of $$x$$ and $$y$$:

\begin{align} \phi_1 &= \dot x - u = 0 \\ \phi_2 &= \dot y - v = 0 \\ \phi_3 &= \dot u - a\cos\theta = 0 \\ \phi_4 &= \dot v - a\sin\theta = 0 \\ \end{align}

(Where $$u$$ and $$v$$ are the $$x$$- and $$y$$-velocities, respectively). We then set up our Lagrange multipliers:

$$F = \lambda_1 \phi_1 + \lambda_2 \phi_2 + \lambda_3 \phi_3 + \lambda_4 \phi_4$$

Compute derivatives:

\begin{align} F_x &= 0 & F_{\dot x} &= \lambda_1 \\ F_y &= 0 & F_{\dot x} &= \lambda_2 \\ F_u &= -\lambda_1 & F_{\dot x} &= \lambda_3 \\ F_v &= -\lambda_2 & F_{\dot x} &= \lambda_4 \\ && F_\theta &= \lambda_3 a\sin\theta - \lambda_4 a\cos\theta \\ \end{align}

And set up the Euler-Lagrange equations:

\begin{align} \dot\lambda_1 &= 0 \\ \dot\lambda_2 &= 0 \\ \dot\lambda_3 &= -\lambda_1 \\ \dot\lambda_4 &= -\lambda_2 \\ \end{align} \\ \lambda_3 a\sin\theta - \lambda_4 a\cos\theta = 0

From these we can see that:

\begin{align} \lambda_1 &= c_1 \\ \lambda_2 &= c_2 \\ \lambda_3 &= c_3 - c_1 t \\ \lambda_4 &= c_4 - c_2 t \\ \end{align} \\ \theta = \arctan\left(\frac{c_4 - c_2 t}{c_3 - c_1 t}\right)

We can show that with the substitutions $$\theta_0=\arctan\left(-\frac{c_1}{c_2}\right)$$, $$t_0 = \frac{c_1 c_3 + c_2 c_4}{c_1^2 + c_2^2}$$, and $$\tau=\frac{c_1 c_4-c_2 c_3}{c_1^2+c_2^2}$$, we have:

$$\theta = \theta_0 + \arctan\left(\frac{t-t_0}{\tau}\right)$$

In fact, if we look at the plots of $$\theta$$ for each leg of the discrete solution, we see that this is indeed the case:

Now we do get a closed-form expression for $$x$$, $$y$$, $$u$$, and $$v$$ by integrating the equations of motion with this arctan control input, however they are complicated enough that finding the optimal parameters will require numerical solution techniques.

I find that the optimal parameters for the two legs are:

\begin{align} t_1 &= 498.925 & t_2 &= 6175.367 \\ \theta_{0,1} &= 0.543648 & \theta_{0,2} &= 3.238171 \\ t_{0,1} &= 256.232 & t_{0,2} &= 3018.254 \\ \tau_1 &= 113.423 & \tau_2 &= 29.827 \\ \end{align}

And here are some plots of the resulting trajectory. Note the second image is scaled down significantly in the vertical, since the course is 100 times longer in y than in x.

• This is amazing. You knocked 20 seconds off my best time (and about 24 from your previous attempt). The differences are pretty interesting. My velocity loop (vx against vy over time) looks like a straightline approximation to yours. – Dr Xorile May 4 at 2:27