Introduction to Mathematical Programming¶
David E. Bernal NeiraDavidson School of Chemical Engineering, Purdue University
Pedro Maciel Xavier
Davidson School of Chemical Engineering, Purdue University
Azain Khalid
Department of Computer Science, Purdue University
Undergraduate Researcher
Environment Setup¶
If you open this notebook in Google Colab, make sure the runtime is set to Julia before running the next cell.
The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/Project.toml, and install the Julia packages used in this notebook.
import Pkg
IN_COLAB = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP")
function detect_quip_repo_dir()
candidates = (pwd(), normpath(pwd(), ".."))
for candidate in candidates
if isfile(joinpath(candidate, "notebooks_jl", "Project.toml"))
return candidate
end
end
error("Could not locate the QuIP repository root from $(pwd()).")
end
QUIP_REPO_DIR = if IN_COLAB
repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))
if !isdir(repo_dir)
run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
end
repo_dir
else
detect_quip_repo_dir()
end
JULIA_NOTEBOOKS_DIR = joinpath(QUIP_REPO_DIR, "notebooks_jl")
Pkg.activate(JULIA_NOTEBOOKS_DIR)
Pkg.instantiate(; io = devnull)
cd(JULIA_NOTEBOOKS_DIR)
Introduction to Mathematical Programming¶
Modeling¶
The solution to optimization problems requires the development of a mathematical model. Here we will model an example given in the lecture and see how an integer program can be solved practically. This example will use as modeling language JuMP. This open-source Julia package provides flexible access to different solvers and a general modeling framework for linear and nonlinear integer programs. The examples solved here will make use of open-source solvers GLPK and CLP/CBC for linear and mixed-integer linear programming, IPOPT for interior point (non)linear programming, BONMIN for convex integer nonlinear programming, and COUENNE for nonconvex (global) integer nonlinear programming.
Problem statement¶
Suppose there is a company that produces two different products, A and B, which can be sold at different values, and per unit, respectively. The company only counts with a single machine with electricity usage of at most 17kW/day. Producing each A and B consumes and , respectively. Besides, the company can only produce at most 2 more units of B than A per day.
Linear Programming¶
This is a valid model, but it would be easier to solve if we had a mathematical representation. Assuming the units produced of A are and of B are we have
# Plotting
using Plots
function plot_lp(;ns::Integer = 1_000, sol = nothing)
# Create empty plot
plt = plot(;
plot_title="LP Feasible Region",
plot_titlevspan=0.1,
)
# Generate the feasible region plot of this problem
x1 = x2 = range(-0.5, 3.5; length = ns)
# Objective: min 5.5x₁ + 2.1x₂
z(x1, x2) = 5.5x1 + 2.1x2
# Constraints
isfeasible(x1, x2) = (x1 >= 0) && # Bound: x1 ≥ 0
(x2 >= 0) && # Bound: x2 ≥ 0
(x2 <= x1 + 2) && # Constraint 1: x₂ ≤ x₁ + 2
(8x1 + 2x2 <= 17) # Constraint 2: 8 x₁ + 2 x₂ ≤ 17
# Heatmap function
f(x1, x2) = ifelse(isfeasible(x1, x2), z(x1, x2), NaN)
# Plot feasible region
heatmap!(
plt, x1, x2, f;
legend=:topright,
title=raw"$ \max ~ z = 5.5 x_1 + 2.1 x_2 $",
xlims=extrema(x1),
ylims=extrema(x2),
xlabel=raw"$ x_1 $",
ylabel=raw"$ x_2 $",
)
# Make plots of constraints
plot!(plt, x1, (x1) -> x1 + 2;
label=raw"$ x_2 \leq x_1 + 2 $",
color=:green,
)
plot!(plt, x1, (x1) -> (17 - 8x1) / 2;
label=raw"$ 8 x_1 + 2 x_2 \leq 17 $",
color=:blue,
)
# Nonnegativitivy constraints
plot!(plt, zeros(ns), x2;
label=raw"$ x_1 \geq 0 $",
color=:purple,
)
plot!(plt, x1, zeros(ns);
label=raw"$ x_2 \geq 0 $",
color=:red,
)
if sol !== nothing
# Optimal solution
scatter!(plt, [sol[1]], [sol[2]];
label=raw"$ (x_1^\ast, x_2^\ast) $",
marker_z=z,
markershape=:star8,
)
end
return plt
end
plot_lp()Loading...
# Optimization modeling
using JuMP
# Define empty model
lp_model = Model()
# Define the variables
@variable(lp_model, x[1:2] >= 0)
# Define the objective function
@objective(lp_model, Max, 5.5x[1] + 2.1x[2])
# Define the constraints
@constraint(lp_model, c1, x[2] <= x[1] + 2)
@constraint(lp_model, c2, 8x[1] + 2x[2] <= 17)
# Print the model
print(lp_model)Loading...
# Linear programming solver
using GLPK
# Here we solve the optimization problem with GLPK
set_optimizer(lp_model, GLPK.Optimizer)
set_silent(lp_model)
optimize!(lp_model)
# Retrieve solution value for "x"
lp_x = value.(x)
# Display solution of the problem
println(solution_summary(lp_model))
println("* x = $lp_x")solution_summary(; result = 1, verbose = false)
├ solver_name : GLPK
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ ├ raw_status : Solution is optimal
│ └ objective_bound : Inf
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : 1.40800e+01
│ └ dual_objective_value : 1.40800e+01
└ Work counters
└ solve_time (sec) : 2.00009e-03
* x = [1.3, 3.3]
We observe that the optimal solution of this problem is , , leading to a profit of 14.08.
# Mixed-integer solver
using Cbc
# With Cbc instead:
set_optimizer(lp_model, Cbc.Optimizer)
set_silent(lp_model)
optimize!(lp_model)
# Retrieve solution value for "x"
lp_x = value.(x)
# Display solution of the problem
println(solution_summary(lp_model))
println("* x = $lp_x")solution_summary(; result = 1, verbose = false)
├ solver_name : COIN Branch-and-Cut (Cbc)
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ ├ raw_status : Cbc_status = finished - check isProvenOptimal or isProvenInfeasible to see if solution found (or check value of best solution)
Cbc_secondaryStatus = unset (status_ will also be -1)
│ └ objective_bound : 1.40800e+01
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : NO_SOLUTION
│ ├ objective_value : 1.40800e+01
│ └ relative_gap : 0.00000e+00
└ Work counters
├ solve_time (sec) : 9.00006e-03
└ node_count : 0
* x = [1.3, 3.3]
Presolve 0 (-2) rows, 0 (-2) columns and 0 (-4) elements
Optimal - objective value 14.08
After Postsolve, objective 14.08, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 14.08 - 0 iterations time 0.012, Presolve 0.01
plot_lp(;sol = lp_x)Loading...
The solvers GLPK and CLP implement the simplex method (with many improvements) by default, but we can also use an interior point method through the solver IPOPT (interior point optimizer). IPOPT is able to solve not only linear but also nonlinear problems.
# Nonlinear solver
using Ipopt
# Using Ipopt for an interior point method
set_optimizer(lp_model, Ipopt.Optimizer)
set_silent(lp_model)
optimize!(lp_model)
# Retrieve solution value for "x"
lp_x = value.(x)
# Display solution of the problem
println(solution_summary(lp_model))
println("* x = $lp_x")
******************************************************************************
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 https://github.com/coin-or/Ipopt
******************************************************************************
solution_summary(; result = 1, verbose = false)
├ solver_name : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : 1.40800e+01
│ └ dual_objective_value : 1.40800e+01
└ Work counters
├ solve_time (sec) : 2.57000e-01
└ barrier_iterations : 8
* x = [1.3000000135344931, 3.300000029213091]
We obtain the same result as previously, but notice that the interior point method reports a solution subject to a certain tolerance, given by its convergence properties when it can get infinitesimally close (but not directly at) the boundary of the feasible region.
Let’s go back to the slides¶
Integer Programming¶
Now let’s consider that only integer units of each product can be produced, namely
function plot_ilp(;ns::Integer = 1_000, sol = nothing)
# Create empty plot
plt = plot(;
plot_title="ILP Feasible Region",
plot_titlevspan=0.1,
)
# Generate the feasible region plot of this problem
x1 = x2 = range(-0.5, 3.5; length = ns)
# Objective: min 5.5x₁ + 2.1x₂
z(x1, x2) = 5.5x1 + 2.1x2
# Constraints
isfeasible(x1, x2) = (x1 >= 0) && # Bound: x1 ≥ 0
(x2 >= 0) && # Bound: x2 ≥ 0
(x2 <= x1 + 2) && # Constraint 1: x₂ ≤ x₁ + 2
(8x1 + 2x2 <= 17) # Constraint 2: 8 x₁ + 2 x₂ ≤ 17
# Optimal Solution
isoptimal(x1, x2) = (sol !== nothing) && (x1 ≈ sol[1] && x2 ≈ sol[2])
# Heatmap functions
f(x1, x2) = ifelse(isfeasible(x1, x2), z(x1, x2), NaN)
g(x1, x2) = ifelse(isfeasible(x1, x2), 1.0, NaN)
# Plot relaxed feasible region
heatmap!(
plt, x1, x2, f;
legend=:topright,
title=raw"$ z = 5.5 x_1 + 2.1 x_2 $",
xlims=extrema(x1),
ylims=extrema(x2),
xlabel=raw"$ x_1 $",
ylabel=raw"$ x_2 $",
)
# Dim relaxed feasible region
heatmap!(
plt, x1, x2, g;
alpha=0.4,
color=:white,
colorbar_entry=false,
xlims=extrema(x1),
ylims=extrema(x2),
)
# Make plots of constraints
plot!(plt, x1, (x1) -> x1 + 2;
label=raw"$ x_2 \leq x_1 + 2 $",
color=:green,
)
plot!(plt, x1, (x1) -> (17 - 8x1) / 2;
label=raw"$ 8 x_1 + 2 x_2 \leq 17 $",
color=:blue,
)
# Nonnegativitivy constraints
plot!(plt, zeros(ns), x2;
label=raw"$ x_1 \geq 0 $",
color=:purple,
)
plot!(plt, x1, zeros(ns);
label=raw"$ x_2 \geq 0 $",
color=:red,
)
# Feasible solutions
xi = []
xj = []
for i = 0:4, j = 0:4
if isfeasible(i, j) && !isoptimal(i, j)
push!(xi, i)
push!(xj, j)
end
end
scatter!(plt, xi, xj;
label=raw"$ (x_1, x_2) \in \mathbb{Z}^2 $",
marker_z=z,
)
if sol !== nothing
# Optimal solution
scatter!(plt, [sol[1]], [sol[2]];
label=raw"$ (x_1^\ast, x_2^\ast) $",
marker_z=z,
markershape=:star8,
)
end
return plt
end
plot_ilp()Loading...
# Copy model
ilp_model = copy(lp_model)
# Retrieve variable reference
x = ilp_model[:x]
# Add integrality contraint: x₁, x₂ ∈ ℤ
set_integer.(x)
# Print the model
print(ilp_model)Loading...
set_optimizer(ilp_model, Cbc.Optimizer)
set_silent(ilp_model)
optimize!(ilp_model)
print(solution_summary(ilp_model))
# Retrieve
ilp_x = value.(x)
# Display solution of the problem
println(solution_summary(ilp_model))
println("* x = $ilp_x")solution_summary(; result = 1, verbose = false)
├ solver_name : COIN Branch-and-Cut (Cbc)
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ ├ raw_status : Cbc_status = finished - check isProvenOptimal or isProvenInfeasible to see if solution found (or check value of best solution)
Cbc_secondaryStatus = search completed with solution
│ └ objective_bound : 1.18000e+01
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : NO_SOLUTION
│ ├ objective_value : 1.18000e+01
│ └ relative_gap : 0.00000e+00
└ Work counters
├ solve_time (sec) : 3.40002e-02
└ node_count : 0solution_summary(; result = 1, verbose = false)
├ solver_name : COIN Branch-and-Cut (Cbc)
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ ├ raw_status : Cbc_status = finished - check isProvenOptimal or isProvenInfeasible to see if solution found (or check value of best solution)
Cbc_secondaryStatus = search completed with solution
│ └ objective_bound : 1.18000e+01
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : NO_SOLUTION
│ ├ objective_value : 1.18000e+01
│ └ relative_gap : 0.00000e+00
└ Work counters
├ solve_time (sec) : 3.40002e-02
└ node_count : 0
* x = [1.0, 3.0]
Here the solution becomes with an objective of 11.8.
plot_ilp(; sol = ilp_x)Loading...
Let’s go back to the slides¶
Enumeration¶
Enumerating all the possible solutions in this problem might be very efficient (there are only 8 feasible solutions), this we only know from the plot. Assuming that we had as upper bounds for the variables 4, the possible solutions would be 16. With a larger number of variables, the enumeration turns to be impractical. For binary variables (we can always “binarize” the integer variables), the number of possible solutions is .
In many other applications, the possible solutions come from permutations of the integer variables (e.g., assignment problems), which grow as with the size of the input.
This growth makes the problem grow out of control relatively quickly.
# Special mathematical functions used later in the notebook
using SpecialFunctions
function plot_growth(n::Integer = 100)
i = 1:n
plt = plot(;
title=raw"Orders of Magnitude for $ n! $ and $ 2^n $",
xlabel=raw"$ n $",
legend=:topleft,
yscale=:log10,
)
plot!(plt, i, (i) -> 2.0^i; label=raw"$ 2^n $", color=:blue)
plot!(plt, i, SpecialFunctions.gamma; label=raw"$ n! $", color=:red)
plot!(plt, i, (i) -> 3.154E16; color=:gray, label = nothing, linestyle = :dash)
annotate!(plt, first(i), 3.154E16, text("ns in a year", :gray, :left, :bottom, 7))
plot!(plt, i, (i) -> 4.3E26; color=:gray, label = nothing, linestyle = :dash)
annotate!(plt, first(i), 4.3E26, text("age of the universe in ns", :gray, :left, :bottom, 7))
plot!(plt, i, (i) -> 6E79; color=:gray, label = nothing, linestyle = :dash)
annotate!(plt, first(i), 6E79, text("atoms in the universe", :gray, :left, :bottom, 7))
return plt
end
plot_growth()Loading...
Let’s go back to the slides¶
Integer convex nonlinear programming¶
The following constraint: “the production of B minus 1, squared, can only be smaller than 2 minus the production of A” can be incorporated in the following convex integer nonlinear program,
function plot_icnlp(; ns::Integer = 1_000, sol = nothing)
# Create empty plot
plt = plot(;
plot_title="INCLP Feasible Region",
plot_titlevspan=0.1,
)
# Generate the feasible region plot of this problem
x1 = x2 = range(-0.5, 3.5; length = ns)
# Objective: min 5.5x₁ + 2.1x₂
z(x1, x2) = 5.5x1 + 2.1x2
# Constraints
isfeasible(x1, x2) = (x1 >= 0) && # Bound: x1 ≥ 0
(x2 >= 0) && # Bound: x2 ≥ 0
(x2 <= x1 + 2) && # Constraint 1: x₂ ≤ x₁ + 2
(8x1 + 2x2 <= 17) && # Constraint 2: 8 x₁ + 2 x₂ ≤ 17
((x2 - 1)^2 <= 2 - x1) # Constraint 3: (x₂ - 1)² ≤ 2 - x₁
# Optimal Solution
isoptimal(x1, x2) = (sol !== nothing) && (x1 ≈ sol[1] && x2 ≈ sol[2])
# Heatmap functions
f(x1, x2) = ifelse(isfeasible(x1, x2), z(x1, x2), NaN)
g(x1, x2) = ifelse(isfeasible(x1, x2), 1.0, NaN)
# Plot relaxed feasible region
heatmap!(
plt, x1, x2, f;
legend=:topright,
title=raw"$ z = 5.5 x_1 + 2.1 x_2 $",
xlims=extrema(x1),
ylims=extrema(x2),
xlabel=raw"$ x_1 $",
ylabel=raw"$ x_2 $",
)
# Dim relaxed feasible region
heatmap!(
plt, x1, x2, g;
alpha=0.4,
color=:white,
colorbar_entry=false,
xlims=extrema(x1),
ylims=extrema(x2),
)
# Make plots of constraints
plot!(plt, x1, (x1) -> x1 + 2;
label=raw"$ x_2 \leq x_1 + 2 $",
color=:green,
)
plot!(plt, x1, (x1) -> (17 - 8x1) / 2;
label=raw"$ 8 x_1 + 2 x_2 \leq 17 $",
color=:blue,
)
plot!(plt, 2 .- (x1 .- 1) .^ 2, x1,;
label=raw"$ (x_2 - 1)^2 \leq 2 - x_1 $",
color=:orange,
)
# Nonnegativitivy constraints
plot!(plt, zeros(ns), x2;
label=raw"$ x_1 \geq 0 $",
color=:purple,
)
plot!(plt, x1, zeros(ns);
label=raw"$ x_2 \geq 0 $",
color=:red,
)
# Feasible solutions
xi = []
xj = []
for i = 0:4, j = 0:4
if isfeasible(i, j) && !isoptimal(i, j)
push!(xi, i)
push!(xj, j)
end
end
scatter!(plt, xi, xj;
label=raw"$ (x_1, x_2) \in \mathbb{Z}^2 $",
marker_z=z,
)
if sol !== nothing
# Optimal solution
scatter!(plt, [sol[1]], [sol[2]];
label=raw"$ (x_1^\ast, x_2^\ast) $",
marker_z=z,
markershape=:star8,
)
end
return plt
end
plot_icnlp()Loading...
# Define the model
icnlp_model = copy(ilp_model)
# Retrieve variable reference
x = icnlp_model[:x]
# (x₂ - 1)² ≤ 2 - x₁
@constraint(icnlp_model, c3, (x[2] - 1)^2 <= 2 - x[1])
# Print the model
print(icnlp_model)Loading...
# BONMIN solver
using AmplNLWriter
using Bonmin_jll
Bonmin_Optimizer() = AmplNLWriter.Optimizer(Bonmin_jll.amplexe);
set_optimizer(icnlp_model, Bonmin_Optimizer)
optimize!(icnlp_model)
icnlp_x = value.(x)
println(solution_summary(icnlp_model))
println("* x = $icnlp_x")Bonmin 1.8.9 using Cbc 2.10.8 and Ipopt 3.14.13
bonmin:
******************************************************************************
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 https://github.com/coin-or/Ipopt
******************************************************************************
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT -12.775 9 0.004
NLP0012I
Num Status Obj It time Location
NLP0014I 1 INFEAS 0.24999981 16 0.007
NLP0014I 2 OPT -12.4125 5 0.003
NLP0014I 3 OPT -9.7000002 6 0.002
NLP0014I 4 OPT -9.7000001 6 0.002
NLP0014I 5 INFEAS 0.24999981 16 0.008
NLP0014I 6 OPT -9.7000001 6 0.001
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT -9.7 0 0
Cbc0004I Integer solution of -9.7 found after 6 iterations and 0 nodes (0.02 seconds)
Cbc0001I Search completed - best objective -9.699999999999999, took 6 iterations and 0 nodes (0.02 seconds)
Cbc0032I Strong branching done 2 times (33 iterations), fathomed 0 nodes and fixed 1 variables
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
"Finished"
solution_summary(; result = 1, verbose = false)
├ solver_name : AmplNLWriter
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : bonmin: Optimal
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : NO_SOLUTION
│ ├ objective_value : 9.70000e+00
│ └ dual_objective_value : 9.70000e+00
└ Work counters
└ solve_time (sec) : 3.57000e-01
* x = [1.0, 2.0]
In this case the optimal solution becomes with an objective of 9.7.
plot_icnlp(; sol = icnlp_x)Loading...
Integer non-convex programming¶
The last constraint “the production of B minus 1 squared can only be greater than the production of A plus one half” can be incorporated in the following convex integer nonlinear program
function plot_incnlp(; ns::Integer = 1_000, sol = nothing)
# Create empty plot
plt = plot(;
plot_title="INCLP Feasible Region",
plot_titlevspan=0.1,
)
# Generate the feasible region plot of this problem
x1 = x2 = range(-0.5, 3.5; length = ns)
# Objective: min 5.5x₁ + 2.1x₂
z(x1, x2) = 5.5x1 + 2.1x2
# Constraints
isfeasible(x1, x2) = (x1 >= 0) && # Bound: x1 ≥ 0
(x2 >= 0) && # Bound: x2 ≥ 0
(x2 <= x1 + 2) && # Constraint 1: x₂ ≤ x₁ + 2
(8x1 + 2x2 <= 17) && # Constraint 2: 8 x₁ + 2 x₂ ≤ 17
((x2 - 1)^2 <= 2 - x1) && # Constraint 3: (x₂ - 1)² ≤ 2 - x₁
((x2 - 1)^2 >= 1/2 + x1) # Constraint 4: (x₂ - 1)² ≤ 1/2 + x₁
# Optimal Solution
isoptimal(x1, x2) = (sol !== nothing) && (x1 ≈ sol[1] && x2 ≈ sol[2])
# Heatmap functions
f(x1, x2) = ifelse(isfeasible(x1, x2), z(x1, x2), NaN)
g(x1, x2) = ifelse(isfeasible(x1, x2), 1.0, NaN)
# Plot relaxed feasible region
heatmap!(
plt, x1, x2, f;
legend=:topright,
title=raw"$ z = 5.5 x_1 + 2.1 x_2 $",
colobar_label=raw"$ z $",
xlims=extrema(x1),
ylims=extrema(x2),
xlabel=raw"$ x_1 $",
ylabel=raw"$ x_2 $",
)
# Dim relaxed feasible region
heatmap!(
plt, x1, x2, g;
alpha=0.4,
color=:white,
colorbar_entry=false,
xlims=extrema(x1),
ylims=extrema(x2),
)
# Make plots of constraints
plot!(plt, x1, (x1) -> x1 + 2;
label=raw"$ x_2 \leq x_1 + 2 $",
color=:green,
)
plot!(plt, x1, (x1) -> (17 - 8x1) / 2;
label=raw"$ 8 x_1 + 2 x_2 \leq 17 $",
color=:blue,
)
plot!(plt, 2 .- (x1 .- 1) .^ 2, x1,;
label=raw"$ (x_2 - 1)^2 \leq 2 - x_1 $",
color=:orange,
)
plot!(plt, -1/2 .+ (x1 .- 1) .^ 2, x1,;
label=raw"$ (x_2 - 1)^2 \leq \frac{1}{2} + x_1 $",
color=:magenta,
)
# Nonnegativitivy constraints
plot!(plt, zeros(ns), x2;
label=raw"$ x_1 \geq 0 $",
color=:purple,
)
plot!(plt, x1, zeros(ns);
label=raw"$ x_2 \geq 0 $",
color=:red,
)
# Feasible solutions
xi = []
xj = []
for i = 0:4, j = 0:4
if isfeasible(i, j) && !isoptimal(i, j)
push!(xi, i)
push!(xj, j)
end
end
scatter!(plt, xi, xj;
label=raw"$ (x_1, x_2) \in \mathbb{Z}^2 $",
marker_z=z,
)
if sol !== nothing
# Optimal solution
scatter!(plt, [round(sol[1])], [round(sol[2])];
label=raw"$ (x_1^\ast, x_2^\ast) $",
marker_z=z,
markershape=:star8,
)
end
return plt
end
plot_incnlp()Loading...
# Define the model
incnlp_model = copy(icnlp_model)
# Retrieve variable reference
x = incnlp_model[:x]
# (x₂ - 1)² ≥ 1/2 + x₁
@constraint(incnlp_model, c4, (x[2] - 1)^2 >= 1/2 + x[1])
# Print the model
print(incnlp_model)Loading...
# Trying to solve the problem with BONMIN we might obtain the optimal solution,
# but we have no guarantees
set_optimizer(incnlp_model, Bonmin_Optimizer)
optimize!(incnlp_model)
# println(solution_summary(incnlp_model)) # Uncomment to print out the solution summary
if result_count(incnlp_model) > 0
incnlp_x = value.(x)
println("* x = $incnlp_x")
endBonmin 1.8.9 using Cbc 2.10.8 and Ipopt 3.14.13
bonmin:
******************************************************************************
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 https://github.com/coin-or/Ipopt
******************************************************************************
Assertion failed: isNlpFeasible(minlp, primalTolerance), file BonHeuristicDiveMIP.cpp, line 133
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT -2.7500001 7 0.003
# COUENNE solver
using Couenne_jll
Couenne_Optimizer() = AmplNLWriter.Optimizer(Couenne_jll.amplexe);
# Trying to solve the problem with global MINLP solver COUENNE
set_optimizer(incnlp_model, Couenne_Optimizer)
optimize!(incnlp_model)
println(solution_summary(incnlp_model))
if result_count(incnlp_model) > 0
incnlp_x = value.(x)
println("* x = $incnlp_x")
endCouenne 0.5.8 -- an Open-Source solver for Mixed Integer Nonlinear Optimization
Mailing list: couenne@list.coin-or.org
Instructions: http://www.coin-or.org/Couenne
couenne:
ANALYSIS TEST: NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT -2.7500001 8 0.003
Couenne: new cutoff value 0.0000000000e+00 (0.013 seconds)
NLP0014I 2 OPT -0 0 0
Loaded instance "C:\Users\azain\AppData\Local\Temp\jl_kYelFW\model.nl"
Constraints: 4
Variables: 2 (2 integer)
Auxiliaries: 5 (4 integer)
Coin0506I Presolve 6 (-2) rows, 3 (-4) columns and 14 (-8) elements
Clp0006I 0 Obj 0 Dual inf 7.599998 (2)
Clp0006I 4 Obj -6.95
Clp0000I Optimal - objective value -6.95
Clp0032I Optimal objective -6.95 - 4 iterations time 0.002, Presolve 0.00
Clp0000I Optimal - objective value -6.95
NLP Heuristic: Couenne: new cutoff value -4.2000000000e+00 (0.014 seconds)
NLP0014I 3 OPT -4.2 0 0
no solution.
Clp0000I Optimal - objective value -6.95
Optimality Based BT: 0 improved bounds
Probing: 0 improved bounds
NLP Heuristic: no solution.
Cbc0013I At root node, 0 cuts changed objective from -6.95 to -4.2 in 2 passes
Cbc0014I Cut generator 0 (Couenne convexifier cuts) - 0 row cuts average 0.0 elements, 1 column cuts (1 active)
Cbc0004I Integer solution of -4.2 found after 2 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective -4.2, took 2 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
"Finished"
Linearization cuts added at root node: 8
Linearization cuts added in total: 8 (separation time: 0s)
Total solve time: 0.002s (0.002s in branch-and-bound)
Lower bound: -4.2
Upper bound: -4.2 (gap: 0.00%)
Branch-and-bound nodes: 0
solution_summary(; result = 1, verbose = false)
├ solver_name : AmplNLWriter
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : couenne: Optimal
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : NO_SOLUTION
│ ├ objective_value : 4.20000e+00
│ └ dual_objective_value : 4.20000e+00
└ Work counters
└ solve_time (sec) : 1.20000e-01
* x = [0.0, 2.0]
In this case the optimal solution becomes with an objective of 4.2.
plot_incnlp(; sol = incnlp_x)Loading...
We are able to solve non-convex MINLP problems. However, the complexity of these problems leads to significant computational challenges that need to be tackled.
Installing Powerful commercial solvers¶
Gurobi¶
Gurobi is one of the most powerful LP and MIP solvers available today. They provide free academic licenses. In order to install the software, visit their Website, create an account (preferably with your academic email), and obtain a license. Once you do that, you can download and use the software.
BARON¶
BARON is one of the most powerful MINLP solvers available today. Students from the University System of Georgia or CMU and UIUC affiliates are eligible for a free license. In order to install the software, visit their Website, create an account (with your academic email), and obtain a license. Once you do that you can download and use the software.