I want to solve a very simple quadratic optimization problem in R where one of the constraints is an equality constraint related to the sum of a vector. I tried to use the quadprog package but I cannot figure out how to make it work. The help function does not provide an example where the one of the constraints equals the sum. But perhaps quadprog is not the right package to use...
My problem is the following. I would like to shrink/expand a vector of values in such a way that the sum is equal to a set value and their values are bound by 0 and a maximum value. As the results need to be as close to the initial values as possible I want to minimize the residual sum of squares. Below I present a toy example. It is easy to see that for example (5,5,6) is in the solution space but not the optimal solution.
x_start <- c(4,3,7)
tot <- 16
# Objective function that needs to be minimized
obj <- function(x){return((x-x_start)^2)}
# Subject to:
# Sum constraint
con1 <- function(x) {return(sum(x)-tot)}
# Bounds 0 <= x <=6
x_max <- 6
x_min <- 0
There are a number of packages in R that can be used for this. See the Optimization Task View on CRAN: https://cran.r-project.org/web/views/Optimization.html
Here are two possibilities.
1) CVXR We can use CVXR for convex optimization.
library(CVXR)
x_start <- c(4,3,7)
tot <- 16
x <- Variable(3)
objective <- Minimize(sum((x - x_start)^2))
constraints <- list(sum(x) == tot, x >= 0, x <= 6)
problem <- Problem(objective, constraints)
soln <- solve(problem)
xval <- soln$getValue(x)
xval
## [,1]
## [1,] 5.500036
## [2,] 4.499964
## [3,] 6.000000
soln$value
## [1] 5.5
soln$status
## [1] "optimal"
2) limSolve Another package that can do this is limSolve. This solve the problem minimize ||Ax-B||^2 over the vector x subject to Ex = F and Gx >= H.
library(limSolve)
lsei(A = diag(3), B = x_start,
E = rep(1, 3), F = tot,
G = rbind(diag(3), -diag(3)), H = c(0, 0, 0, -6, -6, -6))
giving:
$X
[1] 5.5 4.5 6.0
$residualNorm
[1] 7.105427e-15
$solutionNorm
[1] 5.5
$IsError
[1] FALSE
$type
[1] "lsei"