Optimal Allocation Problems - Part I
Authors@R: person(“Annice”, “Najafi”, email = “annicenajafi27@gmail.com”)
The following R code is related to the special case of the problem - “Optimal Allocations In the Construction of k-out-of-n Reliability Systems” - by Derman et al, 1974.
References:
Derman, C., Lieberman, G. J., & Ross, S. M. (1974). Optimal allocations in the construction of k-out-of-n reliability systems. Management Science, 21(3), 241-250.
I am going to start this section by considering the most basic problem which is the situation where we are considering money that gets allocated for an n component system that is getting built. Assume the case where we need to have at least one successful component for the system to work.
Case k=1
Let’s think about this, if we are allocating money simultaneously then we are essentially trying to maximize the probability of one component being successful. That would be \(1-(1-p(x_1))(1-p(x_2))...(1-p(x_n))\) (\(1 - \prod_1^n(1-p(x_i))\)) with \(x_i\) being the money allocated at the ith stage and \(i\) from 1 to n (total number of stages).
The probability of success of the component depends on how much we invest. If we invest 0 dollars, the probability of success is 0 meaning p(0)=0.
Now consider the case where we are unaware of whether a component will be successful or not until it is built. In that case we start with the first component and the probability of success of the component would be p(x_1). If it is unsuccessful then we move on to the next component so its probability would be conditioned on the first component being unsuccessful which is \(p(x_1)(1-p(x_2))\) and so on (\(p(x_1) + \sum_{i=2}^n \prod_{j=1}^{i-1} (1 - p(x_j)) p(x_i)\)).
Below we are going to answer the question of whether investing all of our money into the first task is optimal or dividing our money into all tasks equally and invest in each of them. We take two cases where p(x) is concave and convex and then proceed by increasing the number of components we would need to be successful (k).
Load relevant libraries
library(Deriv)
library(ggplot2)Concave p(x) function
# Define the function for probability of success
P <- function(x) {
alpha <- 0.1
return((exp(alpha * x) - 1) / (exp(alpha * x) + 1))
}
#If you take the log of the equation we talked about above for the
#non-sequential case
alpha <- 0.1
log_expr <- expression(log(1 - (1 - (exp(alpha * x) - 1) / (exp(alpha * x) + 1))))
# and find the first derivative of log(1 - P(x)) with respect to x
first_derivative <- Deriv(log_expr, "x")
# and then the second derivative of log(1 - P(x)) with respect to x
# to check if it is concave or convex - remember that a function is
# concave if its second derivative is negative and convex if positive
second_derivative <- Deriv(first_derivative, "x")
x <- seq(0, 100, 1) #assume the amount we are investing is between 0 and 100
evaluated_second_deriv <- eval(second_derivative)
#Which we see here that the function is actually concave
ggplot(data.frame(x=x, y=P(x)), aes(x=x, y=y))+geom_point(stroke=2,shape=8, color = "#D96098")+
geom_line(size=2, color = "#D96098")+
labs(x="x", y="p(x)")+
ggtitle("p(x) is non-decreasing in x")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))+labs(tag="A")->plt.A## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plt.A
ggplot(data.frame(x=x, y=evaluated_second_deriv), aes(x=x, y=y))+
geom_point(stroke=2,shape=8, color = "#D96098")+
geom_line(size=2, color = "#D96098")+
geom_hline(yintercept=0, linetype="dashed", size=2)+
labs(x="x", y="p\"(x)")+
ggtitle("p(x) is a concave function")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))+labs(tag="B")->plt.B
plt.B
Define our strategies
Again, simulating the success of the outcome is like flipping a coin that has a p(x) probability of getting heads, if p(x) is higher then the probability of getting heads is higher.
#Let's say I am taking the strategy of allocating A/N for each component,
#essentially I am dividing my total money (A) over all N component
#' Simulate problem under strategy 1
#'
#' @param A all of the money we are investing
#' @param n the number of stages
#'
#' @return whether that one task is completed successfully or not
#' @export
simulate_problem_strategy_1 <- function(A, n){
count.success <- 0
for(i in 1:n){
if(rbinom(1, 1, P(A/n)) == 1){
count.success <- count.success + 1
}
}
if(count.success > 1){
return(1)
}else{
return(0)
}
}
#Second strategy is to put all of money down for the first task
#' Simulate problem under strategy 2
#'
#' @param A all of the money we are investing
#' @param n the number of stages
#'
#' @return whether that one task is completed successfully or not
#' @export
simulate_problem_strategy_2 <- function(A, n){
if(rbinom(1, 1, P(A)) == 0){
return(0)
}else{
return(1)
}
}Optimal strategy for a concave p(x)
A <- 10
n <- 100
num_simulations <- 1000
numeric(num_simulations) -> strategy_1.rewards
numeric(num_simulations) -> strategy_2.rewards
for(k in 1:num_simulations){
simulate_problem_strategy_1(A, n) -> strategy_1.rewards[k]
simulate_problem_strategy_2(A, n) -> strategy_2.rewards[k]
}
ggplot(data.frame(name=c("Equally distributed", "All in 1"),
success_rate = c(sum(strategy_1.rewards)/num_simulations, sum(strategy_2.rewards)/num_simulations)),
aes(x = name, y = success_rate))+
geom_bar(stat = "identity", width=0.2, fill= "#D96098")+
labs(x="Strategy", y="Success Rate")+
ggtitle("Optimal solution")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))
Let’s try a case where p(x) is a convex function
Convex p(x) function
# Define the function for probability of success
P <- function(x) {
return(0.01 * sqrt(x))
}
log_expr <- expression(log(1 - 0.01 * sqrt(x)))
# and find the first derivative of log(1 - P(x)) with respect to x
first_derivative <- Deriv(log_expr, "x")
# and then the second derivative of log(1 - P(x)) with respect to x
# to check if it is concave or convex
second_derivative <- Deriv(first_derivative, "x")
x <- seq(0, 100, 1)
evaluated_second_deriv <- eval(second_derivative)
#Which we see here that the function is actually concave
ggplot(data.frame(x=x, y=P(x)), aes(x=x, y=y))+geom_point(stroke=2,shape=8, color = "#D96098")+
geom_line(size=2, color = "#D96098")+
labs(x="x", y="p(x)")+
ggtitle("p(x) is non-decreasing in x")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))+labs(tag="A")->plt.A
plt.A
ggplot(data.frame(x=x, y=evaluated_second_deriv), aes(x=x, y=y))+
geom_point(stroke=2,shape=8, color = "#D96098")+
geom_line(size=2, color = "#D96098")+
geom_hline(yintercept=0, linetype="dashed", size=2)+
labs(x="x", y="p\"(x)")+
ggtitle("p(x) is a convex function")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))+labs(tag="B")->plt.B
plt.B
Optimal strategy for a convex p(x)
A <- 100
n <- 100
num_simulations <- 1000
numeric(num_simulations) -> strategy_1.rewards
numeric(num_simulations) -> strategy_2.rewards
for(k in 1:num_simulations){
simulate_problem_strategy_1(A, n) -> strategy_1.rewards[k]
simulate_problem_strategy_2(A, n) -> strategy_2.rewards[k]
}
ggplot(data.frame(name=c("Equally distributed", "All in 1"),
success_rate = c(sum(strategy_1.rewards)/num_simulations, sum(strategy_2.rewards)/num_simulations)),
aes(x = name, y = success_rate))+
geom_bar(stat = "identity", width=0.2, fill= "#D96098")+
labs(x="Strategy", y="Success Rate")+
ggtitle("Optimal solution")+
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size = 15),
text = element_text(size=18))
k>1
Now let’s increase the number of successful components needed for the system to work to k…
#' Simulate problem under strategy 1
#'
#' @param A all of the money we are investing
#' @param n stages (total components)
#' @param k components we would like to have successfully built
#'
#' @return whether k successful components are built or not
#' @export
simulate_problem_strategy_1 <- function(A, n, k){
count.success <- 0
for(i in 1:n){
if(rbinom(1, 1, P(A/n)) == 1){
count.success <- count.success + 1
}
}
if(count.success >= k){
return(1)
}else{
return(0)
}
}
#Second strategy is to put all of money down for the first k tasks
#' Simulate problem under strategy 2
#'
#' @param A all of the money we are investing
#' @param n stages (total components)
#' @param k components we would like to have successfully built
#'
#' @return whether k successful components are built or not
#' @export
simulate_problem_strategy_2 <- function(A, n, k){
count.success <- 0
for(i in 1:k){
if(rbinom(1, 1, P(A/k)) == 1){
count.success <- count.success + 1
}
}
if(count.success >= k){
return(1)
}else{
return(0)
}
}Simulate problem for k>1 successful components needed
A <- 10000
n <- 50
num_simulations <- 1000
numeric(num_simulations) -> strategy_1.rewards
numeric(num_simulations) -> strategy_2.rewards
storing.strategy.1.sums <- numeric(10)
storing.strategy.2.sums <- numeric(10)
k.index <- 1
for(k in seq(1, n, 2)){
for(p in 1:num_simulations){
simulate_problem_strategy_1(A, n, k) -> strategy_1.rewards[p]
simulate_problem_strategy_2(A, n, k) -> strategy_2.rewards[p]
}
sum(strategy_1.rewards)->storing.strategy.1.sums[k.index]
sum(strategy_2.rewards)->storing.strategy.2.sums[k.index]
k.index <- k.index + 1
}Comparison of strategies
df.k <- data.frame(k = seq(1, 50, 2), strategy_1 = storing.strategy.1.sums/num_simulations, strategy_2 = storing.strategy.2.sums/num_simulations)
ggplot() +
geom_line(data = df.k, aes(x = k, y = strategy_1, colour = "Equally distributed"), size = 2) +
geom_point(data = df.k, aes(x = k, y = strategy_1, colour = "Equally distributed"), shape = 8, stroke = 2) +
geom_line(data = df.k, aes(x = k, y = strategy_2, colour = "All in k"), size = 2) +
geom_point(data = df.k, aes(x = k, y = strategy_2, colour = "All in k"), shape = 8, stroke = 2) +
labs(x = "k (n=50)", y = "Success Rate", colour = "Strategy") +
ggtitle("Optimal solution") +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5, size = 20),
axis.text = element_text(size = 15),
text = element_text(size = 18)
) +
scale_color_manual(values = c("Equally distributed" = "#D96098", "All in k" = "#BEAEE2"))
Interesting! We see that the first strategy (plotted in pink) is optimal even if we increase the number of successful components k when p(x) is convex.
Coming soon…
Part II of these series for more advances problems on optimal allocation.