I’m going to make some assumptions in this answer:
Two moons have masses $m_1$ and $m_2$ and orbit a planet of mass $m_p$, where $m_1, m_2\ll m_p$.
The planet is far enough away from the star that any gravitational/tidal effects from that star are negligible.
There are no other planets capable of destabilizing any moons in orbits reasonably close to the planet (i.e. well within its Hill sphere). This, along with the second assumption, means that we can effectively treat the moon system as a miniature planetary system.
There are two cases to look at: where $m_1 \gg m_2$, and where $m_1\sim m_2$.
1. $m_1 \gg m_2$
In this scenario, we see the possibility of $m_1$ capturing $m_2$ just as Neptune is thought to have captured Triton, involving a three-body collisionless encounter. $m_2$ would originally have been part of a binary moon system of some sort (see my second section) which then interacted with $m_1$; the other binary partner was ejected and $m_2$ became a satellite of $m_1$ (see Agnor & Hamilton (2006)). Assuming that the ejected binary partner had a mass $m_3$, the three bodies would have had to have interacted at a distance
$$r=a\left(\frac{3m_1}{m_2+m_3}\right)^{1/3}\tag{1}$$
where $a$ was the semi-major axis of $m_2$ and $m_3$. There shouldn’t be any issues with applying this model to the moon system.
2. $m_1\sim m_2$
This is a standalone scenario, but I realized that it is also needed to explain the formation of the original binary system in the first setup. It has the advantage that no third body is needed (and thus no other binary system is needed), but it has the disadvantage that a relatively narrow class of initial orbits will permit a successful finish.
Ochai et al. (2014) apply the phenomenon of tidal energy dissipation to the formation of binary planets (here, we apply it to binary moons, because if $m_1~m_2$, it may be more accurate to refer to the system as a binary system, rather than as a satellite and sub-satellite). Given radii of $R_1$ and $R_2$ and a pericenter distance of $q_{12}$, the energy dissipated after each encounter is
$$E=\frac{Gm_1^2}{R_2}\left[\left(\frac{R_2}{q_{12}}\right)^6T_2(\eta_2)+\left(\frac{R_2}{q_{12}}\right)^8T_3(\eta_2)\right]+\frac{Gm_2^2}{R_1}\left[\left(\frac{R_1}{q_{12}}\right)^6T_2(\eta_1)+\left(\frac{R_1}{q_{12}}\right)^8T_3(\eta_1)\right]\tag{2}$$
where $\eta_i\equiv[m_i/(m_1+m_2)]^{½}$ and $T_2$ and $T_3$ are fifth degree polynomial functions with coefficients given in Portegies Zwart & Meinen (1993).
The authors carried out simulations of three gas giant planets orbiting a star (some simulations involved two, while others had the third play a minor role), and found different results for different values of the inner planet’s semi-major axis:
Here, “HJs” stands for “Hot Jupiters”.
It is certainly true that long-term (or even short-term!) stability might be problematic, especially because of tidal interactions with the parent planet. However, many setups will lead to the successful formation of a binary planet.
All of this, of course, assumes that the less massive moon lies within the Hill sphere of the more massive one.