Geometric minimal action method

[oameye]

When I try to optimise the following system Ipopt complains that the system is not differentiable. Is the following intergral expression invalid?

using InfiniteOpt, Ipopt

T = 1
opt = Ipopt.Optimizer   # desired solver
ns = 100;             # 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) + 1*t
yy(t) = 0.3 .* (- xx(t) .^ 2 .+ 1)

@variable(model, du, Infinite(t))
@variable(model, dv, Infinite(t))
@variable(model, u, Infinite(t), start = xx)
@variable(model, v, Infinite(t), start = yy)

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

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

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

optimize!(model)

EXIT: Invalid number in NLP function or derivative detected.

[pulsipher]

This type of error occurs when the objective and/or constraints have a function or derivative that is ill-posed (e.g., divides by zero, takes the square root of a negative number).

In this case, adding appropriate guesses to du and dv should fix your problem. The default is 0 which will lead to problems in this case this the derivative of the square root objective will be ill-posed.

Regarding your formulation, what are du and dv? If these are suppose to be the derivatives of u and v then you should defined them as such.

On a side note, it is preferred that issues are reserved for bugs/feature requests and that usage questions like this one are asked in the forum (i.e., the Discussions tab).

[oameye]

Ha good catch, I will try.

On a side note, it is preferred that issues are reserved for bugs/feature requests and that usage questions like this one are asked in the forum (i.e., the Discussions tab).

My apologies, I wasn't aware. I think you can convert it to a discussion in the sidebar.

[pulsipher]

Thanks for your understanding, I have converted this to a discussion.

[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\to0$ 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 InfiniteOpt.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 attractors (-1, 0) and (1, 0).

If I correct my mistake of the derivative:

using InfiniteOpt, Ipopt, Plots

T = 1
opt = Ipopt.Optimizer   # desired solver
ns = 101;             # 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) + 1*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)

@objective(model, Min, ∫(√(∂(u, t)^2 + ∂(u, t)^2)*√(u^2 + v^2) - (∂(u, t)*u + ∂(u, t)*v), 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)

Ipopt says I have to0 few degrees of freedom:

Exception of type: TOO_FEW_DOF in file "Interfaces/IpIpoptApplication.cpp" at line 662:
 Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!

EXIT: Problem has too few degrees of freedom.

Do I have to introduce the velocity as a separate variable?

[pulsipher]

The key issue here is that with the two ODEs and the initial conditions, the problem is fully determined. Using control language, this problem has two state variables u and v, but no control variables that can be optimized (these provide degrees of freedom i.e., the knob you can turn).

Consider the steady state case where we have:
$$u - u^3 - 10 u v^2 = 0$$
$$-v(1 - u^2) = 0$$
We have two equations and two unknowns, so $u$ and $v$ are fully determined (granted these are nonlinear so there may be multiple solutions or no solutions)

Hence, adding the terminal constraints likely makes the model infeasible.

[oameye]

Thank you very much for the explanation! I have managed to get the problem working 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 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)

Do you think there is a way to also do it for the geometric minimal action method? The main problem is that I need to 2-norm (sqrt) which seems to give instabilities.

This would be a geometric implemetation:

using InfiniteOpt, Ipopt, Plots

T = 1
opt = Ipopt.Optimizer   # desired solver
ns = 101;             # 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, ∫(√((∂(v, t)^2+∂(u, t)^2)*(dv^2+dv^2))  - (∂(v, t)*dv + ∂(u, t)*du), t))
# the sqrt makes it unstable


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

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)

which give the instability error again

EXIT: Invalid number in NLP function or derivative detected.

H

[pulsipher]

I am glad to see things are coming along.

To deal with the square root, a common trick is to introduce an auxiliary variable and square it. Consider the following formulation:

using InfiniteOpt, Ipopt

model = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(model, t in [0, 1], num_supports = 101)
@variable(model, y >= -1, Infinite(t))
@objective(model, Min, integral(sqrt(y), t))

which leads to invalid derivative evaluations. We can reformulate it into:

model = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(model, t in [0, 1], num_supports = 101)
@variable(model, y >= -1, Infinite(t))
@variable(model, aux, Infinite(t)) # auxiliary variable
@objective(model, Min, integral(aux, t))
@constraint(model, aux^2 == y) # aux == sqrt(y) --> aux^2 == y

which avoids the invalid derivative evaluations.

[pulsipher]

P.S. if you wish to also minimize the final time T, then this becomes a minimum final time problem. See #289