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

Specific Domains Optimization (Mathematical)
,
1

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?

2

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
3

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