Nan values when adding constraint to inversion

Specific Domains Optimization (Mathematical)
1

Hi Julia community,

I have been trying to add a constraint to my inversion process:

using Optimization, Optim, OptimizationOptimJL

# Objective function
    function objfun(x, simulation, realsol)
        simulation.model.iceflow.A[]=x[1]*1e-17  #scaling for optimization
        
        iceflow_prob = ODEProblem{false}(Huginn.SIA2D, model.iceflow.H, params.simulation.tspan, tstops=params.solver.tstops, simulation)
        iceflow_sol = solve(iceflow_prob, 
                        params.solver.solver, 
                        callback=cb_MB, 
                        tstops=params.solver.tstops, 
                        reltol=params.solver.reltol, 
                        abstol=params.solver.reltol,
                        save_everystep=params.solver.save_everystep, 
                        progress=params.solver.progress, 
                        progress_steps=params.solver.progress_steps)

        H_obs = realsol[1]
        V_obs = realsol[2]
        
        H_pred = iceflow_sol.u[end]        
        map!(x -> ifelse(x>0.0,x,0.0), H_pred, H_pred)
        
        V_pred = Huginn.avg_surface_V(H_pred, simulation)[3] 

        ofv = 0.0
        if any((s.retcode != :Success for s in iceflow_sol))
            ofv = 1e12
        else
            mask_H = H_obs .!= 0
            mask_V = V_pred .!= 0
            
            ofv_H = mean((H_pred[mask_H] .- H_obs[mask_H]) .^ 2)
            ofv_V = mean((V_obs[mask_V] .- V_pred[mask_V]) .^ 2)

            ofv = x[2] * ofv_H + x[3] * ofv_V
            
        end
        
        println("A = $x")
        println("MSE = $ofv")
        
        return ofv
    end

    # Optimization
    fn(x, p) = objfun(x, p[1], p[2])
    cons(res, x, p) = (res .= x[2] + x[3] - 1.0)
    
    realsol = (glacier.H_glathida, glacier.V[1:end-1, 1:end-1])
    
    # Bounds for parameters
    lower_bound = [0.0085,0.0, 0.0] #A, weight_H, weight_V    
    upper_bound = [8.0, 1.0, 1.0] 
    
    optfun = OptimizationFunction(fn, Optimization.AutoForwardDiff(), cons=cons)
    optprob = OptimizationProblem(optfun, [4.0,0.5,0.5], (simulation, realsol), lb = 
    lower_bound,ub = upper_bound, lcons = [0.0], ucons = [0.0])
    sol = solve(optprob, IPNewton(), x_tol=1.0e-2, f_tol = 1e-1)

The problem is when I add the constraint (x[2]+x[3] = 1.0) as above I get that the parameters become NAN’s:

┌ Warning: First function call produced NaNs. Exiting. Double check that none of the initial conditions, parameters, or timespan values are NaN.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/initdt.jl:250
┌ Warning: Automatic dt set the starting dt as NaN, causing instability. Exiting.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:572
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase ~/.julia/packages/SciMLBase/8XHkk/src/integrator_interface.jl:574
A = [NaN, NaN, NaN]

Note that without the constraint the inversion works perfectly fine.

The only changes I made are the following:

 cons(res, x, p) = (res .= x[2] + x[3] - 1.0)
 optfun = OptimizationFunction(fn, Optimization.AutoForwardDiff(), cons=cons)
 optprob = OptimizationProblem(optfun, [4.0,0.5,0.5], (simulation, realsol), lb = 
 lower_bound,ub = upper_bound, lcons = [0.0], ucons = [0.0])
 sol = solve(optprob, IPNewton(), x_tol=1.0e-2, f_tol = 1e-1) #solver was BFGS() before

Does someone has a clue on how to solve this (or properly apply constraints for my case)?

Any help is appreciated, thanks in advance!

2

You could try adding print statements to your objfun to see where the NaN are coming from and what parameters are being evaluated?

3

Basically I have a forward model that predicts Ice Thickness and surface velocity. These can be compared to observational data for both ice thickness and velocities.

Basically the parameters that I am trying to optimize are a parameter A which is just a constant and weights for MSE calculated for velocity and MSE calculated for ice thickness.

I have printed inbetween and what happens is that for the first iteration the process starts with initial conditions (which are propagated properly) and then for the second iteration the solver proceeds with nan values for all of the parameters. Also without the constraint the optimization runs as expected. So therefore I am pretty sure that the problem lies with my implementation of the constraints. But I can’t seem to figure out what the problem is.