#!/usr/bin/env python3# Copyright 2025, Gurobi Optimization, LLC# Facility location: a company currently ships its product from 5 plants# to 4 warehouses. It is considering closing some plants to reduce# costs. What plant(s) should the company close, in order to minimize# transportation and fixed costs?## Note that this example uses lists instead of dictionaries.  Since# it does not work with sparse data, lists are a reasonable option.## Based on an example from Frontline Systems:#   http://www.solver.com/disfacility.htm# Used with permission.importgurobipyasgpfromgurobipyimportGRB# Warehouse demand in thousands of unitsdemand=[15,18,14,20]# Plant capacity in thousands of unitscapacity=[20,22,17,19,18]# Fixed costs for each plantfixedCosts=[12000,15000,17000,13000,16000]# Transportation costs per thousand unitstransCosts=[[4000,2000,3000,2500,4500],[2500,2600,3400,3000,4000],[1200,1800,2600,4100,3000],[2200,2600,3100,3700,3200],]# Range of plants and warehousesplants=range(len(capacity))warehouses=range(len(demand))# Modelm=gp.Model("facility")# Plant open decision variables: open[p] == 1 if plant p is open.open=m.addVars(plants,vtype=GRB.BINARY,obj=fixedCosts,name="open")# Transportation decision variables: transport[w,p] captures the# optimal quantity to transport to warehouse w from plant ptransport=m.addVars(warehouses,plants,obj=transCosts,name="trans")# You could use Python looping constructs and m.addVar() to create# these decision variables instead.  The following would be equivalent# to the preceding two statements...## open = []# for p in plants:#     open.append(m.addVar(vtype=GRB.BINARY,#                          obj=fixedCosts[p],#                          name="open[%d]" % p))## transport = []# for w in warehouses:#     transport.append([])#     for p in plants:#         transport[w].append(m.addVar(obj=transCosts[w][p],#                                      name="trans[%d,%d]" % (w, p)))# The objective is to minimize the total fixed and variable costsm.ModelSense=GRB.MINIMIZE# Production constraints# Note that the right-hand limit sets the production to zero if the plant# is closedm.addConstrs((transport.sum("*",p)<=capacity[p]*open[p]forpinplants),"Capacity")# Using Python looping constructs, the preceding would be...## for p in plants:#     m.addConstr(sum(transport[w][p] for w in warehouses)#                 <= capacity[p] * open[p], "Capacity[%d]" % p)# Demand constraintsm.addConstrs((transport.sum(w)==demand[w]forwinwarehouses),"Demand")# ... and the preceding would be ...# for w in warehouses:#     m.addConstr(sum(transport[w][p] for p in plants) == demand[w],#                 "Demand[%d]" % w)# Save modelm.write("facilityPY.lp")# Guess at the starting point: close the plant with the highest fixed costs;# open all others# First open all plantsforpinplants:open[p].Start=1.0# Now close the plant with the highest fixed costprint("Initial guess:")maxFixed=max(fixedCosts)forpinplants:iffixedCosts[p]==maxFixed:open[p].Start=0.0print(f"Closing plant{p}")breakprint("")# Use barrier to solve root relaxationm.Params.Method=2# Solvem.optimize()# Print solutionprint(f"\nTOTAL COSTS:{m.ObjVal:g}")print("SOLUTION:")forpinplants:ifopen[p].X>0.99:print(f"Plant{p} open")forwinwarehouses:iftransport[w,p].X>0:print(f"  Transport{transport[w,p].X:g} units to warehouse{w}")else:print(f"Plant{p} closed!")