Variational method for integrability-breaking Richardson-Gaudin models