Plotting sensitivity of variable to parameter

Specific Domains Optimization (Mathematical)
1

Hi! I am using JuMP to model bertran competition in a market with firms and households.

One of the parameters in the model is the marginal cost for firms. I would like to evaluate how the average price of goods varies as this parameter changes.

So my question is whether it is possible to do this natively in JuMP? i.e. to evaluate the sensitivity of model outputs to changes in parameters. Or if the best way is just to put it in a giant for loop and evaluate it that way.

I am relatively new to Julia (1 week) and happy to be told if I am approaching this problem the wrong way. Many thanks!

2

Hi @laurence_jj, welcome to the forum (and to JuMP) :smile:

It’s probably simplest to put in a loop.

Can you share a reproducible example of your model? What solver? Where does the parameter appear?

3

Hi, I am justing the Ipopt solver. The parameter of interest I was trying to examine was μᴹᶜ, here is the code:

using JuMP
using Distributions
using Ipopt
using QuadGK
using Suppressor

# Model parameters
num_firms = 40
θ = 2.44  # relationship utility

μᴹᶜ = 1.62     # average marginal cost
σᴹᶜ = 0.73     # marginal cost variance
μᵣᴹᶜ = 0.0     # incumbent MC advantage

αᵘ = 1.14      # rate sensitivity
u₀ = 1.17

# Define distributions for consumer characteristics
χ_dist = Normal(0, 1)  # Distribution of consumer characteristics (χᵢ)

# Firm-specific base costs (mcⱼ)
incumbent_advantage = [fill(μᵣᴹᶜ, trunc(Int, num_firms/2)); fill(0, trunc(Int, num_firms/2))]


# Function to calculate expected marginal cost for firm j
function expected_MC(j)
    function MC_integrand(χ)
        # Function to calculate marginal cost for firm j and consumer type χ
        function calculate_MC(j, χ)
            mcj = μᴹᶜ + incumbent_advantage[j]  # Firm-specific component
            mci = 4.0 * cdf(Normal(0, 1), χ)  # Consumer-specific component
            return mcj + mci
        end

        return calculate_MC(j, χ) * pdf(χ_dist, χ)
    end
    expected_cost, _ = quadgk(MC_integrand, -4, 4, rtol=1e-4)
    return expected_cost
end

# Function to calculate market share
function market_share(prices, firm_idx)
    function integrand(χ)
        # Vectorized utility calculation
        utilities = -prices .* αᵘ .+ [k <= trunc(Int, num_firms/2) ? θ : 0.0 for k in 1:num_firms]
        
        # Logit choice probability
        exp_uj = exp(utilities[firm_idx])
        exp_sum = sum(exp.(utilities))+ exp(u₀)  # Include outside option in denominator
        choice_prob = exp_uj / exp_sum
        
        return choice_prob * pdf(χ_dist, χ)
    end
    
    share, _ = quadgk(integrand, -4, 4, rtol=1e-4)
    return share
end


# Register the market share function with JuMP
function market_share_jump(x...)
    prices = [value(x[i]) for i in 1:num_firms]
    firm_idx = round(Int, value(x[num_firms + 1]))
    return market_share(prices, firm_idx)
end


#################
#
#  create JuMP model
#
############################

# Create the model
model = Model(Ipopt.Optimizer)

# Decision variables: prices for each firm
@variable(model, p[1:num_firms])

# Initial price guesses using expected marginal costs
expected_MCs = [expected_MC(j) for j in 1:num_firms]
set_start_value.(p, expected_MCs .* 1.5)

# Register the nonlinear function
register(model, :market_share_jump, num_firms + 1, market_share_jump; autodiff=true)

# Constraints
@constraint(model, price_floor[j=1:num_firms], p[j] >= expected_MCs[j])
@constraint(model, price_ceiling[j=1:num_firms], p[j] <= 50.0)

# Profit expressions using actual market share calculation
@NLexpression(model, profit[j=1:num_firms], 
    (p[j] - expected_MCs[j]) * market_share_jump(p..., j))

# Total profit objective
@NLobjective(model, Max, sum(profit[j] for j in 1:num_firms))

# Optimizer settings
set_optimizer_attribute(model, "max_iter", 1000)
set_optimizer_attribute(model, "tol", 1e-4)
set_optimizer_attribute(model, "print_level", 0)

# Solve
optimize!(model)

#################
#
#  analyse results
#
############################


# Get results
optimal_prices = value.(p)

println("Optimization status: ", termination_status(model))
println("\nOptimal prices:")
println("Incumbents (firms 1-15):")
for j in 1:trunc(Int, num_firms/2)

    println("Firm $j: $(round(optimal_prices[j], digits=2)) (Expected MC: $(round(expected_MCs[j], digits=2)))")
end
println("\nEntrants (firms 16-30):")
for j in trunc(Int, (num_firms/2)+1):num_firms
    println("Firm $j: $(round(optimal_prices[j], digits=2)) (Expected MC: $(round(expected_MCs[j], digits=2)))")
end

# Calculate final market shares using optimal prices
final_shares = [market_share(optimal_prices, j) for j in 1:num_firms]

println("\nMarket shares:")
println("Average incumbent share: $(round(mean(final_shares[1:trunc(Int, num_firms/2)])*100, digits=2))%")
println("Average entrant share: $(round(mean(final_shares[trunc(Int, (num_firms/2)+1):num_firms])*100, digits=2))%")

Many thanks!

4

Here’s how I’d write your model:

julia> using JuMP

julia> import Distributions

julia> import Ipopt

julia> import QuadGK

julia> function main(;
           num_firms = 40,
           θ = 2.44,
           μᴹᶜ = 1.62,
           σᴹᶜ = 0.73,
           μᵣᴹᶜ = 0.0,
           αᵘ = 1.14,
           u₀ = 1.17,
       )
           num_firms_half = div(num_firms, 2)
           χ_dist = Distributions.Normal(0, 1)
           incumbent_advantage = [fill(μᵣᴹᶜ, num_firms_half); zeros(num_firms_half)]
           function expected_MC(j)
               expected_cost, _ = QuadGK.quadgk(-4, 4; rtol = 1e-4) do χ
                   mcj = μᴹᶜ + incumbent_advantage[j]
                   mci = 4.0 * Distributions.cdf(Distributions.Normal(0, 1), χ)
                   return (mcj + mci) * Distributions.pdf(χ_dist, χ)
               end
               return expected_cost
           end
           function market_share(prices, firm_idx)
               share, _ = QuadGK.quadgk(-4, 4; rtol = 1e-4) do χ
                   utilities = map(1:num_firms) do k
                       return -αᵘ * prices[k] + ifelse(k <= num_firms_half, θ, 0.0)
                   end
                   exp_uj = exp(utilities[firm_idx])
                   exp_sum = sum(exp, utilities) + exp(u₀)
                   choice_prob = exp_uj / exp_sum        
                   return choice_prob * Distributions.pdf(χ_dist, χ)
               end
               return share
           end
           market_share_jump(j, x...) = market_share(collect(x), round(Int, j))
           expected_MCs = expected_MC.(1:num_firms)
           model = Model(Ipopt.Optimizer)
           set_silent(model)
           @variable(
               model,
               expected_MCs[j] <= p[j in 1:num_firms] <= 50,
               start = 1.5 * expected_MCs[j],
           )
           @operator(model, op_market_share, 1 + num_firms, market_share_jump)
           @expression(
               model,
               profit[j in 1:num_firms], 
               (p[j] - expected_MCs[j]) * op_market_share(j, p...),
           )
           @objective(model, Max, sum(profit))
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return value.(p)
       end
main (generic function with 1 method)

julia> solutions = [x => main(; μᴹᶜ = x) for x in range(1.0, 2.0; length = 11)];

julia> map(solutions) do (x, y)
           return x => extrema(y)
       end
11-element Vector{Pair{Float64, Tuple{Float64, Float64}}}:
 1.0 => (4.35244388352189, 4.352443948268603)
 1.1 => (4.418130951995438, 4.418134349252176)
 1.2 => (4.485512143991204, 4.485512166093286)
 1.3 => (4.554574506199629, 4.554575449707015)
 1.4 => (4.625314534369957, 4.625314701417606)
 1.5 => (4.697712923216194, 4.697713471679339)
 1.6 => (4.771746944243104, 4.771752086201073)
 1.7 => (4.847386725269032, 4.847387100544168)
 1.8 => (4.9245956732298195, 4.924596132256306)
 1.9 => (5.003329940153706, 5.0033315602162824)
 2.0 => (5.08354131761681, 5.083541948581439)

I assume that I got things correct, but I may have made a typo somewhere, etc.

Note that @NL is now legacy syntax. There’s a new nonlinear interface: Nonlinear Modeling · JuMP

5

many thanks, I will check out the new interface!

6

I am finding the new interface won’t let me add the @expressions, which wasn’t an issue with the @NLexpressions:

┌ Warning: The model has been modified since the last call to `optimize!` (or `optimize!` has not been called yet). If you are iteratively querying solution information and modifying a model, query all the results first, then modify the model.
└ @ JuMP C:\XXXXXX
ERROR: OptimizeNotCalled()
Stacktrace:

Updated code:

   model = Model(Ipopt.Optimizer)

    # Decision variables: base prices for all firms
    #@variable(model, p[1:num_firms])
    @variable(model, p_ext[1:N_ent])
    @variable(model, p_int[1:N_inc])
    @variable(model, p_slope[1:N_inc])

    # Initial values
    set_start_value.(p_ext, 3.0 .+ randn(N_ent))  # Random values normally distributed around 5.0
    set_start_value.(p_int, 1 .* randn(N_inc))  # Random values normally distributed around 0.01
    set_start_value.(p_slope, 1 .* randn(N_inc))  # Random values normally distributed around 0.01

    # Bounds for slopes
    @constraint(model, slope_bounds[j=1:N_inc], 0.95 <= p_slope[j] <= 1.1)
    
    # Price bounds (for base prices)
    @constraint(model, price_ceiling_ent[j=1:N_ent], 0.01 <= p_ext[j] <= 10.0)
    @constraint(model, price_ceiling_int[j=1:N_inc], -1.0 <= p_int[j] <= 1.5)

    # Register the modified market share function
    register(model, :market_share_jump, N_inc + N_ent + N_inc + 
        1, market_share_jump; autodiff=true)
    register(model, :floor, 1, floor; autodiff = true)

    # First create the separate profit expressions
    # @NLexpression(model, profit_inc[j=1:N_inc], 
    #     ((p_int[j] + p_slope[j] * expected_MCs[j]) - expected_MCs[j]) * 
    #         market_share_jump(p_int..., p_ext..., p_slope..., j))

    # @NLexpression(model, profit_ext[j=1:N_ent], 
    #     (p_ext[j] - expected_MCs[j+N_inc]) * 
    #         market_share_jump(p_int..., p_ext..., p_slope..., j+N_inc) - cᵃ)


    @expression(model, profit_inc[j=1:N_inc], 
        ((p_int[j] + p_slope[j] * expected_MCs[j]) - expected_MCs[j]) * 
            market_share_jump(p_int..., p_ext..., p_slope..., j))

    @expression(model, profit_ext[j=1:N_ent], 
        (p_ext[j] - expected_MCs[j+N_inc]) * 
            market_share_jump(p_int..., p_ext..., p_slope..., j+N_inc) - cᵃ)

    # Then use sum explicitly in the objective
    @objective(model, Max, 
        sum(profit_inc[j] for j in 1:N_inc) + sum(profit_ext[j] for j in 1:N_ent))

    # Optimizer settings
    set_optimizer_attribute(model, "max_iter", 1000)
    set_optimizer_attribute(model, "tol", 1e-4)
    set_optimizer_attribute(model, "print_level", 0)

    # Solve
    optimize!(model);

    # Get optimal prices
    optimal_p_ext = value.(p_ext)
    optimal_p_int = value.(p_int)
    optimal_slopes = value.(p_slope)
7

Take a careful read of my example and compare it to your original code.

Note that I used @operator instead of register.

You’re hitting the error because the code is now trying to trace your nonlinear expressions, and you calll value inside the market_share_jump function, which is not allowed. (This previously worked because the values of x that get passed in are Float64, and value(x::Float64) is defined to return x.

1 Like
8

thanks!

1 Like