Rank constrained LMIs

This example requires the solver LMIRank (and a semidefinite solver)

A lot of problems, in particular in control theory, can be addressed using rank constraints. Unfortunately, adding a simple rank constraint to a matrix in a standard LMI problem leads to an NP-hard 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 LMIRank can be found in [Orsi et. al.].

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 et. al.]. The problem is to design a dynamic output feedback controller for the system x' = Ax+Bu, y=Cx

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 = set(Bp*(A*X+X*A')*Bp' < -W) + set(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 + set([X eye(4);eye(4) Y] >= 0);   % not needed, see below
F = F + set(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 (i.e. a non-strict constraint, so these appended cones are not affected by the option shift.) Hence, the first constraint in the previous piece of code can (and should) be removed.  Also note the use of a non-strict rank constraint. A strict inequality can be used, but will currently be interpreted as non-strict, hence it should not be used in order to avoid confusion. 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 minimizing 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.

solvesdp(F);

The following call ensures that the solver for the initial guess is SeDuMi, and uses an accuracy recommended by the author of LMIRank.

solvesdp(F,[],sdpsettings('lmirank.solver','sedumi','sedumi.eps',0))

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 immediately see that we effectively have a rank of 6.

eig(double([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(double([X eye(4);eye(4) Y]))

ans =
     6

rank(double([X eye(4);eye(4) Y]),1e-6)

ans =
     6