We present a variational method for approximating the ground state of spin
models close to (Richardson-Gaudin) integrability. This is done by
variationally optimizing eigenstates of integrable Richardson-Gaudin models,
where the toolbox of integrability allows for an efficient evaluation and
minimization of the energy functional. The method is shown to return exact
results for integrable models and improve substantially on perturbation theory
for models close to integrability. For large integrability-breaking
interactions, it is shown how (avoided) level crossings necessitate the use of
excited states of integrable Hamiltonians in order to accurately describe the
ground states of general non-integrable models