In [1]:
using LightGraphs
using TikzGraphs, TikzPictures
using PyCall
include("MIGP.jl")
Out[1]:
solve_MIGP (generic function with 2 methods)

Visualize the branch and bound algorithm on a simple optimization problem

$min$ $Z$

$s.t.$ $Z \geq x_1 + x_2$

$6x_2 - 2x_3 \geq 10$

$4x_1 \geq 3x_3$

$2x_1 \geq x_2$

$x_i \geq 0, integer$ where $i\in$ {1,2,3}

This notebook walks through the implementation of the branch and bound algorithm to add integer constraints to GP-compatible optimization problems. This is a useful extension of the GP optimization framework, and has applications in aerospace system design problems.

The cell below is python code, but shows the declaration of the simple problem described above. The INT_ before the variable name is the declaration that the variable is constrained to be integer; the Julia code below will parse the python model and constrain anv variables where the names start with that flag to be integer.

Note that because geometric programming constraints are only convex under a logarithmic change of variables, zero is not a valid integer value, since log(0) = -$\infty$. Therefore, zero is replaced with a small positive number (10^-30), which has the same practical effect.

For more on how geometric programming constraints are solved, see https://docs.mosek.com/modeling-cookbook/expo.html#geometric-programming.

In [2]:
from gpkit import Model, parse_variables, Vectorize, SignomialsEnabled, Variable, SignomialEquality, units

class Test1(Model):
    """ Test equation

    Variables
    ---------
    INT_x1                  [-]     x1
    INT_x2                  [-]     x2
    INT_x3                  [-]     x3
    Z                   [-]     Z

    """
    def setup(self):
        exec parse_variables(Test1.__doc__)
        ZERO = 1e-30
        x1 = INT_x1
        x2 = INT_x2
        x3 = INT_x3
        constraints = [Z >= x1+x2,
            6*x2 >= 10+ 2*x3,
            4*x1 >= 3*x3,
            2*x1 >= x2,
            x1 >= ZERO,
            x2 >= ZERO,
            x3 >= ZERO

        ]

        self.cost = Z

        return constraints
syntax: extra token "gpkit" after end of expression

The code above is saved in the file mico_test.py, which is imported below.

In [3]:
PyCall.pyversion
@pyimport mico_tests as mt

Lightgraphs and Tikzplots will be used to visualize the running of the algorithm; some of the graph generation setup code is shown below.

In [4]:
G = Graph()
Out[4]:
{0, 0} undirected simple Int64 graph
In [5]:
mutable struct CostedNode
    vars::Array{String}
    vals::Array{Float64}
    vars_expanded::Array{String}
    best_vals::Array{Float64}
    cost::Float64
    level::Int16
    ID::Int16
    parentID::Int16
    label::String

end
In [6]:
function addNode!(Node, G, Node_List)
    add_vertex!(G)
    push!(Node_List,Node)
end
Out[6]:
addNode! (generic function with 1 method)

Solve Relaxed Problem

The algorithm starts by solving the relaxed problem; the optimal non-integer solution of all variables that satisfies the constraints in computed. This is done using an interior point method solver, which solves very quickly for this small problem. However, one of the advantages of this method is that even for large problems with thousands of free variables solution times are still on the order of single seconds.

In [7]:
Node_List = CostedNode[]
var_names, values,cost  = mt.run_Test1(verbosity = 2)
N0 = CostedNode(var_names, values,[],[], cost, 1, 1, 0, "R")
addNode!(N0, G, Node_List)

plot(G,["R"],node_style="draw, rounded corners, fill=blue!10", options="scale=2, font=\\huge\\sf")
Using solver 'mosek'
Solving for 4 variables.
Number of Hessian non-zeros: 7
* Solving exponential optimization problem on dual form. *
* The following log information refers to the solution of the dual problem. *
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : GECO (general convex optimization problem)
  Constraints            : 5               
  Cones                  : 0               
  Scalar variables       : 10              
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Matrix reordering started.
Local matrix reordering started.
Local matrix reordering terminated.
Matrix reordering terminated.
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : GECO (general convex optimization problem)
  Constraints            : 5               
  Cones                  : 0               
  Scalar variables       : 10              
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer  - threads                : 4               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 3
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 7                 conic                  : 0               
Optimizer  - Semi-definite variables: 0                 scalarized             : 0               
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 22                after factor           : 27              
Factor     - dense dim.             : 0                 flops                  : 5.20e+01        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.7e+00  1.2e+02  2.1e+02  0.00e+00   -2.057410035e+02  0.000000000e+00   1.0e+00  0.00  
1   1.5e-01  1.0e+01  1.8e+01  -1.00e+00  -1.723410155e+02  2.366551184e+01   8.8e-02  0.00  
2   2.6e-02  1.7e+00  3.1e+00  -8.19e-01  -2.471272075e+01  6.251800046e+01   1.5e-02  0.00  
3   7.2e-03  4.6e-01  8.8e-01  1.39e+00   -1.598704296e+01  6.965627889e+00   4.3e-03  0.00  
4   1.0e-03  7.3e-02  1.4e-01  9.87e-01   -7.959873841e-01  2.435784190e+00   7.8e-04  0.00  
5   8.3e-05  5.3e-02  1.3e-02  9.84e-01   7.313874228e-01   1.022475284e+00   7.5e-05  0.00  
6   2.8e-06  4.6e-02  1.2e-03  9.74e-01   8.976056091e-01   9.252794099e-01   9.8e-06  0.00  
7   1.4e-07  2.3e-03  7.0e-05  1.01e+00   9.152537750e-01   9.167936709e-01   5.3e-07  0.00  
8   5.3e-10  8.5e-06  4.7e-07  1.00e+00   9.162846690e-01   9.162953032e-01   2.9e-09  0.00  
9   1.8e-12  2.9e-08  2.9e-09  1.00e+00   9.162906978e-01   9.162907654e-01   1.6e-11  0.00  
Optimizer terminated. Time: 0.00    

solsta = 1, prosta = 1

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 9.1629069779e-01    nrm: 1e+00    Viol.  con: 3e-11    var: 0e+00  
  Dual.    obj: 9.1629076537e-01    nrm: 7e+01    Viol.  con: 0e+00    var: 8e-07  
* End solution on dual form. *
Transforming to primal solution.
Solving took 0.011 seconds.
result packing took 3.2% of solve time
solution checking took 1.5e+02% of solve time
Out[7]:

This solution is the base node of the tree which will be searched to determine the optimal solution.

In [8]:
dump(N0)
CostedNode
  vars: Array{String}((4,))
    1: String "Z_Test1"
    2: String "INT_x1_Test1"
    3: String "INT_x3_Test1"
    4: String "INT_x2_Test1"
  vals: Array{Float64}((4,)) [2.5, 0.833333, 1.0e-30, 1.66667]
  vars_expanded: Array{String}((0,))
  best_vals: Array{Float64}((0,)) Float64[]
  cost: Float64 2.4999999147791763
  level: Int16 1
  ID: Int16 1
  parentID: Int16 0
  label: String "R"

Define some functions

These functions are supporting funtions to implement the branch and bound algorithm.

In [9]:
function identify_INT_vars(vars, vals)
    
    #count int variables
    c = 0
    int_vars = String[]
    int_vals = Float64[]
    
    for (i,var) in enumerate(vars)
        key = split(var,"_")[1]
        if key == "INT"
            push!(int_vars,var)
            push!(int_vals,vals[i])
        end
    end
        return int_vars,int_vals
    end
Out[9]:
identify_INT_vars (generic function with 1 method)
In [10]:
function nearest_ints(val)
    rounded = round(val)
    
    if rounded >= val
        return (rounded, rounded-1)
    else
        return (rounded+1,rounded)
    end
end
Out[10]:
nearest_ints (generic function with 1 method)
In [11]:
function isInt(i)
   tol = 1e-3
    return min(i - floor(i), ceil(i)-i) < tol
end
Out[11]:
isInt (generic function with 1 method)
In [12]:
function cleanup_input(varname::Array{String})
    output = String[]
    
    for str in varname
        trim = split(str, '.')[1]
        elem = split(trim, '_')
        s = ""
        for e in elem
            if e == "Test1"
                s = s[1:end-1]
                break
            else
                s = s*String(e)*'_'
            end
        end
        push!(output,s)
    end
    return output
end
Out[12]:
cleanup_input (generic function with 2 methods)
In [13]:
function constrain_variable(Parent::CostedNode, G::Graph, constrained_vars, constrained_vals,g, 
        var2expand::String, current_best)
    
    if isInt(current_best)
        #findall(x->x==2, A)
        j = findall(x->x== var2expand,cleanup_input(Parent.vars))
        child_vars = copy(Parent.vars)
        child_vals = copy(Parent.vals)
        
        deleteat!(child_vars, j)
        deleteat!(child_vals, j)
        N = length(G.fadjlist)
        add_vertex!(G)
        l = var2expand*"="*string(round(current_best))
        best_child = CostedNode(child_vars, child_vals, 
            [Parent.vars_expanded...,var2expand],[Parent.best_vals...,current_best],
            Parent.cost,Parent.level+1, N+1, Parent.ID, l)
        
        add_edge!(G, Parent.ID, best_child.ID)
        worst_child = nothing
        
    else    
        (above, below) = nearest_ints(current_best)
        expanded = [Parent.vars_expanded..., var2expand]
        
        best_above = [Parent.best_vals..., above]
        best_below = [Parent.best_vals..., below]
        
        (vars_above, vals_above, cost_above) = g(expanded, best_above)
        (vars_below, vals_below, cost_below) = g(expanded, best_below)
        
        N = length(G.fadjlist)
        
        labove = var2expand*"="*string(above)
        Nabove = CostedNode(vars_above, vals_above,
            expanded,best_above,
            cost_above, Parent.level+1,N+1, Parent.ID, labove)
        
        lbelow = var2expand*"="*string(below)
        Nbelow = CostedNode(vars_above, vals_above,
            expanded,best_below,
            cost_above, Parent.level+1,N+2, Parent.ID, lbelow)
        
        add_vertex!(G)
        add_vertex!(G)
        add_edge!(G, Parent.ID, Nabove.ID)
        add_edge!(G, Parent.ID, Nbelow.ID)

        if cost_above == nothing && cost_below == nothing
            println("No feasible integer solution")
        elseif cost_above == nothing
            println("Likely infeasible above")
            best_val = below
            best_cost = cost_below[1]
            remaining_vars, remaining_values = identify_INT_vars(vars_below,vals_below)
            best_child = Nbelow
            worst_child = Nabove
        elseif cost_below == nothing
            println("Likely infeasible below")
            best_val = above
            best_cost = cost_above[1]
            remaining_vars, remaining_values = identify_INT_vars(vars_above,vals_above)
            best_child = Nabove
            worst_child = Nbelow

        elseif cost_above[1] <= cost_below[1]
            best_val = above
            best_cost = cost_above[1]
            remaining_vars, remaining_values = identify_INT_vars(vars_above,vals_above)
            best_child = Nabove
            worst_child = Nbelow

        else
            best_val = below
            best_cost = cost_below[1]
            remaining_vars, remaining_values = identify_INT_vars(vars_below,vals_below)
            best_child = Nbelow
            worst_child = Nabove
        end
    end
    
    return (best_child, worst_child)
    
end
Out[13]:
constrain_variable (generic function with 1 method)
In [14]:
Int_Vars, Int_Vals = identify_INT_vars(cleanup_input(N0.vars),N0.vals)
cols = fill("fill=red!10",2*length(Int_Vars)+1)
cols[1] = "fill=blue!10"

function expand!(ParentNode,Graph, Int_Vars, Int_Vals,g)
    
    current_rank = ParentNode.level    #Identify current level of the tree (which keys the variable to expand)
    child_rank = current_rank+1 #Put all children one level below
    Nvars = length(Int_Vars)
    #i = current_rank  #Determine the variable which gets expanded based on the rank
    var2expand   = Int_Vars[1]
    val2expand   = Int_Vals[1]
    
    (best_child, worst_child) = constrain_variable(ParentNode, Graph, Int_Vars, Int_Vals,g, 
        var2expand, val2expand)
    
end
Out[14]:
expand! (generic function with 1 method)
In [15]:
function drawgraph(G, Node_List, cols)
    labels = fill("", length(Node_List))
    node_styles=Dict()
    for Node in Node_List
        e = split(Node.label, '_')
        labels[Node.ID] = e[min(2,length(e))]
    end
    for (i,col) in enumerate(cols[1:length(Node_List)])
        node_styles[i] =col
    end
    plot(G, labels,node_style="draw, rounded corners, fill=blue!10", node_styles=node_styles,
        options="scale=2, font=\\huge\\sf")
end
Out[15]:
drawgraph (generic function with 1 method)

Constrain the first integer variable

The second step in the algorithm is to pick on of the integer variables, and add two new nodes to the tree; one where the variable is assigned the closest integer value greater than the relaxed solution, and once which is the closest integer value smaller than the relaxed solution. The problem is re-solved for each node, and whichever has the higher cost (or is infeasible) is pruned.

In [16]:
g(vars, vals) = mt.run_Test1(vars, vals)
Out[16]:
g (generic function with 1 method)
In [17]:
(Parent, Pruned) = expand!(N0, G, Int_Vars, Int_Vals, g)
push!(Node_List, Parent)
push!(Node_List, Parent)
Node_List[Pruned.ID] = Pruned
cols[Parent.ID] = "fill=blue!10"
Likely infeasible below
Likely infeasible model at ['INT_x1']
 [0.]
Out[17]:
"fill=blue!10"
In [18]:
drawgraph(G, Node_List, cols)
Out[18]:

Expand subsequent nodes

The lowest-cost node is then expanded, and the process repeated for all remaining integer variables

In [69]:
Int_Vars, Int_Vals = identify_INT_vars(cleanup_input(Parent.vars),Parent.vals)
(Parent, Pruned) = expand!(Parent, G, Int_Vars, Int_Vals, g)
if Parent != nothing
    cols[Parent.ID] = "fill=blue!10"
    push!(Node_List, Parent)
else
    println("No valid integer solution")
end
    
if Pruned != nothing
    push!(Node_List, Parent)
    Node_List[Pruned.ID] = Pruned
end
drawgraph(G, Node_List, cols)
Out[69]:

Compare Timing With JuMP

There are many readily available solvers that can easily handle this simple problem. JuMP provides access to many of these, and a convenient formulation.

Note: Mosek is not required to run this model; any solver compatible with JuMP (many of which are open source) can be used. https://www.juliaopt.org/

In [20]:
using JuMP
using Mosek
In [26]:
m = Model(solver = MosekSolver())
@variable(m, 0 <= x1 <= Inf,Int)
@variable(m, 0 <= x2 <= Inf,Int)
@variable(m, 0 <= x3 <= Inf,Int)
@variable(m, 0 <= Z <= Inf)

@objective(m, Min, Z )
@constraint(m, x1+x2 <= Z )
@constraint(m, 10 + 2x3 <= 6x2 )
@constraint(m, 4x1 >= 3x3 )
@constraint(m, 2x1 >= x2 )

print(m)

JuMP_Time = @elapsed status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x1 = ", getvalue(x1))
println("x2 = ", getvalue(x2))
println("x3 = ", getvalue(x3))
Min Z
Subject to
 x1 + x2 - Z ≤ 0
 2 x3 - 6 x2 ≤ -10
 4 x1 - 3 x3 ≥ 0
 2 x1 - x2 ≥ 0
 x1 ≥ 0, integer
 x2 ≥ 0, integer
 x3 ≥ 0, integer
 Z ≥ 0
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : LO (linear optimization problem)
  Constraints            : 4               
  Cones                  : 0               
  Scalar variables       : 4               
  Matrix variables       : 0               
  Integer variables      : 3               

Optimizer started.
Mixed integer optimizer started.
Threads used: 4
Presolve started.
Presolve terminated. Time = 0.00
Presolved problem: 0 variables, 0 constraints, 0 non-zeros
Presolved problem: 0 general integer, 0 binary, 0 continuous
Clique table size: 0
BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME  
0        1        0        0        3.0000000000e+00     3.0000000000e+00     0.00e+00    0.0   
An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.
The relative gap is 0.00e+00(%).
An optimal solution satisfying the absolute gap tolerance of 0.00e+00 has been located.
The absolute gap is 0.00e+00.

Objective of best integer solution : 3.000000000000e+00      
Best objective bound               : 3.000000000000e+00      
Construct solution objective       : Not employed
Construct solution # roundings     : 0
User objective cut value           : 0
Number of cuts generated           : 0
Number of branches                 : 0
Number of relaxations solved       : 1
Number of interior point iterations: 0
Number of simplex iterations       : 0
Time spend presolving the root     : 0.00
Time spend in the heuristic        : 0.00
Time spend in the sub optimizers   : 0.00
  Time spend optimizing the root   : 0.00
Mixed integer optimizer terminated. Time: 0.00

Optimizer terminated. Time: 0.00    


Integer solution solution summary
  Problem status  : PRIMAL_FEASIBLE
  Solution status : INTEGER_OPTIMAL
  Primal.  obj: 3.0000000000e+00    nrm: 1e+01    Viol.  con: 0e+00    var: 0e+00    itg: 0e+00  
Objective value: 3.0
x1 = 1.0
x2 = 2.0
x3 = 0.0
In [37]:
include("MIGP.jl")
@pyimport runGP
In [61]:
#Some output printing

function print_val(var_name, var_list, value_list)
    i = findfirst(x->occursin(var_name,x), var_list)
    println(var_name*": "*string(value_list[i]))
end
Out[61]:
print_val (generic function with 1 method)
In [67]:
pkg_name = "mico_tests"
mod_name = "Test1"
GP_Time = @elapsed out = solve_MIGP(runGP.run_GP, pkg_name, mod_name)

print_val("x1", out[1], out[2])
print_val("x2", out[1], out[2])
print_val("x3", out[1], out[2])
x1: 1.0
x2: 2.0
x3: 1.0000000000000024e-30
Likely infeasible model at ['INT_x2']
 [1.]
Likely infeasible model at ['INT_x2', 'INT_x3', 'INT_x1']
 [2.e+00 1.e-30 1.e-15]
In [68]:
@show JuMP_Time
@show GP_Time
JuMP_Time = 0.001681676
GP_Time = 0.129456586
Out[68]:
0.129456586

If you can use JuMP (or linear programming) - you should!

Solving this problem using a dedicated MILP solver is significantly faster. However, we are interested in a class of problems not representably as a LP, or currently supported by JuMP. The next notebook (Branch_and_Bound_eVTOL) shows the algorithm functioning on a more interesting example, that takes advantage of the non-linear capability of the GP formulation.