using LightGraphs
using TikzGraphs, TikzPictures
using PyCall
include("MIGP.jl")
$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.
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
The code above is saved in the file mico_test.py, which is imported below.
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.
G = Graph()
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
function addNode!(Node, G, Node_List)
add_vertex!(G)
push!(Node_List,Node)
end
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.
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")
This solution is the base node of the tree which will be searched to determine the optimal solution.
dump(N0)
These functions are supporting funtions to implement the branch and bound algorithm.
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
function nearest_ints(val)
rounded = round(val)
if rounded >= val
return (rounded, rounded-1)
else
return (rounded+1,rounded)
end
end
function isInt(i)
tol = 1e-3
return min(i - floor(i), ceil(i)-i) < tol
end
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
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
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
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
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.
g(vars, vals) = mt.run_Test1(vars, vals)
(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"
drawgraph(G, Node_List, cols)
The lowest-cost node is then expanded, and the process repeated for all remaining integer variables
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)
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/
using JuMP
using Mosek
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))
include("MIGP.jl")
@pyimport runGP
#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
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])
@show JuMP_Time
@show GP_Time
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.