Arbitrary precision arithmetic
This tutorial was generated using Literate.jl. Download the source as a .jl
file.
The purpose of this tutorial is to explain how to use a solver which supports arithmetic using a number type other than Float64
.
Required packages
This tutorial uses the following packages:
using JuMP
import CDDLib
import Clarabel
Higher-precision arithmetic
To create a model with a number type other than Float64
, use GenericModel
with an optimizer which supports the same number type:
model = GenericModel{BigFloat}(Clarabel.Optimizer{BigFloat})
A JuMP Model
├ value_type: BigFloat
├ solver: Clarabel
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
The syntax for adding decision variables is the same as a normal JuMP model, except that values are converted to BigFloat
:
@variable(model, -1 <= x[1:2] <= sqrt(big"2"))
2-element Vector{GenericVariableRef{BigFloat}}:
x[1]
x[2]
Note that each x
is now a GenericVariableRef{BigFloat}
, which means that the value of x
in a solution will be a BigFloat
.
The lower and upper bounds of the decision variables are also BigFloat
:
lower_bound(x[1])
-1.0
typeof(lower_bound(x[1]))
BigFloat
upper_bound(x[2])
1.414213562373095048801688724209698078569671875376948073176679737990732478462102
typeof(upper_bound(x[2]))
BigFloat
The syntax for adding constraints is the same as a normal JuMP model, except that coefficients are converted to BigFloat
:
@constraint(model, c, x[1] == big"2" * x[2])
\[ x_{1} - 2.0 x_{2} = 0.0 \]
The function is a GenericAffExpr
with BigFloat
for the coefficient and variable types;
constraint = constraint_object(c)
typeof(constraint.func)
GenericAffExpr{BigFloat, GenericVariableRef{BigFloat}}
and the set is a MOI.EqualTo{BigFloat}
:
typeof(constraint.set)
MathOptInterface.EqualTo{BigFloat}
The syntax for adding and objective is the same as a normal JuMP model, except that coefficients are converted to BigFloat
:
@objective(model, Min, 3x[1]^2 + 2x[2]^2 - x[1] - big"4" * x[2])
\[ 3.0 x_{1}^2 + 2.0 x_{2}^2 - x_{1} - 4.0 x_{2} \]
typeof(objective_function(model))
GenericQuadExpr{BigFloat, GenericVariableRef{BigFloat}}
Here's the model we have built:
print(model)
Min 3.0 x[1]² + 2.0 x[2]² - x[1] - 4.0 x[2]
Subject to
c : x[1] - 2.0 x[2] = 0.0
x[1] ≥ -1.0
x[2] ≥ -1.0
x[1] ≤ 1.414213562373095048801688724209698078569671875376948073176679737990732478462102
x[2] ≤ 1.414213562373095048801688724209698078569671875376948073176679737990732478462102
Let's solve and inspect the solution:
optimize!(model)
@assert is_solved_and_feasible(model; dual = true)
solution_summary(model)
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -6.42857e-01
Dual objective value : -6.42857e-01
* Work counters
Solve time (sec) : 4.24953e-01
Barrier iterations : 5
The value of each decision variable is a BigFloat
:
value.(x)
2-element Vector{BigFloat}:
0.4285714246558161076147072906813123533593766450416896337912086518811186790735189
0.2142857123279078924828007272730108809297577877991360649674411645247653239673801
as well as other solution attributes like the objective value:
objective_value(model)
-0.6428571428571422964607590389935242587959291815638830868454759876473734138856053
and dual solution:
dual(c)
1.571428571977140845343978069015092190548250919787945065022059071052557047888015
This problem has an analytic solution of x = [3//7, 3//14]
. Currently, our solution has an error of approximately 1e-9
:
value.(x) .- [3 // 7, 3 // 14]
2-element Vector{BigFloat}:
-3.915612463813864137890116218069194783529738937637362776690309892355053792476207e-09
-1.957806393231484987012703404784527926486578220746844549760948961746906215408591e-09
But by reducing the tolerances, we can obtain a more accurate solution:
set_attribute(model, "tol_gap_abs", 1e-32)
set_attribute(model, "tol_gap_rel", 1e-32)
optimize!(model)
@assert is_solved_and_feasible(model)
value.(x) .- [3 // 7, 3 // 14]
2-element Vector{BigFloat}:
-4.120732596246374574619292889406407106157605546563218305172773512099467866195165e-32
-7.146646610782677659152301436088423235900252780211057986251367981130623553333357e-32
Rational arithmetic
In addition to higher-precision floating point number types like BigFloat
, JuMP also supports solvers with exact rational arithmetic. One example is CDDLib.jl, which supports the Rational{BigInt}
number type:
model = GenericModel{Rational{BigInt}}(CDDLib.Optimizer{Rational{BigInt}})
A JuMP Model
├ value_type: Rational{BigInt}
├ solver: CDD
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
As before, we can create variables using rational bounds:
@variable(model, 1 // 7 <= x[1:2] <= 2 // 3)
2-element Vector{GenericVariableRef{Rational{BigInt}}}:
x[1]
x[2]
lower_bound(x[1])
1//7
typeof(lower_bound(x[1]))
Rational{BigInt}
As well as constraints:
@constraint(model, c1, (2 // 1) * x[1] + x[2] <= 1)
\[ 2//1 x_{1} + x_{2} \leq 1//1 \]
@constraint(model, c2, x[1] + 3x[2] <= 9 // 4)
\[ x_{1} + 3//1 x_{2} \leq 9//4 \]
and objective functions:
@objective(model, Max, sum(x))
\[ x_{1} + x_{2} \]
Here's the model we have built:
print(model)
Max x[1] + x[2]
Subject to
c1 : 2//1 x[1] + x[2] ≤ 1//1
c2 : x[1] + 3//1 x[2] ≤ 9//4
x[1] ≥ 1//7
x[2] ≥ 1//7
x[1] ≤ 2//3
x[2] ≤ 2//3
Let's solve and inspect the solution:
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)
* Solver : CDD
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"Optimal"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Objective value : 5//6
* Work counters
The optimal values are given in exact rational arithmetic:
value.(x)
2-element Vector{Rational{BigInt}}:
1//6
2//3
objective_value(model)
5//6
value(c2)
13//6