JuMP; Rewriting an @eval @NLexpression with a user defined function in local scope w/ Expr

Specific Domains Optimization (Mathematical)
1

I’m trying to use a user-defined function that takes in several distinct arguments, and then a splatted N-dim vector of parameters as an @NLexpression to be used in a constraint of an optimization problem.

I described the problem in this thread.

I got the following to work:

function my_constraint_i(q_b_i, q_b_arr...)
    q_b= [q_b_arr...]
   ....
    STUFF INVOLVING q_b
    constriant_i_out = __function_of_q_b_vector__
   ....
   return constriant_i_out
end

JuMP.register(m, :my_constraint_i, (N+1), my_constraint_i, autodiff=true)
foc_constraint = @eval @NLexpression(m, [i=1:N], my_constraint_i(i, $(q_b...)))
@NLconstraint(m, [i=1:N, foc_constraint[i] == 0)

where my_constraint_i is a function that takes in an index i, and a splatted N-dim vector of parameters and outputs the i’th constraint (for constraining an optimizer w/ a a set of N nonlinear equations in N unknowns).

However, @eval appears to be defined in global scope, and I’m hoping to run the solver inside a function (this is an inner loop of a bigger optimization problem). Is it possible to rewrite this with Expr?

The following did not work :frowning:

foc_constraint = @NLexpression(m, [i=1:N], Expr(:call, :my_constraint_i, [i, (q_b[j] for j=1:N)...]...))
2

If you post runnable pieces of code, it will be easier for people to help.

2 Likes
3

Any line that contains @eval @NLexpression is pretty suspicious. If you need to do anything beyond the syntax that’s supported by @NLexpression then you should look at the raw expression input format. This format lets provide a julia expression object.

This should basically work (haven’t tested it):

for i in 1:N
    m.addNLconstraint(m, :($(Expr(:call, :my_constraint_i, [i, (q_b[j] for j=1:N)...]...)) == 0))
end
1 Like
4

Here’s a full example of this kind of thing:

https://github.com/tkoolen/MotionCaptureJointCalibration.jl/blob/987b9885e6fb83edf5ff61871b636c0b9bc21d31/src/solve.jl#L107

1 Like
5

Thank you guys very much. Your answers have helped me understand a bit more of what’s going on, but @miles.lubin’s suggestion doesn’t quite work and I’m not sure why (I’m getting a Method error saying “no method matching expr_to_nodedata”).

Below is a runnable example. The goal is to solve the set of nonlinear equations defined as in the pdf screenshot below (Note: the sum over h^k(s)/(1-H^k(s)) is just a number I’m calling H in the code)

image
image1506×338 36.6 KB

Thank you again for the help!

using JuMP
using Ipopt
using Parameters

q_a = [3.2, 1.6]
q_o = [2.5, 2.1]
b = [1000, 700]
c = [900, 800]
sigma_t = [0.5, 0.2]
H = [0.000004, 0.000004]

m = Model(solver=IpoptSolver(print_level=5, tol=1e-8))

@variables m begin
    q_b[i=1:length(q_a)]
end

@NLobjective(m, Min, sum((q_b[i] - q_a[i])^2 for i in 1:length(q_a)))

function foc_constraint_i(q_b_i, gamma, q_b_arr...)
    q_b= [q_b_arr...]

    q_b_i = round(Int64, q_b_i)

    b_min_c = b - c
    sigma_sq = sigma_t .* sigma_t

    scaling_factor = (q_o[q_b_i] / gamma) * H[q_b_i]

    variance_utility = exp(-(gamma^2 / 2) * ( (sigma_sq .* b_min_c)' * b_min_c))
    profit_utility = exp(gamma * (dot(q_b, b_min_c)))

    utility = (profit_utility * variance_utility) * scaling_factor

    lhs = q_b[q_b_i] - utility
    rhs = gamma * (sigma_sq[q_b_i] * b_min_c[q_b_i]) - scaling_factor

    constraint_i = lhs - rhs

    return constraint_i
end

foc_constraint_i(100, 0.000001, q_a...) ##Test that this function works

JuMP.register(m, :foc_constraint_i, (2+length(q_a)), foc_constraint_i, autodiff=true)

for i in 1:length(q_a)
    JuMP.addNLconstraint(m, :($(Expr(:call, :foc_constraint_i, [i, 0.00001, (q_b[j] for j=1:length(q_a))...]...)) == 0))
end

status = solve(m)

q_b_min = getvalue(q_b)
answer = getobjectivevalue(m)
6

I should note that I’ve (eventually) managed to rewrite this without using the function syntax, but I have more complex problems of this sort further into my project and so it would be really helpful to be able to get by with constraints/objectives defined as user-defined functions…

7

The following seems to work instead:

for i in 1:length(q_a)
    JuMP.addNLconstraint(m, :($(Expr(:call, :foc_constraint_i, i, 0.00001, (q_b[j] for j=1:length(q_a))...)) == 0))
end
2 Likes