View source: R/derivedEquations.R
adjointSymb | R Documentation |
Compute adjoint equations of a function symbolically
adjointSymb(f, states = names(f), parameters = NULL, inputs = NULL)
f |
Named vector of type character, the functions |
states |
Character vector of the ODE states for which observations are available |
parameters |
Character vector of the parameters |
inputs |
Character vector of the "variable" input states, i.e. time-dependent parameters (in contrast to the forcings). |
The adjoint equations are computed with respect to the functional
(x, u) -> int( ||x(t) - xD(t)||^2 + ||u(t) - uD(t)||^2, dt),
where x are the states being constrained
by the ODE, u are the inputs and xD and uD indicate the trajectories to be best
possibly approached. When the ODE is linear with respect to u, the attribute inputs
of the returned equations can be used to replace all occurences of u by the corresponding
character in the attribute. This guarantees that the input course is optimal with
respect to the above function.
Named vector of type character with the adjoint equations. The vector has attributes "chi" (integrand of the chisquare functional), "grad" (integrand of the gradient of the chisquare functional), "forcings" (character vector of the forcings necessary for integration of the adjoint equations) and "inputs" (the input expressed as a function of the adjoint variables).
## Not run: ###################################################################### ## Solve an optimal control problem: ###################################################################### library(bvpSolve) # O2 + O <-> O3 # O3 is removed by a variable rate u(t) f <- c( O3 = " build_O3 * O2 * O - decay_O3 * O3 - u * O3", O2 = "-build_O3 * O2 * O + decay_O3 * O3", O = "-build_O3 * O2 * O + decay_O3 * O3" ) # Compute adjoints equations and replace u by optimal input f_a <- adjointSymb(f, states = c("O3"), inputs = "u") inputs <- attr(f_a, "inputs") f_tot <- replaceSymbols("u", inputs, c(f, f_a)) forcings <- attr(f_a, "forcings") # Initialize times, states, parameters times <- seq(0, 15, by = .1) boundary <- data.frame( name = c("O3", "O2", "O", "adjO3", "adjO2", "adjO"), yini = c(0.5, 2, 2.5, NA, NA, NA), yend = c(NA, NA, NA, 0, 0, 0)) pars <- c(build_O3 = .2, decay_O3 = .1, eps = 1) # Generate ODE function func <- funC(f = f_tot, forcings = forcings, jacobian = "full", boundary = boundary, modelname = "example5") # Initialize forcings (the objective) forcData <- data.frame(time = times, name = rep(forcings, each=length(times)), value = rep( c(0.5, 0, 1, 1), each=length(times))) forc <- setForcings(func, forcData) # Solve BVP out <- bvptwpC(x = times, func = func, parms = pars, forcings = forc) # Plot solution par(mfcol=c(1,2)) t <- out[,1] M1 <- out[,2:4] M2 <- with(list(uD = 0, O3 = out[,2], adjO3 = out[,5], eps = 1, weightuD = 1), eval(parse(text=inputs))) matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states") abline(h = .5, lty=2) legend("topright", legend = names(f), lty=1, col=1:3) matplot(t, M2, type="l", lty=1, col=1, xlab="time", ylab="value", main="input u") abline(h = 0, lty=2) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.