  We constrain the scaling relation between optical richness ($\lambda$) and
halo mass ($M$) for a sample of SDSS redMaPPer galaxy clusters within the
context of the {\it Planck} cosmological model. We use a forward modeling
approach where we model the probability distribution of optical richness for a
given mass, $P(\ln \lambda| 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 8,312 SDSS redMaPPer clusters with $20\le
\lambda \le 100$ and $0.10\le z_{\lambda}\le0.33$, and show that the log-normal
distribution model for $P(\lambda|M)$, with four free parameters, well
reproduces the measured abundances and lensing profiles simultaneously. The
constraints are characterized by the mean relation, $\left\langle
\ln{\lambda}\right\rangle(M)=A+B\ln(M/M_{\rm pivot})$, with
$A=3.207^{+0.044}_{-0.046}$ and $B=0.993^{+0.041}_{-0.055}$ (68\%~CL), where
the pivot mass scale $M_{\rm pivot}=3\times 10^{14} h^{-1}M_\odot$, and the
scatter $\sigma_{\mathrm{\ln\lambda}|M}=\sigma_0+q\ln(M/M_{\rm pivot})$ with
$\sigma_0=0.456^{+0.047}_{-0.039}$ and $q=-0.169^{+0.035}_{-0.026}$. We find
that a large scatter in halo masses is required at the lowest richness bins
($20\le \lambda \lesssim 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\le \lambda \le 100$.
