I don't like the notation you're using, so I will use a different one. Consider the ellipses defined below
$\mathcal{E}_1 := \displaystyle\left\{ (x_1, y_1) \in \mathbb{R}^2 : \left(\frac{ x_1 - x_{10}}{a_1}\right)^2 + \left(\frac{ y_1 - y_{10}}{b_1}\right)^2 = 1 \right\}$
$\mathcal{E}_2 := \displaystyle\left\{ (x_2, y_2) \in \mathbb{R}^2 : \left(\frac{ x_2 - x_{20}}{a_2}\right)^2 + \left(\frac{ y_2 - y_{20}}{b_2}\right)^2 = 1 \right\}$
Note that we have 4 unknowns and 2 constraints. Hence, we have 2 degrees of freedom. It would be neater to solve the optimization problem in 2 variables only. But, how?
I introduce functions $\gamma_1, \gamma_2 : [0, 2 \pi] \to \mathbb{R}^2$ defined by
$\gamma_1 (\theta_1) = \left[\begin{array}{c} x_{10} + a_1 \cos(\theta_1)\\ y_{10} + b_1 \sin(\theta_1)\end{array}\right]$
$\gamma_2 (\theta_2) = \left[\begin{array}{c} x_{20} + a_2 \cos(\theta_2)\\ y_{20} + b_2 \sin(\theta_2)\end{array}\right]$
In other words, I have parametrized the ellipses. Since you want to minimize the distance between the elipses, we have the following optimization problem
$\displaystyle\min_{(\theta_1, \theta_2)} \| \gamma_1 (\theta_1) - \gamma_2 (\theta_2)\|_2^2 \quad{} \text{subject to} \quad{} (\theta_1, \theta_2) \in [0, 2 \pi]^2$
which you can solve using partial derivatives. No need for Lagrange multipliers.