Нелинейная оптимизация / программирование с целочисленными переменными в R - PullRequest
5 голосов
/ 05 апреля 2020

Интересно, кто-нибудь может предложить несколько пакетов для решения задачи нелинейной оптимизации, которые могут предоставить целочисленные переменные для оптимального решения? Проблема состоит в том, чтобы минимизировать функцию с ограничением равенства, подчиняющимся некоторым ограничениям нижней и верхней границы.

Я использовал пакет 'nloptr' в R для задачи нелинейной оптимизации, которая хорошо работала, но теперь хотела бы расширить метод, чтобы некоторые переменные были целыми числами. Из моего использования и понимания nloptr до сих пор он может возвращать только непрерывные, а не целочисленные переменные для оптимального решения.

Я полагаю, что эту проблему необходимо решить с помощью смешанного целочисленного нелинейного программирования.

Один из примеров проблемы в форме для nloptr:

min f(x) (x-y)^2/y + (p-q)^2/q
so that (x-y)^2/y + (p-q)^2/q = 10.2

where 
x and p are positive integers not equal to 0 
and 
y and q may or may not be positive integers not equal to 0

Код nloptr для этого в R будет выглядеть так:

library('nloptr')

x1 <- c(50,25,20,15)

fn <- function(x) {
  (((x[1] - x[2])^2)/x[2]) + (((x[3] - x[4])^2)/x[4])
  }

heq <- function(x) {
  fn(x)-10.2
}

lower_limit <- c(0,0,0,0)
upper_limit <- c(67.314, 78, 76.11, 86)


slsqp(x1, fn, lower = lower_limit, upper = upper_limit,  hin = NULL, heq = heq, control = list(xtol_rel = 1e-8, check_derivatives = FALSE))

Это выведет следующее:

$par
[1] 46.74823 29.72770 18.93794 16.22137

$value
[1] 10.2

$iter
[1] 6

$convergence
[1] 4

$message
[1] "NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached."

Это результат, который я ищу, но как указано выше, мне нужны x и p как целые числа.

Я посмотрел на https://cran.r-project.org/web/views/Optimization.html, в котором есть действительно хороший список пакетов для смешанного целочисленного нелинейного программирования, но просто Интересно, имел ли кто-нибудь опыт работы с кем-либо из них и что, по их мнению, может быть наиболее подходящим для решения проблемы, как указано выше.

Был похожий вопрос об этом, опубликованный около 7 лет назад go здесь, но он закончился ссылкой на страницу крана, так что подумал, что стоит спросить снова.

Спасибо очень за ваш вклад.

ура,

Андрей

Ответы [ 5 ]

4 голосов
/ 07 апреля 2020

Вот пример того, как не не работает с CVXR, без более простой целевой функции. Я подозреваю, что проблема не выпуклая, даже с ограничениями, поэтому требуется альтернативный вариант.

#base example from https://cvxr.rbind.io/cvxr_examples/cvxr_gentle-intro/
#install.packages("CVXR")
library(CVXR)

#modified for Stackoverflow integer MIQP ####
#Solves, but terms not normalised by y and q respectively

# Variables minimized over
x <- Variable(1, integer=TRUE)
y <- Variable(1)
p <- Variable(1, integer=TRUE)
q <- Variable(1)

# Problem definition (terms not normalised by y and q respectively)
objective <- Minimize((x - y)^2 + (p - q)^2)
constraints <- list(x >= 0, y >= 0, p >= 0, q >= 0, 
                    x <= 67.314, y <= 78, p <= 76.11, q <= 86)
prob2.1 <- Problem(objective, constraints)

# Problem solution
solution2.1 <- solve(prob2.1)
solution2.1$status
solution2.1$value
solution2.1$getValue(x)
solution2.1$getValue(y)
solution2.1$getValue(p)
solution2.1$getValue(q)


#modified for Stackoverflow integer NLP (not integer) ####
#Does not solve, not convex?

# Variables minimized over
x <- Variable(1)
y <- Variable(1)
p <- Variable(1)
q <- Variable(1)

# Problem definition
objective <- Minimize((x - y)^2/y + (p - q)^2/q)
constraints <- list(x >= 0, y >= 0, p >= 0, q >= 0, 
                    x <= 67.314, y <= 78, p <= 76.11, q <= 86)
prob2.1 <- Problem(objective, constraints)

# Problem solution
solution2.1 <- solve(prob2.1, gp = TRUE)
solution2.1 <- solve(prob2.1, gp = FALSE)
# solution2.1$status
# solution2.1$value
# solution2.1$getValue(x)
# solution2.1$getValue(y)
# solution2.1$getValue(p)
# solution2.1$getValue(q)


#modified for Stackoverflow integer MINLP ####
#Does not solve

# Variables minimized over
x <- Variable(1, integer=TRUE)
y <- Variable(1)
p <- Variable(1, integer=TRUE)
q <- Variable(1)

# Problem definition
objective <- Minimize((x - y)^2/y + (p - q)^2/q)
constraints <- list(x >= 0, y >= 0, p >= 0, q >= 0, 
                    x <= 67.314, y <= 78, p <= 76.11, q <= 86)
prob2.1 <- Problem(objective, constraints)

# Problem solution
solution2.1 <- solve(prob2.1, gp = TRUE)
solution2.1 <- solve(prob2.1, gp = FALSE)
# solution2.1$status
# solution2.1$value
# solution2.1$getValue(x)
# solution2.1$getValue(y)
# solution2.1$getValue(p)
# solution2.1$getValue(q)


#modified for Stackoverflow integer NLP (not integer) ####
#objective multiplied by y*q, Does not solve, not convex?

# Variables minimized over
x <- Variable(1)
y <- Variable(1)
p <- Variable(1)
q <- Variable(1)

# Problem definition
objective <- Minimize((x - y)^2*q + (p - q)^2*y)
constraints <- list(x >= 0, y >= 0, p >= 0, q >= 0, 
                    x <= 67.314, y <= 78, p <= 76.11, q <= 86)
prob2.1 <- Problem(objective, constraints)

# Problem solution
solution2.1 <- solve(prob2.1, gp = TRUE)
solution2.1 <- solve(prob2.1, gp = FALSE)
# solution2.1$status
# solution2.1$value
# solution2.1$getValue(x)
# solution2.1$getValue(y)
# solution2.1$getValue(p)
# solution2.1$getValue(q)
2 голосов
/ 08 апреля 2020

Поскольку эта проблема имеет тип, который трудно решить, любой общий алгоритм не гарантированно подходит для этой конкретной задачи (см. Теорему о бесплатном обеде). Действительно, многие алгоритмы даже не могут сходиться к глобальному оптимуму для сложной задачи. Интересно, что случайный поиск проблемного пространства, по крайней мере, в конечном итоге сойдет, потому что в конечном итоге он ищет все пространство!

tl / dr Попробуйте перечислить проблемное пространство. Например, если ваши четыре переменные являются целыми числами от 0 до 80, есть только ~ 80 ^ 4 = ~ 40 миллионов комбинаций, через которые вы могли бы пройти l oop. Промежуточным вариантом может быть (если только две переменные являются целыми числами) для решения проблемы с помощью методов оптимизации для двух оставшихся переменных с учетом значения для двух целых чисел (может быть, это теперь выпуклая задача?) И l oop через для целочисленные значения.

2 голосов
/ 08 апреля 2020

ROI - вариант для проблем с MINLP. Я считаю, что у него есть доступ к некоторым подходящим решателям. Это позволяет получить доступ к neos (описано в другом ответе на ваш вопрос).

Если вам интересно посмотреть, как выглядит оптимизация ROI, вот LP (пример линейного программирования:

#ROI example https://epub.wu.ac.at/5858/1/ROI_StatReport.pdf
#install.packages("ROI")
library(ROI)
ROI_available_solvers()

ROI_registered_solvers() #has one solver loaded by default

## Get and load "lpsolve" solver
#install.packages("ROI.plugin.lpsolve", repos=c("https://r-forge.r-project.org/src/contrib",
#                                   "http://cran.at.r-project.org"),dependencies=TRUE)
library(ROI.plugin.lpsolve)
ROI_registered_solvers() #Now it is available to use

## Describe model
A <- rbind(c(5, 7, 2), c(3, 2, -9), c(1, 3, 1))
dir <- c("<=", "<=", "<=")
rhs <- c(61, 35, 31)
lp <- OP(objective = L_objective(c(3, 7, -12)),
         constraints = L_constraint(A, dir = dir, rhs = rhs),
         bounds = V_bound(li = 3, ui = 3, lb = -10, ub = 10, nobj = 3),
         maximum = TRUE)

## When you have a model, you can find out which solvers can solve it
ROI_available_solvers(lp)[, c("Package", "Repository")]

## Solve model
lp_sol <- ROI_solve(lp, solver = "lpsolve")
2 голосов
/ 08 апреля 2020

Я попробовал следующий код, используя ваш пример, чтобы попытаться повторить пример nloptr в исходном вопросе:

#base example from https://cvxr.rbind.io/cvxr_examples/cvxr_gentle-intro/
#install.packages("CVXR")
library(CVXR)

#modified for Stackoverflow integer MINLP (MIQP) ####
#Solves

# Variables minimized over
x <- Variable(1, integer=TRUE)
y <- Variable(1)
p <- Variable(1, integer=TRUE)
q <- Variable(1)
z <- Variable(1)

# Problem definition (terms not normalised by y and q respectively)
objective <- Minimize((x - y)^2 + (p - q)^2 -z)
constraints <- list(x <= 67.314, y <= 78, p <= 76.11, q <= 86, z == 10.2)
prob2.1 <- Problem(objective, constraints)

# Problem solution
solution2.1 <- solve(prob2.1)
solution2.1$status
solution2.1$value
solution2.1$getValue(x)
solution2.1$getValue(y)
solution2.1$getValue(p)
solution2.1$getValue(q)
solution2.1$getValue(z)

Однако я получаю это значение -10.19989, когда оно должно быть 0.

> solution2.1$status
[1] "optimal"
> solution2.1$value
[1] -10.19989
> solution2.1$getValue(x)
[1] -1060371
> solution2.1$getValue(y)
[1] -1060371
> solution2.1$getValue(p)
[1] -1517
> solution2.1$getValue(q)
[1] -1517.002
> solution2.1$getValue(z)
[1] 10.2

Я не могу понять, что мне нужно сделать для того, чтобы вышеперечисленное работало, как в примере с nloptr, но при условии, что x и p являются целыми значениями!

Cheers, Andrew

0 голосов
/ 08 апреля 2020

rneos - это пакет, который позволяет получить доступ к neos , бесплатной службе решения с многочисленными алгоритмами, включая некоторые подходящие для задач MINLP (например, BONMIN и Couenne, см. Список здесь ). К сожалению, проблему необходимо отформатировать как модель GAMS или AMPL. Для вас это может означать изучение некоторых базовых c GAMS, и в этом случае, возможно, вы могли бы просто использовать программное обеспечение GAMS , см. Здесь ? Бесплатная версия может быть достаточно для ваших целей. Он может быть запущен как командная строка, поэтому вы можете вызывать его из R. 1011 *

#rneos example
#from p11 of https://www.pfaffikus.de/talks/rif/files/rif2011.pdf

#install.packages("rneos")
library(rneos)
#library(devtools)
#install_github("duncantl/XMLRPC")
library(XMLRPC)
## NEOS: ping
Nping()
## NEOS: listCategories
NlistCategories()
## NEOS: listSolversInCategory
NlistSolversInCategory(category = "lp")
## NEOS: getSolverTemplate
template <- NgetSolverTemplate(category = "lp", solvername = "MOSEK", inputMethod = "GAMS")
template
#gams file below sourced from https://github.com/cran/rneos/blob/master/inst/ExGAMS/TwoStageStochastic.gms
modc <- paste(paste(readLines("TwoStageStochastic.gms"), collapse = "\n"), "\n")
cat(modc)
argslist <- list(model = modc, options = "", wantlog = "", comments = "")
xmls <- CreateXmlString(neosxml = template, cdatalist = argslist)
## NEOS: printQueue
NprintQueue()
## NEOS: submitJob
(test <- NsubmitJob(xmlstring = xmls, user = "rneos", interface = "", id = 0))
## NEOS: getJobStatus
NgetJobStatus(obj = test)
## NEOS: getFinalResults
NgetFinalResults(obj = test)
...