Given the following real input parameters:
$$
\left(x_i,\,y_i\right),
\quad \quad \quad
\left(x_f,\,y_f\right),
\quad \quad \quad
\left(a_x,\,a_y\right),
\quad \quad \quad
v_i > 0
$$
and define the differences:
$$
\Delta x := x_f - x_i\,,
\quad \quad \quad
\Delta y := y_f - y_i
$$
we need to solve the system of equations:
$$
\begin{cases}
\Delta x = v_i\cos\theta\,\Delta t + a_x\,\Delta t^2/2 \\
\Delta y = v_i\sin\theta\,\Delta t + a_y\,\Delta t^2/2 \\
\end{cases}
$$
where solutions $\Delta t > 0$ and $0 \le \theta < 2\pi$ are of interest.
So, rewriting the system of equations:
$$
\begin{cases}
\frac{2\,\Delta x - a_x\,\Delta t^2}{2\,v_i\,\Delta t} = \cos\theta \\
\frac{2\,\Delta y - a_y\,\Delta t^2}{2\,v_i\,\Delta t} = \sin\theta \\
\end{cases}
$$
it follows that:
$$
\left(\frac{2\,\Delta x - a_x\,\Delta t^2}{2\,v_i\,\Delta t}\right)^2 + \left(\frac{2\,\Delta y - a_y\,\Delta t^2}{2\,v_i\,\Delta t}\right)^2 = 1
$$
that is, by reworking:
$$
\left(a_x^2 + a_y^2\right)\Delta t^4 - 4\left(v_i^2 + a_x\,\Delta x + a_y\,\Delta y\right)\Delta t^2 + 4\left(\Delta x^2 + \Delta y^2\right) = 0
$$
where, defined the parameters:
$$
a := a_x^2 + a_y^2\,,
\quad \quad \quad
b := - 4\left(v_i^2 + a_x\,\Delta x + a_y\,\Delta y\right),
\quad \quad \quad
c := 4\left(\Delta x^2 + \Delta y^2\right)
$$
it's clear that there are solutions $\Delta t > 0$ if and only if:
$$
a + c \ne 0
\quad \quad \quad
\text{and}
\quad \quad \quad
b < 0
\quad \quad \quad
\text{and}
\quad \quad \quad
b^2 - 4\,a\,c \ge 0
$$
and in this case it's only one if:
$a = 0$ from which $\Delta t = \sqrt{\frac{c}{-b}}$;
$c = 0$ from which $\Delta t = \sqrt{\frac{-b}{a}}$;
$b^2 - 4\,a\,c = 0$ from which $\Delta t = \sqrt{\frac{-b}{2\,a}}$;
otherwise are two:
$$
\Delta t_1 = \sqrt{\frac{- b - \sqrt{b^2 - 4\,a\,c}}{2\,a}} < \Delta t_2 = \sqrt{\frac{- b + \sqrt{b^2 - 4\,a\,c}}{2\,a}}\,.
$$
In particular, if there are solutions $\Delta t > 0$, the initial velocity vector is trivially determined:
$$
v_{x,i} = v_i\,\cos\theta = \frac{\Delta x}{\Delta t} - \frac{a_x}{2}\,\Delta t\,,
\quad \quad \quad
v_{y,i} = v_i\,\sin\theta = \frac{\Delta y}{\Delta t} - \frac{a_y}{2}\,\Delta t
$$
and consequently also the angle measured with respect to the positive semi-axis of the abscissa:
$$
\theta =
\begin{cases}
\pi + \arctan\left(v_{y,i}/v_{x,i}\right) & \text{if} \; v_{x,i} < 0 \\
3\pi/2 & \text{if} \; v_{x,i} = 0 \, \land \, v_{y,i} < 0 \\
\pi/2 & \text{if} \; v_{x,i} = 0 \, \land \, v_{y,i} > 0 \\
2\pi + \arctan\left(v_{y,i}/v_{x,i}\right) & \text{if} \; v_{x,i} > 0 \, \land v_{y,i} < 0 \\
\arctan\left(v_{y,i}/v_{x,i}\right) & \text{if} \; v_{x,i} > 0 \, \land v_{y,i} \ge 0 \\
\end{cases}
$$
which, if desired, is equivalent to writing:
$$
\theta =
\begin{cases}
2\pi + \text{arctan2}\left(v_{y,i},\,v_{x,\,i}\right) & \text{if} \; v_{y,i} < 0 \\
\text{arctan2}\left(v_{y,i},\,v_{x,\,i}\right) & \text{if} \; v_{y,i} \ge 0 \\
\end{cases}
$$
where $\text{arctan2}(y,\,x)$ associates an angle in the interval $(-\pi,\,\pi]$ for each point $(x,\,y) \ne (0,\,0)$.
Since this is a uniformly accelerated motion, both in the horizontal and vertical directions, the parabolic trajectories will generally be rototranslated in the plane. For example, assuming:
$$
\left(x_i,\,y_i\right) = (1,\,1),
\quad \quad \quad
\left(x_f,\,y_f\right) = (3,\,2),
\quad \quad \quad
\left(a_x,\,a_y\right) = (1,\,-1),
\quad \quad \quad
v_i = 2
$$
the two trajectories solution of the problem are:
where blue is for $\Delta t_1 \approx 1.06$ and $\theta_1 \approx 47.4°$, while red is for $\Delta t_2 \approx 2.98$ and $\theta_2 \approx 114°$.
Finally, it isn't difficult to prove that the area enclosed by the two parabolic arcs is equal to:
$$
\mathcal{A} = \frac{\left|a_x\,\Delta y - a_y\,\Delta x\right|\left(\Delta t_2^2 - \Delta t_1^2\right)}{12} = \frac{\sqrt{15}}{2} \approx 1.94\,.
$$
I hope it's clear enough, good luck! ^_^
Best Answer
(i) Since the unknown vector $\vec v_p $ satisfies an equation of the form $\vec v_p t = \sum_{k=0}^n \vec c_k t^k$ in which each $\vec c_k$ is a known constant vector, you can write $\vec v_p = t^{-1} \sum_{k=0}^n \vec c_k t^k$. Denote the right-hand side as $\vec R(t)$ and treat it as a known parametrized particle path in 3-space.
(ii) Your goal is to find the value of $t>0$ that minimizes the squared length $L^2 = ||\vec R(t)||^2$. If you visualize the trajectory traced by $\vec R(t)$ as a smooth space curve, then it is clear that $L$ is the radius of a sphere centered at the origin that can be drawn tangent to this curve. Imagine temporarily that the space curve is parametrized by its arclength. Then at the critical time where the curve contacts this critical sphere, the velocity vector $\vec T= \vec R'$ is perpendicular to the radius vector $\vec R$. This scalar equation $\vec R \cdot \vec R'=0$ can be solved for the unknown scalar $t$. Thus you now know the critical value of $\vec R$ and of $L$ which is a candidate for the absolute minimum.
(iii) Thus we have reduced the problem to solving a single scalar equation, and that equation (after a bit of massaging) can be re-written as a polynomial in $t$. Newton's method is the standard numerical algorithm for solving such polynomial equations when the degree exceeds two. Or you may have access to a built-in root-solver.