Creating and multiplying a large matrix high order of JuMP variables

Specific Domains Optimization (Mathematical)
1

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?

1 Like
2

See the syntax notes from the documentation:

All expressions must be simple scalar operations. You cannot use dot, matrix-vector products, vector slices, etc. Translate vector operations into explicit sum() operations or use the AffExpr plus auxiliary variable trick described below.

There is no operator overloading provided to build up nonlinear expressions. For example, if x is a JuMP variable, the code 3x will return an AffExpr object that can be used inside of future expressions and linear constraints. However, the code sin(x) is an error. All nonlinear expressions must be inside of macros.

You could try expressing the elements of your matrix by using NLexpression.

3

Miles,

Thank you for the response.

I realize that I can create new @NLexpression() each time I carry out operations, like this:

@NLexpression(mdl, mult[j=1:N],prod(mat[i,j] for i in 1:N))

But, now I need to create several of these expressions with different names, I tried to do this several days ago and just ended up finding a workaround, but it would be useful to know how to do this. Here is my attempt, if I have an array of variable names as:

ss = [Symbol("xy$(int)") for int in 1:Ni];

then I want to assign these names programatically I tried:

x = [Array(Any,Nck[int],) for int in 1:Ni];
for int in 1:Ni
  for j in 1:Nck[int]
    x[int][j] = ts_JuMP[int,j];
  end
  N = length(x[int])
  mat = repmat(x[int],1,N);
  :($(ss[int])) = @NLexpression(mdl,[j=1:N],prod(mat[i,j] for i in 1:N))
end

but I get:

ERROR: MethodError: Cannot `convert` an object of type Array{JuMP.NonlinearExpression,1} to an object of type Symbol
This may have arisen from a call to the constructor Symbol(...),
since type constructors fall back to convert methods.
 in setindex!(::Array{Symbol,1}, ::Array{JuMP.NonlinearExpression,1}, ::Int64) at ./array.jl:415
 in macro expansion; at ./REPL[7]:7 [inlined]
 in anonymous at ./<missing>:?

How do I do this?

@miles.lubin I found this example where you explain how to make variables programatically:
http://julia-programming-language.2336112.n4.nabble.com/Defining-Variables-Programmatically-td6013.html

But I was wondering if there is a better way to do this? i.e. a way that will not lead to the unstable behavior that was talked about.

Also I was thinking that it would be clearer to define the variables in a format such as

x1 = @NLexpression(
x2 = @NLexpression(..
.
.
.
xn = @NLexpression(..

Is this possible? Does this fix the instability issues?
Thanks!

4

In order to solve my problem, I ended up creating many many @NLexpressions when I do a print(mdl) I have over 1000 subexpressions. Can this create an issue? Because, the model is building and then it fails as:

julia> include("Moon_Lander.jl")
WARNING: Method definition stateEquations(JuMP.Model, Array{Any, 2}, Array{Any, 2}, Int64) in module Main at /home/febbo/.julia/v0.5/NLOptControl/test/Moon_Lander.jl:18 overwritten at /home/febbo/.julia/v0.5/NLOptControl/test/Moon_Lander.jl:18.
This is Ipopt version 3.12.1, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:      519
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      530

Total number of variables............................:       48
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       34
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.63e+01 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
WARNING: Problem in step computation; switching to emergency mode.
   1r 0.0000000e+00 1.63e+01 9.99e+02   1.0 0.00e+00  20.0 0.00e+00 0.00e+00R  1
WARNING: Problem in step computation; switching to emergency mode.
Restoration phase is called at point that is almost feasible,
  with constraint violation 0.000000e+00. Abort.
Restoration phase in the restoration phase failed.

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   1.0000000000000000e+01    1.6251899999999999e+01
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.0000000000000000e+01    1.6251899999999999e+01


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 2
Total CPU secs in IPOPT (w/o function evaluations)   =      0.004
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Restoration Failed!
WARNING: Ipopt finished with status Restoration_Failed
WARNING: Not solved to optimality, status: Error

I seem to remember having a similar issue months back and then I think that when I reduced my @NLexpressions it was fixed. In my current case, it would be difficult to do this though.

Also:

  1. I am testing my @NLexpression() with data and they are giving the correct results
  2. I in any division operation, I added +eps() in the denominator to avoid /0

Any thoughts?

5

Unless there’s a bug that I’m not aware of, having 1000 subexpressions on its own is not any sort of issue in terms of correctness of the model.