GloptiPoly 3: moments, optimization and semidefinite programming

Didier Henrion, Jean Bernard Lasserre, Johan Lofberg

What is GloptiPoly ?

Gloptipoly 3 is intended to solve, or at least approximate, the Generalized Problem of Moments (GPM), an infinite-dimensional optimization problem which can be viewed as an extension of the classical problem of moments . From a theoretical viewpoint, the GPM has developments and impact in various areas of mathematics such as algebra, Fourier analysis, functional analysis, operator theory, probability and statistics, to cite a few. In addition, and despite a rather simple and short formulation, the GPM has a large number of important applications in various fields such as optimization, probability, finance, control, signal processing, chemistry, cristallography, tomography, etc. For an account of various methodologies as well as some of potential applications, the interested reader is referred to and the nice collection of papers .

The present version of GloptiPoly 3 can handle moment problems with polynomial data. Many important applications in e.g. optimization, probability, financial economics and optimal control, can be viewed as particular instances of the GPM, and (possibly after some transformation) of the GPM with polynomial data.

The approach is similar to that used in the former version 2 of GloptiPoly . The software allows to build up a hierarchy of semidefinite programming (SDP), or linear matrix inequality (LMI) relaxations of the GPM, whose associated monotone sequence of optimal values converges to the global optimum. For more details on the approach, the interested reader is referred to .

Installation

GloptiPoly 3 is a freeware subject to the General Public Licence (GPL) policy. It can be downloaded at

www.laas.fr/∼\simhenrion/software/gloptipoly3

The package, available as a compressed archive, consists of several m-files and subdirectories, and it contains no binaries. Extracted files are placed in a gloptipoly3 directory that should be declared in the Matlab working path, using e.g. Matlab’s command

GloptiPoly 3 uses by default the semidefinite programming solver SeDuMi , so this package should be properly installed. Other semidefinite solvers can also be used provided they are installed and interfaced through YALMIP .

Getting started

to run interactively the basic example that follows.

Consider the classical problem of minimizing globally the two-dimensional six-hump camel back function

The function has six local minima, two of them being global minima.

Using GloptiPoly 3, this optimization problem can be modeled as a moment problem as follows:

Once the moment problem is modeled, a semidefinite solver can be used to solve it numerically. Here we use SeDuMi which is assumed to be installed and accessible from the Matlab working path:

The flag status = 1 means that the moment problem is solved successfully and that GloptiPoly can extract two globally optimal solutions reaching the objective function obj = -1.0316.

From version 2 to version 3

The major changes incorporated into GloptiPoly when passing from version 2 to 3 can be summarized as follows:

Use of native polynomial objects and object-oriented programming with specific classes for multivariate polynomials, measures, moments, and corresponding overloaded operators. In contrast with version 2, the Symbolic Toolbox for Matlab (gateway to the Maple kernel) is not required anymore to process polynomial data.

Generalized problems of moments featuring several measures with semialgebraic support constraints and linear moment constraints can be processed and solved. Version 2 was limited to moment problems on a unique measure without moment constraints.

Explicit moment substitutions are carried out to reduce the number of variables and constraints.

The moment problems can be solved numerically with any semidefinite solver, provided it is interfaced through YALMIP. In contrast, version 2 used only the solver SeDuMi.

Solving generalized problems of moments

GloptiPoly 3 uses advanced Matlab features for object-oriented programming and overloaded operators. The user should be familiar with the following basic objects.

A multivariate polynomial is an affine combination of monomials, each monomial depending on a set of variables. Variables can be declared in the Matlab working space as follows:

Variables, monomials and polynomials are defined as objects of class mpol.

All standard Matlab operators have been overloaded for mpol objects:

to delete all existing GloptiPoly variables from the Matlab working space.

2 Measures (meas)

Variables can be associated with real-valued measures, and one variable is associated with only one measure. For GloptiPoly, measures are identified with a label, a positive integer. When starting a GloptiPoly session, the default measure has label 1. By default, all created variables are associated with the current measure. Measures can be handled with the class meas as follows:

to delete all existing GloptiPoly measures from the working space. Note that this does not delete existing GloptiPoly variables.

3 Moments (mom)

Linear combinations of moments of a given measure can be manipulated with the mom class as follows:

The notation I[p]d[k] stands for ∫p dμk\int p\>d\mu_{k} where pp is a polynomial of the variables associated with measure dμkd\mu_{k}, and kk is the measure label.

Note that it makes no sense to define moments over several measures, or nonlinear moment expressions:

Note also the distinction between a constant term and the mass of a measure:

Finally, let us mention three equivalent notations to refer to the mass of a measure:

The first command refers explicitly to the measure, the second command is a handy short-cut to refer to a measure via its variables, and the third command refers to GloptiPoly’s labeling of measures.

4 Support constraints (supcon)

Support constraints are modeled by objects of class supcon. The first command means that variable xx must satisfy x3+2x2−x−2=(x−1)(x+1)(x+2)=0x^{3}+2x^{2}-x-2=(x-1)(x+1)(x+2)=0, i.e. measure dμ1(x)d\mu_{1}(x) must be discrete, a linear combination of three Dirac at 11, −1-1 and −2-2. The second command restricts measure dμ2(y)d\mu_{2}(y) within the unit disk.

Note that it makes no sense to define a support constraint on several measures:

5 Moment constraints (momcon)

We can constrain linearly the moments of several measures:

Moment constraints are modeled by objects of class momcon.

For GloptiPoly an objective function to be minimized or maximized is considered as a particular moment constraint:

The latter syntax is a handy short-cut which directly converts an mpol object into an momcon object.

6 Floating point numbers (double)

Variables in a measure can be assigned numerical values:

which is equivalent to enforcing a discrete support for the measure. Here dμ1d\mu_{1} is set to the Dirac at the point 22.

The double operator converts a measure or its variables into a floating point number:

Discrete measure supports consisting of several points can be specified in an array:

7 Moment SDP problems (msdp)

GloptiPoly 3 can manipulate and solve Generalized Problems of Moments (GPM) as defined in :

where measures dμkd\mu_{k} are supported on basic semialgebraic sets

In the above notations, gik(x)g_{ik}(x), hjk(x)h_{jk}(x) are given real polynomials and bjb_{j} are given real constants. The decision variables in the GPM are measures dμk(x)d\mu_{k}(x), and GloptiPoly 3 allows to optimize over them through their moments

where the αk\alpha_{k} are multi-indices.

8 Solving moment problems msol

Once a moment problem is defined, it can be solved numerically with the instruction msol. In the sequel we give several examples of GPMs handled with GloptiPoly 3.

Following , given a multivariate polynomial g0(x)g_{0}(x), the unconstrained optimization problem

can be formulated as a linear moment optimization problem

In Section 3 we already encountered an example of an unconstrained polynomial optimization solved with GloptiPoly 3. Let us revisit this example:

This indicates that the global minimum is attained with a discrete measure supported on two points. The measure can be constructed from the knowledge of its first moments of degree up to 6:

When converting to floating point numbers with the operator double, it is essential to make the distinction between mpol and mom objects:

The first instruction mmon generates a vector of monomials v of class mpol, so the command double(v) calls the convertor @mpol/double which evaluates a polynomial expression on the discrete support of a measure (here two points). The last command double(mom(v)) calls the convertor @mom/double which returns the value of the moments obtained after solving the moment problem.

Note that when inputing moment problems on a unique measure whose mass is not constrained, GloptiPoly assumes by default that the measure has mass one, i.e. that we are seeking a probability measure. Therefore, if g0 is the polynomial defined previously, the two instructions

are equivalent. See also Section 5.3 for handling masses of measures and Section 5.8.2 for more information on mass constraints.

8.2 Constrained minimization

Consider now the constrained polynomial optimization problem

is a basic semialgebraic set described by given polynomials gi(x)g_{i}(x). Following , this (nonconvex polynomial) problem can be formulated as the (convex linear) moment problem

As an example, consider the non-convex quadratic problem of Section 4.4 in :

Each constraint in this problem is interpreted by GloptiPoly 3 as a support constraint on the measure associated with variable xx, see Section 5.4:

The whole problem can be entered as follows:

Since status=0 the moment SDP problem can be solved but it is impossible to detect global optimality. The value obj=-6.0000 is then a lower bound on the global minimum of the quadratic problem.

The measure associated with the problem variables can be retrieved as follows:

Its vector of moments can be built as follows:

These moments are the decision variables of the SDP problem solved with the above msol command. Their numerical values can be retrieved as follows:

The numerical moment matrix can be obtained using the following commands:

As explained in , we can build a hierarchy of nested moment SDP problems, or relaxations, whose solutions converge monotically and asymptotically to the global optimum, under mild technical assumptions. By default the command msdp builds the relaxation of lowest order, equal to half the degree of the highest degree monomial in the polynomial data. An additional input argument can be specified to build higher order relaxations:

We observe that the moment SDP problems feature an increasing number of variables and constraints. They generate a mononotically increasing sequence of lower bounds on the global optimum, which is eventually reached numerically at the fourth relaxation:

8.3 Rational minimization

Minimization of a rational function can also be formulated as a linear moment problem. Given two polynomials g0(x)g_{0}(x) and h0(x)h_{0}(x), consider the rational optimization problem

is a basic semialgebraic set described by given polynomials gi(x)g_{i}(x). Following , the corresponding moment problem is given by

As an example, consider the one-variable rational minimization problem [4, Ex. 2]:

We can solve this problem with GloptiPoly 3 as follows:

8.4 Several measures

GloptiPoly 3 can handle several measures whose moments are linearly related.

For example, consider the GPM arising when solving polynomial optimal control problems as detailed in . We are seeking two occupation measures dμ1(x,u)d\mu_{1}(x,u) and dμ2(x)d\mu_{2}(x) of a state vector x(t)x(t) and input vector u(t)u(t) whose time variation are governed by the differential equation

Given a polynomial test function g(x)g(x) we can relax the dynamics constraint with the moment constraint

linking linearly moments of dμ1d\mu_{1} and dμ2d\mu_{2}. As explained in , a lower bound on the minimum time achievable by any feedback control law u(x)u(x) is then obtained by minimizing the mass of dμ1d\mu_{1} over all possible measures dμ1d\mu_{1}, dμ2d\mu_{2} satisfying the support and moment constraints. The gap between the lower bound and the exact minimum time is narrowed by enlarging the class of test functions gg.

In the following script we solve this moment problem in the case of a double integrator with state and input constraints:

For the initial condition x0=[1  1]x_{0}=[1\>\>1] the exact minimum time is equal to 3.53.5. In Table 1 we report the monotically increasing sequence of lower bounds obtained by solving moment problems with test functions of increasing degrees. We used the above script and the semidefinite solver SeDuMi 1.1R3.

9 Using YALMIP

By default GloptiPoly 3 uses the semidefinite solver SeDuMi for solving numerically SDP moment problems. It is however possible to use any solver interfaced through YALMIP by setting a configuration flag with the mset command:

Parameters for YALMIP, handled with the YALMIP command sdpsettings, can be forwarded to GloptiPoly 3 with the mset command. For example, the following command tells YALMIP to use the SDPT3 solver (instead of SeDuMi) when solving moment problems with GloptiPoly:

10 SeDuMi parameters settings

The default parameters settings of SeDuMi can be altered as follows:

where pars is a structure of parameters consistent with SeDuMi’s format.

11 Exporting moment SDP problems

A moment problem P of class msdp can be converted into SeDuMi’s input format:

The SDP problem can then be solved with SeDuMi as follows:

See for more information on SeDuMi’s input data format.

Similarly, a moment SDP problem can be converted into YALMIP’s input format:

where variable F contains the LMI constraints (YALMIP class lmi), h is the objective function (YALMIP class sdpvar) and y is the vector of moments (YALMIP class sdpvar). The SDP problem can then be solved with any semidefinite solver interfaced through YALMIP as follows:

12 Moment substitutions

By performing explicit moment substitutions it is often possible to reduce significantly the number of variables and constraints in moment SDP problems. Version 2 of GloptiPoly implemented these substitutions for mixed-integer 0-1 problems only . With version 3, these substitutions can be carried out in full generality.

GloptiPoly 3 carries out moment substitutions as soon as the left hand-side of a support or moment equality constraint consists of an isolated monic monomial. Otherwise, no substitution is achieved and the equality constraint is preserved.

For example, consider the AW29AW_{2}^{9} Max-Cut problem studied in [3, §4.7], with variables xix_{i} taking values −1-1 or +1+1 for i=1,…,9i=1,\ldots,9. These integer constraints can be expressed algebraically as xi2=1x^{2}_{i}=1. The following piece of code builds up the third relaxation of this problem:

We see that out of the 5005 moments (corresponding to all the monomials of 9 variables of degree up to 6), only 465 linearly independent moments appear in a reduced moment matrix of dimension 130.

With the following syntax, moment substitutions are not carried out:

Only the mass is substituted, and the remaining 5004 moments linked by 6435 linear equalities (many of which are redundant) now appear explicitly in a full-size moment matrix of dimension 220.

References