Поддержка комплексных чисел

На этой странице описываются комплекснозначные переменные и ограничения, поддерживаемые JuMP. Рабочий пример использования этих функций см. в руководстве Quantum state discrimination.

Комплекснозначные переменные

Комплекснозначная переменная создается посредством ComplexPlane:

julia> model = Model();

julia> @variable(model, x in ComplexPlane())
real(x) + imag(x) im

Обратите внимание, что x является не объектом VariableRef, а аффинным выражением с коэффициентами типа Complex{Float64}:

julia> typeof(x)
GenericAffExpr{ComplexF64, VariableRef}

Среда JuMP создала две внутренние вещественные переменные с именами "real(x)" и "imag(x)":

julia> all_variables(model)
2-element Vector{VariableRef}:
 real(x)
 imag(x)

julia> name.(all_variables(model))
2-element Vector{String}:
 "real(x)"
 "imag(x)"

Чтобы получить аффинное выражение с вещественным значением, представляющее каждую переменную, используйте функции real и imag применительно к x:

julia> typeof(real(x))
AffExpr (alias for GenericAffExpr{Float64, GenericVariableRef{Float64}})

julia> typeof(imag(x))
AffExpr (alias for GenericAffExpr{Float64, GenericVariableRef{Float64}})

Чтобы создать анонимную переменную, используйте именованный аргумент set:

julia> model = Model();

julia> x = @variable(model, set = ComplexPlane())
_[1] + _[2] im

Комплекснозначная переменная и границы начальных значений

Так как у комплекснозначных переменных нет общего порядка, определение границы комплекснозначной переменной неоднозначно. При передаче вещественных или комплексных значений в таких именованных аргументах, как lower_bound, upper_bound и start_value, JuMP применяет вещественные и мнимые части к соответствующим вещественным переменным.

julia> model = Model();

julia> @variable(
           model,
           x in ComplexPlane(),
           lower_bound = 1.0,
           upper_bound = 2.0 + 3.0im,
           start = 4im,
       )
real(x) + imag(x) im

julia> vars = all_variables(model)
2-element Vector{VariableRef}:
 real(x)
 imag(x)

julia> lower_bound.(vars)
2-element Vector{Float64}:
 1.0
 0.0

julia> upper_bound.(vars)
2-element Vector{Float64}:
 2.0
 3.0

julia> start_value.(vars)
2-element Vector{Float64}:
 0.0
 4.0

Вы можете изменить границы и начальные значения, передав imag(x) или real(x) в соответствующую функцию:

julia> set_lower_bound(imag(x), 2)

julia> lower_bound(imag(x))
2.0

julia> delete_upper_bound(real(x))

julia> has_upper_bound(real(x))
false

julia> set_start_value(imag(x), 3)

julia> start_value(imag(x))
3.0

Комплекснозначные ограничения равенства

JuMP преобразует комплекснозначные ограничения равенства в два ограничения с вещественными значениями: одно из них представляет вещественную часть, а другое — мнимую. Таким образом, комплекснозначные ограничения равенства могут решаться любым решателем, который поддерживает тип ограничений с вещественными значениями.

Например:

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x[1:2]);

julia> @constraint(model, (1 + 2im) * x[1] + 3 * x[2] == 4 + 5im)
(1 + 2im) x[1] + 3 x[2] = (4 + 5im)

julia> optimize!(model)

julia> value.(x)
2-element Vector{Float64}:
 2.5
 0.5

эквивалентно

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x[1:2]);

julia> @constraint(model, 1 * x[1] + 3 * x[2] == 4)  # вещественная часть
x[1] + 3 x[2] = 4

julia> @constraint(model, 2 * x[1] == 5)  # мнимая часть
2 x[1] = 5

julia> optimize!(model)

julia> value.(x)
2-element Vector{Float64}:
 2.5
 0.5

Это также верно в отношении комплекснозначных переменных:

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x in ComplexPlane());

julia> @constraint(model, (1 + 2im) * x + 3 * x == 4 + 5im)
(4 + 2im) real(x) + (-2 + 4im) imag(x) = (4 + 5im)

julia> optimize!(model)

julia> value(x)
1.3 + 0.6000000000000001im

что равносильно

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x_real);

julia> @variable(model, x_imag);

julia> @constraint(model, x_real - 2 * x_imag + 3 * x_real == 4)
4 x_real - 2 x_imag = 4

julia> @constraint(model, x_imag + 2 * x_real + 3 * x_imag == 5)
2 x_real + 4 x_imag = 5

julia> optimize!(model)

julia> value(x_real) + value(x_imag) * im
1.3 + 0.6000000000000001im

Эрмитовы положительно полуопределенные конусы

JuMP поддерживает создание эрмитовых матриц.

julia> model = Model();

julia> @variable(model, H[1:3, 1:3] in HermitianPSDCone())
3×3 LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 real(H[1,1])                    …  real(H[1,3]) + imag(H[1,3]) im
 real(H[1,2]) - imag(H[1,2]) im     real(H[2,3]) + imag(H[2,3]) im
 real(H[1,3]) - imag(H[1,3]) im     real(H[3,3])

Среда JuMP создала девять внутренних вещественных переменных решения:

julia> all_variables(model)
9-element Vector{VariableRef}:
 real(H[1,1])
 real(H[1,2])
 real(H[2,2])
 real(H[1,3])
 real(H[2,3])
 real(H[3,3])
 imag(H[1,2])
 imag(H[1,3])
 imag(H[2,3])

и ограничение Vector{VariableRef}-in-MOI.HermitianPositiveSemidefiniteConeTriangle:

julia> num_constraints(model, Vector{VariableRef}, MOI.HermitianPositiveSemidefiniteConeTriangle)
1

Множество MOI.HermitianPositiveSemidefiniteConeTriangle можно эффективно связать мостом с MOI.PositiveSemidefiniteConeTriangle, поэтому оно может быть решено любым решателем, поддерживающим положительно полуопределенные ограничения.

Каждый элемент H является аффинным выражением с коэффициентами типа Complex{Float64}:

julia> typeof(H[1, 1])
GenericAffExpr{ComplexF64, VariableRef}

julia> typeof(H[2, 1])
GenericAffExpr{ComplexF64, VariableRef}

Начальные значения

При задании начального значения будьте внимательны: задавать следует только верхний треугольник вещественных переменных и верхний треугольник мнимых переменных, исключая диагональ:

julia> import LinearAlgebra

julia> function set_hermitian_start(
           H::LinearAlgebra.Hermitian,
           start::LinearAlgebra.Hermitian,
       )
           for j in 1:size(H, 2), i in 1:j
               set_start_value(real(H[i, j]), real(start[i, j]))
               if i < j
                   set_start_value(imag(H[i, j]), imag(start[i, j]))
               end
           end
           return
       end
set_hermitian_start (generic function with 1 method)

julia> H0 = LinearAlgebra.Hermitian(
           [1 (2+3im) (5+6im); (2-3im) 4 (7+8im); (5-6im) (7-8im) 9],
       )
3×3 LinearAlgebra.Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}:
 1+0im  2+3im  5+6im
 2-3im  4+0im  7+8im
 5-6im  7-8im  9+0im

julia> set_hermitian_start(H, H0)

julia> start_value.(all_variables(model))
9-element Vector{Float64}:
 1.0
 2.0
 4.0
 5.0
 7.0
 9.0
 3.0
 6.0
 8.0

Эрмитовы положительно полуопределенные ограничения

HermitianPSDCone также можно использовать в макросе @constraint:

julia> model = Model();

julia> @variable(model, x[1:2])
2-element Vector{VariableRef}:
 x[1]
 x[2]

julia> import LinearAlgebra

julia> H = LinearAlgebra.Hermitian([x[1] 1im; -1im -x[2]])
2×2 LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 x[1]  im
 -im   -x[2]

julia> @constraint(model, H in HermitianPSDCone())
[x[1]  im
 -im   -x[2]] ∈ HermitianPSDCone()

Матрица H в H in HermitianPSDCone() должна быть типа LinearAlgebra.Hermitian. Если матрица имеет другой тип, выдается ошибка build_constraint.