Question about applying a customized branch and bound in Julia and the JuMP package

I am trying to develop a customised branch and bound algorithm in Julia with the JuMP package. The problem only has binary variables.

What I do is relaxing the problem to its linear counterpart and getting the variable values from the linear model. I define a loop to detect the fractional values. At every fractional solution the following should be conducted.
If x is a fractional variable, create two new branches (two new models) one with x >= 1 and one with x =< 0 constraints. I have error with this one and I put it below.
The created models should be stored somewhere (in a container or vector) with appropriate names so I can call them and if they violate the lower bound or upper bound they can be deleted. I’m not sure how to do this.

mp_linear = Model(Gurobi.Optimizer)
@variable(mp_linear, 0 ≤ ψ[1:H] ≤ 1 )
@variable(mp_linear, 0 ≤ χ[1:H , 1:J] ≤ 1  )
@objective(mp_linear, Min, sum( T*zᵖ[h]*ψ[h] for h in 1:H ) + sum(χ[h,j] * zᵐ[h,j] for h in 1:H for j in 1:J) )
@constraint(mp_linear, [h in 1:H , j in 1:J], χ[h,j] ≤ ψ[h])
@constraint(mp_linear, [h in 1:H , j in 1:J], s[h,j] + ψ[h] ≤ χ[h,j] + 1)
@constraint(mp_linear, 335*χ[1,1]+335*χ[2,1] + 335* χ[3,1] + 335* χ[1,2] + 335 *χ[2,2] + 
    335* χ[3,2] + 335* χ[1,3] + 335* χ[2,3] + 335 *χ[3,3]   >= 380.0)

function initiate()
        χᵃ = value.(χ)
        @show χᵃ
        for h in 1:H
            for j in 1:J
                if 0 < χᵃ[h,j] < 1
                    global m1 = copy(mp_linear)
                    @constraint(m1, χ[h,j] ≥ 1)
                    global m0 = copy(mp_linear)
                    @constraint(m0, χ[h,j] ≤ 1)


and this is the error I see

Solved in 7 iterations and 0.00 seconds (0.00 work units)
Optimal objective  2.268656716e+01

User-callback calls 66, time in user-callback 0.00 sec
χᵃ = [0.0 0.0 0.0; 0.0 0.0 0.0; 0.3781094527363184 0.37810945273631835 0.37810945273631835]

Can you please help with this?
Is there a proper example somewhere? I am sure I will have heaaches for the rest of the algorithm lol.

Can you please help with this?

χ belongs to mp_linear so you cannot use it in m1. You need to find the corresponding variable in m1. Moreover, you should use copy_model instead of copy.
So, here is what you should do:

m1, r1 = copy_model(mp_linear)
χ1 = getindex.(r1, χ)
# or
χ1 = m1[:χ]
@constraint(m1, χ1[h,j] ≥ 1)

Is there a proper example somewhere?

GitHub - lanl-ansi/Juniper.jl: A JuMP-based Nonlinear Integer Program Solver implements a Branch & Bound that even supports MINLP (i.e., arbitrary nonlinear constraints). Might be not the most easy to read since it is quite general.

Why? What’s the larger context?

The created models should be stored somewhere

Creating a new model for each branch is the wrong way to do this (the overhead cost is too large). You probably just want to be modifying the bounds of the variables with set_lower_bound and set_upper_bound.

Then, instead of storing the models you can just store the vectors of bounds.

See also GitHub - Wikunia/Bonobo.jl: A general branch and bound framework