This approach is using homogeneous coordinates. Write a point of the plane as $p=(x,y,1)^T$. Then you can come up with symmetric matrices $A$ and $B$ to represent the circle and the ellipse, respectively, in such a way that $p^T\cdot A\cdot p=0$ means the point is on the circle and similarly for $B$. Note that the matrix $A$ will contain the unknown radius as a parameter, or more specifically will contain its square, $r^2$. I'm giving the concrete matrices for your use case a bit further down.
Now consider a linear combination $C=A+\lambda B$. It represents the pencil of conics, i.e. all the conic sections that share the four points of intersection with the two conics $A$ and $B$. For your intuition I suggest you imagine a circle radius where the circle actually has four real intersections with the ellipse, otherwise some of the common points will have complex coordinates which makes intuition harder. So depending on how you choose $\lambda$ you get all the conics through these four common points.
In this pencil of conics there are usually three degenerate cases, where the conic becomes a pair of lines. Computing these can be an important step in intersecting two conics. If you think of your four points of intersection as $a, b, c, d$ then you get the line pairs $(ab,cd), (ac,bd), (ad,bc)$. The way to find these special cases is by requiring $\det C\overset!=0$ which leads to a cubic equation in $\lambda$.
Now you are interested in the special case where the circle touches the ellipse. Which means that two of the points of intersection actually coincide. Which means that two of the degenerate conics are the same. Which means that one of the solutions of your cubic equation in $\lambda$ has multiplicity two. This your can characterized via the discriminant of that cubic polynomial being zero. So that gives you a single equation in $r$. You can find the point of contact in a later step based on that radius.
The above is generic, applicable to any pair of conics with some parametrization, which I assume may be of value to other readers as well. Below are concrete computations for your specific case. Your circle is $x^2+y^2-r^2=0$ so you get
$$A=\begin{pmatrix}1&0&0\\0&1&0\\0&0&r^2\end{pmatrix}$$
Your ellipse is $b^2(x^2-2hx+h^2)+a^2(y^2-2ky+k^2)-a^2b^2=0$ so you get
$$B=\begin{pmatrix}
b^2 & 0 & -b^2h\\
0 & a^2 & -a^2k\\
-b^2h & -a^2k & b^2h^2+a^2k^2-a^2b^2
\end{pmatrix}$$
My computer algebra system tells me that
$$\det(C)=-a^4b^4\lambda^3 + (r^2a^2b^2 - a^4b^2 - a^2b^4 + a^2b^2h^2 + a^2b^2k^2)\lambda^2 + (r^2a^2 + r^2b^2 - a^2b^2 + b^2h^2 + a^2k^2)\lambda + r^2$$
The discriminant of this is of degree 4 in $r^2$. I could paste the generic formula here, but I think you're probably better of substituting concrete numbers before computing the discriminant. Note that this will give you 4 distinct solutions where there circle touches the ellipse. Some of them may have the circle and the ellipse intersect in other points beside the one where they touch.
