  This paper presents a study of the use of numerical simulation and Bayesian
optimisation techniques to investigate the dynamics of celestial systems.
Initially, the study focuses on Lagrange points in restricted three-body
systems where a 2D three-body system simulator is employed to locate the five
Lagrange points. An appropriate loss function is developed to capture the
gravitational stability of the system, and the stability properties of the
different Lagrange points are explored. Additionally, the study investigates
how varying the number of variables for the satellite impacts the search for
the Lagrange points. Finally, the scope of the study is expanded to explore
stable configurations in multi-star systems represented by regular convex
n-gons. In this case, Bayesian optimisation is used to find suitable settings
for the n-gon's radius and the stars' velocity vectors, such that the overall
system is stable.
