[Help Wanted] NL Constraints

[CharlesRSmith44]

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: https://github.com/pulsipher/InfiniteOpt.jl/issues/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:

1. @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 https://github.com/jump-dev/JuMP.jl/issues/2402). JuMP NLP won't be extendable until it is fundamentally redone via https://github.com/jump-dev/MathOptInterface.jl/issues/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, ξ) <= δ)

[CharlesRSmith44]

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.

[CharlesRSmith44]

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.

[CharlesRSmith44]

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.

[CharlesRSmith44]

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.

[CharlesRSmith44]

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.

[CharlesRSmith44]

Okay I'll try to explain more fully.

We want to maximize the expected value of p.x subject to the constraint that p == min(Ax, p_bar) where x is drawn from a beta distribution. The parameters of the beta distribution are fixed in this case. The parameters A and p_bar are allowed to vary.

The complication comes from the fact that we can't just use a minimum function in infinite opt, so we need to use big M to get around it. The big M specification of C == min(A,B) is:

1. C <= A
2. C <= B (this ensures that C is smaller than the smallest element of A and B. However C could still be -1000 if A and B == 0 which isn't right. So, we create and indicator which equals 1 when B >= A and 0 otherwise.
3. C >= A - M*(1-ind)
4. C >= B - M ind Now, if B>= A we need C == A as our solution. This occurs because in that case ind == 1 so C>= A - M(1-ind) = A. If B < A, we need C == B as our solution. This occurs because in that case ind == 0. so C >= B- M*ind = B.

I updated the code your @constraint comment. I think I have the formulation of the big M aspect right in the code.

However, I still get the error "MethodError: no method matching isless(::GenericAffExpr{Float64,GeneralVariableRef}, ::Float64)" when I do the comparison in the @constraint.

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, 2.0]
A = [3.0 10.0; 10.0 3.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
@constraint(m, ind .== (p_bar .>= A*x))
@constraint(m, p .<= A*x)
@constraint(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[x.*p]

# Training Model
@time optimize!(m);

[CharlesRSmith44]

If I create an arbitrary variable w, with the same dimensions as x, the code works just fine. Also, I did need to switch the indicator direction you were correct about that. Good catch. Is there anyway to do this with x the infinite parameter rather than just some vector?

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

N = 2
p_bar = [5.0, 5.0]
A = [3.0 10.0; 10.0 3.0]
w = [0.0; 1.0]

# Setting up model 
m = InfiniteModel(Cbc.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
@constraint(m, ind .== (p_bar .>= A*w))
@constraint(m, p .<= A*w)
@constraint(m, p .<= p_bar)
@constraint(m, p .>= A*w .- M .* (1 .- ind))
@constraint(m, p .>= p_bar .- M.*ind)

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

# Training Model
@time optimize!(m);

'''

[CharlesRSmith44]

Also this issue occurs if x is an infinite variable, not just when x is an infinite parameter. For the problem I'm trying to solve x should be an infinite parameter, but just wondering if that helps you understand it.

[pulsipher]

First, the syntax @constraint(m, ind .== (p_bar .>= A*x)) is incorrect since @constraint is not meant to contain more than one operator, but this contains 2 (.== and .>=). The reason it happens to work with w is that (p_bar .>= A*w) evaluates to an array of booleans since it contains no model variables or parameters. We do not currently provide a native symbolic indicator function which is what I think you are trying to accomplish here. But one can incorporate big-M based reformulations to handle indicator functions via reformulation as I showed before.

Let's return to problem you want to solve. If I understand correctly, your desired formulation is: max 𝔼[p(x)' x] s.t. p(x) = min(Ax, p_bar), ∀ x ∈ X (where X is the co-domain of the beta distribution) where p(x) is a vector of infinite variables, x is a vector of beta random infinite parameters, A is matrix of constants, and p_bar is a vector of constants. This actually doesn't need to be posed as an optimization problem since there are no degrees of freedom (i.e., p(x) is constrained to one possible form by the constraint).

Let's first ignore the last observation and model this as an optimization problem. This actually doesn't require the introduction of big-M constraints to handle the min because none of its inputs are decision variables, it is just a function of infinite parameters x. This means that once we generate samples for x the min function will evaluate to a constant (not containing any decision variables). Thus, we can employ an infinite_parameter_function from InfiniteOpt to represent it (see https://pulsipher.github.io/InfiniteOpt.jl/stable/guide/expression/#Infinite-Parameter-Functions). With this observation, let's formulate the problem in InfiniteOpt!

using InfiniteOpt, JuMP, Clp, Distributions, Random, LinearAlgebra
using JuMP.Containers: @container # convenient for infinite parameter functions

Random.seed!(42) # seed to reproduce same samples for the sake of example

N = 2
p_bar = [5.0, 5.0]
A = [3.0 10.0; 10.0 3.0]

# Initialize the  model 
m = InfiniteModel(Clp.Optimizer)

# Define the random parameters
@infinite_parameter(m, x[i=1:N] in Beta(1.0, 50.0), num_supports = 1000)

# Define the decision variables
@infinite_variable(m, p[i=1:N](x), start = 0)

# Define the infinite parameter function for each element of min(Ax, p_bar)
@container(my_func[i = 1:N], parameter_function(x_sample -> min((A * x_sample)[i], p_bar[i]), x))

# Define the objective
@objective(m, Max, expect(p' * x, x))

# Define the constraints
@constraint(m, p .== my_func)

# Optimize the model (i.e., discretize and solve the model)
optimize!(m)

# Get results
println("Objective value: ", objective_value(m))
println("Values of p(x): ", value.(p))

Notice in the above that p(x) = my_func, so we can equivalently formulate:

using InfiniteOpt, JuMP, Clp, Distributions, Random, LinearAlgebra
using JuMP.Containers: @container # convenient for infinite parameter functions

Random.seed!(42) # seed to reproduce same samples for the sake of example

N = 2
p_bar = [5.0, 5.0]
A = [3.0 10.0; 10.0 3.0]

# Initialize the  model 
m = InfiniteModel(Clp.Optimizer)

# Define the random parameters
@infinite_parameter(m, x[i=1:N] in Beta(1.0, 50.0), num_supports = 1000) 

# Define the infinite parameter function for each element of min(Ax, p_bar)
@container(my_func[i = 1:N], parameter_function(x_sample -> min((A * x_sample)[i], p_bar[i]), x))

# Define the objective
@objective(m, Max, expect(my_func' * x, x))

# Optimize the model (i.e., discretize and solve the model)
optimize!(m)

# Get results
println("Objective value: ", objective_value(m))

Now we readily see that our model contains no decision variables which backs up our original observation that this need not be an optimization problem since it has no degrees of freedom. Thus, we can forego the optimization and just compute the expectation directly using sample average approximation to obtain the same answer.

using Distributions, LinearAlgebra, Random

Random.seed!(42)

N = 2
S = 1000 # number of samples
p_bar = [5.0, 5.0]
A = [3.0 10.0; 10.0 3.0]

# Define the distribution and get the samples 
dist = Beta(1.0, 50.0)
x = hcat([rand(dist, S) for i = 1:N]...)

# Compute the expectation via sample average approximation
my_func(x_sample) = [min((A * x_sample)[i], p_bar[i]) for i = 1:N]
objective = 1 / S * sum(my_func(x[s, :])' * x[s, :]  for s = 1:S) 

# Get results
println("Objective value: ", objective)

I hope that helps.

[CharlesRSmith44]

This is fantastic. Thank you so much. One, hopefully final, question:

What's the syntax if I want to make A a variable in my model, rather than some external variable, but still want the same container function?

i,e.

using InfiniteOpt, JuMP, Clp, Distributions, Random, LinearAlgebra
using JuMP.Containers: @container # convenient for infinite parameter functions

Random.seed!(42) # seed to reproduce same samples for the sake of example

N = 2
p_bar = [5.0, 5.0]
A = [3.0 10.0; 10.0 3.0]

# Initialize the  model 
m = InfiniteModel(Clp.Optimizer)

# Define the random parameters
@infinite_parameter(m, x[i=1:N] in Beta(1.0, 50.0), num_supports = 1000) 

# Define A as a variable in the model: 
@hold_variable(m, A[i=1:N, j=1:N] >= 0)

# Define the infinite parameter function for each element of min(Ax, p_bar)
@container(my_func[i = 1:N], parameter_function(x_sample -> min((A * x_sample)[i], p_bar[i]), x))

[pulsipher]

Glad to be of help.

First, as a clarification the @container macro is just a convenient tool for constructing arrays and has nothing specifically to do with infinite parameter functions.

using JuMP.Containers: @container

# Make a vector `a` where each element is `i + 2`
@container(a[i = 1:10], i + 2)

# The above is equivalent to 
a = Vector{Int}(undef, 10) # preallocate the array 
for i = 1:10
    a[i] = i + 2
end

Now let's consider you problem where the A matrix is considered to be a nonnegative matrix of finite decision variables (i.e., first-stage variables in this context).

max  𝔼_x[min(Ax, p_bar)'x]
 A
s.t. A >= 0

Here, we can no longer use parameter_function in InfiniteOpt to handle min(Ax, p_par) since it now contains decision variables A (infinite parameter functions can only contain infinite parameters and constants). However, the good news is this formulation amounts to a variant of a maximin problem which means we can reformulate the min(Ax, p_par) by introducing a continuous variable z(x) and adding some inequality constraints:

 max   𝔼_x[z(x)'x]
A,z(x)
 s.t.  z(x) <= Ax,    ∀ x ∈ X
       z(x) <= p_bar, ∀ x ∈ X
       A >= 0

This reformulation is exact since the objective wants to make z(x) as large as possible, but this is limited by Ax and p_bar (whichever is smaller). The other piece of good news is that this formulation is linear! Thus, we can readily implement this in InfiniteOpt:

using InfiniteOpt, JuMP, Clp, Distributions, LinearAlgebra

# Define the constants
N = 2
p_bar = [5.0, 5.0]

# Initialize the  model 
m = InfiniteModel(Clp.Optimizer)

# Define the random parameters
@infinite_parameter(m, x[1:N] in Beta(1.0, 50.0), num_supports = 1000) 

# Define A as a matrix of variables: 
@hold_variable(m, A[1:N, 1:N] >= 0)

# Define the reformulation infinite variables z(x)
@infinite_variable(m, z[1:N](x))

# Define the objective
@objective(m, Max, expect(z' * x, x))

# Define the constraints 
@constraint(m, z .<= A * x)
@constraint(m, z .<= p_bar)

# Solve the model and get the results
optimize!(m)
opt_objective = objective_value(m)
opt_A = value.(A)

Hope that helps!

[CharlesRSmith44]

Thank you. This is all very helpful. I'm beginning to understand it. What if we are trying to minimize, so that "this reformulation is exact since the objective wants to make z(x) as large as possible, but this is limited by Ax and p_bar (whichever is smaller)." is not satisfied.

i.e.

using InfiniteOpt, JuMP, Clp, Distributions, LinearAlgebra

# Define the constants
N = 2
p_bar = [5.0, 5.0]

# Initialize the  model 
m = InfiniteModel(Clp.Optimizer)

# Define the random parameters
@infinite_parameter(m, x[1:N] in Beta(1.0, 50.0), num_supports = 1000) 

# Define A as a matrix of variables: 
@hold_variable(m, A[1:N, 1:N] >= 0)

# Define the reformulation infinite variables z(x)
@infinite_variable(m, z[1:N](x))

# Define the objective
@objective(m, Min, expect(z' * x, x))

# Define the constraints 
#@constraint(m, z == min(A*x, p_bar))
@constraint(m, z .<= A * x)
@constraint(m, z .<= p_bar)

# Solve the model and get the results
optimize!(m)
opt_objective = objective_value(m)
opt_A = value.(A)

The optimizer runs, but it says the problem is unbounded as expected. Essentially, is there a way to generate an exact minimum function using constraints involving decision variables?

[pulsipher]

The previous reformulation indeed is no longer exact for a minimization problem. For such a problem, we need to introduce binary variables y(x) with a sufficiently large constant M:

   min     𝔼_x[z(x)'x]
A,z(x),y(x)
   s.t.    z(x) <= Ax,                  ∀ x ∈ X
           z(x) <= p_bar,               ∀ x ∈ X
           z(x) >= Ax - My(x),          ∀ x ∈ X
           z(x) >= p_bar - M(1 - y(x)), ∀ x ∈ X
           y(x) ∈ {0, 1},               ∀ x ∈ X
           A >= 0

So in InfiniteOpt we have:

using InfiniteOpt, JuMP, Cbc, Distributions, LinearAlgebra

# Define the constants
N = 2
p_bar = [5.0, 5.0]
M = 1000 # replace with an appropriate value

# Initialize the  model 
m = InfiniteModel(Cbc.Optimizer) # Gurobi would be much faster

# Define the random parameters
@infinite_parameter(m, x[1:N] in Beta(1.0, 50.0), num_supports = 1000) 

# Define A as a matrix of variables: 
@hold_variable(m, A[1:N, 1:N] >= 0)

# Define the reformulation infinite variables z(x) and y(x)
@infinite_variable(m, z[1:N](x))
@infinite_variable(m, y[1:N](x), Bin)

# Define the objective
@objective(m, Max, expect(z' * x, x))

# Define the constraints 
@constraint(m, z .<= A * x)
@constraint(m, z .<= p_bar)
@constraint(m, z .>= A * x - M * y)
@constraint(m, z .>= p_bar - M * (1 .- y))

# Solve the model and get the results
optimize!(m)
opt_objective = objective_value(m)
opt_A = value.(A)

Note that we now have a MILP and that it will scale exponentially in complexity with num_supports due to the creation of more and more binary variables (so this may become intractable as we use a large amount of supports). Also, this problem will be infeasible if M is too small, but will get harder and harder to solve for excessively large values of M.

All of these reformulations pertain to optimization theory fundamentals and it might be worth it to review some optimization modeling fundamentals such as those covered in the slides from this optimization modeling course: https://laurentlessard.com/teaching/524-intro-to-optimization/#:~:text=This%20course%20is%20an%20introduction,in%20industry%20and%20research%20applications

Thus, for future reference I would recommend first referring to those resources to better focus questions asked here with concerns specific to InfiniteOpt and limit those that pertain to general optimization modeling questions.

I hope that helps and thank you for your interest so far in using InfiniteOpt! We are genuinely very grateful for you and all of our users that are using a tool we have put a lot of time into developing to make a difference.

[pulsipher]

@smitch151 Did the above reformulations work for the problem you're working on?

[CharlesRSmith44]

Yes that was perfect. Thank you!

[pulsipher]

Perfect, I will go ahead and close this issue then. Please reach out again in the future with any other questions you have. Happy modeling!

[pulsipher]

As a follow-up, full nonlinear support beyond what JuMP does was added in v0.5.0 via #161.