De momento puedo ofrecer una solución analítica, bastante engorrosa pero manejable. Espero que alguien pueda dar una solución más sencilla.
Para reducir el número de parámetros, podemos establecer un sistema de coordenadas tal que $\mathbf{p}_1=(\alpha,0)$ , $\mathbf{p}_2=(0,\beta)$ , $\mathbf{v}_1=(0,1)$ y definir entonces $m=v_{2y}/v_{2x}$ .
Empecemos con la ecuación genérica de una sección cónica, en la que podemos poner a $1$ el coeficiente de $x^2$ porque se trata de una elipse: $$ x^2+By^2+Cxy+Dx+Ey+F=0. $$
Tenemos cuatro condiciones: la cónica pasa por $\mathbf{p}_1$ , $\mathbf{p}_2$ y es tangente a $\mathbf{v}_1$ , $\mathbf{v}_2$ . Estas condiciones se traducen en cuatro ecuaciones: $$ \begin{align} \cases{ \alpha^2+\alpha D+F=0\\ \beta^2B+\beta E+F=0\\ 2m\beta B+\beta C+D+mE=0\\ \alpha C+E=0 } \end{align} $$ A partir de ahí, se pueden encontrar expresiones para $B$ , $C$ , $D$ , $E$ en función de $F$ . Conectándolos a la fórmulas para los semiejes $a$ y $b$ podemos entonces calcular la excentricidad $\epsilon=\sqrt{1-b^2/a^2}$ . Cuanto menor sea la excentricidad, más redonda será la elipse, por lo que debemos encontrar el mínimo de $\epsilon$ en función de $F$ .
Computé $d\epsilon/dF$ con Mathematica y descubrí que desaparece para $$ F=-\frac{\alpha ^2 \beta \left(\alpha ^2 \beta +\beta ^3-2 \alpha ^2 \beta m^2+2 \alpha ^3 m\right)}{\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}. $$
Una vez $F$ se pueden calcular los demás coeficientes y hallar la ecuación de la elipse "más redonda": $$ \begin{align} \cases{ B=\frac{\displaystyle\alpha ^2 \left(\beta ^2+\alpha ^2 \left(2 m^2+1\right)+2 \alpha \beta m\right)}{\displaystyle\beta ^2\left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}\\ \\ C=\frac{\displaystyle2 \alpha ^2 m \left(-\alpha ^2+\beta ^2+2 \alpha \beta m\right)}{\displaystyle\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}\\ \\ D=-\frac{\displaystyle2 \alpha ^2 m \left(2 \beta ^3+m \left(\alpha ^3+3 \alpha \beta ^2\right)\right)}{\displaystyle\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}\\ \\ E=-\frac{\displaystyle2 \alpha ^3 m \left(-\alpha ^2+\beta ^2+2 \alpha \beta m\right)}{\displaystyle\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)} } \end{align} $$
A comprobar con GeoGebra confirma que éste es efectivamente el valor de $F$ dando una excentricidad mínima.
EDITAR.
Siguiendo la sugerencia de Rahul, podemos obtener un resultado mucho más sencillo. Elige un sistema de coordenadas tal que $\mathbf{p}_1=(\alpha,0)$ , $\mathbf{p}_2=(-\alpha,0)$ , y definimos $m=v_{1y}/v_{1x}$ , $n=v_{2y}/v_{2x}$ .
Con esas opciones, la ecuación de la elipse se puede escribir como: $$ x^2+By^2-{m+n\over mn}xy+\alpha{m-n\over mn}y-\alpha^2=0. $$ La excentricidad $\epsilon$ es entonces función de $B$ y $d\epsilon/dB=0$ para $$ B=1+\frac{(m+n)^2}{2 m^2 n^2}. $$ En resumen, la ecuación de la elipse más redonda resulta ser: $$ x^2+\left(1+\frac{(m+n)^2}{2 m^2 n^2}\right)y^2-{m+n\over mn}xy+\alpha{m-n\over mn}y-\alpha^2=0. $$
Añadido por Rahul:
Las cosas se simplifican aún más si utilizamos la pendiente del normales , $\mu=-1/m$ y $\nu=-1/n$ : $$ x^2+\left(1+\tfrac12(\mu+\nu)^2\right)y^2+(\mu+\nu)xy+\alpha(\mu-\nu)y-\alpha^2=0. $$
