Максимизация проблемы нелинейных ограничений с использованием R-пакета nloptr

Мне нужно максимизировать целевую функцию для некоторых проблем с помощью пакета R 'nloptr'. Я пробовал базовое правило «Максимизировать f (x) ‹=> Минимизировать -f (x)», но оно не работает. Я не уверен, что не так, либо его использую, либо есть другой способ.

Вот полный пример. Текущее решение - это просто начальный вектор с минимальным целевым значением. Но я должен получить решение, которое максимизирует целевую функцию. Может кто-нибудь, пожалуйста, помогите мне, как его получить. Спасибо!

library(nloptr)
X = log(rbind(c(1.350, 8.100),   
          c(465.000,     423.000),
          c(36.330  ,    119.500),
          c(27.660   ,   115.000), 
          c(1.040     ,  5.500),
          c(11700.000, 50.000),  
          c(2547.000  ,  4603.000),
          c(187.100    , 419.000), 
          c(521.000  ,   655.000), 
          c(10.000    ,  115.000), 
          c(3.300    ,   25.600),  
          c(529.000   ,  680.000), 
          c(207.000  ,   406.000), 
          c(62.000    ,  1320.000),
          c(6654.000 ,   5712.000),
          c(9400.000  ,  70.000),  
          c(6.800      , 179.000), 
          c(35.000   ,   56.000),  
          c(0.120  ,     1.000),   
          c(0.023   ,    0.400),   
          c(2.500  ,     12.100),  
          c(55.500  ,    175.000), 
          c(100.000  ,   157.000), 
          c(52.160   ,   440.000), 
          c(87000.000 ,  154.500), 
          c(0.280  ,     1.900),   
          c(0.122  ,     3.000),   
          c(192.000 ,    180.000)))
n = nrow(X)
q = 0.5
x0 = cbind(8,4)
x01 = x0[1]
x02 = x0[2]
x1 = X[,1]
x2 = X[,2]
pInit = c(0.1614860, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 
          0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 
          0.0000000, 0.0000000, 0.0000000, 0.7124934, 0.0000000, 0.0000000,
          0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 
          0.1260206, 0.0000000, 0.0000000, 0.0000000)

eval_f0 = function(p) {
  obj0 = mean((n * p ) ^ q)
  grad0 = rbind(q * ((n * p) ^ (q - 1))/((mean((n * p ) ^ q))^2))
  return(list("objective" = obj0, "gradient" = grad0))
}

eval_g_eq0 = function(p) {
  sum0 = sum(x1 * p) - x01
  sum1 =  sum(x2 * p) - x02
  sum2 = sum(p) - 1
  constr0 = rbind(sum0, sum1, sum2)
  grad0 = rbind(x1, x2, rep(1,n))
  return(list("constraints" = constr0, "jacobian" = grad0))
}

local_opts <- list( "algorithm" = "NLOPT_LD_AUGLAG",
                    "xtol_rel" = 1.0e-7 )
opts <- list( "algorithm" = "NLOPT_LD_AUGLAG",
              "xtol_rel" = 1.0e-7,
              "maxeval" = 10000,
              "local_opts" = local_opts )

  res1 = nloptr(x0 = c(pInit), 
                eval_f = eval_f0, 
                lb = c(rep(0, n)), 
                ub = c(rep(Inf, n)), 
                eval_g_eq = eval_g_eq0, 
                opts = opts )

  weight = res1$solution
  fval0 = res1$objective

print(list(fval0, weight))

r nonlinear-optimization
person Jafer    schedule 01.11.2018    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Обратите внимание, что градиент (и якобиан) в вашей начальной точке pInit не является конечным, что усложняет эту задачу для любого решателя на основе градиента. Я буду использовать другую начальную точку, немного дальше от границы.

В любом случае, кажется, легче найти максимум с помощью лагранжевого решателя из пакета alabama. С вашими определениями выше до x1 = X[,1]; x2 = X[,2] возможное решение выглядит так:

f1 <- function(p) mean((n * p ) ^ q)
heq1 <- function(p) 
    c(sum(x1 * p) - x01, sum(x2 * p) - x02, sum(p) - 1)

Для простоты мы позволяем программе расчета градиентов и якобианов. Чтобы найти максимум, примените решатель к отрицательному значению целевой функции.

sol <- alabama::auglag(rep(0.1, 28), fn=function(p) -f1(p), heq=heq1)
cat("The maximum value is:", -sol$value, '\n')
## The maximum value is: 0.7085338

Условия равенства выполнены, см.

heq1(sol$par)
## [1] -1.685957e-08  3.721533e-08 -2.935964e-08

и найденное решение

sol$par
##  [1] 0.012186842 0.006640286 0.006706268 0.006418224 0.014501609 0.405618998  
##  [7] 0.003531462 0.005458189 0.005582029 0.005158098 0.008072278 0.005510394  
## [13] 0.005653117 0.002935642 0.003861549 0.123009564 0.004021866 0.009866779  
## [19] 0.024385229 0.027101557 0.011436010 0.006184886 0.007473135 0.004162962  
## [25] 0.245429952 0.019978294 0.010919515 0.008195238

Мне было бы интересно узнать, является ли это разумным решением для вас! Я проверил его по нескольким отправным точкам, и всегда получалось одинаково.

person Hans W.    schedule 01.11.2018
comment
Спасибо за подробный ответ. Указание на проблему очень заметно. Однако установка начальной точки вручную в моем случае не применима. Я привел простой пример, в котором начальная точка была создана на основе теоретической идеи, лежащей в основе этого. Фактически, я должен использовать его в исследовании моделирования, где он будет рассчитываться программой при каждой репликации. - person Jafer; 02.11.2018
comment
Вы должны избегать нулевых компонентов в вашем начальном векторе в любом случае при применении решателя оптимизации на основе градиента - нет никакого способа обойти это. Хорошая новость: добавление небольшой ненулевой суммы к pInit, такой как pInit+0.001, приведет к тому же максимуму. и это можно сделать программно. (Или вы добавляете 0.001 только к тем компонентам, которые равны нулю.) - person Hans W.; 02.11.2018
comment
Есть еще одна проблема: когда мы используем эту процедуру для моделирования, она дает ошибку Ошибка в собственном ($ hessian, симметричный = ИСТИНА, only.values ​​= ИСТИНА): бесконечные или отсутствующие значения в 'x'. Я не знаю, как указать на это и преодолеть это. - person Jafer; 03.11.2018
comment
Вы могли намекнуть, что открыли новый вопрос Ошибка в задаче нелинейной оптимизации: бесконечные или отсутствующие значения для обращения к сообщению об ошибке. - person Hans W.; 06.11.2018