Optimize a vector when constraint coefficients change with decision variables

coolhand 0
Viewed 84

Say I have a state of current scores held in a vector:

current_score <- c(60, 59, 78, 79, 76)

I want to decide what scores to increase. I can increase them to some value (say 95) if I choose that index. Each has an associated cost if chosen.

decision <- rep(1, times=length(current_score))
cost <- c(2792, 5207, 1814, 642, 1037)

I also want to have the weighted score/cost of the entire vector to meet some minimum goal.

goal <- 70

How can I find the optimal indexes to choose?

Here's the function I've tried using nonlinear optimization with optim but it always returns 0. I'm also not sure how I can implement it with binary constraints (decision vector should be either 1 or 0).

fn <- function(decision, current_score, cost, goal)
{
    # if decision = TRUE, score = 95
    # this would provide the constraint coefficients if done in lpSolve
    current_score[which(decision == TRUE)] <- 95

    # calculate weighted score after decision
    wt_score <- sum(current_score*cost)/sum(cost)   
    
    # if decision doesn't meet goal, assign arbitrary large value to reject this choice
    if(wt_score < goal)
    {
        decision_cost <- .Machine$double.xmax
    } else
    {
        decision_cost <- sum(decision * cost)
    } 
    
    return(decision_cost)
}

I would normally set this up as a linear program with lpSolve but the constraint coefficients change along with the decision variables. I'm not sure if I can implement the binary decision constraints with optim. With the approach below, it returns .Machine$double.xmax which is obviously incorrect.

upper <- rep(1, times=length(decision))
lower <- rep(0, times=length(decision))

result <- optim(par=decision, fn=fn, upper=upper, lower=lower, method=c("L-BFGS-B"), current_score = current_score, cost = cost, goal=goal)

For example, decision = c(1,0,0,0,0) provides wt_score = 73.39 and decision_cost = 2792 so it's a better solution.

Edit: FWIW, this is feasible using a genetic algorithm within Excel, but I would like to do it in R

Below it's mentioned this can be set up as a linear problem. The issue I have is how to implement the constraints because they change as the decision variables change. E.g.,

# return a vector that is changed based on the decision variables
constraint_vector <- function(score, x)
{
  score <- score[x==TRUE] <- 95
}

lp_min <- function(score, cost, crvs, goal)
{
  # create a repair vector
  f.obj <- cost
  # not constant, but a function of decision variable x
  f.con <- constraint_vector(score,x)
  f.dir <- c(">=")
  f.rhs <- goal
  
  solution <- lp("in", f.obj, f.con, f.dir, f.rhs, all.bin=TRUE)
}
optimization r
2 Answers

I think a mathematical model can look like:

decision variables: x(i): select score i (binary variable)
                    s(i): final score (continuous variable)
data:               cs(i): current (old) score
                    c(i): cost

min cost = sum(i, c(i)*x(i))
    s(i) = cs(i)*(1-x(i)) + 95*x(i)        (1)
    sum(i, s(i)*c(i)) >= 70*sum(i,c(i))    (2)
    x(i) ∈ {0,1}                           (3)

This is a linear MIP (Mixed-Integer Programming) Model. This can be implemented and solved using any MIP solver. Even LpSolve.

For the programmers who specialize in copy-paste:

current_score <- c(60, 59, 78, 79, 76)  # better name: initial score
cost <- c(2792, 5207, 1814, 642, 1037)  # coefficients
goal <- 70                              # better name: minimum weighted score


library(lpSolve)
n <- length(cost) # number of x(i) variables (and number of s(i) variables)
N <- 2*n # total number of variables (x and s)
M <- 2*n+1 # total number of constraints (binary variables need extra constraints in this interface!!!)

# columns: x first, then s
# rows: constraints (1),(2) and (3) in that order
A <- matrix(0,nrow=M,ncol=N)
for (i in 1:n) {
  A[i,i]       = current_score[i]-95;
  A[i,n+i]     = 1
  A[n+1,n+i]   = cost[i]
  A[n+1+i,i]   = 1
}

rhs <- c(current_score,goal*sum(cost),rep(1,n))
signs <- c(rep("=",n),">=",rep("<=",n))
ints <- 1:n
objs <- c(cost,rep(0,n))
res <- lp("min",objs,A,signs,rhs,T,ints)

res
data.frame(x=res$solution[1:n],s=res$solution[(n+1):(2*n)])

This prints:

> res
Success: the objective function is 2792 

> data.frame(x=res$solution[1:n],s=res$solution[(n+1):(2*n)])
  x  s
1 1 95
2 0 59
3 0 78
4 0 79
5 0 76

Notes:

  1. There are much better tools and solvers for expressing and solving these types of models.
  2. When implementing this in a particular environment, focus on the mathematical model. This defines the model we want to solve. And it is easier to understand and to reason about than a bunch of code. I.e. the mathematical model is worth more than code, and should be all you need. Implementation in a (programming or other) language is precise, detailed work, but conceptually simple. The heavy lifting is already done. However, more than once, I am disappointed in the dexterity of programmers with this. Way too much hand-holding is being asked for. May be an education problem. I don't know.
Erwin Kalvelagen0

if a decision can only be combination of 5 binary choices, there are only 32 combinations to search exhaustively. evaluate all of them explicitly. use this to get started :

do.call(expand.grid,rep(list(0:1),5))

p.s. i believe you mislabled the argument score of your function its never referenced outside of the function definition

Nir Graham0
Related