of Ising-type models, which is only feasible when long-range interactions are neglected. Here, this formalism is extended to include long-range interactions by introducing a modified free energy involving only effective short-range interaction parameters. These parameters can be systematically calculated by using the Gibbs-Bogoliubov variational principle. The technique is illustrated with the calculation of titration curves of homogeneous and heterogeneous polyelectrolytes in a wide range of ionic strengths. The correction of the site protonation free energy (first order correction) is enough to obtain an excellent agreement between theory and Monte Carlo simulations. Corrections to other cluster parameters (higher-order corrections) are also implemented. In general, the correction to a particular parameter represents the average change in the long-range energy when a new interaction is created in the polyelectrolyte. The method presented here represents an improvement in the description of the ionization state of polyelectrolytes that can be relevant in a wide range of areas