These examples illustrate how YALMIP can be used together with the solver VSDP (Verified Semidefinite Programming) to obtain rigorous bounds on achievable objective, and how to use YALMIP, VSDP and INTLAB to compute robust solutions to problems with interval data.

We will first concentrate on the computation of rigorous solutions, i.e. solutions that are verified to be correct, with lower and upper bounds, despite possible errors due to floating point arithmetics. Then we will look at how interval data can be used in YALMIP models.

Make sure that your INTLAB installation is configured and started (you have to run the script startintlab.m to initiate various global variables)

The problem

To illustrate the problem, we use the following demo example from VSDP. This problem is strictly feasible, with optimal objective 1/2.

d = 1e-4;

y = sdpvar(4,1);
Z = [ -y(2)         1/2+0.5*y(1)  -y(3);
      1/2+0.5*y(1)   d            -y(4);
      -y(3)          -y(4)         d];

F   = [Z > 0];
obj = -y(1) - 2*d*y(2);

We begin by solving the problem with SeDuMi.

solvesdp(F,obj,sdpsettings('solver','sedumi'));
format long
double(obj)
ans =
    4.999988671591639e-001

SeDuMi has apparently managed to find a better solution than the optimal!. Clearly something is wrong here. Indeed, if we check feasibility, we see that the solution is infeasible.

>> checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                Type|   Primal residual|   Dual residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   Matrix inequality|      -3.9349e-010|     2.5658e-010|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

This is not too surprising. No SDP solver (except LMILAB which not is recommended for other reasons) is guaranteed to return feasible solutions.

Solving the problem using SDPT3 actually returns a strictly feasible solution, and consequently an objective value larger than 1/2.

solvesdp(F,obj,sdpsettings('solver','sdpt3'));
double(obj)
ans =
    5.000000403019513e-001

checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                Type|   Primal residual|   Dual residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   Matrix inequality|       7.7907e-012|     9.7858e-012|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Based on this result, we know an upper bound on the achievable objective, and we have the corresponding solution.

Verified bounds and feasible solutions

The question is if we can derive both lower and upper bounds, guaranteed to hold also with floating point errors in the computation of the solution, using any SDP solver. This is where VSDP comes into play.

In the default setting, if a user tells YALMIP to use VSDP, then YALMIP will use VSDP to compute a rigorously feasible solution, thus leading to, in YALMIPs notation, an upper bound on the achievable performance.

solvesdp(F,obj,sdpsettings('solver','vsdp'));
double(obj)
    5.000018188263583e-001

checkset(F)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|      Constraint|                Type|   Primal residual|   Dual residual|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Numeric value|   Matrix inequality|       1.9726e-010|     2.5658e-010|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Notice that (assuming you have a standard installation of YALMIP) SeDuMi will be called from VSDP. However, despite the fact that SeDuMi, as we saw above, returns an infeasible solution, VSDP manage to create a feasible solution, and the objective returned is a verified upper bound on the achievable objective, also under floating point errors.

To get a lower bound on the achievable objective, an option has to be turned on, and additional outputs have to be generated.

Note that in YALMIP everything is minimization in standard setting, while the syntax in VSDP is the other way around. Hence, to get a lower bound on our minimization problem, we need to turn on upper bound computation in VSDP. By default, only vsdp.verifiedlower is turned on.

options = sdpsettings('solver','vsdp','vsdp.verifiedupper',1,'savesolveroutput',1);
sol = solvesdp(F,obj,options);
-sol.solveroutput.fU
ans =
    4.999988953605884e-001

This is a verified lower bound for our problem. Our optimal value is somewhere between 4.999988953605884e-001 and 5.000018188263583e-001, consistent with the known optimal value.

Robust semidefinite programming

YALMIP supports exact robust semidefinite programming for problems with interval data via the robust optimization framework. However, exact solutions typically require very large models, since they lead to enumeration in the SDP case.

Almost for free, the interface to VSDP enables approximate robust semidefinite programming for problems with interval data in YALMIP.

Let us consider a robust stability analysis problem. We want to show that the matrix A satisfies a Lyapunov inequality, and find the minimal trace solution

A = [-1 2;1 -3];
P = sdpvar(2,2);

F = [A'*P + P*A < -eye(2)];
obj = trace(P);

solvesdp(F,obj)

However, let us now assume that every element of A is perturbed by a small amount. A robust solution is easily constructed in YALMIP by using the uncertain declaration.

W = sdpvar(2,2,'full');
A = [-1 2;1 -3] + W;

F = [A'*P + P*A < -eye(2)];
F = [F, uncertain(W)];
F = [F, -0.1 < W < 0.1];
obj = trace(P);

solvesdp(F,obj)

Although it works well here, one has to understand that the complexity of robustifying a semidefinite constraint with interval data is exponential. One new constraint has to be introduced for every vertex of the uncertainty polytope. In this case, 4 constraints are generated by YALMIPs robust optimization module.

An alternative is to describe the data using the interval data class in INTLAB, and solve the problem using VSDP. YALMIP handles interval numbers just as standard doubles, hence the following code would compute a solution which is feasible for all instances in the interval data.

W = midrad(zeros(2),0.1*ones(2))

A = [-1 2;1 -3] + W;

F = [A'*P + P*A < -eye(2)];
obj = trace(P);

solvesdp(F,obj)

Notice that the objective value here is larger than when we solved it using the standard robust optimization framework. The reason is that VSDP only solves an outer approximation of the uncertain problem.

Comments

Verified solutions and interval data is (should be...) entirely integrated in the infrastructure of YALMIP. Hence, you can use interval data in problems with convexity-aware nonlinear operators or sum-of-squares models, etc.

As an illustration, the following problem solves a robust sum-of-squares problem.

sdpvar x t s

% Uncertain linear term
p = x^4 + s*x + 1

% Uncertainty
F = [uncertain(s), -0.5 < s < 0.5];

% Lower bound SOS
F = [F, sos(p - t)];

% Solve robust problem
solvesos(F,-t)
double(t)
 

The interval-data version is implemented almost in the same way. Note once again that the solution is conservative.

sdpvar x t

% Uncertainty
s = midrad(0,0.5);

% Uncertain linear term
p = x^4 + s*x + 1

% Lower bound SOS
F = [sos(p - t)];

% Solve robust problem
solvesos(F,-t)
double(t)