This paper reports a search for Higgs boson pair (hh) production in association with a vector boson ($$W\; {\text {o}r}\; Z$$$W\text{o}rZ$) using 139 fb$$^{-1}$$${}^{-1}$of proton–proton collision data at$$\sqrt{s}=13\,\text {TeV}$$$\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 $$$W\to \ell \nu ,Z\to \ell \ell ,\nu \nu$with$$\ell =e, \mu $$$\ell =e,\mu$) and the Higgs bosons each decay into a pair ofb-quarks. It targetsVhhsignals from both non-resonanthhproduction, present in the Standard Model (SM), and resonanthhproduction, 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-resonantVhhproduction 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 resonanceH, in the mass range 260–1000 GeV, that decays intohh, and the other is the production of a heavier neutral pseudoscalar resonanceAthat decays into aZboson andHboson, where theAboson mass is 360–800 GeV and theHboson mass is 260–400 GeV. Constraints are also derived in the parameter space of two-Higgs-doublet models.