diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 3fd1b4e63..c2369d55d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -23,6 +23,7 @@ jobs: - OptimizationCMAEvolutionStrategy - OptimizationEvolutionary - OptimizationGCMAES + - OptimizationIpopt - OptimizationManopt - OptimizationMetaheuristics - OptimizationMOI @@ -65,7 +66,7 @@ jobs: GROUP: ${{ matrix.group }} - uses: julia-actions/julia-processcoverage@v1 with: - directories: src,lib/OptimizationBBO/src,lib/OptimizationCMAEvolutionStrategy/src,lib/OptimizationEvolutionary/src,lib/OptimizationGCMAES/src,lib/OptimizationManopt/src,lib/OptimizationMOI/src,lib/OptimizationMetaheuristics/src,lib/OptimizationMultistartOptimization/src,lib/OptimizationNLopt/src,lib/OptimizationNOMAD/src,lib/OptimizationOptimJL/src,lib/OptimizationOptimisers/src,lib/OptimizationPolyalgorithms/src,lib/OptimizationQuadDIRECT/src,lib/OptimizationSpeedMapping/src + directories: src,lib/OptimizationBBO/src,lib/OptimizationCMAEvolutionStrategy/src,lib/OptimizationEvolutionary/src,lib/OptimizationGCMAES/src,lib/OptimizationIpopt/src,lib/OptimizationManopt/src,lib/OptimizationMOI/src,lib/OptimizationMetaheuristics/src,lib/OptimizationMultistartOptimization/src,lib/OptimizationNLopt/src,lib/OptimizationNOMAD/src,lib/OptimizationOptimJL/src,lib/OptimizationOptimisers/src,lib/OptimizationPolyalgorithms/src,lib/OptimizationQuadDIRECT/src,lib/OptimizationSpeedMapping/src - uses: codecov/codecov-action@v5 with: file: lcov.info diff --git a/lib/OptimizationIpopt/LICENSE b/lib/OptimizationIpopt/LICENSE new file mode 100644 index 000000000..ac2363b14 --- /dev/null +++ b/lib/OptimizationIpopt/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2025 Sebastian Micluța-Câmpeanu and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/lib/OptimizationIpopt/Project.toml b/lib/OptimizationIpopt/Project.toml new file mode 100644 index 000000000..0a3189f7d --- /dev/null +++ b/lib/OptimizationIpopt/Project.toml @@ -0,0 +1,35 @@ +name = "OptimizationIpopt" +uuid = "43fad042-7963-4b32-ab19-e2a4f9a67124" +authors = ["Sebastian Micluța-Câmpeanu and contributors"] +version = "0.1.0" + +[deps] +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" + +[compat] +Ipopt = "1.10.3" +LinearAlgebra = "1.10.0" +ModelingToolkit = "10.20" +Optimization = "4.3.0" +SciMLBase = "2.90.0" +SparseArrays = "1.10.0" +SymbolicIndexingInterface = "0.3.40" +Zygote = "0.7" +julia = "1.10" + +[extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[targets] +test = ["Aqua", "ModelingToolkit", "Random", "ReverseDiff", "Test", "Symbolics", "Zygote"] diff --git a/lib/OptimizationIpopt/src/OptimizationIpopt.jl b/lib/OptimizationIpopt/src/OptimizationIpopt.jl new file mode 100644 index 000000000..73b97e655 --- /dev/null +++ b/lib/OptimizationIpopt/src/OptimizationIpopt.jl @@ -0,0 +1,219 @@ +module OptimizationIpopt + +using Optimization +using Ipopt +using LinearAlgebra +using SparseArrays +using SciMLBase +using SymbolicIndexingInterface + +export IpoptOptimizer + +const DenseOrSparse{T} = Union{Matrix{T}, SparseMatrixCSC{T}} + +struct IpoptOptimizer end + +function SciMLBase.supports_opt_cache_interface(alg::IpoptOptimizer) + true +end + +function SciMLBase.requiresgradient(opt::IpoptOptimizer) + true +end +function SciMLBase.requireshessian(opt::IpoptOptimizer) + true +end +function SciMLBase.requiresconsjac(opt::IpoptOptimizer) + true +end +function SciMLBase.requiresconshess(opt::IpoptOptimizer) + true +end + +function SciMLBase.allowsbounds(opt::IpoptOptimizer) + true +end +function SciMLBase.allowsconstraints(opt::IpoptOptimizer) + true +end + +include("cache.jl") +include("callback.jl") + +function __map_optimizer_args(cache, + opt::IpoptOptimizer; + maxiters::Union{Number, Nothing} = nothing, + maxtime::Union{Number, Nothing} = nothing, + abstol::Union{Number, Nothing} = nothing, + reltol::Union{Number, Nothing} = nothing, + hessian_approximation = "exact", + verbose = false, + progress = false, + callback = nothing, + kwargs...) + jacobian_sparsity = jacobian_structure(cache) + hessian_sparsity = hessian_lagrangian_structure(cache) + + eval_f(x) = eval_objective(cache, x) + eval_grad_f(x, grad_f) = eval_objective_gradient(cache, grad_f, x) + eval_g(x, g) = eval_constraint(cache, g, x) + function eval_jac_g(x, rows, cols, values) + if values === nothing + for i in 1:length(jacobian_sparsity) + rows[i], cols[i] = jacobian_sparsity[i] + end + else + eval_constraint_jacobian(cache, values, x) + end + return + end + function eval_h(x, rows, cols, obj_factor, lambda, values) + if values === nothing + for i in 1:length(hessian_sparsity) + rows[i], cols[i] = hessian_sparsity[i] + end + else + eval_hessian_lagrangian(cache, values, x, obj_factor, lambda) + end + return + end + + lb = isnothing(cache.lb) ? fill(-Inf, cache.n) : cache.lb + ub = isnothing(cache.ub) ? fill(Inf, cache.n) : cache.ub + + prob = Ipopt.CreateIpoptProblem( + cache.n, + lb, + ub, + cache.num_cons, + cache.lcons, + cache.ucons, + length(jacobian_structure(cache)), + length(hessian_lagrangian_structure(cache)), + eval_f, + eval_g, + eval_grad_f, + eval_jac_g, + eval_h + ) + progress_callback = IpoptProgressLogger(cache.progress, cache, prob) + intermediate = (args...) -> progress_callback(args...) + Ipopt.SetIntermediateCallback(prob, intermediate) + + if !isnothing(maxiters) + Ipopt.AddIpoptIntOption(prob, "max_iter", maxiters) + end + if !isnothing(maxtime) + Ipopt.AddIpoptNumOption(prob, "max_cpu_time", maxtime) + end + if !isnothing(reltol) + Ipopt.AddIpoptNumOption(prob, "tol", reltol) + end + if verbose isa Bool + Ipopt.AddIpoptIntOption(prob, "print_level", verbose * 5) + else + Ipopt.AddIpoptIntOption(prob, "print_level", verbose) + end + Ipopt.AddIpoptStrOption(prob, "hessian_approximation", hessian_approximation) + + for kw in pairs(kwargs) + if kw[2] isa Int + Ipopt.AddIpoptIntOption(prob, string(kw[1]), kw[2]) + elseif kw[2] isa Float64 + Ipopt.AddIpoptNumOption(prob, string(kw[1]), kw[2]) + elseif kw[2] isa String + Ipopt.AddIpoptStrOption(prob, string(kw[1]), kw[2]) + else + error("Keyword argument type $(typeof(kw[2])) not recognized") + end + end + + return prob +end + +function map_retcode(solvestat) + status = Ipopt.ApplicationReturnStatus(solvestat) + if status in [ + Ipopt.Solve_Succeeded, + Ipopt.Solved_To_Acceptable_Level, + Ipopt.User_Requested_Stop, + Ipopt.Feasible_Point_Found + ] + return ReturnCode.Success + elseif status in [ + Ipopt.Infeasible_Problem_Detected, + Ipopt.Search_Direction_Becomes_Too_Small, + Ipopt.Diverging_Iterates + ] + return ReturnCode.Infeasible + elseif status == Ipopt.Maximum_Iterations_Exceeded + return ReturnCode.MaxIters + elseif status in [Ipopt.Maximum_CpuTime_Exceeded + Ipopt.Maximum_WallTime_Exceeded] + return ReturnCode.MaxTime + else + return ReturnCode.Failure + end +end + +function SciMLBase.__solve(cache::IpoptCache) + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime) + + opt_setup = __map_optimizer_args(cache, + cache.opt; + abstol = cache.solver_args.abstol, + reltol = cache.solver_args.reltol, + maxiters = maxiters, + maxtime = maxtime, + cache.solver_args...) + + opt_setup.x .= cache.reinit_cache.u0 + + start_time = time() + status = Ipopt.IpoptSolve(opt_setup) + + opt_ret = map_retcode(status) + + if cache.progress + # Set progressbar to 1 to finish it + Base.@logmsg(Base.LogLevel(-1), "", progress=1, _id=:OptimizationIpopt) + end + + minimum = opt_setup.obj_val + minimizer = opt_setup.x + + stats = Optimization.OptimizationStats(; time = time() - start_time, + iterations = cache.iterations, fevals = cache.f_calls, gevals = cache.f_grad_calls) + + finalize(opt_setup) + + return SciMLBase.build_solution(cache, + cache.opt, + minimizer, + minimum; + original = opt_setup, + retcode = opt_ret, + stats = stats) +end + +function SciMLBase.__init(prob::OptimizationProblem, + opt::IpoptOptimizer; + maxiters::Union{Number, Nothing} = nothing, + maxtime::Union{Number, Nothing} = nothing, + abstol::Union{Number, Nothing} = nothing, + reltol::Union{Number, Nothing} = nothing, + kwargs...) + cache = IpoptCache(prob, opt; + maxiters, + maxtime, + abstol, + reltol, + kwargs... + ) + cache.reinit_cache.u0 .= prob.u0 + + return cache +end + +end # OptimizationIpopt diff --git a/lib/OptimizationIpopt/src/cache.jl b/lib/OptimizationIpopt/src/cache.jl new file mode 100644 index 000000000..1eb5ddac1 --- /dev/null +++ b/lib/OptimizationIpopt/src/cache.jl @@ -0,0 +1,328 @@ +mutable struct IpoptCache{T, F <: OptimizationFunction, RC, LB, UB, I, S, + JT <: DenseOrSparse{T}, HT <: DenseOrSparse{T}, + CHT <: DenseOrSparse{T}, CB, O} <: SciMLBase.AbstractOptimizationCache + const f::F + const n::Int + const num_cons::Int + const reinit_cache::RC + const lb::LB + const ub::UB + const int::I + const lcons::Vector{T} + const ucons::Vector{T} + const sense::S + J::JT + H::HT + cons_H::Vector{CHT} + const callback::CB + const progress::Bool + f_calls::Int + f_grad_calls::Int + iterations::Cint + obj_expr::Union{Expr, Nothing} + cons_expr::Union{Vector{Expr}, Nothing} + const opt::O + const solver_args::NamedTuple +end + +function Base.getproperty(cache::IpoptCache, name::Symbol) + if name in fieldnames(OptimizationBase.ReInitCache) + return getfield(cache.reinit_cache, name) + end + return getfield(cache, name) +end +function Base.setproperty!(cache::IpoptCache, name::Symbol, x) + if name in fieldnames(OptimizationBase.ReInitCache) + return setfield!(cache.reinit_cache, name, x) + end + return setfield!(cache, name, x) +end + +function SciMLBase.get_p(sol::SciMLBase.OptimizationSolution{ + T, + N, + uType, + C +}) where {T, N, uType, C <: IpoptCache} + sol.cache.p +end +function SciMLBase.get_observed(sol::SciMLBase.OptimizationSolution{ + T, + N, + uType, + C +}) where {T, N, uType, C <: IpoptCache} + sol.cache.f.observed +end +function SciMLBase.get_syms(sol::SciMLBase.OptimizationSolution{ + T, + N, + uType, + C +}) where {T, N, uType, C <: IpoptCache} + variable_symbols(sol.cache.f) +end +function SciMLBase.get_paramsyms(sol::SciMLBase.OptimizationSolution{ + T, + N, + uType, + C +}) where {T, N, uType, C <: IpoptCache} + parameter_symbols(sol.cache.f) +end + +function IpoptCache(prob, opt; + callback = nothing, + progress = false, + kwargs...) + reinit_cache = OptimizationBase.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit` + + num_cons = prob.ucons === nothing ? 0 : length(prob.ucons) + if prob.f.adtype isa ADTypes.AutoSymbolics || (prob.f.adtype isa ADTypes.AutoSparse && + prob.f.adtype.dense_ad isa ADTypes.AutoSymbolics) + f = Optimization.instantiate_function( + prob.f, reinit_cache, prob.f.adtype, num_cons; + g = true, h = true, cons_j = true, cons_h = true) + else + f = Optimization.instantiate_function( + prob.f, reinit_cache, prob.f.adtype, num_cons; + g = true, h = true, cons_j = true, cons_vjp = true, lag_h = true) + end + T = eltype(prob.u0) + n = length(prob.u0) + + J = if isnothing(f.cons_jac_prototype) + zeros(T, num_cons, n) + else + convert.(T, f.cons_jac_prototype) + end + lagh = !isnothing(f.lag_hess_prototype) + H = if lagh # lag hessian takes precedence + convert.(T, f.lag_hess_prototype) + elseif !isnothing(f.hess_prototype) + convert.(T, f.hess_prototype) + else + zeros(T, n, n) + end + cons_H = if lagh + Matrix{T}[zeros(T, 0, 0) for i in 1:num_cons] # No need to allocate this up if using lag hessian + elseif isnothing(f.cons_hess_prototype) + Matrix{T}[zeros(T, n, n) for i in 1:num_cons] + else + [convert.(T, f.cons_hess_prototype[i]) for i in 1:num_cons] + end + lcons = prob.lcons === nothing ? fill(T(-Inf), num_cons) : prob.lcons + ucons = prob.ucons === nothing ? fill(T(Inf), num_cons) : prob.ucons + + sys = f.sys isa SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing} ? + nothing : f.sys + obj_expr = f.expr + cons_expr = f.cons_expr + + solver_args = NamedTuple(kwargs) + + return IpoptCache( + f, + n, + num_cons, + reinit_cache, + prob.lb, + prob.ub, + prob.int, + lcons, + ucons, + prob.sense, + J, + H, + cons_H, + callback, + progress, + 0, + 0, + Cint(0), + obj_expr, + cons_expr, + opt, + solver_args + ) +end + +function eval_objective(cache::IpoptCache, x) + l = cache.f(x, cache.p) + cache.f_calls += 1 + return cache.sense === Optimization.MaxSense ? -l : l +end + +function eval_constraint(cache::IpoptCache, g, x) + cache.f.cons(g, x) + return +end + +function eval_objective_gradient(cache::IpoptCache, G, x) + if cache.f.grad === nothing + error("Use OptimizationFunction to pass the objective gradient or " * + "automatically generate it with one of the autodiff backends." * + "If you are using the ModelingToolkit symbolic interface, pass the `grad` kwarg set to `true` in `OptimizationProblem`.") + end + cache.f.grad(G, x) + cache.f_grad_calls += 1 + + if cache.sense === Optimization.MaxSense + G .*= -one(eltype(G)) + end + + return +end + +function jacobian_structure(cache::IpoptCache) + if cache.J isa SparseMatrixCSC + rows, cols, _ = findnz(cache.J) + inds = Tuple{Int, Int}[(i, j) for (i, j) in zip(rows, cols)] + else + rows, cols = size(cache.J) + inds = Tuple{Int, Int}[(i, j) for j in 1:cols for i in 1:rows] + end + return inds +end + +function eval_constraint_jacobian(cache::IpoptCache, j, x) + if isempty(j) + return + elseif cache.f.cons_j === nothing + error("Use OptimizationFunction to pass the constraints' jacobian or " * + "automatically generate i with one of the autodiff backends." * + "If you are using the ModelingToolkit symbolic interface, pass the `cons_j` kwarg set to `true` in `OptimizationProblem`.") + end + # Get and cache the Jacobian object here once. `evaluator.J` calls + # `getproperty`, which is expensive because it calls `fieldnames`. + J = cache.J + cache.f.cons_j(J, x) + if J isa SparseMatrixCSC + nnz = nonzeros(J) + @assert length(j) == length(nnz) + for (i, Ji) in zip(eachindex(j), nnz) + j[i] = Ji + end + else + j .= vec(J) + end + return +end + +function hessian_lagrangian_structure(cache::IpoptCache) + lagh = cache.f.lag_h !== nothing + if cache.f.lag_hess_prototype isa SparseMatrixCSC + rows, cols, _ = findnz(cache.f.lag_hess_prototype) + return Tuple{Int, Int}[(i, j) for (i, j) in zip(rows, cols) if i <= j] + end + sparse_obj = cache.H isa SparseMatrixCSC + sparse_constraints = all(H -> H isa SparseMatrixCSC, cache.cons_H) + if !lagh && !sparse_constraints && any(H -> H isa SparseMatrixCSC, cache.cons_H) + # Some constraint hessians are dense and some are sparse! :( + error("Mix of sparse and dense constraint hessians are not supported") + end + N = length(cache.u0) + inds = if sparse_obj + rows, cols, _ = findnz(cache.H) + Tuple{Int, Int}[(i, j) for (i, j) in zip(rows, cols) if i <= j] + else + Tuple{Int, Int}[(row, col) for col in 1:N for row in 1:col] + end + lagh && return inds + if sparse_constraints + for Hi in cache.cons_H + r, c, _ = findnz(Hi) + for (i, j) in zip(r, c) + if i <= j + push!(inds, (i, j)) + end + end + end + elseif !sparse_obj + # Performance optimization. If both are dense, no need to repeat + else + for col in 1:N, row in 1:col + push!(inds, (row, col)) + end + end + return inds +end + +function eval_hessian_lagrangian(cache::IpoptCache{T}, + h, + x, + σ, + μ) where {T} + if cache.f.lag_h !== nothing + cache.f.lag_h(h, x, σ, Vector(μ)) + + if cache.sense === Optimization.MaxSense + h .*= -one(eltype(h)) + end + + return + end + if cache.f.hess === nothing + error("Use OptimizationFunction to pass the objective hessian or " * + "automatically generate it with one of the autodiff backends." * + "If you are using the ModelingToolkit symbolic interface, pass the `hess` kwarg set to `true` in `OptimizationProblem`.") + end + # Get and cache the Hessian object here once. `evaluator.H` calls + # `getproperty`, which is expensive because it calls `fieldnames`. + H = cache.H + fill!(h, zero(T)) + k = 0 + cache.f.hess(H, x) + sparse_objective = H isa SparseMatrixCSC + if sparse_objective + rows, cols, _ = findnz(H) + for (i, j) in zip(rows, cols) + if i <= j + k += 1 + h[k] = σ * H[i, j] + end + end + else + for i in 1:size(H, 1), j in 1:i + k += 1 + h[k] = σ * H[i, j] + end + end + # A count of the number of non-zeros in the objective Hessian is needed if + # the constraints are dense. + nnz_objective = k + if !isempty(μ) && !all(iszero, μ) + if cache.f.cons_h === nothing + error("Use OptimizationFunction to pass the constraints' hessian or " * + "automatically generate it with one of the autodiff backends." * + "If you are using the ModelingToolkit symbolic interface, pass the `cons_h` kwarg set to `true` in `OptimizationProblem`.") + end + cache.f.cons_h(cache.cons_H, x) + for (μi, Hi) in zip(μ, cache.cons_H) + if Hi isa SparseMatrixCSC + rows, cols, _ = findnz(Hi) + for (i, j) in zip(rows, cols) + if i <= j + k += 1 + h[k] += μi * Hi[i, j] + end + end + else + # The constraints are dense. We only store one copy of the + # Hessian, so reset `k` to where it starts. That will be + # `nnz_objective` if the objective is sprase, and `0` otherwise. + k = sparse_objective ? nnz_objective : 0 + for i in 1:size(Hi, 1), j in 1:i + k += 1 + h[k] += μi * Hi[i, j] + end + end + end + end + + if cache.sense === Optimization.MaxSense + h .*= -one(eltype(h)) + end + + return +end diff --git a/lib/OptimizationIpopt/src/callback.jl b/lib/OptimizationIpopt/src/callback.jl new file mode 100644 index 000000000..ab1f28b0e --- /dev/null +++ b/lib/OptimizationIpopt/src/callback.jl @@ -0,0 +1,82 @@ +struct IpoptState + alg_mod::Cint + iter_count::Cint + obj_value::Float64 + inf_pr::Float64 + inf_du::Float64 + mu::Float64 + d_norm::Float64 + regularization_size::Float64 + alpha_du::Float64 + alpha_pr::Float64 + ls_trials::Cint + z_L::Vector{Float64} + z_U::Vector{Float64} + lambda::Vector{Float64} +end + +struct IpoptProgressLogger{C <: IpoptCache, P} + progress::Bool + cache::C + prob::P +end + +function (cb::IpoptProgressLogger)( + alg_mod::Cint, + iter_count::Cint, + obj_value::Float64, + inf_pr::Float64, + inf_du::Float64, + mu::Float64, + d_norm::Float64, + regularization_size::Float64, + alpha_du::Float64, + alpha_pr::Float64, + ls_trials::Cint +) + n = cb.cache.n + m = cb.cache.num_cons + u, z_L, z_U = zeros(n), zeros(n), zeros(n) + g, lambda = zeros(m), zeros(m) + scaled = false + Ipopt.GetIpoptCurrentIterate(cb.prob, scaled, n, u, z_L, z_U, m, g, lambda) + + original = IpoptState( + alg_mod, + iter_count, + obj_value, + inf_pr, + inf_du, + mu, + d_norm, + regularization_size, + alpha_du, + alpha_pr, + ls_trials, + z_L, + z_U, + lambda + ) + + opt_state = Optimization.OptimizationState(; + iter = Int(iter_count), u, objective = obj_value, original) + cb.cache.iterations = iter_count + + if cb.cache.progress + maxiters = cb.cache.solver_args.maxiters + msg = "objective: " * + sprint(show, obj_value, context = :compact => true) + if !isnothing(maxiters) + # we stop at either convergence or max_steps + Base.@logmsg(Base.LogLevel(-1), msg, progress=iter_count / maxiters, + _id=:OptimizationIpopt) + end + end + if !isnothing(cb.cache.callback) + # return `true` to keep going, or `false` to terminate the optimization + # this is the other way around compared to Optimization.jl callbacks + !cb.cache.callback(opt_state, obj_value) + else + true + end +end diff --git a/lib/OptimizationIpopt/test/additional_tests.jl b/lib/OptimizationIpopt/test/additional_tests.jl new file mode 100644 index 000000000..abc15bcde --- /dev/null +++ b/lib/OptimizationIpopt/test/additional_tests.jl @@ -0,0 +1,262 @@ +using Optimization, OptimizationIpopt +using Zygote +using Test +using LinearAlgebra + +# These tests were automatically translated from the Ipopt tests, https://github.com/coin-or/Ipopt +# licensed under Eclipse Public License - v 2.0 +# https://github.com/coin-or/Ipopt/blob/stable/3.14/LICENSE + +@testset "Additional Ipopt Examples" begin + @testset "Simple 2D Example (MyNLP)" begin + # Based on MyNLP example from Ipopt + # minimize -x[1] - x[2] + # s.t. x[2] - x[1]^2 = 0 + # -1 <= x[1] <= 1 + + function simple_objective(x, p) + return -x[1] - x[2] + end + + function simple_constraint(res, x, p) + res[1] = x[2] - x[1]^2 + end + + optfunc = OptimizationFunction(simple_objective, Optimization.AutoZygote(); + cons = simple_constraint) + prob = OptimizationProblem(optfunc, [0.0, 0.0], nothing; + lb = [-1.0, -Inf], + ub = [1.0, Inf], + lcons = [0.0], + ucons = [0.0]) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] ≈ 1.0 atol=1e-6 + @test sol.u[2] ≈ 1.0 atol=1e-6 + @test sol.objective ≈ -2.0 atol=1e-6 + end + + @testset "Luksan-Vlcek Problem 1" begin + # Based on LuksanVlcek1 from Ipopt examples + # Variable dimension problem + + function lv1_objective(x, p) + n = length(x) + obj = 0.0 + for i in 1:(n-1) + obj += 100 * (x[i]^2 - x[i+1])^2 + (x[i] - 1)^2 + end + return obj + end + + function lv1_constraints(res, x, p) + n = length(x) + for i in 1:(n-2) + res[i] = 3 * x[i+1]^3 + 2 * x[i+2] - 5 + sin(x[i+1] - x[i+2]) * sin(x[i+1] + x[i+2]) + + 4 * x[i+1] - x[i] * exp(x[i] - x[i+1]) - 3 + end + end + + # Test with n = 5 + n = 5 + x0 = fill(3.0, n) + + optfunc = OptimizationFunction(lv1_objective, Optimization.AutoZygote(); + cons = lv1_constraints) + prob = OptimizationProblem(optfunc, x0, nothing; + lcons = fill(4.0, n-2), + ucons = fill(6.0, n-2)) + sol = solve(prob, IpoptOptimizer(); maxiters = 1000, reltol=1e-6) + + @test SciMLBase.successful_retcode(sol) + @test length(sol.u) == n + # Check constraints are satisfied + res = zeros(n-2) + lv1_constraints(res, sol.u, nothing) + @test all(4.0 .<= res .<= 6.0) + end + + @testset "Bound Constrained Quadratic" begin + # Simple bound-constrained quadratic problem + # minimize (x-2)^2 + (y-3)^2 + # s.t. 0 <= x <= 1, 0 <= y <= 2 + + quadratic(x, p) = (x[1] - 2)^2 + (x[2] - 3)^2 + + optfunc = OptimizationFunction(quadratic, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, [0.5, 1.0], nothing; + lb = [0.0, 0.0], + ub = [1.0, 2.0]) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] ≈ 1.0 atol=1e-6 + @test sol.u[2] ≈ 2.0 atol=1e-6 + end + + @testset "Nonlinear Least Squares" begin + # Formulate a nonlinear least squares problem + # minimize sum((y[i] - f(x, t[i]))^2) + # where f(x, t) = x[1] * exp(x[2] * t) + + t = [0.0, 0.5, 1.0, 1.5, 2.0] + y_data = [1.0, 1.5, 2.1, 3.2, 4.8] # Some noisy exponential data + + function nls_objective(x, p) + sum_sq = 0.0 + for i in 1:length(t) + pred = x[1] * exp(x[2] * t[i]) + sum_sq += (y_data[i] - pred)^2 + end + return sum_sq + end + + optfunc = OptimizationFunction(nls_objective, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, [1.0, 0.5], nothing) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.objective < 0.1 # Should fit reasonably well + end + + @testset "Mixed Integer-like Problem (Relaxed)" begin + # Solve a problem that would normally have integer constraints + # but relax to continuous + # minimize x^2 + y^2 + # s.t. x + y >= 3.5 + # 0 <= x, y <= 5 + + objective(x, p) = x[1]^2 + x[2]^2 + + function constraint(res, x, p) + res[1] = x[1] + x[2] + end + + optfunc = OptimizationFunction(objective, Optimization.AutoZygote(); + cons = constraint) + prob = OptimizationProblem(optfunc, [2.0, 2.0], nothing; + lb = [0.0, 0.0], + ub = [5.0, 5.0], + lcons = [3.5], + ucons = [Inf]) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] + sol.u[2] ≈ 3.5 atol=1e-6 + @test sol.u[1] ≈ sol.u[2] atol=1e-6 # By symmetry + end + + @testset "Barrier Method Test" begin + # Test problem where barrier method is particularly relevant + # minimize -log(x) - log(y) + # s.t. x + y <= 2 + # x, y > 0 + + function barrier_objective(x, p) + if x[1] <= 0 || x[2] <= 0 + return Inf + end + return -log(x[1]) - log(x[2]) + end + + function barrier_constraint(res, x, p) + res[1] = x[1] + x[2] + end + + optfunc = OptimizationFunction(barrier_objective, Optimization.AutoZygote(); + cons = barrier_constraint) + prob = OptimizationProblem(optfunc, [0.5, 0.5], nothing; + lb = [1e-6, 1e-6], + ub = [Inf, Inf], + lcons = [-Inf], + ucons = [2.0]) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] + sol.u[2] ≈ 2.0 atol=1e-4 + @test sol.u[1] ≈ 1.0 atol=1e-4 + @test sol.u[2] ≈ 1.0 atol=1e-4 + end + + @testset "Large Scale Sparse Problem" begin + # Create a sparse optimization problem + # minimize sum(x[i]^2) + sum((x[i] - x[i+1])^2) + # s.t. x[1] + x[n] >= 1 + + n = 20 + + function sparse_objective(x, p) + obj = sum(x[i]^2 for i in 1:n) + obj += sum((x[i] - x[i+1])^2 for i in 1:(n-1)) + return obj + end + + function sparse_constraint(res, x, p) + res[1] = x[1] + x[n] + end + + optfunc = OptimizationFunction(sparse_objective, Optimization.AutoZygote(); + cons = sparse_constraint) + x0 = fill(0.1, n) + prob = OptimizationProblem(optfunc, x0, nothing; + lcons = [1.0], + ucons = [Inf]) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] + sol.u[n] >= 1.0 - 1e-6 + end +end + +@testset "Different Hessian Approximations" begin + # Test various Hessian approximation methods + + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + + x0 = [0.0, 0.0] + p = [1.0, 100.0] + + @testset "BFGS approximation" begin + optfunc = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + sol = solve(prob, IpoptOptimizer(); + hessian_approximation = "limited-memory") + + @test SciMLBase.successful_retcode(sol) + @test sol.u ≈ [1.0, 1.0] atol=1e-4 + end + + @testset "SR1 approximation" begin + optfunc = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + sol = solve(prob, IpoptOptimizer(); + hessian_approximation = "limited-memory", + limited_memory_update_type = "sr1") + + @test SciMLBase.successful_retcode(sol) + @test sol.u ≈ [1.0, 1.0] atol=1e-4 + end +end + +@testset "Warm Start Tests" begin + # Test warm starting capabilities + + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + + x0 = [0.5, 0.5] + p = [1.0, 100.0] + + optfunc = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + + # First solve + cache = init(prob, IpoptOptimizer()) + sol1 = solve!(cache) + + @test SciMLBase.successful_retcode(sol1) + @test sol1.u ≈ [1.0, 1.0] atol=1e-4 + + # Note: Full warm start testing would require modifying the problem + # and resolving, which needs reinit!/remake functionality +end diff --git a/lib/OptimizationIpopt/test/advanced_features.jl b/lib/OptimizationIpopt/test/advanced_features.jl new file mode 100644 index 000000000..1e34b253d --- /dev/null +++ b/lib/OptimizationIpopt/test/advanced_features.jl @@ -0,0 +1,246 @@ +using Optimization, OptimizationIpopt +using Zygote +using Test +using LinearAlgebra +using SparseArrays + +# These tests were automatically translated from the Ipopt tests, https://github.com/coin-or/Ipopt +# licensed under Eclipse Public License - v 2.0 +# https://github.com/coin-or/Ipopt/blob/stable/3.14/LICENSE + +@testset "Advanced Ipopt Features" begin + + @testset "Custom Tolerances and Options" begin + # Test setting various Ipopt-specific options + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + + x0 = [0.0, 0.0] + p = [1.0, 100.0] + + optfunc = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + + # Test with tight tolerances + sol = solve(prob, IpoptOptimizer(); + reltol = 1e-10, + acceptable_tol = 1e-8, + acceptable_iter = 5) + + @test SciMLBase.successful_retcode(sol) + @test sol.u ≈ [1.0, 1.0] atol=1e-8 + end + + @testset "Constraint Violation Tolerance" begin + # Test problem with different constraint tolerances + function obj(x, p) + return x[1]^2 + x[2]^2 + end + + function cons(res, x, p) + res[1] = x[1] + x[2] - 2.0 + res[2] = x[1]^2 + x[2]^2 - 2.0 + end + + optfunc = OptimizationFunction(obj, Optimization.AutoZygote(); cons = cons) + prob = OptimizationProblem(optfunc, [0.5, 0.5], nothing; + lcons = [0.0, 0.0], + ucons = [0.0, 0.0]) + + sol = solve(prob, IpoptOptimizer(); + constr_viol_tol = 1e-8) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] + sol.u[2] ≈ 2.0 atol=1e-7 + @test sol.u[1]^2 + sol.u[2]^2 ≈ 2.0 atol=1e-7 + end + + @testset "Derivative Test" begin + # Test with derivative checking enabled + function complex_obj(x, p) + return sin(x[1]) * cos(x[2]) + exp(-x[1]^2 - x[2]^2) + end + + optfunc = OptimizationFunction(complex_obj, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, [0.1, 0.1], nothing) + + # Run with derivative test level 1 (first derivatives only) + sol = solve(prob, IpoptOptimizer(); + derivative_test = "first-order", + derivative_test_tol = 1e-4) + + @test SciMLBase.successful_retcode(sol) + end + + @testset "Linear Solver Options" begin + # Test different linear solver options if available + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + + x0 = zeros(10) # Larger problem + p = [1.0, 100.0] + + # Extend Rosenbrock to n dimensions + function rosenbrock_n(x, p) + n = length(x) + sum = 0.0 + for i in 1:2:n-1 + sum += (p[1] - x[i])^2 + p[2] * (x[i+1] - x[i]^2)^2 + end + return sum + end + + optfunc = OptimizationFunction(rosenbrock_n, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + + # Test with different linear solver strategies + sol = solve(prob, IpoptOptimizer(); + linear_solver = "mumps") # or "ma27", "ma57", etc. if available + + @test SciMLBase.successful_retcode(sol) + # Check that odd indices are close to 1 + @test all(isapprox(sol.u[i], 1.0, atol=1e-4) for i in 1:2:length(x0)-1) + end + + @testset "Scaling Options" begin + # Test problem that benefits from scaling + function scaled_obj(x, p) + return 1e6 * x[1]^2 + 1e-6 * x[2]^2 + end + + function scaled_cons(res, x, p) + res[1] = 1e3 * x[1] + 1e-3 * x[2] - 1.0 + end + + optfunc = OptimizationFunction(scaled_obj, Optimization.AutoZygote(); + cons = scaled_cons) + prob = OptimizationProblem(optfunc, [1.0, 1.0], nothing; + lcons = [0.0], + ucons = [0.0]) + + # Solve with automatic scaling + sol = solve(prob, IpoptOptimizer(); + nlp_scaling_method = "gradient-based") + + @test SciMLBase.successful_retcode(sol) + # Check constraint satisfaction + res = zeros(1) + scaled_cons(res, sol.u, nothing) + @test abs(res[1]) < 1e-6 + end + + @testset "Restoration Phase Test" begin + # Problem that might trigger restoration phase + function difficult_obj(x, p) + return x[1]^4 + x[2]^4 + end + + function difficult_cons(res, x, p) + res[1] = x[1]^3 + x[2]^3 - 1.0 + res[2] = x[1]^2 + x[2]^2 - 0.5 + end + + optfunc = OptimizationFunction(difficult_obj, Optimization.AutoZygote(); + cons = difficult_cons) + # Start from an infeasible point + prob = OptimizationProblem(optfunc, [2.0, 2.0], nothing; + lcons = [0.0, 0.0], + ucons = [0.0, 0.0]) + + sol = solve(prob, IpoptOptimizer(); + required_infeasibility_reduction = 0.9) + + if SciMLBase.successful_retcode(sol) + # Check constraint satisfaction if successful + res = zeros(2) + difficult_cons(res, sol.u, nothing) + @test norm(res) < 1e-4 + end + end + + @testset "Mu Strategy Options" begin + # Test different barrier parameter update strategies + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + + x0 = [0.0, 0.0] + p = [1.0, 100.0] + + optfunc = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + + # Test adaptive mu strategy + sol = solve(prob, IpoptOptimizer(); + mu_strategy = "adaptive", + mu_init = 0.1) + + @test SciMLBase.successful_retcode(sol) + @test sol.u ≈ [1.0, 1.0] atol=1e-4 + + # Test monotone mu strategy + sol2 = solve(prob, IpoptOptimizer(); + mu_strategy = "monotone") + + @test SciMLBase.successful_retcode(sol2) + @test sol2.u ≈ [1.0, 1.0] atol=1e-4 + end + + @testset "Fixed Variable Handling" begin + # Test problem with effectively fixed variables + function fixed_var_obj(x, p) + return (x[1] - 1)^2 + (x[2] - 2)^2 + (x[3] - 3)^2 + end + + optfunc = OptimizationFunction(fixed_var_obj, Optimization.AutoZygote()) + # Fix x[2] = 2.0 by setting equal bounds + prob = OptimizationProblem(optfunc, [0.0, 2.0, 0.0], nothing; + lb = [-Inf, 2.0, -Inf], + ub = [Inf, 2.0, Inf]) + + sol = solve(prob, IpoptOptimizer(); + fixed_variable_treatment = "make_parameter") + + @test SciMLBase.successful_retcode(sol) + @test sol.u ≈ [1.0, 2.0, 3.0] atol=1e-6 + end + + @testset "Acceptable Point Termination" begin + # Test reaching an acceptable point rather than optimal + function slow_converge_obj(x, p) + return sum(exp(-10 * (x[i] - i/10)^2) for i in 1:length(x)) + end + + n = 5 + optfunc = OptimizationFunction(slow_converge_obj, + Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, zeros(n), nothing; + sense = Optimization.MaxSense) + + sol = solve(prob, IpoptOptimizer(); + acceptable_tol = 1e-4, + acceptable_iter = 10, + maxiters = 50) + + @test SciMLBase.successful_retcode(sol) + end +end + +@testset "Output and Logging Options" begin + # Test various output options + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + + x0 = [0.0, 0.0] + p = [1.0, 100.0] + + optfunc = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optfunc, x0, p) + + @testset "Verbose levels" begin + for verbose_level in [false, 0, 3, 5] + sol = solve(prob, IpoptOptimizer(); verbose = verbose_level) + @test SciMLBase.successful_retcode(sol) + end + end + + @testset "Timing statistics" begin + sol = solve(prob, IpoptOptimizer(); print_timing_statistics = "yes") + @test SciMLBase.successful_retcode(sol) + end +end diff --git a/lib/OptimizationIpopt/test/problem_types.jl b/lib/OptimizationIpopt/test/problem_types.jl new file mode 100644 index 000000000..9db84f19e --- /dev/null +++ b/lib/OptimizationIpopt/test/problem_types.jl @@ -0,0 +1,328 @@ +using Optimization, OptimizationIpopt +using Zygote +using Test +using LinearAlgebra +using SparseArrays + +# These tests were automatically translated from the Ipopt tests, https://github.com/coin-or/Ipopt +# licensed under Eclipse Public License - v 2.0 +# https://github.com/coin-or/Ipopt/blob/stable/3.14/LICENSE + +@testset "Specific Problem Types" begin + + @testset "Optimal Control Problem" begin + # Discretized optimal control problem + # minimize integral of u^2 subject to dynamics + + N = 20 # number of time steps + dt = 0.1 + + function control_objective(z, p) + # z = [x1, x2, ..., xN, u1, u2, ..., uN-1] + # Minimize control effort + u_start = N + 1 + return sum(z[i]^2 for i in u_start:length(z)) + end + + function dynamics_constraints(res, z, p) + # Enforce dynamics x[i+1] = x[i] + dt * u[i] + for i in 1:N-1 + res[i] = z[i+1] - z[i] - dt * z[N + i] + end + # Initial condition + res[N] = z[1] - 0.0 + # Final condition + res[N+1] = z[N] - 1.0 + end + + n_vars = N + (N-1) # states + controls + n_cons = N + 1 # dynamics + boundary conditions + + optfunc = OptimizationFunction(control_objective, AutoZygote(); + cons = dynamics_constraints) + z0 = zeros(n_vars) + prob = OptimizationProblem(optfunc, z0; + lcons = zeros(n_cons), + ucons = zeros(n_cons)) + + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] ≈ 0.0 atol=1e-6 # Initial state + @test sol.u[N] ≈ 1.0 atol=1e-6 # Final state + end + + @testset "Portfolio Optimization" begin + # Markowitz portfolio optimization with constraints + # minimize risk (variance) subject to return constraint + + n_assets = 5 + # Expected returns (random for example) + μ = [0.05, 0.10, 0.15, 0.08, 0.12] + # Covariance matrix (positive definite) + Σ = [0.05 0.01 0.02 0.01 0.00; + 0.01 0.10 0.03 0.02 0.01; + 0.02 0.03 0.15 0.02 0.03; + 0.01 0.02 0.02 0.08 0.02; + 0.00 0.01 0.03 0.02 0.06] + + target_return = 0.10 + + function portfolio_risk(w, p) + return dot(w, Σ * w) + end + + function portfolio_constraints(res, w, p) + # Sum of weights = 1 + res[1] = sum(w) - 1.0 + # Expected return >= target + res[2] = dot(μ, w) - target_return + end + + optfunc = OptimizationFunction(portfolio_risk, AutoZygote(); + cons = portfolio_constraints) + w0 = fill(1.0/n_assets, n_assets) + prob = OptimizationProblem(optfunc, w0; + lb = zeros(n_assets), # No short selling + ub = ones(n_assets), # No leverage + lcons = [0.0, 0.0], + ucons = [0.0, Inf]) + + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sum(sol.u) ≈ 1.0 atol=1e-6 + @test dot(μ, sol.u) >= target_return - 1e-6 + @test all(sol.u .>= -1e-6) # Non-negative weights + end + + @testset "Geometric Programming" begin + # Geometric program in standard form + # minimize c^T * x subject to geometric constraints + + function geometric_obj(x, p) + # Objective: x1 * x2 * x3 (in log form: log(x1) + log(x2) + log(x3)) + return exp(x[1]) * exp(x[2]) * exp(x[3]) + end + + function geometric_cons(res, x, p) + # Constraint: x1^2 * x2 / x3 <= 1 + # In exponential form: 2*x1 + x2 - x3 <= 0 + res[1] = exp(2*x[1] + x[2] - x[3]) - 1.0 + # Constraint: x1 + x2 + x3 = 1 (in exponential variables) + res[2] = exp(x[1]) + exp(x[2]) + exp(x[3]) - 3.0 + end + + optfunc = OptimizationFunction(geometric_obj, AutoZygote(); + cons = geometric_cons) + x0 = zeros(3) # log variables start at 0 (original variables = 1) + prob = OptimizationProblem(optfunc, x0; + lcons = [-Inf, 0.0], + ucons = [0.0, 0.0]) + + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + # Check constraints + res = zeros(2) + geometric_cons(res, sol.u, nothing) + @test res[1] <= 1e-6 + @test abs(res[2]) <= 1e-6 + end + + @testset "Parameter Estimation" begin + # Nonlinear least squares parameter estimation + # Fit exponential decay model: y = a * exp(-b * t) + c + + # Generate synthetic data + true_params = [2.0, 0.5, 0.1] + t_data = collect(0:0.5:5) + y_data = @. true_params[1] * exp(-true_params[2] * t_data) + true_params[3] + # Add noise + # y_data += 0.03 * randn(length(t_data)) + y_data += [0.05, 0.01, 0.01, 0.025, 0.0001, 0.004, 0.0056, 0.003, 0.0076, 0.012, 0.0023] + + function residual_sum_squares(params, p) + a, b, c = params + residuals = @. y_data - (a * exp(-b * t_data) + c) + return sum(residuals.^2) + end + + optfunc = OptimizationFunction(residual_sum_squares, AutoZygote()) + # Initial guess + params0 = [1.0, 1.0, 0.0] + prob = OptimizationProblem(optfunc, params0; + lb = [0.0, 0.0, -1.0], # a, b > 0 + ub = [10.0, 10.0, 1.0]) + + sol = solve(prob, IpoptOptimizer(), tol=1e-10, acceptable_tol=1e-10) + + @test SciMLBase.successful_retcode(sol) + # Parameters should be close to true values (within noise) + @test sol.u[1] ≈ true_params[1] atol=0.2 + @test sol.u[2] ≈ true_params[2] atol=0.1 + @test sol.u[3] ≈ true_params[3] atol=0.05 + end + + @testset "Network Flow Problem" begin + # Minimum cost flow problem + # Simple network: source -> 2 intermediate nodes -> sink + + # Network structure (4 nodes, 5 edges) + # Node 1: source, Node 4: sink + # Edges: (1,2), (1,3), (2,3), (2,4), (3,4) + + # Edge costs + costs = [2.0, 3.0, 1.0, 4.0, 2.0] + # Edge capacities + capacities = [10.0, 8.0, 5.0, 10.0, 10.0] + # Required flow from source to sink + required_flow = 15.0 + + function flow_cost(flows, p) + return dot(costs, flows) + end + + function flow_constraints(res, flows, p) + # Conservation at node 2: flow in = flow out + res[1] = flows[1] - flows[3] - flows[4] + # Conservation at node 3: flow in = flow out + res[2] = flows[2] + flows[3] - flows[5] + # Total flow from source + res[3] = flows[1] + flows[2] - required_flow + # Total flow to sink + res[4] = flows[4] + flows[5] - required_flow + end + + optfunc = OptimizationFunction(flow_cost, Optimization.AutoZygote(); + cons = flow_constraints) + flows0 = fill(required_flow / 2, 5) + prob = OptimizationProblem(optfunc, flows0, nothing; + lb = zeros(5), + ub = capacities, + lcons = zeros(4), + ucons = zeros(4)) + + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test all(sol.u .>= -1e-6) # Non-negative flows + @test all(sol.u .<= capacities .+ 1e-6) # Capacity constraints + # Check flow conservation + res = zeros(4) + flow_constraints(res, sol.u, nothing) + @test norm(res) < 1e-6 + end + + @testset "Robust Optimization" begin + # Simple robust optimization problem + # minimize worst-case objective over uncertainty set + + function robust_objective(x, p) + # Minimize max_{u in U} (x - u)^T * (x - u) + # where U = {u : ||u||_inf <= 0.5} + # This simplifies to minimizing ||x||^2 + ||x||_1 + return sum(x.^2) + sum(abs.(x)) + end + + function robust_constraints(res, x, p) + # Constraint: sum(x) >= 1 + res[1] = sum(x) - 1.0 + end + + n = 3 + optfunc = OptimizationFunction(robust_objective, Optimization.AutoZygote(); + cons = robust_constraints) + x0 = fill(1.0/n, n) + prob = OptimizationProblem(optfunc, x0, nothing; + lcons = [0.0], + ucons = [Inf]) + + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) + @test sum(sol.u) >= 1.0 - 1e-6 + end + + # @testset "Complementarity Constraint" begin + # # Mathematical program with complementarity constraints (MPCC) + # # Reformulated using smoothing + + # function mpcc_objective(x, p) + # return (x[1] - 1)^2 + (x[2] - 2)^2 + # end + + # function mpcc_constraints(res, x, p) + # # Original complementarity: x[1] * x[2] = 0 + # # Smoothed version: x[1] * x[2] <= epsilon + # ε = 1e-6 + # res[1] = x[1] * x[2] - ε + # # Additional constraint: x[1] + x[2] >= 1 + # res[2] = x[1] + x[2] - 1.0 + # end + + # optfunc = OptimizationFunction(mpcc_objective, Optimization.AutoZygote(); + # cons = mpcc_constraints) + # x0 = [0.5, 0.5] + # prob = OptimizationProblem(optfunc, x0, nothing; + # lb = [0.0, 0.0], + # lcons = [-Inf, 0.0], + # ucons = [0.0, Inf]) + + # sol = solve(prob, IpoptOptimizer()) + + # @test SciMLBase.successful_retcode(sol) + # # Should satisfy approximate complementarity + # @test sol.u[1] * sol.u[2] < 1e-4 + # @test sol.u[1] + sol.u[2] >= 1.0 - 1e-6 + # end +end + +@testset "Stress Tests" begin + @testset "High-dimensional Problem" begin + # Large-scale quadratic program + n = 100 + + # Random positive definite matrix + A = randn(n, n) + Q = A' * A + I + b = randn(n) + + function large_quadratic(x, p) + return 0.5 * dot(x, Q * x) - dot(b, x) + end + + optfunc = OptimizationFunction(large_quadratic, Optimization.AutoZygote()) + x0 = randn(n) + prob = OptimizationProblem(optfunc, x0) + + sol = solve(prob, IpoptOptimizer(); + maxiters = 1000) + + @test SciMLBase.successful_retcode(sol) + # Check optimality: gradient should be near zero + grad = Q * sol.u - b + @test norm(grad) < 1e-4 + end + + @testset "Highly Nonlinear Problem" begin + # Trigonometric test problem + function trig_objective(x, p) + n = length(x) + return sum(sin(x[i])^2 * cos(x[i])^2 + + exp(-abs(x[i] - π/4)) for i in 1:n) + end + + n = 10 + optfunc = OptimizationFunction(trig_objective, Optimization.AutoZygote()) + x0 = randn(n) + prob = OptimizationProblem(optfunc, x0; + lb = fill(-2π, n), + ub = fill(2π, n)) + + sol = solve(prob, IpoptOptimizer(); + hessian_approximation = "limited-memory") + + @test SciMLBase.successful_retcode(sol) + end +end diff --git a/lib/OptimizationIpopt/test/runtests.jl b/lib/OptimizationIpopt/test/runtests.jl new file mode 100644 index 000000000..40f3faf20 --- /dev/null +++ b/lib/OptimizationIpopt/test/runtests.jl @@ -0,0 +1,117 @@ +using Optimization, OptimizationIpopt +using Zygote +using Symbolics +using Test +using SparseArrays +using ModelingToolkit +using ReverseDiff + +rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 +x0 = zeros(2) +_p = [1.0, 100.0] +l1 = rosenbrock(x0, _p) + +optfunc = OptimizationFunction((x, p) -> -rosenbrock(x, p), Optimization.AutoZygote()) +prob = OptimizationProblem(optfunc, x0, _p; sense = Optimization.MaxSense) + +callback = function (state, l) + display(l) + return false +end + +sol = solve(prob, IpoptOptimizer(); callback, hessian_approximation = "exact") +@test SciMLBase.successful_retcode(sol) +@test sol ≈ [1, 1] + +sol = solve(prob, IpoptOptimizer(); callback, hessian_approximation = "limited-memory") +@test SciMLBase.successful_retcode(sol) +@test sol ≈ [1, 1] + +function _test_sparse_derivatives_hs071(backend, optimizer) + function objective(x, ::Any) + return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3] + end + function constraints(res, x, ::Any) + res .= [ + x[1] * x[2] * x[3] * x[4], + x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 + ] + end + prob = OptimizationProblem( + OptimizationFunction(objective, backend; cons = constraints), + [1.0, 5.0, 5.0, 1.0]; + sense = Optimization.MinSense, + lb = [1.0, 1.0, 1.0, 1.0], + ub = [5.0, 5.0, 5.0, 5.0], + lcons = [25.0, 40.0], + ucons = [Inf, 40.0]) + sol = solve(prob, optimizer) + @test isapprox(sol.objective, 17.014017145179164; atol = 1e-6) + x = [1.0, 4.7429996418092970, 3.8211499817883077, 1.3794082897556983] + @test isapprox(sol.u, x; atol = 1e-6) + @test prod(sol.u) >= 25.0 - 1e-6 + @test isapprox(sum(sol.u .^ 2), 40.0; atol = 1e-6) + return +end + +@testset "backends" begin + backends = ( + AutoForwardDiff(), + AutoReverseDiff(), + AutoSparse(AutoForwardDiff()) + ) + for backend in backends + @testset "$backend" begin + _test_sparse_derivatives_hs071(backend, IpoptOptimizer()) + end + end +end + +# Include additional tests based on Ipopt examples +# These tests were automatically translated from the Ipopt tests, https://github.com/coin-or/Ipopt +# licensed under Eclipse Public License - v 2.0 +# https://github.com/coin-or/Ipopt/blob/stable/3.14/LICENSE +include("additional_tests.jl") +include("advanced_features.jl") +include("problem_types.jl") + +@testset "tutorial" begin + rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 + x0 = zeros(2) + _p = [1.0, 1.0] + + cons(res, x, p) = (res .= [x[1]^2 + x[2]^2, x[1] * x[2]]) + + function lagh(res, x, sigma, mu, p) + lH = sigma * [2 + 8(x[1]^2) * p[2]-4(x[2] - (x[1]^2)) * p[2] -4p[2]*x[1] + -4p[2]*x[1] 2p[2]] .+ [2mu[1] mu[2] + mu[2] 2mu[1]] + res .= lH[[1, 3, 4]] + end + lag_hess_prototype = sparse([1 1; 0 1]) + + optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); + cons = cons, lag_h = lagh, lag_hess_prototype) + prob = OptimizationProblem(optprob, x0, _p, lcons = [1.0, 0.5], ucons = [1.0, 0.5]) + sol = solve(prob, IpoptOptimizer()) + + @test SciMLBase.successful_retcode(sol) +end + +@testset "MTK cache" begin + @variables x + @parameters a = 1.0 + @named sys = OptimizationSystem((x - a)^2, [x], [a];) + sys = complete(sys) + prob = OptimizationProblem(sys, [x => 0.0]; grad = true, hess = true) + cache = init(prob, IpoptOptimizer(); verbose = false) + @test cache isa OptimizationIpopt.IpoptCache + sol = solve!(cache) + @test sol.u ≈ [1.0] # ≈ [1] + + @test_broken begin # needs reinit/remake fixes + cache = Optimization.reinit!(cache; p = [2.0]) + sol = solve!(cache) + @test sol.u ≈ [2.0] # ≈ [2] + end +end