  To interpret adaptive-optics observations of (216) Kleopatra, we need to
describe an evolution of multiple moons, orbiting an extremely irregular body
and including their mutual interactions. Such orbits are generally
non-Keplerian and orbital elements are not constants. Consequently, we use a
modified $N$-body integrator, which was significantly extended to include the
multipole expansion of the gravitational field up to the order $\ell = 10$. Its
convergence was verified against the `brute-force' algorithm. We computed the
coefficients $C_{\ell m},S_{\!\ell 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 $\chi^2$
metric then accounts for the absolute astrometry, the relative astrometry (2nd
moon with respect to 1st), angular velocities, and also silhouettes,
constraining the pole orientation. This allowed us to derive the orbital
elements of Kleopatra's two moons. 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, $P_1 = (1.822359\pm0.004156)\,{\rm d}$, $P_2 =
(2.745820\pm0.004820)\,{\rm 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), so that regular librations of the critical
argument are not present. The resulting mass of Kleopatra, $m_1 =
(1.49\pm0.16)\cdot10^{-12}\,M_\odot$ or $2.97\cdot10^{18}\,{\rm kg}$, is
significantly lower than previously thought. An implication explained in the
accompanying paper (Marchis et al.) is that (216) Kleopatra is a critically
rotating body.
