Optimize distance from origin to point on superellipsoid I've been trying to solve this problem with a couple of methods so far (including Lagrange multipliers) but for each method I end up with an unsolvable equation. 
My first approach (Lagrange multipliers) starts off with optimizing $x^2+y^2+z^2$ constrained to $\left(\left|\dfrac{x}{a}\right|^r+\left|\dfrac{y}{b}\right|^r\right)^{\frac{t}{r}}+\left|\dfrac{z}{c}\right|^t-1 = 0$ where $a$, $b$, $c$, $r$, and $t$ are variables that tweak the shape of the superellipsoid. From there, I go from this: 
$$
\begin{bmatrix}
\dfrac{\partial}{\partial x}\left(x^2 + y^2 + z^2 - \lambda\left(\left(\left|\dfrac{x}{a}\right|^r+\left|\dfrac{y}{b}\right|^r\right)^{\frac{t}{r}}+\left|\dfrac{z}{c}\right|^t-1\right)\right)\\
\dfrac{\partial}{\partial y}\left(x^2 + y^2 + z^2 - \lambda\left(\left(\left|\dfrac{x}{a}\right|^r+\left|\dfrac{y}{b}\right|^r\right)^{\frac{t}{r}}+\left|\dfrac{z}{c}\right|^t-1\right)\right)\\
\dfrac{\partial}{\partial z}\left(x^2 + y^2 + z^2 - \lambda\left(\left(\left|\dfrac{x}{a}\right|^r+\left|\dfrac{y}{b}\right|^r\right)^{\frac{t}{r}}+\left|\dfrac{z}{c}\right|^t-1\right)\right)\\
\dfrac{\partial}{\partial\lambda}\left(x^2 + y^2 + z^2 - \lambda\left(\left(\left|\dfrac{x}{a}\right|^r+\left|\dfrac{y}{b}\right|^r\right)^{\frac{t}{r}}+\left|\dfrac{z}{c}\right|^t-1\right)\right)
\end{bmatrix}
=
\begin{bmatrix}
0\\0\\0\\0
\end{bmatrix}
$$
to 
$$
\begin{bmatrix}
2x - \dfrac{t}{r}\lambda\left(\left(\dfrac{x}{a}\right)^r+\left(\dfrac{y}{b}\right)^r\right)^{\dfrac{t-r}{r}}\dfrac{x^{r-1}}{a^r}\\
2y - \dfrac{t}{r}\lambda\left(\left(\dfrac{x}{a}\right)^r+\left(\dfrac{y}{b}\right)^r\right)^{\dfrac{t-r}{r}}\dfrac{y^{r-1}}{b^r}\\
2z-\dfrac{t\lambda}{c^t}z^{t-1}\\
\left(\left(\dfrac{x}{a}\right)^r+\left(\dfrac{y}{b}\right)^r\right)^{\frac{t}{r}}+\left(\dfrac{z}{c}\right)^t-1
\end{bmatrix}
=
\begin{bmatrix}
0\\0\\0\\0
\end{bmatrix}
$$
to
$$
\begin{cases}
2x - \dfrac{t}{r}2z^{2-t}c^tt^{-1}\left(\left(\dfrac{x}{a}\right)^r+\left(\dfrac{y}{b}\right)^r\right)^{\dfrac{t-r}{r}}\dfrac{x^{r-1}}{a^r}=0\\
2y - \dfrac{t}{r}2z^{2-t}c^tt^{-1}\left(\left(\dfrac{x}{a}\right)^r+\left(\dfrac{y}{b}\right)^r\right)^{\dfrac{t-r}{r}}\dfrac{y^{r-1}}{b^r}=0\\
\left(\left(\dfrac{x}{a}\right)^r+\left(\dfrac{y}{b}\right)^r\right)^{\frac{t}{r}}+\left(\dfrac{z}{c}\right)^t=1
\end{cases}
$$
where I'm pretty much stuck, and sympy won't help.
My second attempt is converting the implicit equation of a superellipsoid into its polar form, then using implicit differentiation to get $\dfrac{\partial r}{\partial\theta} = 0$ and $\dfrac{\partial r}{\partial\theta}=0$. I get:
$$
\left(\left(\dfrac{r \sin{\theta} \cos{\phi}}{a}\right)^{n} + \left(\dfrac{r \sin{\phi} \sin{\theta}}{b}\right)^{n}\right)^{\frac{t}{n}} + \left(\dfrac{r \cos{\theta}}{c}\right)^{t} = 1
$$
and sympy gives me:
$$
\begin{cases}
\frac{r \left(\left(\frac{r \cos{\left (\theta \right )}}{c}\right)^{t} \sin^{2}{\left (\theta \right )} - \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n}\right)^{\frac{t}{n}} \cos^{2}{\left (\theta \right )}\right)}{\left(\left(\frac{r \cos{\left (\theta \right )}}{c}\right)^{t} + \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n}\right)^{\frac{t}{n}}\right) \sin{\left (\theta \right )} \cos{\left (\theta \right )}} = 0 \\
\frac{r \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} \sin^{2}{\left (\phi \right )} - \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n} \cos^{2}{\left (\phi \right )}\right) \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n}\right)^{\frac{t}{n}}}{\left(\left(\frac{r \cos{\left (\theta \right )}}{c}\right)^{t} \left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \cos{\left (\theta \right )}}{c}\right)^{t} \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n} + \left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n}\right)^{\frac{t}{n}} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n} \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n}\right)^{\frac{t}{n}}\right) \sin{\left (\phi \right )} \cos{\left (\phi \right )}} = 0
\end{cases}
$$
which boils down to:
$$
\begin{cases}
\left(\frac{r \cos{\left (\theta \right )}}{c}\right)^{t} \sin^{2}{\left (\theta \right )} = \left(\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n}\right)^{\frac{t}{n}} \cos^{2}{\left (\theta \right )} \\
\begin{cases}
\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} \sin^{2}{\left (\phi \right )} = \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n} \cos^{2}{\left (\phi \right )}  \\
or \\
\left(\frac{r \sin{\left (\theta \right )} \cos{\left (\phi \right )}}{a}\right)^{n} + \left(\frac{r \sin{\left (\phi \right )} \sin{\left (\theta \right )}}{b}\right)^{n} = 0
\end{cases}
\end{cases}
$$
at which point sympy won't help once again, and I've no idea how to solve for $\theta$ and $\phi$.
Anyways, has someone got an alternative method for solving this? Or is it just me who can't solve the equations?

Use the change of variable
$$\left(\frac xa\right)^r=\rho^{r/t}\cos^2\theta,
\\\left(\frac yb\right)^r=\rho^{r/t}\sin^2\theta.$$
Then the constraint becomes
$$\rho+\left(\frac zc\right)^t-1=0$$ or
$$z=c(1-\rho)^{1/t}.$$
The function to be maximized is now
$$\rho^{2/t}(a^2\cos^{4/r}\theta+b^2\sin^{4/r}\theta)+c^2(1-\rho)^{2/t}.$$
Differentiating on $\theta$ and canceling,
$$-a^2\sin\theta\cos^{4/r-1}\theta+b^2\cos\theta\sin^{4/r-1}\theta=0$$ or
$$\tan^{4/r-2}\theta=\frac{a^2}{b^2}$$ gives you the solutions in $\theta$.
And differentiating on $\rho$,
$$\rho^{2/t-1}(a^2\cos^{4/r}\theta+b^2\sin^{4/r}\theta)-c^2(1-\rho)^{2/t-1}$$
gives
$$\frac{1-\rho}{\rho}=\left(\frac{a^2\cos^{4/r}\theta+b^2\sin^{4/r}\theta}{c^2}\right)^{t/(2-t)}.$$
