# JuMP and user defined objective function: JuMPArray dimension problem

**URL:** <https://discourse.julialang.org/t/jump-and-user-defined-objective-function-jumparray-dimension-problem/17470>\
**Category:** Optimization (Mathematical)\
**Created:** [November 13, 2018, 1:05pm UTC](https://discourse.julialang.org/t/jump-and-user-defined-objective-function-jumparray-dimension-problem/17470 "2018-11-13T13:05:22Z")\
**Posts on this page:** 1\
**Showing post:** 8

<div class="post-metadata">

**Author:** ![axelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axelle/32/4488_2.png) [@axelle](https://discourse.julialang.org/u/axelle)\
**Post date:** [November 14, 2018, 10:32am UTC](https://discourse.julialang.org/t/jump-and-user-defined-objective-function-jumparray-dimension-problem/17470/8 "2018-11-14T10:32:36Z")

</div>

I cannot seem to get it working yet on a larger version of my model… I now have two problems:

1. I now get a `StackOverflowError:` with nothing more stated by Julia than that. I looked for existing posts on this topic and ounf this one [StackOverflowError - #2 by rdeits](https://discourse.julialang.org/t/stackoverflowerror/7773/2) but which does not seem to apply in my case unless I am missing something (I was able before to solve the whole model before with Optim for some values of the parameter). Why could this error arise?

2. I was getting some `NaN` expressions as my function `bk_obj` is not defined for some values of my two unknown on which the algorithm can stumble upon while searching for a solution. I tried to circumvent it in a rather wild manner by including an if condition in the function itself returning a `0` when a `Nan` is found (as I am looking for a max, i thought it could be ok…). I never solved numerically an optimisation problem even in Matlab or any other language, so I am really unsure if this is the way to go. I would be very grateful for any guidanceon this more general question too and I apologize for my ignorance…

Many thanks.

```julia
using QuadGK
using JuMP
using Clp
using Ipopt

module Para 
    Re=1.3
    Rl=1.1
    E=2.0
    points=2
    lowerB=0.01
    aversion=3
    Gamma=3.0
end

    import ..Para 
    u(c) = (c).^Para.Gamma
    
    function bk_obj(Sb, 
                    cbar,
                    D::AbstractFloat)

        Lf= Para.Rl/(Para.Re- Para.Rl) * D * cbar
        Sf= Para.E - D - Lf
        P_star= P_star= (Lf /Sb)

        θhigh = (Para.Re * (1 - Sb) + Para.Re * (1 - Sb) * (Para.Re / Para.Rl) - cbar * (P_star) * D) / (cbar * D * (Para.Re - (P_star)) )

             if isinf(θhigh) || isnan(θhigh)
            return 0.0
                else
        
        c1_ND = D * cbar
        c1_D = (1 - Sb) + Lf
        c2ND(θ) = (((1 - Sb) - θ* cbar * D) * Para.Rl + Sb * Para.Re + Lf * Para.Rl + Sf * Para.Re ) / (1 - θ)
        c2D(θ) = (1 - Sb) + Lf + ( (Sb + Sf) * Para.Re ) / (1 - θ)
        val1(θ) = u(θ * c1_ND + (1 - θ)* c2ND(θ)) 
        val2(θ) = u(θ * c1_D + (1 - θ)* c2D(θ)) 

        obj = (quadgk(val1, 0.1, θhigh)[1] + quadgk(val2, θhigh, 0.999)[1])
        
            if isinf(obj) || isnan(obj) 
                    return 0.0 
            else return obj
            end #end of if condition ensuring that integral is defined 
        end #end of if condition ensuring that threshold is defined
    end 
        
    function model_creator(D::AbstractFloat) 

        m = Model(solver = IpoptSolver())   
        @variable(m, x[1:2])

        function bk_obj_inner(x1, x2)
        return bk_obj(x1, x2, D)
        end
    
        C2ND_min(x1, x2) = ((1 - x1) - x1* D) * Para.Rl + x1 * Para.Re 
        + (Para.Rl/(Para.Re- Para.Rl) * D * x2) * Para.Rl 
        + (Para.E - D - Para.Rl/(Para.Re- Para.Rl) * D * x2) * Para.Re 
        
        JuMP.register(m, :bk_obj_inner, 2, bk_obj_inner, autodiff=true)
        JuMP.register(m, :C2ND_min, 2, C2ND_min, autodiff=true)
        @NLobjective(m, Max, bk_obj_inner(x[1], x[2]) )
        @NLconstraint(m, 0.0 <= x[1] <= Para.E )
        @NLconstraint(m, 0.001 <= x[2] <= 1.0 )
        @NLconstraint(m, 0 <=C2ND_min(x[1], x[2]))
        JuMP.solve(m)
        return getvalue(x)

    end

    model_creator(1.2)

```

and the error:

```julia

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.8, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 4
Number of nonzeros in Lagrangian Hessian.............: 0

Total number of variables............................: 2
                     variables with only lower bounds: 0
                variables with lower and upper bounds: 0
                     variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 3
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 2
        inequality constraints with only upper bounds: 1

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
   0 -0.0000000e+000 1.00e-003 6.00e-001 0.0 0.00e+000 - 0.00e+000 0.00e+000 0

StackOverflowError:

```

---

_[View the full topic](https://discourse.julialang.org/t/jump-and-user-defined-objective-function-jumparray-dimension-problem/17470)._