原问题

$$ \max_{x,y,z} f(x,y,z)=8xyz\\ s.t.\ \frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1 $$

拉格朗日函数

$$ L(x,y,z,\lambda)=f(x,y,z)+\lambda (\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}-1) $$

对 $L$ 求偏导 $$ \frac{\partial L}{\partial x}=8yz+\frac{2\lambda x}{a^2}=0 \\ \frac{\partial L}{\partial y}=8xz+\frac{2\lambda y}{b^2}=0 \\ \frac{\partial L}{\partial z}=8xy+\frac{2\lambda z}{c^2}=0 \\ $$

import numpy as np
from scipy.optimize import minimize

e = 1e-10
f = lambda x: 8 * x[0] * x[1] * x[2]
cons = (
    {'type': 'eq', 'fun': lambda x: x[0]**2 + x[1]**2 + x[2]**2 - 1},
    {'type': 'ineq', 'fun': lambda x: x[0] - e},
    {'type': 'ineq', 'fun': lambda x: x[1] - e},
    {'type': 'ineq', 'fun': lambda x: x[2] - e},
)

x0 = np.array((1.0, 1.0, 1.0))
res = minimize(f, x0, 
               method='SLSQP',
               constraints=cons)

if res.success:
    print(res.fun)
    print(res.x)
else:
    print(res.message)

1.5396007243645415
[0.57735022 0.57735022 0.57735038]

from sympy import symbols, diff, solve

x1 = symbols('x1')
x2 = symbols('x2')
x3 = symbols('x3')
# a = symbols('a')
# b = symbols('b')
# c = symbols('c')
a, b, c = 1, 1, 1
lbd = symbols('lambda')

L = 8*x1*x2*x3 + lbd * (x1**2/a**2 + x2**2/b**2 + x3**2/c**2 - 1)

difyL_x1 = diff(L, x1)
difyL_x2 = diff(L, x2)
difyL_x3 = diff(L, x3)
difyL_lbd = diff(L, lbd)

# KKT solver
aa = solve([difyL_x1, difyL_x2, difyL_x3, difyL_lbd], [x1, x2, x3, lbd])
aa

[(-1, 0, 0, 0),
 (0, -1, 0, 0),
 (0, 0, -1, 0),
 (0, 0, 1, 0),
 (0, 1, 0, 0),
 (1, 0, 0, 0),
 (-sqrt(3)/3, -sqrt(3)/3, sqrt(3)/3, -4*sqrt(3)/3),
 (-sqrt(3)/3, sqrt(3)/3, -sqrt(3)/3, -4*sqrt(3)/3),
 (sqrt(3)/3, -sqrt(3)/3, -sqrt(3)/3, -4*sqrt(3)/3),
 (sqrt(3)/3, sqrt(3)/3, sqrt(3)/3, -4*sqrt(3)/3),
 (-sqrt(3)/3, -sqrt(3)/3, -sqrt(3)/3, 4*sqrt(3)/3),
 (-sqrt(3)/3, sqrt(3)/3, sqrt(3)/3, 4*sqrt(3)/3),
 (sqrt(3)/3, -sqrt(3)/3, sqrt(3)/3, 4*sqrt(3)/3),
 (sqrt(3)/3, sqrt(3)/3, -sqrt(3)/3, 4*sqrt(3)/3)]