A method to decompose geophysical flows into a balanced flow (defined by its potential vorticity, PV) and an imbalanced component (inertia–gravity waves, IGWs) is introduced. The balanced flow, called the optimal potential vorticity (OPV) balance, is a solution of an IGW-permitting dynamics in which the amount of IGWs is minimal. The residual IGWs are those spontaneously generated by the vortical flow during the numerical integration in which the PV anomaly grows slowly over a sufficiently long ramp period toward a prescribed PV field. The OPV balanced flow is obtained, iteratively, in a cycle of backward and forward integrations where IGWs are removed and PV is restored in every loop. The method is applied to the flow of unsteady vortices in the three-dimensional baroclinic non-hydrostatic dynamics on the $f$-plane and to the single-layer shallow-water dynamics on the sphere. Both applications show that the iterative method converges strongly, after only a few iterations, to the balanced flow.