Solution 1
Write the line we are looking for as: $ax+by+cz = 0$ (in the projective plane). The purpose is to determine the vectors of coefficients $(a,b,c)$. Since the line can be seen as a point in the dual plane, we shall consider the following normalization: $a^2+b^2+c^2 = 1$.
Then you can formulate the problem as follows: Minimize $d(l,p)^2 + d(l,q)^2$ under the constraints: $av_x + bv_y + cv_z = 0$ and $a^2+b^2+c^2 = 1$.
Now: $d(l,p)^2 = (ap_x +bp_y + c)^2/(a^2 + b^2)$ and $d(l,q)^2 = (aq_x +bq_y + c)^2/(a^2 + b^2)$
Form the lagrangian of optimization problem: $L(a,b,c,u,v) = \frac{(ap_x +bp_y + c)^2}{a^2 + b^2} + \frac{(ap_x +bp_y + c)^2}{a^2 + b^2} +u(av_x + bv_y + cv_z) + v(a^2+b^2+c^2-1)$
You need to compute the derivative with respect to $a$, $b$, $c$, $u$, $v$, then multiply by denominator and you get a polynomial system to solve.
Note that the Lagrange multipliers $(u,v)$ are auxiliary unknowns. The second solution avoids introducing these unknowns.
Solution 2:
Take 2 distinct lines $l_1$ and $l_2$ passing through $v$. Then the line we are looking for can be written $l=\alpha l_1 + \beta l_2$.
Express the function to be minimized as a function of $\alpha$ and $\beta$. Take the derivative and you will take a polynomial in $\alpha, \beta$ (it should be homogeneous).