Clustering problem involving multidepots and customers requiring commodities located exclusively in an specific depot

I'm trying to solve a clustering problem that's similiar to a VRP Pickup and Delivery problem with multiple depots and customers.

Each customer demands a commodity that is exclusively found on one depot. I want to generate routes (not even trying to find the optimal order) that fulfill customers' requirements but also minimizing the total cost (distance) of those routes, having a set of Q-capacitated vehicles.

I draw my problem here:

Any help will be deeply appreciated. I'm trying to study how to make these "clusters" happen. Tried VRP Pickup and Delivery problems, Dial-a-Ride problems and so. But I can't find a way to solve this by some heuristic (seems extremely NP-hard for an exact method).

This is a tough problem indeed, but I am not sure about the "extremely NP-hard" part :). All problems which are NP-hard are...very hard.

This looks like a multi-commodity flow problem, one commodity per depot.

It is natural to decompose such a problem as follows. For each customer $$v\in V$$, for each commodity $$k\in K$$ we assume that the demand $$D_{vk}\in \mathbb{R}^+$$ is known. Let $$\Omega_k^v$$ denote all possible paths from depot/commodity $$k\in K$$ to customer $$v \in V$$, and let $$\Omega_k =\bigcup_{v \in V}\Omega_k^v$$ denote the set of all paths from depot $$k$$.

Let $$y_{p}$$ be a binary variable that takes value $$1$$ if and only if path $$p\in \Omega_k^v$$ is used to satisfy the demand of customer $$v\in V$$ for commodity $$k\in K$$. Let $$n_{ij} \in \mathbb{N}$$ be a variable for the number of vehicles of capacity $$Q$$ required on edge $$(i,j)\in A$$, and let $$c_{ij}$$ denote the cost of using one vehicle on edge $$(i,j)$$.

You want to minimize the overall cost :

$$\sum_{(i,j)\in A}c_{ij} n_{ij}$$ subject to:

1. Each customer $$v$$ must be serviced by exactly one path for each commodity: $$\sum_{p\in \Omega_k^v } y_{p} = 1\quad \forall v \in V, \; \forall k\in K$$
2. Capacity constraints must hold: $$\sum_{k \in K}\sum_{v \in V}\sum_{p\in \Omega_k^v|(i,j) \in p} D_{vk}y_{p} \le Q n_{ij}\quad \forall (i,j) \in A$$

Note that there is an exponential number of paths in a network, so you will not be able to define the entire set $$\Omega_k$$. You have two options:

1. Define a subset of such paths, e.g., the $$b$$ best ones (shortest ones) for each $$(k,v)$$. See how the model grows when $$b$$ varies, and how the quality of the solution evolves.
2. Use column generation to define paths dynamically. But since you have $$n_{ij}$$ variables that are integer, you will need to embed this in a branch-and-bound tree. If you are not familiar with branch-and-price, I would go with option $$1$$. In my experience option $$1$$ already gives good results.

Also, note that once you solve this problem, you do not have individual routes per vehicle. You "only" have the routing of each commodity, for each customer, and the number of vehicles that are required on each edge. For individual routes per vehicle, you need to solve another problem after this one : you need to "cover" your network with routes, such that on each edge you have the required number of vehicles. And last but not least, all of this assumes that you may use multiple vehicles to send commodity $$k$$ to customer $$v$$, as you may need to transfer the commodity to another vehicle at one point (in order to minimize $$n_{ij}$$).