0
$\begingroup$

Similar to my question at Stackoverflow I want to calculate curve fitting coefficents of a non-linear function with some nasty inequality constraints. The results from my linked question are quite fine, but I have to adept it to come closer to my "real-world" problem.

The goal is the same, the constraints are different

  • Minimize the residuals of the fitted curve.
  • Take the non-linear inequality constraint during fitting into account

The predicted function of the curve is

$y_p=f(x)=a_1\cdot(1-a_2\cdot x^{a_{3}}\cdot (x-1))$

In my opinion, the cost function could be approached with least squares

$E=\sum(y_m-f(x))^2$

The constraint $h(x)$ can be written as

$\frac{V_2}{(1-a_2\cdot x_2^{a_{3}}\cdot (x_2-1))}-\frac{V_1}{(1-a_2\cdot x_1^{a_{3}}\cdot (x_1-1))}\leq0$

with set values of $V_2, V_1, x_2$ and $x_1$

Sample code in R

library(tidyverse)
library(plotly)
#library(NlcOptim)
#library(Rsolnp)

## define set curve function ----
outerSolve <- function(x, a){

  a[1] * (1 - innerSolve(x, a) * (x - 1))

}

# part of the curve function
innerSolve <- function(x, a){

  a[2] * x ^ a[3]

}

checkConstraint <- function(a){

  checkValue <- Vcon[2] / (1 - innerSolve(xcon[2], a) * (xcon[2] - 1)) -
    Vcon[1] / (1 - innerSolve(xcon[1], a) * (xcon[1] - 1))

  if(checkValue <= 0) {

    return(list(status = TRUE, value = checkValue))

  } else {

    return(list(status = FALSE, value = checkValue))
  }

}

# create sample data
sampleData <- (function(){

  # set x range
  x = seq(50, 800, 5)

  # set coefficents
  setCoef <- c(80, 0.004, -0.3)

  # calculate y with some noise and create dataframe
  set.seed(123)
  tmpDF <- tibble(
    x = x,
    y = outerSolve(x, setCoef) + rnorm(length(x), 0, 0.02)
         )

  p <- plot_ly(x = ~ x, y = ~ y, data = tmpDF, type = 'scatter', mode = 'markers')

  return(list(setCoef = setCoef, raw = tmpDF, plot = p))

})()

# start values for target coefficents
startCoef <- c(max(sampleData$raw$y), 0.0001, -0.1)

# curve fitting with nls
refCurveFitting <- (function(){

  tmpModel = nls(
  y ~ outerSolve(x, c(find.a1, find.a2, find.a3)),
  data = sampleData$raw,
  start = c(
    find.a1 = startCoef[1],
    find.a2 = startCoef[2],
    find.a3 = startCoef[3]
  ),
  trace = TRUE,
  control = nls.control(maxiter = 100)
)

return(list(model = tmpModel, gof = tibble(AIC = AIC(refCurveFitting$model),
                                       BIC = BIC(refCurveFitting$model))
            )
       )

})()

# set constraint values
Vcon <- c(150, 135)
xcon <- c(50, 200)

with the set

  • coefficients $a_2$, $a_3$ (setCoef) and
  • values for $V_2, V_1, x_2$ and $x_1$ (Vcon, xcon)

the constraint is not fulfilled.

checkConstraint(sampleData$setCoef)
$status
[1] FALSE

$value
[1] 1.497501

I have saved the nls model (curve fitting) as reference in refCurveFitting$model. I have created the cost function but I was not able to find a optimiziation routine or package which provides me enough help to guide me through.

minFunction <- function(a){
sampleData$raw %>%
      mutate(predicted.y = outerSolve(x, a),
             squaredRes = (y - predicted.y) ^ 2
      ) %>%
      summarise(sum.squaredRes = sum(squaredRes)) %>%
      pull(sum.squaredRes)
}

What's the best way to proceed?


Edit

Some real data, could be imported into R with

read.table("clipboard", header = TRUE, dec = ",")

"xm" "ym"
63,6    70,66
66,73   70,61
69,75   70,61
72,8    70,6
75,67   70,54
79,46   70,51
82,21   70,46
85,8    70,38
88,9    70,27
92,52   70,2
95,42   70,13
99,1    70,08
102,43  69,99
106,1   69,89
109,49  69,82
113,13  69,71
116,03  69,62
119,78  69,53
123,22  69,48
126,73  69,39
130,4   69,3
133,97  69,22
138,16  69,15
141,66  69,11
144,98  69,04
148,48  68,99
152,2   68,96
155,53  68,89
159,24  68,81
162,88  68,77
166,47  68,71
169,9   68,67
173,17  68,6
176,98  68,57
179,91  68,52
183,46  68,47
187,11  68,42
190,57  68,33
194,12  68,29
197,77  68,27
201,26  68,17
204,82  68,14
208,38  68,1
211,88  68,05
214,97  67,97
218,32  67,92
221,9   67,91
225,61  67,83
229,63  67,76
232,67  67,71
236,24  67,7
239,71  67,64
242,96  67,6
246,84  67,55
250,34  67,5
254,03  67,44
257,55  67,39
260,59  67,34
264,02  67,32
267,7   67,25
270,91  67,2
274,58  67,12
278,25  67,1
281,76  67,06
285,31  67,05
288,84  66,98
291,66  66,95
295,55  66,91
299,09  66,84
302,76  66,77
306,12  66,72
309,68  66,71
312,32  66,73
316,16  66,66
320,19  66,58
323,16  66,54
326,88  66,5
330,42  66,46
333,78  66,37
337,47  66,37
340,99  66,31
344,57  66,28
348,07  66,23
350,96  66,19
354,53  66,16
358,11  66,1
361,71  66,05
365,2   65,99
368,96  65,96
372,41  65,91
375,9   65,87
379,54  65,81
382,38  65,79
385,95  65,75
389,42  65,76
393,72  65,68
396,85  65,58
400,17  65,57
403,78  65,53
407,2   65,47
410,81  65,46
413,94  65,42
418,03  65,4
420,97  65,32
424,64  65,3
428,22  65,23
431,59  65,22
435,06  65,21
438,39  65,14
441,94  65,11
445,38  65,06
448,94  65,07
452,33  65,02
456,04  64,99
459,71  64,92
463,36  64,89
466,6   64,82
470,25  64,81
473,36  64,72
477,38  64,7
480,1   64,72
483,7   64,63
487,33  64,64
490,92  64,58
494,45  64,56

with the found coefficients $a_2$ = 0.00322 and $a_3$ = -0.41534 the values for the constraints are

$V_2$ = 143.5972, $V_1$ = 153.399, $x_2$ = 48.242 and $x_1$ = 338.442

I will go on with a modified Levenberg-Marquardt nonlinear least-squares algorithm.

$\endgroup$
  • 3
    $\begingroup$ Your basic problem is that the data and model are inconsistent with the constraint. Something has to give in a major way. Exactly how do you want to balance out the costs of violating the constraint with the costs of producing a curve that obviously is a bad fit to the data? $\endgroup$ – whuber May 17 at 16:05
  • $\begingroup$ I have slightly adjusted the set constraint values, to get a TRUE status check with checkConstraint(c(78.5294, 0.0028644, -0.299623)). The values are the result after 6 nls iteration steps. Violating the constraints isn't an option, it is some kind of physical backup/rule. The data itself is not that steady and can not cover the physical relation completely. So, if I have a higher residuals as a tradeoff, I will take it. $\endgroup$ – Patrick May 17 at 18:28
  • $\begingroup$ That's a start. But if the constraint is known to hold, then there must be something inherently wrong with the model. You have enough data to explore that and to modify the model accordingly. Is it possible your synthetic example is unrealistic? It would help to see realistic data. $\endgroup$ – whuber May 17 at 19:03
  • $\begingroup$ You could try a transformation into an unconstrained problem. Let $\zeta_i:= 1 - a_2 x_i^{a_3} (x_i - 1)$ for $i=1$, $2$ which as I understand must be positive. You can re-parameterize the model using $[a_1, \zeta_1, \, \zeta_2]$, and the constraint becomes $\zeta_1/\zeta_2 \leqslant V_1 / V_2$. Then you can re-parameterize again so that the three parameters are unconstrained. For instance $\zeta_1 = e^\lambda$ $\zeta_2 = e^\lambda (V_2/V_1 + e^\mu)$ with $\lambda$ and $\mu$ unconstrained. The convergence may become more difficult to get, though. $\endgroup$ – Yves May 21 at 19:18
  • $\begingroup$ Thanks for the input, I will give it a try. I'm pretty happy with the Levenberg-Marquardt approach. After the final iteration step, I'm searching for the lowest residual and a true condition for the constraints. For the sample case the RSS value of the final step is round about 0.2 and 0.8 for a true constraint condition. I'm aware this is not the optimal way of doing it, but in this case I think it's a good tradeoff. $\endgroup$ – Patrick May 22 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.