The ionization state of charged macromolecules in solution is usually determined by the extent of the binding processes. These processes are very sensitive to the ionic strength of the medium, which are of long-range nature. The ionization properties of weak polyelectrolytes can be described by means 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.represents the average change in the long range energy when a new interaction (of the type described by the cluster parameter), is created in the polyelectrolyte. The method presented here represents an improvement in the description of the ionization state of charged macromolecules that can be relevant in a wide range of areas such as biochemistry, environmental chemistry, materials science, etc.