A search for the rare decays $B_s^0 \to \pi^+\pi^-\mu^+\mu^-$ and $B^0 \to \pi^+\pi^-\mu^+\mu^-$ is performed in a data set corresponding to an integrated luminosity of 3.0 fb$^{-1}$ collected by the LHCb detector in proton-proton collisions at centre-of-mass energies of 7 and 8 TeV. Decay candidates with pion pairs that have invariant mass in the range 0.5-1.3 GeV/$c^2$ and with muon pairs that do not originate from a resonance are considered. The first observation of the decay $B_s^0 \to \pi^+\pi^-\mu^+\mu^-$ and the first evidence of the decay $B^0 \to \pi^+\pi^-\mu^+\mu^-$ are obtained and the branching fractions, restricted to the dipion-mass range considered, are measured to be $\mathcal{B}(B_s^0 \to \pi^+\pi^-\mu^+\mu^-)=(8.6\pm 1.5\,({\rm stat}) \pm 0.7\,({\rm syst})\pm 0.7\,({\rm norm}))\times 10^{-8}$ and $\mathcal{B}(B^0 \to \pi^+\pi^-\mu^+\mu^-)=(2.11\pm 0.51\,({\rm stat}) \pm 0.15\,({\rm syst})\pm 0.16\,({\rm norm}) )\times 10^{-8}$, where the third uncertainty is due to the branching fraction of the decay $B^0\to J/\psi(\to \mu^+\mu^-)K^*(890)^0(\to K^+\pi^-)$, used as a normalisation.