CODEWITHSUNDEEP
r
Rich Garrod
Rich Garrod

Reputation: 25

Constrain Optimisation Problems in R

I am trying to set up an optimisation script that will look at a set of models, fit curves to the models and then optimise across them, subject to a few parameters.

Essentially, I have revenue as a function of cost, in a diminishing function, and I have this for multiple portfolios, say 4 or 5. As an input, I have cost and revenue figures, at set increments. What I want to do is fit a curve to the portfolio of the form Revenue=A*cost^B, and then optimise across the different portfolios to find the optimal cost split between each portfolio for a set budget.

The code below (I apologise for the inelegance of it, I'm sure there are MANY improvements to be made!) essentially reads in my data, in this case, a simulation, creates the necessary data frames (this is likely where my inelegance comes in), calculates the necessary variables for the curves for each simulation and produces graphics to check the fitted curve to the data.

My problem is that now I have 5 curves of the form:

revenue = A * Cost ^ B (different A, B and cost for each function)

And I want to know, given the 5 variables, how should I split my cost between them, so I want to optimise the sum of the 5 curves subject to

Cost <= Budget

I know that I need to use constrOptim, but I have spent literally hours banging my head against my desk (literally hours, not literally banging my head...) and I still can't figure out how to set up the function so that it maximises revenue, subject to the cost constraint...

Any help here would be greatly appreciated, this has been bugging me for weeks.

Thanks!

Rich

## clear all previous data

rm(list=ls())
detach()
objects()

library(base)
library(stats)

## read in data

sim<-read.table("input19072011.txt",header=TRUE)
sim2<-data.frame(sim$Wrevenue,sim$Cost)

## identify how many simulations there are - here you can change the 20 to the number of steps but all simulations must have the same number of steps

portfolios<-(length(sim2$sim.Cost)/20)

## create a matrix to input the variables into

a<-rep(1,portfolios)
b<-rep(2,portfolios)
matrix<-data.frame(a,b)

## create dummy vector to hold the revenue predictions

k<-1
j<-20

for(i in 1:portfolios){

test<-sim2[k:j,]

rev9<-test[,1]
cost9<-test[,2]

ds<-data.frame(rev9,cost9)

rhs<-function(cost, b0, b1){
b0 * cost^b1

m<- nls(rev9 ~ rhs(cost9, intercept, power), data = ds, start = list(intercept = 5,power = 1))

matrix[i,1]<-summary(m)$coefficients[1]
matrix[i,2]<-summary(m)$coefficients[2]

k<-k+20
j<-j+20

}

## now there exists a matrix of all of the variables for the curves to optimise

matrix
multiples<-matrix[,1]
powers<-matrix[,2]
coststarts<-rep(0,portfolios)

## check accuracy of curves

k<-1
j<-20

for(i in 1:portfolios){

dev.new()

plot(sim$Wrevenue[k:j])
lines(multiples[i]*(sim$Cost[k:j]^powers[i]))

k<-k+20
j<-j+20

}

Upvotes: 2

Views: 689

Answers (1)

Vincent Zoonekynd
Vincent Zoonekynd

Reputation: 32351

If you want to find the values cost[1],...,cost[5] that maximize revenue[1]+...+revenue[5] subject to the constraints cost[1]+...+cost[5]<=budget (and 0 <= cost[i] <= budget), you can parametrize the set of feasible solutions as follows

cost[1] = s(x[1]) * budget
cost[2] = s(x[2]) * ( budget - cost[1] )
cost[3] = s(x[3]) * ( budget - cost[1] - cost[2])
cost[4] = s(x[4]) * ( budget - cost[1] - cost[2] - cost[3] )
cost[5] = budget - cost[1] - cost[2] - cost[3] - cost[4]

where x[1],...,x[4] are the parameters to find (with no constraints on them) and s is any bijection between the real line R and the segment (0,1).

# Sample data
a <- rlnorm(5)
b <- rlnorm(5)
budget <- rlnorm(1)

# Reparametrization
s <- function(x) exp(x) / ( 1 + exp(x) )
cost <- function(x) {
  cost <- rep(NA,5)
  cost[1] = s(x[1]) * budget
  cost[2] = s(x[2]) * ( budget - cost[1] )
  cost[3] = s(x[3]) * ( budget - cost[1] - cost[2])
  cost[4] = s(x[4]) * ( budget - cost[1] - cost[2] - cost[3] )
  cost[5] = budget - cost[1] - cost[2] - cost[3] - cost[4]
  cost  
}

# Function to maximize
f <- function(x) {
  result <- sum( a * cost(x) ^ b )
  cat( result, "\n" )
  result
}

# Optimization
r <- optim(c(0,0,0,0), f, control=list(fnscale=-1))
cost(r$par)

Upvotes: 5

Related Questions