Adjust graph done in R

I would like help with my graph. I generated a graph, but I would like to insert IdealPoint and SolutionFinal together with the graph. I didn't insert all the code, because it's really big. Unfortunately, I was unable to make executable code, but I believe I don't have to for this case. I would just like someone to help me insert these two (IdealPoint and SolutionFinal) in the graph.

library(rdist)
library(geosphere)
library(tidyverse)

df<-structure(list(Propertie = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19), Latitude = c(-23.8, -23.9, -23.5, -23.4, -23.6,-23.9, -23.2, -23.5, -23.8, -23.7, -23.8, -23.9, -23.4, -23.9, 
                                                                                                -23.9, -23.2, -23.3, -23.7, -23.8), 
                   Longitude = c(-49.1, -49.3,-49.4, -49.7, -49.6, -49.6, -49.6, -49.6, -49.6, -49.6, -49.7,-49.2, -49.5, -49.8, -49.5, -49.3, -49.3, -49.2, -49.5), 
                   Waste = c(526,350, 526, 469, 285, 175, 175, 350, 350, 175, 350, 175, 175, 364,175, 175, 350, 45.5, 54.6)), 
              class = "data.frame", row.names = c(NA, -19L)) 

K<-1
Filter1<-1
Filter2<-1

coordinates<-subset(df,select=c("Latitude","Longitude")) 
d<-distm(coordinates[,2:1]) 
diag(d)<-1000000 
min_distance<-as.matrix(apply(d,MARGIN=2,FUN=min))
limit<-mean(min_distance)+3*sd(min_distance) 

search_vec <- function(mat, vec, dim = 1, tol = 1e-7, fun = all)
  which(apply(mat, dim, function(x) fun((x - vec) > tol)))
indice_exclude<-search_vec(min_distance,limit,fun=any)
if(is_empty(indice_exclude)==FALSE){
  coordinates<-coordinates[-c(indice_exclude),]
  df<-df[-c(indice_exclude),]
  d<-d[-c(indice_exclude),-c(indice_exclude)]}

#Filter 1
Q1<-matrix(quantile(df$Waste, probs = 0.25)) 
Q3<-matrix(quantile(df$Waste, probs = 0.75))
L<-Q1-1.5*(Q3-Q1)
S<-Q3+1.5*(Q3-Q1)
if (Filter1==2){
  d<-subset(d,df$Waste>L[1],df$Waste>L[1])
  df<-subset(df,Waste>L[1])
}
if (Filter1==3){
  d<-subset(d,df$Waste<S[1],df$Waste<S[1])
  df<-subset(df,Waste<S[1])
}
coordinates<-subset(df,select=c("Latitude","Longitude"))
d<-as.dist(d)
fit.average<-hclust(d,method="average") 
npropriedades<-dim(df)[1]

p<-1
clusters<-cutree(fit.average, p) 
nclusters<-matrix(table(clusters))
while (min(nclusters)>1) {
  p<-p+1
  clusters<-cutree(fit.average, p) 
  nclusters<-matrix(table(clusters))}
p<-p-1

mean<-matrix(nrow=p-1,ncol=3)
colnames(mean)<-c("coverage","Range of Coverage","Waste")
l<-1
for(k in 2:p){
  
  clusters<-cutree(fit.average, k) 
  nclusters<-matrix(table(clusters))  
  df$cluster <- clusters 
  
  centro_massa<-matrix(nrow=k,ncol=2)
   for(i in 1:k){
    centro_massa[i,]<-c(mean(subset(df,cluster==i)$Latitude),
                        mean(subset(df,cluster==i)$Longitude))}
  coordinates$cluster<-clusters 
  centro_massa<-cbind(centro_massa,matrix(c(1:k),ncol=1)) 
  
  coverage<-matrix(nrow=k,ncol=1)
  for(i in 1:k){
    coverage[i,]<-max(data.matrix(dist(rbind(subset(coordinates,cluster==i)[,2:1],centro_massa[i,2:1])))[nclusters[i]+1,])}
  
 soma_Waste<-matrix(nrow=k,ncol=1)
  for(i in 1:k){
    soma_Waste[i,]<-sum(subset(df,cluster==i)["Waste"])
  }
  
  mean[l,]<-c(max(coverage),max(coverage)-min(coverage),min(soma_Waste)) 
  l<-l+1
}

q1<-matrix(quantile(mean[,1], probs = 0.25))
q3<-matrix(quantile(mean[,1], probs = 0.75))
s<-q3+1.5*(q3-q1)
if(Filter2==1 & s[1]>0){
  Mean<-subset(mean, mean[,1] <=s[1])
} else{Mean<-mean}

Mean<-Mean[,2:3]

IdealPoint<-as.matrix(t(c(min(Mean[,1]),max(Mean[,2]))))
aux_solution<-as.matrix(dist(rbind(Mean,IdealPoint)))
IdealPoint

distancia_solution<-min(aux_solution[as.matrix(dim(Mean))[1,1]+1,1:as.matrix(dim(Mean))[1,1]])
a<-which(aux_solution[dim(Mean)[1]+1,]==distancia_solution)
SolutionFinal<-Mean[a[1],]
SolutionFinal

plot(Mean[,1],Mean[,2],xlab="Range of Coverage", ylab="Minimal waste")

If you change the last line to

points(SolutionFinal[1],SolutionFinal[2])

Then it does work. But the problem is that it is stacking all the points. So SolutionFinal[1],SolutionFinal[2] is the same as IdealPoint[,1],IdealPoint[,2] and the top-left point in your plot.

You can see them if you jitter the points a bit with jitter

plot(jitter(Mean[,1],1.5),Mean[,2],xlab="Range of Coverage", ylab="Minimal waste")
points(jitter(IdealPoint[,1],1.5),IdealPoint[,2],col="red")
points(jitter(SolutionFinal[1],1.5),jitter(SolutionFinal[2],1.5),col="purple")