Question about Nonlinear expressions in jump

Specific Domains Optimization (Mathematical)
1

Hi, I am trying to use Jump (v 1.23.2) with Julia (1.11.0) with nonlinear constraints where I essential invert a constructed matrix depending on the variable x.
Namely, something like

m = Model(optimizer_with_attributes(Ipopt.Optimizer))
@variable(m, x[1:4*N])
@NLobjective(m, Min, sum(x)^2)
aVal(x::vector) = [x[1], x[2]; x[3], x[4]] # construct 2x2 matrix depending on variable x. Simplified in this example. 
@constraint(m, eqconsc[i=1:N], blah = inv(aVal(x[(i-1)*4+1:4*i]))*x[(i-1)*4+1:4*i]) # compute blah = (A(x))^{-1} x

and it appears I am unable to have a nonlinear expression inside a nonlinear expression (ie inv, a nonlinear expression takes in another nonlinear expression aVal)). Is there anyway around this or am I making a mistake in the syntax?

2

Hi @lukebhan, welcome to the forum :smile:

There are a couple of things going on:

  1. You’re mixing the old @NL nonlinear syntax with the new
  2. aVal(x::vector) = [x[1], x[2]; x[3], x[4]] should be aVal(x::Vector) = [x[1] x[2]; x[3] x[4]]
  3. Is this the error you’re talking about?
julia> aVal(x[(i-1)*4+1:4*I])
2×2 Matrix{VariableRef}:
 x[1]  x[2]
 x[3]  x[4]

julia> inv(aVal(x[(i-1)*4+1:4*I]))
ERROR: MethodError: no method matching NonlinearExpr(::NonlinearExpr)

Closest candidates are:
  GenericNonlinearExpr{V}(::Symbol, ::Vector{Any}) where V<:AbstractVariableRef
   @ JuMP ~/.julia/packages/JuMP/6RAQ9/src/nlp_expr.jl:95
  GenericNonlinearExpr{V}(::Symbol, Any...) where V<:AbstractVariableRef
   @ JuMP ~/.julia/packages/JuMP/6RAQ9/src/nlp_expr.jl:84

Stacktrace:
 [1] oneunit(::Type{NonlinearExpr})
   @ Base ./number.jl:372
 [2] lutype(T::Type)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:210
 [3] lu(A::Matrix{VariableRef}, pivot::LinearAlgebra.RowMaximum; check::Bool)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:300
 [4] lu(A::Matrix{VariableRef}, pivot::LinearAlgebra.RowMaximum)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:299
 [5] lu(A::Matrix{VariableRef})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:299
 [6] inv(A::Matrix{VariableRef})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/dense.jl:917
 [7] top-level scope
   @ REPL[22]:1

julia> oneunit(NonlinearExpr)
ERROR: MethodError: no method matching NonlinearExpr(::NonlinearExpr)

Closest candidates are:
  GenericNonlinearExpr{V}(::Symbol, ::Vector{Any}) where V<:AbstractVariableRef
   @ JuMP ~/.julia/packages/JuMP/6RAQ9/src/nlp_expr.jl:95
  GenericNonlinearExpr{V}(::Symbol, Any...) where V<:AbstractVariableRef
   @ JuMP ~/.julia/packages/JuMP/6RAQ9/src/nlp_expr.jl:84

Stacktrace:
 [1] oneunit(::Type{NonlinearExpr})
   @ Base ./number.jl:372
 [2] top-level scope
   @ REPL[23]:1

This happens because you can’t call inv on a matrix of nonlinear expressions, but we don’t catch this particular error so bad things happen. I’ll see if we can fix this.

There’s an easy fix for your problem. Instead of writing inv(A) * x = blah, do A * blah = x:

N = 3
using JuMP
model = Model()
@variable(model, x[1:4N])
@objective(model, Min, sum(x)^2)
aVal(x::Vector) = [x[1] x[2]; x[3] x[4]]
@variable(model, blah[1:2, 1:N])
for i in 1:N
    @constraint(model, aVal(x[(i-1)*4+1:4*i]) * blah .== x[(i-1)*4+1:4*i])
end

Note that in most cases, whether you use JuMP or regular Julia, you almost never want to call inv(::Matrix). There is almost always a better way of solving the problem.

1 Like
3

Hi @odow
Thanks so much for your help and that is the exact error I am getting! So the issue is not that simple (sorry for not being clear). Namely, my expression is of the form x_{k+1} = A(x_k) * x_k + B(x_k) * x_k + C(x_k) where A, B have different types of inversion calls in them. Is there anyway to handle this? IE, my syntax really is something like:

for i in 1:N
    @constraint(model, x[i*4+1:4*i] = inv(aVal( x[(i-1)*4+1:4*i)])) *  x[(i-1)*4+1:4*i)] +inv(bVal( x[(i-1)*4+1:4*i)]))*  x[(i-1)*4+1:4*i)] + C(x[(i-1)*4+1:4*i)])])
end

I know this is a challenging expression, but I don’t think there is an easy way to fully write it out. Perhaps I can just write out the inverse operation myself since the dimension is static? Any other ideas?

1 Like
4

Add more temporary variables:

@variable(model, a_inv_x[1:2, 1:N])
@variable(model, b_inv_x[1:2, 1:N])
for i in 1:N
    x_i = x[4(i-1).+(1:4)]
    @constraint(model, aVal(x_i) * a_inv_x .== x_i)
    @constraint(model, bVal(x_i) * b_inv_x .== x_i)
    @constraint(model, x[4(i).+(1:4)] .== a_inv_x + b_inv_x + C(x_i))
end
1 Like
5

Ah, I see this is a nice trick, I did not realized. Thank you for your help!

1 Like
6

No worries. Let me know if you have other questions