feasp

Compute solution to given system of LMIs

Syntax

[tmin,xfeas] = feasp(lmisys,options,target)


Description

[tmin,xfeas] = feasp(lmisys,options,target) computes a solution xfeas (if any) of the system of LMIs described by lmisys. The vector xfeas is a particular value of the decision variables for which all LMIs are satisfied.

Given the LMI system

$$N^{T}LxN\le M^{T}R(x)M,$$

(1)

xfeas is computed by solving the auxiliary convex program:

Minimize t subject to N^TL(x)N–M^TR(x)M≤tI.

The global minimum of this program is the scalar value tmin returned as first output argument by feasp. The LMI constraints are feasible if tmin ≤ 0 and strictly feasible if tmin < 0. If the problem is feasible but not strictly feasible, tmin is positive and very small. Some post-analysis may then be required to decide whether xfeas is close enough to feasible.

The optional argument target sets a target value for tmin. The optimization code terminates as soon as a value of t below this target is reached. The default value is target = 0.

Note that xfeas is a solution in terms of the decision variables and not in terms of the matrix variables of the problem. Use dec2mat to derive feasible values of the matrix variables from xfeas.

Control Parameters

The optional argument options gives access to certain control parameters for the optimization algorithm. This five-entry vector is organized as follows:

• options(1) is not used.

• options(2) sets the maximum number of iterations allowed to be performed by the optimization procedure (100 by default).

• options(3) resets the feasibility radius. Setting options(3) to a value R > 0 further constrains the decision vector x = (x₁, . . ., x_N) to lie within the ball

  $${\displaystyle \sum _{i=1}^N{x}_i^2<R^2}$$

  In other words, the Euclidean norm of xfeas should not exceed R. The feasibility radius is a simple means of controlling the magnitude of solutions. Upon termination, feasp displays the f-radius saturation, that is, the norm of the solution as a percentage of the feasibility radius R.

  The default value is R = 109. Setting options(3) to a negative value activates the “flexible bound” mode. In this mode, the feasibility radius is initially set to 108, and increased if necessary during the course of optimization

• options(4) helps speed up termination. When set to an integer value J > 0, the code terminates if t did not decrease by more than one percent in relative terms during the last J iterations. The default value is 10. This parameter trades off speed vs. accuracy. If set to a small value (< 10), the code terminates quickly but without guarantee of accuracy. On the contrary, a large value results in natural convergence at the expense of a possibly large number of iterations.

• options(5) = 1 turns off the trace of execution of the optimization procedure. Resetting options(5) to zero (default value) turns it back on.

Setting option(i) to zero is equivalent to setting the corresponding control parameter to its default value. Consequently, there is no need to redefine the entire vector when changing just one control parameter. To set the maximum number of iterations to 10, for instance, it suffices to type

options=zeros(1,5)       % default value for all parameters 
options(2)=10

Memory Problems

When the least-squares problem solved at each iteration becomes ill conditioned, the feasp solver switches from Cholesky-based to QR-based linear algebra (see Memory Problems for details). Since the QR mode typically requires much more memory, MATLAB^® may run out of memory and display the message

??? Error using ==> feaslv 
Out of memory. Type HELP MEMORY for your options.

You should then ask your system manager to increase your swap space or, if no additional swap space is available, set options(4) = 1. This will prevent switching to QR and feasp will terminate when Cholesky fails due to numerical instabilities.

Examples

Solve System of LMIs

Consider the problem of finding P > I such that:

$$A_1^{T}P+P{A}_1<0,$$

$$A_2^{T}P+P{A}_2<0,$$

$$A_3^{T}P+P{A}_3<0,$$

with data

$$A_1=\left(\begin{array}{cc}-1& 2\\ 1& -3\end{array}\right),A_2=\left(\begin{array}{cc}-0.8& 1.5\\ 1.3& -2.7\end{array}\right),A_3=\left(\begin{array}{cc}-1.4& 0.9\\ 0.7& -2.0\end{array}\right).$$

This problem arises when studying the quadratic stability of the polytope of the matrices, $$Co\left\{A_1,A_2,A_3\right\}$$.

To assess feasibility using feasp, first enter the LMIs.

setlmis([]) 
p = lmivar(1,[2 1]);

A1 = [-1 2;1 -3];
A2 = [-0.8 1.5; 1.3 -2.7];
A3 = [-1.4 0.9;0.7 -2.0];

lmiterm([1 1 1 p],1,A1,'s');     % LMI #1 
lmiterm([2 1 1 p],1,A2,'s');     % LMI #2 
lmiterm([3 1 1 p],1,A3,'s');     % LMI #3 
lmiterm([-4 1 1 p],1,1);         % LMI #4: P 
lmiterm([4 1 1 0],1);            % LMI #4: I 
lmis = getlmis;

Call feasp to a find a feasible decision vector.

[tmin,xfeas] = feasp(lmis);

 Solver for LMI feasibility problems L(x) < R(x)
    This solver minimizes  t  subject to  L(x) < R(x) + t*I
    The best value of t should be negative for feasibility

 Iteration   :    Best value of t so far 
 
     1                        0.972718
     2                        0.870460
     3                       -3.136305

 Result:  best value of t:    -3.136305
          f-radius saturation:  0.000% of R =  1.00e+09
 

The result tmin = -3.1363 means that the problem is feasible. Therefore, the dynamical system $$\underset{}{\overset{\dot{}}{x}}=A\left(t\right)x$$ is quadratically stable for $$A\left(t\right)\in Co\left\{A_1,A_2,A_3\right\}.$$

To obtain a Lyapunov matrix P proving the quadratic stability, use dec2mat.

P = dec2mat(lmis,xfeas,p)

P = 2×2

  270.8553  126.3999
  126.3999  155.1336

It is possible to add further constraints on this feasibility problem. For instance, the following command bounds the Frobenius norm of P by 10 while asking tmin to be less than or equal to –1.

options = [0,0,10,0,0];
[tmin,xfeas] = feasp(lmis,options,-1);

 Solver for LMI feasibility problems L(x) < R(x)
    This solver minimizes  t  subject to  L(x) < R(x) + t*I
    The best value of t should be negative for feasibility

 Iteration   :    Best value of t so far 
 
     1                        0.988505
     2                        0.872239
     3                       -0.476638
     4                       -0.920574
     5                       -0.920574
***                 new lower bound:    -3.726964
     6                       -1.011130
***                 new lower bound:    -1.602398

 Result:  best value of t:    -1.011130
          f-radius saturation:  91.385% of R =  1.00e+01
 

The third entry of options sets the feasibility radius to 10 while the third argument to feasp, -1, sets the target value for tmin. This constraint yields tmin = -1.011 and a matrix P with largest eigenvalue $${\lambda }_{max}\left(P\right)$$ = 8.4653.

P = dec2mat(lmis,xfeas,p);
e = eig(P)

e = 2×1

    3.8875
    8.4653

References

The feasibility solver feasp is based on Nesterov and Nemirovski's Projective Method described in:

Nesterov, Y., and A. Nemirovski, Interior Point Polynomial Methods in Convex Programming: Theory and Applications, SIAM, Philadelphia, 1994.

Nemirovski, A., and P. Gahinet, “The Projective Method for Solving Linear Matrix Inequalities,” Proc. Amer. Contr. Conf., 1994, Baltimore, Maryland, p. 840–844.

The optimization is performed by the C-MEX file feaslv.mex.

Version History

Introduced before R2006a