I'll walk through the derivation of equations capturing the problem. Symbolically, things get pretty ugly (with degree-$14$ polynomials in hundreds of terms), so numerical methods are crucial. But, to show that my approach works, here's a sample configuration:
I'm rearranging and relabeling things a bit.
For a square of side-length $2s$, coordinatize with
$$A=(-s,-s) \qquad B = (s,-s) \qquad C = (s,s) \qquad D = (-s,s) \tag1$$
Let $A'$ and $B'$ be the centers of the circles (offset from $A$ and $B$ instead of OP's $A$ and $D$), with
$$A' := A + (h+k)\,(1,1) \qquad B' := B + (h-k)\,(-1,1) \tag2$$
for some $h$ and $k$. Let $A''$ and $B''$ be the points of tangency of those circles with the ellipse; we'll write:
$$A'' := A' + (h+k)\,(\cos\alpha, \sin\alpha) \qquad B'' := B' + (h-k)\, (-\cos\beta,\sin\beta) \tag3$$
(I've negated $\cos\beta$ in $B''$ for a kind of mirror symmetry with $A''$.)
Finally, let $C'$, $D'$, $T$ be the points of tangency of the square wth the ellipse; with
$$C' := C + (u-v)\,(0,-1) \qquad D' := D + (u+v)\,(0,-1) \qquad T := (t,s) \tag4$$
Now to derive some equations relating the parameters $h$, $k$, $\alpha$, $\beta$, $u$, $v$, $t$ ...
The tangency of the circles requires
$$|A'B'|^2 = (2h)^2 \quad\to\quad s^2+k^2-2hs=0 \quad\to\quad h = \frac{s^2+k^2}{2s} \tag5$$
That was easy! Continuing ...
To get the conic tangent to the square at $C'$, $D'$, $T$, I use a trick of nudging copies of two points to get a total of five, which determine an conic; then I collapse the copies back onto the originals. Algebraically, I define
$$C'' := C' + {0,\varepsilon} \qquad D'' := D' + (0,\varepsilon) \tag6$$
and use this determinant formula shown in this answer for the conic through the five points $C'$, $D'$, $T$, $C''$, $D''$; the result has a factor of $\varepsilon^2$, which I conveniently cancel (because it's non-zero) before then substituting $\varepsilon\to 0$ everywhere else (because it's vanishingly-small), to get
$$\begin{align}
0 &= x^2\,\left(s \left(u^2 + v^2\right) - 2 t u v\right)
\;+\; (y - (s+u))^2 \; s \left(s^2 - t^2\right) \\
&\quad+ 2 x (y - (s+u))\; v \left(s^2 - t^2\right) \\
&\quad- s (s u - t v)^2
\end{align} \tag7$$
Let the tangent vector to the conic at $(x,y)$ be $(\dot{x},\dot{y})$. The differential form of $(9)$ relates these values:
$$\begin{align}
0 &= \dot{x}x\,\left(s \left(u^2 + v^2\right) - 2 t u v\right)
\;+\; \dot{y}(y - (s+u)) \; s \left(s^2 - t^2\right) \\
&\quad+ (\dot{y}x +\dot{x}(y - (s+u)))\; v \left(s^2 - t^2\right)
\end{align} \tag{7'}$$
This must hold when $(x,y)=T$ and $(\dot{x},\dot{y})=(1,0)$, which tells us
$$st \left(u^2 + v^2\right) = u v \left(s^2+t^2\right) \quad\to\quad t = \frac{su}{v},\frac{sv}{u}\tag{8}$$
The solution $t=su/v$ is extraneous (it reduces $(7)$ to a line through $C'$ and $D'$, so we take $t=sv/u$. With this, $(7)$ and $(7')$ become
$$\begin{align}
u^2 x^2 \;+\; s^2 (y - (s+u))^2 \;+\; 2 s v x (y - (s+u)) &\;=\; s^2 (u^2 - v^2) \tag{9} \\[4pt]
\dot{x}\left(u^2 x + s v (y - (s+u))\right)
\;+\; \dot{y}\,s (v x + s (y - (s+u))) &\;=\; 0 \tag{9'}
\end{align}$$
We want these relations to be hold for points $A''$ and $B''$, with corresponding tangent vectors $(-\sin\alpha,\cos\alpha)$ and $(\sin\beta,\cos\beta)$. So, "all we have to do" is substitute (recalling $(5)$) and solve for some chosen four of $k$, $u$, $v$, $\alpha$, $\beta$ in terms of the fifth.
Any way of going about it seems to be a mess. From the first substitutions, we have four polynomials equations (two from $(9)$ and two from $(9')$) that are quartic in $k$, and quadratic in $u$ and $v$, with trigs of $\alpha$ and $\beta$ floating around. Eliminating variables makes things worse and worse.
As for the trigs ... In order that Mathematica can manipulate things as simple polynomials, I used the tangent half-angle substitution
$$a := \tan\frac12\alpha \qquad b := \tan\frac12\beta \tag{10}$$
so that
$$\cos \alpha = \frac{1-a^2}{1+a^2} \qquad \sin\alpha = \frac{2a}{1+a^2} \qquad
\cos\beta = \frac{1-b^2}{1+b^2} \qquad \sin\beta = \frac{2b}{1+b^2} \tag{11}$$
Note that these made anything prettier.
At this point, I resort to numerical methods. The figure above shows a case $s=5$ and $u=-3$, with remaining values taken from the last row of this table of Mathematica's full set of real solutions (and $h$ and $t$ calculated from $(5)$ and $(8)$):
$$\begin{array}{cccc}
k & v & a & b \\
\hline
\phantom{-}44.9302 & \phantom{-}659.359 & -0.975401 & -2.45736 \\
-44.9302 & -659.359 & -2.45736 & -0.975401 \\
\hline
-16.7662 & -120.821 & 0.431426 & -0.883977 \\
\hline
-19.4486 & \phantom{-}49.6271 & -0.767683 & -2.4851 \\
\phantom{-}19.4486 & -49.6271 & -2.4851 & -0.767683 \\
\hline
\phantom{-}25.3198 & \phantom{-}11.4631 & -1.07842 & -2.46041 \\
-25.3198 & -11.4631 & -2.46041 & -1.07842 \\
\hline
-12.7342 & -3.69367 & 0.420587 & -1.33054 \\
\hline
\phantom{-}12.2945 & \phantom{-}0.342764 & -2.35157 & \phantom{-}0.326022 \\
-12.2945 & -0.342764 & \phantom{-}0.326022 & -2.35157 \\
\hline
\phantom{-}0.460327 & \phantom{-}2.46871 & 0.598755 & 1.20824 \\
-0.460327 & -2.46871 & 1.20824 & 0.598755 \\
\hline
\end{array}$$
Clearly, only last two of the above cases has a chance of corresponding to the prescribed configuration of an "inscribed ellipse"; the others put tangent points $C'$ and $D'$ well beyond the square. Generally, it's not even guaranteed that conic $(7)$ is an ellipse. (The condition for ellipse-ness is actually pretty simple. I'll leave that to the reader.) Moreover, when we do have an inscribed ellipse, as desired, the circles may turn out to be internally tangent to it. (I had to go through a few choices of $s$ and $u$ to generate a configuration that looked like OP's original.)
Anyway ... That's about as far as I expect to take this analysis. I hope it's been helpful.
