Constrained Optimisation

February 16, 2025

Unconstrained Optimisation

Consider the problem

\[ \min_{x_1, x_2} f(x_1, x_2) \]

where \(f(x_1, x_2) = x_1^2 + x_2^2 - 3x_1x_2\).

The gradient of this function is

\[ \nabla f(x_1, x_2) = \begin{bmatrix} 2x_1 - 3x_2 \\ 2x_2 - 3x_1 \end{bmatrix} \]

The gradient is zero at the minima of the function. Setting the gradient to zero, we get

\[ \begin{align*} 2x_1 - 3x_2 &= 0 \\ 2x_2 - 3x_1 &= 0 \end{align*} \]

Solving these equations, we get \(x_1 = 0\) and \(x_2 = 0\). The Hessian of the function is

\[ H = \begin{bmatrix} 2 & -3 \\ -3 & 2 \end{bmatrix} \]

The eigenvalues

\[ \begin{align*} \lambda_1 &= 5 \\ \lambda_2 &= -1 \end{align*} \]

# Calculating the eigenvalues of the Hessian

import numpy as np

H = np.array([[2, -3], [-3, 2]])
eigenvalues = np.linalg.eig(H)[0]
print(eigenvalues)

[ 5. -1.]

So Hessian is not positive definite or negative definite at \((0,0)\). Hence, it is not a minima or a maxima. It is a saddle point.

In such a situation, the minima lies on the constraints. Hence, the problem defined above is unsolvable in the current form.

Linear Constraints

Example 0

Lets add bound constraints to the problem.

\[ \begin{align*} \min_{x_1, x_2} &\ x_1^2 + x_2^2 -3x_1 x_2 \\ \text{subject to} &\\ & x_1 \in [-10, 10] \\ & x_2 \in [-10, 10] \end{align*} \]

import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(-10, 10, 100)
x2 = np.linspace(-10, 10, 100)

X1, X2 = np.meshgrid(x1, x2)

Z = X1**2 + X2**2 - 3*X1*X2

plt.contour(X1, X2, Z, 20)
#Change the color map
plt.contourf(X1, X2, Z, 20, cmap='RdGy')
# setting the grid on
plt.grid()
# Plotting the minima in bold red
plt.plot(10, 10, 'bo', markersize=10)
plt.plot(-10, -10, 'bo', markersize=10)
# Show the temperature on the right
plt.colorbar()
# equal aspect ratio
plt.axis('equal')

plt.show()

This is unsatisfactory. We still need an automated method to solve this problem. If we see the above problem, it consists of four inequality constraints in the standard form.

\[ \begin{align*} g_1(x_1, x_2) &= x_1 - 10 \leq 0 \\ g_2(x_1, x_2) &= x_2 - 10 \leq 0 \\ g_3(x_1, x_2) &= -x_1 - 10 \leq 0 \\ g_4(x_1, x_2) &= -x_2 - 10 \leq 0 \end{align*} \]

This is a complex problem. So lets devise simpler problems and devise a method.

Example 1 (Single equality constraint)

\[ \begin{align*} \min_{x_1, x_2} & \ x_1^2 + x_2^2 -3x_1 x_2 \\ \text{subject to} &\\ & x_1 + 2x_2 -8 = 0 \end{align*} \]

Lets solve it graphically.

import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(-10, 10, 100)
x2 = np.linspace(-10, 10, 100)

X1, X2 = np.meshgrid(x1, x2)

Z = X1**2 + X2**2 - 3*X1*X2

plt.contour(X1, X2, Z, 20)
# setting the grid on
plt.grid()
# equal aspect ratio
plt.axis('equal')
# Show the temperature on the right
plt.colorbar()

G = X1 + 2*X2 - 8

plt.contour(X1, X2, G, levels=[0], colors='b')

# Showing the direction of the gradient of Z along the constraint

for i in range(-10,10,1):
    x = i
    y = (8 - x)/2
    plt.quiver(x, y, 2*x - 3*y, 2*y - 3*x, color='blue')


plt.show()

Now lets solve it algebraically.

First we define a Lagrangian, which is

\[ L(x_1, x_2, \lambda) = f(x_1, x_2) + \lambda g(x_1, x_2) \]

The gradient of the Lagrangian is

\[ \nabla L(x_1, x_2, \lambda) = \begin{bmatrix} 2x_1 - 3x_2 + \lambda \\ 2x_2 - 3x_1 + 2\lambda \\ x_1 + 2x_2 - 8 \end{bmatrix} \]

Setting the gradient to zero, we get

\[ \begin{align*} 2x_1 - 3x_2 + \lambda &= 0 \\ 2x_2 - 3x_1 + 2\lambda &= 0 \\ x_1 + 2x_2 - 8 &= 0 \end{align*} \]

Solving these equations we get

# Solving the above equations

from sympy import symbols, Eq, solve

x1, x2, lambda_ = symbols('x1 x2 lambda')

eq1 = Eq(2*x1 - 3*x2 + lambda_, 0)
eq2 = Eq(2*x2 - 3*x1 + 2*lambda_, 0)
eq3 = Eq(x1 + 2*x2 - 8, 0)

solution = solve((eq1, eq2, eq3), (x1, x2, lambda_))
print(solution)

{lambda: 20/11, x1: 32/11, x2: 28/11}

Exampe 2 (Single inequality constraint)

\[ \begin{align*} \min_{x_1, x_2} &\ x_1^2 + x_2^2 - 3x_1 x_2 \\ \text{subject to} & \\ & x_1 + 2x_2 -8 \leq 0 \end{align*} \]

Lets solve it graphically.

import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(-10, 10, 100)
x2 = np.linspace(-10, 10, 100)

X1, X2 = np.meshgrid(x1, x2)

Z = X1**2 + X2**2 - 3*X1*X2

plt.contour(X1, X2, Z, 20)
# setting the grid on
plt.grid()
# equal aspect ratio
plt.axis('equal')
# Show the temperature on the right
plt.colorbar()

G = X1 + 2*X2 - 8

plt.contour(X1, X2, G, levels=[0], colors='b')

# Showing the feasible region

plt.fill_between(x1, (8 - x1)/2, -10, color='gray', alpha=0.5)

# Showing the direction of the gradient of Z along the constraint

for i in range(-10,10,1):
    x = i
    y = (8 - x)/2
    plt.quiver(x, y, 2*x - 3*y, 2*y - 3*x, color='blue')


plt.show()

In this case, the Lagrangian is

\[ L(x_1, x_2, \lambda, s) = f(x_1, x_2) + \lambda (g(x_1, x_2) + s^2) \]

The gradient of the Lagrangian is

\[ \nabla L(x_1, x_2, \lambda, s) = \begin{bmatrix} 2x_1 - 3x_2 + \lambda \\ 2x_2 - 3x_1 + 2\lambda \\ x_1 + 2x_2 - 8 + s^2 \\ 2\lambda s\end{bmatrix} \]

Setting the gradient to zero, we get

\[ \begin{align*} 2x_1 - 3x_2 + \lambda &= 0 \\ 2x_2 - 3x_1 + 2\lambda &= 0 \\ x_1 + 2x_2 - 8 + s^2 &= 0 \\ 2\lambda s &= 0 \end{align*} \]

Here we have switching conditions.

Case 1: \(s = 0\) Case 2: \(\lambda = 0\) Case 3: \(s \neq 0\) and \(\lambda \neq 0\)

All the solutions can be obtained by solving the above equations.

# Solving the above equations

from sympy import symbols, Eq, solve

x1, x2, lambda_, s = symbols('x1 x2 lambda s')

eq1 = Eq(2*x1 - 3*x2 + lambda_, 0)
eq2 = Eq(2*x2 - 3*x1 + 2*lambda_, 0)
eq3 = Eq(x1 + 2*x2 - 8 + s**2, 0)
eq4 = Eq(2*lambda_*s, 0)

solution = solve((eq1, eq2, eq3, eq4), (x1, x2, lambda_, s))
print(solution)

[(0, 0, 0, -2*sqrt(2)), (0, 0, 0, 2*sqrt(2)), (32/11, 28/11, 20/11, 0)]

So which of these three is the correct solution? Lets check the feasibility of the solutions.

# Checking the feasibility of the solutions

for sol in solution:
    x1 = sol[0]
    x2 = sol[1]
    if x1 + 2*x2 - 8 <= 0:
        print(sol)

(0, 0, 0, -2*sqrt(2))
(0, 0, 0, 2*sqrt(2))
(32/11, 28/11, 20/11, 0)

All are feasible. So which is the minimum? Lets check the value of the function at these points.

# Checking the value of the function at these points

for sol in solution:
    x1 = sol[0]
    x2 = sol[1]
    print(x1**2 + x2**2 - 3*x1*x2)

0
0
-80/11

So clearly the third point is the minimum.

Example 3 (Single inequality constraint)

\[ \begin{align*} \min_{x_1, x_2} &\ x_1^2 + x_2^2 - 3x_1 x_2 \\ \text{subject to} & \\ & x_1 + 2x_2 -8 \geq 0 \end{align*} \]

Lets solve it graphically.

import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(-10, 10, 100)
x2 = np.linspace(-10, 10, 100)

X1, X2 = np.meshgrid(x1, x2)

Z = X1**2 + X2**2 - 3*X1*X2

plt.contour(X1, X2, Z, 20)
# setting the grid on
plt.grid()
# equal aspect ratio
plt.axis('equal')
# Show the temperature on the right
plt.colorbar()

G = - X1 - 2*X2 + 8

plt.contour(X1, X2, G, levels=[0], colors='b')

# Showing the feasible region

plt.fill_between(x1, (8 - x1)/2, 10, color='gray', alpha=0.5)

# Showing the direction of the gradient of Z along the constraint

for i in range(-10,10,1):
    x = i
    y = (8 - x)/2
    plt.quiver(x, y, 2*x - 3*y, 2*y - 3*x, color='blue')


plt.show()

In this case, the Lagrangian is

\[ L(x_1, x_2, \lambda, s) = f(x_1, x_2) + \lambda (g(x_1, x_2) - s^2) \]

The gradient of the Lagrangian is

\[ \nabla L(x_1, x_2, \lambda, s) = \begin{bmatrix} 2x_1 - 3x_2 - \lambda \\ 2x_2 - 3x_1 - 2\lambda \\ -x_1 - 2x_2 + 8 + s^2 \\ 2\lambda s \end{bmatrix} \]

Setting the gradient to zero, we get

\[ \begin{align*} 2x_1 - 3x_2 - \lambda &= 0 \\ 2x_2 - 3x_1 - 2\lambda &= 0 \\ -x_1 - 2x_2 + 8 + s^2 &= 0 \\ 2\lambda s &= 0 \end{align*} \]

Here we have switching conditions.

Case 1: \(s = 0\) Case 2: \(\lambda = 0\) Case 3: \(s = 0\) and \(\lambda = 0\)

All the solutions can be obtained by solving the above equations.

# Solving the above equations

from sympy import symbols, Eq, solve

x1, x2, lambda_, s = symbols('x1 x2 lambda s')

eq1 = Eq(2*x1 - 3*x2 - lambda_, 0)
eq2 = Eq(2*x2 - 3*x1 - 2*lambda_, 0)
eq3 = Eq(-x1 - 2*x2 + 8 + s**2, 0)
eq4 = Eq(2*lambda_*s, 0)

solution = solve((eq1, eq2, eq3, eq4), (x1, x2, lambda_, s))
print(solution)

[(32/11, 28/11, -20/11, 0), (0, 0, 0, -2*sqrt(2)*I), (0, 0, 0, 2*sqrt(2)*I)]

So which of these three is the correct solution? Lets check the feasibility of the solutions.

# Checking the feasibility of the solutions

for sol in solution:
    x1 = sol[0]
    x2 = sol[1]
    if -x1 - 2*x2 + 8 >= 0:
        print(sol)

(32/11, 28/11, -20/11, 0)
(0, 0, 0, -2*sqrt(2)*I)
(0, 0, 0, 2*sqrt(2)*I)

All are feasible. So which is the minimum? Lets check the value of the function at these points.

# Checking the value of the function at these points

for sol in solution:
    x1 = sol[0]
    x2 = sol[1]
    print(x1**2 + x2**2 - 3*x1*x2)

-80/11
0
0

Non-linear Constraints

Example 3

\[ \begin{align*} \min_{x_1, x_2} & x_1^2 + x_2^2 - 3x_1 x_2 \\ \text{subject to} & \\ & x_1, x_2 \in [-10, 10] \end{align*} \]

The Lagrangian of the problem is

\[ L(x_1, x_2, \lambda_1, \lambda_2, \lambda_3, \lambda_4, s_1, s_2, s_3, s_4) = f(x_1, x_2) + \lambda_1 (x_1 - 10 + s_1^2) + \lambda_2 (x_2 - 10 +s_2^2) +\lambda_3 (-x_1 - 10 + s_3^2) + \lambda_4 (-x_2 - 10 + s_4^2) \]

The gradient of the Lagrangian is

\[ \nabla L(x_1, x_2, \lambda_1, \lambda_2, \lambda_3, \lambda_4, s_1, s_2, s_3, s_4) = \begin{bmatrix} 2x_1 - 3x_2 + \lambda_1 - \lambda_3 \\ 2x_2 - 3x_1 + \lambda_2 - \lambda_4 \\ x_1 - 10 + s_1^2 \\ x_2 - 10 + s_2^2 \\ -x_1 - 10 + s_3^2 \\ -x_2 - 10 + s_4^2 \\ 2\lambda_1 s_1 \\ 2\lambda_2 s_2 \\ 2\lambda_3 s_3 \\ 2\lambda_4 s_4 \end{bmatrix} \]

Setting these equations to zero, we get

# Solving the above equations

from sympy import symbols, Eq, solve

x1, x2, lambda1, lambda2, lambda3, lambda4, s1, s2, s3, s4 = symbols('x1 x2 lambda1 lambda2 lambda3 lambda4 s1 s2 s3 s4')

eq1 = Eq(2*x1 - 3*x2 + lambda1 - lambda3, 0)
eq2 = Eq(2*x2 - 3*x1 + lambda2 - lambda4, 0)
eq3 = Eq(x1 - 10 + s1**2, 0)
eq4 = Eq(x2 - 10 + s2**2, 0)
eq5 = Eq(-x1 - 10 + s3**2, 0)
eq6 = Eq(-x2 - 10 + s4**2, 0)
eq7 = Eq(2*lambda1*s1, 0)
eq8 = Eq(2*lambda2*s2, 0)
eq9 = Eq(2*lambda3*s3, 0)
eq10 = Eq(2*lambda4*s4, 0)

solution = solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10), (x1, x2, lambda1, lambda2, lambda3, lambda4, s1, s2, s3, s4))

print(solution)

[(10, -10, -50, 0, 0, -50, 0, -2*sqrt(5), -2*sqrt(5), 0), (10, -10, -50, 0, 0, -50, 0, -2*sqrt(5), 2*sqrt(5), 0), (10, -10, -50, 0, 0, -50, 0, 2*sqrt(5), -2*sqrt(5), 0), (10, -10, -50, 0, 0, -50, 0, 2*sqrt(5), 2*sqrt(5), 0), (-10, 10, 0, -50, -50, 0, -2*sqrt(5), 0, 0, -2*sqrt(5)), (-10, 10, 0, -50, -50, 0, -2*sqrt(5), 0, 0, 2*sqrt(5)), (-10, 10, 0, -50, -50, 0, 2*sqrt(5), 0, 0, -2*sqrt(5)), (-10, 10, 0, -50, -50, 0, 2*sqrt(5), 0, 0, 2*sqrt(5)), (0, 0, 0, 0, 0, 0, -sqrt(10), -sqrt(10), -sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), -sqrt(10), -sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), -sqrt(10), sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), -sqrt(10), sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), sqrt(10), -sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), sqrt(10), -sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), sqrt(10), sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, -sqrt(10), sqrt(10), sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), -sqrt(10), -sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), -sqrt(10), -sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), -sqrt(10), sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), -sqrt(10), sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), sqrt(10), -sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), sqrt(10), -sqrt(10), sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), sqrt(10), sqrt(10), -sqrt(10)), (0, 0, 0, 0, 0, 0, sqrt(10), sqrt(10), sqrt(10), sqrt(10)), (-15, -10, 0, 0, 0, 25, -5, -2*sqrt(5), -sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, -5, -2*sqrt(5), sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, -5, 2*sqrt(5), -sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, -5, 2*sqrt(5), sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, 5, -2*sqrt(5), -sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, 5, -2*sqrt(5), sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, 5, 2*sqrt(5), -sqrt(5)*I, 0), (-15, -10, 0, 0, 0, 25, 5, 2*sqrt(5), sqrt(5)*I, 0), (-10, -10, 0, 0, 10, 10, -2*sqrt(5), -2*sqrt(5), 0, 0), (-10, -10, 0, 0, 10, 10, -2*sqrt(5), 2*sqrt(5), 0, 0), (-10, -10, 0, 0, 10, 10, 2*sqrt(5), -2*sqrt(5), 0, 0), (-10, -10, 0, 0, 10, 10, 2*sqrt(5), 2*sqrt(5), 0, 0), (-10, -15, 0, 0, 25, 0, -2*sqrt(5), -5, 0, -sqrt(5)*I), (-10, -15, 0, 0, 25, 0, -2*sqrt(5), -5, 0, sqrt(5)*I), (-10, -15, 0, 0, 25, 0, -2*sqrt(5), 5, 0, -sqrt(5)*I), (-10, -15, 0, 0, 25, 0, -2*sqrt(5), 5, 0, sqrt(5)*I), (-10, -15, 0, 0, 25, 0, 2*sqrt(5), -5, 0, -sqrt(5)*I), (-10, -15, 0, 0, 25, 0, 2*sqrt(5), -5, 0, sqrt(5)*I), (-10, -15, 0, 0, 25, 0, 2*sqrt(5), 5, 0, -sqrt(5)*I), (-10, -15, 0, 0, 25, 0, 2*sqrt(5), 5, 0, sqrt(5)*I), (15, 10, 0, 25, 0, 0, -sqrt(5)*I, 0, -5, -2*sqrt(5)), (15, 10, 0, 25, 0, 0, -sqrt(5)*I, 0, -5, 2*sqrt(5)), (15, 10, 0, 25, 0, 0, -sqrt(5)*I, 0, 5, -2*sqrt(5)), (15, 10, 0, 25, 0, 0, -sqrt(5)*I, 0, 5, 2*sqrt(5)), (15, 10, 0, 25, 0, 0, sqrt(5)*I, 0, -5, -2*sqrt(5)), (15, 10, 0, 25, 0, 0, sqrt(5)*I, 0, -5, 2*sqrt(5)), (15, 10, 0, 25, 0, 0, sqrt(5)*I, 0, 5, -2*sqrt(5)), (15, 10, 0, 25, 0, 0, sqrt(5)*I, 0, 5, 2*sqrt(5)), (10, 10, 10, 10, 0, 0, 0, 0, -2*sqrt(5), -2*sqrt(5)), (10, 10, 10, 10, 0, 0, 0, 0, -2*sqrt(5), 2*sqrt(5)), (10, 10, 10, 10, 0, 0, 0, 0, 2*sqrt(5), -2*sqrt(5)), (10, 10, 10, 10, 0, 0, 0, 0, 2*sqrt(5), 2*sqrt(5)), (10, 15, 25, 0, 0, 0, 0, -sqrt(5)*I, -2*sqrt(5), -5), (10, 15, 25, 0, 0, 0, 0, -sqrt(5)*I, -2*sqrt(5), 5), (10, 15, 25, 0, 0, 0, 0, -sqrt(5)*I, 2*sqrt(5), -5), (10, 15, 25, 0, 0, 0, 0, -sqrt(5)*I, 2*sqrt(5), 5), (10, 15, 25, 0, 0, 0, 0, sqrt(5)*I, -2*sqrt(5), -5), (10, 15, 25, 0, 0, 0, 0, sqrt(5)*I, -2*sqrt(5), 5), (10, 15, 25, 0, 0, 0, 0, sqrt(5)*I, 2*sqrt(5), -5), (10, 15, 25, 0, 0, 0, 0, sqrt(5)*I, 2*sqrt(5), 5)]

So the optimal solutions are the points that satisfy the following conditions.

\[ \lambda_i > 0 \quad \text{and} \quad \lambda_i s_i = 0 \]

# Identifying the optimal solutions

for sol in solution:
    x1 = sol[0]
    x2 = sol[1]
    lambda1 = sol[2]
    lambda2 = sol[3]
    lambda3 = sol[4]
    lambda4 = sol[5]
    s1 = sol[6]
    s2 = sol[7]
    s3 = sol[8]
    s4 = sol[9]
    if lambda1 >= 0 and \
       lambda2 >= 0 and \
       lambda3 >= 0 and \
       lambda4 >= 0 and \
       s1.is_real and \
       s2.is_real and \
       s3.is_real and \
       s4.is_real:
    
        print(x1, x2, lambda1, lambda2, lambda3, lambda4, s1, s2, s3, s4)

0 0 0 0 0 0 -sqrt(10) -sqrt(10) -sqrt(10) -sqrt(10)
0 0 0 0 0 0 -sqrt(10) -sqrt(10) -sqrt(10) sqrt(10)
0 0 0 0 0 0 -sqrt(10) -sqrt(10) sqrt(10) -sqrt(10)
0 0 0 0 0 0 -sqrt(10) -sqrt(10) sqrt(10) sqrt(10)
0 0 0 0 0 0 -sqrt(10) sqrt(10) -sqrt(10) -sqrt(10)
0 0 0 0 0 0 -sqrt(10) sqrt(10) -sqrt(10) sqrt(10)
0 0 0 0 0 0 -sqrt(10) sqrt(10) sqrt(10) -sqrt(10)
0 0 0 0 0 0 -sqrt(10) sqrt(10) sqrt(10) sqrt(10)
0 0 0 0 0 0 sqrt(10) -sqrt(10) -sqrt(10) -sqrt(10)
0 0 0 0 0 0 sqrt(10) -sqrt(10) -sqrt(10) sqrt(10)
0 0 0 0 0 0 sqrt(10) -sqrt(10) sqrt(10) -sqrt(10)
0 0 0 0 0 0 sqrt(10) -sqrt(10) sqrt(10) sqrt(10)
0 0 0 0 0 0 sqrt(10) sqrt(10) -sqrt(10) -sqrt(10)
0 0 0 0 0 0 sqrt(10) sqrt(10) -sqrt(10) sqrt(10)
0 0 0 0 0 0 sqrt(10) sqrt(10) sqrt(10) -sqrt(10)
0 0 0 0 0 0 sqrt(10) sqrt(10) sqrt(10) sqrt(10)
-10 -10 0 0 10 10 -2*sqrt(5) -2*sqrt(5) 0 0
-10 -10 0 0 10 10 -2*sqrt(5) 2*sqrt(5) 0 0
-10 -10 0 0 10 10 2*sqrt(5) -2*sqrt(5) 0 0
-10 -10 0 0 10 10 2*sqrt(5) 2*sqrt(5) 0 0
10 10 10 10 0 0 0 0 -2*sqrt(5) -2*sqrt(5)
10 10 10 10 0 0 0 0 -2*sqrt(5) 2*sqrt(5)
10 10 10 10 0 0 0 0 2*sqrt(5) -2*sqrt(5)
10 10 10 10 0 0 0 0 2*sqrt(5) 2*sqrt(5)

Hence the optimal solution is \((10, 10)\) and \((-10, -10)\).

Example 4

\[ \begin{align*} \min_{x_1, x_2} \ f(x_1, x_2) \\ \text{subject to} \\ x_1^2 + x_2^2 - 5 \leq 0 \end{align*} \]

Lets solve it graphically.

import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(-10, 10, 100)
x2 = np.linspace(-10, 10, 100)

X1, X2 = np.meshgrid(x1, x2)

Z = X1**2 + X2**2 - 3*X1*X2

plt.contour(X1, X2, Z, 20)
#Change the color map
plt.contourf(X1, X2, Z, 20, cmap='RdGy')
# setting the grid on
plt.grid()
# equal aspect ratio
plt.axis('equal')
# Show the temperature on the right
plt.colorbar()

G = X1**2 + X2**2 - 64

plt.contour(X1, X2, G, levels=[0], colors='b')

# Plotting the gradient field of G
for i in range(0, 360, 18):
    x = 8 * np.cos(np.deg2rad(i))
    y = 8 * np.sin(np.deg2rad(i))
    plt.quiver(x, y, 2 * x, 2 * y, color='blue', scale=50)

plt.show()

The Lagrangian of the problem is

\[ L(x_1, x_2, \lambda, s) = f(x_1, x_2) + \lambda (g(x_1, x_2) + s^2) \]

The gradient of the Lagrangian is

\[ \nabla L(x_1, x_2, \lambda, s) = \begin{bmatrix} 2x_1 - 3x_2 + \lambda \\ 2x_2 - 3x_1 + 2\lambda \\ x_1^2 + x_2^2 - 5 + s^2 \\ 2\lambda s \end{bmatrix} \]

Setting the gradient to zero, we get

\[ \begin{align*} 2x_1 - 3x_2 + \lambda &= 0 \\ 2x_2 - 3x_1 + 2\lambda &= 0 \\ x_1^2 + x_2^2 - 5 + s^2 &= 0 \\ 2\lambda s &= 0 \end{align*} \]

# Solving the above equations

from sympy import symbols, Eq, solve

x1, x2, lambda_, s = symbols('x1 x2 lambda s')

eq1 = Eq(2*x1 - 3*x2 + lambda_, 0)
eq2 = Eq(2*x2 - 3*x1 + 2*lambda_, 0)
eq3 = Eq(x1**2 + x2**2 - 5 + s**2, 0)
eq4 = Eq(2*lambda_*s, 0)

solution = solve((eq1, eq2, eq3, eq4), (x1, x2, lambda_, s))
print(solution)

[(0, 0, 0, -sqrt(5)), (0, 0, 0, sqrt(5)), (-8*sqrt(565)/113, -7*sqrt(565)/113, -5*sqrt(565)/113, 0), (8*sqrt(565)/113, 7*sqrt(565)/113, 5*sqrt(565)/113, 0)]

Applying the KKT conditions,

# Checking if lambda's are positive and s is real

for sol in solution:
    x1 = sol[0]
    x2 = sol[1]
    lambda_ = sol[2]
    s = sol[3]
    if lambda_ >= 0 and s.is_real:
        print(sol)

(0, 0, 0, -sqrt(5))
(0, 0, 0, sqrt(5))
(8*sqrt(565)/113, 7*sqrt(565)/113, 5*sqrt(565)/113, 0)