0

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
Michiel van Dijk
  • 157
  • 1
  • 9

2 Answers2

1

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"
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
1

Try to express your quadratic problem (i.e. the objective and constraints) in matrix and vector form. For the sum constraint, you can use the fact that the sum of vectors is the dot product with the vector of ones, so the corresponding row on the matrix of constraints is c(1,1,1) and the rhs for that constraint is tot.

According to the documentation, the meq parameter in solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE) specifies the number of constraints that are equality constraints,

the first meq constraints are treated as equality constraints

so change this parameter to 1 (or however many equality constraints you need) and make the first column of Amat a vector of ones.

Juan Carlos Ramirez
  • 2,054
  • 1
  • 7
  • 22