Optimization Problems

Today, we want to review our work with using derivatives to solve optimization problems and extend this to solving problems involving more than one variable; introducing the language and notation of partial derivatives.

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
from mpl_toolkits import mplot3d

PROBLEM 1: Economics and Profit

Starting from the relationship \(Profit = Revenue - Cost\); let’s consider an example.

• \(C(x) = \text{cost of advertisements} = 900 + 400x - x^2\)

900 dollars to get setup

print cost of 400\(x\)

some volume savings \(x^2\)

• \(R(x) = \text{revenue due to x advertisements} = 600x - 6x^2\)

sales 600 per ad

subtracting 6\(x^2\) for diminishing returns

• profit = revenue - cost: \(P(x) = R(x) - C(x)\)

To find max: \(P'(x) = 0\)!

PROBLEM

A wire four feet long is cut in two pieces. One piece forms a circle of radius \(r\), the other forms a square of side \(x\). Choose \(r\) to minimize the sum of their areas. Then choose \(r\) to maximize.

PROBLEM

A fixed wall makes one side of a rectangle. We have 200 feet of fence for the other three sides. Maximize the area A in 4 steps:

1. Draw a picture of the situation.
2. Select one unknown quantity as \(x\) (but not A!).
3. Find all other quantities in terms of \(x\).
4. Solve \(dA/dx =0\) and check endpoints.

A Second Approach

THE SECOND DERIVATIVE TEST

\[f''(x) > 0 \quad g''(x) < 0\]

def f(x): return x**2 - 1
def g(x): return -x**2 + 1

x = np.linspace(-2,2,100)
plt.figure()
plt.plot(x, f(x), label = '$f(x)$')
plt.plot(x, g(x), label = '$g(x)$')
plt.axhline(color = 'black')
plt.legend(loc = 'upper left')
plt.title('$f(x) = x^2 - 1$\n$g(x) = -x^2 + 1$', loc = 'left')
plt.savefig('images/maxsec.png')

Extending the work

\[f(x) = ax^3 + b^2x^2 + srx + 12b\]

• \(\frac{df}{da} = x^3\)
• \(\frac{df}{dx} = 3ax^2 + 2b^2x + sr\)
• \(\frac{df}{ds} = rx\)

Curves vs. Surfaces

• As functions we have:

\[y = f(x) \quad vs. \quad z = f(x,y)\]

• With first derivatives we have:

\[\frac{df}{dx} \quad \quad \frac{\partial f}{\partial x} \text{and} \frac{\partial f}{\partial y}\]

Derivative of \(f\) with respect to \(x\) as opposed to the partial derivative of f with respect to \(x\) or \(y\).

A Picture of 3D

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
#plt.savefig('images/3d_axes.png')

#example function z
def f(x, y):
    return (x ** 2 + y ** 2)

#now, we define both our x and y domain and range
x = np.linspace(-6, 6, 30)
y = np.linspace(-6, 6, 30)

#meshgrid creates a single object of all our x, y pairs
#and then we can substitute this directly into our function
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
ax.contour3D(X, Y, Z, 50, cmap='binary')
plt.savefig('images/surface3d.png')

ax.view_init(70, 35)
fig

Function in 3D

A familiar version of our quadratic in 3d would be:

\[f(x, y) = x^2 + y^2\]

Slicing Dimensions

We can consider slices of our surface to return our 2D perspective. This is essentially the meaning of a partial derivative.

fig = plt.figure(figsize=plt.figaspect(0.5))

ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
cset = ax.contour(X, Y, Z, zdir='z')
cset = ax.contour(X, Y, Z, zdir='x')
cset = ax.contour(X, Y, Z, zdir='y')
ax.plot_surface(X, Y, f(X, Y), rstride=8, cstride=8, alpha=0.1)

ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
cset = ax.contour(X, Y, Z, zdir='z', offset = 0)
cset = ax.contour(X, Y, Z, zdir='x', offset = -6)
cset = ax.contour(X, Y, Z, zdir='y', offset = 6)
ax.plot_surface(X, Y, f(X, Y), rstride=8, cstride=8, alpha=0.1)
#plt.savefig('images/3dslicin.png')

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x12197e080>

Analytic Solutions

\[f(x,y) = x^2 + y^2\]

\[\frac{\partial f}{\partial x} = 2x \quad \rightarrow \quad 2x = 0\]

\[\frac{\partial f}{\partial y} = 2y \quad \rightarrow \quad 2y = 0\]

Seems our minimum is located at the point \((0, 0, f(0, 0))\).

PROBLEMS

1. Consider \(f(x,y) = 11x^2 -2xy + 2y^2 + 3y\). Find its critical point, and discuss whether this is a maximum or a minimum value.
2. Consider \(f(x,y) = \displaystyle \frac{x - y}{2x^2 + 8y^2 + 3}\). Use SymPy to find the partial derivatives and critical points. Argue based on the plot whether these are maximum or minimum values.

plt.figure()
plt.contour(X, Y, f(X, Y))
plt.plot(0, 0, 'x', markersize = 12, label = 'minimum')
plt.legend(frameon = 'False')
#plt.savefig('images/topdown.png')

<matplotlib.legend.Legend at 0x121a11f28>

The Second Derivative in 3D

\(f_{xx}, f_{yy},\) and \(f_{xy}\)

• We put these together using the idea of a discriminant, definied as follows:

\[D(a,b) = f_{xx}(a,b)f_{yy}(a,b) - f_{xy}(a,b)^2\]

and we will determine maximums and minimums as follows:

• If \(D > 0\); local max or min
• If \(D < 0\); saddle point
• If \(D = 0\); inconclusive

Example

\[f(x,y) = (x^2 + y^2)e^{-x}\]

import sympy as sy

x, y = sy.symbols('x y')

def f(x, y): return (x**2 + y**2)*sy.exp(-x)

fx = sy.diff(f(x, y), x)
fy = sy.diff(f(x, y), y)

eq1 = sy.Eq(fx, 0)
eq2 = sy.Eq(fy, 0)

sy.solve(eq1)

[{x: -sqrt(-y**2 + 1) + 1}, {x: sqrt(-y**2 + 1) + 1}]

sy.solve(eq2)

[{y: 0}]

sy.solve(eq1)[0][x].subs(y, 0)

0

sy.solve(eq1)[1][x].subs(y, 0)

2

fxx = sy.diff(fx, x)
fyy = sy.diff(fy, y)
fxy = sy.diff(fx, y)

D = fxx * fyy - fxy**2

sy.pprint(sy.factor(D))

  ⎛ 2          2    ⎞  -2⋅x
2⋅⎝x  - 4⋅x - y  + 2⎠⋅ℯ

D.subs(x, 0).subs(y, 0)

4

D.subs(x, 2).subs(y, 0)

-4*exp(-4)

def f(x, y):
    return (x**2 + y**2)*np.exp(-x)

#now, we define both our x and y domain and range
x = np.linspace(-1, 3, 30)
y = np.linspace(-2,2, 30)

#meshgrid creates a single object of all our x, y pairs
#and then we can substitute this directly into our function
X, Y = np.meshgrid(x, y)

#set up the figure
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');

#plot the surface
#ax.plot_surface(X, Y, f(X, Y), cmap = 'viridis')

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z'); ax.plot_surface(X, Y, f(X, Y), cmap = 'viridis')
ax.view_init(2, 85)
plt.savefig('images/soln3d.png')

for i in range(2):
    for j in range(2):
        print(i, j, i**2 + j**2)

0 0 0
0 1 1
1 0 1
1 1 2

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');


for i in range(10):
    for j in range(10):
        #print(i, j, i**2 + j**2)
        ax.scatter3D(i, j, i**2 + j**2)

plt.savefig('images/nestloop.png')

The Nested Loop

for i in range(10):
    for j in range(10):
        #print(i, j, i**2 + j**2)
        ax.scatter3D(i, j, i**2 + j**2)