feasp

Compute solution to given system of LMIs

Syntax

[tmin,xfeas] = feasp(lmisys)

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

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

Description

[tmin,xfeas] = feasp(lmisys) computes a solution, if any exists, of a system of LMIs and returns a vector xfeas of particular values of the decision variables for which all LMIs in the system are satisfied.

Given an LMI system,

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

feasp computes xfeas 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. The LMI constraints are feasible if tmin ≤ 0 and strictly feasible if tmin < 0.

[tmin,xfeas] = feasp(lmisys,options) specifies additional options for the algorithm such as maximum iterations, termination conditions, and constraints on the magnitude of decision variables.

[tmin,xfeas] = feasp(lmisys,options,target) specifies a target value for tmin. The optimization stops when t falls below this value.

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).$$

Calculating the quadratic stability of the polytope of the matrices $$Co\left\{A_1,A_2,A_3\right\}$$ requires finding a solution to this system of LMIs.

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 tmin result -3.1363 means that the problem is feasible. Therefore, the dynamic 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

You can add further constraints on this feasibility problem. For instance, the following command bounds the Frobenius norm of P by 10 while specifying that tmin is 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 a tmin result of -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

Input Arguments

lmisys — LMI system to solve

array

LMI system to solve, specified as an array that you obtain using getlmis.

target — Target value for tmin

0 (default) | scalar

Target value for tmin, specified as a scalar. The optimization terminates when t falls below this target.

options — Control parameters for optimization

[0 0 0 0 0] (default) | five-element vector

Control parameters for optimization, specified as a five-element vector, where:

• options(1) is ignored.

• options(2) sets the maximum number of iterations before the optimization terminates. If options(2) = 0, then feasp defaults to a maximum of 100 iterations.

• options(3) sets the feasibility radius. The feasibility radius R is a constraint on the decision vector x = (x₁, . . ., x_N), such that it lies within the ball given by

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

  In other words, feasp constrains xfeas such that its Euclidean norm does 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.

  If options(3) = 0, the feasibility radius defaults to R = 109. You can change this value as follows:

  • options(3) > 0 — Set the feasibility radius to options(3).

  • options(3) < 0 — Activate flexible bound mode. In this mode, feasp starts with a feasibility radius of 108 and increases this value if necessary during the optimization.

• options(4) sets a tolerance for termination. The optimization terminates if t does not decrease by more than 1 percent (in relative terms) during the last J iterations. When options(4) = 0, J defaults to 10 iterations. Reduce this value to terminate the optimization faster at the possible cost of some accuracy. Increase the value to allow more iterations for convergence at the expense of a possibly longer computation.

• options(5) controls the trace of execution of the optimization procedure. When options(5) = 0, then feasp turns the trace on by default. Set options(5) = 1 to turns the trace off.

Setting option(i) to zero is equivalent to setting the corresponding control parameter to its default value. Consequently, you do not need to redefine the entire vector when changing just one control parameter. For instance, use the following code for options to set the maximum number of iterations to 20, while keeping all other options at their default values.

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

Output Arguments

tmin — Feasibility parameter

scalar

Feasibility parameter, returned as a scalar value. The LMI constraints are feasible if tmin ≤ 0 and strictly feasible if tmin < 0. If tmin is positive and very small, the problem is feasible but not strictly feasible. Some post-analysis might then be required to decide whether xfeas is close enough to feasible.

xfeas — Solution

vector

Solution to the feasibility problem, returned as a vector. If a solution exists, xfeas contains values of the decision variables that satisfy the LMI system. xfeas is in terms of the decision variables, not in terms of the matrix variables of the problem. Use dec2mat to derive feasible values of the matrix variables from xfeas.

Tips

• 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 Tips for details). Since the QR mode typically requires much more memory, MATLAB^® might run out of memory and display the following message.

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

  If you see this message, increase your swap space. If no additional swap space is available, set options(4) = 1. Doing so prevents switching to QR and causes feasp to terminate when Cholesky-based linear algebra fails due to numerical instabilities.

Algorithms

The feasibility solver of feasp is based on Nesterov and Nemirovskii's projective method described in [1] and [2].

References

[1] Nesterov, Yurii, and Arkadii Nemirovskii. Interior-Point Polynomial Algorithms in Convex Programming Society for Industrial and Applied Mathematics, 1994. https://doi.org/10.1137/1.9781611970791

[2] Nemirovskii, A., and P. Gahinet. “The Projective Method for Solving Linear Matrix Inequalities.” In Proceedings of 1994 American Control Conference - ACC ’94, 1:840–44. Baltimore, MD, USA: IEEE, 1994. https://doi.org/10.1109/ACC.1994.751861.

Version History

Introduced before R2006a