# Optimization with JuMP: an affine expression enters a nonlinear constraint (dynamic discrete choice)

I am trying to use JuMP to solve a dynamic programming problem, but I cannot get MPEC (with Ipopt) to work.

I am solving a prediction problem. There are 2 products, “safe” and “risky”. The safe product always yields utility 0 and the risky product yields utility `β x + ξ`, where `x` is a (time-varying) observable characteristic, `β` a regression coefficient and `ξ` a (time-varying) utility shock. In each period, if the decision maker opts for the safe product, they get utility 0 and exit the market. If the decision maker opts for the risky one, then with probability `q` (time-varying) they get it and exit the market, and with probability `1-q` they fail and advance to the next period.

I want to minimize the sum of squared utility shocks `ξ`, subject to the constraint that the model-implied market shares `s(β, ξ)` match the data market shares `s0`. Here is a minimal (non-working) example:

``````using JuMP, Ipopt, Random
# generate fake data
N = 100  # number of periods
β = 1.0  # fundamental parameter
x = randn(N)  # exog. regressor
ξ0 = 0.25 * randn(N)  # good-specific iid shock; in fact E[ξ] ≜ 0
u = x * β .+ ξ0  # mean utility of good 1
q = rand(N)  # probability of winning when choosing good 1 (good 2 always has q = 1, u = 0)
v, V = zeros(N), zeros(N)  # initialize choice-specific value function (1) and value function
for i=1:N
v = q .* u + (1 .- q) .* vcat(V[2:end], 0.0)  # solve cvf(1) for all t
V = log.(1 .+ exp.(v))  # solve emax for all t
end

s0 = exp.(v) ./(1 .+ exp.(v))  # generate observed market shares

# code up MPEC
function obj(ξ...)
sum(ξ[i] .^ 2 for i=1:N)
end

function getMeanUtilities(β, ξ::Vararg{TT, NN} where NN) where {TT <: Union{Real, VariableRef}}
# recover mean utilities
u = Vector{Any}(undef, N)

for i=1:N
u[i] = x[i] * β + ξ[i]
end
u
end

function getShares(u::Vararg{TT, NN} where NN) where {TT <: Union{Real, VariableRef}}

# # initialize and recover cvf
v = Vector{TT}(undef, N)
V = Vector{TT}(undef, N)
s = Vector{TT}(undef, N)

for i=1:N
v[i] = q[i] * u[i]
V[i] = log(1 + exp(v[i]))
end

for t=1:N
for i=1:N
v[i] = q[i] * u[i] + (1-q[i]) * (i == N ? 0 : V[i+1])
V[i] = log(1 + exp(v[i]))
end
end

for i=1:N
s[i] = exp(v[i])/(1 + exp(v[i]))
end
s
end

# now define model
mdl = Model(Ipopt.Optimizer)
@variables(mdl, begin
β
ξ[1:N]
u[1:N]
end
)
register(mdl, :obj, N, obj; autodiff=true)
register(mdl, :getMeanUtilities, N+1, getMeanUtilities;
autodiff=true)
register(mdl, :getShares, N, getShares; autodiff=true)
@constraint(mdl, [i=1:N], getMeanUtilities(β, ξ...)[i] == u[i])
@NLconstraint(mdl, [i=1:N], getShares(u...)[i] == s0[i])
``````

The error I get when trying to define the nonlinear constraint is the following:

``````MethodError: Cannot `convert` an object of type Array{GenericAffExpr{Float64,VariableRef},1} to an object of type VariableRef
``````

How should I proceed?

User-defined functions cannot return a vector output.

This example suggests a work-around: https://jump.dev/JuMP.jl/stable/examples/nlp_tricks/#User-defined-functions-with-vector-outputs.

If you have lots of functions, you will need to use the raw-expression syntax: Registering an array of functions and using them in constraints and objectives - #2 by odow (more docs: https://jump.dev/JuMP.jl/stable/nlp/#Raw-expression-input).

2 Likes

Thanks @odow. Here’s a now-working example, cribbing the `Memoize` function from https://jump.dev/JuMP.jl/stable/examples/nlp_tricks/#User-defined-functions-with-vector-outputs:

``````using JuMP, Ipopt, Random
# generate fake data
N = 100  # number of periods
β0 = 1.0  # regression parameter
x = randn(MersenneTwister(1234), N)  # exog. regressor
ξ0 = 0.25 * randn(MersenneTwister(2345), N)  # good-specific iid shock; in fact E[ξ] ≜ 0
u0 = x * β0 .+ ξ0  # mean utility of good 1
q = rand(MersenneTwister(3456), N)  # probability of winning when choosing good 1 (good 2 always has q = 1, u = 0)
v, V = zeros(N), zeros(N)  # initialize cvf(1), emax(1)
for i=1:N
v = q .* u0 + (1 .- q) .* vcat(V[2:end], 0.0)  # solve cvf(1) for all t
V = log.(1 .+ exp.(v))  # solve emax for all t
end

s0 = exp.(v) ./(1 .+ exp.(v))  # generate "observed" market shares

# code up MPEC
function obj(ξ...)
sum(ξ[i] .^ 2 for i=1:N)
end

function getMeanUtilities(β, ξ::Vararg{TT, NN} where NN) where {TT <: Union{Real, VariableRef}}
# recover mean utilities
u = Vector{TT}(undef, N)

for i=1:N
u[i] = x[i] * β + ξ[i]
end
u
end

function getShares(u::Vararg{TT, NN} where NN) where {TT <: Union{Real, VariableRef}}

# # initialize and recover cvf
v = Vector{TT}(undef, N)
V = Vector{TT}(undef, N)
s = Vector{TT}(undef, N)

for i=1:N
v[i] = q[i] * u[i]
V[i] = log(1 + exp(v[i]))
end

tol = 1e-14
dist = 1e3
v_ = similar(v)

while dist > tol
v_ = copy(v)
for i=1:N
v[i] = q[i] * u[i] + (1-q[i]) * (i == N ? 0.0 : V[i+1])
V[i] = log(1 + exp(v[i]))
end
dist = norm(v_ - v)
end

for i=1:N
s[i] = exp(v[i])/(1 + exp(v[i]))
end
s
end

# now define model
"""
memoize(foo::Function, n_outputs::Int)

Take a function `foo` and return a vector of length `n_outputs`, where each
element is a function that returns the `i`'th output of `foo`.

To avoid duplication of work, cache the most-recent evaluations of `foo`.
Because `foo_i` is auto-differentiated with ForwardDiff, our cache needs to
work when `x` is a `Float64` and a `ForwardDiff.Dual`.
"""
function memoize(foo::Function, n_outputs::Int)
last_x, last_f = nothing, nothing
last_dx, last_dfdx = nothing, nothing
function foo_i(i, x::T...) where {T <: Real}
if T == Float64
if x != last_x
last_x, last_f = x, foo(x...)
end
return last_f[i]::T
else
if x != last_dx
last_dx, last_dfdx = x, foo(x...)
end
return last_dfdx[i]::T
end
end
return [(x...) -> foo_i(i, x...) for i = 1:n_outputs]
end

# memoize mean utility and share functions
memoized_u = memoize(getMeanUtilities, N)
memoized_s = memoize(getShares, N)

# run model
function runModel(u_functions, s_functions)
model = Model(Ipopt.Optimizer)
@variables(model, begin
β
ξ[1:N]
u[1:N]
end
)
for (i, f) in enumerate(u_functions)
f_sym = Symbol("f_\$(i)")
register(model, f_sym, N+1, f; autodiff = true)
add_NL_constraint(model, :(\$(u[i]) == \$(f_sym)(\$(β), \$(ξ...))))
end

for (i, s) in enumerate(s_functions)
s_sym = Symbol("s_\$(i)")
register(model, s_sym, N, s; autodiff = true)
add_NL_constraint(model, :(\$(s0[i]) == \$(s_sym)(\$(u...))))
end

@constraint(model, β <= 3.0)
@constraint(model, β >= -3.0)

register(model, :obj, N, obj; autodiff=true)
@NLobjective(model, Min, obj(ξ...))
set_optimizer_attribute(model, "max_cpu_time", 420.0)

optimize!(model)
@show value(β)
# return model
end

b = runModel(memoized_u, memoized_s)

``````
1 Like