Let us try Ross Millikan's suggestion, cf. his comment above: Let us minimize the (higher-order) functional
$$\tag{1} S~:=~ \frac{1}{2}\int_{t_i}^{t_f} \!dt~ a^2, \qquad a~\equiv~\dot{v},\qquad v~\equiv~\dot{x}, $$
with fixed endpoints
$$\tag{2} x(t_i)~=~x_i , \qquad v(t_i)~=~v_i , \qquad x(t_f)~=~x_f , \qquad v(t_f)~=~v_f , $$
and with fixed initial and final times, $t_i$ and $t_f$. Let $j\equiv\dot{a}$ denote the jerk. An infinitesimal variation of the functional $S$ yields
$$ \delta S~=~\int_{t_i}^{t_f} \! dt~ a \delta a~\stackrel{\text{int. by parts}}{=}~ \int_{t_i}^{t_f} \!dt~\frac{dj}{dt} \delta x +\left[ a \underbrace{\delta v}_{=0} - j \underbrace{\delta x}_{=0}\right]^{t_f}_{t_i} $$ $$\tag{3} ~\stackrel{(2)}{=}~\int_{t_i}^{t_f} \!dt~\frac{dj}{dt} \delta x. $$
The higher-order Euler-Lagrange equation becomes
$$\tag{4} \frac{dj}{dt}~=~0, $$
cf. the fundamental lemma of calculus of variations. In other words, the jerk
$$\tag{5} j~=~A$$
should be a constant of motion; the acceleration
$$\tag{6} a(t)~=~At+B$$
should be an affine function in $t$; the velocity
$$\tag{7} v(t)~=~\frac{A}{2} t^2+ Bt +C $$
is a quadratic polynomial in $t$; while the position
$$\tag{8} x(t)~=~\frac{A}{6} t^3+ \frac{B}{2}t^2+Ct +D$$
is a cubic polynomial in $t$. Thus we arrive at the (cubic) spline fitting problem mentioned in Floris' answer: Four linear equations (2) with four unknowns $(A,B,C,D)$. To solve them, notice the similarity between the two formulas
$$ \tag{9} \frac{v_i+v_f}{2}~\stackrel{(2)+(7)}{=}~\frac{A}{4}(t_i^2+t_f^2)+B\frac{t_i+t_f}{2}+C $$
and
$$ \tag{10} \frac{\Delta x}{\Delta t}~\stackrel{(2)+(8)}{=}~\frac{A}{6}(t_i^2+t_it_f+t_f^2)+B\frac{t_i+t_f}{2}+C. $$
The difference of (9) and (10) depends only on one unknown $A$! Here we have defined
$$\tag{11} \Delta t ~:=~ t_f-t_i, \qquad \Delta x ~:=~ x_f-x_i\quad\text{and} \quad\Delta v ~:=~ v_f-v_i.$$
The jerk and average acceleration become
$$ \tag{12} j~\stackrel{(5)}{=}~A~\stackrel{(9)+(10)}{=}~6\frac{v_i+v_f}{(\Delta t)^2}-12\frac{\Delta x}{(\Delta t)^3}$$
and
$$ \tag{13} \langle a\rangle ~:=~\frac{\int_{t_i}^{t_f} \! dt~ a}{\Delta t} ~\stackrel{(6)}{=}~A\frac{t_i+t_f}{2}+B~\stackrel{(2)+(7)}{=}~\frac{\Delta v}{\Delta t}, $$
respectively. Similarly, the remaining unknowns $B$, $C$, and $D$ can now easily be determined uniquely in terms of the boundary conditions (2).
Best Answer
Sticking to one dimension for simplicity, at a constant acceleration, $a$, the distance travelled in a time $t$ is simply:
$$ s = v_0 t + \frac{1}{2}at^2 $$
So the average velocity, $v_{av} = s/t$, is:
$$ v_{av} = v_0 + \frac{1}{2}at $$
But acceleration $\times$ time is just the change in velocity i.e. $at = v - v_0$ so:
$$ v_{av} = v_0 + \frac{1}{2}(v - v_0) = \frac{v_0 + v}{2} $$