- Generated on for SCIP Doxygen Documentation by
1.15.0
To find globally optimal solutions for (mixed-integer) nonlinear optimization problems, SCIP relies on a linear relaxations to estimate the gap between the objective function value of a best known feasible solution and the optimal value of the problem. To construct such a relaxation, linear under- or overestimators of nonlinear functions are often used. In the following, we illustrate how to extend the techniques available in SCIP by using the nonlinear handler plugin. This plugin type allows close interaction with the code that constructs the linear relaxation and performs variable bound tightening for nonlinear constraints.
As an exemplary application from computational chemistry, consider the solution of the Lennard-Jones Cluster problem. Given a fixed number of particles, the task is to position them in 3-dimensional space such that the sum of the Lennard-Jones potentials between pairs of particles is minimized. The Lennard-Jones potential combines a repulsive and an attractive potential that depend on the Euclidean distance of two particles. A widely used formula is
\[4\epsilon\left[\left(\frac{\sigma}{d}\right)^{12} - \left(\frac{\sigma}{d}\right)^{6}\right],\]
where \(d\) is the distance between the particles, and \(\epsilon\), \(\sigma>0\) are parameters. The minimization of the Lennard-Jones potential for \(\epsilon=\sigma=1\) is a popular test problem for global optimization.
The problem for \(N\) particles can be formulated in SCIP as follows:
\begin{equation}\begin{aligned} \min\; & 4\sum_{i=1}^N \sum_{j=i+1}^N p_{ij} \\ \text{s.t.}\; & p_{ij} \geq (\Vert x^i - x^j\Vert_2^2)^{-6} - (\Vert x^i - x^j\Vert_2^2)^{-3} && \forall i,j\in\{1,\ldots,N\}, i < j, \\ & x^i \in [-9,9]^3 && i\in\{1,\ldots,N\}. \end{aligned} \end{equation}
Bounding the particles' positions by the box \([-9,9]^3\) is an artificial constraint that is required to ensure convergence of SCIP's algorithm.
To solve this problem, SCIP needs to compute linear underestimators for the function \(f(r) := r^{-6} - r^{-3}\), where \(r\) is a variable that SCIP introduces to stand for the distance of two particles ( \(\Vert x^i - x^j\Vert_2^2\)). In the standard approach, \(r^{-6}\) is recognized as convex and \(-r^{-3}\) as concave. A convex underestimator for \(f(r)\) is then derived by replacing \(-r^{-3}\) by the secant on the function between current lower and upper bounds on \(r\). While this gives a valid underestimator, it is far from being tight, which results in a weak dual bound.
A tighter relaxation could be obtained if \(r^{-6}\) and \(-r^{-3}\) were not underestimated separately, but the convex envelope of \(f(r)\) were used instead. This is achieved by means of the nonlinear handler in this example. It implements the following central callbacks:
Given an algebraic expression, check whether it is of the form \(r^{-6} - r^{-3} - p\). Here, \(r\) and \(p\) could be any expressions, but will correspond to a quadratic expression \(\Vert x^i - x^j\Vert_2^2\) and a variable expression \(p_{ij}\), respectively, for some particles \(i,j\) in this example. If an expression of this structure is found, the callback informs the handler for nonlinear constraints that it will provide underestimators and participate in bound tightening. To produce underestimators that can be added to the LP relaxation, the nonlinear handler requires variables that stand for \(r\) and \(p\). SCIP ensures that these are made available.
In the actual implementation, the handler checks for \(a r^{-6} - a r^{-3} + b p\) for some coefficients \(a\in\{-1,1\}\) and \(b\neq 0\), because SCIP's handler for nonlinear constraints may multiply \(p_{ij} \geq r_{ij}^{-6} - r_{ij}^{-3}\) by \(-1\). If \(a<0\), the handler will provide overestimators instead of underestimators. For the remainder of this documentation, we will assume \(a=1\) and \(b=-1\), though.
At a node of the branch-and-bound tree, let \((r',p')\) be the value of the expressions \(r\) and \(p\) at the current node and \(\ell\) and \(u\) be the bounds on the expression \(r\) (SCIP computes these via interval arithmetics). The nonlinear handler needs to compute a linear underestimator for \(f(r)-p\) that is as tight as possible at \((r',p')\). Note that \(f(r)\) attains its minimum at \(r_{\min}=\sqrt[3]{2}\approx 1.2599\) and has an inflection point at \(r_{\inf}=\sqrt[3]{\frac{7}{2}}\approx 1.51829\). The function is convex for \(0<r\leq r_{\inf}\) and concave for \(r\geq r_{\inf}\).
\[ f'(\tilde r) = \frac{f(u) - f(\tilde r)}{u - \tilde r}. \]
This equation is solved by using Newton's Method to find a root of \(f(u) - f(\tilde r) - f'(\tilde r) ( u - \tilde r)\), using \(\frac{1}{2}(r_{\min}+r_{\inf})\) as starting pointing. If \(\tilde r < \ell\), then the secant between \(\ell\) and \(u\) can be used as linear underestimator.Given bounds on \(r\) and \(p\), compute the range of \(f(r)-p\). The minimum of \(f(r)\) is attained either at \(\sqrt[3]{2}\), the lower, or the upper bound on \(r\). The maximum of \(f(r)\) is attained at the lower or upper bound on \(r\), if \(r>0\).
Given an interval \([\underline{g},\overline{g}]\) for \(f(r)-p\) and bounds \([\underline{p},\overline{p}]\) for \(p\), compute bounds on \(r\). The code first considers the equation \(z^2 - z\in [\underline{g},\overline{g}] + [\underline{p},\overline{p}]\) and finds an interval \([\underline{z},\overline{z}]\) such that \([\underline{z},\overline{z}]^2-[\underline{z},\overline{z}]\supseteq [\underline{g},\overline{g}] + [\underline{p},\overline{p}]\). Next, bounds on \(r\) are given by \([\underline{z},\overline{z}]^{-\frac{1}{3}}\).
See the Install file
1.15.0