Solving system of polynomial equations stemming from a constrained optimization problem

Question

Constrained optimization (maximization and minimization) problem with two equality constraints.

$$\begin{array}{ll} \text{extremize} & x y z\\ \text{subject to} & x^2 + y^2 + z^2 = 9\\ & x + y - z = 3\end{array}$$

The Lagrangian, $\mathcal{L} : \mathbb{R}^5 \to \mathbb{R}$, is

$$\mathcal{L}(x,y,z,\mu_1,\mu_2) = xyz - \mu_1(x^2 + y^2 + z^2 - 9) - \mu_2(x + y - z - 3)$$

The first-order (necessary) conditions (FOCs) are

1. $yz - 2\mu_1x - \mu_2 = 0$;
2. $xz - 2\mu_1y - \mu_2 = 0$;
3. $xy - 2\mu_1z + \mu_2 = 0$;
4. $x^2 + y^2 + z^2 = 9$;
5. $x + y - z = 3$.

Question: How can this system of equations best be solved? In other words, how can I find all candidate constrained maximizers or minimizers?

When I tried to solve the system (by first solving (1),(2),(3) for $\mu_2$ and consequently the resulting equation for $\mu_1$, or by first solving (1),(2),(3) for $\mu_1$ and consequently the resulting equation for $\mu_2$)
I arrived at $x = y = -z$, but in that case we only reach $(1,1,1)$ (by (5)), which violates the other equality constraint (4), so these approaches apparentely do not lead to any solutions, but there definitely are (several) solutions (as I've checked on WolframAlpha).

How to find these?

Answer 1

You can reduce the problem in the beginning to: $$\text{optimize} \ xy(x+y-3) \ \text{subject to} \\ \ x^2+y^2+(x+y-3)^2=9.$$ Hence: $$L(x,y,\lambda)=x^2y+xy^2-3xy+\lambda(9-x^2-y^2-(x+y-3)^2).$$ FOC: $$\begin{cases} L_x=2xy+y^2-3y-2x\lambda-2\lambda(x+y-3)=0\\ L_y=x^2+2xy-3x-2y\lambda-2\lambda(x+y-3)=0\\ L_{\lambda}=9-x^2-y^2-(x+y-3)^2=0\\ \end{cases}.$$ $(2)-(1)$: $$x^2-y^2-3(x-y)+2\lambda(x-y)=0 \Rightarrow (x-y)(x+y-3+2\lambda)=0.$$ Case 1: $x-y=0 \Rightarrow y=x$, which is plugged to $(3)$: $$9-2x^2-(2x-3)^2=0 \Rightarrow x_1=0,x_2=2.$$ Can you continue?

Answer 2

Solving systems of polynomial equations is easier with symbolic computation. Using SymPy:

from sympy import *

x, y, z, mu1, mu2 = symbols('x y z mu1 mu2', real=True)

# define Lagrangian
L = (x * y * z) - mu1 * (x**2 + y**2 + z**2 - 9) - mu2 * (x + y - z - 3)

# compute partial derivatives
partials = [ diff(L,var) for var in [x, y, z, mu1, mu2] ]

# find where partial derivatives vanish
print solve_poly_system(partials, x, y, z, mu1, mu2)

This script outputs the following list of solutions.

[(-1, 2, -2, 1, -2), (0, 0, -3, 0, 0), (0, 3, 0, 0, 0), (2, -1, -2, 1, -2), (2, 2, 1, 1, -2), (3, 0, 0, 0, 0)]