[Help Wanted] NL Constraints

[charles-r-smith]

Hello, I'm just checking to see whether using @NLconstraints/@NLoptimize is possible in infinite opt? If so, is there an example in any of the documentation? I didn't find any.

The following thread: #16 suggests that it was in the cards at some point, but unsure if it was ever fully developed.

The code I'm trying to write includes the following constraints:

@constraint(m, p == min(p_bar, max(zeros(N), transpose(A)*p + c - x))) 

which causes an error presumably because of the use of min/max. The error is: MethodError: no method matching isless(::GenericQuadExpr{Float64,GeneralVariableRef}, ::Float64)
Closest candidates are:
isless(!Matched::Float64, ::Float64) at float.jl:465
isless(!Matched::Missing, ::Any) at missing.jl:87
isless(!Matched::AbstractFloat, ::AbstractFloat) at operators.jl:156
...

The second constraint is:

prob(x,c,w) = 1 ./(1 .+exp.(-1000*(x.*c-w))) # approximates indicator function of x*c>=w
JuMP.register(m, :prob, 3, prob)
for i = 1:N
    @NLconstraint(m, delta[i] = prob(x[i],c[i],w[i]))
end

This gives the error: MethodError: no method matching register(::InfiniteModel, ::Symbol, ::Int64, ::typeof(prob)). So I assume there's no ability for users to define functions in constraints.

The other way I try to define this constraints:

for i = 1:N
    @NLconstraint(m, delta[i] == 1 ./(1 .+exp.(-1000*(x[i].*c[i]-w[i])))) #ensuring prob x*c >= w = delta
end

This gives the error: MethodError: no method matching _init_NLP(::InfiniteModel)

Thank you,

[pulsipher]

Hi there!

InfiniteOpt.jl does not currently support general nonlinear expressions (e.g., NLconstraints) like JuMP. The details of what expressions are accepted are given on https://pulsipher.github.io/InfiniteOpt.jl/dev/guide/expression/. This limitation is due to JuMP's nonlinear interface being nonextendable for packages like this one (see jump-dev/JuMP.jl#2402). JuMP NLP won't be extendable until it is fundamentally redone via jump-dev/MathOptInterface.jl#846 (maybe in 2-3 years).

However, in accordance with #16, I do plan on implementing our own NLP interface from the ground up (aiming at sometime this year) since JuMP is a ways out, but this will require a significant amount of work.

Regarding your above code:

• Your first constraint would need to be a @NLconstraint in JuMP anyway since @constraint only supports linear and quadratic expressions (not min and max). However, the max operator can be reformulated via traditional linear programming reformulations (i.e., introducing an auxiliary variable with inequality constraints).
• Your indicator function could be approximated via big-M constraints (i.e., using binary variables) and thus avoid the need for nonlinear expressions.

If we want the constraint Prob(f(ξ) <= 0) <= δ (a chance constraint where f denotes the RHS function) we can equivalently represent it Expectation[indicator_(f(ξ) <= 0)(ξ)] <= δ. This can in turn be represented by big-M constraints (introducing an infinite variable y(ξ) in {0,1} for the indicator function and a sufficiently large scalar constant M). For example, if f(ξ) := w(ξ) - c * x(ξ) then in InfiniteOpt.jl we can do the following:

using JuMP, InfiniteOpt, Distributions
c = 42 # whatever it is 
δ = 0.95 # whatever it is
M = 1000 # a sufficient upper bound to relax the constraints
m = InfiniteModel()
@infinite_parameter(m, ξ in Normal(0, 1), num_supports = 1000)
@infinite_variable(m, w(ξ))
@infinite_variable(m, x(ξ))
@infinite_variable(m, y(ξ), Bin) # for the indicator function
@constraint(m, big_m_constr, w - c * x <= y * M)
@constraint(m, probabilistic_constr, expect(1 - y, ξ) <= δ)

[charles-r-smith]

How would you do the first constraint with the min/max. I tried defining boundary conditions, but that kept erroring (I think it's because I can only define boundary conditions with respect to the infinite parameter, not with respect to other variables).

for i = 1:N
   @BDconstraint(m, ((transpose(A)*p+c-x)[i] in [-1000000.0,-1e-16]), p[i] == 0) # transpose(A) * p + c -x <= 0 
   @BDconstraint(m, ((transpose(A)*p+c-x)[i] in [0.0, p_bar[i]]), p[i] == (transpose(A)*p+c-x)[i])
 @BDconstraint(m, ((transpose(A)*p+c-x)[i] in [p_bar[i]+1e-9, 1000000.0]), p[i] == p_bar[i])
end

For the probability aspect, when I incorporate why you suggested, I get the following error:

MathOptInterface.UnsupportedConstraint{MathOptInterface.SingleVariable,MathOptInterface.ZeroOne}: MathOptInterface.SingleVariable-in-MathOptInterface.ZeroOne constraint is not supported by the model.

Any thoughts? Thank you very much!

[pulsipher]

With respect to the use of @BDconstraint it corresponds to bounded constraints which are constraints over some sub-domain of the overall infinite domain which is parameterized by infinite parameters. Thus, the second argument specifies the sub-domain in terms of the infinite parameter(s) and the third argument specifies the constraint that will be enforced over that sub-domain. For example, consider the following time valued illustration:

using JuMP, InfiniteModel
m = InfiniteModel()
@infinite_parameter(m, t in [0, 10]) # overall time domain (infinite domain) is the interval [0, 10]
@infinite_variable(m, x(t)) # define some time variable

@constraint(m, x <= 42) # this will be enforced over the entirety of the infinite domain (i.e., t in [0, 10])
@BDconstraint(m, t in [0, 10], x <= 42) # equivalent to above with @constraint since [0, 10] is the full domain
@BDconstraint(m, t == 0, x ==5) # this will only be enforced over the subdomain t in [0, 0] (a subset of [0, 10])
@BDconstraint(m, t in [2, 5], x^2 <= 9) # this constraint will only be enforced over the subdomain t in [2, 5]

Thus, @BDconstraint exists to limit the scope of a particular constraint with respect to its infinite domain dependencies.

With respect to your min-max constraint, my intuition is that it can probably be reformulated into a more convenient form when the entire optimization problem is considered. I am not exactly sure what overall formulation you are tackling, but I am guessing there might be a way to drop the min operator and then only have to deal with the max operator which is straightforward to reformulate if the objective is minimizing. I don't quite follow the reformulation you are attempting above, but the error is due to the nature of what bounded constraints are in accordance with my conversation above. If you can get it in a form with only the max operator then we can reformulate it as follows:

using JuMP, InfiniteOpt, Distributions
m = InfiniteModel()
@infinite_parameter(m, ξ in Normal(0, 1), num_supports = 1000)
@infinite_variable(m, x(ξ)) 

# Enforce the constraint max(42, x) <= 80
@infinite_variable(m, z(ξ)) # auxiliary variable to represent max(42, x)
@constraint(m, z >= 42) # helps characterize the auxiliary variable
@constraint(m, z >= x) # helps characterize the auxiliary variable
@constraint(m, z <= 80) # defines the constraint we wanted with the max operator

Finally, the MOI error is due to your choice of optimizer (solver). This is because the big-M constraint reformulation uses binary variables and you'll need to use an appropriate mixed-integer solver. Depending on the rest of the model, you could probably use Gurobi, Cplex, or Cbc.

[charles-r-smith]

I see what you are saying. I guess my problem is a little more involved than the example and for that reason I'm unsure how to connect the two. Mainly because my constraint is w == max(42,x) where w is a variable in my model rather than a constant. I don't see how you can define that using the auxiliary variable method.

[pulsipher]

That is not a problem at all. To formulate the constraint w(ξ) = max(42, x(ξ)) we do the following:

using JuMP, InfiniteOpt, Distributions
m = InfiniteModel()
@infinite_parameter(m, ξ in Normal(0, 1), num_supports = 1000)
@infinite_variable(m, x(ξ)) 
@infinite_variable(m, w(ξ))

# Enforce the equivalent set of constraints
@infinite_variable(m, z(ξ)) # auxiliary variable to represent max(42, x)
@constraint(m, z >= 42) # helps characterize the auxiliary variable
@constraint(m, z >= x) # helps characterize the auxiliary variable
@constraint(m, w == z) # defines the constraint we wanted with the max operator

Since w = z, this can be further simplified to:

using JuMP, InfiniteOpt, Distributions
m = InfiniteModel()
@infinite_parameter(m, ξ in Normal(0, 1), num_supports = 1000)
@infinite_variable(m, x(ξ)) 
@infinite_variable(m, w(ξ))

# Enforce the equivalent set of constraints
@constraint(m, w >= 42)
@constraint(m, w >= x)

Note that I am making some assumptions about the objective function (namely that it is minimizing). These types of reformations stem from considering min-max problems in linear programming. For more information refer to https://math.stackexchange.com/questions/2446606/linear-programming-set-a-variable-the-max-between-two-another-variables and http://apmonitor.com/me575/index.php/Main/MiniMax.

[charles-r-smith]

using JuMP, InfiniteOpt, Distributions
m = InfiniteModel()
@infinite_parameter(m, ξ in Normal(0, 1), num_supports = 1000)
@infinite_variable(m, x(ξ)) 
@infinite_variable(m, w(ξ))

# Enforce the equivalent set of constraints
@infinite_variable(m, z(ξ)) # auxiliary variable to represent max(42, x)
@constraint(m, z >= 42) # helps characterize the auxiliary variable
@constraint(m, z >= x) # helps characterize the auxiliary variable
@constraint(m, w == z) # defines the constraint we wanted with the max operator

There are many values of z that satisfy the constraints z >= 42 and z >= x. One of them is z = max(42,x). However, if x is 50, the value z = 60 would satisfy those constraints. To me this doesn't ensure that z = max(42,x), just that z is greater than both 42 and x. Am I missing something?

[pulsipher]

This reformulation depends on the objective function. If the objective will tend for z to be minimized then this reformulation is exact. Otherwise, it can be reformulated by introducing a binary variable. Again, this is all better explained in detail here https://math.stackexchange.com/questions/2446606/linear-programming-set-a-variable-the-max-between-two-another-variables (you can also refer to online linear programming course materials like https://laurentlessard.com/teaching/524-intro-to-optimization/#:~:text=This%20course%20is%20an%20introduction,in%20industry%20and%20research%20applications.). Without knowing your full formulation, I cannot determine which reformulation will suit your situation.

[charles-r-smith]

Ah, this makes so much more sense. Thank you.

[pulsipher]

@smitch151 If there are no further questions, would you mind me closing this issue?

[pulsipher]

Seeing that the immediate issue is resolved, I will close this issue. Further discussion about the development of the NLP interface can be made in #16.

[charles-r-smith]

I'm trying to revisit in a simple case. Let's say I want to have the constraint p == min(A*x, p_bar).

From my understanding of the Big M stuff you described, I need four constraints:

1. p <= A*x
2. p <= p_bar
3. p >= Ax - Mind
4. p >= p_bar - M*ind

where ind = 1 when p_bar < A*x.

When I try to code this up I get the following error: "MethodError: no method matching isless(::GenericAffExpr{Float64,GeneralVariableRef}, ::Float64)"

How can I generate this indicator to do this comparison.

using Pkg
Pkg.add(Pkg.PackageSpec(;name="InfiniteOpt", version="0.3.2"))

using LinearAlgebra, JuMP, Ipopt, Distributions, Random
using InfiniteOpt

N = 2
p_bar = [5.0, 5.0]
A = [1.0 1.0; 2.0 2.0]

# Setting up model 
m = InfiniteModel(Ipopt.Optimizer);

# Parameter everything depends upon (shocks to each node)
@infinite_parameter(m, x[i=1:N] in Beta(1.0,50.0), num_supports = 10000) #parameter of shocks

# setting up variables, starting positino, and constraints
@infinite_variable(m, p[i=1:N](x), start = 0) # clearing vector. 

# setting up constraint
#@constraint(m, p == min(A*x, p_bar))
M = 100
@infinite_variable(m, ind[i=1:N](x), Bin) #ind indicates if p_bar>x
ind .== (p_bar .>= A * x)
@constraint(m, p .<= A*x)
@constriant(m, p .<= p_bar)
@constraint(m, p .>= A* x - M.*ind)
@constraint(m, p .>= p_bar - M.*(1-ind))


@objective(m, Min, expect(sum(x.*p),x)) # min/max E[(ci * xi + p_bar i - pi(x)))]

# Training Model
@time optimize!(m);

Once again, thank you.

[pulsipher]

The error is occurring since the comparison syntax p_bar .>= A * x is only valid inside @constraint to define constraints as per typical JuMP syntax.

Your reformulation doesn't appear to be correct. Would you please algebraically define the optimization problem you are trying to solve before any reformulation is applied to it and also label each symbol as a constant, infinite parameter, or infinite variable (perhaps this could be provided via a short attachment). Once, I better understand the full problem you are trying to solve, I should be better able to help you. This sounds like an interesting application class.