JuMP and user defined objective function: JuMPArray dimension problem

Specific Domains Optimization (Mathematical)
8

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 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.

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:


******************************************************************************
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: