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!