An important use of calculus is the solving of maximum-minimum problems, that is, fidning the absolute maximum of minimum value of some varying quantity and that point at which that maximum of minimum occurs.

There’s an extensive treatment of optimization on Scipy’s lecture notes.

The general strategy for solving these problems involves translating the problem into an equation in one variable. Then one can use derivatives to find out critical points and evaluate whether these points are maximum or minimum values over a (closed or open) interval.

Let’s take a textbook example. A stereo manufacturer determines that in order to sell \(x\) units of a new stereo, the price per unit, in dollars, must be

\[p(x) = 1000 - x\] The manufacturer also determines that the total cost of producing \(x\) units is given by

\[C(x) = 3000 + 20x\] (a) find the total revenue (b) find the total profit (c) how many units must the company produce and sell to maximize profit? (d) what is the maximum profit? (e) what price per unit gives the maximum profit?

The total revenue is the price times the units sold.

from sympy import *

x = symbols("x")
revenue = expand((1000 - x) * x)
revenue
## -x**2 + 1000*x

\[R(x) = 1000x - x^2\]

The total profit is the total revenue minus the cost.

cost = 3000 + (20 * x)
profit = expand(revenue - cost)
profit
## -x**2 + 980*x - 3000

\[P(x) = -x^2 + 980x - 3000\]

To find how many units are needed to maximize profits, we need to find the first derivative and find its critical values (where the derivative is zero), we can use second derivative test to figure out if we have a minimum or a maximum.

The rule is that when a function’s slope is zero at \(x\), and the second derivative at \(x\) is:

  • less than 0, it is a local maximum
  • greater than 0, it is a local minimum
  • equal to 0, then the test fails (there may be other ways of finding out though)
deriv = diff(profit, x)
deriv
## 980 - 2*x
solve(deriv, x)
## [490]
dderiv = diff(profit, x, 2)
dderiv
## -2

\[P^\prime(x) = -2x + 980\] \[P^\prime(0) = 490\] \[P^{\prime\prime}(x) = -2\] The maximum profit is given by

profit.evalf(subs={x:490})
## 237100.000000000

The price per unit needed to make the maximum profit is

price = 1000 - x
ppunit = price.evalf(subs={x:490})
ppunit
## 510.000000000000

Scipy methods

scipy has a handful of methods to optimize functions.

As an example, suppose you derived that your profit function is:

\[P(x) = -x^2 + 980x + 20x\]

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.optimize import fmin


def f(x):
    return -x ** 2 + 980 * x - 3000


# fmin finds the minimum, we flip the sign to find the max
max_x = fmin(lambda x: -f(x), 0) 
## Optimization terminated successfully.
##          Current function value: -237100.000000
##          Iterations: 43
##          Function evaluations: 86
print(f"Maximum value occurs at {round(max_x.item(), 2)} units sold")
## Maximum value occurs at 490.0 units sold
print(f"Maximum value is {round(f(max_x.item()), 2)} USD")
## Maximum value is 237100.0 USD
x_values = np.linspace(0, 1000)
y_values = f(x_values)

fig, ax = plt.subplots()
ax.plot(x_values, y_values)
ax.plot(max_x.item(), f(max_x.item()), color="red", marker="x")
plt.title("profit maximum")