撰寫一個函式來優化非線性函式nleqslv

我有一個我想要求解的非線性方程組，但我發現很難撰寫函式來優化nleqslv.

這是我想在數學上做的事情；我想最小化：

所以 a 值是常數，我想搜索最小化這個總和的 x 值。

問題是對我所擁有的價值的大量限制。在數學上，按順序是：

因此，如果 k<N，xs的累積和必須始終大于或等于 0 且小于或等于 scm

但是，所有 x 的累積總和必須等于 0：

最后，每個 xi 都可以是負數（第一個除外）并且受最小值和最大值的約束，這些值是它之前所有 x 的累積總和的函式：

我在 R 中設定了假值，以解決這個問題的簡單版本：

a <- c(20, 34, 22, 27)
scm <- 9300
finj <- function(x, inj_max){
  (1/(10 * x 1)) * inj_max
}

fwit <- function(x, wit_max){
  log( x  1) * wit_max
}
INJ <- 4650
WIT <- 4650

這些函式將最后一個約束轉換為：

fn <- function(x){
  
  #(0 <= 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 <> length of x!

我收到這個錯誤是合乎邏輯的，因為我有很多約束，所以我不知道如何解決這個優化問題，或者我如何重寫約束來避免這個錯誤。

uj5u.com熱心網友回復：

我不認為你想要nleqslv，這是針對非線性方程組的。您正在嘗試最小化具有多個引數的單個函式。optim從基礎 R 應該作業。

至于約束，每個引數都以最小值和最大值為界，但界限是順序相關的，這使得它有點棘手。一種方法是將輸入順序轉換為允許的空間。這允許函式接受任何實際值作為輸入，因為它會自動轉換它們以滿足約束。我用于pnorm轉換。

另一件要考慮的事情是該問題具有N - 1自由度，因為 sum(x)必須為 0。處理該問題的方法是僅將N - 1引數傳遞給要優化的函式，然后設定x[N]為-sum(x[-N])。

下面是一些使用“假值”的示例代碼：

scm <- 9300
INJ <- 4650
WIT <- 4650

a <- c(20, 34, 22, 27)

fT <- function(xT) {
  # transforms the input values xT into values that meet the problem constraints
  x <- numeric(length(xT)   1)
  mini <- 0             # the minimum for parameter 1
  maxi <- min(scm, INJ) # the maximum for parameter 1
  x[1] <- (maxi - mini)*(pnorm(xT[1]))   mini # transform xT[1] to a value between mini and maxi
  xcumsum <- x[1]
  
  for (i in 2:length(xT)) {
    mini <- max(-xcumsum, -WIT*log(xcumsum   1))     # calculate the minimum for parameter i
    maxi <- min(scm - xcumsum, INJ/(10*xcumsum   1)) # calculate the maximum for parameter i
    x[i] <- (maxi - mini)*(pnorm(xT[i]))   mini      # transform xT[i] to a value between mini and maxi
    xcumsum <- xcumsum   x[i]
  }
  
  x[i   1] <- -xcumsum
  return(x)
}

fn <- function(xT) {
  return(sum(a*fT(xT)))
}

# optimize fn using a vector of N - 1 zeros as the initial guess
> 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 <- fT(optim(numeric(length(a) - 1), fn)$par)
> x
[1]  4650 -4650  4650 -4650

傳遞x[1:3]給fn給出了最小函式值：

> fn(head(x, -1))
[1] -88350

通過fn手動計算來檢查：

> sum(a*x)
[1] -88350

UPDATE1：從關于收斂到次優區域最小值的評論中，我嘗試了optim適用于此處的各種可用方法，并且“L-BFGS-B”方法確實找到了這種情況的全域，但很難說如果它通常會收斂到全域最小值：

cm <- 4650
INJ <- 4650
WIT <- 4650
a <- c(20, 19, 22, 27)

> optim(numeric(length(a) - 1), fn)$value
[1] -32550
> optim(numeric(length(a) - 1), fn, method = "BFGS")$value
[1] -32550
> optim(numeric(length(a) - 1), fn, method = "CG")$value
[1] -32550
> optim(numeric(length(a) - 1), fn, method = "L-BFGS-B")$value
[1] -37200
> optim(numeric(length(a) - 1), fn, method = "SANN")$value
[1] -32550.1

UPDATE2：針對上述代碼如何處理N - 1引數的問題，我會指出幾點：

1. 一個N - 1長度向量被傳遞給optim（見numeric(length(a) - 1)）
2. fT accepts a vector (xT) and outputs a vector of length length(xT) 1 (see x <- numeric(length(xT) 1) and x[i 1] <- -xcumsum)
3. I never created the object N, but i = N - 1 once the for loop completes, so x[i 1] <- -xcumsum affects x the same as N <- length(a); x[N] <- -sum(x[-N]) since xcumsum is a lagging cumulative sum

標籤：r 优化 非线性优化

上一篇：根據R中的日期和小時以15分鐘的間隔聚合資料

下一篇：我的plot()代碼是否有問題使其無法編織？