Moment based relaxation of polynomials programs
YALMIP comes with a built-in module for polynomial programming using
moment relaxations. This can be used for finding lower bounds on constrained
polynomial programs (inequalities and equalities, element-wise and in the
semidefinite cone), and to extract the corresponding optimizers. The implementation is entirely based on high-level
YALMIP code, and can be somewhat inefficient for large problems (the
inefficiency would then show in the setup of the problem, not actually
solving the semidefinite resulting program). For larger problems, you might
want to check out the dedicated software package
Gloptipoly.
For the underlying theory of moment relaxations, the reader is referred to
[Lasserre].
The following code calculates a lower bound on a concave quadratic
optmization problem.
sdpvar x1 x2 x3
obj = -2*x1+x2-x3;
F = set(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>0);
F = F + set(4-(x1+x2+x3)>0);
F = F + set(6-(3*x2+x3)>0);
F = F + set(x1>0);
F = F + set(2-x1>0);
F = F + set(x2>0);
F = F + set(x3>0);
F = F + set(3-x3>0);
solvemoment(F,obj);
double(obj)
ans =
-6.0000
|
Notice that YALMIP does not recover variables by default, a fact showing up in the
difference between lifted variables and actual nonlinear variables (lifted
variables are the variables used in the semidefinite relaxation to model
nonlinear variables) The lifted variables can be obtained by using the
command relaxdouble . The quadratic
constraint above is satisfied in the lifted variables, but not in the true
variables, as the following code illustrates.
relaxdouble(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24)
ans =
23.2648
double(x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24)
ans =
-2.0000
|
An tighter relaxation can be obtained by using a higher order relaxation (the
lowest possible is used if it is not specified).
solvemoment(F,obj,[],2);
double(obj)
ans =
-5.6593
|
The obtained bound can be used iteratively to improve the bound by adding
dynamically generated cuts.
solvemoment(F+set(obj>double(obj)),obj,[],2);
double(obj)
ans =
-5.3870
solvemoment(F+set(obj>double(obj)),obj,[],2);
double(obj)
ans =
-5.1270
|
The known true minima, -4, is found in the fourth order relaxation.
solvemoment(F,obj,[],4);
double(obj)
ans =
-4.0000
|
The true global minima is however not recovered with the lifted variables, as we can see if we
check the current solution (still violates the nonlinear constraint).
checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint | Type| Primal residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| Element-wise| -0.88573|
| #2| Numeric value| Element-wise| 1.834|
| #3| Numeric value| Element-wise| 5.668|
| #4| Numeric value| Element-wise| 1.834|
| #5| Numeric value| Element-wise| 0.16599|
| #6| Numeric value| Element-wise| 2.0873e-006|
| #7| Numeric value| Element-wise| 0.33198|
| #8| Numeric value| Element-wise| 2.668|
+++++++++++++++++++++++++++++++++++++++++++++++++++
|
Extracting solutions
To extract a globally optimal solution, we need two output
arguments. The first output is a diagnostic structure (standard solution
structure from the semidefinite solver), the second output is the
(hopefully) extracted globally optimal solutions and the third output is
a data structure containing all data needed to extract variables.
[sol,x] = solvemoment(F,obj,[],4);
x{1}
ans =
0.5000
0.0000
3.0000
x{2}
ans =
2.0000
-0.0000
0.0000
|
We can easily check that these are globally optimal solutions since
they reach the lower bound -4 and satisfy the constraints (up to
numerical precision).
setsdpvar([x1;x2;x3],x{1});
double(obj)
ans =
-4.0000
checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Type| Primal residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| Element-wise (quadratic)| -1.1034e-006|
| #2| Numeric value| Element-wise| 0.5|
| #3| Numeric value| Element-wise| 3|
| #4| Numeric value| Element-wise| 0.5|
| #5| Numeric value| Element-wise| 1.5|
| #6| Numeric value| Element-wise| 5.956e-007|
| #7| Numeric value| Element-wise| 3|
| #8| Numeric value| Element-wise| 4.6084e-007|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
Polynomial semidefinite constraints
Nonlinear semidefinite constraints can be
added using exactly the same notation and syntax. The following example is taken from
[D. Henrion, J. B. Lasserre].
sdpvar x1 x2
obj = -x1^2-x2^2;
F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
[sol,x] = solvemoment(F,obj,[],2);
setsdpvar([x1;x2],x{1});
double(obj)
ans =
-4.00003
checkset(F)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Type| Primal residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| LMI (quadratic)| -0.00034633|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
Advanced features
A number of advanced features are available. We just briefly introduce
these here by a quick example where we refine the extracted solution
using a couple of Newton steps on an algebraic systems defining the global
solutions given the optimal moment matrices, and change the numerical tolerance in the extraction
algorithm. Finally, we calculate some different global solutions using
the optimal moment matrices. Please check the moment relaxation example in
yalmipdemo for details.
sdpvar x1 x2
obj = -x1^2-x2^2;
F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
ops = sdpsettings('moment.refine',5','moment.creftol',1e-8);
[sol,data] = solvemoment(F,obj,ops);
xopt1 = extractsolution(data,sdpsettings('moment.refine',0));
xopt2 = extractsolution(data,sdpsettings('moment.refine',1));
xopt3 = extractsolution(data,sdpsettings('moment.refine',10));
xopt4 = extractsolution(data,sdpsettings('moment.creftol',1e-3,'moment.refine',5));
|
The moment relaxation solver can also be called using a more standard
YALMIP notation, by simply defining the solver as 'moment' .
sdpvar x1 x2
obj = -x1^2-x2^2;
F = set([1-4*x1*x2 x1;x1 4-x1^2-x2^2]);
sol = solvesdp(F,obj,sdpsettings('solver','moment','moment.order',2));
setsdpvar(sol.momentdata.x,sol.xoptimal{1});
double(obj)
ans =
-4.00003
checkset(F)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Type| Primal residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Numeric value| LMI (quadratic)| -0.00034633|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|