Estimate a region of attraction

[Marco Ghibaudi]

Dear professor,

I'm trying to estimate a region of attraction of the origin for a nonlinear discrete-time system with two equilibria: (0,0) and (2,4), in particular through some simulations I've seen that the (possible) region of attraction is quite small, on the other hand the result of my eta given by my code is a very large number (the kind of region of attraction that I'm trying to estimate is a level set of the founded Lyapunov function).

Code example:

close all

clear all

clc

% sdp variables

x1 = sdpvar(1,1);

x2 = sdpvar(1,1);

coeff_V = sdpvar(1,14);

eta = sdpvar(1,1);

% parameters

gamma = 1e-3;

% Lyapunov function structure (polynomial of degree 4)

V = 0;

k = 0;

for i = 0:4

   for j = 0:4

      if (i+j >= 1) & (i+j <= 4)

        k = k+1;

        V = V+coeff_V(k)*x1^j*x2^i;

      end

   end

end

% dynamics of the system

x1next = 1/4*x1*x2;

x2next = x1+1/2*x2;

% estimate of the Lyapunov function

F = [];

F = [F; sos(V-gamma*[x1 x2]*[x1 x2]')];

F = [F; sos(-(replace(V,[x1 x2],[x1next x2next])-V+gamma*[x1 x2]*[x1 x2]'))];

F = [F; sos(-(V-eta))];

solvesos(F,[],[],[coeff_V eta]);

[Johan Löfberg]

first, [V, V_coeff] = polynomial([x1;x2],4) simplifies your code...

However, the problem is trivially not feasible. You are asking for V to be a polynomial larger than a quadratic, and at the same time globally smaller than eta. Obviously not possible