# Import lpSolve package
library(lpSolve)
## Warning: package 'lpSolve' was built under R version 4.1.1
# Set coefficients of the objective function
f.obj <- c(5, 7)
 
# Set matrix corresponding to coefficients of constraints by rows
# Do not consider the non-negative constraint; it is automatically assumed
f.con <- matrix(c(1, 0,
                  2, 3,
                  1, 1), nrow = 3, byrow = TRUE)

# Set unequality signs
f.dir <- c("<=",
           "<=",
           "<=")

# Set right hand side coefficients
f.rhs <- c(16,
           19,
           8)

# Final value (z)
lp("max", f.obj, f.con, f.dir, f.rhs)
## Success: the objective function is 46
# Variables final values
lp("max", f.obj, f.con, f.dir, f.rhs)$solution
## [1] 5 3
# Sensitivities
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$sens.coef.from
## [1] 4.666667 5.000000
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$sens.coef.to
## [1] 7.0 7.5
# Dual Values (first dual of the constraints and then dual of the variables)
# Duals of the constraints and variables are mixed
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals
## [1] 0 2 1 0 0
# Duals lower and upper limits
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals.from
## [1] -1.000000e+30  1.600000e+01  6.333333e+00 -1.000000e+30 -1.000000e+30
lp("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals.to
## [1] 1.0e+30 2.4e+01 9.5e+00 1.0e+30 1.0e+30

Daftar Pustaka : https://towardsdatascience.com/linear-programming-in-r-444e9c199280