This paper reports a search for Higgs boson pair (hh) production in association with a vector boson ($W\; {\text {o}r}\; Z$) using 139 fb$^{-1}$ of proton–proton collision data at $\sqrt{s}=13\,\text {TeV}$ recorded with the ATLAS detector at the Large Hadron Collider. The search is performed in final states in which the vector boson decays leptonically ($W\rightarrow \ell \nu ,\, Z\rightarrow \ell \ell ,\nu \nu $ with $\ell =e, \mu $) and the Higgs bosons each decay into a pair of b-quarks. It targets Vhh signals from both non-resonant hh production, present in the Standard Model (SM), and resonant hh production, as predicted in some SM extensions. A 95% confidence-level upper limit of 183 (87) times the SM cross-section is observed (expected) for non-resonant Vhh production when assuming the kinematics are as expected in the SM. Constraints are also placed on Higgs boson coupling modifiers. For the resonant search, upper limits on the production cross-sections are derived for two specific models: one is the production of a vector boson along with a neutral heavy scalar resonance H, in the mass range 260–1000 GeV, that decays into hh, and the other is the production of a heavier neutral pseudoscalar resonance A that decays into a Z boson and H boson, where the A boson mass is 360–800 GeV and the H boson mass is 260–400 GeV. Constraints are also derived in the parameter space of two-Higgs-doublet models.