Rank constrained semidefinite programming problems

This example requires the solver LMIRANK (and a semidefinite solver). That solver is practically obsolete so the coding below is at risk of being outdated.

Rank constrained problems

A lot of problems, in particular in control theory, can be written using rank constraints on symmetric postive semidefinite matrices. Unfortunately, adding a simple rank constraint to a matrix in a standard SDP problem leads to an intractable optimization problem. Nevertheless, with a reasonable initial guess, a local search can be efficient for finding a low-rank solution in some cases. The solver LMIRANK performs such a local search and is interfaced by YALMIP. Details on the algorithm in LMIRank can be found in Orsi 2006.

A common alternative approach to rank-constrained problems is to use various relaxations. On such example is to penalize the nuclear norm of the matrix, which is readily available in the norm operator.

Dynamic controller design using rank constraints

To illustrate rank constraints in YALMIP, and how to use the solver LMIRANK, we solve one of the examples in Orsi 2006. The problem is to design a dynamic output feedback controller for the system \(\dot{x} = Ax+Bu, y=Cx\).

Define data for the dynamical system.

A = [0.2 0 1 0;0 0.2 0 1;-1 1 0.2 0;1 -1 0 0.2];
B = [0 0 1 0]';
C = [0 1 0 0];

Define matrices X and Y (essentially modeling a Lyapunov matrix and its inverse)

X = sdpvar(4,4);
Y = sdpvar(4,4);

Without going into the theoretical background, the following LMI constraints are necessary for the existence of a controller.

Bp = null(B')';
Cp = null(C)';
W  = eye(3)*1e-6;
F = [Bp*(A*X+X*A')*Bp' <= -W, Cp*(Y*A+A'*Y)*Cp' <= -W];

The rank constraint enter when we want to couple the matrices X and Y, and ensure the existence of a dynamic controller with at most 2 states. Once again stating the equations without any theoretical justification, the resulting constraints are

F = [F, [X eye(4);eye(4) Y] >= 0];   % not needed, see below
F = [F, rank([X eye(4);eye(4) Y]) <= 4+2];

In the current implementation in YALMIP, a rank constraint automatically appends a positive semidefiniteness constraint on the involved matrix. Hence, the first constraint in the previous piece of code can (and should) be removed. At the moment, the rank operator is implemented rather straightforwardly using a nonlinear operator, but requires some dedicated code internally (read hack) to work when actually calling a solver. Hence, it should be emphasized that the current implementation is experimental and can be subject to changes in order to make the operator better integrated.

At this point, we are ready to solve the problem. The solver LMIRANK needs an initial guess, and this can be supplied either by the user (by initializing variables and using the ‘usex0’ option in standard fashion) or automatically by YALMIP. YALMIP computes an initial guess by removing the rank constraint and then minimizes the trace of the rank constrained matrix (or sum of traces of all rank constrained matrices).

Note that LMIRANK does not support objective functions, but only solves feasibility problems. If you want to optimize an objective function you need to, e.g., perform a bisection.

The following call will automatically select an SDP solver to find the initial guess, solve the initial problem, and then call LMIRANK (if installed) to search for a low-rank solution.

optimize(F);

The following call ensures that the solver for the initial guess is SeDuMi.

optimize(F,[],sdpsettings('lmirank.solver','sedumi'))

Note that these computations only give a feasible pair X and Y. To actually compute the controller (called reconstruction) requires some additional steps.

If we compute the eigenvalues of our rank constrained matrix, we see that we effectively have a rank of 6.

eig(value([X eye(4);eye(4) Y]))
ans =
   -0.0000
   -0.0000
    1.6780
    2.5342
    2.7754
    2.9674
    6.2289
    6.8932

indeed, this is the rank reported by YALMIP.

rank([X eye(4);eye(4) Y])
Linear scalar (real, derived, 1 variables, current value : 6)

Note that rank computations are inherently hard from a numerical point of view, hence it may easily happen that the reported rank is higher than the desired rank, even though LMIRank claimed feasibility. In such cases, a more detailed analysis using tolerances might be necessary (YALMIP always uses the default tolerance.)

rank(value([X eye(4);eye(4) Y]))
ans =
     6
rank(value([X eye(4);eye(4) Y]),1e-6)
ans =
     6