[R-sig-dyn-mod] Problem with dede's lagvalue (deSolve)
Nicola Gambaro
n|co|@ @end|ng |rom g@mb@ro@co@uk
Tue Apr 12 16:53:34 CEST 2022
Dear Thomas,
Thank you for your reply and apologies for the delay. I wanted to take the time to warmly thank you for your very informative and prompt replies to this question and another recent one from me.
I think I have now understood what the problem is!
Kind regards,
Nicola
> On 15 Mar 2022, at 12:00, r-sig-dynamic-models-request using r-project.org wrote:
>
> Send R-sig-dynamic-models mailing list submissions to
> r-sig-dynamic-models using r-project.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> or, via email, send a message with subject or body 'help' to
> r-sig-dynamic-models-request using r-project.org
>
> You can reach the person managing the list at
> r-sig-dynamic-models-owner using r-project.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-dynamic-models digest..."
>
>
> Today's Topics:
>
> 1. Modeling event with non-instantaneous duration
> (Rafael Ayala Hernandez)
> 2. Re: Problem with dede's lagvalue (deSolve) (Thomas Petzoldt)
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Tue, 15 Mar 2022 03:16:07 +0000
> From: Rafael Ayala Hernandez <Rafael.AyalaHernandez using oist.jp>
> To: "r-sig-dynamic-models using r-project.org"
> <r-sig-dynamic-models using r-project.org>
> Subject: [R-sig-dyn-mod] Modeling event with non-instantaneous
> duration
> Message-ID: <85F4496B-E06A-458E-AA76-0F0BFD9EE064 using oist.jp>
> Content-Type: text/plain; charset="us-ascii"
>
> Dear all,
>
> I am currently using the deSolve package in my asteRisk package to perform orbit propagation of satellites.
>
> I am now aiming to model orbital maneuvers, which can be seen as sudden, transient changes in acceleration due to engine bursts.
>
> I have been able to model instantaneous bursts as events with pre-defined times. However, for a more accurate representation of such maneuvers, it would be required to model them as changes in acceleration that are present for a given period of time (corresponding to the time during which engines are fired).
>
> I would therefore like to ask if there is any recommended way of modeling this type of non-instantaneous events in ODE systems with the deSolve package or other related packages?
>
> Thanks a lot in advance
>
> Best wishes,
>
> Rafa
>
>
>
> ------------------------------
>
> Message: 2
> Date: Tue, 15 Mar 2022 08:10:59 +0100
> From: Thomas Petzoldt <thomas.petzoldt using tu-dresden.de>
> To: <r-sig-dynamic-models using r-project.org>
> Cc: <nicola using gambaro.co.uk>
> Subject: Re: [R-sig-dyn-mod] Problem with dede's lagvalue (deSolve)
> Message-ID: <0762cefd-18d3-a147-5daf-ec06b967bb12 using tu-dresden.de>
> Content-Type: text/plain; charset="utf-8"; Format="flowed"
>
> Hi Nicola,
>
> the error indicates that the history buffer is exhausted. You can
> increase the buffer by setting mxhist:
>
> y_out_2 <- dede(y = yinit, times = times_sd, func = supplydrive,
> parms = parms_sd, method = "lsoda",
> control = list(mxhist = 1e6))
>
> This is found on the ` ?dede help page.
>
> Unfortunately, this solves only part of the problem in your case,
> because then the solver seems to get stuck due to numerical problems.
> Some of the states grow exponentially to extreme values.
>
> To see this, you can put a print or cat in the model function:
>
> cat("time=", t, ", slag=", slag, ", llag=", llag, ",
> taulag=", taulag, "\n")
>
> ... or you can try a brute-force simulation with another solver, e.g.
> "vode" and a modified absolute tolerance (either a large value or zero).
> The relative tolerance remains unchanged.
>
> times_sd <- seq(0, 500, by = 1)
> y_out_2 <- dede(y = yinit, times = times_sd, func = supplydrive,
> parms = parms_sd, method = "vode",
> control = list(mxhist = 1e6), atol = 0)
>
> plot(y_out_2)
>
> The output contains then at least part of the solution and the plot
> shows where the it ran into trouble. Setting the time step to a larger
> value (e.g. by=1) speeds up simulation, because a too small time step
> just limits internal step size adaptation.
>
> Thomas
>
>
>
> On 10.03.2022 at 19:06 Nicola Gambaro wrote:
>> Dear R dynamics community
>>
>> I am trying to solve a dynamic model based on DDEs with deSolve’s dede solver. However, I am encountering problems understanding and solving errors related to the lagvalue function. Here’s the code:
>>
>> ----------------------------------------------------
>> library(deSolve)
>> supplydrive <- function(t, y, parms) {
>> with(as.list(c(y, parms)),{
>> #y[1] = demand
>> #y[2] = price
>> #y[3] = supply
>> #y[4] = extraction rate
>> #y[5] = production stock
>> #y[6] = manufacturing stock
>> #y[7] = in-use stock
>> #y[8] = waste management stock
>>
>> #define delays
>> #Extraction rate // s
>> slag <- t - s
>> if (slag <= 0) {
>> e_slag <- y[4]
>> } else {
>> e_slag <- lagvalue(slag, 4)
>> }
>> #Stocks // l
>> llag <- t - l
>> if (llag <= 0) {
>> s3_llag <- y[7]
>> s2_llag <- y[6]
>> } else {
>> s3_llag <- lagvalue(llag, 7)
>> s2_llag <- lagvalue(llag, 6)
>> }
>> #Price // tau
>> taulag <- t - tau
>> if (taulag <= 0) {
>> pricelag <- y[2]
>> dpricelag = 0
>> } else {
>> pricelag <- lagvalue(taulag, 2)
>> dpricelag <- lagderiv(taulag, 2)
>> }
>>
>> #waste and dissipation flows
>> w1 = eta*y[4] + gamma*(1-eta)*e_slag#delay s<---
>> w4 = theta*(alpha4*s3_llag + alpha7*s2_llag)#delay l
>> d1 = lambda*y[5]
>> d2 = lambda*y[6]
>> d3 = lambda*y[7]
>> d4 = lambda*y[8]
>> totald = d1+d2+d3+d4
>>
>> #supply, demand, price
>> dy1 = r*y[1]
>> dy2 = c*y[2]*((y[1]-y[3])/y[1])
>> dy3 = y[4] - (w1 + w4 + totald)
>>
>> #extraction rate
>> if (dpricelag > 0)
>> dy4 = y[4] * (dpricelag/pricelag) #delay tau
>> else
>> dy4 = z*dy1
>>
>> #stocks
>> dy5 = y[4] - alpha1*y[5] + alpha2*y[6] + alpha5*y[8] - w1 - d1
>> dy6 = alpha1*y[5] - alpha2*y[6] - alpha3*y[6] + alpha6*y[8] - alpha7*y[6] - d2
>> dy7 = alpha3*y[6] - alpha4*y[7] - d3
>> dy8 = alpha4*y[7] - alpha5*y[8] - alpha6*y[8] + alpha7*y[6] - w4 - d4
>>
>> list(c(dy1, dy2, dy3, dy4, dy5, dy6, dy7, dy8))
>> })
>> }
>>
>> yinit <- c(1, 1, 1, 1, 1, 1, 1, 1)
>> times_sd <- seq(0, 200, by = 0.1)
>> parms_sd <- c(alpha = c(0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1),
>> tau = 120, r = 0.1, c = 0.8, s=2, l=1,
>> eta=0.2, gamma=0.3, theta=0.5,
>> lambda=0.01, z=0.2)
>> y_out_2 <- dede(y = yinit, times = times_sd, func = supplydrive, parms = parms_sd, method = "lsoda”)
>>
>> —————————————————
>>
>> This runs fine. However, if the times sequence length is increased to, say,
>>
>> times_sd <- seq(0, 500, by = 0.1) ,
>>
>> then I get the following error
>>
>> Error in lagvalue(taulag, 2) :
>> illegal input in lagvalue - lag, 355.98, too large, at time = 475.98
>>
>> However, this shouldn’t be the case because the conditional statements in the function should make sure that the lagtime is not to a time before the beginning of the simulation (t=0). Is this a memory issue? How could I go about solving this issue?
>>
>> Many thanks,
>> Nicola
>> _______________________________________________
>> R-sig-dynamic-models mailing list
>> R-sig-dynamic-models using r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>
> --
> Dr. Thomas Petzoldt
> https://tu-dresden.de/Members/thomas.petzoldt
>
>
>
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>
>
> ------------------------------
>
> End of R-sig-dynamic-models Digest, Vol 168, Issue 1
> ****************************************************
More information about the R-sig-dynamic-models
mailing list