Array disjunctions

[timholy]

I am hoping to perform optimization on two matrices MW and MH constrained by their product, MW*MH = D where D is diagonal. If I define rtD = sqrt.(sqrt.(diag(D))) then the following work:

    @variables(model, begin
        -rtD[i]*rtD[j] <= MW[i=1:r, j=1:r] <= rtD[i]*rtD[j]
        -rtD[i]*rtD[j] <= MH[i=1:r, j=1:r] <= rtD[i]*rtD[j]
    end)
    @constraints(model, begin
        # McCormick envelopes for MW * MH .== D; these might need splitting for the disjunction
        [i=1:r, j=1:r, k=1:r], sum(-rtD[i]*rtD[k]*MH[k,j] - MW[i,k]*rtD[k]*rtD[j] - rtD[i]*rtD[j]*rtD[k]^2 for k=1:r) <= D[i,j]
        [i=1:r, j=1:r, k=1:r], sum( rtD[i]*rtD[k]*MH[k,j] + MW[i,k]*rtD[k]*rtD[j] - rtD[i]*rtD[j]*rtD[k]^2 for k=1:r) <= D[i,j]
        [i=1:r, j=1:r, k=1:r], sum( rtD[i]*rtD[k]*MH[k,j] - MW[i,k]*rtD[k]*rtD[j] + rtD[i]*rtD[j]*rtD[k]^2 for k=1:r) >= D[i,j]
        [i=1:r, j=1:r, k=1:r], sum( MW[i,k]*rtD[k]*rtD[j] - rtD[i]*rtD[k]*MH[k,j] + rtD[i]*rtD[j]*rtD[k]^2 for k=1:r) >= D[i,j]
    end)

but adding

    @disjunction(model,
        -rtD[i]*rtD[j] <= MW[i=1:r, j=1:r] <= 0,
                     0 <= MW[i=1:r, j=1:r] <= rtD[i]*rtD[j],
        reformulation = :big_m,
        name = :MW_disjunction)

results in

julia> include("testdata.jl")
ERROR: LoadError: UndefVarError: `i` not defined
Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/JuMP/KiENn/src/macros.jl:2246 [inlined]
 [2] macro expansion
   @ ~/.julia/packages/DisjunctiveProgramming/RK89U/src/macros.jl:69 [inlined]
 [3] envelope_wsparse(U::Matrix{Float64}, D::Matrix{Float64})
   @ Main ~/data/DaeWoo/SMF/opt.jl:55
 [4] envelope_wsparse(U::Matrix{Float64}, D::Diagonal{Float64, Vector{Float64}})
   @ Main ~/data/DaeWoo/SMF/opt.jl:69
 [5] top-level scope
   @ ~/data/DaeWoo/SMF/testdata.jl:14
 [6] include(fname::String)
   @ Base.MainInclude ./client.jl:478
 [7] top-level scope
   @ REPL[13]:1
in expression starting at /home/tim/data/DaeWoo/SMF/testdata.jl:14

caused by: Indexing an `Array` with keyword arguments is not supported.

Indexing with keyword arguments _is_ supported for `Container.DenseAxisArray`
(and `Container.SparseAxisArray`).

Force the container type by passing `container = DenseAxisArray` to any of
the JuMP macros. For example:

julia> model = Model();

julia> @variable(model, x[i=1:3], container = DenseAxisArray)
1-dimensional DenseAxisArray{VariableRef,1,...} with index sets:
    Dimension 1, Base.OneTo(3)
And data, a 3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]

julia> x[i=2]
x[2]

Stacktrace:
 [1] getindex(x::Matrix{VariableRef}; kwargs::Base.Pairs{Symbol, UnitRange{Int64}, Tuple{Symbol, Symbol}, NamedTuple{(:i, :j), Tuple{UnitRange{Int64}, UnitRange{Int64}}}})
   @ JuMP ~/.julia/packages/JuMP/KiENn/src/JuMP.jl:1038
 [2] macro expansion
   @ ~/.julia/packages/MutableArithmetics/Zy1nV/src/rewrite.jl:322 [inlined]
 [3] macro expansion
   @ ~/.julia/packages/JuMP/KiENn/src/macros.jl:508 [inlined]
 [4] macro expansion
   @ ~/.julia/packages/DisjunctiveProgramming/RK89U/src/macros.jl:66 [inlined]
...

[hdavid16]

Hi @timholy, are you trying to define the disjunction for every i, j? If so, you currently need to add the disjunction inside a loop for each index. For your reference, we are working on a rewrite that will allow indexing disjunctions as is done with the constraint macro in JuMP.

[hdavid16]

@timholy , the new version 0.4.0 will allow indexing in the @disjunct macro (to create arrays of disjunctions). Let me know if there is still interest in this and I can provide some sample code.