MATH06, Univariate optimization
Back to the previous page｜Optimization
List of posts to read before reading this article
Contents
- Method for root finding
- Implementation with sympy for simple case
- Implementation with scipy for more realistic problems
Method for root finding
Bracketing methods
Bisection method
False position (regula falsi)
Interpolation
Iterative methods
Newton’s method (and similar derivative-based methods)
Secant method
Steffensen’s method
Inverse interpolation
Combinations of methods
Brent’s method
Roots of polynomials
Finding one root
Finding roots in pairs
Finding all roots at once
Exclusion and enclosure methods
Real-root isolation
Finding multiple roots of polynomials
Implementation with sympy for simple case
The classic optimization problem: Minimize the area of a cylinder with unit volume
$$the\ objective\ function\ :\ f ([r,h]) = 2πr^{2} + 2πrh$$
$$the\ equality\ constraint\ :\ g([r,h]) = πr^{2}h − 1 = 0$$
import sympy
sympy.init_printing()
r, h = sympy.symbols("r, h")
Area = 2 * sympy.pi * r**2 + 2 * sympy.pi * r * h
Volume = sympy.pi * r**2 * h
h_r = sympy.solve(Volume - 1)[0]
Area_r = Area.subs(h_r)
rsol = sympy.solve(Area_r.diff(r))[0]
rsol
OUTPUT : \(\frac{2^{\frac{2}{3}}}{2 \sqrt[3]{\pi}}\)
Numerical value
_.evalf()
OUTPUT : 0.541926070139289
Verification
Area_r.diff(r, 2).subs(r, rsol)
OUTPUT : 12π
Area_r.subs(r, rsol)
OUTPUT : \(3 \sqrt[3]{2} \sqrt[3]{\pi}\)
_.evalf()
OUTPUT : \(5.53581044593209\)
Implementation with scipy for more realistic problems
$$the\ objective\ function\ :\ f (r) = 2πr^{2} + \frac{2}{r}$$
Method1
from scipy import optimize
import numpy as np
# object function
def f(r):
return 2 * np.pi * r**2 + 2 / r
# optimization
r_min = optimize.brent(f, brack=(0.1, 4))
r_min, f(r_min)
OUTPUT : \(\left ( 0.5419260772557135, \quad 5.535810445932086\right )\)
SUPPLEMENT
brack keyword argument specify a starting interval for the algorithm
Method2
from scipy import optimize
import numpy as np
# object function
def f(r):
return 2 * np.pi * r**2 + 2 / r
# optimization
optimize.minimize_scalar(f, bracket=(0.1, 4))
OUTPUT :
fun: 5.535810445932086
nfev: 19
nit: 15
success: True
x: 0.5419260772557135
Visualization
import matplotlib.pyplot as plt
import numpy as np
# main graph
r = np.linspace(0.1,2,100)
y = 2*np.pi*r**2 + 2/r
plt.plot(r,y)
plt.ylim([0,30])
# optimization point
plt.plot(0.5419260772557135, 5.535810445932086, marker='*', ms=15, mec='r')
plt.annotate("Optimization point", fontsize=14, family="serif", xy=(0.5419260772557135, 5.535810445932086), xycoords="data", xytext=(+20, +50), textcoords="offset points", arrowprops=dict(arrowstyle="->", connectionstyle="arc3, rad=.5"))
plt.show()
List of posts followed by this article
Reference