  A key characteristic of the anomalous sub-solution equation is that the
solution exhibits algebraic decay rate over long time intervals, which is often
refered to the Mittag-Leffler type stability. For a class of power nonlinear
sub-diffusion models with variable coefficients, we prove that their solutions
have Mittag-Leffler stability when the source functions satisfy natural decay
assumptions. That is the solutions have the decay rate
$\|u(t)\|_{L^{s}(\Omega)}=O\left( t^{-(\alpha+\beta)/\gamma} \right)$ as
$t\rightarrow\infty$, where $\alpha$, $\gamma$ are positive constants,
$\beta\in(-\alpha,\infty)$ and $s\in (1,\infty)$. Then we develop the structure
preserving algorithm for this type of model. For the complete
monotonicity-preserving ($\mathcal{CM}$-preserving) schemes developed by Li and
Wang (Commun. Math. Sci., 19(5):1301-1336, 2021), we prove that they satisfy
the discrete comparison principle for time fractional differential equations
with variable coefficients. Then, by carefully constructing the fine the
discrete supsolution and subsolution, we obtain the long time optimal decay
rate of the numerical solution $\|u_{n}\|_{L^{s}(\Omega)}=O\left(
t_n^{-(\alpha+\beta)/\gamma} \right)$ as $t_{n}\rightarrow\infty$, which is
fully agree with the theoretical solution. Finally, we validated the analysis
results through numerical experiments.
