matlab feasp,Compute solution to given system of LMIs

Consider the problem of finding P > I such that:

A1TP+PA1<0,

A2TP+PA2<0,

A3TP+PA3<0,

with data

A1=(-121-3),A2=(-0.81.51.3-2.7),A3=(-1.40.90.7-2.0).

This problem arises when studying the quadratic stability of the polytope of the matrices, Co{A1,A2,A3}.

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 x˙=A(t)x is quadratically stable for A(t)∈Co{A1,A2,A3}.

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 λmax(P) = 8.4653.

P = dec2mat(lmis,xfeas,p);

e = eig(P)

e = 2×1

3.8875

8.4653