I’m trying to solve a two-step optimization problem (to fit a complex econometric model) that requires solving a [fairly straight forward] quadratic program (as a function of the parameters of the model) at every iteration.
Is there any issue with initializing a new JuMP model to solve the quadratic program at each step? Would a different approach (perhaps a direct quad-prog solver or something) be preferable?
A simple example of what I want to do (for intuition) might be something like this:
# Quad Prog Function function solveQadProg( H, f, ceq, beq) m = Model(solver=IpoptSolver(print_level=0, tol=1e-12)) @variables m begin b[i=1:length(f)] >= 0 end @objective(m, Min, sum(H[i]*b[i]^2 + f[i]*b[i] for i=1:length(f) )) @constraint(m, sum( b[i] * ceq[i] for i=1:length(ceq) ) == beq) status = JuMP.solve(m) b_min = getvalue(b) return b_min end ### Main Model #user M_main = Model(solver=IpoptSolver()) @variables M_main begin f_var[i=1:N] beq_var end ## This would require solveQadProg to be registered properly, etc. and wouldn't work as is, but it gets the idea across @expression(M_main, b_quadprog[1:N], @solveQadProgm, [H_data, f_var, ceq_data, beq]) @objective(M_main, sum(someData[i] - b_quadprog[i] for i=1:N)^2 / N) status = JuMP,solve(M_main)