using drake for region of attraction analysis

Question

I am having fun with pydrake for the region of attraction (RoA) analysis. There is one type of non-linear system involving trigonometry that got me stuck. A similar problem has been posted here.

As Russ points out in that post, there are 2 solutions: 1) shifting to another state-space; 2) using Taylor expansion. I am applying the 2nd approach to the pendulum case to get a feeling. The system dynamics are given in here. Here is what I have:

1. Linearize the system, and find an LQR, whose cost-to-go function will help to propose a Lyapunov candidate:

# Prepare
builder = DiagramBuilder()
pendulum = builder.AddSystem(PendulumPlant())
context = pendulum.CreateDefaultContext()
pendulum.get_input_port(0).FixValue(context, [0])
context.SetContinuousState([np.pi, 0])

# linearize the system at [pi, 0]
linearized_pendulum = Linearize(pendulum, context)

# LQR
Q = np.diag((10., 1.))
R = [1.]
(K, S) = LinearQuadraticRegulator(linearized_pendulum.A(), linearized_pendulum.B(), Q, R)

2. Set up pendulum's non-linear dynamics with the control input designed from LQR

# constants
pi = np.pi
g = 9.81
m = 1.
l = 0.5
d = 0.1

# non-linear dyamics
prog = MathematicalProgram()
x = prog.NewIndeterminates(2, "x")
x1 = x[0]
x2 = x[1]
Tsin = -(x1-pi) + (x1-pi)**3/6  # 3rd order taylor expansion of sin(x) around the fix point
r = [pi, 0]  # reference
u = -K.dot(x)[0] # control input
u2 = -K.dot(x-r)[0] # control input with reference
fn = [x2, (u-d*x2-Tsin*m*g*l)/(m*l*l)]

3. RoA analysis with the s-procedure, fix a Lyapunov function and find the Lagrange multiplier.

# cost-to-go of LQR as Lyapunov candidate
V = x.dot(S.dot(x)) 
Vdot = Jacobian([V], x).dot(fn)[0]

# Define the Lagrange multiplier.
lambda_ = prog.NewSosPolynomial(Variables(x), 2)[0].ToExpression()

# Optimization setup
rho = 2.
prog.AddSosConstraint(-Vdot + lambda_*(V-rho))
prog.AddSosConstraint(V)
prog.AddSosConstraint(lambda_)

# Print result
result = Solve(prog)
if (result.is_success()):
    print(f"Verified that {str(V)} < {rho} is in the region of attraction.")
    #print(Polynomial(result.GetSolution(lambda_)))
else:
    print("failed")

No matter how small I set up for $rho$, s-procedure always returns failed. However, from here, it is known that this example has at least an RoA of 2.

I don't know is this a numerical issue or my formulation is wrong. If you could shed some light on this, it would be much appreciated. Thank you!

Answer 1

Your Lyapunov function candidate isn't correct. Currently you use

V = x.dot(S.dot(x)) 

But I think you should use

V = (x - r).dot(S.dot(x - r))

BTW, as an anecdote (so take this with a grain of salt). I had success when I work on a more complicated problem with approach 1 (namely create a new variable xbar = x - r, and then rewrite the whole dynamics using xbar, and treat xbar as the indeterminates). In my problem when I tried approach 2 it gave me numerical error. So if approach 2 doesn't work on your problem, you can try approach 1.