In :
using LightGraphs
using TikzGraphs, TikzPictures
using PyCall
include("MIGP.jl")

Out:
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}

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

In :
PyCall.pyversion
@pyimport mico_tests as mt

In :
G = Graph()

Out:
{0, 0} undirected simple Int64 graph
In :
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 :
function addNode!(Node, G, Node_List)
push!(Node_List,Node)
end

Out:
addNode! (generic function with 1 method)

## Solve Relaxed Problem¶

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

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  - 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:
In :
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¶

In :
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,"_")
if key == "INT"
push!(int_vars,var)
push!(int_vals,vals[i])
end
end
return int_vars,int_vals
end

Out:
identify_INT_vars (generic function with 1 method)
In :
function nearest_ints(val)
rounded = round(val)

if rounded >= val
return (rounded, rounded-1)
else
return (rounded+1,rounded)
end
end

Out:
nearest_ints (generic function with 1 method)
In :
function isInt(i)
tol = 1e-3
return min(i - floor(i), ceil(i)-i) < tol
end

Out:
isInt (generic function with 1 method)
In :
function cleanup_input(varname::Array{String})
output = String[]

for str in varname
trim = split(str, '.')
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:
cleanup_input (generic function with 2 methods)
In :
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)
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)

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)

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)

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
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
remaining_vars, remaining_values = identify_INT_vars(vars_above,vals_above)
best_child = Nabove
worst_child = Nbelow

elseif cost_above <= cost_below
best_val = above
best_cost = cost_above
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
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:
constrain_variable (generic function with 1 method)
In :
Int_Vars, Int_Vals = identify_INT_vars(cleanup_input(N0.vars),N0.vals)
cols = fill("fill=red!10",2*length(Int_Vars)+1)
cols = "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
val2expand   = Int_Vals

(best_child, worst_child) = constrain_variable(ParentNode, Graph, Int_Vars, Int_Vals,g,
var2expand, val2expand)

end

Out:
expand! (generic function with 1 method)
In :
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:
drawgraph (generic function with 1 method)

## Constrain the first integer variable¶

In :
g(vars, vals) = mt.run_Test1(vars, vals)

Out:
g (generic function with 1 method)
In :
(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:
"fill=blue!10"
In :
drawgraph(G, Node_List, cols)

Out:

## Expand subsequent nodes¶

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

# Compare Timing With JuMP¶

In :
using JuMP
using Mosek

In :
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.
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 :
include("MIGP.jl")
@pyimport runGP

In :
#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:
print_val (generic function with 1 method)
In :
pkg_name = "mico_tests"
mod_name = "Test1"
GP_Time = @elapsed out = solve_MIGP(runGP.run_GP, pkg_name, mod_name)

print_val("x1", out, out)
print_val("x2", out, out)
print_val("x3", out, out)

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]


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 :
@show JuMP_Time
@show GP_Time

JuMP_Time = 0.001681676
GP_Time = 0.129456586

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

In [ ]: