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 at r-project.org
You can reach the person managing the list at
r-sig-dynamic-models-owner at 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 at oist.jp>
To: "r-sig-dynamic-models at r-project.org"
<r-sig-dynamic-models at r-project.org>
Subject: [R-sig-dyn-mod] Modeling event with non-instantaneous
duration
Message-ID: <85F4496B-E06A-458E-AA76-0F0BFD9EE064 at 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 at tu-dresden.de>
To: <r-sig-dynamic-models at r-project.org>
Cc: <nicola at gambaro.co.uk>
Subject: Re: [R-sig-dyn-mod] Problem with dede's lagvalue (deSolve)
Message-ID: <0762cefd-18d3-a147-5daf-ec06b967bb12 at 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 at 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 at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
------------------------------
End of R-sig-dynamic-models Digest, Vol 168, Issue 1
****************************************************