Use type-stable nlpmodels

[odow]

Copied from Optimization example. I'll make a PR once CI is merged

#!/usr/bin/env julia
###### AC-OPF using ADNLPModels ######
#
# implementation reference: https://juliasmoothoptimizers.github.io/ADNLPModels.jl/stable/tutorial/
# other AD libraries can be considered: https://juliasmoothoptimizers.github.io/ADNLPModels.jl/stable/
#

import PowerModels
import Symbolics
import ADNLPModels
import ConcreteStructs
import NLPModelsIpopt

ConcreteStructs.@concrete struct DataRepresentation
    data
    ref
    var_lookup
    var_init
    var_lb
    var_ub
    ref_gen_idxs
    lookup_pg
    lookup_qg
    lookup_va
    lookup_vm
    lookup_lij
    lookup_p_lij
    lookup_q_lij
    cost_arrs
    f_bus
    t_bus
    ref_bus_idxs
    ref_buses_idxs
    ref_bus_gens
    ref_bus_arcs
    ref_branch_idxs
    ref_arcs_from
    ref_arcs_to
    p_idxmap
    q_idxmap
    bus_pd
    bus_qd
    bus_gs
    bus_bs
    br_g
    br_b
    br_tr
    br_ti
    br_ttm
    br_g_fr
    br_b_fr
    br_g_to
    br_b_to
end

function load_and_setup_data(file_name)
    data = PowerModels.parse_file(file_name)
    PowerModels.standardize_cost_terms!(data, order=2)
    PowerModels.calc_thermal_limits!(data)
    ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

    # Some data munging to type-stable forms

    var_lookup = Dict{String,Int}()

    var_init = Float64[]
    var_lb = Float64[]
    var_ub = Float64[]

    var_idx = 1
    for (i,bus) in ref[:bus]
        push!(var_init, 0.0) #va
        push!(var_lb, -Inf)
        push!(var_ub, Inf)
        var_lookup["va_$(i)"] = var_idx
        var_idx += 1

        push!(var_init, 1.0) #vm
        push!(var_lb, bus["vmin"])
        push!(var_ub, bus["vmax"])
        var_lookup["vm_$(i)"] = var_idx
        var_idx += 1
    end

    for (i,gen) in ref[:gen]
        push!(var_init, 0.0) #pg
        push!(var_lb, gen["pmin"])
        push!(var_ub, gen["pmax"])
        var_lookup["pg_$(i)"] = var_idx
        var_idx += 1

        push!(var_init, 0.0) #qg
        push!(var_lb, gen["qmin"])
        push!(var_ub, gen["qmax"])
        var_lookup["qg_$(i)"] = var_idx
        var_idx += 1
    end

    for (l,i,j) in ref[:arcs]
        branch = ref[:branch][l]

        push!(var_init, 0.0) #p
        push!(var_lb, -branch["rate_a"])
        push!(var_ub,  branch["rate_a"])
        var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
        var_idx += 1

        push!(var_init, 0.0) #q
        push!(var_lb, -branch["rate_a"])
        push!(var_ub,  branch["rate_a"])
        var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
        var_idx += 1
    end

    @assert var_idx == length(var_init)+1

    ref_gen_idxs = [i for i in keys(ref[:gen])]
    lookup_pg = Dict{Int,Int}()
    lookup_qg = Dict{Int,Int}()
    lookup_va = Dict{Int,Int}()
    lookup_vm = Dict{Int,Int}()
    lookup_lij = Tuple{Int,Int,Int}[]
    lookup_p_lij = Int[]
    lookup_q_lij = Int[]
    cost_arrs = Dict{Int,Vector{Float64}}()

    for (i,gen) in ref[:gen]
        lookup_pg[i] = var_lookup["pg_$(i)"]
        lookup_qg[i] = var_lookup["qg_$(i)"]
        cost_arrs[i] = gen["cost"]
    end

    for (i,bus) in ref[:bus]
        lookup_va[i] = var_lookup["va_$(i)"]
        lookup_vm[i] = var_lookup["vm_$(i)"]
    end

    for (l,i,j) in ref[:arcs]
        push!(lookup_lij, (l,i,j))
        push!(lookup_p_lij,var_lookup["p_$(l)_$(i)_$(j)"])
        push!(lookup_q_lij,var_lookup["q_$(l)_$(i)_$(j)"])
    end

    f_bus = Dict{Int,Int}()
    t_bus = Dict{Int,Int}()

    for (l,branch) in ref[:branch]
        f_bus[l] = branch["f_bus"]
        t_bus[l] = branch["t_bus"]
    end

    ref_bus_idxs = [i for i in keys(ref[:bus])]
    ref_buses_idxs = [i for i in keys(ref[:ref_buses])]
    ref_bus_gens = ref[:bus_gens]
    ref_bus_arcs = ref[:bus_arcs]
    ref_branch_idxs = [i for i in keys(ref[:branch])]
    ref_arcs_from = ref[:arcs_from]
    ref_arcs_to = ref[:arcs_to]

    p_idxmap = Dict(lookup_lij[i] => lookup_p_lij[i] for i in 1:length(lookup_lij))
    q_idxmap = Dict(lookup_lij[i] => lookup_q_lij[i] for i in 1:length(lookup_lij))

    bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])

    bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])

    for (i,bus) in ref[:bus]
        if length(ref[:bus_loads][i]) > 0
            bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
            bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
        end

        if length(ref[:bus_shunts][i]) > 0
            bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
            bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
        end
    end


    br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_ttm = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])

    for (i,branch) in ref[:branch]
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)

        br_g[i] = g
        br_b[i] = b

        br_tr[i] = tr
        br_ti[i] = ti
        br_ttm[i] = tr^2 + ti^2

        br_g_fr[i] = branch["g_fr"]
        br_b_fr[i] = branch["b_fr"]
        br_g_to[i] = branch["g_to"]
        br_b_to[i] = branch["b_to"]
    end

    return DataRepresentation(
        data,
        ref,
        var_lookup,
        var_init,
        var_lb,
        var_ub,
        ref_gen_idxs,
        lookup_pg,
        lookup_qg,
        lookup_va,
        lookup_vm,
        lookup_lij,
        lookup_p_lij,
        lookup_q_lij,
        cost_arrs,
        f_bus,
        t_bus,
        ref_bus_idxs,
        ref_buses_idxs,
        ref_bus_gens,
        ref_bus_arcs,
        ref_branch_idxs,
        ref_arcs_from,
        ref_arcs_to,
        p_idxmap,
        q_idxmap,
        bus_pd,
        bus_qd,
        bus_gs,
        bus_bs,
        br_g,
        br_b,
        br_tr,
        br_ti,
        br_ttm,
        br_g_fr,
        br_b_fr,
        br_g_to,
        br_b_to)
end

function build_opf_optimization_prob(dataset)
    (;data,
    ref,
    var_lookup,
    var_init,
    var_lb,
    var_ub,
    ref_gen_idxs,
    lookup_pg,
    lookup_qg,
    lookup_va,
    lookup_vm,
    lookup_lij,
    lookup_p_lij,
    lookup_q_lij,
    cost_arrs,
    f_bus,
    t_bus,
    ref_bus_idxs,
    ref_buses_idxs,
    ref_bus_gens,
    ref_bus_arcs,
    ref_branch_idxs,
    ref_arcs_from,
    ref_arcs_to,
    p_idxmap,
    q_idxmap,
    bus_pd,
    bus_qd,
    bus_gs,
    bus_bs,
    br_g,
    br_b,
    br_tr,
    br_ti,
    br_ttm,
    br_g_fr,
    br_b_fr,
    br_g_to,
    br_b_to) = dataset

    function opf_objective(x)
        cost = 0.0
        for i in ref_gen_idxs
            pg = x[lookup_pg[i]]
            _cost_arr = cost_arrs[i]
            cost += _cost_arr[1]*pg^2 + _cost_arr[2]*pg + _cost_arr[3]
        end
        return cost
    end

    function opf_constraints!(ret, x)
        offsetidx = 0

        # va_con
        for (reti,i) in enumerate(ref_buses_idxs)
            ret[reti + offsetidx] = x[lookup_va[i]]
        end

        offsetidx += length(ref_buses_idxs)

        #     @constraint(model,
        #         sum(p[a] for a in ref[:bus_arcs][i]) ==
        #         sum(pg[g] for g in ref_bus_gens[i]) -
        #         sum(load["pd"] for load in bus_loads) -
        #         sum(shunt["gs"] for shunt in bus_shunts)*x[lookup_vm[i]]^2
        #     )
        for (reti,i) in enumerate(ref_bus_idxs)
            ret[reti + offsetidx] = sum(x[lookup_pg[j]] for j in ref_bus_gens[i]; init=0.0) -
                bus_pd[i] -
                bus_gs[i]*x[lookup_vm[i]]^2 -
                sum(x[p_idxmap[a]] for a in ref_bus_arcs[i])
        end

        offsetidx += length(ref_bus_idxs)

        #     @constraint(model,
        #         sum(q[a] for a in ref[:bus_arcs][i]) ==
        #         sum(x[lookup_qg[g]] for g in ref_bus_gens[i]) -
        #         sum(load["qd"] for load in bus_loads) +
        #         sum(shunt["bs"] for shunt in bus_shunts)*x[lookup_vm[i]]^2
        #     )
        for (reti,i) in enumerate(ref_bus_idxs)
            ret[reti + offsetidx] = sum(x[lookup_qg[j]] for j in ref_bus_gens[i]; init=0.0) -
                bus_qd[i] +
                bus_bs[i]*x[lookup_vm[i]]^2 -
                sum(x[q_idxmap[a]] for a in ref_bus_arcs[i])
        end

        offsetidx += length(ref_bus_idxs)

        # @NLconstraint(model, p_fr ==  (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
        # power_flow_p_from_con =
        for (reti,(l,i,j)) in enumerate(ref_arcs_from)
            ret[reti + offsetidx] = (br_g[l]+br_g_fr[l])/br_ttm[l]*x[lookup_vm[f_bus[l]]]^2 +
                (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*cos(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) +
                (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*sin(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) -
                x[p_idxmap[(l,i,j)]]
        end

        offsetidx += length(ref_arcs_from)

        # @NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
        # power_flow_p_to_con
        for (reti,(l,i,j)) in enumerate(ref_arcs_to)
            ret[reti + offsetidx] = (br_g[l]+br_g_to[l])*x[lookup_vm[t_bus[l]]]^2 +
                (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*cos(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) +
                (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*sin(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) -
                x[p_idxmap[(l,i,j)]]
        end

        offsetidx += length(ref_arcs_to)

        # @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
        # power_flow_q_from_con
        for (reti,(l,i,j)) in enumerate(ref_arcs_from)
            ret[reti + offsetidx] = -(br_b[l]+br_b_fr[l])/br_ttm[l]*x[lookup_vm[f_bus[l]]]^2 -
                (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*cos(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) +
                (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*sin(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) -
                x[q_idxmap[(l,i,j)]]
        end

        offsetidx += length(ref_arcs_from)

        # @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
        # power_flow_q_to_con
        for (reti,(l,i,j)) in enumerate(ref_arcs_to)
            ret[reti + offsetidx] = -(br_b[l]+br_b_to[l])*x[lookup_vm[t_bus[l]]]^2 -
                (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*cos(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) +
                (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*sin(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) -
                x[q_idxmap[(l,i,j)]]
        end

        offsetidx += length(ref_arcs_to)

        # @constraint(model, va_fr - va_to <= branch["angmax"])
        # @constraint(model, va_fr - va_to >= branch["angmin"])
        # power_flow_vad_con
        for (reti,(l,i,j)) in enumerate(ref_arcs_from)
            ret[reti + offsetidx] = x[lookup_va[f_bus[l]]] - x[lookup_va[t_bus[l]]]
        end

        offsetidx += length(ref_arcs_from)

        # @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
        # power_flow_mva_from_con
        for (reti,(l,i,j)) in enumerate(ref_arcs_from)
            ret[reti + offsetidx] = x[p_idxmap[(l,i,j)]]^2 + x[q_idxmap[(l,i,j)]]^2
        end

        offsetidx += length(ref_arcs_from)

        # @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
        # power_flow_mva_to_con
        for (reti,(l,i,j)) in enumerate(ref_arcs_to)
            ret[reti + offsetidx] = x[p_idxmap[(l,i,j)]]^2 + x[q_idxmap[(l,i,j)]]^2
        end

        offsetidx += length(ref_arcs_to)

        @assert offsetidx == length(ret)
        return ret
    end

    con_lbs = Float64[]
    con_ubs = Float64[]

    #@constraint(model, va[i] == 0)
    for (i,bus) in ref[:ref_buses]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_balance_p_con
    for (i,bus) in ref[:bus]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_balance_q_con
    for (i,bus) in ref[:bus]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_p_from_con
    for (l,i,j) in ref[:arcs_from]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_p_to_con
    for (l,i,j) in ref[:arcs_to]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_q_from_con
    for (l,i,j) in ref[:arcs_from]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_q_to_con
    for (l,i,j) in ref[:arcs_to]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_vad_con
    for (l,i,j) in ref[:arcs_from]
        branch = ref[:branch][l]
        push!(con_lbs, branch["angmin"])
        push!(con_ubs, branch["angmax"])
    end

    #power_flow_mva_from_con
    for (l,i,j) in ref[:arcs_from]
        branch = ref[:branch][l]
        push!(con_lbs, -Inf)
        push!(con_ubs, branch["rate_a"]^2)
    end

    #power_flow_mva_to_con
    for (l,i,j) in ref[:arcs_to]
        branch = ref[:branch][l]
        push!(con_lbs, -Inf)
        push!(con_ubs, branch["rate_a"]^2)
    end

    nlp = ADNLPModels.ADNLPModel!(opf_objective, var_init, var_lb, var_ub, opf_constraints!, con_lbs, con_ubs, backend = :optimized)
    return (
        nlp = nlp,
        con_lbs = con_lbs,
        con_ubs = con_ubs,
    )
end

function solve_opf(file_name)
    start = time()
    dataset = load_and_setup_data(file_name);
    data_load_time = time() - start
    model = build_opf_optimization_prob(dataset)
    nlp = model.nlp
    model_build_time = time() - start - data_load_time
    output = NLPModelsIpopt.ipopt(nlp)
    cost = output.objective
    feasible = (output.primal_feas <= 1e-6)
    total_time = time() - start
    solve_time = total_time - model_build_time - data_load_time
    model_variables = length(dataset.var_init)
    model_constraints = length(model.con_lbs)
    println("variables: $(model_variables), $(length(dataset.var_lb)), $(length(dataset.var_ub))")
    println("constraints: $(model_constraints), $(length(model.con_lbs)), $(length(model.con_ubs))")
    println("")
    println("\033[1mSummary\033[0m")
    println("   case........: $(file_name)")
    println("   variables...: $(model_variables)")
    println("   constraints.: $(model_constraints)")
    println("   feasible....: $(feasible)")
    println("   cost........: $(round(Int, cost))")
    println("   total time..: $(total_time)")
    println("     data time.: $(data_load_time)")
    println("     build time: $(model_build_time)")
    println("     solve time: $(solve_time)")
    println("")

    return Dict(
        "case" => file_name,
        "variables" => model_variables,
        "constraints" => model_constraints,
        "feasible" => feasible,
        "cost" => cost,
        "time_total" => total_time,
        "time_data" => data_load_time,
        "time_build" => model_build_time,
        "time_solve" => solve_time,
        "solution" => Dict(k => output.solution[v] for (k, v) in dataset.var_lookup),
    )
end

if isinteractive() == false
    solve_opf("$(@__DIR__)/../data/pglib_opf_case5_pjm.m")
end

[ccoffrin]

I am curious about how much this impacts the results on NLPModels. It can live in the variants for testing, but we should defer to the developers of NLPModels, if they want this to be the default implementation or not.

[odow]

on the 118 case, it went from

Dict{String, Any} with 10 entries:
  "cost"        => 97213.6
  "variables"   => 1088
  "constraints" => 1539
  "case"        => "data/pglib_opf_case118_ieee.m"
  "time_total"  => 3.76066
  "time_build"  => 1.11021
  "solution"    => Dict("p_102_66_65"=>-1.66211, "p_132_84_85"=>-0…
  "time_solve"  => 2.60632
  "time_data"   => 0.0441308
  "feasible"    => true

to

Dict{String, Any} with 10 entries:
  "cost"        => 97213.6
  "variables"   => 1088
  "constraints" => 1539
  "case"        => "data/pglib_opf_case118_ieee.m"
  "time_total"  => 2.30815
  "time_build"  => 1.22802
  "solution"    => Dict("p_102_66_65"=>-1.66211, "p_132_84_85"=>-0…
  "time_solve"  => 1.03217
  "time_data"   => 0.0479541
  "feasible"    => true

so the solve is much faster

[odow]

The 793 case went from

Summary
   case........: data/pglib_opf_case793_goc.m
   variables...: 5432
   constraints.: 7978
   feasible....: true
   cost........: 260198
   total time..: 46.25211310386658
     data time.: 0.1509079933166504
     build time: 27.710966110229492
     solve time: 18.390239000320435

to

Summary
   case........: data/pglib_opf_case793_goc.m
   variables...: 5432
   constraints.: 7978
   feasible....: true
   cost........: 260198
   total time..: 39.35335683822632
     data time.: 0.1566319465637207
     build time: 29.251823902130127
     solve time: 9.94490098953247

So again, the solve time is 2X faster. But still only a small reduction in total runtime from 46 seconds to 39 seconds.

[odow]

@tmigot @amontoison thoughts?

[odow]

See #51

[tmigot]

Hi @ccoffrin @odow ! Thank you for maintaining this and the feedback.
I noticed that ADNLPModels had new failures with the last update, but I didn't have the time to investigate.

In general, I think type-stable functions behave better in Julia, so it helps the auto-diff.
Now, the slower build time is maybe because the function tree becomes more complicated !? Although this is just a wild guess.

Overall, I am fine with both solutions, it depends more on the benchmark philosophy.
Either we compare one formulation for all modeling tools, or we also compete on what is the best formulation for each modeling tool. What I don't like about the latter is that it opens the door to people always saying that the benchmark will be biais because you didn't find the best formulation, but it will still be informative.

[ccoffrin]

@tmigot thank you for following up! I have completed a detailed study of the two models for your consideration. The key decisions I would like you to provide are,

1. What version of the NLPModels implementation would you like to be the default?
2. Would you like any other version(s) to sick around as a non-default "variant" of NLPModels?

While you consider these choices, let me note that, my original inception of "rosetta-opf" was to be an educational tool for transfer learning in the Julia optimization ecosystem (like rosetta code is for programming languages). It should help folks understand how to model across the different frameworks. Over time it has become well known as a performance benchmark, however I still say this is not the primary objective of the project. For example, the Symbolic AD variant of the JuMP model is more performant on AC-OPF than the default one. However, I have not made it the default implementation here because it is not the best default AD choice for NLP modeling in JuMP.

With that said, here are the two versions of NLPModels we now have,

NLPModels - this is the current implementation
NLPModels-CS - type-stable implementation using ConcreteStructs

Based on the table below the NLPModels-CS version is consistently faster (around 10%-25% in the large size limit).

Case	Vars	Cons	NLPModels	NLPModels-CS	NLPM/NLPM-CS
case3_lmbd	24	28	5.10e-02	2.89e-02	1.77
case5_pjm	44	53	1.58e-01	1.15e-01	1.37
case14_ieee	118	169	2.73e-01	1.99e-01	1.37
case24_ieee_rts	266	315	7.72e-01	5.81e-01	1.33
case30_ieee	236	348	7.03e-01	5.00e-01	1.40
case30_as	236	348	6.32e-01	5.69e-01	1.11
case39_epri	282	401	6.59e-01	4.05e-01	1.63
case57_ieee	448	675	1.05e+00	1.09e+00	0.97
case60_c	518	737	1.84e+00	9.22e-01	2.00
case73_ieee_rts	824	987	2.56e+00	1.53e+00	1.68
case89_pegase	1042	1649	6.96e+00	3.84e+00	1.81
case118_ieee	1088	1539	5.36e+00	2.45e+00	2.19
case162_ieee_dtc	1484	2313	6.69e+00	4.50e+00	1.49
case179_goc	1468	2200	9.79e+00	4.63e+00	2.11
case197_snem	1608	2397	6.88e+00	4.47e+00	1.54
case200_activ	1456	2116	5.47e+00	3.57e+00	1.53
case240_pserc	2558	3617	5.16e+01	1.82e+01	2.83
case300_ieee	2382	3478	1.22e+01	7.64e+00	1.59
case500_goc	4254	6097	3.01e+01	1.84e+01	1.64
case588_sdet	4110	5979	2.31e+01	1.56e+01	1.48
case793_goc	5432	7978	3.80e+01	2.64e+01	1.44
case1354_pegase	11192	16646	1.61e+02	1.12e+02	1.44
case1803_snem	15246	23172	3.37e+02	2.03e+02	1.66
case1888_rte	14480	21494	5.32e+02	N.D.	N.D.
case1951_rte	15018	22075	3.30e+02	2.26e+02	1.46
case2000_goc	19008	29432	4.05e+02	3.30e+02	1.23
case2312_goc	17128	25716	3.14e+02	2.16e+02	1.45
case2383wp_k	17004	25039	2.76e+02	1.95e+02	1.42
case2736sp_k	19088	28356	3.04e+02	2.48e+02	1.23
case2737sop_k	18988	28358	2.89e+02	2.40e+02	1.20
case2742_goc	24540	38196	1.12e+03	7.05e+02	1.58
case2746wp_k	19520	28446	3.13e+02	2.45e+02	1.28
case2746wop_k	19582	28642	N.D.	2.56e+02	N.D.
case2848_rte	21822	32129	6.41e+02	4.26e+02	1.51
case2853_sdet	23028	33154	8.73e+02	5.33e+02	1.64
case2868_rte	22090	32393	7.12e+02	4.80e+02	1.48
case2869_pegase	25086	37813	8.32e+02	6.72e+02	1.24
case3012wp_k	21082	31029	4.36e+02	3.44e+02	1.27
case3022_goc	23238	34990	7.99e+02	5.53e+02	1.45
case3120sp_k	21608	32092	4.61e+02	3.69e+02	1.25
case3375wp_k	24350	35876	6.49e+02	5.33e+02	1.22
case3970_goc	35270	54428	1.97e+03	1.53e+03	1.28
case4020_goc	36696	56957	2.07e+03	1.62e+03	1.28
case4601_goc	38814	59596	2.30e+03	1.95e+03	1.18
case4619_goc	42532	66289	2.47e+03	2.30e+03	1.07
case4661_sdet	34758	51302	1.51e+03	1.24e+03	1.21
case4837_goc	41398	64030	2.24e+03	2.05e+03	1.09
case4917_goc	37872	56917	2.12e+03	1.68e+03	1.26
case5658_epigrids	48552	74821	3.41e+03	3.04e+03	1.12
case6468_rte	49734	75937	4.18e+03	3.36e+03	1.24
case6470_rte	50482	75976	3.70e+03	3.24e+03	1.14
case6495_rte	50426	76124	4.59e+03	3.70e+03	1.24
case6515_rte	50546	76290	4.20e+03	3.40e+03	1.23
case7336_epigrids	62116	95306	5.35e+03	5.04e+03	1.06
case8387_pegase	78748	118702	N.D.	N.D.	N.D.
case9241_pegase	85568	130826	N.D.	N.D.	N.D.
case9591_goc	83572	130588	N.D.	N.D.	N.D.
case10000_goc	76804	112352	N.D.	N.D.	N.D.
case10192_epigrids	89850	139456	N.D.	N.D.	N.D.
case10480_goc	96750	150874	N.D.	N.D.	N.D.
case13659_pegase	117370	170588	N.D.	N.D.	N.D.
case19402_goc	179562	281733	N.D.	N.D.	N.D.
case20758_epigrids	179236	274918	N.D.	N.D.	N.D.
case24464_goc	203374	313641	N.D.	N.D.	N.D.
case30000_goc	208624	307752	N.D.	N.D.	N.D.
case78484_epigrids	674562	1039062	N.D.	N.D.	N.D.

[tmigot]

Thanks again @ccoffrin for the amount of work this represents.

This is very interesting. Do you have the logs of what failed in case1888_rte? and in case8387_pegase? I am suspecting you get an "out-of-memory" at some point !? So the conclusion seems to be that having the type-stable model as a default seems good.

Back to your primary question:

1. What version of the NLPModels implementation would you like to be the default?
2. Would you like any other version(s) to sick around as a non-default "variant" of NLPModels?
   The simple answer is: I don't really have any other suggestion :).

The reason being that this benchmark uses Ipopt so essentially it requires 3 derivatives from the model:

• the gradient;
• the "sparse" Jacobian and Hessian (of the objective + of the lagrangian).

We are not (at the moment) investigating more ways to compute the sparse Jacobian/Hessian with automatic differentiation (even though one issue JuliaSmoothOptimizers/ADNLPModels.jl#204 is directly linked to this topic), but more interested in ways to compute matrix-vector products directly with AD, i.e. without having to compute/evaluate the whole matrix.
We do however give the possibility for the user to manually provide the sparsity pattern, which can save a ton of time for specific applications (like discretized PDE-constrained problems, where the sparsity pattern is essentially known).