Практические руководства
Как решить задачу оптимизации ограничений в R
Оптимизация нелинейных ограничений в R с использованием nloptr
Вступление
Часто в исследованиях в области физических наук мы сталкиваемся с трудной проблемой оптимизации функции (называемой объективной), которая должна удовлетворять ряду ограничений - линейным или нелинейным равенствам и неравенствам. Оптимизаторы обычно также должны придерживаться верхней и нижней границы. Недавно я работал над аналогичной проблемой в квантовой информатике (QIS), где я попытался оптимизировать нелинейную функцию на основе нескольких ограничений, продиктованных правилами физики и математики. Заинтересованные читатели могут найти мою работу Оптимизация созвездия для когерентных состояний с фазовой манипуляцией с приемником смещения для максимизации взаимной информации, в которой я оптимизировал взаимную информацию для четырехфазной манипуляции (QPSK) на основе набора ограничений.
Преследуя проблему нелинейной оптимизации с набором ограничений, я обнаружил, что не все процедуры оптимизации созданы одинаково. На разных языках доступно несколько библиотек, таких как python (scipy.optimize), Matlab (fmincon), C ++ (robotim, nlopt), и R (nloptr). Хотя мой список процедур оптимизации не является исчерпывающим, некоторые из них более надежны, чем другие, некоторые обеспечивают более быстрое выполнение, чем другие, а некоторые имеют лучшую документацию. Основываясь на нескольких ключевых факторах, я считаю, что nloptr, реализованный на языке R, наиболее подходит для нелинейной оптимизации. nloptr использует nlopt, реализованный на C ++, в качестве серверной части. В результате он обеспечивает элегантность языка R и скорость C ++. Процедура оптимизации выполняется быстро за доли секунды даже с допуском порядка 10e-15.
Задача нелинейной оптимизации
Общая задача нелинейной оптимизации обычно имеет вид
где f - целевая функция, g определяет набор ограничений неравенства, h - набор ограничений равенства. xL и xU - нижняя и верхняя границы соответственно. В литературе представлено несколько алгоритмов оптимизации. Например, MMA (Метод движущихся асимптот) ¹ поддерживает произвольные ограничения нелинейного неравенства, (COBYLA) Оптимизация с ограничениями посредством линейного приближения², (ORIG_DRIECT) алгоритм DIRECT³. Алгоритмы оптимизации, которые также поддерживают ограничения нелинейного равенства, включают ISRES (улучшенная стратегия эволюции стохастического ранжирования) ⁴, (AUGLAG) Augmented Lagrangian Algorithm⁵. Полный список таких методов можно найти на справочной странице nlopt C ++ по адресу https://nlopt.readthedocs.io/en/latest/NLopt_Reference/. nloptr использует один из этих методов для решения данной задачи оптимизации.
«MMA (Метод движущихся асимптот) поддерживает произвольные ограничения нелинейного неравенства, (COBYLA) Оптимизация с ограничениями посредством линейного приближения, (ORIG_DRIECT) алгоритм DIRECT. Алгоритмы оптимизации, которые также поддерживают ограничения нелинейного равенства, включают ISRES (улучшенную стратегию эволюции стохастического ранжирования), (AUGLAG) расширенный лагранжев алгоритм ».
В оставшейся части статьи я приведу несколько примеров решения проблемы оптимизации ограничений с помощью R. Лично я использую R Studio, объединяющую компилятор и редактор R. R Studio также предоставляет инструмент knitr, который отлично подходит для написания документации или статей со встроенным кодом, который также может генерировать исходный латексный код и файл pdf. Большая часть представленного здесь примера была изменена из наборов тестов, используемых для проверки функций в пакете nloptr R.
Установка и загрузка библиотеки
Установка nloptr в R довольно проста.
install.packages(“nloptr”) library(‘nloptr’)
Пример 1: Оптимизация с явным градиентом
В первом примере мы минимизируем функцию Rosenbrock Banana
чей градиент задается
Однако не все алгоритмы в nlopt требуют явного градиента, как мы увидим в следующих примерах. Давайте сначала определим целевую функцию и ее градиент:
eval_f <- function(x)
{
return ( 100 * (x[2] - x[1] * x[1])^2 + (1 - x[1])^2 )
}
eval_grad_f <- function(x) {
return( c( -400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]),
200 * (x[2] - x[1] * x[1]) ) )
}
Также нам потребуются начальные значения оптимизаторов:
x0 <- c( -1.2, 1 )
Перед тем, как запустить процедуру минимизации, нам нужно указать, какой алгоритм мы будем использовать. Это можно сделать следующим образом:
opts <- list("algorithm"="NLOPT_LD_LBFGS",
"xtol_rel"=1.0e-8)
Здесь мы будем использовать алгоритм L-BFGS. Теперь мы готовы запустить процедуру оптимизации.
# solve Rosenbrock Banana function res <- nloptr( x0=x0, eval_f=eval_f, eval_grad_f=eval_grad_f, opts=opts)
Мы можем увидеть результат, набрав
print(res)
Функция оптимизирована на (1,1), что является основной истиной. Проверьте себя, запустив код в R Studio.
Пример 2: Минимизация с ограничением неравенства без градиентов
Проблема минимизации:
с a1 = 2, b1 = 0, a2 = −1 и b2 = 1. Мы переупорядочиваем ограничения, чтобы они имели вид g (x) ≤ 0:
Сначала определите целевую функцию
# objective function
eval_f0 <- function( x, a, b ){
return( sqrt(x[2]) )
}
и ограничения
# constraint function
eval_g0 <- function( x, a, b ) {
return( (a*x[1] + b)^3 - x[2] )
}
Определить параметры
# define parameters a <- c(2,-1) b <- c(0, 1)
Теперь решите, используя NLOPT_LN_COBYLA без информации о градиенте
# Solve using NLOPT_LN_COBYLA without gradient information
res1 <- nloptr( x0=c(1.234,5.678),
eval_f=eval_f0,
lb = c(-Inf,0),
ub = c(Inf,Inf),
eval_g_ineq = eval_g0,
opts = list("algorithm"="NLOPT_LN_COBYLA",
"xtol_rel"=1.0e-8),
a = a,
b = b )
print( res1 )
Пример 3: Минимизация с ограничениями равенства и неравенства без градиентов
Мы хотим решить следующую задачу оптимизации ограничений
с учетом ограничений
В этом примере оптимальное решение достигается при (1.00000000, 4.74299963, 3.82114998, 1.37940829).
# Objective Function
eval_f <- function(x)
{
return (x[1]*x[4]*(x[1] +x[2] + x[3] ) + x[3] )
}
# Inequality constraints
eval_g_ineq <- function(x)
{
return (25 - x[1]*x[2]*x[3]*x[4])
}
# Equality constraints
eval_g_eq <- function(x)
{
return ( x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 - 40 )
}
# Lower and upper bounds
lb <- c(1,1,1,1)
ub <- c(5,5,5,5)
#initial values
x0 <- c(1,5,5,1)
Нам также нужно будет определить параметры для оптимизации.
# Set optimization options. local_opts <- list( "algorithm" = "NLOPT_LD_MMA", "xtol_rel" = 1.0e-15 ) opts <- list( "algorithm"= "NLOPT_GN_ISRES", "xtol_rel"= 1.0e-15, "maxeval"= 160000, "local_opts" = local_opts, "print_level" = 0 )
Мы используем NL_OPT_LD_MMA для локальной оптимизации и NL_OPT_GN_ISRES для общей оптимизации. Вы можете установить очень низкий допуск, чтобы получить наилучший результат. Количество итераций задается с помощью maxeval. Установка действительно низкого допуска или очень большого количества итераций может привести к наилучшему приближению за счет увеличения времени вычислений. Наконец, оптимизируем
res <- nloptr ( x0 = x0,
eval_f = eval_f,
lb = lb,
ub = ub,
eval_g_ineq = eval_g_ineq,
eval_g_eq = eval_g_eq,
opts = opts
)
print(res)
В моем случае результат оказался (1 4,768461 3,78758 1,384204), что довольно близко к истине.
Пример 4: Минимизация с множественными ограничениями неравенства без градиентов
Наша целевая функция в этом случае
при условии
с оценками переменных как
Для этой задачи оптимальное решение достигается на (1, 1). Теперь напишем код.
# Objective function
eval_f <- function(x)
{
return ( x[1]^2 + x[2]^2 )
}
# Inequality constraints
eval_g_ineq <- function (x) {
constr <- c(1 - x[1] - x[2],
1 - x[1]^2 - x[2]^2,
9 - 9*x[1]^2 - x[2]^2,
x[2] - x[1]^2,
x[1] - x[2]^2)
return (constr)
}
# Lower and upper bounds
lb <- c(-50, -50)
ub <- c(50, 50)
# Initial values
x0 <- c(3, 1)
Наконец, определите параметры для nloptr следующим образом:
opts <- list( "algorithm" = "NLOPT_GN_ISRES", "xtol_rel" = 1.0e-15, "maxeval"= 160000, "tol_constraints_ineq" = rep( 1.0e-10, 5 ))
а затем выполните оптимизацию
res <- nloptr(
x0 = x0,
eval_f = eval_f,
lb = lb,
ub = ub,
eval_g_ineq = eval_g_ineq,
opts = opts )
print(res)
При указанном допуске и количестве итераций мне удалось достичь оптимального решения (1, 1).
Хотя я не приводил примеры с множественными ограничениями равенства, они очень похожи на Пример 4. Однако не забудьте выбрать алгоритм оптимизации как NLOPT_GN_ISRES.
Ссылки
- К. Сванберг, Метод движущихся асимптот - новый метод структурной оптимизации, Международный журнал численных методов в инженерии, 1987, 24, 359–373.
- Пауэлл М.Дж. (1994) Достижения в области оптимизации и численного анализа, Kluwer Academic, Dordrecht, Метод прямой поисковой оптимизации, моделирующий целевые функции и функции ограничений с помощью линейной интерполяции, стр. 51–67.
- Https://people.cs.vt.edu/~ltw/lecture_notes/HPCS08tut.pdf
- Томас П. Рунарссон и Синь Яо, «Стохастическое ранжирование для эволюционной оптимизации с ограничениями», IEEE Trans. Эволюционные вычисления, т. 4 (№ 3), с. 284–294 (2000).
- Https://www.him.uni-bonn.de/fileadmin/him/Section6_HIM_v1.pdf
- Https://people.kth.se/~krille/mmagcmma.pdf
- Https://cran.r-project.org/web/packages/nloptr/vignettes/nloptr.pdf
- Http://faculty.cas.usf.edu/jkwilde/mathcamp/Constrained_Optimization.pdf
- Https://nlopt.readthedocs.io/en/latest/
Если эта статья принесет вам пользу, пожалуйста, используйте следующие цитаты, чтобы ссылаться на мою работу:
Rahul Bhadani. Nonlinear Optimization in R using nlopt. arXiv preprint arXiv:2101.02912, 2021.
or
@article{bhadani2021nonlinear, title={Nonlinear Optimization in R using nlopt}, author={Rahul Bhadani}, year={2021}, eprint={2101.02912}, archivePrefix={arXiv}, primaryClass={math.OC}, journal={arXiv preprint arXiv:2101.02912}, }
