infiniteopt / InfiniteOpt.jl

An intuitive modeling interface for infinite-dimensional optimization problems.
https://infiniteopt.github.io/InfiniteOpt.jl/stable
MIT License
251 stars 17 forks source link

Geometric minimal action method #353

Closed oameye closed 1 month ago

oameye commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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.

oameye commented 1 month ago

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). image

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?