Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion docs/src/numerical.md
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,10 @@ ModelAnalyzer.Numerical.SmallObjectiveCoefficient
ModelAnalyzer.Numerical.LargeObjectiveCoefficient
ModelAnalyzer.Numerical.SmallObjectiveQuadraticCoefficient
ModelAnalyzer.Numerical.LargeObjectiveQuadraticCoefficient
ModelAnalyzer.Numerical.NonconvexQuadraticConstraint
ModelAnalyzer.Numerical.SmallMatrixQuadraticCoefficient
ModelAnalyzer.Numerical.LargeMatrixQuadraticCoefficient
ModelAnalyzer.Numerical.NonconvexQuadraticObjective
ModelAnalyzer.Numerical.NonconvexQuadraticConstraint
```

These issues are saved in the data structure that is returned from the `ModelAnalyzer.analyze` function:
Expand Down
153 changes: 13 additions & 140 deletions src/feasibility.jl
Original file line number Diff line number Diff line change
Expand Up @@ -815,7 +815,8 @@ function _analyze_complementarity!(model, data)
obj = JuMP.constraint_object(con)
func = obj.func
set = obj.set
func_val = JuMP.value(x -> data.primal_point[x], func) - _set_value(set)
func_val =
JuMP.value.(x -> data.primal_point[x], func) - _set_value(set)
comp_val = MOI.Utilities.set_dot(func_val, data.dual_point[con], set)
if abs(comp_val) > data.atol
push!(data.complementarity, ComplemetarityViolation(con, comp_val))
Expand All @@ -824,9 +825,14 @@ function _analyze_complementarity!(model, data)
return
end

function _set_value(set::MOI.AbstractScalarSet)
return 0.0
end
# not needed because it would have stoped in dualization before
# function _set_value(set::MOI.AbstractScalarSet)
# return 0.0
# end
# function _set_value(set::MOI.Interval)
# error("Interval sets are not supported.")
# return (set.lower, set.upper)
# end

function _set_value(set::MOI.AbstractVectorSet)
return zeros(MOI.dimension(set))
Expand All @@ -844,11 +850,6 @@ function _set_value(set::MOI.EqualTo)
return set.value
end

function _set_value(set::MOI.Interval)
error("Interval sets are not supported.")
return (set.lower, set.upper)
end

function _analyze_objectives!(
model::JuMP.GenericModel{T},
dual_model,
Expand Down Expand Up @@ -915,13 +916,9 @@ function _analyze_objectives!(
end

###

# unsafe as is its checked upstream
function _last_primal_solution(model::JuMP.GenericModel)
if !JuMP.has_values(model)
error(
"No primal solution is available. You must provide a point at " *
"which to check feasibility.",
)
end
return Dict(v => JuMP.value(v) for v in JuMP.all_variables(model))
end

Expand Down Expand Up @@ -1138,137 +1135,13 @@ function _fix_ret(
return ret
end

function _add_with_resize!(vec, val, i)
if i > length(vec)
resize!(vec, i)
end
return vec[i] = val
end

"""
dual_feasibility_report(
point::Function,
model::GenericModel{T};
atol::T = zero(T),
skip_missing::Bool = false,
) where {T}

A form of `dual_feasibility_report` where a function is passed as the first
argument instead of a dictionary as the second argument.

## Example

```jldoctest
julia> model = Model();

julia> @variable(model, 0.5 <= x <= 1, start = 1.3); TODO

julia> dual_feasibility_report(model) do v
return dual_start_value(v)
end
Dict{Any, Float64} with 1 entry:
x ≤ 1 => 0.3 TODO
```
"""
# probablye remove this method
function dual_feasibility_report(
point::Function,
model::JuMP.GenericModel{T};
atol::T = zero(T),
skip_missing::Bool = false,
) where {T}
if JuMP.num_nonlinear_constraints(model) > 0
error(
"Nonlinear constraints are not supported. " *
"Use `dual_feasibility_report` instead.",
)
end
if !skip_missing
constraint_list = JuMP.all_constraints(
model;
include_variable_in_set_constraints = false,
)
for c in constraint_list
if !haskey(point, c)
error(
"point does not contain a dual for constraint $c. Provide " *
"a dual, or pass `skip_missing = true`.",
)
end
end
end
dual_model = _dualize2(model)
map = dual_model.ext[:dualization_primal_dual_map].primal_con_dual_var

dual_var_primal_con = _reverse_primal_con_dual_var_map(map)

function dual_point(jump_dual_var::JuMP.GenericVariableRef{T})
# v is a variable in the dual jump model
# we need the associated cosntraint in the primal jump model
moi_dual_var = JuMP.index(jump_dual_var)
moi_primal_con, i = dual_var_primal_con[moi_dual_var]
jump_primal_con = JuMP.constraint_ref_with_index(model, moi_primal_con)
pre_point = point(jump_primal_con)
if ismissing(pre_point)
if !skip_missing
error(
"point does not contain a dual for constraint $jump_primal_con. Provide " *
"a dual, or pass `skip_missing = true`.",
)
else
return missing
end
end
return point(jump_primal_con)[i]
end

dual_con_to_violation = JuMP.primal_feasibility_report(
dual_point,
dual_model;
atol = atol,
skip_missing = skip_missing,
)

# some dual model constraints are associated with primal model variables (primal_con_dual_var)
# if variable is free
primal_var_dual_con =
dual_model.ext[:dualization_primal_dual_map].primal_var_dual_con
# if variable is bounded
primal_convar_dual_con =
dual_model.ext[:dualization_primal_dual_map].constrained_var_dual
# other dual model constraints (bounds) are associated with primal model constraints (non-bounds)
primal_con_dual_con =
dual_model.ext[:dualization_primal_dual_map].primal_con_dual_con

dual_con_primal_all = _build_dual_con_primal_all(
primal_var_dual_con,
primal_convar_dual_con,
primal_con_dual_con,
)

ret = _fix_ret(dual_con_to_violation, model, dual_con_primal_all)

return ret
end

function _reverse_primal_con_dual_var_map(primal_con_dual_var)
dual_var_primal_con =
Dict{MOI.VariableIndex,Tuple{MOI.ConstraintIndex,Int}}()
for (moi_con, vec_vars) in primal_con_dual_var
for (i, moi_var) in enumerate(vec_vars)
dual_var_primal_con[moi_var] = (moi_con, i)
end
end
return dual_var_primal_con
end

function _dualize2(
model::JuMP.GenericModel,
optimizer_constructor = nothing;
kwargs...,
)
mode = JuMP.mode(model)
if mode != JuMP.AUTOMATIC
if mode == JuMP.MANUAL
error("Dualization does not support solvers in $(mode) mode")
end
dual_model = JuMP.Model()
Expand Down
17 changes: 10 additions & 7 deletions src/infeasibility.jl
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,9 @@ function ModelAnalyzer.analyze(
end
iis = iis_elastic_filter(model, optimizer)
# for now, only one iis is computed
push!(out.iis, IrreducibleInfeasibleSubset(iis))
if iis !== nothing
push!(out.iis, IrreducibleInfeasibleSubset(iis))
end

return out
end
Expand All @@ -278,7 +280,7 @@ function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer)
println(
"iis resolver cannot continue because model is found to be $(status) by the solver",
)
return
return nothing
end

model, reference_map = JuMP.copy_model(original_model)
Expand Down Expand Up @@ -315,19 +317,20 @@ function iis_elastic_filter(original_model::JuMP.GenericModel, optimizer)
var = collect(keys(func.terms))
coef1 = func.terms[var[1]]
coef2 = func.terms[var[2]]
if JuMP.value(var1) > tolerance && JuMP.value(var2) > tolerance
if JuMP.value(var[1]) > tolerance &&
JuMP.value(var[2]) > tolerance
error("IIS failed due numerical instability")
elseif JuMP.value(var[1]) > tolerance
has_lower = JuMP.has_lower_bound(var[1])
JuMP.fix(var[1], 0.0; force = true)
# or delete(model, var1)
# or delete(model, var[1])
delete!(constraint_to_affine, con)
constraint_to_affine[con] = coef2 * var[2]
push!(de_elastisized, (con, var[1], has_lower))
elseif JuMP.value(var[2]) > tolerance
has_lower = JuMP.has_lower_bound(var[2])
JuMP.fix(var[2], 0.0; force = true)
# or delete(model, var2)
# or delete(model, var[2])
delete!(constraint_to_affine, con)
constraint_to_affine[con] = coef1 * var[1]
push!(de_elastisized, (con, var[2], has_lower))
Expand Down Expand Up @@ -570,7 +573,7 @@ function ModelAnalyzer._summarize(
end

function ModelAnalyzer._summarize(io::IO, issue::IrreducibleInfeasibleSubset)
return print(io, "IIS: ", join(map(issue.constraint, _name), ", "))
return print(io, "IIS: ", join(map(_name, issue.constraint), ", "))
end

function ModelAnalyzer._verbose_summarize(
Expand Down Expand Up @@ -629,7 +632,7 @@ function ModelAnalyzer._verbose_summarize(
return print(
io,
"Irreducible Infeasible Subset: ",
join(map(issue.constraint, _name), ", "),
join(map(_name, issue.constraint), ", "),
)
end

Expand Down
35 changes: 19 additions & 16 deletions src/intervals.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,28 +66,31 @@ function Base.iszero(a::Interval)
return iszero(a.hi) && iszero(a.lo)
end

Base.:+(a::Interval) = a
Base.:-(a::Interval) = Interval(-a.hi, -a.lo)
# this code is only used for interval += scalar_coef * interval
# so only bivariate + and * are needes

function Base.:+(a::Interval{T}, b::T) where {T<:Real}
return Interval(a.lo + b, a.hi + b)
end
Base.:+(b::T, a::Interval{T}) where {T<:Real} = a + b
# Base.:+(a::Interval) = a
# Base.:-(a::Interval) = Interval(-a.hi, -a.lo)

function Base.:-(a::Interval{T}, b::T) where {T<:Real}
return Interval(a.lo - b, a.hi - b)
end
function Base.:-(b::T, a::Interval{T}) where {T<:Real}
return Interval(b - a.hi, b - a.lo)
end
# function Base.:+(a::Interval{T}, b::T) where {T<:Real}
# return Interval(a.lo + b, a.hi + b)
# end
# Base.:+(b::T, a::Interval{T}) where {T<:Real} = a + b

# function Base.:-(a::Interval{T}, b::T) where {T<:Real}
# return Interval(a.lo - b, a.hi - b)
# end
# function Base.:-(b::T, a::Interval{T}) where {T<:Real}
# return Interval(b - a.hi, b - a.lo)
# end

function Base.:+(a::Interval{T}, b::Interval{T}) where {T<:Real}
return Interval(a.lo + b.lo, a.hi + b.hi)
end

function Base.:-(a::Interval{T}, b::Interval{T}) where {T<:Real}
return Interval(a.lo - b.hi, a.hi - b.lo)
end
# function Base.:-(a::Interval{T}, b::Interval{T}) where {T<:Real}
# return Interval(a.lo - b.hi, a.hi - b.lo)
# end

## Multiplication
function Base.:*(x::T, a::Interval{T}) where {T<:Real}
Expand All @@ -99,6 +102,6 @@ function Base.:*(x::T, a::Interval{T}) where {T<:Real}
end
end

Base.:*(a::Interval{T}, x::T) where {T<:Real} = x * a
# Base.:*(a::Interval{T}, x::T) where {T<:Real} = x * a

Base.convert(::Type{Interval{T}}, x::T) where {T<:Real} = Interval(x, x)
Loading