关于矩阵的 n 次方的 Julia 约束

Question

我正在使用 Julia JuMP 求解器，并且正在尝试编写一个涉及矩阵幂的约束。

我的变量是一个方阵M和一个整数n。 M 的大小为 s。我们有s>n。

M 是一个方阵，我们可以将其视为具有 s 个顶点的图的邻接矩阵。我想断言这个图有 n 个顶点，它们是长度为 n 的循环的一部分。

所以我添加了 JuMP 约束

@constraint(model, tr(M^n)==n)

当我尝试运行求解器时，出现此错误，这似乎是由M^n 引起的：

  MethodError: zero(::Type{Union{}}) is ambiguous. Candidates:
  zero(::Type{var"#s828"} where var"#s828"<:AbstractIrrational) in Base at irrationals.jl:149
  zero(::Union{Type{P}, P}) where P<:Dates.Period in Dates at /opt/julia-1.5.2/share/julia/stdlib/v1.5/Dates/src/periods.jl:53
  zero(::Type{T}) where T<:Number in Base at number.jl:242
  zero(F::Type{var"#s280"} where var"#s280"<:Union{MathOptInterface.ScalarAffineFunction{T}, MathOptInterface.ScalarQuadraticFunction{T}}) where T in MathOptInterface.Utilities at /home/n7student/.julia/packages/MathOptInterface/VjkSQ/src/Utilities/functions.jl:1987
  zero(a::Type{SA}) where SA<:StaticArrays.StaticArray in StaticArrays at /home/n7student/.julia/packages/StaticArrays/LJQEe/src/linalg.jl:88
  zero(::Type{V}) where V<:AbstractVariableRef in JuMP at /home/n7student/.julia/packages/JuMP/qhoVb/src/JuMP.jl:734
Possible fix, define
  zero(::Type{Union{}})

Stacktrace:
 [1] mutable_operate!(::typeof(zero), ::Array{Union{},2}) at /home/n7student/.julia/packages/MutableArithmetics/0tlz5/src/linear_algebra.jl:202
 [2] mutable_operate_to! at /home/n7student/.julia/packages/MutableArithmetics/0tlz5/src/linear_algebra.jl:207 [inlined]
 [3] _mul! at /home/n7student/.julia/packages/MutableArithmetics/0tlz5/src/dispatch.jl:63 [inlined]
 [4] mul! at /home/n7student/.julia/packages/MutableArithmetics/0tlz5/src/dispatch.jl:68 [inlined]
 [5] * at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:153 [inlined]
 [6] power_by_squaring(::Array{VariableRef,2}, ::Int64) at ./intfuncs.jl:245
 [7] ^ at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/dense.jl:408 [inlined]
 [8] macro expansion at ./none:0 [inlined]
 [9] literal_pow(::typeof(^), ::Array{VariableRef,2}, ::Val{2}) at ./none:0
 [10] macro expansion at /home/n7student/.julia/packages/MutableArithmetics/0tlz5/src/rewrite.jl:227 [inlined]
 [11] macro expansion at /home/n7student/.julia/packages/JuMP/qhoVb/src/macros.jl:440 [inlined]
 [12] top-level scope at ./In[89]:55
 [13] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091
 [14] execute_code(::String, ::String) at /home/n7student/.julia/packages/IJulia/a1SNk/src/execute_request.jl:27
 [15] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/n7student/.julia/packages/IJulia/a1SNk/src/execute_request.jl:86
 [16] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [17] invokelatest at ./essentials.jl:709 [inlined]
 [18] eventloop(::ZMQ.Socket) at /home/n7student/.julia/packages/IJulia/a1SNk/src/eventloop.jl:8
 [19] (::IJulia.var"#15#18")() at ./task.jl:356

Answer 1

您只能在@constraint 中形成线性或二次约束。 正常的解决方法是定义一系列约束

model = Model()
@variable(model, x[i=1:2, j=1:2, k=1:3])
@constraint(model, [k=1:2], x[:, :, k+1] .== x[:, :, k] * x[:, :, 1])

或者，使用非线性用户定义函数，例如：

using JuMP, LinearAlgebra, Ipopt
s, n = 3, 2
f(x...) = LinearAlgebra.tr(reshape(collect(x), (s, s))^n)
model = Model(Ipopt.Optimizer)
register(model, :f, s^2, f; autodiff = true)
@variable(model, x[1:s, 1:s])
@NLconstraint(model, f(x...) == n)
optimize!(model)