Consider the following optimisation model:

I = 1:5; J = 1:5; K = 1:5;

m = Model(HiGHS.Optimizer)
@variable(m, A[i in I, j in J, k in K] >= 0)
@variable(m, y[i in I, k in K] >= 0)

@constraint(m, sum(collect(A[:,:,1]),dims = 2) .== y[:,1])

In the model, I have a 5x5x5 array variable A, and a 5x5 array variable y. I want to have an equation constraint, where the vector obtained by adding the column vectors of the matrix A[:,:,1] together is (entry-wise) equal to the column vector p[:,1].

However, the above code gives an error:

MethodError: no method matching _broadcasted_type(::Base.Broadcast.ArrayStyle{JuMP.Containers.DenseAxisArray}, ::Base.HasShape{2}, ::Type{AffExpr})

How can I fix that?

CodePudding user response：

Apply collect to the other side of the equation:

julia> @constraint(m, vec(sum(collect(A[:,:,1]),dims = 2)) .== collect(y[:,1]))
5-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 A[1,1,1]   A[1,2,1]   A[1,3,1]   A[1,4,1]   A[1,5,1] - y[1,1] == 0.0
 A[2,1,1]   A[2,2,1]   A[2,3,1]   A[2,4,1]   A[2,5,1] - y[2,1] == 0.0
....

EDIT

note that your JuMP experience will be much more enjoyable if you do:

I = 5; J = 5; K = 5;
@variable(m, A[1:I, 1:J, 1:K] >= 0)
@variable(m, y[1:I, 1:K] >= 0)

Now you do not have dense arrays and can do what you initially wanted:

julia> @constraint(m, sum(A[:,:,1],dims = 2) .== y[:,1])
5×1 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 A[1,1,1]   A[1,2,1]   A[1,3,1]   A[1,4,1]   A[1,5,1] - y[1,1] == 0.0
 A[2,1,1]   A[2,2,1]   A[2,3,1]   A[2,4,1]   A[2,5,1] - y[2,1] == 0.0
...