User defined function inside JuMP Objective resulting in MethodError

Specific Domains Optimization (Mathematical)
, , , ,
1

Hi there,

I’m trying to formulate an optimization model that maximizes profit while dealing with the transportation of raw supplies and production of goods. The profit calculation is lengthy and depends upon both decision variables. So, instead of including it in the Objective statement, I defined a function that takes as inputs the decision variables and performs the calculations necessary for profit. However, I receive the following error:
MethodError: Cannot convert an object of type AffExpr to an object of type Float64

Here is my MWE:

using CSV, DataFrames, JuMP, Gurobi

# Load input data
df = DataFrame(a = rand(5),b = rand(5))

# Number of municipalities
n = nrow(df)

# Manure generated
S = df[!, 1]

# Initialize model
model = Model(Gurobi.Optimizer)

# Load parameters
param = 5

# Decision Variables
@variable(model, x[1:n, 1:2], Bin)  # Biomethane (1) and Electricity (2) production in muni i
@variable(model, y[1:n,1:n] >= 0)  # Manure transported (kg/day) from muni i to j

# Here is a function to calculate the network wide profit from the production of Biomethane and/or Electricity.
# Profit is defined as revenue from Biomethane and Electricity minus infrastructure and transportation costs.
# This function is designed to be called in the objective function.

function network_profit(df::DataFrame, y::Matrix{VariableRef}, x::Matrix{VariableRef})
    # Some logic to calculate profit based on manure amount and decision (either biomethane or electricity)
    # Return calculated profit
    
    n = nrow(df)
    
    # initialize manure column ("ArgumentError: Cannot assign to non-existent column")
    df[!, :manure_supply_kg_day] = df[!,1]  
   
    for i in 1:n
        for j in 1:n
            if i != j # Avoid subtracting or adding when the source and destination are the same
                # Subtract resources sent from municipality 'i' to 'j'
                df[i, :manure_supply_kg_day] -= y[i,j]

                # Add resources received by municipality 'i' from 'j'
                df[i, :manure_supply_kg_day] += y[j,i]
            end
        end
    end
   
    # digester size
    df[!, :digester_size_m3] = df[!, :manure_supply_kg_day] ./ 1000 .* 17
    # 1000...manure density? ; 17...digestion time?

    # digester capital cost
    df[!, :digester_capital_cost_USD] = param .* (df[!, :digester_size_m3] ./ param) .^ param

    # annualized digester capital cost (USD/yr)
    df[!, :annualized_digester_capital_cost_USD_yr] = df[!, :digester_capital_cost_USD] .* (param / 100)

    # digester O&M (USD/yr)
    df[!, :digester_OM_USD_yr] = df[!, :digester_capital_cost_USD] .* param

    # biogas production (m3/day)
    df[!, :biogas_production_m3_day] = df[!, :manure_supply_kg_day] .* df[!, 2]
  
    # Apply muni specific CH4 yield to the transported manure
    # Add methane for incoming manure and subtract methane for outgoing manure,
    # while taking into account the specific CH4 yield of the municipality where the manure originates.
    
    # initialize methane column ("ArgumentError: cannot assign to a non-existent column")
    df[!, :methane_production_m3_day] = zeros(nrow(df))
        
    for i in 1:n
        for j in 1:n
            if i != j # Avoid subtracting or adding when the source and destination are the same
                # Subtract resources sent from municipality 'i' to 'j'
                df[i, :methane_production_m3_day] -= y[i,j] * df[i, 2]

                # Add resources received by municipality 'i' from 'j'
                df[i, :methane_production_m3_day] += y[j,i] * df[i, 2]
            end
        end
    end

    # Biomethane calculations

    # daily biomethane production (m3/day)
    df[!, :daily_biomethane_production_m3_day] = df[!, :methane_production_m3_day] ./ 0.94
    # 0.94

    # hourly biomethane production (m3/hr)
    df[!, :hourly_biomethane_production_m3_hr] = df[!, :daily_biomethane_production_m3_day] ./ 24

    # biomethane capital cost (USD)
    df[!, :biomethane_capital_cost_USD] = df[!, :hourly_biomethane_production_m3_hr] .* param

    # daily biomethane revenue (USD/day)
    df[!, :daily_biomethane_revenue_USD_day] = df[!, :daily_biomethane_production_m3_day] .* param

    # annual biomethane revenue (USD/day)
    df[!, :annual_biomethane_revenue_USD_yr] = df[!, :daily_biomethane_revenue_USD_day] .* 365

    # annualized biomethane capital cost (USD/yr)
    df[!, :annualized_biomethane_capital_cost_USD_yr] = df[!, :biomethane_capital_cost_USD] .* (param / 100)

    # biomethane O&M (USD/yr)
    df[!, :biomethane_OM_USD_yr] = df[!, :daily_biomethane_production_m3_day] .* param .* 365

    # annual biomethane cost
    df[!, :annual_biomethane_cost] = df[!, :annualized_biomethane_capital_cost_USD_yr] .+ df[!, :biomethane_OM_USD_yr]

    # annual biomethane profit
    df[!, :annual_biomethane_profit] = df[!, :annual_biomethane_revenue_USD_yr] .- df[!, :annual_biomethane_cost]
    
    # Electricity calculations

    # daily electricity generated (MWh/day)
    df[!, :daily_electricity_generated_MWh_day] = df[!, :methane_production_m3_day] .* param .* 0.0105

    # daily electricity revenue (USD/day)
    df[!, :daily_electricity_revenue_USD_day] = df[!, :daily_electricity_generated_MWh_day] .* param

    # annual electricity revenue (USD/yr)
    df[!, :annual_electricity_revenue_USD_yr] = df[!, :daily_electricity_revenue_USD_day] .* 365

    # electricity capital cost (USD)
    df[!, :electricity_capital_cost_USD] = param .* (df[!, :daily_electricity_generated_MWh_day] ./ 24)

    # annualized electricity capital cost (USD/yr)
    df[!, :annualized_electricity_capital_cost_USD_yr] = df[!, :electricity_capital_cost_USD] .* (param / 100)

    # electricity O&M (USD/yr)
    df[!, :electricity_OM_USD_yr] = df[!, :electricity_capital_cost_USD] .* (param / 100)

    # annual electricity cost
    df[!, :annual_electricity_cost] = df[!, :annualized_electricity_capital_cost_USD_yr] .+ df[!, :electricity_OM_USD_yr]

    # annual electricity profit
    df[!, :annual_electricity_profit] = df[!, :annual_electricity_revenue_USD_yr] .- df[!, :annual_electricity_cost]
    
    # Aggregate calculations

    # infrastructure cost ($ annualized)
    df[!, :infrastructure_cost_USD_annualized] = df[!, :annualized_digester_capital_cost_USD_yr] .+ df[!, :digester_OM_USD_yr] .+ df[!, :annualized_biomethane_capital_cost_USD_yr] .+ df[!, :biomethane_OM_USD_yr] .+ df[!, :annualized_electricity_capital_cost_USD_yr] .+ df[!, :electricity_OM_USD_yr]

    # annual energy revenue (USD/yr)
    df[!, :annual_energy_revenue_USD_yr] = df[!, :annual_biomethane_revenue_USD_yr] .+ df[!, :annual_electricity_revenue_USD_yr]

    # biogas GHG savings
    df[!, :annual_biogas_GHG_savings] = df[!, :methane_production_m3_day] .* param .* param .* 365

    # manure collection emissions (kg/yr)
    df[!, :manure_collection_emissions_kg_yr] = df[!, 1] .* (param .+ param) .* 365

    # annual GHG offset (kg/yr)
    df[!, :annual_GHG_savings_kg_yr] = df[!, :annual_biogas_GHG_savings] .- df[!, :manure_collection_emissions_kg_yr]

    # annual GHG offset (MM kg/yr)
    df[!, :annual_GHG_savings_MM_kg_yr] = df[!, :annual_GHG_savings_kg_yr] ./ 1000000
    
    # Aggregate calculations for Optimization
   
    # Initialize values
    biometh_rev = 0
    elect_rev = 0
    biometh_inf = 0
    elect_inf = 0

    biometh_rev = sum(df[!, :annual_biomethane_revenue_USD_yr] .* x[:,1])
    elect_rev = sum(df[!, :annual_electricity_revenue_USD_yr] .* x[:,2])
    
    df[!, :elect_infrastructure_cost_USD_annualized] = df[!, :annualized_digester_capital_cost_USD_yr] .+ df[!, :digester_OM_USD_yr] .+ df[!, :annualized_electricity_capital_cost_USD_yr] .+ df[!, :electricity_OM_USD_yr]  
    df[!, :biometh_infrastructure_cost_USD_annualized] = df[!, :annualized_digester_capital_cost_USD_yr] .+ df[!, :digester_OM_USD_yr] .+ df[!, :annualized_biomethane_capital_cost_USD_yr] .+ df[!, :biomethane_OM_USD_yr]
    
    biometh_inf = sum(df[!, :biometh_infrastructure_cost_USD_annualized] .* x[:,1])
    elect_inf = sum(df[!, :elect_infrastructure_cost_USD_annualized] .* x[:,2])

    network_profit = biometh_rev + elect_rev - biometh_inf - elect_inf
    
    return network_profit
end


# Objective
@objective(model, Max, network_profit(df, y, x))

# Constraints

# Material Balance: do not export more manure than available in each muni
for i=1:n
    @constraint(model, sum(y[i,j] for j=1:n) <= S[i])
end

# Processing: each muni can process either Biomethane (1) or Electricity (2)
for j=1:n
    @constraint(model, sum(x[i,1] + x[i,2] for i=1:n) <= 1)
end

optimize!(model)
2

Welcome to the Julia discourse forum.

Can you narrow down from the error where it’s coming from in the code. If I had to guess, it’s probably from here?:

3

Yes, I believe the error is because I’m trying to perform calculations in a df (expecting Float64) using JuMP variables (AffExpr). I can see why this is problematic, but I can’t understand how to perform the calculations for my Objective (profit) without involving the decision variables and my original df.

4

Hi @philds, welcome to the forum :smile:

You’re correct that the error be because the column expects a Float64.

One fix is to use something like df[!, :manure_supply_kg_day] = AffExpr.(df[!,1]). But I think you should consider a much larger change in your approach.

For example:

using DataFrames, JuMP, Gurobi

function network_profit(
    df::DataFrame,
    y::Matrix{VariableRef},
    x::Matrix{VariableRef},
    i,
)
    param = 5
    n = nrow(df)
    manure_supply_kg_day =
        df[i, 1] + sum(y[j, i] - y[i, j] for j in 1:n if i != j)
    digester_size_m3 =
        manure_supply_kg_day / 1_000 * 17
    digester_capital_cost_USD =
        param * (digester_size_m3 / param)^param
    annualized_digester_capital_cost_USD_yr =
        digester_capital_cost_USD * (param / 100)
    digester_OM_USD_yr =
        digester_capital_cost_USD * param
    biogas_production_m3_day =
        manure_supply_kg_day * df[i, 2]
    methane_production_m3_day =
        sum(y[j,i] * df[i, 2] - y[i,j] * df[i, 2] for j in 1:n if i != j)
    daily_biomethane_production_m3_day =
        methane_production_m3_day / 0.94
    hourly_biomethane_production_m3_hr =
        daily_biomethane_production_m3_day / 24
    biomethane_capital_cost_USD =
        hourly_biomethane_production_m3_hr * param
    daily_biomethane_revenue_USD_day =
        daily_biomethane_production_m3_day * param
    annual_biomethane_revenue_USD_yr =
        daily_biomethane_revenue_USD_day * 365
    annualized_biomethane_capital_cost_USD_yr =
        biomethane_capital_cost_USD * (param / 100)
    biomethane_OM_USD_yr =
        daily_biomethane_production_m3_day * param * 365
    annual_biomethane_cost =
        annualized_biomethane_capital_cost_USD_yr + biomethane_OM_USD_yr
    annual_biomethane_profit =
        annual_biomethane_revenue_USD_yr - annual_biomethane_cost
    daily_electricity_generated_MWh_day =
        methane_production_m3_day * param * 0.0105
    daily_electricity_revenue_USD_day =
        daily_electricity_generated_MWh_day * param
    annual_electricity_revenue_USD_yr =
        daily_electricity_revenue_USD_day * 365
    electricity_capital_cost_USD =
        param * (daily_electricity_generated_MWh_day / 24)
    annualized_electricity_capital_cost_USD_yr =
        electricity_capital_cost_USD * (param / 100)
    electricity_OM_USD_yr =
        electricity_capital_cost_USD * (param / 100)
    annual_electricity_cost =
        annualized_electricity_capital_cost_USD_yr + electricity_OM_USD_yr
    annual_electricity_profit =
        annual_electricity_revenue_USD_yr - annual_electricity_cost
    infrastructure_cost_USD_annualized =
        annualized_digester_capital_cost_USD_yr +
        digester_OM_USD_yr +
        annualized_biomethane_capital_cost_USD_yr +
        biomethane_OM_USD_yr +
        annualized_electricity_capital_cost_USD_yr +
        electricity_OM_USD_yr
    annual_energy_revenue_USD_yr =
        annual_biomethane_revenue_USD_yr + annual_electricity_revenue_USD_yr
    annual_biogas_GHG_savings = methane_production_m3_day * param * param * 365
    manure_collection_emissions_kg_yr = df[i, 1] * (param + param) * 365
    annual_GHG_savings_kg_yr =
        annual_biogas_GHG_savings - manure_collection_emissions_kg_yr
    annual_GHG_savings_MM_kg_yr = annual_GHG_savings_kg_yr / 1_000_000
    biometh_rev = annual_biomethane_revenue_USD_yr * x[i, 1]
    elect_rev = annual_electricity_revenue_USD_yr * x[i, 2]
    elect_infrastructure_cost_USD_annualized =
        annualized_digester_capital_cost_USD_yr +
        digester_OM_USD_yr +
        annualized_electricity_capital_cost_USD_yr +
        electricity_OM_USD_yr
    biometh_infrastructure_cost_USD_annualized =
        annualized_digester_capital_cost_USD_yr +
        digester_OM_USD_yr +
        annualized_biomethane_capital_cost_USD_yr +
        biomethane_OM_USD_yr
    biometh_inf = biometh_infrastructure_cost_USD_annualized * x[i, 1]
    elect_inf = elect_infrastructure_cost_USD_annualized * x[i, 2]
    network_profit = biometh_rev + elect_rev - biometh_inf - elect_inf
    return network_profit
end

df = DataFrame(a = rand(5), b = rand(5))
n = nrow(df)
model = Model(Gurobi.Optimizer)
@variable(model, x[1:n, 1:2], Bin)
@variable(model, y[1:n, 1:n] >= 0)
@objective(model, Max, sum(network_profit(df, y, x, i) for i in 1:n))
@constraint(model, [i in 1:n], sum(y[i, :]) <= df[i, 1])
@constraint(model, [i in 1:n], sum(x[i, :]) <= 1)
optimize!(model)

But you have a problem with:

    # digester capital cost
    df[!, :digester_capital_cost_USD] = param .* (df[!, :digester_size_m3] ./ param) .^ param

here digester_size_m3 is an affine function of y, and Gurobi supports only linear or quadratic problems, so it cannot solve problems with y^5.

Was this formulation intended?

2 Likes
5

Thanks for the suggestion @odow! That is a much-improved approach and it worked on my full problem. Also, thanks for pointing out the non-linearity. I did not realize that would cause issues with Gurobi. I will look for a different solver. Thank you!

6

If you have nonlinearity and binary variables, this is a very hard problem to solve.

Before looking for other solvers, I’d ask: is the nonlinearity really necessary? Could you use a quadratic instead? Gurobi can find optimal solutions to mixed-integer quadratic programs.

2 Likes
7

Good point. I introduced the power to model economies of scale from a reference cost ( cost2 = cost1 * (size2/size1)^ 0.85 ). But perhaps I could find a different way to achieve this, such as with a piecewise function.

1 Like