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)}.$$