geometric minimal action method using optimal control

[oameye]

I would like to compute the maximum likelihood path (MLP) of an underdamped stochastic ordinary system. Suppose we SDE of the form:

x^{\prime}(t)=f(x(t))+\varepsilon \eta(t) .

with $\eta$ additive Wiener Process. The geometric minimal action method states that in the limit of $\varepsilon\ll1$ the MLP is given by optimizing the following action integral:

\hat{S}_T(x)=\int_0^1\left\{\left|x^{\prime}\right||f(x)|-(x^{\prime}, f(x))\right\} \ d s

where $\left(\cdot,\cdot\right)$ represents the dot product. Is there a way to encode this in OptimalControl.jl?

A typical example would be a double well problem such as the Maier Stein system given as:

\begin{align}
    \dot{u} &= u-u^3 - 10*u*v^2\\
    \dot{v} &= -(1+u^2)*v
\end{align}

where one wants to know the optimal path between the two attractor (-1, 0) and (1, 0).

[jbcaillau]

hi @oameye do you retain the SDE feature in what you want to solve? we only treat ODEs.

[oameye]

No effectively, the path is for the ODE, i.e., $\varepsilon=0$. But from a physics perspective you would have needed the noise. If you take the average over all the transitions you would have taken the MLP.

[oameye]

I guess my question comes down to if it is possible in OptimalControl.jl to implement the integral above.

[jbcaillau]

If it is an optimal control problem on an ODE, possibly with an integral cost, then yes. Check this example : https://control-toolbox.org/OptimalControl.jl/stable/tutorial-basic-example.html

@oameye Can you write the integral cost for the Maier Stein system above?

[oameye]

Thank you for the reaction! I have managed to get the problem working in Jump for a slightly different "problem", i.e., ~~geometric~~ minimal action method. The integral to minimise is then of the form:

\hat{S}_T(x)=\int_0^T\left\{\left|x^{\prime}-f(x)\right|^2\right\} \ d s

This disadvantage is that know one has to additional minimise $T$. The following code for jump works:

using InfiniteOpt, Ipopt, Plots

T = 50
opt = Ipopt.Optimizer   # desired solver
ns = 501;             # number of points in the time grid
model = InfiniteModel(optimizer_with_attributes(opt));

@infinite_parameter(model, t in [0, T], num_supports = ns)

xx(t) = (-1*(1-t/T) + 1*t/T)
xxx = xx.(supports(t))
yy(t) = 0.3 .* (- xx(t) .^ 2 .+ 1)
yyy = yy.(supports(t))
scatter(xxx, yyy)

@variable(model, u, Infinite(t), start = xx)
@variable(model, v, Infinite(t), start = yy)
du = u - u^3 - 10*u*v^2
dv = -(1 - u^2)*v

@objective(model, Min, ∫((∂(u, t) - du)^2 + (∂(v, t) - dv)^2, t))

@constraint(model, u(0) == -1)
@constraint(model, v(0) == 0)
@constraint(model, u(T) == 1)
@constraint(model, v(T) == 0)

# @constraint(model, ∂(u, t) == u - u^3 - 10*u*v^2)
# @constraint(model, ∂(v, t) == -(1 - u^2)*v )

optimize!(model)

u_opt = value.(u)
v_opt = value.(v)
plot(u_opt, v_opt)
xlims!(-1.1, 1.1)
ylims!(-0.1,0.5)

However, I am not so you if I can implement this in the OptimalControl.jl interface.

• Can I define the time-derivative of the control parameter (the path defined by u and v)?
• My problem doesn't have a state variable (both the variables u and v should be changed). Can I define a problem without a state variable?

[jbcaillau]

hi @oameye ; JuMP / InfiniteOpt code looks nice. As far as I can tell from what I guess from math / code you wrote, your problem looks like:

\begin{align}
& \int_0^T |u(t) - f(x(t))|^2\,\mathrm{d}t \to \min\\
& x'(t) = u(t)
\end{align}

with free final time $T$ (not a problem - I do not see $T$ optimised in what you wrote, though) and some boundary conditions on $x$ at $t=0$ and $t=T$. Correct?

[oameye]

Yeah that is correct :)

I do not see T optimised in what you wrote, though?

Yeah I didn't know how to implement a double optimisation problem, i.e.,

min_{T>0}min_{x}\int_0^T\left\{\left|x^{\prime}-f(x)\right|^2\right\} \ d s 

See technically I would have to run above InfiniteOpt code for different T.

[jbcaillau]

@oameye look like you're want to solve sth like (the math should be clear):

T = 50

@def action begin
    t ∈ [0, T], time
    x ∈ R², state
    u ∈ R², control
    x(0) == [-1, 0]
    x(T) == [1, 0]
    ẋ(t) == u(t)
    ∫( sum( (u(t) - f(x(t))).^2) ) → min
end

f(u, v) = [u - u^3 - 10*u*v^2,  -(1 - u^2)*v]
f(x) = f(x...)

Note that no init is provided, but you can have a look at the doc:

If the final time T is free, just write (but is it well posed? everything seems to be rescalable whatever T)

@def action begin
    T ∈ R, variable
    t ∈ [0, T], time
    x ∈ R², state
    u ∈ R², control
    x(0) == [-1, 0]
    x(T) == [1, 0]
    ẋ(t) == u(t)
    ∫( sum( (u(t) - f(x(t))).^2) ) → min
end

[jbcaillau]

using OptimalControl
using NLPModelsIpopt
using Plots

T = 50

@def action begin
    t ∈ [0, T], time
    x ∈ R², state
    u ∈ R², control
    x(0) == [-1, 0]
    x(T) == [1, 0]
    ẋ(t) == u(t)
    ∫( sum((u(t) - f(x(t))).^2) ) → min
end

f(u, v) = [u - u^3 - 10*u*v^2,  -(1 - u^2)*v]
f(x) = f(x...)

x1(t) = -(1 - t / T) + t / T
x2(t) = 0.3(-x1(t)^2 + 1)
x(t) = [x1(t), x2(t)]
u(t) = f(x(t))
init = (state=x, control=u)

sol = solve(action; init=init, grid_size=50)
sol = solve(action; init=sol, grid_size=1000) # grid refinement

plot(sol)

julia> sol.objective
0.24942662751291103

[oameye]

Hey @jbcaillau,

Thank you very much for figuring this out for me. It works really well.

state =  sol.state.(sol.times)
plot(first.(state), last.(state))

If the final time T is free, just write (but is it well posed? everything seems to be rescalable whatever T)

T is not free. In general it should be minimsed too, to find the maximum likehood path between to stable states in an overdamped autonomous system in the infinitesimal noise limit.

If you guys think it can be useful example as a utilty of the package, I am happy to write up an example page for it in the documentation :)

[jbcaillau]

hi @oameye

If you guys think it can be useful example as a utilty of the package, I am happy to write up an example page for it in the documentation :)

would be great, please do. we use Documenter, you can find examples in the Applications folder. See for instance https://control-toolbox.org/calculus_of_variations/dev

[ocots]

Hi @oameye

If you are still interested in making an example for the documentation, you can follow this set up tutorial : https://github.com/control-toolbox/CTApp.jl/discussions/9

Any feedback is welcome!

[oameye]

I am! I was a bit busy last weeks. I will try to write something next week 😁

[jbcaillau]

@oameye feel free to come back to us! closing the issue.