Let there be given a $2n$-dimenional real symplectic manifold $(M,\omega)$ with a globally defined real function $H:M\times[t_i,t_f] \to \mathbb{R}$, which we will call the Hamiltonian. The time evolution is governed by Hamilton's (or equivalently Liouville's) equations of motion. Here $t\in[t_i,t_f]$ is time.
On one hand, there is the notion of complete integrability, aka. Liouville integrability, or sometimes just called integrability. This means that there exist $n$ independent globally defined real functions
$$I_i, \qquad i\in\{1, \ldots, n\},$$
(which we will call action variables), that pairwise Poisson commute,
$$ \{I_i,I_j\}_{PB}~=~0, \qquad i,j\in\{1, \ldots, n\}.$$
On the other hand, given a fixed point $x_{(0)}\in M$, under mild regularity assumptions, there always exists locally (in a sufficiently small open Darboux$^1$ neighborhood of $x_{(0)}$) an $n$-parameter complete solution for Hamilton's principal function
$$S(q^1, \ldots, q^n; I_1, \ldots,I_n; t)$$
to the Hamilton-Jacobi equation, where
$$I_i, \qquad i\in\{1, \ldots, n\},$$
are integration constants. This leads to a local version of property 1.
The main point is that the global property 1 is rare, while the local property 2 is generic.
--
$^1$ A Darboux neighborhood here means a neighborhood where there exists a set of canonical coordinates aka. Darboux coordinates $(q^1, \ldots, q^n;p_1, \ldots, p_n)$, cf. Darboux' Theorem.
1) A constant of motion
$f(z,t)$ is a (globally defined, smooth) function $f:M\times [t_i,t_f] \to
\mathbb{R}$ of the dynamical variables $z\in M$ and time $t\in[t_i,t_f]$,
such that the map $$[t_i,t_f]~\ni ~t~~\mapsto~~f(\gamma(t),t)~\in~ \mathbb{R}$$
doesn't depend on time for every solution curve $z=\gamma(t)$ to the equations of motion of the system.
An integral of motion/first integral
is a constant of motion $f(z)$ that doesn't depend explicitly on time.
2) In the following let us for simplicity restrict to the case where the system is a finite-dimensional autonomous$^1$ Hamiltonian system with Hamiltonian $H:M \to \mathbb{R}$ on a $2N$-dimensional symplectic manifold $(M,\omega)$.
Such system is called (Liouville/completely) integrable if
there exist $N$ functionally independent$^2$, Poisson-commuting, globally defined functions $I_1, \ldots, I_N: M\to \mathbb{R}$, so that the Hamiltonian $H$ is a function of $I_1, \ldots, I_N$, only.
Such integrable system is called maximally superintegrable if
there additionally exist $N-1$ globally defined integrals of motion $I_{N+1}, \ldots, I_{2N-1}: M\to \mathbb{R}$, so that the combined set $(I_{1}, \ldots, I_{2N-1})$ is functionally independent.
It follows from Caratheodory-Jacobi-Lie theorem that every finite-dimensional autonomous Hamiltonian system on a symplectic manifold $(M,\omega)$ is locally maximally superintegrable in sufficiently small local neighborhoods around any point of $M$ (apart from critical points of the Hamiltonian).
The main point is that (global) integrability is rare, while local integrability is generic.
--
$^1$ An autonomous Hamiltonian system means that neither the Hamiltonian $H$ nor the symplectic two-form $\omega$ depend explicitly on time $t$.
$^2$ Outside differential geometry $N$ functions $I_1, \ldots, I_N$ are called functionally independent if $$\forall F:~~ \left[z\mapsto F(I_1(z), \ldots, I_N(z)) \text{ is the zero-function} \right]~~\Rightarrow~~ F \text{ is the zero-function}.$$
However within differential geometry, which is the conventional framework for dynamical systems, $N$ functions $I_1, \ldots, I_N$ are called functionally independent if $\mathrm{d}I_1\wedge \ldots\wedge \mathrm{d}I_N\neq 0$ is nowhere vanishing. Equivalently, the rectangular matrix $$\left(\frac{\partial I_k}{\partial z^K}\right)_{1\leq k\leq N, 1\leq K\leq 2N}$$ has maximal rank in all points $z$. If only $\mathrm{d}I_1\wedge \ldots\wedge \mathrm{d}I_N\neq 0$ holds a.e., then one should strictly speaking strip the symplectic manifold $M$ of these singular orbits.
Best Answer
I think one must distinguish between chaotic "Hamiltonian" systems and chaotic dissipative systems. In the latter, the phase space volume is not conserved, so it is much more difficult to find "integrals of motion" because Liouville's theorem is broken. Remember, a quantity "A" is an integral of motion if
$\frac{dA}{dt} = \frac{\partial A}{\partial t} + \{A,H\} = 0$, where $H$ is the Hamiltonian. For dissipative chaotic systems, you can't even write down $H$, so it is difficult to see how one could generally find integrals/constants of motion of the system.
However, there is an important class of systems that show up in cosmology for example where you have Hamiltonian "chaos", where essentially the trajectories of the system exhibit all of the properties of chaos: sensitive dependence on initial conditions, diverging trajectories over time, but, the system still has attractors: a famous example is the dynamics of a closed anisotropic universe / Bianchi IX, in shameless self-promotion here: https://arxiv.org/pdf/1311.0389.pdf (in particular, see Page 27) This has of course led to wide debates for years in the cosmology community of whether this is "Really" chaos, since, in principle, the trajectories are predictable, but, I hope this answers your question.
Further, with respect to your Billiards problem / the famous Hadamard billiards, as you can see it is the same as the diagram on Page 27. Therefore, the billiard problem is also an example of Hamiltonian / deterministic / non-dissipative chaos. The phase space has an asymptotic attractor. This hopefully demonstrates that integrals of motion such as the one you found above ($E$ is the total energy of the system, and in this case, is the Hamiltonian, $H$) are only really possible if one can write down a Hamiltonian by Liouville's theorem.