There is a somewhat convoluted if direct resolution by accept-reject. First, a simple differentiation shows that the pdf of the distribution is
$$f(x)=(a+bx^p)\exp\left\{-ax-\frac{b}{p+1}x^{p+1}\right\}$$
Second, since
$$f(x)=ae^{-ax}\underbrace{e^{-bx^{p+1}/(p+1)}}_{\le 1}+bx^pe^{-bx^{p+1}/(p+1)}\underbrace{e^{-ax}}_{\le 1}$$
we have the upper bound
$$f(x)\le g(x)=ae^{-ax}+bx^pe^{-bx^{p+1}/(p+1)}$$
Third, considering the second term in $g$, take the change of variable $\xi=x^{p+1}$, i.e., $x=\xi^{1/(p+1)}$. Then$$\dfrac{\text{d}x}{\text{d}\xi}=\dfrac{1}{p+1}\xi^{\frac{1}{p+1}-1}=\dfrac{1}{p+1}\xi^{\frac{-p}{p+1}}$$
is the Jacobian of the change of variable. If $X$ has a density of the form $\kappa bx^pe^{-bx^{p+1}/(p+1)}$ where $\kappa$ is the normalising constant, then $\Xi=X^{1/(p+1)}$ has the density
$$\kappa b\xi^{\frac{p}{p+1}}e^{-b\xi/(p+1)}\,\dfrac{1}{p+1}\xi^{\frac{-p}{p+1}}=\kappa \dfrac{b}{p+1}e^{-b\xi/(p+1)}$$
which means that (i) $\Xi$ is distributed as an Exponential $\mathcal{E}(b/(p+1))$ variate and (ii) the constant $\kappa$ is equal to one. Therefore, $g(x)$ ends up being equal to the equally weighted mixture of an Exponential $\mathcal{E}(a)$ distribution and the $1/(p+1)$-th power of an Exponential $\mathcal{E}(b/(p+1))$ distribution, modulo a missing multiplicative constant of $2$ to account for the weights:
$$f(x)\le g(x)=2\left(\frac{1}{2} ae^{-ax}+\frac{1}{2} bx^pe^{-bx^{p+1}/(p+1)}\right)$$
And $g$ is straightforward to simulate as a mixture.
An R rendering of the accept-reject algorithm is thus
simuF <- function(a,b,p){
reepeat=TRUE
while (reepeat){
if (runif(1)<.5) x=rexp(1,a) else
x=rexp(1,b/(p+1))^(1/(p+1))
reepeat=(runif(1)>(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/
(a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1))))}
return(x)}
and for an n-sample:
simuF <- function(n,a,b,p){
sampl=NULL
while (length(sampl)<n){
x=u=sample(0:1,n,rep=TRUE)
x[u==0]=rexp(sum(u==0),b/(p+1))^(1/(p+1))
x[u==1]=rexp(sum(u==1),a)
sampl=c(sampl,x[runif(n)<(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/
(a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1)))])
}
return(sampl[1:n])}
Here is an illustration for a=1, b=2, p=3:
