This is just a try to solve:
The parametric equations of an ellipse at the origin is
$$\begin{align*}\vec{c_0(t)} &= \begin{bmatrix}x(t)\\y(t)\end{bmatrix}\\&= \begin{bmatrix} a ~\cos(t)\\b ~\sin(t)\end{bmatrix}\end{align*} $$
Rotate the ellipse $\theta_z$ about the $z$ axis and get this
$$\begin{align*}\vec{c_1(t)} &=R_{\theta_z}~\vec{c_0}(t)\\&= \begin{bmatrix}\cos(\theta_z)&&-\sin(\theta_z)\\\sin(\theta_z)&&\cos(\theta_z)\end{bmatrix}\vec{c_0}(t)\\&= \begin{bmatrix}\cos(\theta_z)&&-\sin(\theta_z)\\\sin(\theta_z)&&\cos(\theta_z)\end{bmatrix}\begin{bmatrix}x(t)\\y(t)\end{bmatrix}\\&= \begin{bmatrix}x(t)~\cos(\theta_z)- y(t)~\sin(\theta_z) \\ x(t)~\sin(\theta_z)+ y(t)~\cos(\theta_z)\end{bmatrix}\\&= \begin{bmatrix}a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z) \\ a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)\end{bmatrix} \end{align*}$$
Translate the ellipse by a vector $\vec{v}_{T}=(\Delta x ~ ~ \Delta y)^T $ to some other place other than the origin
$$\begin{align*}\vec{c_2}(t) &= \vec{c_1}(t) + \vec{v}_T \\&=\begin{bmatrix}a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z) \\ a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)\end{bmatrix}+\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix}\\&= \begin{bmatrix}a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z)+\Delta x \\ a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)+\Delta y\end{bmatrix} \end{align*}$$
We could have also first translated and then rotated to get
$$\begin{align*}\vec{c_{2'}}(t)&=R_{\theta_z}~(\vec{c_0}(t)+ \vec{v}_T)\\&=\begin{bmatrix}\cos(\theta_z)&&-\sin(\theta_z)\\\sin(\theta_z)&&\cos(\theta_z)\end{bmatrix} \left( \begin{bmatrix} a ~\cos(t)\\b ~\sin(t)\end{bmatrix}+\begin{bmatrix}\Delta x\\\Delta y \end{bmatrix}\right)\\&= \begin{bmatrix}\cos(\theta_z)&&-\sin(\theta_z)\\\sin(\theta_z)&&\cos(\theta_z)\end{bmatrix} \begin{bmatrix} a ~\cos(t)+\Delta x\\b ~\sin(t)+\Delta y\end{bmatrix}\\&=\begin{bmatrix} a~\cos(\theta_z) ~\cos(t)-b~\sin(\theta_z)~\sin(t)+\Delta x ~\cos(\theta_z)-\Delta y ~\sin(\theta_z)\\a~\sin(\theta_z) ~\cos(t)+b~\cos(\theta_z)~\sin(t)+\Delta x~ \sin(\theta_z)+\Delta y ~\cos(\theta_z)\end{bmatrix}\end{align*}$$
From this expression we can calculate the position of a point on the circumference of an ellipse for a parameter value $t$. Given a point $\vec{P} = (x_0 ~ ~ y_0)^T$ in the plane, on the ellipse, inside oder outside the ellipse, we can calculate the distance $d$ between this point $\vec{P}$ and "a general point" on the ellipse like this:
$$\begin{align*}d(\vec{c_2}(t),\vec{P})&=f(t)\\&=\left\|\vec{c_2}(t)-\vec{P}\right\|\\&=\left\|\begin{bmatrix}a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z)+\Delta x \\ a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)+\Delta y\end{bmatrix}-\begin{bmatrix} x_0\\y_0\end{bmatrix}\right\|\\&=\left\|\begin{bmatrix}a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z)+\Delta x-x_0 \\ a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)+\Delta y-y_0\end{bmatrix}\right\|\\&=\sqrt{(a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z)+\Delta x-x_0 )^2+(a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)+\Delta y-y_0)^2}\end{align*}$$
To find out the shortest distance $d$ between $\vec{P}$ and the ellipse we will have to find the value of $t$ that minimizes the value of the following expression (under the square root above):
$$(a ~\cos(\theta_z)~\cos(t)- b ~\sin(\theta_z)~\sin(t)+\Delta x-x_0 )^2+(a~\sin(\theta_z) ~\cos(t)+ b ~\cos(\theta_z)~\sin(t)+\Delta y-y_0)^2$$
This isn't a very easy task, but if we find such a $t$ then we are done, because we can plug it in $\vec{c_{2}}(t)$ or $\vec{c_{2'}}(t)$ to get our point on the ellipse, which has the shortest distance from $\vec{P}$. Maybe this is not the best approach to find out $t$, look what I found here. Maybe we should calculate the tangent:
$$\vec{c'_2}(t) = \begin{bmatrix}-a ~\sin(t)~\cos(\theta_z)- b ~\cos(t)~\sin(\theta_z)\\ -a ~\sin(t)~\sin(\theta_z)+ b ~\cos(t)~\cos(\theta_z)\end{bmatrix} $$
and calculate the vector $\vec{l}(t)$ from the point $\vec{P}$ to a general point on the ellipse, we have it already :
$$\vec{l}(t) = \vec{c_2}(t)-\vec{P} = \begin{bmatrix}a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z)+\Delta x-x_0 \\ a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)+\Delta y-y_0\end{bmatrix} $$
Our $t$ can also be found with the condition that the vectors $\vec{l}(t)$ and $\vec{c'_2}(t)$ have to be perpendicular to each other, their dot product has to vanisch. This means: the $t$, that makes the following expression is equal to $0$ is the good old $t$ that we've been looking for:
$$(a ~\cos(t)~\cos(\theta_z)- b ~\sin(t)~\sin(\theta_z)+\Delta x-x_0)(-a ~\sin(t)~\cos(\theta_z)- b ~\cos(t)~\sin(\theta_z))+(a ~\cos(t)~\sin(\theta_z)+ b ~\sin(t)~\cos(\theta_z)+\Delta y-y_0)(-a ~\sin(t)~\sin(\theta_z)+ b ~\cos(t)~\cos(\theta_z))=0 $$
This is just exactly the same condition (Yes, it is!) we get, if we differntiate the above expression (under the square root) and set the result $=0$ to find the min. value! The question remains the same: For which $t$ do we get the minimal value under the square root or for which $t$ does this expression above (=the derivative of the expression under the square root) become $=0$?
Assume $a_1 = a\cos(\theta_z), b_1 = b\sin(\theta_z), a_2 = a\sin(\theta_z), b_2 = b\cos(\theta_z), c_1 = \Delta x - x_0$ and $c_2 = \Delta y - y_0$. Again assume $q_1 =-a_1~c_1-a_2~c_2 ,$ $q_2= -b_1~c_1+b_2~c_2,$ $q_3=(-a_1^2+b_1^2-a_2^2+b_2^2)/2$ and $q_4=-a_1~b_1+a_2~b_2$ and use the geometric identities $\cos^2(t)-\sin^2(t)=\cos(2t)$ and $\sin(t)\cos(t)=\sin(2t)/2$, then the equation becomes $$ q_1 \cdot \sin(t)+q_2 \cdot \cos(t)+ q_3 \cdot \sin(2t) + q_4 \cdot \cos(2t) =0$$
Finding a closed formula is cumbersome, brother! After some research on the topic, this equation can be solved numerically with the newton method to get the desired $t=t^*$. Then the distance will be $d_{min} = \left\|\vec{c}_2(t^*)-\vec{P}\right\|$.
Also have a look at this answer (from Doctor Rob, The Math Forum)
There may be a general formula for this, but if so, it is so
complicated that no one would ever write it out explicitly. I would
find the distance in the following way.
The point Q(x_0,y_0) on the ellipse whose distance from the given point
P(X,Y) is least at a point such that the tangent to the ellipse at Q
is perpendicular to the line PQ. The slope of the tangent at Q can be
found by implicitly differentiating the equation of the ellipse and
solving for dy/dx, then substituting in x0 and y0. Then the line
perpendicular to this tangent has slope -1/(dy/dx), and it passes
through Q, so the point-slope form will give you its equation in terms
of A, B, C, D, E, F, x0, and y0. Since the point P lies on this line,
that gives you an equation that x0 and y0 must satisfy, in addition to
the equation of the ellipse. By eliminating y0 from these two
equations, you get a quartic equation in x0. This may have zero, two,
or four real solutions. Each one will give a corresponding value of
y0.
Then for these solutions, the distance PQ should be computed, and the
smallest of them chosen as the answer. For an ellipse, the largest one
will give you the point farthest from P.
In theory, you can solve the quartic in terms of radicals, but the
result of trying to do that is generally an unmanageably complicated
formula. Once you have the quartic, if you cannot factor it into
linear and quadratic factors, you might as well solve it numerically.
Note that this solution does not depend on the curve being an ellipse,
just a quadratic equation in x and y. It will give you the answer for
parabolas, hyperbolas, and degenerate conics, too.
As an example, take the point P(3,8) and the ellipse
x^2 + x*y + y^2 + 2*x + 4*y - 9 = 0.
Then
2*x + x*(dy/dx) + y + 2*y*(dy/dx) + 2 + 4*(dy/dx) = 0,
dy/dx = -(2*x+y+2)/(x+2*y+4),
and the slope of the perpendicular to the tangent line at Q(x0,y0) is
m = (x0+2*y0+4)/(2*x0+y0+2),
and that line is
y - y0 = [(x0+2*y0+4)/(2*x0+y0+2)]*(x-x0).
Since P(3,8) is on that line, we get the equation
8 - y0 = [(x0+2*y0+4)/(2*x0+y0+2)]*(3-x0),
4 + 17*x0 + x0^2 = y0^2,
and since Q is on the curve, we get the equation
x0^2 + x0*y0 + y0^2 + 2*x0 + 4*y0 - 9 = 0.
Eliminating y0 from these two equations, we end up with the quartic
equation
3*x0^4 + 51*x0^3 + 185*x0^2 - 494*x0 - 39 = 0.
This polynomial doesn't factor, so we solve numerically. There
are two real roots, x0 = -0.0767858556 and x0 = 1.7998308063,
approximately, and the corresponding values of y are
y0 = 1.6433309230 and y0 = -6.1511392960, approximately. The first
point has distance 7.0621422356 and the second 14.2019417499 from P,
so the minimum distance from P to the ellipse is 7.0621422356.
Another approach is to use a Lagrange Multiplier. Set
f(x,y,z) = (x-X)^2 + (y-Y)^2 + z*(A*x^2+B*x*y+C*y^2+D*x+E*y+F).
Then set the partials of f with respect to x, y, and z equal to zero,
and solve for x, y, and z. That will give you x = x0 and y = y0. The
resulting equations will be the same, because eliminating z from the
first two equations will give you the equation of the line PQ we got
above, and the last equation will be the equation of the ellipse.