触发根函数时更改 R 的 deSolve 中的参数值
Changing value of parameters in R's deSolve when root function is triggered
短版
根据评论中的建议，我留下了我在此处发现的问题的简短摘要版本。原文完整的解释可以在下面找到。
总而言之，我目前正在使用 deSolve 包，并且有兴趣实现响应根函数的事件。我知道我可以使用此类事件来引入模型状态变量的突然变化，但我还想修改模型的参数值以响应此类根函数。
这可能吗？
长版
我已经在 R 中实现了一个轨道数值传播器（一个计算 space 卫星位置的函数，给定初始位置和速度状态）。问题被表述为一组 6 个常微分方程 (X , 位置和速度的 Y 和 Z 分量）。我实现的模型计算任意给定时间的加速度，然后我使用 deSolve 包执行积分并计算轨迹。
进行此类计算时必须确定的一个关键参数是参考系的中心，参考系的中心通常位于对卫星施加最显着引力影响的天体的质心。这是因为，尽管原则上可以使用任意参考系进行积分和计算轨迹，但实际上只有当坐标中心位于施加主要引力影响的天体上时，我们才能获得合理的结果（即地球轨道卫星的地球，月球轨道卫星的月球等等），如讨论的那样in this SE question。
最初，我的实现使用了一个恒定的坐标中心，由用户提供或根据不同主要天体的影响范围自动确定。
但是，这不适用于行星际 space 任务的建模，因为施加主要引力影响的天体在轨迹期间会发生变化。一个很好的例子是阿波罗任务，卫星从地球轨道开始，然后移动到月球轨道。
我已经设法检测到中央天体何时发生这种变化，并且 return 它作为积分器结果的一部分。然而，为了实际执行正确的建模，在检测到这些变化时需要更改集成期间使用的中心体。这个“换中心体”的过程涉及两个任务（注意只是坐标中心的移动，不涉及旋转）：
- 将卫星坐标减去要作为新坐标中心的天体坐标（这样，卫星坐标即为新天体坐标）
- 修改指定中心天体的参数值，该参数传递给计算加速度的函数，它是提供给定义 ODE 模型的函数的参数列表的元素之一。
我相信任务 1 可以通过使用 root-activated 事件轻松解决。为此，我在为此目的专门创建的环境中定义了一个变量，用于存储自动计算的天体的值，该天体在积分器的每次迭代中发挥主要的引力影响。在新的迭代中，将计算一个新值，并将其与先前的值进行比较。如果相同，则什么也不会发生。但如果不同，根函数将 return 0，触发事件函数。然后，事件函数将 return 位置减去新中心天体的坐标。
但是，我不确定如何执行任务 2，因为这将涉及更改提供给 ODE 模型的初始参数之一。对此有任何想法将不胜感激！ （要么是我方法的延续，要么是完全不同的方法）。
我留下了相关代码的简化版本。
我的主函数叫做hpop，是用户级函数，传递初始状态向量和其他参数。它看起来像这样：
hpop <- function(position, velocity, dateTime, times, satelliteMass, dragArea,
radiationArea, dragCoefficient, radiationCoefficient,
earthSphericalHarmonicsDegree = 130, solidEarthTides=TRUE,
oceanTides=TRUE, moonSphericalHarmonicsDegree = 150,
centralBody="Earth", ...) {
extraArgs <- list(...)
## This is the environment used to hold the variable that keeps track of what
## the central body should be at each iteration
propagationGlobalVars <- new.env()
## This is the initial value of such variable, which is the user-provided central body, by default Earth
propagationGlobalVars$latestCentralBody <- centralBody
## This is the initial state, composed by 6 variables (X, Y and Z components of position and velocity)
initial_state <- c(position, velocity)
## This is the list of parameters required for trajectory calculation
## There is quite a few, but the 2 most relevant for this question are the last 2
## centralBody is the body that will be used as the center of coordinates, and
## globalVarsEnv is the environment that will be containing the variable that keeps track of what the central body should be
parameters = list(
dateTime = dateTime,
solarArea = radiationArea,
satelliteMass = satelliteMass,
satelliteArea = dragArea,
Cr = radiationCoefficient,
Cd = dragCoefficient,
earthSPHDegree = earthSphericalHarmonicsDegree,
moonSPHDegree = moonSphericalHarmonicsDegree,
SETcorrections = solidEarthTides,
OTcorrections = oceanTides,
centralBody = centralBody,
globalVarsEnv = propagationGlobalVars)
## This calles function ode from the deSolve package, passing the previously defined initial state,
## integration times and the function defining the ode model (code below)
integration_results <- ode(y=initial_state, times=times, func=odeModel,
parms=parameters, method="radau", rtol=1e-13,
atol=1e-16, hini=0.01, ...)
numeric_results <- integration_results[, 1:7]
central_bodies <- names(centralBodiesNum[integration_results[, 8]])
output <- cbind(as.data.frame(numeric_results), central_bodies)
colnames(output) <- c("time", "X", "Y", "Z", "dX", "dY", "dZ", "Central body")
return(output)
}
ODE 模型代码的简化版本，称为 odeModel，我将其作为 func 参数传递给 deSolve 包中的函数 ode，它是：
odeModel <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
state_vector <- state
## Function accel calculates the acceleration and velocity at time t. It returns a list with
## two elements. The first is a numeric vector with the X, Y and Z components of velocity,
## and the X, Y and Z components of acceleration, in this order.
## The second element is the celestial body that, given the position at that iteration,
## exerts the main gravitational influence on the satellite
results <- accel(t, state_vector, dateTime, solarArea, satelliteMass,
satelliteArea, Cr, Cd, earthSPHDegree, SETcorrections,
OTcorrections, moonSPHDegree, centralBody)
centralBody <- results[[2]]
## Now we can compare the central body with that of the previous iteration,
## and assign the new value to the tracking variable
if(centralBody != globalVarsEnv$latestCentralBody) {
message(strwrap(paste("A transition from the sphere of influence of ",
globalVarsEnv$latestCentralBody, " to that of ",
centralBody, " has been detected.", sep=""), initial="", prefix="\n"))
assign("latestCentralBody", centralBody, envir = globalVarsEnv)
## here I would also trigger the root function, to perform Task 1 of the two
## tasks required to change the central body as described above. And also,
## I would need a way to change the value of centralBody in the parameters argument, which is the main issue of this question
}
acceleration <- results[[1]]
dx <- acceleration[1, 1]
dy <- acceleration[1, 2]
dz <- acceleration[1, 3]
d2x <- acceleration[2, 1]
d2y <- acceleration[2, 2]
d2z <- acceleration[2, 3]
## This is the return value of the odeModel function. As specified in the
## documentation for the func argument of the ode function, the first element
## of the list are the values of the derivatives of the model, and the
## second one can be any other variable. Note that since such other variables
## must be numeric, I actually access a named vector of numbers, and then convert
## back to the proper names when outputting final results to the users.
## It is important to provide the central body of each output step so that
## we know what the center of coordinates are at each step
list(c(dx, dy, dz, d2x, d2y, d2z),
centralBodiesNum[centralBody])
})
}
要根据根函数更改参数，可以使用额外的状态变量（下面的 y3），该变量在模型函数中具有导数零并且只能通过事件更改。修改教程示例 (Example3) 中的弹跳球示例，我们得到：
library(deSolve)
ball <- function(t, y, p) {
dy1 <- y[2]
dy2 <- y[3]
dy3 <- 0 # emulates a parameter, derivative = 0 but can be changed by event
list(c(dy1, dy2, dy3))
}
## gravity is essentially a parameter
yini <- c(height = 0, velocity = 10, gravity = -9.8)
rootfunc <- function(t, y, p){
return (y[1])
}
eventfunc <- function(t, y, p) {
y[1] <- 0
y[2] <- -0.9 * y[2]
y[3] <- 0.5 * y[3] # 0.5 just for technical demonstration
return(y)
}
times <- seq(from = 0, to = 20, by = 0.01)
out <- ode(times = times, y = yini, func = ball,
parms = NULL, rootfun = rootfunc,
events = list(func = eventfunc, root = TRUE))
plot(out)
这当然也适用于不止一个参数。它还可以通过强制或查找表进行扩展。
短版
根据评论中的建议，我留下了我在此处发现的问题的简短摘要版本。原文完整的解释可以在下面找到。
总而言之，我目前正在使用 deSolve 包，并且有兴趣实现响应根函数的事件。我知道我可以使用此类事件来引入模型状态变量的突然变化，但我还想修改模型的参数值以响应此类根函数。
这可能吗？
长版
我已经在 R 中实现了一个轨道数值传播器（一个计算 space 卫星位置的函数，给定初始位置和速度状态）。问题被表述为一组 6 个常微分方程 (X , 位置和速度的 Y 和 Z 分量）。我实现的模型计算任意给定时间的加速度，然后我使用 deSolve 包执行积分并计算轨迹。
进行此类计算时必须确定的一个关键参数是参考系的中心，参考系的中心通常位于对卫星施加最显着引力影响的天体的质心。这是因为，尽管原则上可以使用任意参考系进行积分和计算轨迹，但实际上只有当坐标中心位于施加主要引力影响的天体上时，我们才能获得合理的结果（即地球轨道卫星的地球，月球轨道卫星的月球等等），如讨论的那样in this SE question。
最初，我的实现使用了一个恒定的坐标中心，由用户提供或根据不同主要天体的影响范围自动确定。
但是，这不适用于行星际 space 任务的建模，因为施加主要引力影响的天体在轨迹期间会发生变化。一个很好的例子是阿波罗任务，卫星从地球轨道开始，然后移动到月球轨道。
我已经设法检测到中央天体何时发生这种变化，并且 return 它作为积分器结果的一部分。然而，为了实际执行正确的建模，在检测到这些变化时需要更改集成期间使用的中心体。这个“换中心体”的过程涉及两个任务（注意只是坐标中心的移动，不涉及旋转）：
- 将卫星坐标减去要作为新坐标中心的天体坐标（这样，卫星坐标即为新天体坐标）
- 修改指定中心天体的参数值，该参数传递给计算加速度的函数，它是提供给定义 ODE 模型的函数的参数列表的元素之一。
我相信任务 1 可以通过使用 root-activated 事件轻松解决。为此，我在为此目的专门创建的环境中定义了一个变量，用于存储自动计算的天体的值，该天体在积分器的每次迭代中发挥主要的引力影响。在新的迭代中，将计算一个新值，并将其与先前的值进行比较。如果相同，则什么也不会发生。但如果不同，根函数将 return 0，触发事件函数。然后，事件函数将 return 位置减去新中心天体的坐标。
但是，我不确定如何执行任务 2，因为这将涉及更改提供给 ODE 模型的初始参数之一。对此有任何想法将不胜感激！ （要么是我方法的延续，要么是完全不同的方法）。
我留下了相关代码的简化版本。
我的主函数叫做hpop，是用户级函数，传递初始状态向量和其他参数。它看起来像这样：
hpop <- function(position, velocity, dateTime, times, satelliteMass, dragArea,
radiationArea, dragCoefficient, radiationCoefficient,
earthSphericalHarmonicsDegree = 130, solidEarthTides=TRUE,
oceanTides=TRUE, moonSphericalHarmonicsDegree = 150,
centralBody="Earth", ...) {
extraArgs <- list(...)
## This is the environment used to hold the variable that keeps track of what
## the central body should be at each iteration
propagationGlobalVars <- new.env()
## This is the initial value of such variable, which is the user-provided central body, by default Earth
propagationGlobalVars$latestCentralBody <- centralBody
## This is the initial state, composed by 6 variables (X, Y and Z components of position and velocity)
initial_state <- c(position, velocity)
## This is the list of parameters required for trajectory calculation
## There is quite a few, but the 2 most relevant for this question are the last 2
## centralBody is the body that will be used as the center of coordinates, and
## globalVarsEnv is the environment that will be containing the variable that keeps track of what the central body should be
parameters = list(
dateTime = dateTime,
solarArea = radiationArea,
satelliteMass = satelliteMass,
satelliteArea = dragArea,
Cr = radiationCoefficient,
Cd = dragCoefficient,
earthSPHDegree = earthSphericalHarmonicsDegree,
moonSPHDegree = moonSphericalHarmonicsDegree,
SETcorrections = solidEarthTides,
OTcorrections = oceanTides,
centralBody = centralBody,
globalVarsEnv = propagationGlobalVars)
## This calles function ode from the deSolve package, passing the previously defined initial state,
## integration times and the function defining the ode model (code below)
integration_results <- ode(y=initial_state, times=times, func=odeModel,
parms=parameters, method="radau", rtol=1e-13,
atol=1e-16, hini=0.01, ...)
numeric_results <- integration_results[, 1:7]
central_bodies <- names(centralBodiesNum[integration_results[, 8]])
output <- cbind(as.data.frame(numeric_results), central_bodies)
colnames(output) <- c("time", "X", "Y", "Z", "dX", "dY", "dZ", "Central body")
return(output)
}
ODE 模型代码的简化版本，称为 odeModel，我将其作为 func 参数传递给 deSolve 包中的函数 ode，它是：
odeModel <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
state_vector <- state
## Function accel calculates the acceleration and velocity at time t. It returns a list with
## two elements. The first is a numeric vector with the X, Y and Z components of velocity,
## and the X, Y and Z components of acceleration, in this order.
## The second element is the celestial body that, given the position at that iteration,
## exerts the main gravitational influence on the satellite
results <- accel(t, state_vector, dateTime, solarArea, satelliteMass,
satelliteArea, Cr, Cd, earthSPHDegree, SETcorrections,
OTcorrections, moonSPHDegree, centralBody)
centralBody <- results[[2]]
## Now we can compare the central body with that of the previous iteration,
## and assign the new value to the tracking variable
if(centralBody != globalVarsEnv$latestCentralBody) {
message(strwrap(paste("A transition from the sphere of influence of ",
globalVarsEnv$latestCentralBody, " to that of ",
centralBody, " has been detected.", sep=""), initial="", prefix="\n"))
assign("latestCentralBody", centralBody, envir = globalVarsEnv)
## here I would also trigger the root function, to perform Task 1 of the two
## tasks required to change the central body as described above. And also,
## I would need a way to change the value of centralBody in the parameters argument, which is the main issue of this question
}
acceleration <- results[[1]]
dx <- acceleration[1, 1]
dy <- acceleration[1, 2]
dz <- acceleration[1, 3]
d2x <- acceleration[2, 1]
d2y <- acceleration[2, 2]
d2z <- acceleration[2, 3]
## This is the return value of the odeModel function. As specified in the
## documentation for the func argument of the ode function, the first element
## of the list are the values of the derivatives of the model, and the
## second one can be any other variable. Note that since such other variables
## must be numeric, I actually access a named vector of numbers, and then convert
## back to the proper names when outputting final results to the users.
## It is important to provide the central body of each output step so that
## we know what the center of coordinates are at each step
list(c(dx, dy, dz, d2x, d2y, d2z),
centralBodiesNum[centralBody])
})
}
要根据根函数更改参数，可以使用额外的状态变量（下面的 y3），该变量在模型函数中具有导数零并且只能通过事件更改。修改教程示例 (Example3) 中的弹跳球示例，我们得到：
library(deSolve)
ball <- function(t, y, p) {
dy1 <- y[2]
dy2 <- y[3]
dy3 <- 0 # emulates a parameter, derivative = 0 but can be changed by event
list(c(dy1, dy2, dy3))
}
## gravity is essentially a parameter
yini <- c(height = 0, velocity = 10, gravity = -9.8)
rootfunc <- function(t, y, p){
return (y[1])
}
eventfunc <- function(t, y, p) {
y[1] <- 0
y[2] <- -0.9 * y[2]
y[3] <- 0.5 * y[3] # 0.5 just for technical demonstration
return(y)
}
times <- seq(from = 0, to = 20, by = 0.01)
out <- ode(times = times, y = yini, func = ball,
parms = NULL, rootfun = rootfunc,
events = list(func = eventfunc, root = TRUE))
plot(out)
这当然也适用于不止一个参数。它还可以通过强制或查找表进行扩展。