Optimisation using the constraint

Question

I have following data

Rand_Data = data.frame(x = c(  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,
        11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,
        22.,  23.,  24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,
        33.,  34.,  35.,  36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,
        44.,  45.,  46.,  47.,  48.,  49.,  50.,  51.,  52.,  53.,  54.,
        55.,  56.,  57.,  58.,  59.,  60.,  61.,  62.,  63.,  64.,  65.,
        66.,  67.,  68.,  69.,  70.,  71.,  72.,  73.,  74.,  75.,  76.,
        77.,  78.,  79.,  80.,  81.,  82.,  83.,  84.,  85.,  86.,  87.,
        88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,  96.,  97.,  98.,
        99., 100., 101., 102., 103., 104., 105., 106., 107., 108., 109.,
       110., 111., 112., 113., 114., 115., 116., 117., 118., 119., 120.,
       121., 122., 123., 124., 125., 126., 127., 128., 129., 130., 131.,
       132., 133., 134., 135., 136., 137., 138., 139., 140., 141., 142.,
       143., 144., 145., 146., 147., 148., 149., 150., 151., 152., 153.,
       154., 155., 156., 157., 158., 159., 160., 161., 162., 163., 164.,
       165., 166., 167., 168., 169., 170., 171., 172., 173., 174., 175.,
       176., 177., 178., 179., 180., 181., 182., 183., 184., 185., 186.,
       187., 188., 189., 190., 191., 192., 193., 194., 195., 196., 197.,
       198., 199), Val = c( 0.00000000e+00,  7.75383362e-02, -1.34275505e+00, -1.65497205e+00,
       -1.18564171e+00,  1.77660878e+00,  1.84593088e+00,  1.71963694e+00,
        8.21964340e-01,  1.74041137e+00,  2.64637861e+00,  2.06965750e+00,
        1.02828336e+00, -5.00395684e-01, -6.93668766e-02,  4.22030196e-01,
        1.16198422e-01, -5.11722287e-04, -6.57944961e-01, -3.74434142e-01,
       -1.15827213e+00, -1.45830766e+00, -4.90458625e-01,  2.47042268e-02,
        6.93327973e-01,  1.30938478e+00, -4.06142808e-01, -4.61168289e-01,
       -6.72137769e-01, -4.84899922e-01, -8.29915365e-01, -2.41828942e+00,
       -1.20765939e+00,  5.76353406e-01,  2.61955788e+00,  2.98795810e+00,
        4.16904195e+00,  2.93858016e+00,  5.35527961e+00,  5.11840962e+00,
        5.51236637e+00,  5.88671081e+00,  5.38415506e+00,  6.32405264e+00,
        5.65838016e+00,  5.07960498e+00,  5.70407219e+00,  7.29978260e+00,
        6.08456290e+00,  7.43051089e+00,  8.55807191e+00,  8.52981505e+00,
        8.23329180e+00,  9.10450577e+00,  9.85288281e+00,  9.62989249e+00,
        8.69069599e+00,  8.34052238e+00,  7.19647733e+00,  7.24917949e+00,
        6.63422915e+00,  4.95911711e+00,  3.49654704e+00,  1.85980401e+00,
        3.08957065e+00,  4.24267262e+00,  2.38508285e+00,  1.02089838e+00,
        2.19670308e-01,  1.53858340e+00,  2.41653182e+00,  3.69861411e+00,
        2.87607870e+00,  3.21819643e+00,  2.07359072e+00,  1.13000323e+00,
        1.17047772e+00,  1.55913045e+00,  8.41514076e-01,  2.23618634e+00,
        2.99905073e+00,  3.21259672e+00,  4.39274819e+00,  1.99515000e+00,
        1.16175606e+00,  1.70758890e+00,  9.74160444e-01, -3.03140114e-02,
        1.17659213e+00,  3.07232323e-01,  1.08034381e+00,  1.29199568e+00,
        2.45607387e-01, -5.21021566e-01, -4.34639402e-01, -9.60127469e-01,
       -1.03258482e+00, -1.43325172e+00, -2.28217823e+00, -2.09923712e+00,
       -2.06051134e+00, -2.19767637e+00, -2.16617856e+00, -2.45811083e+00,
       -2.79865309e+00, -3.70529459e+00, -2.77461502e+00, -3.27216972e+00,
       -4.79844244e+00, -3.42557001e+00, -4.14314358e+00, -5.25738457e+00,
       -3.49827107e+00, -3.55643468e+00, -5.45626633e+00, -5.33106305e+00,
       -5.25332322e+00, -4.94845946e+00, -5.04795700e+00, -4.75058961e+00,
       -3.37815392e+00, -3.24099974e+00, -4.26774262e+00, -4.35293178e+00,
       -3.70409130e+00, -4.85062599e+00, -4.46039428e+00, -4.73719944e+00,
       -5.53392228e+00, -6.15049133e+00, -6.10977286e+00, -4.87411131e+00,
       -4.49530898e+00, -4.66802205e+00, -6.98105638e+00, -5.77006061e+00,
       -6.27152149e+00, -6.07877879e+00, -5.09850276e+00, -5.68772276e+00,
       -5.16446602e+00, -3.58115497e+00, -3.21209072e+00, -3.61179089e+00,
       -4.40163561e+00, -3.14328137e+00, -3.80444247e+00, -3.91338929e+00,
       -3.98395412e+00, -4.71086150e+00, -2.98308294e+00, -2.51608443e+00,
       -2.86848586e+00, -3.27653288e+00, -3.17684302e+00, -2.88103208e+00,
       -4.90737736e+00, -4.27541967e+00, -2.92691377e+00, -2.19142665e+00,
       -3.08185258e+00, -2.54393704e+00, -2.42485142e+00, -1.90176762e+00,
       -6.52372476e-01, -1.04169430e+00, -5.79717797e-01, -4.03696696e-01,
        5.92495881e-01, -2.14694942e+00, -1.42679299e+00, -9.81665751e-01,
       -3.51830293e+00, -2.24672044e+00, -2.23578918e+00, -9.44689326e-01,
       -1.55913248e+00, -6.17618042e-01, -9.56082794e-01, -1.38326771e+00,
       -1.55108670e+00, -5.74061887e-01,  1.89141995e-01,  8.11432327e-01,
        1.99574115e+00,  1.72686176e+00,  7.17859558e-01,  5.03518805e-01,
        1.91256417e+00,  9.59411484e-01,  5.28613488e-01,  1.17086649e+00,
        1.88075943e+00,  2.41999514e+00,  2.51115680e+00,  4.77976432e+00,
        6.14323607e+00,  6.42564314e+00,  4.82006996e+00,  4.63386067e+00))

And, below is my function (for simplicity, I just consider a Cubic spline fitting)

MyFN = function(a, b, c, d) return(a * Rand_Data$x^3 + b * Rand_Data$x^2  + c * Rand_Data$x + d)

Now I want to estimate the parameters a, b, c, d that minimises the sum of absolute difference between Rand_Data$Val and the fitted value with below constraints

1. Fitted value for the 100th position from MyFN should be less than 2 i.e. MyFN[100] < 2
2. The difference between sum of fitted value and sum of observed MyFN$Val should be less that 0.001

I wonder how can I setup optim() function so that I can estimate a,b,c,d

Here, my data and function and just for example, so I am looking for some generic approach to optimally estimate the parameters.

Answer 1

Since your post Optimisation using the constraint giving very different results from different packages is closed, I think we can continue the investigation in this post and take a look how two approaches behavior differently.

---

Below is the code with your updated example data, including the verification and visualizations as well

set.seed(1)
A <- 1.34
B <- 0.5673
C <- 6.356
D <- -1.234
x <- seq(0.5, 20, length.out = 500)
y <- A + B * x + C * x^2 + D * log(x) + runif(500, 0, 3)
Rand_Data <- data.frame(x = x, Val = y)


library(pracma)
library(nloptr)

eps <- 1e-4
X <- cbind(1, x, x^2, log(x))
f <- function(theta) {
  sum(abs(X %*% theta - y))
}
A <- X[100, , drop = FALSE]
b <- 120 - eps
hin <- function(theta) {
  abs(sum(X %*% theta) - sum(y)) - 1e-3 + eps
}

# use fmincon
p <- fmincon(rep(0, 4), f, A = A, b = b, hin = hin)
sol1 <- p$par

# use nloptr
Hx <- function(theta) {
  X[100, , drop = FALSE] %*% theta - (120 - eps)
}
q <- nloptr(rep(0, 4), f, eval_g_ineq = hin, eval_g_eq = Hx, opts = list("algorithm" = "NLOPT_LN_COBYLA", "xtol_rel" = 1.0e-8))
sol2 <- q$solution


# verify the solution
res <- transform(
  Rand_Data,
  Val_est1 = X %*% sol1,
  Val_est2 = X %*% sol2
)
check_report <- with(
  res,
  list(
    fmincon = list(
      sol = sol1,
      cost = f(sol1),
      valid = abs(sum(Val) - sum(Val_est1)) < 1e-3 & Val_est1[100] < 200
    ),
    nloptr = list(
      sol = sol2,
      cost = f(sol2),
      valid = abs(sum(Val) - sum(Val_est2)) < 1e-3 & Val_est2[100] < 200
    )
  )
)

# plot
plot(x, log(res$Val), type = "p")
lines(x, log(res$Val_est1), col = "red")
lines(x, log(res$Val_est2), col = "blue")
legend(15, 4,
  legend = c("val", "val_est1 (fmincon)", "val_est2 (nloptr)"),
  col = c("black", "red", "blue"),
  pch = c(1, NA, NA),
  lty = c(NA, 1, 1)
)

When you check the report, i.e., check_report, you will see that nloptr failed to meet the requirements

> check_report
$fmincon
$fmincon$sol
[1]   0.2215526   6.9993834   6.1700956 -19.3708030

$fmincon$cost
[1] 1328.184

$fmincon$valid
[1] TRUE


$nloptr
$nloptr$sol
[1]  0.3655464 -0.9045971  6.4630681  0.1513615

$nloptr$cost
[1] 2426.848

$nloptr$valid
[1] FALSE

and the estimation curves are shown as below