using JuMP, Ipopt

function minimize(
    f::Function,
    A::Matrix,
    b::Vector,
    optimizer::JuMP.OptimizerFactory
)
    M, N = size(A)
    @assert length(b) == M
    model = Model(optimizer)
    @variable(model, x[1:N] >= 0)
    @constraint(model, A * x .== b)
    JuMP.register(model, :f, N, (x...) -> f(collect(x)), autodiff = true)
    @NLobjective(model, Min, f(x...))
    optimize!(model)
    return objective_value(model), value.(x)
end


#  min x[1]^2
# s.t. x[1] - x[2] == 2
#      x >= 0
z, x = minimize(x -> x[1]^2, [1 -1], [2], with_optimizer(Ipopt.Optimizer))