Transportation problem & branch-and-bound method

I have to solve the transportation problem with these limitations (A is the Supply row, B is the Demand column, matrix has transportation prices):

I have found the solution and the code works correctly, because I solved the task and the answers are the same. In brackets is the volume for transportation. For example, in element (1, 2) we have 10 tonnes to transport with the price $8 per one tonne = $80 in total.

However, there is a limitation in my task. I need to use 6-tonne truck for transportation. For example, in found solution in element (1, 2) I have to deliver 10 tonnes with the price $8 per tonne. It means $80 in total, but with my limitation it means that I need to book 2 trucks and to pay $120.

So, taking into account this limitation and using branch-and-bound method I need to find optimal solution. How can I do it?

import gurobipy as gp
from gurobipy import GRB

# Define data
supply = [17, 8, 10, 9]  # Supply nodes
demand = [6, 15, 7, 8, 8]  # Demand nodes
cost = [
    [10, 8, 5, 9, 16],
    [4, 3, 4, 11, 12],
    [5, 10, 29, 7, 6],
    [9, 2, 4, 1, 3]
]  # Cost of transportation from supply i to demand j

m = gp.Model()

# Define decision variables
flow = m.addVars(len(supply), len(demand), lb=0, vtype=GRB.INTEGER, name="flow")

# Define supply constraints
for i in range(len(supply)):
    m.addConstr(gp.quicksum(flow[i, j] for j in range(len(demand))) <= supply[i], name=f"supply_{i}")

# Define demand constraints
for j in range(len(demand)):
    m.addConstr(gp.quicksum(flow[i, j] for i in range(len(supply))) >= demand[j], name=f"demand_{j}")

# Define objective function
m.setObjective(gp.quicksum(flow[i, j] * cost[i][j] for i in range(len(supply)) for j in range(len(demand))), sense=GRB.MINIMIZE)

# Solve model
m.optimize()

# Print results
if m.status == GRB.OPTIMAL:
    print(f"Optimal solution found with objective value {m.objVal:.2f}")
    for i in range(len(supply)):
        for j in range(len(demand)):
            if flow[i, j].x > 0:
                print(f"Flow from supply node {i+1} to demand node {j+1}: {flow[i, j].x:.0f}")
else:
    print("No solution found or problem is infeasible.")

Assign the truck count to an integral variable and tie the price to the truck count. When the solver is given this, it is able to find an optimal solution that does not use more than one truck per route.

Also, worry less about the particular branch-and-bound algorithm; just pick the solver category (MILP) and let it do its work.

import pandas as pd
import pulp

truck_capacity = 6
suppliers = pd.RangeIndex(name='supplier', stop=4)
consumers = pd.RangeIndex(name='consumer', stop=5)

supply = pd.Series(
    name='supply', index=suppliers, data=(17, 8, 10, 9),
)
demand = pd.Series(
    name='demand', index=consumers, data=(6, 15, 7, 8, 8),
)
price_per_tonne = pd.DataFrame(
    index=suppliers, columns=consumers,
    data=(
        (10,  8,  5,  9, 16),
        ( 4,  3,  4, 11, 12),
        ( 5, 10, 29,  7,  6),
        ( 9,  2,  4,  1,  3),
    ),
).stack()
price_per_tonne.name = 'price'

flow = pd.DataFrame(
    index=suppliers, columns=consumers,
    data=pulp.LpVariable.matrix(
        name='flow_s%d_c%d', cat=pulp.LpContinuous, lowBound=0,
        indices=(suppliers, consumers),
    ),
).stack()
flow.name = 'flow'

trucks = pd.DataFrame(
    index=suppliers, columns=consumers,
    data=pulp.LpVariable.matrix(
        name='trucks_s%d_c%d', cat=pulp.LpInteger, lowBound=0,
        indices=(suppliers, consumers),
    )
).stack()
trucks.name = 'trucks'

price = truck_capacity * pulp.lpDot(price_per_tonne, trucks)
prob = pulp.LpProblem(name='transportation', sense=pulp.LpMinimize)
prob.setObjective(price)

# The flow must not exceed the supply
for supplier, group in flow.groupby('supplier'):
    prob.addConstraint(
        name=f'flow_supply_s{supplier}',
        constraint=pulp.lpSum(group) <= supply[supplier],
    )

# The flow must exactly meet the demand
for consumer, group in flow.groupby('consumer'):
    prob.addConstraint(
        name=f'flow_demand_c{consumer}',
        constraint=pulp.lpSum(group) == demand[consumer],
    )

# The capacity must be able to carry the flow
for (supplier, consumer), truck_flow in flow.items():
    prob.addConstraint(
        name=f'capacity_s{supplier}_c{consumer}',
        constraint=truck_flow <= trucks[(supplier, consumer)] * truck_capacity
    )

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal

print('Flow:')
flow = flow.apply(pulp.value).unstack(level='consumer')
print(flow)
print()

print('Trucks:')
trucks = trucks.apply(pulp.value).unstack(level='consumer')
print(trucks)
print()

transportation:
MINIMIZE
60*trucks_s0_c0 + 48*trucks_s0_c1 + 30*trucks_s0_c2 + 54*trucks_s0_c3 + 96*trucks_s0_c4 + 24*trucks_s1_c0 + 18*trucks_s1_c1 + 24*trucks_s1_c2 + 66*trucks_s1_c3 + 72*trucks_s1_c4 + 30*trucks_s2_c0 + 60*trucks_s2_c1 + 174*trucks_s2_c2 + 42*trucks_s2_c3 + 36*trucks_s2_c4 + 54*trucks_s3_c0 + 12*trucks_s3_c1 + 24*trucks_s3_c2 + 6*trucks_s3_c3 + 18*trucks_s3_c4 + 0
SUBJECT TO
flow_supply_s0: flow_s0_c0 + flow_s0_c1 + flow_s0_c2 + flow_s0_c3 + flow_s0_c4
 <= 17
...
flow_demand_c0: flow_s0_c0 + flow_s1_c0 + flow_s2_c0 + flow_s3_c0 = 6
...
capacity_s0_c0: flow_s0_c0 - 6 trucks_s0_c0 <= 0
...
Result - Optimal solution found

Objective value:                276.00000000
Enumerated nodes:               0
Total iterations:               33
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.01

Flow:
consumer    0    1    2    3    4
supplier                         
0         0.0  6.0  5.0  6.0  0.0
1         0.0  6.0  2.0  0.0  0.0
2         6.0  0.0  0.0  0.0  4.0
3         0.0  3.0  0.0  2.0  4.0

Trucks:
consumer    0    1    2    3    4
supplier                         
0         0.0  1.0  1.0  1.0  0.0
1         0.0  1.0  1.0  0.0  0.0
2         1.0  0.0  0.0  0.0  1.0
3         0.0  1.0  0.0  1.0  1.0