I am creating a package NLOptControl.jl that uses JuMP and when the final time is a design variable (`tf_var`

), the array DMatrix needs to be recalculated each iteration.

I modified the functionality that calculates the DMatrix so that it is a JuMP expression, but realized that the result will be a matrix with terms `tf_var^N`

where `N`

is the discretization mesh. Anyways, there may be no need to get bogged down with the specifics of why I am trying to do the so, let me provided a simple example:

```
using JuMP, Ipopt
Ni=2; Nck=[10,5];
τ = [zeros(Float64,Nck[int],) for int in 1:Ni];
τ[1] = [-1.0,-0.927484,-0.763842,-0.525646,-0.236234,0.0760592,0.380665,0.647767,0.851225,0.971175];
τ[2] = [-1.0,-0.72048,-0.167181,0.446314,0.885792];
mdl = Model(solver = IpoptSolver());
@variable(mdl, tf_var)
ts_JuMP = [Array(Any, Nck[int],) for int in 1:Ni];
for int in 1:Ni
ts_JuMP[int] = tf_var/2*τ[int] + tf_var/2;
end
x = ts_JuMP[1];
N = length(x)
mat = repmat(x,1,N);
```

which gives:

```
julia> mat
10×10 Array{Any,2}:
0 0 0 0 … 0 0 0
0.03625800000000001 tf_var 0.03625800000000001 tf_var 0.03625800000000001 tf_var 0.03625800000000001 tf_var 0.03625800000000001 tf_var 0.03625800000000001 tf_var 0.03625800000000001 tf_var
0.11807899999999999 tf_var 0.11807899999999999 tf_var 0.11807899999999999 tf_var 0.11807899999999999 tf_var 0.11807899999999999 tf_var 0.11807899999999999 tf_var 0.11807899999999999 tf_var
0.23717700000000003 tf_var 0.23717700000000003 tf_var 0.23717700000000003 tf_var 0.23717700000000003 tf_var 0.23717700000000003 tf_var 0.23717700000000003 tf_var 0.23717700000000003 tf_var
0.381883 tf_var 0.381883 tf_var 0.381883 tf_var 0.381883 tf_var 0.381883 tf_var 0.381883 tf_var 0.381883 tf_var
0.5380296 tf_var 0.5380296 tf_var 0.5380296 tf_var 0.5380296 tf_var … 0.5380296 tf_var 0.5380296 tf_var 0.5380296 tf_var
0.6903325 tf_var 0.6903325 tf_var 0.6903325 tf_var 0.6903325 tf_var 0.6903325 tf_var 0.6903325 tf_var 0.6903325 tf_var
0.8238835 tf_var 0.8238835 tf_var 0.8238835 tf_var 0.8238835 tf_var 0.8238835 tf_var 0.8238835 tf_var 0.8238835 tf_var
0.9256125 tf_var 0.9256125 tf_var 0.9256125 tf_var 0.9256125 tf_var 0.9256125 tf_var 0.9256125 tf_var 0.9256125 tf_var
0.9855875000000001 tf_var 0.9855875000000001 tf_var 0.9855875000000001 tf_var 0.9855875000000001 tf_var 0.9855875000000001 tf_var 0.9855875000000001 tf_var 0.9855875000000001 tf_var
```

Then I need to multiply all of the elements in the rows as

```
prod(mat,2)
```

But, I get the error:

```
julia> prod(mat,2)
ERROR: Cannot multiply a quadratic expression by an aff. expression
in *(::JuMP.GenericQuadExpr{Float64,JuMP.Variable}, ::JuMP.GenericAffExpr{Float64,JuMP.Variable}) at /home/febbo/.julia/v0.5/JuMP/src/operators.jl:217
in macro expansion at ./reduce.jl:104 [inlined]
in macro expansion at ./simdloop.jl:73 [inlined]
in mapreduce_impl(::Base.#identity, ::Base.#*, ::Array{Any,2}, ::Int64, ::Int64, ::Int64) at ./reduce.jl:102
in _mapreduce(::Base.#identity, ::Base.#*, ::Base.LinearFast, ::Array{Any,2}) at ./reduce.jl:162
in _reducedim_init(::Base.#identity, ::Base.#*, ::Base.#one, ::Base.#prod, ::Array{Any,2}, ::Int64) at ./reducedim.jl:111
in mapreducedim at ./reducedim.jl:253 [inlined]
in prod at ./reducedim.jl:303 [inlined]
in prod(::Array{Any,2}, ::Int64) at ./reducedim.jl:305
```

So, the above did not work and I tried:

```
@NLexpression(mdl, ts_JuMP[int=1:Ni,j=1:Nck[int]], tf_var/2*τ[int][j] + tf_var/2)
x = [Array(Any,Nck[int],) for int in 1:Ni];
int=1
for j in 1:Nck[int]
x[int][j] = ts_JuMP[int,j];
end
N = length(x[int])
mat = repmat(x[int],1,N);
```

but, I get this error:

```
julia> prod(mat,2)
ERROR: MethodError: no method matching *(::JuMP.NonlinearExpression, ::JuMP.NonlinearExpression)
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...) at operators.jl:138
*{T<:Union{JuMP.AbstractJuMPScalar,JuMP.GenericNormExpr{2,Float64,JuMP.Variable},JuMP.GenericNorm{P,Float64,JuMP.Variable},JuMP.NonlinearExpression},S}(::Array{S,N}, ::T<:Union{JuMP.AbstractJuMPScalar,JuMP.GenericNormExpr{2,Float64,JuMP.Variable},JuMP.GenericNorm{P,Float64,JuMP.Variable},JuMP.NonlinearExpression}) at /home/febbo/.julia/v0.5/JuMP/src/operators.jl:552
*(::SparseMatrixCSC{Tv,Ti<:Integer}, ::Union{JuMP.AbstractJuMPScalar,JuMP.GenericNormExpr{2,Float64,JuMP.Variable},JuMP.GenericNorm{P,Float64,JuMP.Variable},JuMP.NonlinearExpression}) at /home/febbo/.julia/v0.5/JuMP/src/operators.jl:567
...
in mapreduce_impl(::Base.#identity, ::Base.#*, ::Array{Any,2}, ::Int64, ::Int64, ::Int64) at ./reduce.jl:101
in _mapreduce(::Base.#identity, ::Base.#*, ::Base.LinearFast, ::Array{Any,2}) at ./reduce.jl:162
in _reducedim_init(::Base.#identity, ::Base.#*, ::Base.#one, ::Base.#prod, ::Array{Any,2}, ::Int64) at ./reducedim.jl:111
in mapreducedim at ./reducedim.jl:253 [inlined]
in prod at ./reducedim.jl:303 [inlined]
in prod(::Array{Any,2}, ::Int64) at ./reducedim.jl:305
```

So, any ideas how I can create an array of matrices that include high order JuMP variables?