Creating a macro for defining differential equation constraints in NLOptControl.jl

Specific Domains Optimization (Mathematical)
1

@miles.lubin please correct me if I am wrong, but from what I understand from the docs and our previous discussions, JuMP.jl does not allow the user to programmatically create variables and use ReverseDiffSparse.jl.

To improve the readability and reduce boiler-plate code, I have been working on a macro to allow the user to define the Differential Equation constraints easily. Without this macro, the code looks like this. Where the only thing that the user needs to modify are the expressions inside the @NLexpression() .

I can use the macro for simple linear problems like in this example, but if I try to use it with the macro, then I cannot solve nonlinear problems like the first example.

It fails because I need to create the nonlinear expression inside of the @NLexpression macro. I have tried many things but I cannot seem to get it to work. Am I missing something? Any suggestions?

Thanks!

2

It’s hard to follow the links to understand what you’re trying to do. Could you give a self-contained example?

3

@miles.lubin sorry for the delay, I was on a camping trip.

Basically, what I am trying to do is create a macro that can create a function for the user. Currently, the user would have to write the function as:

function Brachistochrone{T<:Any}(n::NLOpt,x::Array{T,2},u::Array{T,2})
  if n.s.integrationMethod==:tm; L=size(x)[1]; else L=size(x)[1]-1; end
  dx = Array(Any,L,n.numStates)
  dx[:,1]=@NLexpression(n.mdl, [j=1:L], x[j,3]*sin(u[j,1]) );
  dx[:,2]=@NLexpression(n.mdl, [j=1:L], x[j,3]*cos(u[j,1]) );
  dx[:,3]=@NLexpression(n.mdl, [j=1:L], 9.81*cos(u[j,1]) );
  return dx
end

where most of the above is boiler-plate code and the only thing that changes is the third argument that is in the @NLexpression() macro. For instance, if the user was to use my macro, the above function could be defined as:

@DiffEq(Brachistochrone,[x[j,3]*sin(u[j,1]);x[j,3]*cos(u[j,1]);9.81*cos(u[j,1])])

Do, you understand what I am trying to do now? If not, I will try to come up with a more contained example.

Thanks!

4

In principle you could “just” define

@DiffEq(Brachistochrone,[x[j,3]*sin(u[j,1]);x[j,3]*cos(u[j,1]);9.81*cos(u[j,1])])

to generate the code

function Brachistochrone{T<:Any}(n::NLOpt,x::Array{T,2},u::Array{T,2})
  if n.s.integrationMethod==:tm; L=size(x)[1]; else L=size(x)[1]-1; end
  dx = Array(Any,L,n.numStates)
  dx[:,1]=@NLexpression(n.mdl, [j=1:L], x[j,3]*sin(u[j,1]) );
  dx[:,2]=@NLexpression(n.mdl, [j=1:L], x[j,3]*cos(u[j,1]) );
  dx[:,3]=@NLexpression(n.mdl, [j=1:L], 9.81*cos(u[j,1]) );
  return dx
end

JuMP doesn’t currently have any ability to make multidimensional output easy to generate.

5

Is the @DiffEq macro capable of defining differential equation constraints, that get automatically discretized to algaebric constraints? Does it work for linear and nonlinear programs/ODEs? What about PDEs? If so, is there a similar macro for functionals that automatically approximates the integral to an algaebric expression?

6

@mohamed82008 the @DiffEq macro is currently under development. Eventually it will automatically approximate the integral to nonlinear ODEs. There is no plan to do PDEs yet.

2 Likes
7

@miles.lubin I am still struggling to figure this out. The issues is not generating the function above, it happens when I try to actually use it. I get the error:

ERROR: sin is not defined for type Variable. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.
 in macro expansion at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:226 [inlined]
 in macro expansion at /home/febbo/.julia/v0.5/JuMP/src/macros.jl:1193 [inlined]
 in Brachistochrone(::NLOptControl.NLOpt, ::Array{JuMP.Variable,2}, ::Array{JuMP.Variable,2}) at /home/febbo/.julia/v0.5/NLOptControl/src/macros.jl:20
 in OCPdef!(::NLOptControl.NLOpt) at /home/febbo/.julia/v0.5/NLOptControl/src/setup.jl:256
 in #configure!#51(::Array{Any,1}, ::Function, ::NLOptControl.NLOpt) at /home/febbo/.julia/v0.5/NLOptControl/src/setup.jl:379
 in (::NLOptControl.#kw##configure!)(::Array{Any,1}, ::NLOptControl.#configure!, ::NLOptControl.NLOpt) at ./<missing>:0
 in anonymous at ./<missing>:?

I think that this is happening because I am evaluating the expression (eq) that I pass the macro inside the function that I am creating. So, then when I run it it is like I am trying to run:

julia> n.r.x[1,3]*sin(n.r.u[1,1])
ERROR: sin is not defined for type Variable. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.
in sin(::JuMP.Variable) at /home/febbo/.julia/v0.5/JuMP/src/operators.jl:630

So, I am trying to define a nonlinear expression outside of the @NLexpression macro, which does not work. Do you have any ideas how I can do this? I have several other ways, but since nothing is working, I won’t go into it.

8

Don’t try to evaluate the nonlinear expression. You’ll have to generate a call to the macro @NLexpression and put the expression in the body of that call.

9

@miles.lubin thank you for your help with this, unfortunately I still cannot seem to get it running, but I have created a self-contained example that should be easier to work with:

using JuMP
type R
  eq
  mdl::JuMP.Model
end
function R()
  R(Any,JuMP.Model())
end

macro defODE(r,eq)
  r=esc(r);
  quote
      r.eq=$(Meta.quot(eq))
      function $(esc(:DE))(r,x)
         @NLexpression($r.mdl,$r.eq)
      end
  end
end

r=R();
@defODE(r,2+sin(x));
@variable(r.mdl,x);
expr=DE(r,x)

Where the output from the above code is:

ERROR: MethodError: no method matching parseNLExpr_runtime(::JuMP.Model, ::Expr, ::Array{ReverseDiffSparse.NodeData,1}, ::Int64, ::Array{Float64,1})
Closest candidates are:
  parseNLExpr_runtime(::JuMP.Model, ::Number, ::Any, ::Any, ::Any) at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:196
  parseNLExpr_runtime(::JuMP.Model, ::JuMP.Variable, ::Any, ::Any, ::Any) at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:202
  parseNLExpr_runtime(::JuMP.Model, ::JuMP.NonlinearExpression, ::Any, ::Any, ::Any) at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:208
  ...
 in macro expansion at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:226 [inlined]
 in macro expansion at /home/febbo/.julia/v0.5/JuMP/src/macros.jl:1193 [inlined]
 in DE(::R, ::JuMP.Variable) at ./REPL[166]:6
10

@miles.lubin I know that you are probably very busy with the meetup, but I would greatly appreciate any help that you can spare with this, I finally got it working for for the scalar case:

  function test(n,exp,x)
    @eval begin
      x1=$x[1,1]
      x1=eval(x1);
      eq=@NLexpression($n,$exp)
    end
    return eq
  end

the above can be tested with this made-up example:

  using JuMP, Ipopt
  n = Model(solver=IpoptSolver(print_level=0))
  @variable(n,x[1:4,1:4])
  exp=test(n,:(sin(x1)),x)

The issue that I am having is that I am dealing with arrays of expressions, so I tried to extend the above to an N-d case:

  function test(n,exp_arr,x)
    code=quote
      x1=$x[1,1]
      x1=eval(x1);
      @NLexpression($n,$exp_arr[1])
      end
    return eval(code)
  end

which can be tested with:

  using JuMP, Ipopt
  n = Model(solver=IpoptSolver(print_level=0))
  @variable(n,x[1:4,1:4])
  exp_arr=[:(sin(x1)),:(sin(x1))]
  exp=test(n,exp_arr,x)

which fails with:

julia> exp=test(n,exp_arr,x)
ERROR: MethodError: no method matching parseNLExpr_runtime(::JuMP.Model, ::Expr, ::Array{ReverseDiffSparse.NodeData,1}, ::Int64, ::Array{Float64,1})
Closest candidates are:
  parseNLExpr_runtime(::JuMP.Model, ::Number, ::Any, ::Any, ::Any) at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:196
  parseNLExpr_runtime(::JuMP.Model, ::JuMP.Variable, ::Any, ::Any, ::Any) at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:202
  parseNLExpr_runtime(::JuMP.Model, ::JuMP.NonlinearExpression, ::Any, ::Any, ::Any) at /home/febbo/.julia/v0.5/JuMP/src/parsenlp.jl:208
  ...
 in eval(::Module, ::Any) at ./boot.jl:234
 in test(::JuMP.Model, ::Array{Expr,1}, ::Array{JuMP.Variable,2}) at ./REPL[9]:7

The difference is that between an array indexed at a spot and a the item at that spot in a quoted expression can be seen with this simple example:

julia> A=[:(2+3),:(4),:(9-8)];
julia> t=A[1];
julia> eval(quote 
       @show isequal($A[1],$t)
       @show $A[1]
       @show $t
       end)
isequal((Any[:(2 + 3),4,:(9 - 8)])[1],2 + 3) = false
(Any[:(2 + 3),4,:(9 - 8)])[1] = :(2 + 3)
2 + 3 = 5
5

I have tried all kinds of things (which I will not go into), and I cannot get things working. Again, any help would be very nice. Thanks

11

@miles.lubin I ended up just calling the function for each expression that I have, and it works. But it would be more efficient if I could just call the function once

12

This is not a JuMP-specific issue, it’s a question about code generation and macro hygiene.