如何使用 R 最小化具有线性等式约束的非线性目标函数？

Question

最小化
f(z) = sum_(t=2)^IJ (Z_t - Z_t-1)^2

受约束
sum_(j=1)^J (Z_(i-1)J+j+k ) = y^f_i, i = 1,....., I-1.

此优化从财政年度数据系列 (y^f_i) 中找出季度值，然后将这些季度值相加得出年度值。 I 是在系列间隔中考虑的日历年数 (2

如何使用 R 解决这个问题？

我尝试编写代码的方式如下：

library(NlcOptim)
library(readxl)
Calendarization <- read_excel("C:/Users/HP/Desktop/Calendarization.xlsx")
View(Calendarization)

y<-Calendarization$`wholesale price`

objfun = function(z){
  return(sum(z[t] - lag(z[t], k=1))^2)
}

for (t in 2:156){
  objfun
} -> objfun

p0<-0:39

Aeq<-sum(z[((i-1)*4)+j+2])

for (j in 1:4){
  for (i in 1:39){
  Aeq
}->Aeq
}

Beq<- y[i]

x=p0

solnl(x, objfun=objfun, Aeq=Aeq, Beq=Beq)

这是我的数据：

year       wholesale price
1970-1971   0.99
1971-1972   1.32
1972-1973   20.9
1973-1974   2.83
1974-1975   5.78
1975-1976   3.38
1976-1977   3.02
1977-1978   2.88
1978-1979   4.08
1979-1980   5.4
1980-1981   4.51
1981-1982   5.91
1982-1983   6.42
1983-1984   7.07
1984-1985   7.68
1985-1986   8.04
1986-1987   9.62
1987-1988   10.05
1988-1989   9.81
1989-1990   9.6
1990-1991   10.59
1991-1992   11.08
1992-1993   9.42
1993-1994   9.6
1994-1995   12.28
1995-1996   12.58
1996-1997   10.87
1997-1998   12.09
1998-1999   13.66
1999-2000   12.28
2000-2001   11.75
2001-2002   11.49
2002-2003   13.08
2003-2004   13.43
2004-2005   15.06
2005-2006   16.5
2006-2007   18.48
2007-2008   24.74
2008-2009   26.69

Answer 1

这个表述似乎有问题。只看图像中的公式，i=I-1=39-1=38 的最后一个约束求和 z 个元素 (38-1)*4 + 1 + 6、(38-1)*4 + 2 + 6、(38-1)*4 + 3 + 6 和 (38-1)*4 + 4 + 6 这是元素 155 156 157 158 但 z从 1 到 4*39，因此只有 156 个元素。此外，并非所有 z 值都参与约束。

鉴于所引用的问题，让我们将问题更改为有意义的问题，并假设我们要最小化 z 的 I*J 元素的连续差的平方和，这取决于 z 求和的前 4 个元素到 Cal[1, 2]，接下来的 4 个元素求和到 Cal[2, 2]，依此类推，直到 z 的最后 4 个元素求和到 Cal[39, 2]。在这种情况下，我们可以使用所示的克罗内克积将 Aeq 约束矩阵写为块对角矩阵。我们忽略 K。（Cal 在最后的注释中重复显示。）

library(NlcOptim)

I = 39; J = 4
objfun <- function(x) sum(diff(x)^2)
Aeq <- diag(I) %x% matrix(1, 1, J)
Beq <- Cal[, 2]
st <- rep(1, I*J)
res <- solnl(st, objfun, Aeq = Aeq, Beq = Beq)

给予

> res
List of 6
 $ par    : num [1:156, 1] 0.576 0.445 0.182 -0.213 -0.739 ...
 $ fn     : num 25.5
 $ counts : num [1, 1:2] 19332 124
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "nfval" "ngval"
 $ lambda :List of 3
  ..$ lower: num [1:156, 1] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ upper: num [1:156, 1] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ eqlin: num [1:39] 0.263 1.486 2.427 1.662 0.68 ...
 $ grad   : num [1:156, 1] 0.263 0.263 0.263 0.263 -1.486 ...
 $ hessian: num [1:156, 1:156] 1.65775 -1.10379 0.62425 -0.09085 0.00878 ...

注意

Lines <- "year       wholesale price
1970-1971   0.99
1971-1972   1.32
1972-1973   20.9
1973-1974   2.83
1974-1975   5.78
1975-1976   3.38
1976-1977   3.02
1977-1978   2.88
1978-1979   4.08
1979-1980   5.4
1980-1981   4.51
1981-1982   5.91
1982-1983   6.42
1983-1984   7.07
1984-1985   7.68
1985-1986   8.04
1986-1987   9.62
1987-1988   10.05
1988-1989   9.81
1989-1990   9.6
1990-1991   10.59
1991-1992   11.08
1992-1993   9.42
1993-1994   9.6
1994-1995   12.28
1995-1996   12.58
1996-1997   10.87
1997-1998   12.09
1998-1999   13.66
1999-2000   12.28
2000-2001   11.75
2001-2002   11.49
2002-2003   13.08
2003-2004   13.43
2004-2005   15.06
2005-2006   16.5
2006-2007   18.48
2007-2008   24.74
2008-2009   26.69"

Cal <- read.table(text = Lines, skip = 1, col.names = c("year", "wholesale price"),
  check.names = FALSE, strip.white = TRUE)