Adding new polygon attribute column based on location withing another polygon in R

I have here an example with two layers – one containing two large polygons Pol_large, and second containing 4 small polygons Pol_small. For Pol_small
I want now to find the location of each of the small polygons within the large polygons creating a new attribute column location and transporting the attributes from the Pol_large$id column. It would look following:

Pol_small$location
[1] "A1" "A1" "A2" "A2"

Not sure which funciton in sf or sp would be best to do it.

library(sp)

create_poly <- function(coordin,id_number){
  coordinates = matrix(c(coordin[1], coordin[3],
                         coordin[1], coordin[4],
                         coordin[2], coordin[4],
                         coordin[2], coordin[3],
                         coordin[1], coordin[3]), 
                       ncol = 2, byrow = TRUE)
  P_1 = Polygon(coordinates)
  Ps_1 = SpatialPolygons(list(Polygons(list(P_1), ID = id_number)), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  ID_1 <- sapply(slot(Ps_1, "polygons"), function(x) slot(x, "ID"))
  df <- data.frame(rep(0, length(ID_1)), row.names=ID_1)
  Ps1_df <- SpatialPolygonsDataFrame(Ps_1, df)
  Ps1_df
}

# Create shapefile with polygon boundaries of interest:
Pol_1 <- create_poly(c(50,100,10,50),1)
Pol_2 <- create_poly(c(110,170,0,50),1)
Pol_3 <- create_poly(c(55,65,40,45),1)
Pol_4 <- create_poly(c(70,90,30,50),1)
Pol_5 <- create_poly(c(145,160,20,40),1)
Pol_6 <- create_poly(c(120,135,10,30),1)

# Join the polygons into one layer and assign attributes:
Pol_large <- do.call(rbind, list(Pol_1, Pol_2))
Pol_large$id <- c("A1","A2")
Pol_small <- do.call(rbind, list(Pol_3, Pol_4, Pol_5, Pol_6))
Pol_small$id <- c("b1","b2","b3","b4")

# Plot
plot(Pol_large, axes = TRUE)
plot(Pol_small, add=T, col = "blue")

Is this solution optimal? For some reason, in my real large example, this does not seem to work exactly as it should – specifically some smaller polygons with boundaries touching the edges of the large polygons do not get any assignment, they return integer(0)

library(sf)    
a <- unlist(st_within(st_as_sf(Pol_small),st_as_sf(Pol_large)))

for (i in 1:length(a)){
  g <- a[i]
  Pol_small$location[i] <- Pol_large$id[g]
}