## Artificial intelligence for automated and robotized Welding
## Import libraries
library(knitr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
### Number of vertices or Weld spots that need to be welded
n <- 10
## Boundary of our Euclidean Space
max_x <- 500
max_y <- 500
## Generate welding spots using vision sensors
set.seed(6000)
Weld_Spot <- data.frame(id = 1:n, x = runif(n, max = max_x), y = runif(n, max = max_y))
ggplot(Weld_Spot, aes(x, y)) +
geom_point()
## Generate Weld spot dimensions from CAD for precision and accuracy
CAD_File <- as.matrix(stats::dist(select(Weld_Spot, x, y), diag = TRUE, upper = TRUE))
dist_fun <- function(i, j) {
vapply(seq_along(i), function(k) CAD_File[i[k], j[k]], numeric(1L))
}
## Model formulation
library(ompr)
model <- MIPModel() %>%
# we create a variable that is 1 iff we travel from Vertex spot i to j
add_variable(x[i, j], i = 1:n, j = 1:n,
type = "integer", lb = 0, ub = 1) %>%
# a helper variable for the MTZ formulation of the TSP
add_variable(u[i], i = 1:n, lb = 1, ub = n) %>%
# minimize travel distance
set_objective(sum_expr(dist_fun(i, j) * x[i, j], i = 1:n, j = 1:n), "min") %>%
# you cannot go to the same vertex
set_bounds(x[i, i], ub = 0, i = 1:n) %>%
# leave each vertex. in this case each vertex is a transient state
add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n) %>%
#
# visit each vertex
add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 1:n) %>%
# ensure no subtours (arc constraints)
add_constraint(u[i] >= 2, i = 2:n) %>%
add_constraint(u[i] - u[j] + 1 <= (n - 1) * (1 - x[i, j]), i = 2:n, j = 2:n)
model
## Mixed integer linear optimization problem
## Variables:
## Continuous: 10
## Integer: 100
## Binary: 0
## Model sense: minimize
## Constraints: 110
## Mixed integer linear optimization Algorithm
## Variables:
## Continuous: 10
## Integer: 100
## Binary: 0
## Model sense: minimize
## Constraints: 110
## Model is solved using solver library called GLPK
library(ompr.roi)
library(ROI.plugin.glpk)
result <- solve_model(model, with_ROI(solver = "glpk", verbose = TRUE))
## <SOLVER MSG> ----
## GLPK Simplex Optimizer, v4.47
## 110 rows, 110 columns, 434 non-zeros
## 0: obj = 0.000000000e+000 infeas = 2.900e+001 (20)
## * 33: obj = 2.430553712e+003 infeas = 0.000e+000 (1)
## * 71: obj = 1.034476001e+003 infeas = 7.401e-017 (1)
## OPTIMAL SOLUTION FOUND
## GLPK Integer Optimizer, v4.47
## 110 rows, 110 columns, 434 non-zeros
## 100 integer variables, 90 of which are binary
## Integer optimization begins...
## + 71: mip = not found yet >= -inf (1; 0)
## + 315: >>>>> 1.773467222e+003 >= 1.061876113e+003 40.1% (27; 1)
## + 545: >>>>> 1.564760792e+003 >= 1.072548402e+003 31.5% (46; 6)
## + 1031: >>>>> 1.504476380e+003 >= 1.116118613e+003 25.8% (80; 30)
## + 1211: >>>>> 1.452500623e+003 >= 1.121648647e+003 22.8% (86; 42)
## + 2573: >>>>> 1.399547164e+003 >= 1.252722578e+003 10.5% (140; 115)
## + 3454: >>>>> 1.392498047e+003 >= 1.324965760e+003 4.8% (81; 274)
## + 3949: mip = 1.392498047e+003 >= tree is empty 0.0% (0; 567)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----
## To extract the solution we can use get_solution method that will return a data.frame
## which we can further be used with tidyverse packages.
solution <- get_solution(result, x[i, j]) %>%
filter(value > 0)
kable(head(solution, 10))
| x |
6 |
1 |
1 |
| x |
1 |
2 |
1 |
| x |
7 |
3 |
1 |
| x |
9 |
4 |
1 |
| x |
4 |
5 |
1 |
| x |
3 |
6 |
1 |
| x |
5 |
7 |
1 |
| x |
10 |
8 |
1 |
| x |
8 |
9 |
1 |
| x |
2 |
10 |
1 |
## Now we need to link back the indexes in our model with the actual Spots.
paths <- select(solution, i, j) %>%
rename(from = i, to = j) %>%
mutate(trip_id = row_number()) %>%
tidyr::gather(property, idx_val, from:to) %>%
mutate(idx_val = as.integer(idx_val)) %>%
inner_join(Weld_Spot, by = c("idx_val" = "id"))
kable(head(arrange(paths, trip_id),40))
| 1 |
from |
6 |
338.42244 |
421.57400 |
| 1 |
to |
1 |
352.90650 |
378.00491 |
| 2 |
from |
1 |
352.90650 |
378.00491 |
| 2 |
to |
2 |
497.95978 |
404.12478 |
| 3 |
from |
7 |
206.26002 |
268.15366 |
| 3 |
to |
3 |
205.19469 |
331.14037 |
| 4 |
from |
9 |
58.50193 |
98.34538 |
| 4 |
to |
4 |
96.04935 |
174.41220 |
| 5 |
from |
4 |
96.04935 |
174.41220 |
| 5 |
to |
5 |
19.47843 |
261.90337 |
| 6 |
from |
3 |
205.19469 |
331.14037 |
| 6 |
to |
6 |
338.42244 |
421.57400 |
| 7 |
from |
5 |
19.47843 |
261.90337 |
| 7 |
to |
7 |
206.26002 |
268.15366 |
| 8 |
from |
10 |
390.72731 |
313.50804 |
| 8 |
to |
8 |
182.11250 |
70.67535 |
| 9 |
from |
8 |
182.11250 |
70.67535 |
| 9 |
to |
9 |
58.50193 |
98.34538 |
| 10 |
from |
2 |
497.95978 |
404.12478 |
| 10 |
to |
10 |
390.72731 |
313.50804 |
## Plot Weld Spot and show shortest path planning
ggplot(Weld_Spot, aes(x, y)) +
geom_point() +
geom_line(data = paths, aes(group = trip_id)) +
ggtitle(paste0(" Robotized Welding Path Planning for part 6000: Optimal Spot Welding path
generated by vision sensors and CAD_File "))
## Generate welding spots using vision sensors
set.seed(7000)
Weld_Spot <- data.frame(id = 1:n, x = runif(n, max = max_x), y = runif(n, max = max_y))
ggplot(Weld_Spot, aes(x, y)) +
geom_point()
## Generate welding spots using vision sensors
set.seed(9000)
Weld_Spot <- data.frame(id = 1:n, x = runif(n, max = max_x), y = runif(n, max = max_y))
ggplot(Weld_Spot, aes(x, y)) +
geom_point()
