Issue in simple Swing Equation - need help and some suggestions

[juciwdjared]

First, really thanks for sharing this. However I find something that seems need improvements.
So I tried on Swing Equation, which should give me a similar result as using this method (Phase Portrait Plotter on 2D phase plane
[Version 1.0.3]
(https://www.mathworks.com/matlabcentral/fileexchange/110785-phase-portrait-plotter-on-2d-phase-plane#version_history_tab) (370 KB) by [Yu Zhang]

The codes I tried using your method are attached below (note some \times "*" operators are missing):

clf;
eq = [sin(pi/6);cos(pi/6);0];
dev = [1;1;4];
ang_ind = [1,2];
pow_sys = SOStab(eq,dev,ang_ind);
pow_sys.dynamics = [pow_sys.x(3)*pow_sys.x(2);-pow_sys.x(3)pow_sys.x(1);
0.5 - pow_sys.x(1) - 0.2pow_sys.x(3)
];

pow_sys.SoS_out(6,14,0.1);
vopow_sys = value(pow_sys.vcoef_outer);
wopow_sys = value(pow_sys.wcoef_outer);
pow_sys.plot_roa([1,2],3,'outer',0,"\theta","\omega");

The result looks like this, partially correct:

My comments are:

1. It seems your method is sensitive to the settings of D, T, and error (varepsilon). I tried difeerent settings but none of them give expected results and some of them seems not correct. Although maybe I'm wrong, could someone help on this?
2. You should remove the requirement on the initial acknowledgement of equilibrium point values. Refer to the method above I referred to.
3. The good thing of your method is it does NOT limited to the degree of 2 of ODE fcn.

[matteo-tacchi]

Thank you for your interest in the toolbox. Indeed, most solvers have difficulties in solving accurately the SDPs that come from the moment-SoS hierarchy, hence some numerical artifacts in the obtained results. That being said, I will answer your points one by one:

1. Indeed, the method is sensitive to all of its parameters, and identifying good settings is an open question so far; a quick comment on the "partially correct" plot that you shared: although it might sound unsatisfying, this often happens with the method (regardless of the implementation), and the quick fix just consists in ignoring all connected components that do not contain the original equilibrium point; indeed, the current method uses monomials as the basis of polynomials, and in practice (at reasonable degrees) this choice produces very rigid polynomials whose behavior is globally fixed by their local values; here, the solver puts all its power to produce the approximation containing the equilibrium point, at the price of numerical inaccuracy in other regions. Also, remember that the plots have no meaning outside of the admissible set X, as polynomials are only constrained within this set. One can compare it with Lyapunov functions: outside of their domain, these functions are totally free, and they only decrease along trajectory within a neighbourhood of the equilibrium point; the same occurs here, so that you can ignore the red lines at the bottom of the figure.
2. We wanted to remove this requirement, but it is needed both in the system's recasting (enforcing that the equilibrium point is 0) and in the plotting of the RoA approximation. Having 0 as the equilibrium is a must in polynomial optimization with monomials (because they vanish only at this point), and without this assumption numerical accuracy would be even worse than in the current setting.
3. Indeed, this method can handle a lot of dynamical systems, from polynomial systems with degree higher than 2 to trigonometric-polynomial systems; the next steps will be to include even more nonlinearities such as saturations.