GitHub

Minimal Action Method using Optimal Control

The Minimal Action Method is a numerical technique for finding the most probable transition pathway between stable states in stochastic dynamical systems. It achieves this by minimizing an action functional that represents the path's deviation from the deterministic dynamics, effectively identifying the path of least resistance through the system's landscape. This tutorial demonstrates how to implement MAM as an optimal control problem.

Required Packages

using OptimalControl
using NLPModelsIpopt
using Plots, Printf

Problem Setup

We'll consider a 2D system with a double-well flow, called the Maier-Stein model. It is a famous benchmark problem as it exhibits non-gradient dynamics with two stable equilibrium points at (-1,0) and (1,0), connected by a non-trivial transition path. The system's deterministic dynamics are given by:

# Define the vector field
f(u, v) = [u - u^3 - 10*u*v^2,  -(1 - u^2)*v]
f(x) = f(x...)

Optimal Control Formulation

The minimal action path minimizes the deviation from the deterministic dynamics:

function ocp(T)
    action = @def begin
        t ∈ [0, T], time
        x ∈ R², state
        u ∈ R², control
        x(0) == [-1, 0]    # Starting point (left well)
        x(T) == [1, 0]     # End point (right well)
        ẋ(t) == u(t)       # Path dynamics
        ∫( sum((u(t) - f(x(t))).^2) ) → min  # Minimize deviation from deterministic flow
    end
    return action
end

Initial Guess

We provide an initial guess for the path using a simple interpolation:

# Time horizon
T = 50

# Linear interpolation for x₁
x1(t) = -(1 - t/T) + t/T

# Parabolic guess for x₂
x2(t) = 0.3(-x1(t)^2 + 1)
x(t) = [x1(t), x2(t)]
u(t) = f(x(t))

# Initial guess
init = (state=x, control=u)

Solving the Problem

We solve the problem in two steps for better accuracy:

# First solve with coarse grid
sol = solve(ocp(T); init=init, grid_size=50)

# Refine solution with finer grid
sol = solve(ocp(T); init=sol, grid_size=1000)

# Objective value
objective(sol)
0.2494266208039852

Visualizing Results

Let's plot the solution trajectory and phase space:

plot(sol)
Example block output
# Phase space plot
MLP = state(sol).(time_grid(sol))
scatter(first.(MLP), last.(MLP),
        title="Minimal Action Path",
        xlabel="u",
        ylabel="v",
        label="Transition path")
Example block output

The resulting path shows the most likely transition between the two stable states given a transient time $T=50$, minimizing the action functional while respecting the system's dynamics.

Minimize with respect to T

To find the maximum likelihood path, we also need to minimize the transient time T. Hence, we perform a discrete continuation over the parameter T by solving the optimal control problem over a continuous range of final times T, using each solution to initialize the next problem.

objectives = []
Ts = range(1,100,100)
sol = solve(ocp(Ts[1]); display=false, init=init, grid_size=200)
println(" Time   Objective     Iterations")
for T=Ts
    global sol = solve(ocp(T); display=false, init=sol, grid_size=1000, tol=1e-8)
    @printf("%6.2f  %9.6e  %d\n", T, objective(sol), iterations(sol))
    push!(objectives, objective(sol))
end
 Time   Objective     Iterations
  1.00  4.076020e+00  1
  2.00  1.653530e+00  52
  3.00  9.192103e-01  7
  4.00  6.108596e-01  6
  5.00  4.576636e-01  6
  6.00  3.744821e-01  6
  7.00  3.269073e-01  6
  8.00  2.987667e-01  6
  9.00  2.817032e-01  7
 10.00  2.711311e-01  7
 11.00  2.644388e-01  7
 12.00  2.601047e-01  8
 13.00  2.572286e-01  8
 14.00  2.552711e-01  10
 15.00  2.539046e-01  10
 16.00  2.529270e-01  12
 17.00  2.522113e-01  14
 18.00  2.516759e-01  14
 19.00  2.512676e-01  12
 20.00  2.509506e-01  12
 21.00  2.507006e-01  19
 22.00  2.505006e-01  35
 23.00  2.503386e-01  28
 24.00  2.502058e-01  30
 25.00  2.500959e-01  33
 26.00  2.500040e-01  31
 27.00  2.499266e-01  35
 28.00  2.498608e-01  39
 29.00  2.498046e-01  42
 30.00  2.497562e-01  40
 31.00  2.497144e-01  48
 32.00  2.496780e-01  45
 33.00  2.496462e-01  44
 34.00  2.496183e-01  52
 35.00  2.495937e-01  57
 36.00  2.495719e-01  50
 37.00  2.495526e-01  59
 38.00  2.495354e-01  62
 39.00  2.495201e-01  64
 40.00  2.495064e-01  68
 41.00  2.494942e-01  65
 42.00  2.494832e-01  67
 43.00  2.494733e-01  75
 44.00  2.494644e-01  72
 45.00  2.494563e-01  77
 46.00  2.494491e-01  75
 47.00  2.494426e-01  58
 48.00  2.494367e-01  80
 49.00  2.494314e-01  92
 50.00  2.494266e-01  80
 51.00  2.494223e-01  90
 52.00  2.494184e-01  87
 53.00  2.494150e-01  94
 54.00  2.494119e-01  97
 55.00  2.494091e-01  100
 56.00  2.494066e-01  73
 57.00  2.494044e-01  107
 58.00  2.494025e-01  112
 59.00  2.494008e-01  109
 60.00  2.494048e-01  12
 61.00  2.493981e-01  186
 62.00  2.493970e-01  120
 63.00  2.493962e-01  116
 64.00  2.493969e-01  18
 65.00  2.493949e-01  210
 66.00  2.493945e-01  123
 67.00  2.493943e-01  106
 68.00  2.493942e-01  60
 69.00  2.493952e-01  22
 70.00  2.493943e-01  223
 71.00  2.493945e-01  133
 72.00  2.493949e-01  153
 73.00  2.493954e-01  122
 74.00  2.493960e-01  79
 75.00  2.493966e-01  196
 76.00  2.493973e-01  155
 77.00  2.493983e-01  67
 78.00  2.493990e-01  194
 79.00  2.494000e-01  160
 80.00  2.494014e-01  62
 81.00  2.494022e-01  230
 82.00  2.494034e-01  163
 83.00  2.494046e-01  93
 84.00  2.494066e-01  21
 85.00  2.494074e-01  270
 86.00  2.494088e-01  186
 87.00  2.494104e-01  178
 88.00  2.494121e-01  50
 89.00  2.494136e-01  255
 90.00  2.494152e-01  174
 91.00  2.494170e-01  189
 92.00  2.494193e-01  31
 93.00  2.494212e-01  111
 94.00  2.494225e-01  341
 95.00  2.494249e-01  25
 96.00  2.494264e-01  370
 97.00  2.494285e-01  208
 98.00  2.494306e-01  207
 99.00  2.494331e-01  32
100.00  2.494356e-01  27
T_min = Ts[argmin(objectives)]
plt1 = scatter(Ts, log10.(objectives), xlabel="Time", label="Objective (log10)")
vline!(plt1, [T_min], label="Minimum", z_order=:back)
plt2 = scatter(Ts[40:100], log10.(objectives[40:100]), xlabel="Time", label="Objective (log10)")
vline!(plt2, [T_min], label="Minimum", z_order=:back)
plot(plt1, plt2, layout=(2,1), size=(800,800))
Example block output