Results are presented from a search for physics beyond the standard model in proton-proton collisions at s = 13 TeV in channels with two Higgs bosons, each decaying via the process H → bb¯, and large missing transverse momentum. The search uses a data sample corresponding to an integrated luminosity of 137 fb⁻¹ collected by the CMS experiment at the CERN LHC. The search is motivated by models of supersymmetry that predict the production of neutralinos, the neutral partners of the electroweak gauge and Higgs bosons. The observed event yields in the signal regions are found to be consistent with the standard model background expectations. The results are interpreted using simplified models of supersymmetry. For the electroweak production of nearly mass-degenerate higgsinos, each of whose decay chains yields a neutralino (χ~⁰₁) that in turn decays to a massless goldstino and a Higgs boson, (χ~⁰₁) masses in the range 175 to 1025 GeV are excluded at 95% confidence level. For the strong production of gluino pairs decaying via a slightly lighter (χ~⁰₂) to H and a light (χ~⁰₁), gluino masses below 2330 GeV are excluded.