Part 1: Overview and Mathematical Formulation¶
A civil engineering company wants to decide on the projects that they should do. Their objective is to minimize the environmental impact of their projects while making enough profit to keep the company running.
They have a portfolio of 6 possible projects to invest in, where A, B , and C are new infrastructure projects (so-called type 1), and D, E, F are refurbishment projects (so-called type 2).
The environmental impact of each project is given by $I_i$ where $i \in [1,(...),6]$ is the index of the project. $I_i=[90,45,78,123,48,60]$
The profit of each project is given by $P_i$ where $i\in [1,(...),6]$ is the index of the project: $P_i=[120,65,99,110,33,99]$
The company is operating with the following constraints, please formulate the mathematical program that allows solving the problem:
- The company wants to do 3 out of the 6 projects
- the projects of type 2 must be at least as many as the ones of type 1
- the profit of all projects together must be greater or equal than $250$ ($\beta$)
Task 1: Writing the mathematical formulation
Write down every formulation and constraint that is relevant to solve this optimization problem.
Solution $$Min(Z) \sum_{i=1}^{6} x_i I_i $$
subject to
$$ \sum_{i=1}^{6} x_i = 3 $$$$ \sum_{i=4}^{6} x_i \ge \sum_{i=1}^{3} x_i $$$$ \sum_{i=1}^{6} x_i P_i \ge β $$$$x\in(1,0), i\in(1,....,6)$$Extra Challenge (Task 6) $$ \sum_{i=1}^{6} x_i I_i \le x_1 γ + (1-x_1) M $$
where M is a sufficiently big number such as 9999
Task 2: Setting up the software
We'll continue using Gurobi this week, which you've set up in last week's PA. We'll use some other special packages as well. Therefore, be sure to use the special conda environment created for this week.
import gurobipy as gp
Task 3: Setting up the problem
Define any variables you might need to setup your model.
# YOUR_CODE_HERE
# SOLUTION
# Project data
I = [90, 45, 78, 123, 48, 60] # Environmental impact
P = [120, 65, 99, 110, 33, 99] # Profit
# Minimum required profit
beta = 250
M = 100000
# Number of projects and types
num_projects = len(I)
num_type1_projects = 3
num_type2_projects = num_projects - num_type1_projects
Part 2: Create model with Gurobi¶
Remember that examples of using Gurobi to create and optimize a model are provided in the online textbook, and generally consist of the following steps (the first instantiates a class and the rest are executed as methods of the class):
- Define the model (instantiate the class)
- Define variables
- Define objective function
- Add constraints
- Optimize the model
Remember, you can always ask for help to understand a function of gurobi
help(gurobipy.model.addVars)
Task 4: Create the Gurobi model
Create the Gurobi model, set your decision variables, your function and your constrains. Take a look at the book for an example implementation in Python if you don't know where to start.
Solution.
Each line of code does the following:
- Create a Gurobi model
- Add binary variables
- Objective function: Minimize environmental impact
- Constraint: Select exactly 3 projects
- Constraint: Number of type 2 projects must be at least as many as type 1 projects selected
- Constraint: Minimum profit requirement
- Optimize the model
# YOUR_CODE_HERE
# SOLUTION
model = gp.Model("Project_Selection")
x = model.addVars(num_projects, vtype=gp.GRB.BINARY, name="x")
model.setObjective(sum(I[i] * x[i] for i in range(num_projects)),
gp.GRB.MINIMIZE)
model.addConstr(x.sum() == 3, "Select_Projects")
model.addConstr((sum(x[i] for i in range(num_type2_projects, num_projects))
- sum(x[i] for i in range(num_type1_projects)) >= 0),
"Type_Constraint")
model.addConstr(sum(P[i] * x[i] for i in range(num_projects)) >= beta,
"Minimum_Profit")
model.optimize()
Set parameter Username
Set parameter LicenseID to value 2588551
Academic license - for non-commercial use only - expires 2025-11-21
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (win64 - Windows 11.0 (22631.2))
Set parameter LicenseID to value 2588551
Academic license - for non-commercial use only - expires 2025-11-21
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (win64 - Windows 11.0 (22631.2))
CPU model: 13th Gen Intel(R) Core(TM) i7-1365U, instruction set [SSE2|AVX|AVX2]
Thread count: 10 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 3 rows, 6 columns and 18 nonzeros
Model fingerprint: 0xa9bc0e77
Variable types: 0 continuous, 6 integer (6 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+02]
Objective range [5e+01, 1e+02]
Bounds range [1e+00, 1e+00]
RHS range [3e+00, 3e+02]
Presolve time: 0.00s
Presolved: 3 rows, 6 columns, 17 nonzeros
Variable types: 0 continuous, 6 integer (6 binary)
Found heuristic solution: objective 228.0000000
Found heuristic solution: objective 198.0000000
Root relaxation: objective 1.855909e+02, 3 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 cutoff 0 198.00000 198.00000 0.00% - 0s
Explored 1 nodes (3 simplex iterations) in 0.03 seconds (0.00 work units)
Thread count was 12 (of 12 available processors)
Solution count 2: 198 228
Optimal solution found (tolerance 1.00e-04)
Best objective 1.980000000000e+02, best bound 1.980000000000e+02, gap 0.0000%
Task 5: Display your results
Display the model in a good way to interpret and print the solution of the optimization.
# YOUR_CODE_HERE
# SOLUTION
print("Model structure:")
# see the model that you have built in a nice why to interpret
model.display()
if model.status == gp.GRB.OPTIMAL:
print("Optimal Solution:")
for i in range(num_projects):
if x[i].x > 0.9:
print(f"Project {i+1}: Selected")
else:
print("No optimal solution found.")
print("Optimal Objective function Value", model.objVal)
Model structure: Minimize 90.0 x[0] + 45.0 x[1] + 78.0 x[2] + 123.0 x[3] + 48.0 x[4] + 60.0 x[5] Subject To Select_Projects: x[0] + x[1] + x[2] + x[3] + x[4] + x[5] = 3 Type_Constraint: -1.0 x[0] + -1.0 x[1] + -1.0 x[2] + x[3] + x[4] + x[5] >= 0 Minimum_Profit: 120.0 x[0] + 65.0 x[1] + 99.0 x[2] + 110.0 x[3] + 33.0 x[4] + 99.0 x[5] >= 250 Binaries ['x[0]', 'x[1]', 'x[2]', 'x[3]', 'x[4]', 'x[5]'] Optimal Solution: Project 1: Selected Project 5: Selected Project 6: Selected Optimal Objective function Value 198.0
C:\Users\jdding\AppData\Local\Temp\ipykernel_19904\2324086517.py:6: DeprecationWarning: Model.display() is deprecated model.display()
Task 6: Additional constraint
Solve the model with an additional constraint: if project 1 is done then the impact of all projects together should be lower than $\gamma$ with $\gamma=130$.
Paste your model previous model, and call it model2 to keep the results separated and add the new constraint. Then run your second model.
# YOUR_CODE_HERE
# SOLUTION
model2 = gp.Model("Project_Selection")
x = model2.addVars(num_projects, vtype=gp.GRB.BINARY, name="x")
model2.setObjective(sum(I[i] * x[i] for i in range(num_projects)),
gp.GRB.MINIMIZE)
model2.addConstr(x.sum() == 3, "Select_Projects")
model2.addConstr((sum(x[i] for i in range(num_type2_projects, num_projects))
- sum(x[i] for i in range(num_type1_projects)) >= 0),
"Type_Constraint")
model2.addConstr(sum(P[i] * x[i] for i in range(num_projects)) >= beta,
"Minimum_Profit")
# added constraint
gamma = 130
model2.addConstr((sum(I[i] * x[i] for i in range(num_projects))
<= gamma * x[0]+ M * (1 - x[0])),
"Impact_Constraint")
<gurobi.Constr *Awaiting Model Update*>
# YOUR_CODE_HERE
# SOLUTION
model2.optimize()
print("Model structure:")
# see the model that you have built in a nice why to interpret
model2.display()
# Display the solution
if model2.status == gp.GRB.OPTIMAL:
print("Optimal Solution:")
for i in range(num_projects):
if x[i].x > 0.9:
print(f"Project {i+1}: Selected")
else:
print("No optimal solution found.")
print("Optimal Objective function Value", model2.objVal)
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (win64 - Windows 11.0 (22631.2)) CPU model: 13th Gen Intel(R) Core(TM) i7-1365U, instruction set [SSE2|AVX|AVX2] Thread count: 10 physical cores, 12 logical processors, using up to 12 threads Optimize a model with 4 rows, 6 columns and 24 nonzeros Model fingerprint: 0x37d09666 Variable types: 0 continuous, 6 integer (6 binary) Coefficient statistics: Matrix range [1e+00, 1e+05] Objective range [5e+01, 1e+02] Bounds range [1e+00, 1e+00] RHS range [3e+00, 1e+05] Presolve removed 4 rows and 6 columns Presolve time: 0.00s Presolve: All rows and columns removed Explored 0 nodes (0 simplex iterations) in 0.02 seconds (0.00 work units) Thread count was 1 (of 12 available processors) Solution count 1: 228 Optimal solution found (tolerance 1.00e-04) Best objective 2.280000000000e+02, best bound 2.280000000000e+02, gap 0.0000% Model structure: Minimize 90.0 x[0] + 45.0 x[1] + 78.0 x[2] + 123.0 x[3] + 48.0 x[4] + 60.0 x[5] Subject To Select_Projects: x[0] + x[1] + x[2] + x[3] + x[4] + x[5] = 3 Type_Constraint: -1.0 x[0] + -1.0 x[1] + -1.0 x[2] + x[3] + x[4] + x[5] >= 0 Minimum_Profit: 120.0 x[0] + 65.0 x[1] + 99.0 x[2] + 110.0 x[3] + 33.0 x[4] + 99.0 x[5] >= 250 Impact_Constraint: 99960.0 x[0] + 45.0 x[1] + 78.0 x[2] + 123.0 x[3] + 48.0 x[4] + 60.0 x[5] <= 100000 Binaries ['x[0]', 'x[1]', 'x[2]', 'x[3]', 'x[4]', 'x[5]'] Optimal Solution: Project 2: Selected Project 4: Selected Project 6: Selected Optimal Objective function Value 228.0
C:\Users\jdding\AppData\Local\Temp\ipykernel_19904\2231829295.py:8: DeprecationWarning: Model.display() is deprecated model2.display()
