Common Lyapunov matrix from T–S fuzzy LMI differs from published result (scaling issue?)

I am implementing a Takagi–Sugeno (T–S) fuzzy model and computing a common Lyapunov matrix using the MATLAB LMI Toolbox.

The system matrix contains one nonlinear term Z1​, which I replace by its minimum and maximum values to generate the fuzzy rules.

%System matrix
syms Z1
A_general = [
    -2200/79, -2000/79;
    9223372036854775808/4620909390464243, ...
    (381128907510930132896351001478811136*Z1)/ ...
    125959558087561075514773639947125
];
%Premise variable bounds
X_MG10(2) =198.3325; 
x2min = -89.4;   % Last point at which LMI is feasible
x2max =  89.4;

f1max = 1/(x2min + X_MG10(2));
f1min = 1/(x2max + X_MG10(2));

Zmin = f1min;
Zmax = f1max;
%T–S fuzzy rule generation
premiseIdx = 2;
numPremise = length(premiseIdx);
numRules = 2^numPremise;

A_rule = cell(numRules,1);
ruleIndices = dec2bin(0:numRules-1) - '0';

for iRule = 1:numRules
    A_local = A_general;
    if ruleIndices(iRule) == 0
        A_local = subs(A_local, Z1, Zmin);
    else
        A_local = subs(A_local, Z1, Zmax);
    end
    A_rule{iRule} = double(A_local);
end
%Common Lyapunov LMI formulation
N = numel(A_rule);
n = size(A_rule{1},1);

setlmis([])

p = lmivar(1,[n 1]);

for k = 1:N
    lmiterm([k 1 1 p],1,A_rule{k},'s');
end

pp = 1e-6;
S = newlmi;
lmiterm([-S 1 1 p],1,1);

lmiterm([ S 1 1 0], -pp*eye(n));

lmis = getlmis;
[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.056084
     2                        0.012870
     3                        0.012870
     4                        0.012870
     5                    2.404692e-03
     6                    3.445975e-05
     7                       -0.032602

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

P = dec2mat(lmis,xfeas,p)

P = 2×2

1.0e+03 * 3.5077 0.0488 0.0488 0.0445

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Result reported in the reference paper

The paper reports the following common Lyapunov matrix:

P =[103.1 1.43;

1.43 1.31]

My questions

1. Is the difference between the two matrices simply due to scaling / non-uniqueness of LMI feasibility solutions?
2. Should I add normalization or optimization constraints (e.g., trace(P) minimization or I≤P≤ aI to obtain a Lyapunov matrix comparable to the published result?

Any clarification or guidance would be appreciated.

Answer 1

[Sam Chak]

Hi @Manish Kumar

Q1: Is the difference between the two matrices simply due to scaling / non-uniqueness of LMI feasibility solutions?

The positive-definite matrix "P" computed by MATLAB's LMI solver differs significantly from the "P" referenced in the paper. The LMI solver should consistently produce the same result for identical inputs and matrices. While we do not know exactly how the authors generated their "P", I believe they may have used truncated values and chose not to display lengthy numbers, such as 9223372036854775808, 4620909390464243, 381128907510930132896351001478811136, and 125959558087561075514773639947125, due to space constraints. Although LMI solvers are deterministic, they can be sensitive to small changes in input values, and the numerical stability of the solution is contingent upon the conditioning of the problem.

Q2: Should I add normalization or optimization constraints (e.g., trace(P) minimization or I≤P≤ aI to obtain a Lyapunov matrix comparable to the published result?

Nope, I don't advise you to insist on replicating a matrix "P" similar to the published result for two reasons. First, the proposed LMI approach is intended to demonstrate the stability of your designed T–K fuzzy system, which I have verified in the results below. This should be sufficient to convince any experts familiar with fuzzy control systems and stability principles. There is little value in replicating the "P" from the published result unless you aim to provide a counterexample to the authors.

Secondly, I used the same approach to verify the authors' "P", and the results were inconsistent. One state matrix "A" led to a negative definite matrix (which implies stability), while the other did not.

% System matrix
syms Z1
A_general = [
            -2200/79, -2000/79;
             9223372036854775808/4620909390464243, (381128907510930132896351001478811136*Z1)/125959558087561075514773639947125
            ];

% Premise variable bounds
X_MG10(2) = 198.3325; 
x2min     = -89.4;   % Last point at which LMI is feasible
x2max     =  89.4;

f1max     = 1/(x2min + X_MG10(2));
f1min     = 1/(x2max + X_MG10(2));

Zmin      = f1min;
Zmax      = f1max;

% T–S fuzzy rule generation
premiseIdx = 2;
numPremise = length(premiseIdx);
numRules   = 2^numPremise;

A_rule = cell(numRules,1);
ruleIndices = dec2bin(0:numRules-1) - '0';

for iRule = 1:numRules
    A_local = A_general;
    if ruleIndices(iRule) == 0
        A_local = subs(A_local, Z1, Zmin);
    else
        A_local = subs(A_local, Z1, Zmax);
    end
    A_rule{iRule} = double(A_local);
end

% Common Lyapunov LMI formulation
N = numel(A_rule);
n = size(A_rule{1},1);

setlmis([])

p = lmivar(1,[n 1]);

for k = 1:N
    lmiterm([k 1 1 p], 1, A_rule{k}, 's');
end

S   = newlmi;
lmiterm([-S 1 1 p], 1, 1);

pp  = 1e-6;
lmiterm([ S 1 1 0], -pp*eye(n));

lmis = getlmis;
[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.056084
     2                        0.012870
     3                        0.012870
     4                        0.012870
     5                    2.404692e-03
     6                    3.445975e-05
     7                       -0.032602

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

P   = dec2mat(lmis, xfeas, p)

P = 2×2

1.0e+03 * 3.5077 0.0488 0.0488 0.0445

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Verify the LMI result:

A1  = A_rule{1}

A1 = 2×2

1.0e+03 * -0.0278 -0.0253 1.9960 0.0105

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

eig(A1)     % Manually check if A1 is stable

ans = 
   1.0e+02 *

  -0.0867 + 2.2397i
  -0.0867 - 2.2397i

Q1  = A1'*P + P*A1

Q1 = 2×2

1.0e+03 * -0.4798 -0.8460 -0.8460 -1.5361

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

eig(Q1)     % Check if Q1 is negative-definite

ans = 2×1

1.0e+03 * -2.0053 -0.0106

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

A2  = A_rule{2}

A2 = 2×2

1.0e+03 * -0.0278 -0.0253 1.9960 0.0278

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

eig(A2)     % Manually check if A2 is stable

ans = 
   1.0e+02 *

  -0.0004 + 2.2307i
  -0.0004 - 2.2307i

Q2  = A2'*P + P*A2

Q2 = 2×2

-479.7840 -3.3512 -3.3512 -0.2494

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

eig(Q2)     % Check if Q2 is negative-definite

ans = 2×1

-479.8074 -0.2260

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Verify the stability using the "P" from the paper:

P   = [103.1  1.43;
         1.43 1.31];

A1  = A_rule{1};
Q1  = A1'*P + P*A1;
eig(Q1)     % Check if Q1 is negative-definite

ans = 2×1

-60.1737 -18.3750

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

A2  = A_rule{2};
Q2  = A2'*P + P*A2;
eig(Q2)     % Check if Q2 is negative-definite

ans = 2×1

-34.2908 0.9655

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Manish Kumar]

Hi @Sam Chak,

Thank you for the clarification. I have two follow-up questions regarding numerical parameters in the LMI formulation.

1) Choice of lower bound in the positive-definiteness constraint

In my code, I enforce positive definiteness of the Lyapunov matrix as:

% lmiterm([-S 1 1 p],1,1e-3); % P >= 1e-3*I (works)

lmiterm([-S 1 1 p],1,1); % P >= I (also works)

Both choices lead to feasibility for my 2×2 system, but result in differently scaled Lyapunov matrices.

Is there a recommended or standard way to select this lower bound (e.g., problem-dependent scaling, numerical conditioning, or simply a small ε such as 1e-6)?

2) Choice of feasibility threshold on tmin

After solving with feasp, I accept feasibility when:

if tmin < -1e-6

Is 1e-6 a standard numerical threshold for declaring LMI feasibility, or should this be chosen relative to system dimension, matrix scaling, or solver accuracy?

Any practical guidelines on selecting these values would be very helpful.

[Sam Chak]

Hi @Manish Kumar

Q3. Choice of lower bound in the positive-definiteness constraint

This feature provides flexibility, allowing the user to "control" the design. I typically initiate the search with the following:

lmiterm([-S 1 1 p], 1, 1); % to find -P <= I, which literally implies P >= I

lmiterm([ S 1 1 0], 1); % I

This approach lets the LMI solver automatically search for feasible solutions, which are often satisfactory in my applications. However, if the positive-definite matrix P is excessively large or small compared to the state matrix A (to the extent that it may raise concerns among reviewers), I will adjust the weight.

Q4. Choice of feasibility threshold on tmin

Regarding the feasibility threshold tmin, as long as tmin < 0, the solutions are considered strictly feasible. It is important to note that having tmin < -1e-6 is not a standard threshold for determining feasibility. The LMI solver is programmed to find the optimal value of tmin​ within a number of iterations. However, if tmin​ is an extremely small negative value, reviewers may question the numerical stability of the solutions. In your case, since there are only two state matrices A, we can manually verify the feasibility of the positive-definite matrix P.