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 Kleopatras 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 Kleopatras 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.822359pm0.004156),{rm d}$, $P_2 = (2.745820pm0.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.49pm0.16)cdot10^{-12},M_odot$ or $2.97cdot10^{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.