Блог

У меня есть система нелинейных уравнений, которую я хочу решить, но мне трудно написать функцию для оптимизации nleqslv .

Вот что я хочу сделать математически; я хочу свести к минимуму:

Таким образом, значения a являются константами, и я хочу найти значения x, которые минимизируют эту сумму.

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

Таким образом, совокупная сумма xs всегда должна быть больше или равна 0 и меньше или равна scm, если klt;N

Однако совокупная сумма всех xs вместе должна быть равна 0:

Наконец, каждый из xi может быть отрицательным (кроме первого) и связан минимальным и максимальным значением, и эти значения являются функцией совокупной суммы всех xs до него:

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

 a lt;- c(20, 34, 22, 27) scm lt;- 9300 finj lt;- function(x, inj_max){  (1/(10 * x 1)) * inj_max }  fwit lt;- function(x, wit_max){  log( x  1) * wit_max } INJ lt;- 4650 WIT lt;- 4650  

Эти функции преобразуют последнее ограничение в:

 fn lt;- function(x){    #(0 lt;= x1) can't be expressed - so I put it x1   0.0000001  c(  x[1] - scm,  x[1]   0.0000001,  x[1]   x[2] - scm,  x[1]   x[2]   0.0000001,  x[1]   x[2]   x[3] - scm,  x[1]   x[2]   x[3]   0.0000001,  x[1]   x[2]   x[3]   x[4],  #now start inj/wit constraints  x[2] * (10 * x[1]  1) - INJ,  x[2]   log(x[1]  1) * WIT   0.0000001,  x[3] * (10 * (x[1]   x[2])  1) - INJ,  x[3]   log(x[1]   x[2]  1) * WIT   0.0000001,  x[4] * (10 * (x[1]   x[2]   x[3])  1) - INJ,  x[4]   log(x[1]   x[2]   x[3]   1) * WIT  0.0000001  )   }  nleqslv(c(4650, -4650, 4650, -4650), fn)  

Я тоже написал это function и попытался решить эту проблему, но получил ошибку:

 Error in nleqslv(c(4650, -4650, 4650, -4650), fn) :   Length of fn result lt;gt; length of x!  

Логично, что я получаю эту ошибку, так как у меня так много ограничений, поэтому я не знаю, как я могу решить эту проблему оптимизации или как я могу переписать ограничения, чтобы не было этой ошибки.

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

Что касается ограничений, то каждый параметр ограничен минимумом и максимумом, но границы зависят последовательно, что делает его немного сложнее. Один из подходов заключается в последовательном преобразовании входных данных в их допустимое пространство. Это позволяет функции принимать любые реальные значения в качестве входных данных, поскольку она автоматически преобразует их в соответствии с ограничениями. Я использовал pnorm для трансформации.

Другая вещь, которую следует учитывать, заключается в том, что проблема имеет N - 1 степени свободы, так sum(x) как должна быть равна 0. Способ справиться с этим состоит в том, чтобы передать N - 1 только параметры функции, которая должна быть оптимизирована, а затем настроена x[N] на -sum(x[-N]) выполнение .

Вот несколько примеров кода, использующего ваши «поддельные значения».:

 scm lt;- 9300 INJ lt;- 4650 WIT lt;- 4650  a lt;- c(20, 34, 22, 27)  fT lt;- function(xT) {  # transforms the input values xT into values that meet the problem constraints  x lt;- numeric(length(xT)   1)  mini lt;- 0 # the minimum for parameter 1  maxi lt;- min(scm, INJ) # the maximum for parameter 1  x[1] lt;- (maxi - mini)*(pnorm(xT[1]))   mini # transform xT[1] to a value between mini and maxi  xcumsum lt;- x[1]    for (i in 2:length(xT)) {  mini lt;- max(-xcumsum, -WIT*log(xcumsum   1)) # calculate the minimum for parameter i  maxi lt;- min(scm - xcumsum, INJ/(10*xcumsum   1)) # calculate the maximum for parameter i  x[i] lt;- (maxi - mini)*(pnorm(xT[i]))   mini # transform xT[i] to a value between mini and maxi  xcumsum lt;- xcumsum   x[i]  }    x[i   1] lt;- -xcumsum  return(x) }  fn lt;- function(xT) {  return(sum(a*fT(xT))) }  # optimize fn using a vector of N - 1 zeros as the initial guess gt; optim(numeric(length(a) - 1), fn) $par [1] 17.3 -23.2 9.2  $value [1] -88350  $counts function gradient   32 NA   $convergence [1] 0  $message NULL  

Возвращаемые значения optim$par являются преобразованными значениями. Инвертировать преобразование с помощью fT :

 x lt;- fT(optim(numeric(length(a) - 1), fn)$par) gt; x [1] 4650 -4650 4650 -4650  

Переход x[1:3] к fn дает минимизированное значение функции:

 gt; fn(head(x, -1)) [1] -88350  

Что подтверждается вычислением fn вручную:

 gt; sum(a*x) [1] -88350  

ОБНОВЛЕНИЕ 1: Из комментариев о сходимости к неоптимальному локальному минимуму я попробовал различные доступные методы optim , которые применимы здесь, и метод «L-BFGS-B» действительно находит глобальный для этого случая, но трудно сказать, будет ли он в целом сходиться к глобальному минимуму:

 cm lt;- 4650 INJ lt;- 4650 WIT lt;- 4650 a lt;- c(20, 19, 22, 27)  gt; optim(numeric(length(a) - 1), fn)$value [1] -32550 gt; optim(numeric(length(a) - 1), fn, method = "BFGS")$value [1] -32550 gt; optim(numeric(length(a) - 1), fn, method = "CG")$value [1] -32550 gt; optim(numeric(length(a) - 1), fn, method = "L-BFGS-B")$value [1] -37200 gt; optim(numeric(length(a) - 1), fn, method = "SANN")$value [1] -32550.1  

ОБНОВЛЕНИЕ 2: В ответ на вопрос о том, как приведенный выше код обрабатывает N - 1 параметры, я укажу на несколько вещей:

1. N - 1 вектор длины передается optim (см. numeric(length(a) - 1) )
2. fT принимает вектор ( xT ) и выводит вектор длины length(xT) 1 (см. x lt;- numeric(length(xT) 1) и x[i 1] lt;- -xcumsum )
3. Я никогда не создавал объект N , но i = N - 1 как только for цикл завершится, это x[i 1] lt;- -xcumsum повлияет x на то же, N lt;- length(a); x[N] lt;- -sum(x[-N]) что и с тех пор xcumsum , — это запаздывающая совокупная сумма