constrained optimization in R setting up constraints

Question

I have been trying to solve a constrained optimization problem in R using constrOptim() (my first time) but am struggling to set up the constraints for my problem.

The problem is pretty straight forward and i can set up the function ok but am a bit at a loss about passing the constraints in.

e.g. problem i've defined is (am going to start with N fixed at 1000 say so i just want to solve for X ultimately i'd like to choose both N and X that max profit):

so i can set up the function as:

fun <- function(x, N, a, c, s) {   ## a profit function
    x1 <- x[1]
    x2 <- x[2]
    x3 <- x[3]    
    a1 <- a[1]
    a2 <- a[2]
    a3 <- a[3]
    c1 <- c[1]
    c2 <- c[2]
    c3 <- c[3]
    s1 <- s[1]
    s2 <- s[2]
    s3 <- s[3]  
    ((N*x1*a1*s1)-(N*x1*c1))+((N*x2*a2*s2)-(N*x2*c2))+((N*x3*a3*s3)-(N*x3*c3))
}

The constraints i need to implement are that:

x1>=0.03
x1<=0.7
x2>=0.03
x2<=0.7
x3>=0.03
x2<=0.7
x1+x2+x3=1

The X here represents buckets into which i need to optimally allocate N, so x1=pecent of N to place in bucket 1 etc. with each bucket having at least 3% but no more than 70%.

Any help much appreciated...

e.g. here is an example i used to test the function does what i want:

    fun <- function(x, N, a, c, s) {   ## profit function
        x1 <- x[1]
        x2 <- x[2]
        x3 <- x[3]    
        a1 <- a[1]
        a2 <- a[2]
        a3 <- a[3]
        c1 <- c[1]
        c2 <- c[2]
        c3 <- c[3]
        s1 <- s[1]
        s2 <- s[2]
        s3 <- s[3]  
        ((N*x1*a1*s1)-(N*x1*c1))+((N*x2*a2*s2)-(N*x2*c2))+((N*x3*a3*s3)-(N*x3*c3))
    };

    x <-matrix(c(0.5,0.25,0.25));

    a <-matrix(c(0.2,0.15,0.1));

    s <-matrix(c(100,75,50));

    c <-matrix(c(10,8,7));

    N <- 1000;

fun(x,N,a,c,s);

So... x[1:3] are variables, while a[1:3], s[1:3] and c[1:3] are given values ? Could you elaborate the formulation you need ? It is not clear to me... — digEmAll
– digEmAll, Commented Dec 20, 2012 at 16:52
yes exactly - a[1:3], s[1:3], c[1:3] represent historic average activation rate, spend, acquisition cost..i'm trying to choose the x[1:3] to maximize the function — andrewm4894
– andrewm4894, Commented Dec 20, 2012 at 16:56

agstudy · Accepted Answer · 2012-12-20 17:47:54Z

2

You can use The lpSolveAPI package.

## problem constants
a <- c(0.2, 0.15, 0.1)
s <- c(100, 75, 50)
c <- c(10, 8, 7)
N <- 1000


## Problem formulation
# x1          >= 0.03
# x1          <= 0.7
#     x2      >= 0.03
#     x2      <= 0.7
#          x3 >= 0.03
# x1 +x2 + x3  = 1
#N*(c1- a1*s1)* x1 + (a2*s2 - c2)* x2 + (a3*s3-  c3)* x3

library(lpSolveAPI)
my.lp <- make.lp(6, 3)

The best way to build a model in lp solve is columnwise;

#constraints by columns
set.column(my.lp, 1, c(1, 1, 0, 0, 1, 1))
set.column(my.lp, 2, c(0, 0, 1, 1, 0, 1))
set.column(my.lp, 3, c(0, 0, 0, 0, 1, 1))
#the objective function ,since we need to max I set negtive max(f) = -min(f)
set.objfn (my.lp, -N*c(c[1]- a[1]*s[1], a[2]*s[2] - c[2],a[3]*s[3]-  c[3]))
set.rhs(my.lp, c(rep(c(0.03,0.7),2),0.03,1))
#constraint types 
set.constr.type(my.lp, c(rep(c(">=","<="), 2),">=","="))

take a look at my model

my.lp
Model name: 

 Model name: 
             C1     C2     C3          
Minimize  10000  -3250   2000          
R1            1      0      0  >=  0.03
R2            1      0      0  <=   0.7
R3            0      1      0  >=  0.03
R4            0      1      0  <=   0.7
R5            1      0      1  >=  0.03
R6            1      1      1   =     1
Kind        Std    Std    Std          
Type       Real   Real   Real          
Upper       Inf    Inf    Inf          
Lower         0      0      0      
 solve(my.lp)

[1] 0    ## sucess :)

get.objective(my.lp)
[1] -1435
get.constraints(my.lp)
[1] 0.70 0.70 0.03 0.03 0.97 1.00
## the decisions variables
get.variables(my.lp)
[1] 0.03 0.70 0.27

edited Dec 20, 2012 at 17:47

answered Dec 20, 2012 at 17:12

agstudy

122k18 gold badges205 silver badges265 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

andrewm4894 Over a year ago

that's great thanks, will take a look to see if i can figure it out. Is N missing from the specification you gave (e.g. i just had N=1000 to start simple). also how can i find the optimum values for x[1:3] after i run solve(my.lp)...? many thanks

agstudy Over a year ago

@user1919374 Yes The N is missing but you can add it. But the more important is to see if you have constraints on the type of your decisions variables? are all integer ,binary for example?

agstudy Over a year ago

@user1919374 I add N and the decisions variables. But you need to integrate the decisions variables types, e.g sset.type(my.lp, c(1,2,3), "integer").

andrewm4894 · Accepted Answer · 2012-12-21 11:11:39Z

1

Hi Just in case of use to anyone i also found an answer as below:

First of all, your objective function can be written a lot more concisely using vector operations:

> my_obj_coeffs <- function(N,a,c,s) N*(a*s-c)


> fun <- function(x,N,a,c,s) sum(my_obj_coeffs(N,a,c,s) * x)

You're trying to solve a linear program, so you can use solve it using the simplex algorithm. There's a lightweight implementation of it in the 'boot' package.

> library(boot)

> solution <- function(obj) simplex(obj, diag(3), rep(0.7,3), diag(3), rep(0.03,3), rep(1,3), 1, maxi=TRUE)

Then for the example parameters you used, you can call that solution function:

> a <- c(0.2,0.15,0.1)
> s <- c(100,75,50)
> c <- c(10,8,7)
> N <- 1000
> solution(my_obj_coeffs(N,a,c,s))

Linear Programming Results

Call : simplex(a = obj(N, a, s, c), A1 = diag(3), b1 = rep(0.7, 3), 
    A2 = diag(3), b2 = rep(0.03, 3), A3 = matrix(1, 1, 3), b3 = 1, 
    maxi = TRUE)

Maximization Problem with Objective Function Coefficients
      [,1]
[1,] 10000
[2,]  3250
[3,] -2000
attr(,"names")
[1] "x1" "x2" "x3"


Optimal solution has the following values
  x1   x2   x3 
0.70 0.27 0.03 
The optimal value of the objective  function is 7817.5.

answered Dec 21, 2012 at 11:11

andrewm4894

1,5584 gold badges17 silver badges40 bronze badges

1 Comment

andrewm4894 Over a year ago

i should also say the problem as i have formulated it always has a simple solution which is to portion as much as possible into the 'best' buckets subject to the constraints - which is pretty obvious in hindsight!

Collectives™ on Stack Overflow

constrained optimization in R setting up constraints

2 Answers 2

3 Comments

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related