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)
}
CodePudding user response:
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
CodePudding user response:
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