Skip to content

Commit e75fff3

Browse files
odowodow
authored
[Utilities] add PenaltyRelaxation (#1995)
1 parent 4eb55f2 commit e75fff3

File tree

5 files changed

+726
-1
lines changed

5 files changed

+726
-1
lines changed

‎docs/src/submodules/Utilities/overview.md

+61
Original file line numberDiff line numberDiff line change
@@ -315,6 +315,67 @@ $$ \begin{aligned}
315315
In IJulia, calling `print` or ending a cell with
316316
[`Utilities.latex_formulation`](@ref) will render the model in LaTeX.
317317

318+
## Utilities.PenaltyRelaxation
319+
320+
Pass [`Utilities.PenaltyRelaxation`](@ref) to [`modify`](@ref) to relax the
321+
problem by adding penalized slack variables to the constraints. This is helpful
322+
when debugging sources of infeasible models.
323+
324+
```jldoctest
325+
julia> model = MOI.Utilities.Model{Float64}();
326+
327+
julia> x = MOI.add_variable(model);
328+
329+
julia> MOI.set(model, MOI.VariableName(), x, "x")
330+
331+
julia> c = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(2.0));
332+
333+
julia> map = MOI.modify(model, MOI.Utilities.PenaltyRelaxation(Dict(c => 2.0)));
334+
335+
julia> print(model)
336+
Minimize ScalarAffineFunction{Float64}:
337+
0.0 + 2.0 v[2]
338+
339+
Subject to:
340+
341+
ScalarAffineFunction{Float64}-in-LessThan{Float64}
342+
0.0 + 1.0 x - 1.0 v[2] <= 2.0
343+
344+
VariableIndex-in-GreaterThan{Float64}
345+
v[2] >= 0.0
346+
347+
julia> map[c]
348+
MathOptInterface.ScalarAffineFunction{Float64}(MathOptInterface.ScalarAffineTerm{Float64}[MathOptInterface.ScalarAffineTerm{Float64}(1.0, MathOptInterface.VariableIndex(2))], 0.0)
349+
```
350+
351+
You can also modify a single constraint using [`Utilities.ScalarPenaltyRelaxation`](@ref):
352+
```jldoctest
353+
julia> model = MOI.Utilities.Model{Float64}();
354+
355+
julia> x = MOI.add_variable(model);
356+
357+
julia> MOI.set(model, MOI.VariableName(), x, "x")
358+
359+
julia> c = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(2.0));
360+
361+
julia> f = MOI.modify(model, c, MOI.Utilities.ScalarPenaltyRelaxation(2.0));
362+
363+
julia> print(model)
364+
Minimize ScalarAffineFunction{Float64}:
365+
0.0 + 2.0 v[2]
366+
367+
Subject to:
368+
369+
ScalarAffineFunction{Float64}-in-LessThan{Float64}
370+
0.0 + 1.0 x - 1.0 v[2] <= 2.0
371+
372+
VariableIndex-in-GreaterThan{Float64}
373+
v[2] >= 0.0
374+
375+
julia> f
376+
MathOptInterface.ScalarAffineFunction{Float64}(MathOptInterface.ScalarAffineTerm{Float64}[MathOptInterface.ScalarAffineTerm{Float64}(1.0, MathOptInterface.VariableIndex(2))], 0.0)
377+
```
378+
318379
## Utilities.MatrixOfConstraints
319380

320381
The constraints of [`Utilities.Model`](@ref) are stored as a vector of tuples

‎docs/src/submodules/Utilities/reference.md

+7
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,13 @@ Utilities.identity_index_map
9090
Utilities.ModelFilter
9191
```
9292

93+
## Penalty relaxation
94+
95+
```@docs
96+
Utilities.PenaltyRelaxation
97+
Utilities.ScalarPenaltyRelaxation
98+
```
99+
93100
## MatrixOfConstraints
94101

95102
```@docs

‎src/Utilities/Utilities.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ include("mockoptimizer.jl")
7777
include("cachingoptimizer.jl")
7878
include("universalfallback.jl")
7979
include("print.jl")
80-
80+
include("penalty_relaxation.jl")
8181
include("lazy_iterators.jl")
8282

8383
include("precompile.jl")

‎src/Utilities/penalty_relaxation.jl

+296
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,296 @@
1+
# Copyright (c) 2017: Miles Lubin and contributors
2+
# Copyright (c) 2017: Google Inc.
3+
#
4+
# Use of this source code is governed by an MIT-style license that can be found
5+
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
6+
7+
"""
8+
ScalarPenaltyRelaxation(penalty::T) where {T}
9+
10+
A problem modifier that, when passed to [`MOI.modify`](@ref), destructively
11+
modifies the constraint in-place to create a penalized relaxation of the
12+
constraint.
13+
14+
!!! warning
15+
This is a destructive routine that modifies the constraint in-place. If you
16+
don't want to modify the original model, use `JuMP.copy_model` to create a
17+
copy before calling [`MOI.modify`](@ref).
18+
19+
## Reformulation
20+
21+
The penalty relaxation modifies constraints of the form ``f(x) \\in S`` into
22+
``f(x) + y - z \\in S``, where ``y, z \\ge 0``, and then it introduces a penalty
23+
term into the objective of ``a \\times (y + z)`` (if minimizing, else ``-a``),
24+
where ``a`` is `penalty`
25+
26+
When `S` is [`MOI.LessThan`](@ref) or [`MOI.GreaterThan`](@ref), we omit `y` or
27+
`z` respectively as a performance optimization.
28+
29+
## Return value
30+
31+
`MOI.modify(model, ci, ScalarPenaltyRelaxation(penalty))` returns `y + z` as a
32+
[`MOI.ScalarAffineFunction`](@ref). In an optimal solution, query the value of
33+
this function to compute the violation of the constraint.
34+
35+
## Examples
36+
37+
```jldoctest; setup=:(import MathOptInterface; const MOI = MathOptInterface)
38+
julia> model = MOI.Utilities.Model{Float64}();
39+
40+
julia> x = MOI.add_variable(model);
41+
42+
julia> c = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(2.0));
43+
44+
julia> f = MOI.modify(model, c, MOI.Utilities.ScalarPenaltyRelaxation(2.0));
45+
46+
julia> print(model)
47+
Minimize ScalarAffineFunction{Float64}:
48+
0.0 + 2.0 v[2]
49+
50+
Subject to:
51+
52+
ScalarAffineFunction{Float64}-in-LessThan{Float64}
53+
0.0 + 1.0 v[1] - 1.0 v[2] <= 2.0
54+
55+
VariableIndex-in-GreaterThan{Float64}
56+
v[2] >= 0.0
57+
58+
julia> f isa MOI.ScalarAffineFunction{Float64}
59+
true
60+
```
61+
"""
62+
struct ScalarPenaltyRelaxation{T} # <: MOI.AbstractFunctionModification
63+
# We don't make this a subtype of AbstractFunctionModification to avoid some
64+
# ambiguities with generic methods in Utilities and Bridges. The underlying
65+
# reason is that these reformulations can be written using the high-level
66+
# MOI API, so we don't need special handling for bridges and utilities.
67+
penalty::T
68+
end
69+
70+
function _change_sense_to_min_if_necessary(
71+
::Type{T},
72+
model::MOI.ModelLike,
73+
) where {T}
74+
sense = MOI.get(model, MOI.ObjectiveSense())
75+
if sense != MOI.FEASIBILITY_SENSE
76+
return sense
77+
end
78+
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
79+
f = zero(MOI.ScalarAffineFunction{T})
80+
MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
81+
return MOI.MIN_SENSE
82+
end
83+
84+
function MOI.modify(
85+
model::MOI.ModelLike,
86+
ci::MOI.ConstraintIndex{F,<:MOI.AbstractScalarSet},
87+
relax::ScalarPenaltyRelaxation{T},
88+
) where {T,F<:Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}}
89+
sense = _change_sense_to_min_if_necessary(T, model)
90+
y = MOI.add_variable(model)
91+
z = MOI.add_variable(model)
92+
MOI.add_constraint(model, y, MOI.GreaterThan(zero(T)))
93+
MOI.add_constraint(model, z, MOI.GreaterThan(zero(T)))
94+
MOI.modify(model, ci, MOI.ScalarCoefficientChange(y, one(T)))
95+
MOI.modify(model, ci, MOI.ScalarCoefficientChange(z, -one(T)))
96+
scale = sense == MOI.MIN_SENSE ? one(T) : -one(T)
97+
a = scale * relax.penalty
98+
O = MOI.get(model, MOI.ObjectiveFunctionType())
99+
obj = MOI.ObjectiveFunction{O}()
100+
MOI.modify(model, obj, MOI.ScalarCoefficientChange(y, a))
101+
MOI.modify(model, obj, MOI.ScalarCoefficientChange(z, a))
102+
return one(T) * y + one(T) * z
103+
end
104+
105+
function MOI.modify(
106+
model::MOI.ModelLike,
107+
ci::MOI.ConstraintIndex{F,MOI.GreaterThan{T}},
108+
relax::ScalarPenaltyRelaxation{T},
109+
) where {T,F<:Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}}
110+
sense = _change_sense_to_min_if_necessary(T, model)
111+
# Performance optimization: we don't need the z relaxation variable.
112+
y = MOI.add_variable(model)
113+
MOI.add_constraint(model, y, MOI.GreaterThan(zero(T)))
114+
MOI.modify(model, ci, MOI.ScalarCoefficientChange(y, one(T)))
115+
scale = sense == MOI.MIN_SENSE ? one(T) : -one(T)
116+
a = scale * relax.penalty
117+
O = MOI.get(model, MOI.ObjectiveFunctionType())
118+
obj = MOI.ObjectiveFunction{O}()
119+
MOI.modify(model, obj, MOI.ScalarCoefficientChange(y, a))
120+
return one(T) * y
121+
end
122+
123+
function MOI.modify(
124+
model::MOI.ModelLike,
125+
ci::MOI.ConstraintIndex{F,MOI.LessThan{T}},
126+
relax::ScalarPenaltyRelaxation{T},
127+
) where {T,F<:Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}}
128+
sense = _change_sense_to_min_if_necessary(T, model)
129+
# Performance optimization: we don't need the y relaxation variable.
130+
z = MOI.add_variable(model)
131+
MOI.add_constraint(model, z, MOI.GreaterThan(zero(T)))
132+
MOI.modify(model, ci, MOI.ScalarCoefficientChange(z, -one(T)))
133+
scale = sense == MOI.MIN_SENSE ? one(T) : -one(T)
134+
a = scale * relax.penalty
135+
O = MOI.get(model, MOI.ObjectiveFunctionType())
136+
obj = MOI.ObjectiveFunction{O}()
137+
MOI.modify(model, obj, MOI.ScalarCoefficientChange(z, a))
138+
return one(T) * z
139+
end
140+
141+
"""
142+
PenaltyRelaxation(
143+
penalties = Dict{MOI.ConstraintIndex,Float64}();
144+
default::Union{Nothing,T} = 1.0,
145+
)
146+
147+
A problem modifier that, when passed to [`MOI.modify`](@ref), destructively
148+
modifies the model in-place to create a penalized relaxation of the constraints.
149+
150+
!!! warning
151+
This is a destructive routine that modifies the model in-place. If you don't
152+
want to modify the original model, use `JuMP.copy_model` to create a copy
153+
before calling [`MOI.modify`](@ref).
154+
155+
## Reformulation
156+
157+
See [`Utilities.ScalarPenaltyRelaxation`](@ref) for details of the
158+
reformulation.
159+
160+
For each constraint `ci`, the penalty passed to [`Utilities.ScalarPenaltyRelaxation`](@ref)
161+
is `get(penalties, ci, default)`. If the value is `nothing`, because `ci` does
162+
not exist in `penalties` and `default = nothing`, then the constraint is
163+
skipped.
164+
165+
## Return value
166+
167+
`MOI.modify(model, PenaltyRelaxation())` returns a
168+
`Dict{MOI.ConstraintIndex,MOI.ScalarAffineFunction}` that maps each constraint
169+
index to the corresponding `y + z` as a [`MOI.ScalarAffineFunction`](@ref). In
170+
an optimal solution, query the value of these functions to compute the violation
171+
of each constraint.
172+
173+
## Relax a subset of constraints
174+
175+
To relax a subset of constraints, pass a `penalties` dictionary and set
176+
`default = nothing`.
177+
178+
## Supported constraint types
179+
180+
The penalty relaxation is currently limited to modifying
181+
[`MOI.ScalarAffineFunction`](@ref) and [`MOI.ScalarQuadraticFunction`](@ref)
182+
constraints in the linear sets [`MOI.LessThan`](@ref), [`MOI.GreaterThan`](@ref),
183+
[`MOI.EqualTo`](@ref) and [`MOI.Interval`](@ref).
184+
185+
It does not include variable bound or integrality constraints, because these
186+
cannot be modified in-place.
187+
188+
To modify variable bounds, rewrite them as linear constraints.
189+
190+
## Examples
191+
192+
```jldoctest; setup=:(import MathOptInterface; const MOI = MathOptInterface)
193+
julia> model = MOI.Utilities.Model{Float64}();
194+
195+
julia> x = MOI.add_variable(model);
196+
197+
julia> c = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(2.0));
198+
199+
julia> map = MOI.modify(model, MOI.Utilities.PenaltyRelaxation(default = 2.0));
200+
201+
julia> print(model)
202+
Minimize ScalarAffineFunction{Float64}:
203+
0.0 + 2.0 v[2]
204+
205+
Subject to:
206+
207+
ScalarAffineFunction{Float64}-in-LessThan{Float64}
208+
0.0 + 1.0 v[1] - 1.0 v[2] <= 2.0
209+
210+
VariableIndex-in-GreaterThan{Float64}
211+
v[2] >= 0.0
212+
213+
julia> map[c] isa MOI.ScalarAffineFunction{Float64}
214+
true
215+
```
216+
217+
```jldoctest; setup=:(import MathOptInterface; const MOI = MathOptInterface)
218+
julia> model = MOI.Utilities.Model{Float64}();
219+
220+
julia> x = MOI.add_variable(model);
221+
222+
julia> c = MOI.add_constraint(model, 1.0 * x, MOI.LessThan(2.0));
223+
224+
julia> map = MOI.modify(model, MOI.Utilities.PenaltyRelaxation(Dict(c => 3.0)));
225+
226+
julia> print(model)
227+
Minimize ScalarAffineFunction{Float64}:
228+
0.0 + 3.0 v[2]
229+
230+
Subject to:
231+
232+
ScalarAffineFunction{Float64}-in-LessThan{Float64}
233+
0.0 + 1.0 v[1] - 1.0 v[2] <= 2.0
234+
235+
VariableIndex-in-GreaterThan{Float64}
236+
v[2] >= 0.0
237+
238+
julia> map[c] isa MOI.ScalarAffineFunction{Float64}
239+
true
240+
```
241+
"""
242+
mutable struct PenaltyRelaxation{T}
243+
default::Union{Nothing,T}
244+
penalties::Dict{MOI.ConstraintIndex,T}
245+
246+
function PenaltyRelaxation(
247+
p::Dict{MOI.ConstraintIndex,T};
248+
default::Union{Nothing,T} = one(T),
249+
) where {T}
250+
return new{T}(default, p)
251+
end
252+
end
253+
254+
function PenaltyRelaxation(; kwargs...)
255+
return PenaltyRelaxation(Dict{MOI.ConstraintIndex,Float64}(); kwargs...)
256+
end
257+
258+
function PenaltyRelaxation(
259+
d::Dict{<:MOI.ConstraintIndex,T};
260+
kwargs...,
261+
) where {T}
262+
return PenaltyRelaxation(convert(Dict{MOI.ConstraintIndex,T}, d); kwargs...)
263+
end
264+
265+
function MOI.modify(model::MOI.ModelLike, relax::PenaltyRelaxation{T}) where {T}
266+
map = Dict{MOI.ConstraintIndex,MOI.ScalarAffineFunction{T}}()
267+
for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
268+
_modify_penalty_relaxation(map, model, relax, F, S)
269+
end
270+
return map
271+
end
272+
273+
function _modify_penalty_relaxation(
274+
map::Dict{MOI.ConstraintIndex,MOI.ScalarAffineFunction{T}},
275+
model::MOI.ModelLike,
276+
relax::PenaltyRelaxation,
277+
::Type{F},
278+
::Type{S},
279+
) where {T,F,S}
280+
for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
281+
penalty = get(relax.penalties, ci, relax.default)
282+
if penalty === nothing
283+
continue
284+
end
285+
try
286+
map[ci] = MOI.modify(model, ci, ScalarPenaltyRelaxation(penalty))
287+
catch err
288+
if err isa MethodError && err.f == MOI.modify
289+
@warn("Skipping PenaltyRelaxation for ConstraintIndex{$F,$S}")
290+
return
291+
end
292+
rethrow(err)
293+
end
294+
end
295+
return
296+
end

0 commit comments

Comments
 (0)