Initialization of model variables appears to not work in this nonlinear example

[tomdstone]

For context, I'm modeling a specific MRI sequence and trying to find optimal parameters. Basically, I can set one of the variables (called TR) fixed and then outside of JuMP do a grid search to optimize it, but when I insert it into my optimization I am unable to optimize as even though I provide inbounds parameters for initialization, the model appears to not be using some of them. I have gotten this exact same issue with both MadNLP.jl and Ipopt.jl

First, some auxiliary functions and structs (this is the relevant stuff copied over from a larger project so its a little gross, sorry):

using Pkg
Pkg.activate(temp=true)
Pkg.add(["JuMP", "MadNLP"])
# Pkg.add(["JuMP", "Ipopt"])

using JuMP
using MadNLP
# using Ipopt

const γ = 2.675e8 # rad / s / T

## Structs used

Base.@kwdef struct Scanner
    B0=3.
    B1=1e-5
    Gmax=300e-3
    Smax=730
    ADC_Δt=1e-6
end

Base.@kwdef struct ScalingParams
    propG=1.0
    propSlew=1.0
end

Base.@kwdef struct TissueParams
    T₁=0.8
    T₂=50e-3
    D=1e-10
end

Base.@kwdef struct dwSSFPParams
    α
    G
    S
    τ
    TR
end

## Functions used

b_buxton(dwp::dwSSFPParams) = (γ * dwp.G * dwp.τ)^2 * dwp.TR

β_buxton(dwp::dwSSFPParams) = (γ * dwp.G * dwp.τ)^2 * dwp.τ

function M⁻(;
    tp::TissueParams,
    dwp::dwSSFPParams,
    b = b_buxton(dwp),
    β = β_buxton(dwp),
    bD = b * tp.D,
    βD = β * tp.D,
    M₀ = 3.,
)
    A₁ = exp(-bD)
    A₂ = exp(-βD)
    A₂3rd = exp(-βD / 3) # A₂ ^ (1/3)

    cα = cosd(dwp.α)
    sα = sind(dwp.α)

    E₁ = exp(-dwp.TR / tp.T₁)
    E₂ = exp(-dwp.TR / tp.T₂)

    s = E₂ * A₁ * A₂3rd^(-4) * (1 - E₁ * cα) + E₂ * A₂3rd^(-1) * (cα - 1)
    r = 1 - E₁ * cα + E₂^2 * A₁ * A₂3rd * (cα - E₁)
    K = (1 - E₁ * A₁ * cα - E₂^2 * A₁^2 * A₂3rd^(-2) * (E₁ * A₁ - cα)) / (E₂ * A₁ * A₂3rd^(-4) * (1 + cα) * (1 - E₁ * A₁))
    F₁ = K - sqrt(K^2 - A₂^2)

    return - M₀ * (1 - E₁) * E₂ * A₂3rd^(-2) * (F₁ - E₂ * A₁ * A₂3rd^2) * sα / (r - F₁ * s)
end

function b_eff(dwp_diff::dwSSFPParams; tp::TissueParams, M₀::Real = 3.)
    dwp_non  = dwSSFPParams(α = dwp_diff.α, G = 0., S = 0., τ = 0., TR = dwp_diff.TR)

    return -1 / tp.D * (log(M⁻(; tp, dwp = dwp_diff, M₀)) - log(M⁻(; tp, dwp = dwp_non, M₀)))
end

function η(S, dwp::dwSSFPParams, readout_time)
    return S * sqrt(readout_time / dwp.TR)
end

Now, I set up the model

b = 1e9
readout_time = 40e-3
sp = ScalingParams()
scanner = Scanner()
tp = TissueParams()
rf_time = 6e-3
init_ = (α = 10., τ = 2e-3, S = 0.95 * (scanner.Smax * sp.propSlew), ξ = 1e-4, TR = 60e-3)

dead_time = 2e-3

model = Model(MadNLP.Optimizer)

@variable(model, 1. <= α <= 179.)
@variable(model, 0. <= τ)
@variable(model, 0. <= S <= scanner.Smax * sp.propSlew)
@variable(model, 0. <= ξ)
@variable(model, 0. <= TR)

@constraint(model, model[:ξ] <= model[:τ])
@constraint(model, model[:ξ] * model[:S] <= scanner.Gmax * sp.propG)
@constraint(model, model[:τ] + model[:ξ] + rf_time + readout_time + dead_time <= model[:TR])
@constraint(model, b == b_eff(
    dwSSFPParams(model[:α], model[:ξ] * model[:S], model[:S], model[:τ], model[:TR]);
    tp = tp,
    M₀ = scanner.B0
))

set_start_value(model[:α], init_.α)
set_start_value(model[:τ], init_.τ)
set_start_value(model[:S], init_.S)
set_start_value(model[:ξ], min(0.95 * scanner.Gmax * sp.propG / init_.S, init_.ξ))
set_start_value(model[:TR], init_.TR)

@variable(model, 0 <= signal)
@constraint(model, model[:signal] == M⁻(;
    tp = tp,
    dwp = dwSSFPParams(model[:α], model[:ξ] * model[:S], model[:S], model[:τ], model[:TR]),
    M₀ = scanner.B0
))

@objective(
    model,
    Max,
    η(
        model[:signal],
        dwSSFPParams(model[:α], model[:ξ] * model[:S], model[:S], model[:τ], model[:TR]),
        readout_time
    )
)

# checking that starting values are good

dwp_init = dwSSFPParams(
    model[:α] |> start_value,
    (model[:ξ] |> start_value) * (model[:S] |> start_value),
    model[:S] |> start_value,
    model[:τ] |> start_value,
    model[:TR] |> start_value
)

signal_init = M⁻(;
    tp = tp,
    dwp = dwp_init,
    M₀ = scanner.B0
)

η(
    signal_init,
    dwp_init,
    readout_time
)

b_eff(
    dwp_init;
    tp = tp,
    M₀ = scanner.B0
) # starts near, but not at, correct b value

actually optimizing

optimize!(model)

results = (
    α = value(model[:α]),
    τ = value(model[:τ]),
    ξ = value(model[:ξ]),
    G = value(model[:ξ]) * value(model[:S]),
    S = value(model[:S]),
    TR = value(model[:TR]),
    η = objective_value(model),
    signal = value(model[:signal]),
    b = value(b_eff(
        dwSSFPParams(model[:α], model[:ξ] * model[:S], model[:S], model[:τ], model[:TR]);
        tp = tp,
        M₀ = scanner.B0
    )),
    status = termination_status(model)
)

note that S * ξ > scanner.Gmax * sp.propG, is significantly too large, which makes K=Inf and thus F₁ becomes NaN in the M⁻ equation.
Also note that ξ, τ, and S are not the values I initially set them to, while α and TR are.

# the values that it has result in NaNs, thus optimization fails

dwp_res = dwSSFPParams(results.α, results.G, results.S, results.τ, results.TR)

M⁻(;
    tp = tp,
    dwp = dwp_res,
    M₀ = scanner.B0
)

b_eff(
    dwp_res;
    tp = tp,
    M₀ = scanner.B0
)

η(
    M⁻(;
        tp = tp,
        dwp = dwp_res,
        M₀ = scanner.B0
    ),
    dwp_res,
    readout_time
)

But when I fix the TR, we get the following:

TR = 60e-3
b = 1e9
scanner = Scanner()
tp = TissueParams()
rf_time = 6e-3
init_ = (α = 10., τ = 2e-3, S = 100., ξ = 1e-4)

model2 = Model(MadNLP.Optimizer)
# model2 = Model(Ipopt.Optimizer)

@variable(model2, 1. <= α <= 179.)
@variable(model2, 0. <= τ <= TR) # needed the <= TR/2 for convergence stability
@variable(model2, 0. <= S <= scanner.Smax)
@variable(model2, 0. <= ξ <= TR/2) # needed the <= TR/2 for convergence stability

@constraint(model2, model2[:ξ] <= model2[:τ])
@constraint(model2, model2[:τ] + model2[:ξ] + rf_time <= TR)
@constraint(model2, model2[:ξ] * model2[:S] <= scanner.Gmax)

set_start_value(model2[:S], init_.S)
set_start_value(model2[:α], init_.α)
set_start_value(model2[:τ], init_.τ)
set_start_value(model2[:ξ], init_.ξ)

@constraint(model2, b == b_eff(
    dwSSFPParams(model2[:α], model2[:ξ] * model2[:S], model2[:S], model2[:τ], TR);
    tp = tp,
    M₀ = scanner.B0
))

@variable(model2, 0 <= signal)
@constraint(model2, model2[:signal] == M⁻(;
    tp = tp,
    dwp = dwSSFPParams(model2[:α], model2[:ξ] * model2[:S], model2[:S], model2[:τ], TR),
    M₀ = scanner.B0
))

@objective(
    model2,
    Max,
    η(
        model2[:signal],
        dwSSFPParams(model2[:α], model2[:ξ] * model2[:S], model2[:S], model2[:τ], TR),
        TR - rf_time - model2[:τ] - model2[:ξ]
    )
)

dwp_init2 = dwSSFPParams(
    model2[:α] |> start_value,
    (model2[:ξ] |> start_value) * (model2[:S] |> start_value),
    model2[:S] |> start_value,
    model2[:τ] |> start_value,
    TR
)

signal_init2 = M⁻(;
    tp = tp,
    dwp = dwp_init2,
    M₀ = scanner.B0
)

η(signal_init2, dwp_init2, readout_time)

b_eff(dwp_init2; tp = tp,M₀ = scanner.B0)

optimize!(model2)

results2 = (
    α = value(model2[:α]),
    τ = value(model2[:τ]),
    ξ = value(model2[:ξ]),
    G = value(model2[:ξ]) * value(model2[:S]),
    S = value(model2[:S]),
    η = objective_value(model2),
    signal = value(model2[:signal]),
    status = termination_status(model2)
)

All these numbers are reasonable:

dwp_res2 = dwSSFPParams(results2.α, results2.G, results2.S, results2.τ, TR)

b_eff(dwp_res2; tp, M₀=scanner.B0)

signal_res2 = M⁻(;
    tp = tp,
    dwp = dwp_res2,
    M₀ = scanner.B0
)

η(
    signal_res2,
    dwp_res2,
    TR - (results2.τ + results2.ξ + rf_time)
)

Inspecting the results tuples, it seems like the first optimizer starts at a place that is not in the domain even though I gave it a spot in the domain. What could be going on here?

Mac M4 Max, Julia version 1.12.6

[odow]

Inspecting the results tuples, it seems like the first optimizer starts at a place that is not in the domain even though I gave it a spot in the domain. What could be going on here?

Your first solve failed:

julia> optimize!(model)
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.

Number of nonzeros in equality constraint Jacobian...:       11
Number of nonzeros in inequality constraint Jacobian.:        7
Number of nonzeros in Lagrangian Hessian.............:       34

The Jacobian for the equality constraints contains an invalid number

Number of Iterations....: 0

Number of objective function evaluations             = 0
Number of objective gradient evaluations             = 0
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 1
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 0.000

EXIT: Invalid number in NLP function or derivative detected.

It is not valid to query value of a variable without checking what the primal status is. In your case:

julia> primal_status(model)
INFEASIBLE_POINT::ResultStatusCode = 3

There is no guarantee that the returned value will be the initial starting point.

[tomdstone]

Is there a way to determine what the most recent failure point was? If so, then I can check if perhaps the derivative at that point is the issue. A cursory google search suggests that what is returned by value is the most recent point the solver was working with. So might it be possible to infer that, given no objective or gradient evaluations were performed, no steps had been taken yet and what is returned by value is in fact the initial point it evaluated the function at? It seems extremely unlikely that the value 0.001 would exactly come up for multiple variables if it had actually taken a step using the gradient (especially considering the same thing happened with multiple solvers).

[odow]

Is there a way to determine what the most recent failure point was?

Not easily, because the error happens before Ipopt logs the first iterate (so you can't use a callback https://github.com/jump-dev/Ipopt.jl#solver-specific-callbacks).

A cursory google search suggests that what is returned by value is the most recent point the solver was working with

This is not true. We provide no guarantee on what solution will be returned by the solver. If the solver fails, it might return the starting point, a current feasible point, no solution, or a completely garbled solution that is an infeasible point.

So might it be possible to infer that, given no objective or gradient evaluations were performed, no steps had been taken yet and what is returned by value is in fact the initial point it evaluated the function at? It seems extremely unlikely that the value 0.001 would exactly come up for multiple variables if it had actually taken a step using the gradient (especially considering the same thing happened with multiple solvers).

Ipopt has a number of options that perturb the problem data, including variable bounds and starting values. One example is https://coin-or.github.io/Ipopt/OPTIONS.html#OPT_bound_push

If you set:

set_attribute(model, "bound_push", 1e-5)
optimize!(model)

Ipopt can successfully start a few iterates, but it still gets stuck somewhere and gives up.

In general, I don't find it useful to analyse the particular details of solver iterates. Spend more time writing a nice model. For non-convex problems, good starting points are essential.

[odow]

We've previously discussed building tools to debug nonlinear problems: #3664

[tomdstone]

Thanks! That's very helpful.

[tomdstone]

I'm marking it as closed because it seems that the intended behavior is to perturb the problem data, and I should really rewrite my problem not using SI units as the scales vary wildly across variables.