Complex number support
JuMP has a limited range of support for complex-valued variables and constraints. Since no current solvers natively support complex-valued variables and constraints, JuMP reformulates all complex expressions into their real and imaginary parts.
Complex-valued variables
Create a complex-valued variable using ComplexPlane
:
julia> model = Model();
julia> @variable(model, x in ComplexPlane())
real(x) + (0.0 + 1.0im) imag(x)
Note that x
is not a VariableRef
; instead, it is an affine expression with Complex{Float64}
-valued coefficients:
julia> typeof(x)
GenericAffExpr{ComplexF64, VariableRef}
Behind the scenes, JuMP has created two real-valued variables, with names "real(x)"
and "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)"
Use the real
and imag
functions on x
to return a real-valued affine expression representing each variable:
julia> typeof(real(x))
AffExpr (alias for GenericAffExpr{Float64, VariableRef})
julia> typeof(imag(x))
AffExpr (alias for GenericAffExpr{Float64, VariableRef})
To create an anonymous variable, use the set
keyword argument:
julia> model = Model();
julia> x = @variable(model, set = ComplexPlane())
_[1] + (0.0 + 1.0im) _[2]
Complex-valued variable bounds
Because complex-valued variables lack a total ordering, the definition of a variable bound for a complex-valued variable is ambiguous. If you pass a real- or complex-valued argument to keywords such as lower_bound
, upper_bound
, and start_value
, JuMP will apply the real and imaginary parts to the associated real-valued variables.
julia> model = Model();
julia> @variable(
model,
x in ComplexPlane(),
lower_bound = 1.0,
upper_bound = 2.0 + 3.0im,
start = 4im,
)
real(x) + (0.0 + 1.0im) imag(x)
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
Complex-valued equality constraints
JuMP reformulates complex-valued equality constraints into two real-valued constraints: one representing the real part, and one representing the imaginary part. Thus, complex-valued equality constraints can be solved any solver that supports the real-valued constraint type.
For example:
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.0 + 2.0im) x[1] + (3.0 + 0.0im) x[2] = 4.0 + 5.0im
julia> optimize!(model)
julia> value.(x)
2-element Vector{Float64}:
2.5
0.5
is equivalent to
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) # real component
x[1] + 3 x[2] = 4.0
julia> @constraint(model, 2 * x[1] == 5) # imag component
2 x[1] = 5.0
julia> optimize!(model)
julia> value.(x)
2-element Vector{Float64}:
2.5
0.5
This also applies if the variables are complex-valued:
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.0 + 2.0im) real(x) + (-2.0 + 4.0im) imag(x) = 4.0 + 5.0im
julia> optimize!(model)
julia> value(x)
1.3 + 0.6000000000000001im
which is equivalent to
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.0
julia> @constraint(model, x_imag + 2 * x_real + 3 * x_imag == 5)
2 x_real + 4 x_imag = 5.0
julia> optimize!(model)
julia> value(x_real) + value(x_imag) * im
1.3 + 0.6000000000000001im
Hermitian PSD Cones
JuMP supports creating matrices where are Hermitian.
julia> model = Model();
julia> @variable(model, H[1:3, 1:3] in HermitianPSDCone())
3×3 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
real(H[1,1]) … real(H[1,3]) + (0.0 + 1.0im) imag(H[1,3])
real(H[1,2]) + (-0.0 - 1.0im) imag(H[1,2]) real(H[2,3]) + (0.0 + 1.0im) imag(H[2,3])
real(H[1,3]) + (-0.0 - 1.0im) imag(H[1,3]) real(H[3,3])
Behind the scenes, JuMP has created nine real-valued decision variables:
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])
and a Vector{VariableRef}-in-MOI.HermitianPositiveSemidefiniteConeTriangle
constraint:
julia> num_constraints(model, Vector{VariableRef}, MOI.HermitianPositiveSemidefiniteConeTriangle)
1
The MOI.HermitianPositiveSemidefiniteConeTriangle
set can be efficiently bridged to MOI.PositiveSemidefiniteConeTriangle
, so it can be solved by any solver that supports PSD constraints.
Each element of H
is an affine expression with Complex{Float64}
-valued coefficients:
julia> typeof(H[1, 1])
GenericAffExpr{ComplexF64, VariableRef}
julia> typeof(H[2, 1])
GenericAffExpr{ComplexF64, VariableRef}
Hermitian PSD constraints
The HermitianPSDCone
can also be used in the @constraint
macro:
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] (0.0 + 1.0im)
(0.0 - 1.0im) (-1.0 - 0.0im) x[2]
julia> @constraint(model, H in HermitianPSDCone())
[x[1] (0.0 + 1.0im);
(0.0 - 1.0im) (-1.0 + 0.0im) x[2]] ∈ HermitianPSDCone()