Hello! I was searching to transition some of my deterministic ILP models to a stochastic version.

I started investigating InfiniteOpt.jl and made a simple stochastic flavor of the JuMP LP example.

I just modeled a parameter of the constraint to be a Normal around the initial value, i.e. `ξ ~ Normal(6,1)`

instead of scalar `ξ = 6`

. I was expecting the result to be the same since the expectation `E(ξ) = 6`

and for this problem I think the Value of the Stochastic Solution (VSS) should be zero. However InfiniteOpt.jl returns somewhat random results.

Do I need to somehow increase the granularity for the transition of an infinite problem description to the finite version? If yes, how can I do it ? I tried playing around with the `num_supports`

argument but didn’t have any luck.

MWE:

The deterministic initial LP program:

```
using JuMP, Graphs, InfiniteOpt, Distributions, HiGHS
# Deterministic problem according to https://jump.dev/JuMP.jl/stable/tutorials/getting_started/getting_started_with_JuMP/#An-example
function modeldet()
model = Model(HiGHS.Optimizer)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
ξ = 6
@constraint(model, c1, ξ * x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
@objective(model, Min, 12x + 20y)
optimize!(model)
return model
end
```

and the InfiniteOpt.jl version:

```
# Introduce some uncertainty in the parameters using probability density functions
function modelinf()
model = InfiniteModel(HiGHS.Optimizer)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
@infinite_parameter(model, ξ ~ Normal(6,1))
@constraint(model, c1, ξ * x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
@objective(model, Min, 12x + 20y)
optimize!(model)
return model
end
```

Plotting the two problem constraints gives a clear intuition where the optimum decision variable should be. The green constraint is stochastic and parametrized by `ξ`

.

## Makie code for plotting

```
using GLMakie
function plotprob()
f = Figure()
a = Axis(f[1,1])
for ξ in rand(Normal(6,1), 50)
plot!(a, -100:100, x -> (100 - ξ * x)/8, color = :green, linewidth=1)
end
plot!(a, -100:100, x -> (120 - 7x)/12, color = :blue, linewidth=5)
# objective value
# for ob in 0:50:500
# plot!(a, 0:30, x -> (ob - 12x)/20, color = :black)
# end
return f
end
```

The deterministic `JuMP`

program finds the the optimum solution:

```
julia> md = modeldet();
Running HiGHS 1.5.3 [date: 1970-01-01, git hash: 45a127b78]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
2 rows, 2 cols, 4 nonzeros
2 rows, 2 cols, 4 nonzeros
Presolve : Reductions: rows 2(-0); columns 2(-0); elements 4(-0) - Not reduced
Problem not reduced by presolve: solving the LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Pr: 2(220) 0s
2 2.0500000000e+02 Pr: 0(0) 0s
Model status : Optimal
Simplex iterations: 2
Objective value : 2.0500000000e+02
HiGHS run time : 0.00
julia> value(md[:x])
15.000000000000005
```

But the InfiniteOpt.jl struggles to find the same

```
julia> mdi = modelinf();
Running HiGHS 1.5.3 [date: 1970-01-01, git hash: 45a127b78]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
2 rows, 2 cols, 4 nonzeros
2 rows, 2 cols, 4 nonzeros
Presolve : Reductions: rows 2(-9); columns 2(-0); elements 4(-18)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Pr: 2(200) 0s
2 2.7496874512e+02 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model status : Optimal
Simplex iterations: 2
Objective value : 2.7496874512e+02
HiGHS run time : 0.00
julia> value(mdi[:x])
17.91406209370758
```

I also get high variability, which appear to be consistently above the optimum

```
julia> scatter([value(modelinf()[:x]) for _ in 1:40])
```

So my question is how can I write the stochastic program in order to get a reliable solution ?

### Feeding in samples of a Random Variable instead of the distribution ?

A bit of irrelevant but I was also wondering whether I can feed into the model samples instead of a well defined distribution. The idea is that I am estimating the parameters using MCMC Bayesian modeling, thus I end up having a chain of values that describe the parameter.

Ofc, I could try to simplify the chain into some kind of discrete distribution, but it feels like going from finite samples to an infinite PDF just to go back to finite problem (through InfiniteOpt.jl transition, i.e. finite samples) might be a waste of accuracy (?)