Julia JuMP StackOverflow 在定义线性矩阵不等式约束时出错
Julia JuMP StackOverflow Error when defining a constraint on linear matrix inequalities
我受到 this post 的启发,从 MATLAB 的 LMI 工具箱切换到 Julia 实现,我尝试将 post 中建议的脚本改编为我试图解决的问题。但是,我没有做到这一点,在尝试使用 JuMP 解决以下 LMI 问题时,一直导致 Julia 崩溃或给我一个 WhosebugError
,
Lyapunov stability of hybrid systems
对于未知对称矩阵Ui,Wi、M(和Pi),其中Ui,Wi大于或等于0。
我已经 运行 逐行设置模型,发现它在我定义最后一个约束时发生:
@SDconstraint(model, Ai'*Pi + Pi*Ai + Ei'*Ui*Ei ⪯ 0)
我在 JuMP 上看到了其他 post 的相关信息,它们通常是变量定义不当的结果,但这些知识并没有让我弄清楚这个问题。这是产生错误的代码的最小形式,
using JuMP
using MosekTools
using LinearAlgebra
E1p = [1 0; 0 1];
A1 = [-1 10; -100 -1];
Ei = E1p
Ai = A1
n=2
model = Model(Mosek.Optimizer)
@variable(model, Pi[i=1:n, j=1:n], Symmetric)
@variable(model, Ui[i=1:n, j=1:n], Symmetric)
@variable(model, Wi[i=1:n, j=1:n], Symmetric)
@variable(model, M[i=1:n, j=1:n], Symmetric)
@SDconstraint(model, Ui ⪰ 0)
@SDconstraint(model, Wi ⪰ 0)
@SDconstraint(model, Ei'*M*Ei ⪰ 0)
@SDconstraint(model, Pi - Ei'*Wi*Ei ⪰ 0)
@SDconstraint(model, Ai'*Pi + Pi*Ai + Ei'*Ui*Ei ⪯ 0)
optimize!(model)
Pi = value.(Pi)
重复很长时间的错误的第一行是:
in mutable_operate! at MutableArithmetics/bPWR4/src/linear_algebra.jl:132
这看起来像是一个错误。
这是一个 MWE
using JuMP
model = Model()
E = [1 1; 1 1]
@variable(model, P[1:size(E, 1), 1:size(E, 2)], Symmetric)
julia> @SDconstraint(model, E' * P + P * E + E' * P * E <= 0)
ERROR: WhosebugError:
Stacktrace:
[1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64}) (repeats 79984 times)
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/linear_algebra.jl:132
如果我放弃第一个学期,我会得到
julia> @SDconstraint(model, P * E + E' * P * E <= 0)
ERROR: MethodError: no method matching *(::Matrix{AffExpr})
Closest candidates are:
*(::StridedMatrix{T} where T, ::LinearAlgebra.AbstractQ) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:680
*(::StridedMatrix{T} where T, ::LinearAlgebra.Adjoint{var"#s832", var"#s831"} where {var"#s832", var"#s831"<:LinearAlgebra.AbstractQ}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:720
*(::StridedVecOrMat{T} where T, ::LinearAlgebra.Adjoint{var"#s832", var"#s831"} where {var"#s832", var"#s831"<:LinearAlgebra.LQPackedQ}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lq.jl:254
...
Stacktrace:
[1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::Matrix{AffExpr}) (repeats 2 times)
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/linear_algebra.jl:132
[2] operate_fallback!(::MutableArithmetics.IsMutable, ::Function, ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/interface.jl:330
[3] operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/rewrite.jl:80
[4] macro expansion
@ ~/.julia/packages/MutableArithmetics/bPWR4/src/rewrite.jl:276 [inlined]
[5] macro expansion
@ ~/.julia/packages/JuMP/y5vgk/src/macros.jl:447 [inlined]
[6] top-level scope
@ REPL[15]:1
所以这可能与https://github.com/jump-dev/MutableArithmetics.jl/issues/84有关。
我会开一个问题。 (编辑:https://github.com/jump-dev/MutableArithmetics.jl/issues/86)
作为解决方法,首先将其拆分为一个表达式:
julia> @expression(model, ex, E' * P + P * E + E' * P * E)
2×2 Matrix{AffExpr}:
3 P[1,1] + 4 P[1,2] + P[2,2] 4 P[1,2] + 2 P[2,2] + 2 P[1,1]
2 P[1,1] + 4 P[1,2] + 2 P[2,2] 4 P[1,2] + 3 P[2,2] + P[1,1]
julia> @SDconstraint(model, ex <= 0)
[-3 P[1,1] - 4 P[1,2] - P[2,2] -2 P[1,1] - 4 P[1,2] - 2 P[2,2];
-2 P[1,1] - 4 P[1,2] - 2 P[2,2] -P[1,1] - 4 P[1,2] - 3 P[2,2]] ∈ PSDCone()
我受到 this post 的启发,从 MATLAB 的 LMI 工具箱切换到 Julia 实现,我尝试将 post 中建议的脚本改编为我试图解决的问题。但是,我没有做到这一点,在尝试使用 JuMP 解决以下 LMI 问题时,一直导致 Julia 崩溃或给我一个 WhosebugError
,
Lyapunov stability of hybrid systems
对于未知对称矩阵Ui,Wi、M(和Pi),其中Ui,Wi大于或等于0。
我已经 运行 逐行设置模型,发现它在我定义最后一个约束时发生:
@SDconstraint(model, Ai'*Pi + Pi*Ai + Ei'*Ui*Ei ⪯ 0)
我在 JuMP 上看到了其他 post 的相关信息,它们通常是变量定义不当的结果,但这些知识并没有让我弄清楚这个问题。这是产生错误的代码的最小形式,
using JuMP
using MosekTools
using LinearAlgebra
E1p = [1 0; 0 1];
A1 = [-1 10; -100 -1];
Ei = E1p
Ai = A1
n=2
model = Model(Mosek.Optimizer)
@variable(model, Pi[i=1:n, j=1:n], Symmetric)
@variable(model, Ui[i=1:n, j=1:n], Symmetric)
@variable(model, Wi[i=1:n, j=1:n], Symmetric)
@variable(model, M[i=1:n, j=1:n], Symmetric)
@SDconstraint(model, Ui ⪰ 0)
@SDconstraint(model, Wi ⪰ 0)
@SDconstraint(model, Ei'*M*Ei ⪰ 0)
@SDconstraint(model, Pi - Ei'*Wi*Ei ⪰ 0)
@SDconstraint(model, Ai'*Pi + Pi*Ai + Ei'*Ui*Ei ⪯ 0)
optimize!(model)
Pi = value.(Pi)
重复很长时间的错误的第一行是:
in mutable_operate! at MutableArithmetics/bPWR4/src/linear_algebra.jl:132
这看起来像是一个错误。
这是一个 MWE
using JuMP
model = Model()
E = [1 1; 1 1]
@variable(model, P[1:size(E, 1), 1:size(E, 2)], Symmetric)
julia> @SDconstraint(model, E' * P + P * E + E' * P * E <= 0)
ERROR: WhosebugError:
Stacktrace:
[1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64}) (repeats 79984 times)
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/linear_algebra.jl:132
如果我放弃第一个学期,我会得到
julia> @SDconstraint(model, P * E + E' * P * E <= 0)
ERROR: MethodError: no method matching *(::Matrix{AffExpr})
Closest candidates are:
*(::StridedMatrix{T} where T, ::LinearAlgebra.AbstractQ) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:680
*(::StridedMatrix{T} where T, ::LinearAlgebra.Adjoint{var"#s832", var"#s831"} where {var"#s832", var"#s831"<:LinearAlgebra.AbstractQ}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:720
*(::StridedVecOrMat{T} where T, ::LinearAlgebra.Adjoint{var"#s832", var"#s831"} where {var"#s832", var"#s831"<:LinearAlgebra.LQPackedQ}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lq.jl:254
...
Stacktrace:
[1] mutable_operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::Matrix{AffExpr}) (repeats 2 times)
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/linear_algebra.jl:132
[2] operate_fallback!(::MutableArithmetics.IsMutable, ::Function, ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/interface.jl:330
[3] operate!(::typeof(MutableArithmetics.sub_mul), ::Matrix{AffExpr}, ::LinearAlgebra.Adjoint{Int64, Matrix{Int64}}, ::LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}, ::Matrix{Int64})
@ MutableArithmetics ~/.julia/packages/MutableArithmetics/bPWR4/src/rewrite.jl:80
[4] macro expansion
@ ~/.julia/packages/MutableArithmetics/bPWR4/src/rewrite.jl:276 [inlined]
[5] macro expansion
@ ~/.julia/packages/JuMP/y5vgk/src/macros.jl:447 [inlined]
[6] top-level scope
@ REPL[15]:1
所以这可能与https://github.com/jump-dev/MutableArithmetics.jl/issues/84有关。
我会开一个问题。 (编辑:https://github.com/jump-dev/MutableArithmetics.jl/issues/86)
作为解决方法,首先将其拆分为一个表达式:
julia> @expression(model, ex, E' * P + P * E + E' * P * E)
2×2 Matrix{AffExpr}:
3 P[1,1] + 4 P[1,2] + P[2,2] 4 P[1,2] + 2 P[2,2] + 2 P[1,1]
2 P[1,1] + 4 P[1,2] + 2 P[2,2] 4 P[1,2] + 3 P[2,2] + P[1,1]
julia> @SDconstraint(model, ex <= 0)
[-3 P[1,1] - 4 P[1,2] - P[2,2] -2 P[1,1] - 4 P[1,2] - 2 P[2,2];
-2 P[1,1] - 4 P[1,2] - 2 P[2,2] -P[1,1] - 4 P[1,2] - 3 P[2,2]] ∈ PSDCone()