constrained optimization of a complicated function

Question

I need to optimize the following to find the maximum value for r1:

ad = 0.95*M_D + 0.28*G_D + 0.43*S_D + 2.25*Q_D
as = 0.017*M_A + 0.0064*G_A + 0.0076*S_A + 0.034*Q_A
ccb = 0.0093*M_CC+ 0.0028*G_CC + 0.0042*S_CC + 0.0186*Q_CC
ccd = 0.0223*M_CD + 0.0056*G_CD + 0.0078*S_CD + 0.0446*Q_CD
apb = 1.28*M_P + 2.56*Q_P 
r1=(1+ccb*(1+ccd))*ad*as*100/(130-apb)

subject to the following constraints:

0 <= M_CD <= M_CC <= M_A <= M_D <= M_P <= 9
0 <= G_CD <= G_CC <= G_A <= G_D <= 9
0 <= S_CD <= S_CC <= S_A <= S_D <= 9
0 <= Q_CD <= Q_CC <= Q_A <= Q_D <= Q_P <= 3

The approach I've tried before doesn't work very well and I'm hoping to find a better solution.

[Task View "Optimization"](http://cran.r-project.org/web/views/Optimization.html) — Roland, Jun 07 '13 at 07:36
Is there an error here? `r1=(1+ccb*(1+ccd))*ad*ad*100/(130)` no `as` and no `apb` — user1609452, Jun 07 '13 at 07:46
as it stands the solution is the upper bound for all variables — user1609452, Jun 07 '13 at 08:18
Can you provide some background here? I find it hard to believe there's a real-world problem with that many input variables and that many semi-independent constraints (especially since, as already noted, the maximum is obviously at max(all inputs) ). PS @Roland, great reference page but sadly not useful to a naive OP. — Carl Witthoft, Jun 07 '13 at 11:50
Yeah sorry user it's `r1=(1+ccb*(1+ccd))*ad*as*100/(130-apb)` I made some this mistake when simplifying my problem — Wicelo, Jun 07 '13 at 17:55
Again maximising ccb, ccd, ad, as, apb would seem to be the solution leading to the same solution upper bound for all variables. — user1609452, Jun 07 '13 at 18:27

user1609452 · Answer 1 · 2013-06-07T18:37:23.547

Once the problem is stated correctly you maybe able to firstly map the parameters to lower and upper bounds of [0,1]. You can then implement the inequalities inside your function and optimise using an algorithm which accepts basic lower and upper bound constraints. nlminb could be used but the vignette suggests the algorithm used may not be the best.

UPDATE:

With OP revised function

dumFun <- function(p){
    p[1] -> M_CD; p[2] -> M_CC; p[3] -> M_A; p[4] -> M_D; p[5] -> M_P;
    M_P <- 9*M_P; M_D <- M_P*M_D; M_A <- M_D*M_A; M_CC <- M_A*M_CC; M_CD <- M_CC*M_CD; 
    p[6] -> G_CD; p[7] -> G_CC; p[8] -> G_A; p[9] -> G_D;
    G_D <- 9*G_D; G_A <- G_D*G_A; G_CC <- G_A*G_CC; G_CD <- G_CC*G_CD; 
    p[10] -> S_CD; p[11] -> S_CC; p[12] -> S_A; p[13] -> S_D;
    S_D <- 9*S_D; S_A <- S_D*S_A; S_CC <- S_A*S_CC; S_CD <- S_CC*S_CD; 
    p[14] -> Q_CD; p[15] -> Q_CC; p[16] -> Q_A; p[17] -> Q_D; p[18] -> Q_P;
    Q_P <- 3*Q_P; Q_D <- Q_P*Q_D; Q_A <- Q_D*Q_A; Q_CC <- Q_A*Q_CC; Q_CD <- Q_CC*Q_CD; 

    ad = 0.95*M_D + 0.28*G_D + 0.43*S_D + 2.25*Q_D
    as = 0.017*M_A + 0.0064*G_A + 0.0076*S_A + 0.034*Q_A
    ccb = 0.0093*M_CC+ 0.0028*G_CC + 0.0042*S_CC + 0.0186*Q_CC
    ccd = 0.0223*M_CD + 0.0056*G_CD + 0.0078*S_CD + 0.0446*Q_CD
    apb = 1.28*M_P + 2.56*Q_P 
    r1=(1+ccb*(1+ccd))*ad*as*100/(130-apb)
    -r1
}
require(minqa)
p <- rep(.1, 18)
bobyqa(p, dumFun, lower = rep(0, length(p)), upper = rep(1, length(p)))
parameter estimates: 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 
objective: -9.65605526502482 
number of function evaluations: 97 
>

I think I understand what bobyqa does but I don't understand how you do to make constraints function of other variables. What do you mean by mapping parameters to [0,1] and why are you linking paramaters by multiplying them in dumfun function ? — Wicelo, Jun 09 '13 at 05:51
for example if `0 < x,y,z < 1` and redefine `y1 = x*y` and z1 = y*z then `0 < z1 < y1 < x < 1`. — user1609452, Jun 09 '13 at 08:38
Sorry I'm not sure to understand, the goal is to get the value of paramaters which must be integers from 0 to 9. — Wicelo, Jun 12 '13 at 18:34
The fact that they are integers is rather important and something you failed to mention. — user1609452, Jun 12 '13 at 18:36
sorry I mentionned it on my previous posts. On this thread too but a moredator suppressed this information. — Wicelo, Jun 12 '13 at 18:42

score 0 · Accepted Answer · answered Jun 12 '13 at 18:36

I finally solved my problem not with vectorization but with C. My program containing 14nested loops executes 100 to 1000 times faster with C than with R ! Which is sad because I didn't learn anything new from that and it prooves that R can be pretty useless on some problems but what can we do.

constrained optimization of a complicated function

2 Answers2