### Optimization Algorithms

There is a Python module for the Steepest Descent algorithm and one for the Newton algorithm.  Since the Newton algorithm uses Eigenmatrices (also called Eigenspaces) for calculating the inverse of the Hessian matrix, I provide a detailed module for Eigenmatrices.

I also provide two modules to demonstrate the use of the Python optimization packages Scipy and GEKKO for solving the optimization problem with two inequality constraints:  max F subject to G1 <= c1 and G2 <= c2, where F is the objective function to be maximized, G1 and G2 are the constraint functions, and c1 and c2 are the constraint constants.  For learning about optimization algorithms, I recommend Antoniou and Lu (2021).

##### Python Code

To receive the Python code for this topic, use the Contact form to request the zip file and I will e-mail it to you.

The optimization algorithm modules use the Sympy, Numpy, Matplotlib, Math, and GEKKO packages as shown in the table below.  You should install these in your Python installation.

##### References

Antoniou, A. and Lu, W-S.  (2021).  Practical Optimization:  Algorithms and Engineering Applications.  Springer.

Steepest Descent

Eigenmatrices

Newton

Inequality Constraints-Scipy

Inequality Constraints-GEKKO ​Purpose

This module implements the Steepest Descent algorithm.  This algorithm requires gradients, and in this module you learn the Python for gradients as well as other matrix and vector calculus methods.

New Python

New Python includes a Python function to encapsulate the calculations for finding the new direction of search; the array() function to create a Numpy array; array functions reshape(), transpose(), and dot(); and the Math package.

Import:  Import the Math package.  We will use the sqrt() function.

Tracing:  For iterative programs like this, it is very informative to use tracing to print intermediate results.  The “tracing = True” statement is near the beginning of the module.

Python function:  Since there are a lot of calculations in the interation to convergence, we take the bulk of that and put it into a function called NewDirection.  The code for a function is indented and the end of a function happens when the code is not indented.  Also note that a function usually requires parameters (here f, g, xk, and ak) and returns results (here ak, dk).

We use Sympy for symbolic derivations and Numpy for numerical calculations.  Looking ahead to Step 2, you can see that we use Sympy to derive the symbolic gradient (Jacobian) from the symbolic function, resulting in a gradient in the form of a Sympy Matrix.  This is passed as one of the parameters to the NewDirection function where we substitute the current x values and then convert the Sympy Matrix to a Numpy array.

Important:  When creating a Numpy array, you must always set the data type using the dtype parameter.  Here we specify the data type as float instead of integer, for example.  Also note that since the mathematics uses column vectors, we need to reshape the Numpy array.

Step 2:  We use the Jacobian() method of a Sympy Matrix since a gradient is a one-dimensional Jacobian.

Step 3:  Initialize a Numpy array to store the column vector points xk on the path.

Step 4:  Optimization algorithms are all iterative searches for a local solution.  To implement the iteration, you use a Python while loop that stops when the algorithm converges.

##### Steepest Descent

Figure produced by Steepest Descent module:

Blue:  Level curves for the Himmelblau function.

Red:   Path taken by the Steepest Descent algorithm to the local minimum.

The Himmelblau function has four minimum values and the minimum found by the algorithm depends on the starting value that you provide. ##### Eigenmatrices

Purpose

The Newton algorithm requires calculation of the inverse of the Hessian, and the most efficient method for this is to use the eigenvalue and eigenvector matrices.  Therefore, in preparation for the next module for the Newton algorithm, I demonstrate how to calculate eigenvalues and eigenvectors and how to multiply matrices.

New Python

New Python includes the hessian(), pprint(), eigh(), and diag() functions.

Step 2:  Use the hessian() method of a Sympy Matrix since a gradient is a one-dimensional Jacobian.  To print a Sympy matrix, use the “pretty print” function pprint().

Step 3:  The Numpy linalg subpackage provides the eig() and eigh() functions for creating eigenvalues and eigenvectors.  Since eigh() is designed for symmetric matrices, this is what we use.  Numpy provides the diag(array) function that creates a matrix using the array values for the diagonal—just what we need for the eigenvalue matrix.

Steps 4-6:  Using the transpose() and dot() functions, use the formula for diagonalization of H.  The output will be in scientific notation unless you suppress it in favour of floating point, which we do on line 50.

Purpose

We implement the Newton algorithm by making some simple changes to the Steepest Descent module.

Code Changes

There is no new Python.  This module uses what you learned in the modules for Steepest Descent and Eigenmatrices.  Here is a description of changes made to the Steepest Descent module.

NewDirection Function:

Most of the changes are in the NewDirection() function.

• Change the parameters (line 23).

• Keep the code for the gradient (up to line 31).

• Replace the remaining code with code from the Eigenmatrices module.

Step 1:  All the same.

Step 2:  Add the code for deriving the Hessian that was used in the Eigenmatrices module.

Step 3:  This is very similar, and differs only by using dk, not adk.

Step 4:  Again, this is very similar and differs only by using dk, not adk.

Step 5:  All the same except for the plot title.

Newton

Figure produced by Newton module:

Blue:  Level curves for the Himmelblau function.

Red:   Path taken by the Newton algorithm to the local minimum.

The Himmelblau function has four minimum values and the minimum found by the algoritm depends on the starting value that you provide. ##### Inequality Constraints - Scipy

Purpose

In this module, you are introduced to the minimize() function in Scipy.  This function provides you with access to unconstrained optimization solvers such as CG (conjugate gradient) and BFGS (Newton) and to constrained optimization solvers such as COBYLA and SLSQP (which we use in this module).

New Python

Import:  Since all we need from Scipy is the minimize function, it makes sense to import only this function.

Step 1:  Since this is a maximization problem, we minimize the negative of the objective function.  Also, since minimize() requires a list of variables, we use x = [x1, x2].

Step 2:  We need to give the algorithm starting values.  The minimize() function requires the starting values in a Numpy array.  A Numpy array is very much like a Python list, but the array uses less memory and is processed faster so it is better for computation.  We define the array using zeros and then assign the values.

Step 3:  Convergence is always improved if we give an algorithm a range of uncertainty, in this case referred to as “bounds.”  Minimize() wants these bounds in a tuple.

Step 4:  Tell minimize() that the constraints are both inequality constraints.

Step 5:  Tell minimize() to use the SLSQP method and to display the iteration details.

##### Inequality Constraints - GEKKO

Purpose

In this module, you are introduced to GEKKO, which provides an interface to some sophisticated optimization algorithms including APOPT, BPOPT, and IPOPT.

New Python

Import:  You import the entire GEKKO package as shown.

Step 1:  Create an object for your model.  In subsequent steps you will attach starting values, bounds, and equations to this model object using the methods Var(), Obj(), and Equation().

Step 2:  The Var() method  provides the starting values and bounds.

Step 3:  The objective function is defined using the Obj() method and the constraints are defined using the Equation() method.

Step 4:  As you can tell so far, this is a very clean and simple way to do things.  There is no need for you to calculate the optimal value of the objective function since that is obtained by the objfcnval option.