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)
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)