How can I minimize nonlinear objective function with linear equality constraint with R?

Viewed 83

enter image description here

Minimize
f(z) = sum_(t=2)^IJ (Z_t - Z_t-1)^2

Subject to constraint
sum_(j=1)^J (Z_(i-1)J+j+k ) = y^f_i, i = 1,....., I-1.

This optimization find out quarterly values from fiscal year data series (y^f_i) and then sum those quarterly value to find out annual value. I is the number of calendar years considered in the series interval (2<=I<=N). J is quarter value and K is the number of period of calendar year i which are in fiscal year i-1. In my case, I = 39, J = 4, K = 2

How can I solve this problem using R?

The way I tried to write code is provided below:

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)

Here is the data I have:

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

1 Answers

There seems to be something wrong with the formulation. Just looking at the formulas in the image, the last constraint for i=I-1=39-1=38 sums z elements (38-1)*4 + 1 + 6, (38-1)*4 + 2 + 6, (38-1)*4 + 3 + 6 and (38-1)*4 + 4 + 6 which is elements 155 156 157 158 but z goes from 1 to 4*39 and so has only 156 elements. Furthermore not all the z values participate in a constraint.

Given the problems cited let us change the problem to o ne that makes sense and assume we want to minimize the sum of the squares of the successive differences of the I*J elements of z subject to the first 4 elements of z summing to Cal[1, 2], the next 4 summing to Cal[2, 2], and so on up to the last 4 elements of z summing to Cal[39, 2]. In that case we can write the Aeq constraint matrix as a block diagonal matrix using the kronecker product shown. We ignore K. (Cal is shown reproducibly in the Note at the end.)

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)

giving

> 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 ...

Note

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)
Related