Here is your simulation, it was really quick.

===

#special case when r = 1

r<- 1

#Number of players (or some large amount to calculate the proportion of players accurately)

b <- 10000

#Used to store the results

S <- double(b)

for(i in 1:b){

#a will track the number of successes within an attempt

a <- 0

#n is the number of trials within an attempt

n <-0

p <- 0.001

while(a < r){

#For every iteration, increase the number of trials by 1

n <- n + 1

p <- min(c(0.001+0.001*n,0.2))

#was a success or not

y <- rbinom(1,1,p)

if(y==1){a = a+1}

}

#record the number of attempts

S[i] <- n

#print every 100 iterations

if(i%%1000 == 0){print(i)}

}

#histogram

hist(S,freq=F,xlab="n",,main=paste("Proportion of people and their required number of trials for p = ",round(p,6), "and r = ", r))

#plot mean

abline(v = mean(S), col = "red",lwd = 2)

#plot 2.5% and 97.5% quantiles for the 95% middle

abline(v=quantile(S,c(0.025,0.975)),col="purple", lwd = 2)

===

The probability formula is 0.001 + 0.001*n (where n is the number of attempts, 0.001 is the decimal representation of 0.1%). I used min(0.001 + 0.001*n,0.2) because when n gets large enough and surpasses 0.2, min(formula,0.2) = 0.2 (which gives the cap you desired). The cap is hit when 0.001 + 0.001*n = 0.2 or when n = 199

The average number of trials is about 39, 95% of people will fall between 6 and 85. It's hard to work this out theoretically (unlike the previous cases) because the trials are no longer independent.