TDM 20200: Project 8 – Optimization with Python and SciPy

Project Objectives

Optimization is a core concept in data science and engineering. In this project, you will learn how to formulate optimization problems in code and solve them using scipy.optimize. Along the way, you will build intuition about what optimization is doing, explore real computational methods, and compare algorithmic approaches.

Learning Objectives

Understand what an optimization problem is
Define objective functions and constraints in code
Understand different optimization algorithms in SciPy
Compare methods and interpret their results
Apply optimization to real scenarios

If AI is used in any cases, such as for debugging, research, etc., we now require that you submit a link to the entire chat history. For example, if you used ChatGPT, there is an “Share” option in the conversation sidebar. Click on “Create Link” and please add the shareable link as a part of your citation.

The project template in the Examples Book now has a “Link to AI Chat History” section; please have this included in all your projects. If you did not use any AI tools, you may write “None”.

We allow using AI for learning purposes; however, all submitted materials (code, comments, and explanations) must all be your own work and in your own words. No content or ideas should be directly applied or copy pasted to your projects. Please refer to the-examples-book.com/projects/spring2026/syllabus#guidance-on-generative-ai. Failing to follow these guidelines is considered as academic dishonesty.

Questions

Question 1 (2 points)

Optimization means finding the best solution under a set of rules. In math, this is often described as minimizing (or maximizing) a function. In this question, you will start with the simplest possible case: minimizing a function of one variable. So, you will write such a function, visualize it, then use SciPy’s optimizer to find its minimum point automatically. We begin by importing our standard data libraries.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

Now let’s define an objective function ($f(x)$). The SciPy minimize function will use this function to calculate values of the objective function as it searches for the minimum. You simply need to take the input and return the output.

In this one-variable problem, the input is a single number ($x$). However, for problems with multiple variables, the input will be an array of numbers. For example, if we have a function of $x$ and $y$, the input will be $[x, y]$, or ($v$), where $v[0] = x$ and $v[1] = y$.

def objective_function(x):
    """
    A simple 1D function to minimize.

    Example: x^2 + 2x + 1 is a simple quadratic function with a minimum at $x = -1$.
    """
    return  x**2 + 2*x + 1 # TODO: put your own function here

Once your function is defined, it is important to visualize it before optimizing. You can use online graphing tools like Desmos to plot your function, or you can keep it all in Python using matplotlib.

For optimization problems, it is important to have a decent first guess at the minimum, also known as the initial guess. This is used to help the optimizer find the minimum, and for complex problems it can greatly improve the speed and reliability of the solution.

Below is some simple matplotlib code to plot your function from x = -2 to x = 6. Feel free to change the range to something that makes sense for your function.

x_vals = np.linspace(-2, 6, 300)
y_vals = objective_function(x_vals)

plt.figure(figsize=(10,5))
plt.plot(x_vals, y_vals, linewidth=2)
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("Visualization of your objective function")
plt.grid(alpha=0.3)
plt.show()

For a single variable problem, it should be pretty clear where the minimum or a minima is.

Before we run the optimizer, we need to choose which algorithm to use. SciPy’s minimize function supports many methods, each with different strengths. Some are gradient-based, some are derivative-free, and some handle constraints or bounds. Here’s a short overview:

Common optimization methods in SciPy
Method	Good for	Constraints	Notes
`BFGS`	Smooth, unconstrained problems	No	Gradient-based; efficient when derivatives are smooth
`Nelder-Mead`	Noisy or derivative-unavailable situations	No	Simplex method; slower but robust
`Powell`	Derivative-free	No	Works fairly well without gradients
`L-BFGS-B`	Smooth with bounds	Bounds only	Like BFGS but supports box constraints
`SLSQP`	Problems with constraints	Yes	Frequently used for constrained optimization

For this question, we will use BFGS. You will soon have a chance to experiment with others.

The following cell runs the optimizer using BFGS starting from an initial guess. Try different initial guesses and see how the result changes.

initial_guess = 0.0  # You can change this
result = minimize(objective_function, initial_guess, method="BFGS")

print(f"Minimum value: {result.fun:.6f}")
print(f"At x = {result.x[0]:.6f}")
print(f"Success? {result.success}")

if not result.success:
    print("Message:", result.message)

Deliverables

1.1 Define a 1D objective function.
1.2 Plot it and build intuition about its minimum.
1.3 Use SciPy’s optimizer with BFGS to find the minimum.
1.4 Briefly reflect on how the numeric result compares to the plot.
1.5 Does the minimum value and location match what you saw in the plot?
1.6 What happens if you change initial_guess?

Question 2 (2 points)

In many real-life problems, we must respect rules or limits. In optimization, these are called constraints. For example, “x must be positive” or “x must be no greater than 5.”

SciPy’s minimize function supports constraints through its constraints parameter. In this question, you will define an objective similar to before, but you will also constrain the search to a feasible region.

Let’s start by defining a simple objective function of one variable. There is an example function in the code below. Define your own objective function. It should look like something you might have used in Question 1, but you’re welcome to write a new one.

def constrained_objective(x):
    return  x[0]**2 + 2*x[0] + 1 # TODO: define a 1D objective using x[0]

Now let’s discuss constraints. In SciPy, an inequality constraint is a function defined so that:

constraint_fun() >= 0

For example:

To require x >= 0: lambda x: x[0]
To require x ⇐ 5: lambda x: 5 - x[0]
To require x >= 10: lambda x: x[0] - 10

On the contrary, an equality constraint is a function defined so that:

constraint_fun() == 0

For example:

To require x == 10: lambda x: x[0] - 10
In a multivariable problem, to require x + y == 10: lambda x, y: x[0] + y[0] - 10 (in SciPy, the syntac should be: v[0] + v[1])

(You can also have more complex constraints, known as nonlinear constraints.)

Now, please define at least two constraints below. Think about whether you want a bounded interval like [0,5], or something more interesting. These inequality constraints should be added to a constraints list, similar to the example below:

# Example: constraints for 0 <= x <= 5
constraints = [
    {"type": "ineq", "fun": lambda x: x[0]},       # x >= 0
    {"type": "ineq", "fun": lambda x: 5 - x[0]},  # x <= 5
]

Now define your own constraints (you may use the example above or choose different bounds):

constraints = [
    # TODO: define at least 2 inequality constraint dictionaries
]

Once your objective and constraints are defined, run the optimizer using the SLSQP method, which supports constraints.

initial_guess = [1.0]  # try a value inside the feasible region
result = minimize(constrained_objective, initial_guess, method="SLSQP", constraints=constraints)

print(f"Optimal x = {result.x[0]:.6f}")
print(f"Objective = {result.fun:.6f}")
print("Constraint values at optimum:")
for c in constraints:
    print(" ", c["fun"](result.x))

if not result.success:
    print("Message:", result.message)

Now create a quick plot of your objective function over a range of x values, and mark the feasible region and the optimum. You only need a simple plot; the goal is to verify that the solution makes sense visually.

Deliverables

2.1 Define a constrained objective.
2.2 Create at least two constraints.
2.3 Use SciPy’s optimizer with constraints.
2.4 Print and plot the solution.
2.5 How does the constraints change where the minimum occurred?
2.6 Could the unconstrained minimum be infeasible?

Question 3 (2 points)

In real-world optimization, you often have more than one variable. For example, a company might want to choose how many units of different products to make to maximize its profit, subject to resource limits. In this question, you will work with a small real-world production planning problem, and use optimization to find the best solution.

A factory makes two products:

Product A yields $10 profit per unit, requires 2 hours labor and 1 unit material per unit produced,
Product B yields $15 profit per unit, requires 1 hour labor and 3 units material per unit produced.

Available resources:

100 hours labor,
120 units material.

We want to maximize profit but subject to the labor and material limits.

SciPy’s optimizer only minimizes. In order to maximize a function, you want to minimize the negative of the function. For example, to maximize $f(x) = x^2$, you would minimize $-f(x)$.

Start by defining the negative profit function.

def neg_profit(x):
    # x = [units of A, units of B]
    return  # TODO: return negative total profit

Next, define the constraint functions:

Total labor used should be ≤ 100,
Total material used should be ≤ 120,
Units of A and B must be ≥ 0.

Write these as inequality constraints for minimize(…).

constraints = [
    # TODO: labor constraint (ineq)
    # TODO: material constraint (ineq)
    # TODO: non-negativity for product A
    # TODO: non-negativity for product B
]

Now solve the optimizer using the SLSQP method.

initial_guess = [10, 10]  # a feasible guess, feel free to change this
result = minimize(neg_profit, initial_guess, method="SLSQP", constraints=constraints)

A_opt, B_opt = result.x
max_profit = # take the negative of the result.fun, as we minimized negative profit, so the maximum is the negative of that value

print("Optimization successful?", result.success)
print(f"Units of Product A: {A_opt:.4f}")
print(f"Units of Product B: {B_opt:.4f}")
print(f"Maximum profit: ${max_profit:.2f}")

# Print resource usage
labor_used =  # TODO: compute
material_used =  # TODO: compute
print(f"Labor used: {labor_used:.2f} / 100")
print(f"Material used: {material_used:.2f} / 120")

Deliverables

3.1 Define the negative profit objective.
3.2 Add all resource constraints.
3.3 Run the optimizer and print the optimal plan
3.4 Check that constraints are satisfied.
3.5 Why did we minimize negative profit?
3.6 What would happen if we removed the non-negativity constraints?

Question 4 (2 points)

So far you have optimized functions of one and two variables. Now let’s take another 2-variable example but also visualize it. For functions of two variables, contour plots provide an intuitive picture of the optimization landscape. Imagine a topographic map where lines represent levels of equal height.

Below is a template for a 2-variable objective, where one array (v) is passed in, and the function returns a single value.

def objective2d(v):
    x, y = v[0], v[1]
    return  # TODO: a smooth 2D function with a clear minimum, e.g. (1, 2): f = (x-1)^2 + (y-2)^2

Next, build a contour plot of the function over a grid of x and y values. On the contour plot, mark the optimizer’s solution with a red star.

Since many of you probably have not used contour plots or meshgrid, the following code is provided to help you get started.

The np.meshgrid function creates a grid of $x$ and $y$ values, and the objective2d function is called for each point in the grid. The result is a 2D array of function values. np.linspace is used to create the grid of $x$ and $y$ values from $x = -2$ to $x = 4$, and $y = -2$ to $y = 4$. Feel free to change these ranges to something that makes sense for your function.

# Build grid
x = np.linspace(-2, 4, 200)
y = np.linspace(-2, 4, 200)
X, Y = np.meshgrid(x, y)

Z = objective2d([X, Y])

plt.figure(figsize=(10,8))
cs = plt.contour(X, Y, Z, levels=20)
plt.clabel(cs, inline=True, fontsize=8)
plt.scatter(result.x[0], result.x[1], color="red", s=150, marker="*", label="Optimum")
plt.xlabel("x")
plt.ylabel("y")
plt.title("2D Objective Contour Plot")
plt.legend()
plt.show()

Now solve this problem using SciPy. For this question there are no constraints unless you choose to add them.

initial_guess = [1, 1]
result = minimize(objective2d, initial_guess, method="BFGS")

print(f"Optimal x = {result.x[0]:.6f}, y = {result.x[1]:.6f}")
print(f"Optimal value = {result.fun:.6f}")

Deliverables

4.1 Define a 2D objective function.
4.2 Solve it with SciPy.
4.3 Create a contour plot and mark the optimum.
4.4 Short written reflection.
4.5 How does the contour plot help you "see" where the minimum is?
4.6 Does the optimizer’s solution match the contour intuition?

Question 5 (2 points)

Optimization is not just about math, it appears in many areas of data science. One such area is model fitting: finding parameters of a model that make it best match observed data. In this question, you will perform a simple form of curve fitting through optimization. First, we will simulate some noisy data from a linear model. For example:

y = ax + b + noise

The noise can be added using np.random.normal(mean, std, size) which represents size number of random numbers from a normal distribution. In this problem, we will just use a mean of 0 and a standard deviation of 1 (standard normal distribution).

Your job is to estimate the best-fit parameters ($a$) and ($b$) by minimizing the sum of squared errors (SSE).

Below is the data generation code, please read through and run it to get your data.

real_a = # define some real value for a
real_b = # define some real value for b
np.random.seed(42)
x_data = np.linspace(0, 10, 20)
y_true = real_a * x_data + real_b
y_data = y_true + np.random.normal(0, 1, len(x_data))

plt.scatter(x_data, y_data)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Noisy data")
plt.show()

The formula for SSE is:

$SSE = sum((y_{data} - y_{pred})^2)$

Now define the model function and SSE objective. The goal is to find parameters [a, b] that make the line best fit the data in least squares sense.

def model(params, x):
    a, b = params[0], params[1]
    return a * x + b

def sse(params):
    y_pred = # From our input params and x_data, calculate the predicted y values
    return # return the SSE

Estimate the parameters by minimizing the SSE. Use method="BFGS".

initial_guess = # make an initial guess for a and b, that is different from your real_a and real_b
result = # run the optimizer with minimize

a_opt, b_opt = result.x
print(f"Estimated a = {a_opt:.4f}, b = {b_opt:.4f}")
print(f"SSE = {result.fun:.4f}")

Finally, we can plot the fitted line on top of the data.

x_plot = np.linspace(0, 10, 100)
plt.scatter(x_data, y_data, label="Data")
plt.plot(x_plot, model([a_opt, b_opt], x_plot), label=f"Fit: y={a_opt:.2f}x+{b_opt:.2f}")
plt.plot(x_plot, y_true, linestyle="--", alpha=0.6, label="True line")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Curve Fitting via Optimization")
plt.legend()
plt.show()

Deliverables

5.1 Estimate parameters by minimizing SSE.
5.2 Plot the fitted model on top of the data.
5.3 Write a short reflection on the connection to regression
5.4 How does the optimizer’s solution match the true values?
- What would happen if you changed the initial guess?

Submitting your Work

Once you have completed the questions, save your Jupyter notebook. You can then download the notebook and submit it to Gradescope.

Items to submit

firstname_lastname_project8.ipynb

It is necessary to document your work, with comments about each solution. All of your work needs to be your own work, with citations to any source that you used. Please make sure that your work is your own work, and that any outside sources (people, internet pages, generative AI, etc.) are cited properly in the project template.

You must double check your .ipynb after submitting it in gradescope. A very common mistake is to assume that your .ipynb file has been rendered properly and contains your code, markdown, and code output even though it may not.

Please take the time to double check your work. See here for instructions on how to double check this.

You will not receive full credit if your .ipynb file does not contain all of the information you expect it to, or if it does not render properly in Gradescope. Please ask a TA if you need help with this.