Aims. To interpret adaptive-optics observations of (216) Kleopatra, we need to describe an evolution of multiple moons orbiting an extremely irregular body and include their mutual interactions. Such orbits are generally non-Keplerian and orbital elements are not constants. Methods. Consequently, we used a modified N-body integrator, which was significantly extended to include the multipole expansion of the gravitational field up to the order ℓ = 10. Its convergence was verified against the 'brute-force' algorithm. We computed the coefficients Cℓ m, Sℓ m for Kleopatra's shape, assuming a constant bulk density. For Solar System applications, it was also necessary to implement a variable distance and geometry of observations. Our χ2 metric then accounts for the absolute astrometry, the relative astrometry (second moon with respect to the first), angular velocities, and silhouettes, constraining the pole orientation. This allowed us to derive the orbital elements of Kleopatra's two moons. Results. Using both archival astrometric data and new VLT/SPHERE observations (ESO LP 199.C-0074), we were able to identify the true periods of the moons, P1 = (1.822359 ± 0.004156) d, P2 = (2.745820 ± 0.004820) d. They orbit very close to the 3:2 mean-motion resonance, but their osculating eccentricities are too small compared to other perturbations (multipole, mutual), meaning that regular librations of the critical argument are not present. The resulting mass of Kleopatra, m1 = (1.49 ± 0.16) × 10-12 M· or 2.97 × 1018 kg, is significantly lower than previously thought. An implication explained in the accompanying paper is that (216) Kleopatra is a critically rotating body.