2012-12-20 158 views
2

我一直在尝试解决R中使用constrOptim()(我的第一次)的约束优化问题,但我正在努力为我的问题设置约束。约束优化R建立约束

的问题是非常简单的,我可以设置功能确定,但​​在大约经过约束损失感到有点。

例如我已经定义的问题是(我打算以N固定在1000开始说,所以我只想解决X最终我想选择N和X的最大利润):

所以我可以设置功能为:

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)) 
} 

我需要实现的约束是:

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

的X这里表示成我需要最优地分配N个区段,所以N个X1 = pecent放置在桶1等,每个桶具有至少3%但不超过70%。

非常感谢任何帮助...

例如,这里是我用来测试功能的一个例子,我想要什么:

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); 
+0

所以...'X [1:3],'是变量,而一个'[1:3]','S [1:3]'和'C [1:3]'中给出值?你能否详细说明你需要的配方?我不清楚... – digEmAll

+0

是的 - a [1:3],s [1:3],c [1:3]代表历史平均激活率,支出,购置成本..我试图选择x [1:3]来最大限度地发挥 – andrewm4894

回答

2

您可以使用lpSolveAPI包。

## 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) 

打造LP模型解决的最好办法是纵列;

#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),">=","=")) 

看看我的模型

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 
+0

的功能,非常感谢,请看看我能否弄明白。你给出的规范中缺少N(例如,我刚刚有N = 1000,开始简单)。在我运行解决(my.lp)之后,我怎样才能找到x [1:3]的最佳值...??非常感谢 – andrewm4894

+0

@ user1919374是的N丢失,但你可以添加它。但更重要的是看你是否对决策变量的类型有所限制?是全部整数,例如二进制? – agstudy

+0

全部整数..... – andrewm4894

1

嗨只是在使用任何人的情况下,我也找到了答案如下:

首先,你的目标函数可以写成更简洁地使用矢量操作:

> 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) 

你正在尝试解决一个线性程序,所以你可以使用simplex算法来解决它。 'boot'包中有一个轻量级的实现。

> 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. 
+0

我也应该说这个问题,因为我已经制定了它总是有一个简单的解决方案,是尽可能分成受限制的'最佳'桶 - 事后很明显,这是很明显的! – andrewm4894