We introduce the grand-reaction method for coarse-grained simulations of acid-base equilibria in a system coupled to a reservoir at a given pH and concentration of added salt. It can be viewed as an extension of the constant-pH method and the reaction ensemble, combining explicit simulations of reactions within the system, and grand-canonical exchange of particles with the reservoir. Unlike the previously introduced methods, the grand-reaction method is applicable to acid-base equilibria in the whole pH range because it avoids known artifacts. However, the method is more general, and can be used for simulations of any reactive system coupled to a reservoir of a known composition. To demonstrate the advantages of the grand-reaction method, we simulated a model system: A solution of weak polyelectrolytes in equilibrium with a buffer solution. By carefully accounting for the exchange of all constituents, the method ensures that all chemical potentials are equal in the system and in the multi-component reservoir. Thus, the grand-reaction method is able to predict non-monotonic swelling of weak polyelectrolytes as a function of pH, that has been known from mean-field predictions and from experiments but has never been observed in coarse-grained simulations. Finally, we outline possible extensions and further generalizations of the method, and provide a set of guidelines to enable safe usage of the method by a broad community of users.<br><br>