# Nan values when adding constraint to inversion

**URL:** <https://discourse.julialang.org/t/nan-values-when-adding-constraint-to-inversion/109922>\
**Category:** Optimization (Mathematical)\
**Created:** [February 8, 2024, 11:11am UTC](https://discourse.julialang.org/t/nan-values-when-adding-constraint-to-inversion/109922 "2024-02-08T11:11:48Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![vivekag7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivekag7/32/202869_2.png) [@vivekag7](https://discourse.julialang.org/u/vivekag7)\
**Post date:** [February 8, 2024, 11:11am UTC](https://discourse.julialang.org/t/nan-values-when-adding-constraint-to-inversion/109922/1 "2024-02-08T11:11:48Z")

</div>

Hi Julia community,

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

```julia
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:

```julia
┌ 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:

```julia
 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!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 8, 2024, 10:38pm UTC](https://discourse.julialang.org/t/nan-values-when-adding-constraint-to-inversion/109922/2 "2024-02-08T22:38:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![vivekag7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vivekag7/32/202869_2.png) [@vivekag7](https://discourse.julialang.org/u/vivekag7)\
**Post date:** [February 8, 2024, 11:20pm UTC](https://discourse.julialang.org/t/nan-values-when-adding-constraint-to-inversion/109922/3 "2024-02-08T23:20:29Z")

</div>

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.