Fastest curve from $p_0$ to $p_1$

Disclaimer: this is only a partial answer to the question as formulated, mainly because it doesn't correctly take into account the condition that the end-point $\vec{p_1}$ is prescribed as well. Please read the bottom lines. Furthermore it is assumed that the problem is two-dimensional.

It is well known that a straight line is the shortest path between two points, but it can only be realized with constant velcities, e.g. as given by: \begin{eqnarray*} x(t) = v_x.t + s_x \\ y(t) = v_y.t + s_y \end{eqnarray*} $\vec{s} = (s_x,s_y)$ is an initial position, $\vec{v} = (v_x,v_y)$ is a (constant) velocity, $t$ is time.
Replace the linear relationship by a quadratic one. It's somewhat reasonable to suppose that it will give us the most efficient solution with non-constant velocities, as required by the OP. \begin{eqnarray*} x(t) = \frac{1}{2} a_x.t^2 + v_x.t + s_x \\ y(t) = \frac{1}{2} a_y.t^2 + v_y.t + s_y \end{eqnarray*} Here $\vec{a} = (a_x,a_y)$ is the constant acceleration. The resulting curve is a quadratic spline.
The begin- and end-points are given. It's easy to see, with $T$ as the end-time: $$ \vec{p_0} = \left[ \begin{array}{c} s_x \\ s_y \end{array} \right] $$ $$ \vec{p_1} = \left[ \begin{array}{c} \frac{1}{2} a_x.T^2 + v_x.T + s_x \\ \frac{1}{2} a_y.T^2 + v_y.T + s_y \end{array} \right] $$ The begin- and end-velocities are given. It's easy to see, with $T$ as the end-time: $$ \vec{v_0} = \vec{v}(0) = \left[ \begin{array}{c} v_x \\ v_y \end{array} \right] $$ $$ \vec{v_1} = \vec{v}(T) = \left[ \begin{array}{c} a_x.T + v_x \\ a_y.T + v_y \end{array} \right] $$ From the last two equations we find: $$ \frac{v_x(T) - v_x(0)}{a_x} = T \quad ; \quad \frac{v_y(T) - v_y(0)}{a_y} = T $$ And we know that $\sqrt{a_x^2 + a_y^2} = 1$. So let's define a constant and normed acceleration: $$ \left[ \begin{array}{c} a_x \\ a_y \end{array} \right] = \left[ \begin{array}{c} v_x(T) - v_x(0) \\ v_y(T) - v_y(0) \end{array} \right] / \sqrt{(v_x(T) - v_x(0))^2 + (v_y(T) - v_y(0))^2} $$ Or $$\vec{a} = \frac{\vec{v_1} - \vec{v_0}}{\left|\vec{v_1} - \vec{v_0}\right|} $$ Everything is defined now. For example the time $T$ is easy to calculate from: $$ T = \frac{v_x(T) - v_x(0)}{a_x} = \sqrt{(v_x(T) - v_x(0))^2 + (v_y(T) - v_y(0))^2} = \left|\vec{v_1} - \vec{v_0}\right| $$ But I'm a bit worried about the "normed acceleration". Let's assume instead that: $$ \vec{a} = \frac{\vec{v_1} - \vec{v_0}}{\left|\vec{v_1} - \vec{v_0}\right|} \times A $$ Where $A$ (a scalar) has the dimension of an acceleration. Then: $$ T = \frac{\left|\vec{v_1} - \vec{v_0}\right|}{A} $$ has the dimension of time. Which sounds better to me, as a physicist :-)

As has been said above, this is a partial answer to the question. The reason is that we can calculate the end-position now, instead of imposing it: $$ \vec{p_1} = \vec{p_0} + \frac{\left|\vec{v_1} - \vec{v_0}\right|}{A} \frac{\vec{v_1}+\vec{v_0}}{2} $$ One reason that I've undeleted this post again is that there seems to be no other possibility than the quadratic spline if one wants to maintain a constant optimal acceleration from the beginning to the end. More information about quadratic splines is found in: Splines (PDF).