Арифметика произвольной точности

Цель данного руководства — объяснить, как использовать решатель, поддерживающий арифметические операции с числовым типом, отличным от Float64.

В этом руководстве используются следующие пакеты.

using JuMP
import CDDLib
import Clarabel

Арифметика более высокой точности

Чтобы создать модель с числовым типом, отличным от Float64, используйте GenericModel с оптимизатором, который поддерживает этот числовой тип:

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

Синтаксис добавления переменных решения такой же, как и в обычной модели JuMP, но значения преобразуются в BigFloat:

@variable(model, -1 <= x[1:2] <= sqrt(big"2"))

2-element Vector{GenericVariableRef{BigFloat}}:
 x[1]
 x[2]

Обратите внимание, что каждая переменная x теперь является GenericVariableRef{BigFloat}, то есть значение x в решении будет иметь тип BigFloat.

Нижние и верхние границы переменных решения также имеют тип 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

Синтаксис добавления ограничений такой же, как и в обычной модели JuMP, но коэффициенты преобразуются в BigFloat:

@constraint(model, c, x[1] == big"2" * x[2])

c : x[1] - 2.0 x[2] = 0.0

Функция представляет собой GenericAffExpr с типом BigFloat для коэффициентов и переменных:

constraint = constraint_object(c)
typeof(constraint.func)

GenericAffExpr{BigFloat, GenericVariableRef{BigFloat}}

а множество представляет собой MOI.EqualTo{BigFloat}:

typeof(constraint.set)

MathOptInterface.EqualTo{BigFloat}

Синтаксис добавления целевой функции такой же, как и в обычной модели JuMP, но коэффициенты преобразуются в BigFloat:

@objective(model, Min, 3x[1]^2 + 2x[2]^2 - x[1] - big"4" * x[2])

3.0 x[1]² + 2.0 x[2]² - x[1] - 4.0 x[2]

typeof(objective_function(model))

GenericQuadExpr{BigFloat, GenericVariableRef{BigFloat}}

Вот построенная модель:

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

Давайте получим и проверим решение:

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.65813e+00
  Barrier iterations : 5

Значение каждой переменной решения имеет тип BigFloat:

value.(x)

2-element Vector{BigFloat}:
 0.4285714246558161076147072906813123533593766450416896337912086518811186790735189
 0.2142857123279078924828007272730108809297577877991360649674411645247653239673801

Этот же тип имеют и другие атрибуты решения, например целевое значение:

objective_value(model)

-0.6428571428571422964607590389935242587959291815638830868454759876473734138856053

и двойственное решение:

dual(c)

1.571428571977140845343978069015092190548250919787945065022059071052557047888015

Эта задача имеет аналитическое решение x = [3//7, 3//14]. В настоящее время решение имеет погрешность примерно 1e-9:

value.(x) .- [3 // 7, 3 // 14]

2-element Vector{BigFloat}:
 -3.915612463813864137890116218069194783529738937637362776690309892355053792476207e-09
 -1.957806393231484987012703404784527926486578220746844549760948961746906215408591e-09

Однако, уменьшив допуски, можно получить более точное решение:

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

Арифметические операции с рациональными числами

Помимо типов чисел с плавающей запятой высокой точности, таких как BigFloat, JuMP также поддерживает решатели с точными арифметическими операциями с рациональными числами. Одним из примеров является библиотека CDDLib.jl, которая поддерживает числовой тип Rational{BigInt}:

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

Как и ранее, мы можем создать переменные с рациональными границами:

@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}

А также ограничения:

@constraint(model, c1, (2 // 1) * x[1] + x[2] <= 1)

c1 : 2//1 x[1] + x[2] ≤ 1//1

@constraint(model, c2, x[1] + 3x[2] <= 9 // 4)

c2 : x[1] + 3//1 x[2] ≤ 9//4

и целевые функции:

@objective(model, Max, sum(x))

x[1] + x[2]

Вот построенная модель:

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

Давайте получим и проверим решение:

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

Оптимальные значения даны в виде точных рациональных чисел:

value.(x)

2-element Vector{Rational{BigInt}}:
 1//6
 2//3

objective_value(model)

5//6

value(c2)

13//6