WIP: add a replacement for Symbolics.jl

[odow]

x-ref https://github.com/jump-dev/MathOptInterface.jl/issues/2533

[ccoffrin]

Can we also update the package deps in this PR?

[odow]

I did in #38. I'll tag a release

[ccoffrin]

Is this still required? https://github.com/lanl-ansi/MathOptSymbolicAD.jl/blob/od/symbolic/Project.toml#L14

[odow]

Hopefully by the time I finish this PR, no. It's still a WIP

[codecov[bot]]

Welcome to Codecov :tada:

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

Thanks for integrating Codecov - We've got you covered :open_umbrella:

[odow]

import MathOptInterface as MOI
import MathOptSymbolicAD

function _expr_to_symbolics(expr::MathOptSymbolicAD._Node)
    if expr.head == MathOptSymbolicAD._OP_COEFFICIENT
        return MOI.VariableIndex(-expr.index)
    elseif expr.head == MathOptSymbolicAD._OP_VARIABLE
        return MOI.VariableIndex(expr.index)
    elseif expr.head == MathOptSymbolicAD._OP_OPERATION
        return MOI.ScalarNonlinearFunction(
            expr.operation,
            Any[_expr_to_symbolics(c) for c in expr.children],
        )
    else
        @assert expr.head == MathOptSymbolicAD._OP_INTEGER_EXPONENT
        return MOI.ScalarNonlinearFunction(
            :^,
            Any[_expr_to_symbolics(expr.children[1]), expr.index],
        )
    end
end

x = MOI.VariableIndex.(1:4)
model = MOI.Nonlinear.Model()
for i in 1:3
    MOI.Nonlinear.add_constraint(model, :($(x[4]) * log($(x[i]) + $i)), MOI.LessThan(2.0))
end

variable_order = Dict(xi.value => i for (i, xi) in enumerate(x))
functions = map(values(model.constraints)) do c
    expr = MathOptSymbolicAD._to_expr(model, c.expression, variable_order, Any[])
    return MathOptSymbolicAD._Function(model, expr)
end

f = first(functions)
g = _expr_to_symbolics(f.expr)
x, ∇f, H, ∇²f = MathOptSymbolicAD.gradient_and_hessian(x -> x.value > 0, g)

julia> x
2-element Vector{MathOptInterface.VariableIndex}:
 MOI.VariableIndex(1)
 MOI.VariableIndex(2)

julia> ∇f
2-element Vector{Any}:
 log(+(MOI.VariableIndex(2), MOI.VariableIndex(-1)))
 *(MOI.VariableIndex(1), /((1), +(MOI.VariableIndex(2), MOI.VariableIndex(-1))))

julia> H
2-element Vector{Tuple{Int64, Int64}}:
 (1, 2)
 (2, 2)

julia> ∇²f
2-element Vector{Any}:
 /((1), +(MOI.VariableIndex(2), MOI.VariableIndex(-1)))
 *(MOI.VariableIndex(1), /((-1), ^(+(MOI.VariableIndex(2), MOI.VariableIndex(-1)), (2))))

[odow]

using JuMP
import PowerModels
import MathOptSymbolicAD
import Ipopt
function main(case, backend)
    pm = PowerModels.instantiate_model(
        case,
        PowerModels.ACPPowerModel,
        PowerModels.build_opf,
    );
    model = pm.model
    set_optimizer(model, Ipopt.Optimizer)
    set_attribute(
        model,
        MOI.AutomaticDifferentiationBackend(),
        backend,
    )
    optimize!(model)
    return solution_summary(model)
end
a = main("test/data/pglib_opf_case13659_pegase.m", MathOptSymbolicAD.SymbolicAD.SymbolicMode())
b = main("test/data/pglib_opf_case13659_pegase.m", MathOptSymbolicAD.DefaultBackend())
c = main("test/data/pglib_opf_case13659_pegase.m", MOI.Nonlinear.SparseReverseMode())

julia> a
* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 8.94805e+06
  Dual objective value : 8.83438e+06

* Work counters
  Solve time (sec)   : 6.70950e+01
  Barrier iterations : 66


julia> b
* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 8.94805e+06
  Dual objective value : 8.83438e+06

* Work counters
  Solve time (sec)   : 8.47096e+01
  Barrier iterations : 66


julia> c
* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 8.94805e+06
  Dual objective value : 8.83438e+06

* Work counters
  Solve time (sec)   : 7.91290e+01
  Barrier iterations : 66

[odow]

Even this now works:

julia> using JuMP, Ipopt

julia> import MathOptSymbolicAD: SymbolicAD

julia> begin
           model = Model(Ipopt.Optimizer)
           @variable(model, x)
           @operator(model, op_foo, 1, x -> x^2)
           @objective(model, Min, op_foo(x - 0.5))
           attr = MOI.AutomaticDifferentiationBackend()
           set_attribute(model, attr, SymbolicAD.SymbolicMode())
           optimize!(model)
           assert_is_solved_and_feasible(model)
           @show value(x)
       end;
This is Ipopt version 3.14.17, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        1

Total number of variables............................:        1
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.5000000e-01 0.00e+00 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 0.00e+00 0.00e+00  -1.7 5.00e-01    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 0.018

EXIT: Optimal Solution Found.
value(x) = 0.5