# MARKETING OPTIMIZATION R-SCRIPT
# Bob Agnew, raagnew1@gmail.com, www.raagnew.com
# The algorithm imbedded in this script optimizes the assignment of offers to
# prospects, subject to stipulated offer quantity constraints and the limit of
# at most one offer per prospect.  This demo problem, with a million
# prospects and three offers, solves quickly on my home PC.  This script can be
# adapted for a different number of offers and alternative equality or inequality
# constraints, as long as they can be linearized, although the dual must be
# formulated for solution.  Inequality constraints would require a different
# R solver (nlminb) which incorporates lower bounds.
#
# Offer vector contains the optimal offers; zero indicates no offer to a prospect
#
date() # Start time stamp
# Number of Prospects
n <- 1000000
# Initialize prospect offer and profit at zero, i.e., no offer.
Profit <- Offer <- rep(0,n)
# Simulate Prospect Profits for Three Offers
# Different profit vectors will be generated for each execution, which has the
# advantage of showing that solution is robust across various instances.
p1 <- 100*runif(n,0,1) - 10
p2 <- .6*p1 + .4*(60*runif(n,0,1) - 6)
p3 <- .4*p1 + .6*(40*runif(n,0,1) - 4)
# Ordinarily, offer profit "scores" would be input from an external text file using the
# scan function, i.e., inp <- scan("infile", list(0,0,0)); p1 <- inp[[1]];
# p2 <- inp[[2]]; p3 <- inp[[3]].  You can check this out by saving an instance of the
# simulated vectors via write(rbind(p1,p2,p3), file = "infile", ncolumns = 3) and then
# use a different version of this script with scan rather than simulation. 
# Mean Offer Profits
c(mean(p1),mean(p2),mean(p3))
# Stipulated Offer Quantities; 400000 get no offer
b <- c(300000,200000,100000)
# Dual Function
dual <- function(u) {
d1 <- p1 - u[1]; d2 <- p2 - u[2]; d3 <- p3 - u[3]
v <- pmax(d1,d2,d3,rep(0,n))
y <- b%*%u + sum(v)
y
}
# Dual Minimization using Nonlinear Minimization Function
out <- nlm(dual,c(0,0,0))
# Dual Minimization Output
out
# Keep dual minimum and estimates.
mindual <- out$minimum
u <- out$estimate
d1 <- p1 - u[1]; d2 <- p2 - u[2]; d3 <- p3 - u[3]
v <- pmax(d1,d2,d3)
ord <- order(v, decreasing = TRUE)
s <- c(0,0,0)
h <- 100000000
for (j in 1:sum(b)) {
k <- s < b
m <- ord[j]
i <- which.max(c(k[1]*d1[m]-(1-k[1])*h,k[2]*d2[m]-(1-k[2])*h,k[3]*d3[m]-(1-k[3])*h))
Offer[m] <- i
Profit[m] <- c(p1[m],p2[m],p3[m])[i]
s[i] <- s[i] + 1
}
# Dual Minimum (Upper Bound)
mindual
# Primal Optimum (Maximum Profit)
sum(Profit)
# Note that we are extremely close, in some instances right on.
# Offer Quantities
s
# Offer Profits
c(sum(Profit[Offer == 1]),sum(Profit[Offer == 2]),sum(Profit[Offer == 3]))
# Offers to first 100 prospects
Offer[1:100]
# In actual application, you would export a text file containing the optimal
# prospect offers, i.e., write(Offer, file = "outfile", ncolumns = 1).
date() # End time stamp