In this chapter we show how to solve an optimization problem for many scenarios, where each scenario has slightly different input data. We also consider the generation of random scenarios.
5.1 Scenario analysis
Consider again the Duplo example with objective function coefficients s_1 and s_2:
\begin{array}{llll}
\text{maximize}& s_1\cdot x_1 + s_2 \cdot x_2 && \text{(profit)} \\
\text{subject to} &x_1 + 2x_2 \leq 6 &&\text{(big bricks)} \\
&2x_1 + 2x_2 \leq 8 &&\text{(small bricks)} \\
&x_1, x_2 \geq 0
\end{array}
with decision variables
x_1: Number of chairs
x_2: Number of tables
We can solve the Duplo example for different objective function coefficients, with the results as specified in the table below.
Scenario
s_1
s_2
Optimal solution
1
15
20
(2,2)
2
9
25
(0,3)
3
20
19
(4,0)
4
16
21
(2,2)
To solve the Duplo problem repeatedly in Python using a different scenario every time, we can use a for-loop. We specify the different values of s_1 and s_2 in two separate lists, where s1[i] is the value of s_1 in scenario i+1 (recall Python starts indexing at zero), and similarly s2[i] the value of s_2 in scenario i+1.
from gurobipy import Model, GRB# Different scenarioss1 = [15, 9, 20, 16]s2 = [20, 25, 19, 21]for i inrange(0,4):# Initialize Gurobi model. duplo = Model(name='Duplo problem')# Suppress Gurobi Optimizer output duplo.setParam('OutputFlag',0)# Declare the two decision variables. x1 = duplo.addVar(name='chairs') x2 = duplo.addVar(name='tables')# Specify the objective function in terms of s1[i], s2[i] duplo.setObjective(s1[i]*x1 + s2[i]*x2, sense=GRB.MAXIMIZE)# Add the resource constraints on bricks. duplo.addConstr(x1 +2*x2 <=6, name='big-bricks') duplo.addConstr(2*x1 +2*x2 <=8, name='small-bricks')# Optimize the duplo duplo.optimize()# Print solutionprint("Optimal solution for scenario", i+1, "is [", x1.X, ",", x2.X, "]")
Set parameter Username
Set parameter LicenseID to value 2715549
Academic license - for non-commercial use only - expires 2026-09-29
Optimal solution for scenario 1 is [ 2.0 , 2.0 ]
Optimal solution for scenario 2 is [ 0.0 , 3.0 ]
Optimal solution for scenario 3 is [ 4.0 , 0.0 ]
Optimal solution for scenario 4 is [ 2.0 , 2.0 ]
A couple of remarks are in place for the code above. First of all, note that all the lines below the line for i in range(0,4): are indented because we want the model to be constructed, optimized and printed in every iteration.
We use the command duplo.setParam('OutputFlag',0) to suppress the output that Gurobi usually generates when we optimize a linear optimization problem. We do this to avoid that this output is printed in every interation of the for-loop.
Also note that the objective function is defined as s1[i]*x1 + s2[i]*x2 meaning that in iteration i, the objective function uses the values s1[i] and s2[i] (which are therefore different in every iteration as they depend on i).
In the print() statement we print a combination of text and the values i+1, x1.X and x2.X where the latter two are the values of the optimal solution in scenario i+1.
5.2 Random scenarios
In the previous section, we chose the values of the objective function coefficients in every scenario ourselves. These can also be generated randomly using random number generation functions in Python contained in the NumPy package, which contains a lot of the mathematical functionality that Python has to offer.
This package is typically imported under the aliasnp, which is shorthand notation so that we do not have to type numpy everytime.
import numpy as np
We treat two functions for generating random numbers. The first one can be used to generate fractional random numbers.
np.random.rand(n) : Generates n random numbers from the interval [0,1].
The generated numbers are stored in a list. Technically speaking, they are stored in a NumPy array, which you can think of as a “mathematical” list (but this is not important for us).
import numpy as np# Generate five random (fractional) numbers from [0,1]x = np.random.rand(5)# Print the numbersprint(x)
Note that if you want to generate a (fractional) number from an interval [a,b], you can do this by generating a number x from the interval [0,1], multiplying this number by (b-a) and then adding a to the result.
In other words, if x is a random number from [0,1], then the number a + (b-a) \cdot x is a random number from [a,b].
import numpy as np# Interval valuesa =3b =7# Generate random (fractional) number from [0,1]x = np.random.rand(1)# Random number from interval [a,b] = [3,7]y = a + (b-a)*x# Random number from interval [a,b]print(y)
[5.14693509]
Note that, in the above code snippet, even though we only generate one number, it is still given to us in a list. If you only want to get the number (not in a list) you should index the list at position 0. This is important if you want to use such a number for defining the objective function or a constraint.
print(y[0])
5.146935088585712
You can also generate multiple numbers from an interval [a,b] if x = [x_0, x_1, \dots] is a list of numbers generated from [0,1]. If you multiply x with (b-a), then all numbers in the list are multiplied by (b-a), i.e., (b-a)*x = [(b-a)*x_0, (b-a)*x_1,\dots]. Similarly, if we add a to (b-a)*x, then a is added to every number in the list (b-a)*x, i.e., we get a + (b-a)*x = [a + (b-a)*x_0, a + (b-a)*x_1, \dots].
import numpy as np# Interval valuesa =3b =7# Generate five random (fractional) numbers from [0,1]x = np.random.rand(5)# Five random numbers from interval [a,b] = [3,7]y = a + (b-a)*x# Print the numbersprint(y)
The second function we discuss can generate integer-valued numbers from a specified range of integers:
np.random.randint(low=a, high=b, size=n) : Generates n random integers from the set \{a, a+1, \dots, b-1\}.
import numpy as np# Generate five random integers from the set {2,3,4,5,6,7,8,9}y = np.random.randint(low=2, high=10, size=5)# Print the numbersprint(y)
[9 7 6 5 6]
We can use these random numbers to generate random scenarios and then solve the Duplo example for those scenarios.
from gurobipy import Model, GRB# Different 20 random scenarioss1 = np.random.randint(low=15, high=25, size=20) # Numbers from {15,16,...,24}s2 = np.random.randint(low=20, high=35, size=20) # Numbers from {20,12,...,34}for i inrange(0,20):# Initialize Gurobi model. duplo = Model(name='Duplo problem')# Suppress Gurobi Optimizer output duplo.setParam('OutputFlag',0)# Declare the two decision variables. x1 = duplo.addVar(name='chairs') x2 = duplo.addVar(name='tables')# Specify the objective function in terms of s1[i], s2[i] duplo.setObjective(s1[i]*x1 + s2[i]*x2, sense=GRB.MAXIMIZE)# Add the resource constraints on bricks. duplo.addConstr(x1 +2*x2 <=6, name='big-bricks') duplo.addConstr(2*x1 +2*x2 <=8, name='small-bricks')# Optimize the duplo duplo.optimize()# Print solutionprint("Optimal solution for scenario", i+1, "is [", x1.X, ",", x2.X, "]")
Optimal solution for scenario 1 is [ 2.0 , 2.0 ]
Optimal solution for scenario 2 is [ 4.0 , 0.0 ]
Optimal solution for scenario 3 is [ 2.0 , 2.0 ]
Optimal solution for scenario 4 is [ 2.0 , 2.0 ]
Optimal solution for scenario 5 is [ 2.0 , 2.0 ]
Optimal solution for scenario 6 is [ 2.0 , 2.0 ]
Optimal solution for scenario 7 is [ 2.0 , 2.0 ]
Optimal solution for scenario 8 is [ 2.0 , 2.0 ]
Optimal solution for scenario 9 is [ 2.0 , 2.0 ]
Optimal solution for scenario 10 is [ 2.0 , 2.0 ]
Optimal solution for scenario 11 is [ 2.0 , 2.0 ]
Optimal solution for scenario 12 is [ 2.0 , 2.0 ]
Optimal solution for scenario 13 is [ 2.0 , 2.0 ]
Optimal solution for scenario 14 is [ 2.0 , 2.0 ]
Optimal solution for scenario 15 is [ 2.0 , 2.0 ]
Optimal solution for scenario 16 is [ 2.0 , 2.0 ]
Optimal solution for scenario 17 is [ 2.0 , 2.0 ]
Optimal solution for scenario 18 is [ 2.0 , 2.0 ]
Optimal solution for scenario 19 is [ 2.0 , 2.0 ]
Optimal solution for scenario 20 is [ 2.0 , 2.0 ]
Note that the solution (2,2) appears very often, meaning it is a stable solution in the sense that even though the objective function coefficients change, the solution remains optimal.
You can also generate the random numbers inside of the for-loop. Note that in this case you should use s1[0] and s2[0] in the objective function. As discussed before, this is because s1 = np.random.randint(low=15, high=25, size=1) is a list of length 1, and so to get the number in this list, you have to index s1 at index 0.
Note that in every iteration, first two random numbers s1 and s2 are generated, which are then later used to define the objective function. Because new numbers are generated at the start of every iteration, we are solving a different problem every time.
from gurobipy import Model, GRBfor i inrange(0,20):# Generate random coefficients s1 = np.random.randint(low=15, high=25, size=1) s2 = np.random.randint(low=15, high=25, size=1)# Initialize Gurobi model. duplo = Model(name='Duplo problem')# Suppress Gurobi Optimizer output duplo.setParam('OutputFlag',0)# Declare the two decision variables. x1 = duplo.addVar(name='chairs') x2 = duplo.addVar(name='tables')# Specify the objective function in terms of s1[0], s2[0] duplo.setObjective(s1[0]*x1 + s2[0]*x2, sense=GRB.MAXIMIZE)# Add the resource constraints on bricks. duplo.addConstr(x1 +2*x2 <=6, name='big-bricks') duplo.addConstr(2*x1 +2*x2 <=8, name='small-bricks')# Optimize the duplo duplo.optimize()# Print solutionprint("Optimal solution for scenario", i+1, "is [", x1.X, ",", x2.X, "]")
Optimal solution for scenario 1 is [ 4.0 , 0.0 ]
Optimal solution for scenario 2 is [ 4.0 , 0.0 ]
Optimal solution for scenario 3 is [ 4.0 , 0.0 ]
Optimal solution for scenario 4 is [ 2.0 , 2.0 ]
Optimal solution for scenario 5 is [ 2.0 , 2.0 ]
Optimal solution for scenario 6 is [ 2.0 , 2.0 ]
Optimal solution for scenario 7 is [ 2.0 , 2.0 ]
Optimal solution for scenario 8 is [ 2.0 , 2.0 ]
Optimal solution for scenario 9 is [ 2.0 , 2.0 ]
Optimal solution for scenario 10 is [ 4.0 , 0.0 ]
Optimal solution for scenario 11 is [ 4.0 , 0.0 ]
Optimal solution for scenario 12 is [ 4.0 , 0.0 ]
Optimal solution for scenario 13 is [ 2.0 , 2.0 ]
Optimal solution for scenario 14 is [ 2.0 , 2.0 ]
Optimal solution for scenario 15 is [ 2.0 , 2.0 ]
Optimal solution for scenario 16 is [ 2.0 , 2.0 ]
Optimal solution for scenario 17 is [ 2.0 , 2.0 ]
Optimal solution for scenario 18 is [ 4.0 , 0.0 ]
Optimal solution for scenario 19 is [ 2.0 , 2.0 ]
Optimal solution for scenario 20 is [ 4.0 , 0.0 ]