如何最小化具有 R 线性等式约束的非线性 objective 函数?
How can I minimize nonlinear objective function with linear equality constraint with R?
最小化
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<=I<=N)。 J 是季度值,K 是第 i 个日历年在 i-1 会计年度中的期间数。
在我的例子中,我 = 39,J = 4,K = 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
表述似乎有问题。只看图中的公式,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
个元素的连续差异的平方和,该差异受前 4 个元素影响z 的元素求和为 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)
最小化
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<=I<=N)。 J 是季度值,K 是第 i 个日历年在 i-1 会计年度中的期间数。 在我的例子中,我 = 39,J = 4,K = 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
表述似乎有问题。只看图中的公式,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
个元素的连续差异的平方和,该差异受前 4 个元素影响z 的元素求和为 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)