Ослабление билинейного члена в целевой функции - с использованием конвертов CVXR и Маккормика

У меня есть целевая функция

$целевая функция$

с некоторыми ограничениями

$ограничения$

который я хочу минимизировать.

Я хочу использовать пакет R CVXR и конверты McCormick.

Давайте проверим код:

library(CVXR)  # if necessary
x <- Variable(1)
y <- Variable(1)
w <- Variable(1)
objective <- Minimize(5*x^2 + 14*w + 10*y^2 -76*x -108*y +292)
constraints <- list(x >= 0, x <= 100,
                    y >= 0, y <= 100,
                    x+2*y <= 10, 
                    x+y <= 6,
                    w >= 0, w >= 100*x + 100*y - 10000,  # constraints according to McCormick envelopes
                    w <= 100*y, w <= 100*x)  # constraints according to McCormick envelopes
prob_OF <- Problem(objective, constraints)
solution_OF <- solve(prob_OF)
solution_OF$value
## -125.0667
solution_OF$getValue(x)                  
## 2.933333
solution_OF$getValue(y)
## 3.066667
solution_OF$getValue(w)
## 1.000135e-30

Но, к сожалению, есть один недостаток:

значения для $x-value$ и $y-value$ являются приблизительными значениями, потому что McCormick - это ослабление ограничений.
These are the true values:
- $true_xy$

Чтобы преодолеть этот недостаток, вы можете использовать этот подход (см. ответ josliber), где вы разбиваете интервал на одну из переменных. Как и josliber, я разделился на переменную x...

В следующей части возникают мои проблемы, но сначала взгляните на мой код:

library(CVXR)  # if necessary

# the two variables for the objective function
x <- Variable(1)
y <- Variable(1)

# five variables for w - substitute of the bilinear part 
w1 <- Variable(1)
w2 <- Variable(1)
w3 <- Variable(1)
w4 <- Variable(1)
w5 <- Variable(1)

# five boolean variables to know which range of the x-axis is used
z1 <- Variable(1, boolean=TRUE)
z2 <- Variable(1, boolean=TRUE)
z3 <- Variable(1, boolean=TRUE)
z4 <- Variable(1, boolean=TRUE)
z5 <- Variable(1, boolean=TRUE)

# objective function with the five linear substitutes for xy
objective <- Minimize(5*x^2 + 14*w1 + 14*w2 + 14*w3 + 14*w4 + 14*w5 + 10*y^2 -76*x -108*y +292)

# for convenience - to build up the constraints easier
x_Range=matrix(rep(seq(from=0,to=100,by=20),each=2)[-c(1,12)],nrow=2,byrow=FALSE)
##     [,1] [,2] [,3] [,4] [,5]
##[1,]    0   20   40   60   80
##[2,]   20   40   60   80  100
y_Range=matrix(c(0,100),nrow=2,byrow=FALSE)
##     [,1]
##[1,]    0
##[2,]  100

# FIRST alternative of the constraints
constraints <- list(x >= 0, x <= 100,
                    y >= 0, y <= 100,
                    x+2*y <= 10,  
                    x+y <= 6,
                    ## w1
                    w1 <= x_Range[2,1]*y + x1*y_Range[1,1] - x_Range[2,1]*y_Range[1,1]*z1,  #w1 <= xU1*y + x1*yL - xU1*yL*z1,
                    w1 <= x1*y_Range[2,1] + x_Range[1,1]*y - x_Range[1,1]*y_Range[2,1]*z1,  #w1 <= x1*yU + xL1*y - xL1*yU*z1,
                    x1 >= x_Range[1,1]*z1, x1 <= x_Range[2,1]*z1,                           #xL1*z1 <= x1 <= xU1*z1,
                    ## w2
                    w2 <= x_Range[2,2]*y + x2*y_Range[1,1] - x_Range[2,2]*y_Range[1,1]*z2,  #w2 <= xU2*y + x2*yL - xU2*yL*z2,
                    w2 <= x2*y_Range[2,1] + x_Range[1,2]*y - x_Range[1,2]*y_Range[2,1]*z2,  #w2 <= x2*yU + xL2*y - xL2*yU*z2,
                    x2 >= x_Range[1,2]*z2, x2 <= x_Range[2,2]*z2,                           #xL2*z2 <= x2 <= xU2*z2,
                    ## w3
                    w3 <= x_Range[2,3]*y + x3*y_Range[1,1] - x_Range[2,3]*y_Range[1,1]*z3,  #w3 <= xU3*y + x3*yL - xU3*yL*z3,
                    w3 <= x3*y_Range[2,1] + x_Range[1,3]*y - x_Range[1,3]*y_Range[2,1]*z3,  #w3 <= x3*yU + xL3*y - xL3*yU*z3,
                    x3 >= x_Range[1,3]*z3, x3 <= x_Range[2,3]*z3,                           #xL3*z3 <= x3 <= xU3*z3,
                    ## w4
                    w4 <= x_Range[2,4]*y + x4*y_Range[1,1] - x_Range[2,4]*y_Range[1,1]*z4,  #w4 <= xU4*y + x4*yL - xU4*yL*z4,
                    w4 <= x4*y_Range[2,1] + x_Range[1,4]*y - x_Range[1,4]*y_Range[2,1]*z4,  #w4 <= x4*yU + xL4*y - xL4*yU*z4,
                    x4 >= x_Range[1,4]*z4, x4 <= x_Range[2,4]*z4,                           #xL4*z4 <= x4 <= xU4*z4,
                    ## w5
                    w5 <= x_Range[2,5]*y + x5*y_Range[1,1] - x_Range[2,5]*y_Range[1,1]*z5,  #w5 <= xU5*y + x5*yL - xU5*yL*z5,
                    w5 <= x5*y_Range[2,1] + x_Range[1,5]*y - x_Range[1,5]*y_Range[2,1]*z5,  #w5 <= x5*yU + xL5*y - xL5*yU*z5,
                    x5 >= x_Range[1,5]*z5, x5 <= x_Range[2,5]*z5                            #xL5*z5 <= x5 <= xU5*z5
                    )

# SECOND alternative of the constraints
constraints <- list(x >= 0, x <= 100,
                    y >= 0, y <= 100,
                    x+2*y <= 10,  
                    x+y <= 6, 
                    ## w1
                    w1 >= x_Range[1,1]*y + x1*y_Range[1,1] - x_Range[1,1]*y_Range[1,1]*z1,  #w1 >= xL1*y + x1*yL - xL1*yL*z1,
                    w1 >= x_Range[2,1]*y + x1*y_Range[2,1] - x_Range[2,1]*y_Range[2,1]*z1,  #w1 >= xU1*y + x1*yU - xU1*yU*z1,
                    w1 <= x_Range[2,1]*y + x1*y_Range[1,1] - x_Range[2,1]*y_Range[1,1]*z1,  #w1 <= xU1*y + x1*yL - xU1*yL*z1,
                    w1 <= x1*y_Range[2,1] + x_Range[1,1]*y - x_Range[1,1]*y_Range[2,1]*z1,  #w1 <= x1*yU + xL1*y - xL1*yU*z1,
                    x1 >= x_Range[1,1]*z1, x1 <= x_Range[2,1]*z1,                           #xL1*z1 <= x1 <= xU1*z1,
                    ## w2
                    w2 >= x_Range[1,2]*y + x2*y_Range[1,1] - x_Range[1,2]*y_Range[1,1]*z2,  #w2 >= xL2*y + x2*yL - xL2*yL*z2,
                    w2 >= x_Range[2,2]*y + x2*y_Range[2,1] - x_Range[2,2]*y_Range[2,1]*z2,  #w2 >= xU2*y + x2*yU - xU2*yU*z2,
                    w2 <= x_Range[2,2]*y + x2*y_Range[1,1] - x_Range[2,2]*y_Range[1,1]*z2,  #w2 <= xU2*y + x2*yL - xU2*yL*z2,
                    w2 <= x2*y_Range[2,1] + x_Range[1,2]*y - x_Range[1,2]*y_Range[2,1]*z2,  #w2 <= x2*yU + xL2*y - xL2*yU*z2,
                    x2 >= x_Range[1,2]*z2, x2 <= x_Range[2,2]*z2,                           #xL2*z2 <= x2 <= xU2*z2,
                    ## w3
                    w3 >= x_Range[1,3]*y + x3*y_Range[1,1] - x_Range[1,3]*y_Range[1,1]*z3,  #w3 >= xL3*y + x3*yL - xL3*yL*z3,
                    w3 >= x_Range[2,3]*y + x3*y_Range[2,1] - x_Range[2,3]*y_Range[2,1]*z3,  #w3 >= xU3*y + x3*yU - xU3*yU*z3,
                    w3 <= x_Range[2,3]*y + x3*y_Range[1,1] - x_Range[2,3]*y_Range[1,1]*z3,  #w3 <= xU3*y + x3*yL - xU3*yL*z3,
                    w3 <= x3*y_Range[2,1] + x_Range[1,3]*y - x_Range[1,3]*y_Range[2,1]*z3,  #w3 <= x3*yU + xL3*y - xL3*yU*z3,
                    x3 >= x_Range[1,3]*z3, x3 <= x_Range[2,3]*z3,                           #xL3*z3 <= x3 <= xU3*z3,
                    ## w4
                    w4 >= x_Range[1,4]*y + x4*y_Range[1,1] - x_Range[1,4]*y_Range[1,1]*z4,  #w4 >= xL4*y + x4*yL - xL4*yL*z4,
                    w4 >= x_Range[2,4]*y + x4*y_Range[2,1] - x_Range[2,4]*y_Range[2,1]*z4,  #w4 >= xU4*y + x4*yU - xU4*yU*z4,
                    w4 <= x_Range[2,4]*y + x4*y_Range[1,1] - x_Range[2,4]*y_Range[1,1]*z4,  #w4 <= xU4*y + x4*yL - xU4*yL*z4,
                    w4 <= x4*y_Range[2,1] + x_Range[1,4]*y - x_Range[1,4]*y_Range[2,1]*z4,  #w4 <= x4*yU + xL4*y - xL4*yU*z4,
                    x4 >= x_Range[1,4]*z4, x4 <= x_Range[2,4]*z4,                           #xL4*z4 <= x4 <= xU4*z4,
                    ## w5
                    w5 >= x_Range[1,5]*y + x5*y_Range[1,1] - x_Range[1,5]*y_Range[1,1]*z5,  #w5 >= xL5*y + x5*yL - xL5*yL*z5,
                    w5 >= x_Range[2,5]*y + x5*y_Range[2,1] - x_Range[2,5]*y_Range[2,1]*z5,  #w5 >= xU5*y + x5*yU - xU5*yU*z5,
                    w5 <= x_Range[2,5]*y + x5*y_Range[1,1] - x_Range[2,5]*y_Range[1,1]*z5,  #w5 <= xU5*y + x5*yL - xU5*yL*z5,
                    w5 <= x5*y_Range[2,1] + x_Range[1,5]*y - x_Range[1,5]*y_Range[2,1]*z5,  #w5 <= x5*yU + xL5*y - xL5*yU*z5,
                    x5 >= x_Range[1,5]*z5, x5 <= x_Range[2,5]*z5                            #xL5*z5 <= x5 <= xU5*z5
                    )

# Now the final calculations
prob_OF <- Problem(objective, constraints)
solution_OF <- solve(prob_OF)

# results
### of FIRST constraint alternative
solution_OF$value
## NA
solution_OF$getValue(x1)
## NA
solution_OF$getValue(x2)
## NA
solution_OF$getValue(x3)
## NA
solution_OF$getValue(x4)
## NA
solution_OF$getValue(x5)
## NA

# results
### of SECOND constraint alternative
solution_OF$value
## 15.99776
solution_OF$getValue(x1)
## 4.206985
solution_OF$getValue(x2)
## 28.49989
solution_OF$getValue(x3)
## 49.34965
solution_OF$getValue(x4)
## 69.59733
solution_OF$getValue(x5)
## 89.71508

Вопросы:

Why is josliber only using two of the four inequalities?
- Therefore also my attempt of the second constraint alternative in the code

Can you please help me to solve my issue?!
- I think I have to use second constraint alternative but here the values of solution_OF$getValue(x'ses) are so high that like josliber mentioned to sum them up won't get me the expected result of x=2.
- В обоих вариантах есть два ограничения x+2*y <= 10 и x+y <= 6. Должен ли я преобразовать их из-за того, что я разделил диапазон x на пять частей?
Можно ли также разделить более чем на одну переменную? Ссылаясь на ответ от josliber.

tueftla 17.06.2021 источник

comment

К вашему сведению: есть глобальные решатели, которые делают это автоматически. Я решаю множество невыпуклых задач и никогда не записываю эти конверты Маккормика. - Erwin Kalvelagen 17.06.2021

comment

@Erwin Kalvelagen: И не могли бы вы назвать такой решатель для конвертов Маккормика - лучше всего применимый в среде CVXR. Я новичок в этой теме, поэтому пишу. За каждое удобство я благодарен. - tueftla 18.06.2021

comment

CVXR предназначен только для выпуклых моделей (CVX означает ConVeX). - Erwin Kalvelagen 18.06.2021

comment

Не CVXR, но с использованием optim в базе R:

g <- function(z, x = z[1], y = z[2]) 5*x^2 + 14*x*y + 10*y^2 -76*x -108*y +292 + max(x+2*y-10, 0) + max(x+y-6, 0) + max(-x, 0) + max(-y, 0); ans <- optim(c(1, 1), g); ans$par

- G. Grothendieck 21.06.2021

comment

@ГРАММ. Гротендик: Большое спасибо за то, что показали формулировку в optim. Но меня действительно интересует использование конвертов Маккормика, таких как josliber пояснил в своем ответе. Причина, по которой я это делаю, заключается в том, что я потерпел неудачу в реализации... - tueftla 22.06.2021

Ответы (1)

arrow_upward
3
arrow_downward

Мы можем сделать задачу выпуклой:

min 5z^2+(1/5)y^2-76x-108*y+292
    z = x+(7/5)y
    x+2y <= 10
    x+y <= 6

Это можно ввести в CVXR и решить с помощью любого решателя QP. (Конечно, проверьте мою математику).

Erwin Kalvelagen 18.06.2021

comment

Это также отличный совет переформулировать уравнение, чтобы получить выпуклую задачу. Но я был весьма заинтересован в правильном использовании конвертов Маккормика - хотя вы прокомментировали выше, что вы никогда не будете программировать эти конверты Маккормика вручную. Здесь josliber показал процедуру, но я не смог ее реализовать. Поэтому я и задал этот вопрос... - tueftla; 21.06.2021

Ослабление билинейного члена в целевой функции - с использованием конвертов CVXR и Маккормика

Ответы (1)

Вопросы по теме