Lesson #7 - May 20, 2022

The examples are adapted from:

R. Johansson, "Numerical Python: Scientific Computing and Data Science Applications With Numpy, Scipy and Matplotlib", Apress, 2018

Optimization

Univariate Optimization

Problem: minimize the area of a cylinder with unit volume.

We'll use Brent methd, so we need bracketing values:

Note that we use brent method to fine minima (info here), while we use brentq to find the roots.

Unconstrained Multivariate Optimization

Find the minimum of:

$$f(x_0, x_1) = (x_0-1)^4 +5(x_1-1)^2 - 2x_0x_1$$

The gradient is:

$$ \nabla f = \left[ \begin{array}{c} -2 x_1 + 4(x_0-1)^3 \\ -2x_0 + 10 x_1 -10 \end{array}\right]$$

The Hessian is:

$$ H(f) = \left[ \begin{array}{cc} 12(x_0-1)^2 & -2 \\ -2 & 10 \end{array} \right]$$

We use fmin_ncg which implement Newton's method, info here. Note that the we need an initial guess.

If it is not possible to give values for the gradient and/or the hessia, we can use fmin_bfgs (a quasi-Newton method which estimates numerically the gradient and the Hessian, see here) or fmin_cg (a variant of the steepest descent method, see here). In both methods we need an initial guess.

Sometimes, the fucntion has several local minima, so it is sifficult to chosse an appropriate initial guess, i.e. one for which a global, and non local minimum is found.

Evaluate the function in $x = [3.37653057, 1.48365789]$ and in $x = [-0.5,1.5]$ to verify we did not found the global minimum.

In this case a brute force approach can be usefull in order to find an estimate of a global minimum.

Instead of using directly yhe methods in optimize, we can use a general interface, minimize, selecting the method as an argument.

Note that in this case the optimal point is in attribute x!

Nonlinear least square problems

We saw least square approximation, but law undelying our data depends nonlinearly on some parameters, such as

$$ f(x) = b_0 + b_1 e^{-b_2\,x}$$

we have to use a different approach.

Define the fucntion difference between the data and the law given the values of parameters

We use leastsq in order to estimate the parameters $b_0, b_1, b_2$.

There is as specific tool to fit the values of parameters without having to define the difference function: curve_fit (see here).

Constrained optimization

The variables can be subjected to bounds or constraints.

Find the minimum of

$$f(x_0, x_1) = (x_0 -1)^2 + (x_1-1)^2$$

with the following bound on $x$:

$$x_0 \in [2, 3],\quad x_1 \in [0, 2]$$

The corresponding unconstrained optimization problem gives:

The value of $x_0$ for the miminum is outside the bounds.

One method for constrained optimization is L-BFGS-B (see here), which requires the bounds of the variables.

Slighlty different is the case of problem with constraints.

Find the mimimum of

$$f(x_0, x_1) = (x_0 -1)^2 + (x_1-1)^2$$

such that

$$g(x) = x_1 - 1.75 - (x_0-0.75)^4 \ge 0$$

The constraints are given by a list of dictionaries, where the key is the type of constarint(equality, inequality) and the value is the function defining the constraint. The inequality is always "not less than" type.

The method used isSLSQP (see here).

Linear programming

If both the function to minimize and the constraints are linear, we have a linear programming problem.

Minimize

$$f(x) = -x_0 + 2x_1 -3x_2$$

with the following constraints

$$x_0+x_1 \le 1$$$$-x_0+3x_1 \le 2$$$$-x_1 + x_2 \le 3$$

In this case, the problem and the constraints can be recast in matrix form.

Minimize

$$f(x) = c^T x$$

with the following constraints: $$A x \le b$$

where

$$ x = \left[ \begin{array}{r}x_0\\x_1\\x_2\end{array}\right], c = \left[ \begin{array}{r}-1\\2\\-3\end{array}\right], A = \left[ \begin{array}{ccc} 1 & 1 & 0 \\ -1 & 3 & 0 \\ 0 & -1 & 1 \end{array}\right], b = \left[ \begin{array}{r}1 \\ 2 \\ 3 \end{array}\right]$$

The method used islinprog (see here), which implement a variation of the (classical) simplex method.

Solving ODEs

We want to solve an ordinary differential equation, ODE. Ordinary means that we have only one independent variable, say $x$. The standard explicit form with only a function unkwon $y = y(x)$ is

$$ \frac{d^ny}{dx^n} = f\left( x, y, \frac{dy}{dx}, \ldots, \frac{d^{n-1}y}{dx^{n-1}} \right)$$

The initial condition (the values of $y$, $\frac{dy}{dx}$, $\ldots$, $\frac{d^{n-1}y}{dx^{n-1}}$ at $x=0$) are known, which are named initial value conditions.

An example with a first-order ODE is the following

$$ y'(x) = f(x, y) = x + y(x)^2$$

with initial condition

$$ y(0) = 0 $$

for $$x \in [0, 1.9] $$

Note that while the function $f$ is in closed form, a sample of $y$ at different values of $x$, $x_k$, is returned.

We use the method odeint of module integrate (see here), which implement the LSODA method.

Note that $y'$ gives the slope of the tangent, and therefore in the figure above the solution is always tangent to the arrow field with components $\left( \cos\left(\tan^{-1}f(x,y)\right), \sin\left(\tan^{-1}f(x,y)\right) \right)$, which is called direction field.

An example of a system with two ordinary differential equations is the Lotka-Volterra predator and prey model. Given

$$x(t) = \left[ \begin{array}{r} x_0(t) \\ x_1(t) \end{array} \right], x'(t) = \left[ \begin{array}{r} x'_0(t) \\ x'_1(t) \end{array} \right]$$

we have

$$ x'(t) = f(t,x) = \left\{ \begin{array}{l} ax_0(t) - bx_0x_1(t) \\ cx_0x_1(t) - dx_1(t) \end{array} \right.$$

Where $x_0$ is the number of prey and $x_1$ is the number of predators, $t$ is the independent variable time (which is not directly involved in the definition of $f$!). $a, b, c, d$ are the parameters of the model.

As an example, suppose that at $t=0$ there are 600 preys and 400 predators.

Note that the results is an array where the first dimension is the time and the second dimension the sample of the solution ($x_0$ and $x_1$).

Also in this case is possible to plot the direction field, this time anyway $f$ gives directly the components of the arrows.

The presence of a closed loop means that it is a stable (periodic) solution.

In the case we have a system of ODE of higher order, we can transform it anyway to a system of first order ODES.

For example, in the case of two masses, $m_0$ and $m_1$, connected by springs with stiffness $k_0$ and $k_1$ and dampers with $c_0$ and $c_1$ in free vibration

$$ \left\{ \begin{array}{l} m_0 x_0''(t) +c_0 x_0' + k_0 x_0 - k_1 (x_1-x_0) = 0 \\ m_1 x_1''(t) -c_1 x_1' + k_1 (x_1-x_0) = 0 \end{array}\right.$$

Defining

$y_0 = x_0$, $y_1 = x_0'$, $y_2 = x_1$, $y_3 = x_1'$

we have

$$ \left\{ \begin{array}{l} y_0' = y_1 \\ m_0 y_1' +c_0 y_1 + k_0 y_0 - k_1 (y_2-y_0) = 0 \\ y_2' = y_3 \\ m_1 y_3' -c_1 y_3 + k_1 (y_2-y_0) = 0 \end{array} \right.$$

and finally

$$ \left\{ \begin{array}{l} y_0' = y_1\\ y_1' = -\frac{1}{m_0}\left(c_0 y_1 + k_0 y_0 - k_1 (y_2-y_0) \right)\\ y_2' = y_3\\ y_3' = -\frac{1}{m_1}\left( -c_1 y_3 + k_1 (y_2-y_0) \right)\end{array} \right.$$

In this case we use the general interface to the integrators for ODE, ode, and therefore the function $f$ is defined in slighly different way (the independent variable is the first argument).

We use the same method seen before, LSODA.