Parameter manipulation while solving ODE

New to Julia
#1

Hi,

I am trying to solve an optimal control problem on the following form of Lorenz system:

function lorenz(du, u, p, t)
    theta = p
    du[1] = 10.0 * (u[2] - u[1]) + theta[1]
    du[2] = u[1] * (28.0 - (u[3])) - u[2] + theta[2]
    du[3] = u[1] * u[2] - (8/3) * (u[3]) + theta[3]
end

tspan = (0.0, T)             # time span for nonlinear problem
prob = ODEProblem(lorenz,u0,tspan,theta)                 #ODE problem definition
sol = solve(prob, RK4(),saveat = Δt, dt=Δt)          #solving the NL problem

where (theta) is the control parameter. In the previous code, theta is assumed to be constant throughout the time horizon T. However, in an optimal control problem, the control action (theta)
should be updated at each time step. Is there a way to do that without having to manually do the RK4 algorithm?

1 Like
#2

The differential equations library allows you to define parameters of any type, so functions should also be fine afaik.
Assuming your control parameter only depends on quantities that are known at each step of the integration (i.e time and the state), the following should do the trick:

function controlActionTheta(du,u,t)
   Theta = (maximum(u),t,0) #or whatever you want to define there is
   return Theta
end

function lorenz(du, u, p, t)
    theta = p(du,u,t)
    du[1] = 10.0 * (u[2] - u[1]) + theta[1]
    du[2] = u[1] * (28.0 - (u[3])) - u[2] + theta[2]
    du[3] = u[1] * u[2] - (8/3) * (u[3]) + theta[3]
end

prob = ODEProblem(lorenz,u0,tspan,controlActionTheta)

In case you want to do something different, you might also want to checkout the callback library of DifferentialEquations.jl:
https://diffeq.sciml.ai/stable/features/callback_library/

3 Likes
#3

Indeed, as @Salmon points out, the parameters can also be functions. This is used frequently, see also the Example 3 Ordinary Differential Equations · DifferentialEquations.jl.

3 Likes
#4

Thanks a lot for your prompt response!

1 Like