
Table of contents
- Introduction
-
Implementation 2.1 Unconstrained optimization 2.2 Bounds 2.3 Linear constraints 2.4 Nonlinear constraints 2.5 Applying different constraint types together
- Conclusions
1. Introduction
Optimization is the process of picking the best elements from a set of potential candidates to reach a specific goal.
We perform a lot of optimization tasks in our everyday life: finding the shortest or fastest route to reach a destination, preparing a to-do list with daily assignments ordered by priority, buying groceries.
We can describe such problems starting with the definition of an objective function f(x)
.
Let us imagine that we are organizing a journey to another city, and we are trying assess a suitable departure time. In this example, the objective function f(x)
is the duration of the trip as function of the departure time x
.
We can formulate an optimization problem as the identification of the minimum or maximum value of the objective function. In our example, we want to determine the departure time that will minimize the duration of the trip:

In other scenarios, we may want to maximize f(x)
. For instance, when the objective represents a likelihood or a return of investment. Nevertheless, maximizing a function is equivalent to minimizing its negative. Therefore, one may focus on minimization problems alone:

In real-world applications, we may need to apply constraints to our optimization problem. For example, we may want to find the fastest route, but we are unwilling to pay tolls, or travel at night. We define constrained optimization as the process of minimizing the objective function under some logical conditions that may reflect:
- real-world limitations;
- the physical meaning of the input variables;
- contextual circumstances.
In this post, we share an optimization example using [SciPy](https://scipy.org/)
, a popular Python library for scientific computing. In particular, we explore the most common constraint types: bounds, linear and nonlinear constraints.
2. Implementation
2.1 Unconstrained optimization
We start from a simple unconstrained optimization problem, and add constraints to the input variables later on.
Import the needed libraries:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, Bounds, LinearConstraint, NonlinearConstraint
Imagine the following multivariable objective function:

Its gradient with respect to x₀
and x₁
is

def f(x):
'''Objective function'''
return 0.5*x[0]**2 + 2.5*x[1]**2 + 4 * x[0] * np.sin(np.pi * x[0]) + 5
def df(x):
'''Gradient of the objective function'''
return np.array([x[0] + 4 * np.sin(np.pi * x[0]) + 4 * np.pi * x[0] * np.cos(np.pi * x[0]), 5*x[1]])
Let us generate data and observe the function values for x₀, x₁ ∈ [-1, 1]
:
# Generate data
X0, X1 = np.meshgrid(np.linspace(-1, 1, 100), np.linspace(-1, 1, 100))
Z = f(np.stack([X0, X1]))
# Plot
fig = plt.figure(figsize=(30, 20))
# First subplot
ax = fig.add_subplot(1, 3, 1, projection='3d')
ax.contour3D(X0, X1, Z, 70, cmap='plasma')
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
# Second subplot
ax = fig.add_subplot(1, 3, 2, projection='3d')
ax.contour3D(X0, X1, Z, 70, cmap='plasma')
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.yaxis.set_ticklabels([])
ax.view_init(0, 80)
# Third subplot
ax = fig.add_subplot(1, 3, 3, projection='3d')
ax.contour3D(X0, X1, Z, 70, cmap='plasma')
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.zaxis.set_ticklabels([])
ax.view_init(89, -90);

In our example, the objective function is non-convex and possesses several minima.
This implies that, depending on the starting point, the problem may converge to a different minimizer.
We can solve the optimization problem by leveraging the helpful [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)
function as follows:
# Starting point
x_start = np.array([0.5, 0.5])
# Optimization
result = minimize(
f, x_start, method='trust-constr', jac=df)
result.x.round(3)
Notably, we are applying the trust-constr
method. It allows to optimize a function subject to constraints. More information on the method is available in the package documentation and in "Trust-region methods" (Conn, Gould and Toint; 2000).
The above code snippet returns the found minimizer:
array([-0., 0.])
Let us plot it:
# Minimum from unconstrained optimization
min_x0, min_x1 = np.meshgrid(result.x[0], result.x[1])
min_z = f(np.stack([min_x0, min_x1]))
# Plot
fig = plt.figure(figsize=(30, 20))
# First subplot
ax = fig.add_subplot(1, 3, 1, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
# Second subplot
ax = fig.add_subplot(1, 3, 2, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.yaxis.set_ticklabels([])
ax.view_init(0, 80)
# Third subplot
ax = fig.add_subplot(1, 3, 3, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.zaxis.set_ticklabels([])
ax.view_init(89, -90);

We can now experiment the addition of constraints.
2.2 Bounds
Let us consider our previous example about finding the fastest trip between two cities, with the departure time as input variable. We may expect to find more or less traffic depending on the hour of the day. By aiming at minimizing the duration of the trip, a model may also suggest, for instance, to travel at night.
Although this might result in the shortest trip, we may prefer to travel during daytime. So, we could ask the model to find the shortest trip considering departure times ranging from 7.00 AM to 6.00 PM only.
This is where bounds come in. Bounds are simply equality or inequality constraints on input variables. They allow to evaluate the objective function solely between the specified ranges.
In our case, suppose to have the following acceptable values for x₀
and x₁
:

We can easily pass these values to the [Bounds](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.Bounds.html)
object and perform a new optimization experiment as follows:
lim = [0.25, 0.30, 0.75, 0.8]
bounds = Bounds([lim[0], lim[1]], # [min x0, min x1]
[lim[2], lim[3]]) # [max x0, max x1]
result = minimize(
f, x_start, method='trust-constr', jac=df, bounds=bounds)
result.x.round(3)
The optimization task leads now to a different solution, as the previous point array([0 , 0])
does not fall in the feasible region:
array([0.25, 0.301])
We can finally plot the new minimum and the feasible region, and observe the area in which f(x)
was evaluated:
# Feasible region (bounds)
X0_bound, X1_bound = np.meshgrid(np.linspace(lim[0], lim[2], 20), np.linspace(lim[1], lim[3], 20))
Z_bound = f(np.stack([X0_bound, X1_bound]))
# New minimum within bounds
min_x0_bounds, min_x1_bounds = np.meshgrid(result.x[0], result.x[1])
min_z_bounds = f(np.stack([min_x0_bounds, min_x0_bounds]))
# Plot
fig = plt.figure(figsize=(30, 20))
# First subplot
ax = fig.add_subplot(1, 3, 1, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.scatter(min_x0_bounds, min_x1_bounds, min_z_bounds, marker='o', color='blue', linewidth=10)
ax.plot_surface(X0_bound, X1_bound, Z_bound, color='black', alpha=0.6)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
# Second subplot
ax = fig.add_subplot(1, 3, 2, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.scatter(min_x0_bounds, min_x1_bounds, min_z_bounds, marker='o', color='blue', linewidth=10)
ax.plot_surface(X0_bound, X1_bound, Z_bound, color='black', alpha=0.6)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.yaxis.set_ticklabels([])
ax.view_init(0, 80)
# Third subplot
ax = fig.add_subplot(1, 3, 3, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.scatter(min_x0_bounds, min_x1_bounds, min_z_bounds, marker='o', color='blue', linewidth=10)
ax.plot_surface(X0_bound, X1_bound, Z_bound, color='black', alpha=0.6)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.zaxis.set_ticklabels([])
ax.view_init(89, -90);

2.3 Linear constraints
Linear constraints define linear relations between the optimization variables. For example, let us imagine x₀
and x₁
to be subject to

We can easily rewrite these conditions as a linear system, and pass them to the [LinearConstraint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.LinearConstraint.html)
object before running the optimization task:

linear_constraint = LinearConstraint(
[[1, 0], [1, -1]], [0.25, -np.inf], [np.inf, 0.1])
result = minimize(
f, x_start, method='trust-constr', jac=df, constraints=linear_constraint)
result.x.round(3)
The new solution is
array([0.25, 0.15])
The feasible region of f(x)
corresponds to a portion of the space delimited by an intersection of hyperplanes. Let us plot these boundaries:
# Linear constraints: first hyperplane
X0_lin_1 = np.repeat(0.25, 20)
X1_lin_1, Z_lin_1 = np.meshgrid(np.linspace(-1, 1, 20), np.linspace(4, 10, 20))
# Linear constraints: second hyperplane
X1_lin_2 = np.linspace(-1, 1, 20)
X0_lin_2 = X1_lin_2 + 0.1
# New minimum with linear constraints
min_x0_lin_constr, min_x1_lin_constr = np.meshgrid(result.x[0], result.x[1])
min_z_lin_constr = f(np.stack([min_x0_lin_constr, min_x0_lin_constr]))
# Plot
fig = plt.figure(figsize=(30, 20))
# First subplot
ax = fig.add_subplot(1, 3, 1, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.scatter(min_x0_lin_constr, min_x1_lin_constr, min_z_lin_constr, marker='o', color='blue', linewidth=10)
ax.plot_surface(X0_lin_1, X1_lin_1, Z_lin_1, color='green', alpha=0.3)
ax.plot_surface(X0_lin_2, X1_lin_2, Z_lin_1, color='yellow', alpha=0.2)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
# Second subplot
ax = fig.add_subplot(1, 3, 2, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.scatter(min_x0_lin_constr, min_x1_lin_constr, min_z_lin_constr, marker='o', color='blue', linewidth=10)
ax.plot_surface(X0_lin_1, X1_lin_1, Z_lin_1, color='green', alpha=0.2)
ax.plot_surface(X0_lin_2, X1_lin_2, Z_lin_1, color='yellow', alpha=0.2)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.yaxis.set_ticklabels([])
ax.view_init(0, 80)
# Third subplot
ax = fig.add_subplot(1, 3, 3, projection='3d')
ax.contour3D(X0, X1, Z, 40, cmap='plasma')
ax.scatter(min_x0, min_x1, min_z, marker='o', color='red', linewidth=10)
ax.scatter(min_x0_lin_constr, min_x1_lin_constr, min_z_lin_constr, marker='o', color='blue', linewidth=10)
ax.plot_surface(X0_lin_1, X1_lin_1, Z_lin_1, color='green', alpha=1)
ax.plot_surface(X0_lin_2, X1_lin_2, Z_lin_1, color='yellow', alpha=1)
ax.set_xlabel('$x_{0}$')
ax.set_ylabel('$x_{1}$')
ax.set_zlabel('$f(x)$')
ax.axes.zaxis.set_ticklabels([])
ax.view_init(89, -90);

2.4 Nonlinear constraints
We can also explore the objective function within the region defined by nonlinear constraints using the [NonlinearConstraint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.NonlinearConstraint.html)
object. Assume that x₀
and x₁
are subject to:

We optimize f(x)
as follows:
non_linear_eq= lambda x: x[0]**2 + x[1]**2
non_linear_constr = NonlinearConstraint(
non_linear_eq, 0.2, np.inf)
result = minimize(
f, np.array([0.5, 1]), method='trust-constr', jac=df, constraints=non_linear_constr)
result.x.round(3)
array([-0., 0.447])
Similarly to the previous examples, we could observe the objective and the found minimizer given the current constraint:

2.5 Applying different constraint types together
We can combine bounds as well as linear and nonlinear constraints as follows:
result = minimize(
f,
x_start,
method='trust-constr',
jac=df,
bounds=bounds,
constraints=[linear_constraint, non_linear_constr]
)
result.x.round(3)
array([0.25, 0.371])
We remark that not all optimization methods support bounds and/or constraints. Additional information can be found in the package documentation.
3. Conclusions
In this post, we explored different types of optimization constraints. In particular, we shared practical Python examples using the SciPy
library. The examples come with plots that allow to visually inspect the different constraints.
Related posts:
Survival Analysis: Optimize the Partial Likelihood of the Cox Model