# MARKETING OPTIMIZATION R-SCRIPT
# Bob Agnew, raagnew@comcast.net, http://home.comcast.net/~raagnew
# 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.  Remarkably, this demo problem, with a million
# prospects and three offers, solves in about a minute on my home PC.  Previously,
# I solved it using SAS PROC NLP on a big server and the time required was much
# greater. Ref:
# http://home.comcast.net/~raagnew/Downloads/Dual_Solution_Marketing_Optimization.pdf
# 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. It is also possible
# to allow more than one offer per prospect.  R for Windows or Mac can be downloaded
# free from www.r-project.org. This script can be opened, highlighted, and run from
# the R console.  It can also be run essentially in batch using the source function,
# e.g., source("c:/Optimization/R/mktg_opt.R").  
# Number of Prospects
n <- 1000000
# Initialize prospect offer and profit at zero, i.e., no offer.
Profit <- Offer <- array(0,dim = 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
b <- c(300000,200000,100000)
# Dual Function
dual <- function(u) {
d1 <- p1 - u[1]; d2 <- p2 - u[2]; d3 <- p3 - u[3]
v <- (d1>=0&d1>=d2&d1>=d3)*d1 + (d2>=0&d2>d1&d2>=d3)*d2 + (d3>=0&d3>d1&d3>d2)*d3
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 <- (d1>=d2&d1>=d3)*d1 + (d2>d1&d2>=d3)*d2 + (d3>d1&d3>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]
e1<-k[1]*d1[m]-(1-k[1])*h;e2<-k[2]*d2[m]-(1-k[2])*h;e3<-k[3]*d3[m]-(1-k[3])*h
e <- c(e1,e2,e3)
i <- order(e)[3]
Offer[m] <- i
p <- c(p1[m],p2[m],p3[m])
Profit[m] <- p[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]))
# In actual application, you would export a text file containing the optimal
# prospect offers, i.e., write(Offer, file = "outfile", ncolumns = 1).