Errors in using nls function in R

Question

I have been trying to use the nls function in R on a temperature model. I have to estimate 2 parameters called s0 and Tch.

I have imported a set of data from Excel, in which I have 4 columns of 9998 rows : Date, Treelle, Tretrelevee, and Tcons. I am trying to understand how a heating controller in a building. I want to calculate the temperature (Tdep) tanks to a temperature setpoint (Tcons) and the return tempeature (Tretrelevee) and compare it with the measured temperature (Tdepreelle).

Here is a sample of my data :

   donneesbis <- structure(list(Tcons = c(51.6099999999998, 53.1362499999993, 
                                   45.4742499999993, 44.0543749999995, 57.8999999999997, 49.7168333333327, 
                                   45.7269999999998, 46.0214583333331, 53.1855833333333, 46.9418333333326, 
                                   52.2359166666666, 45.3108333333328, 44.2624999999995, 44.1653749999997, 
                                   47.3951666666659), Tret.relevee = c(43.764166666666, 44.8750000000001, 
                                                                       37.9016666666671, 36.2000000000002, 46.5000000000003, 41.514999999999, 
                                                                       39.9133333333328, 39.6116666666664, 40.2916666666656, 41.8299999999993, 
                                                                       41.1224999999992, 36.5, 36.3250000000002, 43.8758333333327, 38.2700000000001
                                   ), Treelle.30sec = c(47.4999999999989, 48.6250000000003, 40.7016666666667, 
                                                        38.6125000000001, 50.5999999999987, 44.7766666666672, 42.6600000000008, 
                                                        42.5116666666667, 43.7849999999998, 44.6649999999998, 44.1483333333336, 
                                                        38.8999999999999, 38.7250000000001, 47.9008333333323, 41.1833333333334
                                   )), .Names = c("Tcons", "Tret.relevee", "Treelle.30sec"), row.names = c(8478L, 
                                                                                                           7231L, 6122L, 3466L, 9721L, 5064L, 1857L, 2348L, 1052L, 4575L, 
                                                                                                           1352L, 3653L, 3496L, 8654L, 6429L), class = "data.frame")

Here are my functions. They are recursive, which makes it more complicated :

Tdep(t) = sigma(t-1) (Tch - Tretrelevee(t)) + Tretrelevee(t)

where sigma(t) = ((100-2*s0)/32)1.25(Tcons(t)-Tdep(t)) + s0.

In addition, sigma should not be 100-s0.

It is clearer written like this, but in the code I only wrote one function : Tdep(t) = (((100-2*s0)/32)1.25(Tcons(t-1)-Tdep(t-1)) + s0)*(Tch - Tretrelevee(t)) + Tretrelevee(t)

The parameters I want to estimate are s0 and Tch.

Here is the code :

Tdep <-  function(s0,Tch,tps) {
   if (tps == 1)     {return (37.9)}
   else if ((((100-2*s0)/32)*1.25*(donneesbis$Tcons[tps-1]-Tdep(s0,Tch,tps-1))+s0) <= s0) {return ((s0/100)*(Tch - donneesbis$Tret.relevee[tps])+donneesbis$Tret.relevee[tps])}
   else if ((((100-2*s0)/32)*1.25*(donneesbis$Tcons[tps-1]-Tdep(s0,Tch,tps-1))+s0) >= 100 - s0) {return (((100-s0)/100)*(Tch - donneesbis$Tret.relevee[tps])+donneesbis$Tret.relevee[tps])}
   else {return (((((100-2*s0)/32)*1.25*(donneesbis$Tcons[tps-1]-Tdep(s0,Tch,tps-1))+s0)/100)*(Tch - donneesbis$Tret.relevee[tps])+donneesbis$Tret.relevee[tps])}}

With this formula, I can calculate each value one by one, but I cannot calculate a serie of values :

> x <- 1:10
> Tdep(15,60,x)
[1] 37.9
Warning message:
In if (tps == 1) { :
  the condition has length > 1 and only the first element will be used

Though I have :

> Tdep(15,60,15)

[1] 41.5295

Therefore when I use nls :

x <- 1:15 
nls(donneesbis$Treelle.30sec~ Tdepbis(s0,Tch0,x), data = donneesbis, start = list(s1=15,Tch1=60))

I get this error in return :

Error in qr(.swts * attr(rhs, "gradient")) : 
  dims [product 2] do not match the length of object [9998]

And this warning message :

In addition: Warning messages:
1: In if (tps == 1) { :

Thanks to research on this website, I also tried with the "ifelse" function :

Tdepbis <-  function(s0,Tch,tps) {
  ifelse(tps <= 1,37.9,{if ((((100-2*s0)/32)*1.25*(donneesbis$Tcons[tps-1]-    Tdepbis(s0,Tch,tps-1))+s0) <= s0) {return ((s0/100)*(Tch - donneesbis$Tret.relevÃ.e[tps])+donneesbis$Tret.relevÃ.e[tps])}
  else if ((((100-2*s0)/32)*1.25*(donneesbis$Tcons[tps-1]-Tdepbis(s0,Tch,tps-1))+s0) >= 100 - s0) {return (((100-s0)/100)*(Tch - donneesbis$Tret.relevÃ.e[tps])+donneesbis$Tret.relevÃ.e[tps])}
  else {return (((((100-2*s0)/32)*1.25*(donneesbis$Tcons[tps-1]-Tdepbis(s0,Tch,tps-1))+s0)/100)*(Tch - donneesbis$Tret.relevÃ.e[tps])+donneesbis$Tret.relevÃ.e[tps])}
  })
}

And, using nls the same way as before, I get this error message :

Error in donneesbis$Tcons[tps - 1] : 
  only 0's may be mixed with negative subscripts

I have trouble understanding what the first error is, since I have checked that Treelle and Tdep had the same length, and I can't solve the second error.

I could not fine any solution that couls solve my problem on forum or in R help.

Would someone have any idea ?

Thank you !

Answer 1

You need an efficient function. Thus, I would implement it in Rcpp (check if the formulas are correct):

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
NumericVector temps(const NumericVector Tcons,
                    const NumericVector Tretrelevee,
                    const double s0,
                    const double Tch,
                    const double Tdep0) {
  int n = Tcons.size();
  NumericVector res(n);
  double sigma;
  res(0) = Tdep0;
  for (int t = 1; t < n; t++) {
    sigma = (((100-2*s0)/32)*1.25*(Tcons(t-1)-res(t-1)) + s0);
    if (sigma < s0) sigma = s0;
    if (sigma > (100 - s0)) sigma = (100 - s0);
    res(t) = sigma * (Tch - Tretrelevee(t)) + Tretrelevee(t);
  }

  return res;
}

(Follow these instructions if you have never used it. Briefly, you need to install Rtools if you use Windows and then open a C++ file in RStudio. Copy this code into it and source it.)

Then you can do this in R:

fit <- nls(Treelle.30sec ~ temps(Tcons, Tret.relevee, s0, Tch, 37.9), 
      data = donneesbis, start = list(s0 = 1, Tch = 1))
summary(fit)  
#Formula: Treelle.30sec ~ temps(Tcons, Tret.relevee, s0, Tch, 37.9)
#
#Parameters:
#    Estimate Std. Error t value Pr(>|t|)    
#s0  100.1629     0.2307 434.093   <2e-16 ***
#Tch  21.4532    27.1098   0.791    0.443    
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 2.673 on 13 degrees of freedom
#
#Number of iterations to convergence: 4 
#Achieved convergence tolerance: 1.533e-08

plot(donneesbis$Treelle.30sec)
lines(predict(fit), col = "red")

No idea why you chose 37.9 as the starting temperature. 47.5 seems like the obvious choice.