Accurate and efficient modeling of binary black holes (BBHs) is crucial for the detection of gravitational waves (GWs) emitted by them. Closed-form solutions to these systems when they are in the initial inspiral phase are highly sought after and have been worked out by many groups in the post-Newtonian (PN) approximation. Most of these solutions are valid only in certain limits (small eccentricity, no spins, equal mass, etc.). Establishing the integrable nature of PN BBHs opens up the possibility of constructing closed-form solutions since integrability precludes chaos and guarantees the existence of action-angle variables. In this talk, I will discuss our series of efforts in establishing the integrable nature of the most general BBH system (arbitrary masses, spins, and eccentricity) as per the Liouville-Arnold theorem at 2PN order. I will also discuss our recently derived action-angle-based solution for these systems at 1.5PN order.