More of a clarification than an Issue...
Page 4 of the paper states that we initialize u to the LQR solution to help the NN get in the neighborhood of a reasonable linear control function. I am only familiar with LQR's in which we formulate a cost function in the form of the Riccati equation, and at each time step optimize for the lowest cost control(s). How does one generate the a linear control function from a traditional LQR optimization solution? Thank you

Hi, Chang, YaChien,
I'm recently learning your project. It's wonderful, but I'm so curious about the lie derivative that you implemented on the code, why you only do a diagonal operation when code are computing the lie derivative. It looks like different from the formula in the paper. I used the formal formula of lie derivative for learning, but there is no valid funciton proved by dreal4 even if the cost function is very small. Can you give me some suggestion?
Hope you can reply, thank you very much!

Could you provide some details/references about how to find the sum-of-squares polynomials for the inverted pendulum example in the plotting.ipynb file? Some existing code will be highly appreciated.

Thanks in advance.

@YaChienChang
In my case, there are six variables, which leads to the calculation time of CheckStatisfiability (condition, config) is too long to be applied. config has been shrunk to 0.1. Are there any strategies to shorten the computing time?

Hi! @YaChienChang
I was checking your amazing work and approach! However, when plotting the level sets of the Lyapunov functions, why do you substract some values from the NN and LQR approach?

ax.contour(X, Y, V-0.11, 0, linewidths=2, colors='k')
ax.contour(X, Y, V_lqr-2.6, 0, linewidths=2, colors='m',linestyles='--')
ax.contour(X, Y, V_SOS, 0, linewidths=2, colors='b')

Here, why do you substract V-0.11 or V_lqr-2.6?

I'd really appreciate if you can give me some light in here! Thank you!

Thank you for providing us a valuable coding.

I want to ask, how did you find out the specific Lyapunov risk term?

I am not sure the source of the term.

Thanks!

Hi Ya-Chien Chang,
I tried to verify another system with your code, but I couldn't find a ROA and Lyapunov function. I don't know if there is a problem in my choice of the parameters.
Take the inverted pendulum case for example, x is (-6,6), ball_lb is 0.5, ball_ub is 6, and the parameters of Lyapunov_risk.
Any suggestion on this? Thank you very much!

Vielen Dank für interessante Arbeit.
Ich habe versucht den Code für ein ganz simples System testen.
Es läuft aber findet kein Lyapunov funktion für das System und ich kann nicht vestehen woran das liegt?
Irgendwelche Idee?

Hi, I am trying to run your repository but I get the following error: "RuntimeError: expected scalar type Float but found Double" when trying to run the block of Learning and Falsification on both the inverted pendulum and the path following notebooks.
I am not sure how to proceed and I was wondering if you have faced this before.
More details about the error:

0 Lyapunov Risk= 73.07649993896484
===========Verifying==========
Not a Lyapunov function. Found counterexample:
x1 : [-5.607981400644867165, -5.607370072610038392]
x2 : [0.9742785792574937265, 0.974595727232512421]

RuntimeError Traceback (most recent call last)
in
13
14 while i < max_iters and not valid:
---> 15 V_candidate, u = model(x)
16 X0,u0 = model(x_0)
17 f = f_value(x,u)

~/anaconda3/envs/nnlyap/lib/python3.6/site-packages/torch/nn/modules/module.py in _call_impl(self, *input, **kwargs)
887 result = self._slow_forward(*input, **kwargs)
888 else:
--> 889 result = self.forward(*input, **kwargs)
890 for hook in itertools.chain(
891 _global_forward_hooks.values(),

in forward(self, x)
11 def forward(self,x):
12 sigmoid = torch.nn.Tanh()
---> 13 h_1 = sigmoid(self.layer1(x))
14 out = sigmoid(self.layer2(h_1))
15 u = self.control(x)

~/anaconda3/envs/nnlyap/lib/python3.6/site-packages/torch/nn/modules/module.py in _call_impl(self, *input, **kwargs)
887 result = self._slow_forward(*input, **kwargs)
888 else:
--> 889 result = self.forward(*input, **kwargs)
890 for hook in itertools.chain(
891 _global_forward_hooks.values(),

~/anaconda3/envs/nnlyap/lib/python3.6/site-packages/torch/nn/modules/linear.py in forward(self, input)
92
93 def forward(self, input: Tensor) -> Tensor:
---> 94 return F.linear(input, self.weight, self.bias)
95
96 def extra_repr(self) -> str:

~/anaconda3/envs/nnlyap/lib/python3.6/site-packages/torch/nn/functional.py in linear(input, weight, bias)
1751 if has_torch_function_variadic(input, weight):
1752 return handle_torch_function(linear, (input, weight), input, weight, bias=bias)
-> 1753 return torch._C._nn.linear(input, weight, bias)
1754
1755

RuntimeError: expected scalar type Float but found Double

Thank you for your interesting work. We have conducted a follow-up study on your research, where we solved for a controller and certified the Lyapunov conditions for arbitrary nonlinear dynamical systems in the discrete-time setting: https://github.com/jlwu002/nlc_discrete (NeurIPS 2023), paper.

We noticed that the controller used in the paper (and code) for Path Following dynamics does not have a bias term, meaning that the controller $u = 0$ at the equilibrium point $(d_e, \theta_e) = 0$. This does not seem feasible for the path-following setting: when the distance error and angle error are both zero (at equilibrium point), we still need the steering angle to ensure the object is following the unit circle target trajectory.

Hello,
I tried to calculate the LQR solution for the linearized model arround the origin, but I got different values from what you have in the published code, [-23.58639732, -5.31421063], and when I try to run your code with the LQR solution that I have found, the search for the Lyapunov function, which guarantees the stability of the system, does not converge to a valid function.
I also publish the matlab code that I used to find the LQR solutions, using MATLAB:

%%%%%%%%%%%
G = 9.81;
l = 0.5;
m = 0.15;
b = 0.1;

d = -b/(m*l*l);
c = G/l

A=[ 0, 1;
      c, d ]
B=[ 0;
  1/(m*l*l)
      ]
R=eye
Q=eye(2,2)
[K,S,P] = lqr(A,B,Q,R)

%%%%%%%%%%%%

the LQR gain is [1.9773 0.9762]

So can you please tell me how did you get the LQR values?

Thank you

The x0 in the paper should be the stable point, or target point, why in the code it is:
x_0 = torch.zeros([1, 2])?

Hi Ya-Chien,

I tried to run your inverted pendulum case on my own laptop and plot the Lyapunov function and its Lie derivative as shown below:

There are some regions that Lie derivative of V is positive. Also I found it is pretty weird but I am not sure which part I did it incorrectly, running the same code from Jupiter notebook on my laptop does not show Lyapunov conditions are satisfied from CheckLyapunov using dReal. Any suggestion on this? Thank you!