ODE with constants that vary through time

Modeling
1

Hi Stanimals and ODE experts,

I have a toy ODE problem that has two equations with time varying constants. I’ve coded up the model where the constants for the first 3 time points change and the rest are 0. This is just an example to show that I can kind of hack it to work, however, I’m going to have much longer time spans and lots of changing constants. You can think of the constants measuring a known exogenous input at the given time that disturbs the system as written. Is there a way to get this to work in the ideal case I present below?

functions{
  real[] decay(real t, real [] y, real [] theta, 
                       real [] x_r, int[] x_i){
                         
    real dydt[2];
    real theta2;
    
    if(t <= 1){
      theta2 = theta[2] + theta[3] * x_r[1];
      dydt[1] = y[2] * x_r[1] - theta[1] * y[1]; 
      dydt[2] = -theta2 * y[2]  + (1 - x_i[1]) * theta[1] * (1 - y[2]) ; 
    } else if (t <= 2 && t > 1){
      theta2 = theta[2] + theta[3] * x_r[2];
      dydt[1] = y[2] * x_r[2] - theta[1] * y[1]; 
      dydt[2] = -theta2 * y[2]  + (1 - x_i[2]) * theta[1] * (1 - y[2]) ; 
    } else if (t <= 3 && t > 2){
      theta2 = theta[2] + theta[3] * x_r[3];
      dydt[1] = y[2] * x_r[3] - theta[1] * y[1]; 
      dydt[2] = -theta2 * y[2]  + (1 - x_i[3]) * theta[1] * (1 - y[2]) ; 
    } else {
      theta2 = theta[2];
      dydt[1] = - theta[1] * y[1]; 
      dydt[2] = -theta2 * y[2]  + theta[1] * (1 - y[2]) ; 
    }
   
    return dydt;
  }
}
data {
  real theta[3];
  int N_t;
  real times[N_t];  
  real C0[2];
  real r[N_t];
  int x[N_t];
}
parameters{
}
transformed parameters{
  real C[N_t, 2];
  C = integrate_ode_rk45(decay, C0, 0, times, theta, r, x);
}
model{
}
generated quantities{
  real stock[N_t];
  real quality[N_t];

  stock = C[, 1];
  quality = C[, 2];
}

Ideally, I would just be able to call the time-varying constants using t or some iterator that is indexed off of the number of time points. Something like (I know this doesn’t work since t is a real and indexes must be int):

functions{
  real[] decay(real t, real [] y, real [] theta, 
                       real [] x_r, int[] x_i){
                         
    real dydt[2];
    real theta2;
  
      theta2 = theta[2] + theta[3] * x_r[t];
      dydt[1] = y[2] * x_r[t] - theta[1] * y[1]; 
      dydt[2] = -theta2 * y[2]  + (1 - x_i[t]) * theta[1] * (1 - y[2]) ; 
    }
   
    return dydt;
  }
}

To run the first ODE sample with dummy data you can use:

library(rstan)
library(data.table)

r <- c(10, 40, 35, rep(0, 12))
input_data <- list(
  theta = c(stock_decay = 0.05, quality_wearout = 0.03, repetition_wearout = 0.01),
  N_t = length(r), 
  times = seq(from = 1, to = length(r)),
  C0 = c(0, 1),
  r = log(1 + r),
  x = (r != 0)*1
  )

m_test <- stan_model(file = 'de_test.stan')

samples <- sampling(m_test, data=input_data, 
                    algorithm = "Fixed_param", 
                    iter = 1, chain = 1)

X <- extract(samples)
out <- data.frame(
  time = input_data$times,
  stock = c(X[["stock"]]),
  quality = c(X[["quality"]])
  )

data <- as.data.table(out)

With output:

    time    stock   quality
 1:    1 2.276421 0.9474520
 2:    2 5.483019 0.8859394
 3:    3 8.211380 0.8294966
 4:    4 7.810937 0.8137738
 5:    5 7.429992 0.7992599
 6:    6 7.067626 0.7858618
 7:    7 6.722934 0.7734944
 8:    8 6.395054 0.7620780
 9:    9 6.083164 0.7515390
10:   10 5.786483 0.7418099
11:   11 5.504273 0.7328291
12:   12 5.235827 0.7245390
13:   13 4.980474 0.7168863
14:   14 4.737573 0.7098217
15:   15 4.506519 0.7033003

Thanks!

2 Likes
2

This looks like a job for 1d interpolation, which isn’t currently easy to do. Check this out though: Forced ODEs, a start for a case study?.

There might be problems with the ODE solver changing constants discontinuously like that, so you should probably use the interpolation stuff in that forced ODEs thread.

1 Like
3

In my StanCon 2018 submission I ended up implementing a custom leapfrog integrator precisely to have time-varying interpolated input (I used splines). For my case it worked fine, but my ODEs were very much not stiff. For more challenging ODEs doing the interpolation within the ODE call would be a better approach. Since leapfrog is easy to implement, I have first tested whether it integrates nicely in R (contrasting to the deSolve package) which lets you iterate quickly (Euler’s rule didn’t cut it for me).

Hope that helps!

4 Likes
4

Thanks, I used the code provided as inspiration and simplified it tremendously for my problem. I get the same answer as before without all the if-elses. I’ll have to check out if it still works to solve the “inverse” problem.

functions{
    int find_interval_elem(real x, vector exo) {
      int N = num_elements(exo);
      real right = N - x;
      int iter = 0;
      
      if(N == 1) return(1);
      else if(0 == right) return(N);
      else {
        while(right > 1) {
            right += -1; 
            iter += 1;
        }
      }
      return( N - iter );
    }
    
  real[] decay(real t, real [] y, real [] theta, real [] x_r, int[] x_i){                  
    real dydt[2];
    int d = find_interval_elem(t, to_vector(x_r));
    real theta2 = theta[2] + theta[3] * x_r[d];
    int ind = d > 0 ? 1 : 0;

    dydt[1] = y[2] * x_r[d] - theta[1] * y[1]; 
    dydt[2] = -theta2 * y[2]  + (1 - ind) * theta[1] * (1 - y[2]) ; 
   
    return dydt;
  }
}