  We develop a methodology for the computation of the $K\to \ell\nu_\ell
\ell'^+ \ell'^-$ decay width using lattice QCD and present an exploratory study
here. We use a scalar function method to account for the momentum dependence of
the decay amplitude and adopt the infinite volume reconstruction method to
reduce the systematic errors such as the temporal truncation effects and the
finite-volume effects. We then perform a four-body phase-space integral to
obtain the decay width. The only remaining technical problem is the possible
power-law finite-volume effects associated with the process of $K\to\pi\pi
\ell\nu_\ell\to \ell\nu_\ell \ell'^+ \ell'^-$, where the intermediate state
involves multiple hadrons. In this work, we use a gauge ensemble of twisted
mass fermion with a pion mass $m_\pi=352$ MeV and a nearly-physical kaon mass.
At this kinematics, the $\pi\pi$ in the intermediate state cannot be on shell
simultaneously as $2m_\pi>m_K$ and the finite-volume effects associated with
$\pi\pi$ state are exponentially suppressed. Using the developed methods
mentioned above, we calculate the branching ratios for four channels of $K\to
\ell\nu_\ell\ell'^+ \ell'^-$, and obtain the results comparable to the
experimental measurements and ChPT predictions. Our work demonstrates the
capability of lattice QCD to improve Standard Model prediction in $K\to
\ell\nu_\ell \ell'^+ \ell'^-$ decay width.
