We constrain the scaling relation between optical richness (λ) and halo mass (M) for a sample of Sloan Digital Sky Survey (SDSS) red-sequence Matched-filter Probabilistic Percolation (redMaPPer) galaxy clusters within the context of the Planck cosmological model. We use a forward modeling approach where we model the probability distribution of optical richness for a given mass, P(In λ/M). To model the abundance and the stacked lensing profiles, we use an emulator specifically built to interpolate the halo mass function and the stacked lensing profile for an arbitrary set of halo mass and redshift, which is calibrated based on a suite of high-resolution N-body simulations. We apply our method to 8312 SDSS redMaPPer clusters with 20 ≤ λ ≤ 100 and 0.10 ≤ z λ ≤ 0.33 and show that the lognormal distribution model for P(λ/M), with four free parameters, well reproduces the measured abundances and lensing profiles simultaneously. The constraints are characterized by the mean relation, , with and (68% CL), where the pivot mass scale M pivot = 3 ×1014 h -1 M o, and the scatter with and . We find that a large scatter in halo masses is required at the lowest-richness bins (20 ≤ λ ≲ 30) in order to reproduce the measurements. Without such a large scatter, the model prediction for the lensing profiles tends to overestimate the measured amplitudes. This might imply a possible contamination of intrinsically low-richness clusters due to the projection effects. Such a low-mass halo contribution is significantly reduced when applying our method to the sample of 30 ≤ λ ≤ 100.