SIP Extension

News.md

@@ -30,3 +30,35 @@
     - Parametric interval methods and the Implicit optimizer have been move to a separate package (to be tagged shortly.)
     - JIT compilation time has been reduced substantially.
     - Support for silent tag and time limits.
+
+## v0.4.0
+  - 6/7/2020: [**EAGO v0.4.0 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.4.0).
+      - Support for new MOI/JuMP `RawParameter` input and a number of new attributes.
+      - Separates McCormick and ReverseMcCormick libraries (now [McCormick.jl](https://github.com/PSORLab/McCormick.jl) and [ReverseMcCormick.jl](https://github.com/PSORLab/ReverseMcCormick.jl))
+        from main package.  McCormick.jl is reexported.
+      - Relaxation calculations now return NaN values on a domain violation.
+      - Tolerance based validation of cuts has been added to generate numerically safe cuts.
+      - Significantly simplify internal codebase for `EAGO.Optimizer` (no changes to API): fully decouples input problem specifications from the formulation used internally, stack only stores variables that are branched on, and a number of internal rearrangements to clearly delineate different routines.
+      - Add problem classification preprocessing that throws to simpler routines if LP problem types are detected (enables future support for SOCP, MILP, MISOCP, and Convex forms).
+      - Fix multiple bugs and add more transparent error codes.
+    - 06/17/2020: [**EAGO v0.4.1 has been tagged**](https://github.com/PSORLab/EAGO.jl/commit/9c1bcf024a19840a0ac49c8c6da13619a5f3845f#comments) Contains minor bug releases.
+    - 08/29/2020: [**EAGO v0.4.2 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.4.2) Support for Julia v1.5.
+
+## v0.5.0
+- 11/18/2020: [**EAGO v0.5.0 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.5.0)
+    - Introduces the `register_eago_operators!(m::JuMP.Model)` which can be used
+        to register all nonstandard nonlinear terms used in EAGO in any JuMP model.
+      - Introduces `positive`, `negative`, `lower_bnd`, `upper_bnd`, and `bnd`
+        functions which can be used to enforce bounds on intermediate terms in
+        nonlinear expressions (EAGO.Optimizer only).
+      - Adds envelopes: `abs2`, `sinpi`, `cospi`, `fma`, `cbrt`.
+      - Adds envelopes and functions: `xlogx`
+      - Adds envelopes of special functions: `erf`, `erfc`, `erfinv`, `erfcinv`
+      - Adds envelopes of activation functions: `relu`, `gelu`, `elu`, `selu`, `swish1`,`sigmoid`, `softsign`, `softplus`,`bisigmoid`, `pentanh`, `leaky_relu`, `param_relu`.
+      - Error messages in `sip_explicit` have been made more transparent.
+      - Fixes some issues with documentation image rendering and links.
+      - Drops appveyor CI and Travis CI in favor of Github Actions.
+- 11/18/2020 [**EAGO v0.5.1 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.5.1)
+     - Support for Julia ~1 (with limited functionality for Julia 1.0, 1.1).
+- 11/18/2020 **EAGO v0.5.2 has been tagged**
+     - Fix user specified branching variables.

Project.toml

@@ -2,7 +2,7 @@ name = "EAGO"
 uuid = "bb8be931-2a91-5aca-9f87-79e1cb69959a"
 
 authors = ["Matthew Wilhelm <[email protected]>"]
-version = "0.5.2"
+version = "0.6.0"
 
 [deps]
 Cassette = "7057c7e9-c182-5462-911a-8362d720325c"

README.md

@@ -99,36 +99,14 @@ As a global optimization platform, EAGO's solvers can be used to find solutions
 The EAGO package has numerous features: a solver accessible from JuMP/MathOptInterface, domain reduction routines, McCormick relaxations, and specialized non-convex semi-infinite program solvers. A full description of all EAGO features is available in the [**documentation website**](https://psorlab.github.io/EAGO.jl/dev/). A series of example have been provided in the form of Jupyter notebooks in the separate [**EAGO-notebooks**](https://github.com/PSORLab/EAGO-notebooks) repository
 
 ## Recent News
-
-- 6/7/2020: [**EAGO v0.4.0 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.4.0).
-  - Support for new MOI/JuMP `RawParameter` input and a number of new attributes.
-  - Separates McCormick and ReverseMcCormick libraries (now [McCormick.jl](https://github.com/PSORLab/McCormick.jl) and [ReverseMcCormick.jl](https://github.com/PSORLab/ReverseMcCormick.jl))
-    from main package.  McCormick.jl is reexported.
-  - Relaxation calculations now return NaN values on a domain violation.
-  - Tolerance based validation of cuts has been added to generate numerically safe cuts.
-  - Significantly simplify internal codebase for `EAGO.Optimizer` (no changes to API): fully decouples input problem specifications from the formulation used internally, stack only stores variables that are branched on, and a number of internal rearrangements to clearly delineate different routines.
-  - Add problem classification preprocessing that throws to simpler routines if LP problem types are detected (enables future support for SOCP, MILP, MISOCP, and Convex forms).
-  - Fix multiple bugs and add more transparent error codes.
-
-- 06/17/2020: [**EAGO v0.4.1 has been tagged**](https://github.com/PSORLab/EAGO.jl/commit/9c1bcf024a19840a0ac49c8c6da13619a5f3845f#comments) Contains minor bug releases.
-- 08/29/2020: [**EAGO v0.4.2 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.4.2) Support for Julia v1.5.
-- 11/18/2020: [**EAGO v0.5.0 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.5.0)
-  - Introduces the `register_eago_operators!(m::JuMP.Model)` which can be used
-    to register all nonstandard nonlinear terms used in EAGO in any JuMP model.
-  - Introduces `positive`, `negative`, `lower_bnd`, `upper_bnd`, and `bnd`
-    functions which can be used to enforce bounds on intermediate terms in
-    nonlinear expressions (EAGO.Optimizer only).
-  - Adds envelopes: `abs2`, `sinpi`, `cospi`, `fma`, `cbrt`.
-  - Adds envelopes and functions: `xlogx`
-  - Adds envelopes of special functions: `erf`, `erfc`, `erfinv`, `erfcinv`
-  - Adds envelopes of activation functions: `relu`, `gelu`, `elu`, `selu`, `swish1`,
-                                             `sigmoid`, `softsign`, `softplus`,
-                                             `bisigmoid`, `pentanh`, `leaky_relu`,
-                                             `param_relu`.
-  - Error messages in `sip_explicit` have been made more transparent.
-  - Fixes some issues with documentation image rendering and links.
-  - Drops appveyor CI and Travis CI in favor of Github Actions.
-- 11/18/2020 [**EAGO v0.5.1 has been tagged**](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.5.1) Support for Julia ~1 (with limited functionality for Julia 1.0, 1.1).
+- 2/5/2021: [EAGO v0.6.0 has been tagged](https://github.com/PSORLab/EAGO.jl/releases/tag/v0.6.0).
+  - License changed from CC BY-NC-SA 4.0 to MIT
+  - Fix deprecated Ipopt constructor
+  - Fix discrepancy between the returned objective value and the objective evaluated at the solution.
+  - Dramatically decrease allocates and first-run performance of SIP routines.
+  - Add two algorithms for the modifications variants to SIPres detailed in Djelassi, H. and Mitsos A. 2017.
+  - Fix objective interval fallback function.
+  - New SIP interface with extendable subroutines.
 
 For a full list of EAGO release news, see click [**here**](https://github.com/PSORLab/EAGO.jl/releases)
 

src/EAGO.jl

@@ -131,6 +131,10 @@ module EAGO
     include("eago_script/script.jl")
 
     # routines for solving SIPs
-    export SIP_Options, SIP_Result, explicit_sip_solve
-    include("eago_semiinfinite/semi_infinite.jl")
+    export SIP_Result, sip_solve, SIPRes, SIPResRev, SIPHybrid
+    export build_model, set_tolerance_inner!, set_tolerance!, get_disc_set,
+           sip_llp!, sip_bnd!, sip_res!, get_sip_optimizer, check_convergence,
+           LowerLevel1, LowerLevel2, LowerLevel3, LowerProblem, UpperProblem,
+           ResProblem
+    include("eago_semiinfinite/semiinfinite.jl")
 end

src/eago_optimizer/functions/nonlinear/forward_pass.jl

@@ -483,12 +483,12 @@ function forward_power!(k::Int64, children_arr::Vector{Int64}, children_idx::Uni
         num2 = 0.0
     end
 
-    # is output a number (by closure of the reals)?
+    # is output a number (by closure of the reals)
     output_is_number = arg1_is_number && arg2_is_number
     numvalued[k] = output_is_number
 
     # x^1 = x
-    if num2 === 1.0
+    if arg2_is_number && (num2 == 1.0)
         if arg1_is_number
              numberstorage[k] = num1
         else
@@ -497,7 +497,7 @@ function forward_power!(k::Int64, children_arr::Vector{Int64}, children_idx::Uni
         return nothing
 
     # x^0 = 1
-    elseif num2 === 0.0
+    elseif arg2_is_number && (num2 == 0.0)
         if arg1_is_number
              numberstorage[k] = 1.0
         else

src/eago_semiinfinite/interface/classify_sip.jl

@@ -0,0 +1,82 @@
+const CR = JuMP.ConstraintRef
+
+function classify_sip!(mSIP::SIPModel, c::AbstractConstraint)
+    error("Function classify_sip! not defined for a constraint of type = $(typeof(c))")
+end
+
+function update_status(status::SIPCons, term_status::SIPCons)
+    if status == SIPNOTSET
+        return term_status
+    elseif (status == DECISION) && (term_status == DECISION)
+        return DECISION
+    elseif (status == UNCERTAIN) && (term_status == UNCERTAIN)
+        return UNCERTAIN
+    end
+    return SEMIINFINITE
+end
+
+classify_sip(mSIP::SIPModel, v::VariableRef) = mSIP.p[v] ? UNCERTAIN : DECISION
+function classify_sip(mSIP::SIPModel, expr::GenericAffExpr{T,VariableRef}) where T
+    status = SIPNOTSET
+    for term in linear_terms(expr)
+        status = update_status(status, classify_sip(mSIP, term[1]))
+        if status == SEMIINFINITE
+            break
+        end
+    end
+    return status
+end
+function classify_sip(mSIP::SIPModel, expr::GenericQuadExpr{T,VariableRef}) where T
+    status = SIPNOTSET
+    for term in quad_terms(expr)
+        status = update_status(status, classify_sip(mSIP, term[2]))
+        status = update_status(status, classify_sip(mSIP, term[3]))
+        if status == SEMIINFINITE
+            break
+        end
+    end
+    for term in linear_terms(expr)
+        status = update_status(status, classify_sip(mSIP, term[1]))
+        if status == SEMIINFINITE
+            break
+        end
+    end
+    return status
+end
+
+for typ in (:(Vector{GenericQuadExpr{T,VariableRef}}),
+            :(Vector{GenericAffExpr{T,VariableRef}}),
+            :(Vector{VariableRef}))
+    @eval function classify_sip(mSIP::SIPModel, expr::Vector{$typ{T,VariableRef}}) where T
+        status = SIPNOTSET
+        for ex in expr
+            status = update_status(status, classify_sip(mSIP, ex))
+            if status == SEMIINFINITE
+                break
+            end
+        end
+        return status
+    end
+end
+
+function classify_sip!(mSIP::SIPModel, cr::ScalarConstraint{F,S}) where {F<:JuMP.AbstractJuMPScalar, S<:JuMP.AbstractScalarSet}
+    mSIP.constraint_type[cr] = classify_sip(mSIP, constraint_object(cr).func)
+    nothing
+end
+
+function classify_sip!(mSIP::SIPModel, cr::CR{M,C,S}) where {M<:JuMP.AbstractModel, C ,S<:JuMP.AbstractShape}
+    mSIP.constraint_type[cr] = classify_sip(constraint_object(cr))
+    nothing
+end
+
+function classify_sip!(mSIP::SIPModel)
+    constraint_types = list_of_constraint_types(mSIP.m)
+    for (func_typ, set_typ) in constraint_types
+        constraint_refs = all_constraints(mSIP.m, func_typ, set_typ)
+        for cr in constraint_refs
+            classify_sip!(mSIP, cr)
+        end
+    end
+    # TODO: ADD CLASSIFICATION OF NONLINEAR EQUATIONS
+    return nothing
+end

src/eago_semiinfinite/interface/scratch_sheet.jl

@@ -0,0 +1,104 @@
+
+#=
+=#
+
+
+
+#=
+Basic idea behind syntax
+mSIP = SIPModel()
+
+@decision_variable(mSIP, 0.0 <= x <= 1.0)
+@decision_variable(mSIP, z, Bin)
+@uncertain_variable(mSIP, p, 0.0 <= x <= 1.0)
+
+
+@constraint(mSIP, ...)
+@objective(mSIP, Min, ...)
+
+@NLconstraint(mSIP, ...)
+@NLobjective(mSIP, Min, ...)
+
+optimize!()
+=#
+
+using JuMP
+
+import JuMP.object_dictionary
+
+@enum(SIPCons, DECISION, UNCERTAIN, SEMIINFINITE, SIPNOTSET)
+
+"""
+"""
+Base.@kwdef mutable struct SIPModel <: JuMP.AbstractModel
+    m::JuMP.Model                                   = JuMP.Model()
+    p::Dict{JuMP.VariableRef,Bool}                  = Dict{JuMP.VariableRef,Bool}()
+    constr_type::Dict{JuMP.ConstraintRef, SIPCons}  = Dict{JuMP.ConstraintRef, SIPCons}()
+    nl_constraint_type::Dict{Int, SIPCons}          = Dict{Int, SIPCons}()
+    nl_expression_type::Dict{Int, SIPCons}          = Dict{Int, SIPCons}()
+end
+
+#=
+One option is to fully extend JuMP... That requires alot of function definitions
+and the only new functionality we actually want to add is classifying constraints
+based on is it a decision variable or is it a
+=#
+macro decision_variable(args...)
+    esc(quote
+        mSIP = $args[1]
+        inputs_2f = $(args[2:end]...)
+        @show mSIP
+        @show inputs_2f
+        #inputs = $(esc(args))
+        #@show inputs
+        #@show typeof(inputs)
+        vi = @variable($(args...))
+        #model = $inputs[1]
+        #for vi in variable_indices
+        #    model.p[vi] = false
+        #end
+    #    return vi
+        end)
+end
+
+mSIP = SIPModel()
+@decision_variable(mSIP, x)
+
+#=
+using MathOptInterface
+
+m = Model()
+@variable(m, 0 <= y <= 1)
+@variable(m, x)
+@variable(m, a, Bin)
+@variable(m, q[i=1:2])
+@constraint(m, y^2 + y + x<= 0)
+@constraint(m, [x, y-1, y-2] in SecondOrderCone())
+@NLconstraint(m, sin(x) + cos(y) <= 0.0)
+@constraint(m, 2x - 1 ⟂ x)
+@constraint(m, q in SOS2([3,5]))
+@constraint(m, a => {x + y <= 1})
+@SDconstraint(m, [x 2x; 3x 4x] >= ones(2, 2))
+
+A = [1 2; 3 4]
+b = [5,6]
+@constraint(m, con, A * x .== b)
+
+
+list = list_of_constraint_types(m)
+
+cons0 = all_constraints(m, VariableRef, MOI.LessThan{Float64})
+cons1 = all_constraints(m, VariableRef, MOI.GreaterThan{Float64})
+cons2 = all_constraints(m, GenericQuadExpr{Float64,VariableRef}, MOI.LessThan{Float64})
+cons2 = all_constraints(m, GenericQuadExpr{Float64,VariableRef}, MOI.LessThan{Float64})
+
+cons0_1 = cons0[1]
+cons1_1 = cons1[1]
+cons2_1 = cons2[1]
+cons3_1 = cons2[1]
+
+out = constraint_object(cons2_1).func
+l_terms = linear_terms(out)
+q_terms = quad_terms(out)
+#all_consts = all_constraints(m, list[1][1], list[1][1])
+=#

src/eago_semiinfinite/interface/sip_model.jl

@@ -0,0 +1,73 @@
+#=
+
+Macro interface is inspired by InfiniteOpt.jl. Many of the code structures are
+reused. Transitioning this over to a fully InfiniteOpt.jl in the future is desirable
+but most of our use cases involve nonlinear functions so that InfiniteOpt.jl
+nonlinear support is gating.
+
+The general idea here is pretty simple. When we add a variable to the model
+the user specifies whether or not it is a decision variable or an uncertain
+variable we then classify constraints. See https://github.com/PSORLab/EAGO.jl/issues/61
+for a discussion of the roadmap here.
+=#
+
+@enum(SIPCons, DECISION, UNCERTAIN, SEMIINFINITE, SIPNOTSET)
+
+
+"""
+"""
+Base.@kwdef mutable struct SIPData
+    p::Dict{JuMP.VariableRef,Bool}                  = Dict{JuMP.VariableRef,Bool}()
+    constr_type::Dict{JuMP.ConstraintRef, SIPCons}  = Dict{JuMP.ConstraintRef, SIPCons}()
+    nl_constraint_type::Dict{Int, SIPCons}          = Dict{Int, SIPCons}()
+    nl_expression_type::Dict{Int, SIPCons}          = Dict{Int, SIPCons}()
+end
+
+function initialize_sip_data(m::JuMP.Model)
+    m.ext[:sip] = SIPData()
+end
+
+function sip_optimizehook(m::JuMP.Model; kwargs...)
+    data = _getsipdata(m)::SIPData
+    if uncertain_variable_num(data) == 0.0
+        ret = JuMP.optimize!(m::JuMP.Model, ignore_optimize_hook = true, kwargs...)
+    else
+        # ret  = model_sip_solve(m)
+    end
+    return ret
+end
+
+function enable_semiinfinite(m::JuMP.Model)
+    haskey(m.ext, :sip) && error("Model already has semiinfinite programs enabled")
+    initialize_sip_data(m)
+    JuMP.set_optimize_hook(m, sip_optimizehook)
+    return nothing
+end
+
+function ModelWithSIP(args...; kwargs...)
+    m = JuMP.Model(args...; kwargs...)
+    enable_semiinfinite(m)
+    return m
+end
+
+"""
+    @uncertain_variable
+
+Add an *anonymous* variable to `m::SIPModel` described by the keyword
+arguments `kw_args` and returns the variable.
+
+    @uncertain_variable(m, expr, args..., kw_args...)
+
+This macro simply denotes the variable as belonging to the set of descision
+variables and then performs `@variable(m, expr, args..., kw_args...)` to
+add the variable to the basic JuMP model for storage. All JuMP syntax is
+supported.
+"""
+macro uncertain_variable(args...)
+    vi = @variable(args...)
+    model = esc(args[1])
+    for vi in variable_indices
+        model.p[vi] = true
+    end
+    return vi
+end