Examples
Example problem 1
Solve the following nonlinear, nonconvex program with pyipm.py:
minimize (1 + x - x3)2 + (1.61 - [x2 + y2]0.5)2 (1.61 + [x2 + y2]0.5)2
x, y
subject to y - 5⁄3 x + 9⁄20 ≤ 0 and y ≥ 0
where a nonconvex objective function is minimized subject to linear inequality constraints. (Note: a more aesthetically pleasing rendering of the problem definition may be found at the top of the image below.)
To begin with, relevant python packages need to be imported into the script namespace
import numpy as np
import theano
import theano.tensor as T
import pyipm
corresponding to NumPy, Theano, the Theano submodule tensor, and pyipm. The variable T can be used to access the theano.tensor submodule. Then, a Theano vector variable needs to be defined:
x = T.vector('x_dev')
which is assigned to x and named 'x_dev'. All that has been done so far has not been problem specific and it is likely that one will need do these above steps in their own application.
For this problem, the objective function can be written symbolically as
f = (1.0 + x[0] - x[0]**3)**2 + (1.61 - T.sqrt(T.sum(x**2)))**2 * (1.61 + T.sqrt(T.sum(x**2)))**2
where T.sum() and T.sqrt() are the Theano equivalents of the sum and square root functions from NumPy. Note that x[0] corresponds to x and x[1] corresponds to y in the problem statement. The inequality constraints can be defined as follows
ci = T.zeros((2,))
ci = T.set_subtensor(ci[0], 1.6667*x[0] - x[1] - 0.45)
ci = T.set_subtensor(ci[1], x[1])
where T.zeros() constructs an array of zeros much like in NumPy and T.set_subtensor() performs assignment to elements of ci since assignments cannot be done through indexing in Theano.
An initialization for x and y must be provided as the starting point for the optimizer. This can be set as one likes, such as a random initialization or some definite point based on knowledge of the problem. Since pyipm is an infeasible interior-point method, the initialization can be infeasible (i.e. the initialization can violate the inequality constraints) and the solver should be able to reach the feasible region (provided the constraints satisfy a constraint qualification; which they do). To demonstrate this, the infeasible initialization
x0 = np.array([-0.5, -0.25])
was chosen. The solver class, pyipm.IPM, is then instantiated by
problem = pyipm.IPM(x0=x0, x_dev=x, f=f, ci=ci)
and then can be run by calling the class member pyipm.IPM.solve() like so
x, s, lda, fval, kkt = problem.solve()
where x is the solution to the problem, s are the slack variables, lda are the Lagrange multipliers, fval is the objective function value at x, and kkt are the Karush-Kuhn-Tucker (KKT) conditions.
Here is the script as a whole when all the pieces above are put together:
import numpy as np
import theano
import theano.tensor as T
import pyipm
# Theano tensor variable
x = T.vector('x_dev')
# objective function
f = (1.0 + x[0] - x[0]**3)**2 + (1.61 - T.sqrt(T.sum(x**2)))**2 * (1.61 + T.sqrt(T.sum(x**2)))**2
# inequality constraints
ci = T.zeros((2,))
ci = T.set_subtensor(ci[0], 1.6667*x[0] - x[1] - 0.45)
ci = T.set_subtensor(ci[1], x[1])
# variable initialization
x0 = np.array([-0.5, -0.25])
# solver class instantiation
problem = pyipm.IPM(x0=x0, x_dev=x, f=f, ci=ci)
# solve the problem
x, s, lda, fval, kkt = problem.solve()
When this script is run, it will quickly converge to a feasible local minimum (which is the global minimum in this case). The course of the optimization is plotted in the following figure.
