GitHub

Minimal Action Method using Optimal Control

The Minimal Action Method (MAM) is used to find the maximum likelihood transition paths between stable states in dynamical systems. 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:

mysqrt(x) = sqrt(x + 1e-1)
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
sol.objective
0.24942662645678104

Visualizing Results

Let's plot the solution trajectory and phase space:

plot(sol)
Example block output
# Phase space plot
MLP = sol.state.(sol.time_grid)
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=50)
println(" Time   Objective     Iterations")
for T=Ts
    global sol = solve(ocp(T); display=false, init=sol, grid_size=1000)
    @printf("%6.2f  %9.6e  %d\n", T, sol.objective, sol.iterations)
    push!(objectives, sol.objective)
end
 Time   Objective     Iterations
  1.00  4.076020e+00  2
  2.00  1.653532e+00  22
  3.00  9.192122e-01  34
  4.00  6.108613e-01  13
  5.00  4.576650e-01  5
  6.00  3.744832e-01  6
  7.00  3.269081e-01  6
  8.00  2.987674e-01  6
  9.00  2.817036e-01  6
 10.00  2.711315e-01  7
 11.00  2.644391e-01  7
 12.00  2.601049e-01  7
 13.00  2.572288e-01  8
 14.00  2.552712e-01  8
 15.00  2.539047e-01  8
 16.00  2.529271e-01  9
 17.00  2.522114e-01  10
 18.00  2.516760e-01  10
 19.00  2.512676e-01  12
 20.00  2.509506e-01  13
 21.00  2.507006e-01  12
 22.00  2.505007e-01  12
 23.00  2.503386e-01  12
 24.00  2.502058e-01  12
 25.00  2.500959e-01  11
 26.00  2.500040e-01  48
 27.00  2.499266e-01  200
 28.00  2.498608e-01  18
 29.00  2.498046e-01  15
 30.00  2.497563e-01  15
 31.00  2.497144e-01  15
 32.00  2.496780e-01  13
 33.00  2.496462e-01  13
 34.00  2.496183e-01  13
 35.00  2.495937e-01  13
 36.00  2.495719e-01  13
 37.00  2.495526e-01  13
 38.00  2.495354e-01  25
 39.00  2.495201e-01  29
 40.00  2.495064e-01  28
 41.00  2.494942e-01  17
 42.00  2.494832e-01  20
 43.00  2.494733e-01  19
 44.00  2.494644e-01  33
 45.00  2.494563e-01  22
 46.00  2.494491e-01  42
 47.00  2.494426e-01  34
 48.00  2.494367e-01  49
 49.00  2.494314e-01  20
 50.00  2.494266e-01  32
 51.00  2.494223e-01  50
 52.00  2.494185e-01  34
 53.00  2.494150e-01  66
 54.00  2.494119e-01  33
 55.00  2.494091e-01  72
 56.00  2.494066e-01  26
 57.00  2.494044e-01  23
 58.00  2.494022e-01  38
 59.00  2.494008e-01  73
 60.00  2.493993e-01  26
 61.00  2.493981e-01  31
 62.00  2.493970e-01  36
 63.00  2.493962e-01  63
 64.00  2.493955e-01  37
 65.00  2.493949e-01  52
 66.00  2.493945e-01  40
 67.00  2.493943e-01  49
 68.00  2.493942e-01  43
 69.00  2.493942e-01  78
 70.00  2.493943e-01  52
 71.00  2.493945e-01  99
 72.00  2.493949e-01  44
 73.00  2.493954e-01  54
 74.00  2.493959e-01  48
 75.00  2.493966e-01  47
 76.00  2.493973e-01  42
 77.00  2.493981e-01  35
 78.00  2.493990e-01  60
 79.00  2.494000e-01  33
 80.00  2.494011e-01  31
 81.00  2.494022e-01  33
 82.00  2.494034e-01  51
 83.00  2.494046e-01  29
 84.00  2.494060e-01  75
 85.00  2.494074e-01  37
 86.00  2.494088e-01  64
 87.00  2.494104e-01  54
 88.00  2.494119e-01  34
 89.00  2.494136e-01  57
 90.00  2.494152e-01  27
 91.00  2.494170e-01  63
 92.00  2.494188e-01  56
 93.00  2.494206e-01  66
 94.00  2.494225e-01  36
 95.00  2.494245e-01  40
 96.00  2.494264e-01  61
 97.00  2.494285e-01  71
 98.00  2.494306e-01  52
 99.00  2.494327e-01  62
100.00  2.494349e-01  78
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[20:100], log10.(objectives[20: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