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).
