Cost matrixes are discussed in the book: "Assignment Problems" by Rainer Burkard, Mauro Dell'Amico, Silvano Martello (on pages 73, and 200, etc.). Yes, some parts would be unconnected (most likely) in such a large problem.
Seperation and calculation of the cost matrixes is only the start of the problem. Solving using traditional methods, for so many hubs and nodes, using CPLEX (for example) is going to run into out-of-memory problems.
Solution of these problems is complex (nearly NP-hard at that size; but you'll run into time and memory constraints, as opposed to simply not being able to formulate a solution). I can recommend a couple of papers (with 100's of equations) that address the hub-location-allocation problem (HLAP).
"Bi-objective optimization of multi-server intermodal hub-location-allocation problem in congested systems: modeling and solution" by Mahdi Rashidi Kahag, Seyed Taghi Akhavan Niaki, Mehdi Seifbarghy and Sina Zabihi, from the Journal of Industrial Engineering International (June 2019, Volume 15, Issue 2, pp 221–248):
A new multi-objective intermodal hub-location-allocation problem is modeled in this paper in which both the origin and the destination hub facilities are modeled as an M/M/m queuing system. The problem is being formulated as a constrained bi-objective optimization model to minimize the total costs as well as minimizing the total system time. A small-size problem is solved on the GAMS software to validate the accuracy of the proposed model. As the problem becomes strictly NP-hard, an MOIWO algorithm with an efficient chromosome structure and a fuzzy dominance method is proposed to solve large-scale problems.
In this paper, a new bi-objective intermodal hub-location-allocation problem within queuing framework was investigated. The developed model of this problem intended to simultaneously minimize the total network costs and the total time spent in the network. The solution of a small-size problem using the GAMS software was presented to show the accuracy of the developed formulation. As the problem belongs to the class of NP-Hard problems with conflicting objectives, an MOIWO algorithm with a new chromosome structure was utilized to solve the problem. Moreover, a fuzzy ranking method was used in MOIWO to prioritize the solutions obtained in each iteration.".
I kept that quote short as I prefer the second paper.
In the paper: "Managing facility disruption in hub-and-spoke networks: formulations and efficient solution methods" by Nader Azizi, an open access publication of the: Annals of Operations Research (January 2019, Volume 272, Issue 1–2, pp 159–185), in the first section, he offers an introduction to his paper: "The focus of this paper is on the Uncapacitated Single Allocation p-Hub Location problem where the number of hub facilities to open is predetermined. This problem is also known as the Single Allocation p-Hub Median problem.".
He details some of the history of solving this problem:
"Early research on hub-and-spoke systems focused on developing efficient mathematical formulations and solution techniques for this paradigm. O’Kelly (1986, 1987) presented the first mathematical model for the single allocation p-hub median problem. Another pioneer of the hub location research, Campbell (1994), developed a linear integer formulation for the problem. Based on the idea of multicommodity flow, Ernst and Krishnamoorthy (1996) proposed a new set of formulations for both single and multiple allocation cases. A widely used formulation in the literature for single and multiple allocation p-hub median problem is the work of Skorin-Kapov et al. (1996). The MIP models proposed by Skorin-Kapov et al. (1996) yield to tight linear relaxation.
Over the years, a number of solution techniques e.g., approximation and exact methods have been developed and successfully applied to solve instances of the hub location-allocation problem. Examples of such methods include Simulated Annealing (Ernst and Krishnamoorthy 1999), Genetic Algorithm (Kratica et al. 2007), Hybrid GA and Tabu Search (Abdinnour-Helm 1998), Lagrangian Relaxation (Elhedhli and Wu 2010; Contreras et al. 2009), Benders Decomposition (Contreras et al. 2012), Branch and Bound (Ernst and Krishnamoorthy 1998), and Branch and Cut (Yaman and Carello 2005) among others.
More recent studies focus on the extension of the classical version and aim to develop more realistic hub location-allocation problems such as models with congestion (e.g., Elhedhli and Wu 2010; Azizi et al. 2017), flow dependent economies of scale (Camargo et al. 2009; Campbell et al. 2005), and competitive hub location models (Eiselt and Marianov 2009). For an overview of hub location problem the reader may refer to Campbell and O’Kelly (2012) and Contreras (2015).".
"In short, the contribution of the paper include:
(i) An investigation into a reliable hub-and-spoke system incorporating heterogeneous probability of hub failure
(ii) A non-linear model purposely designed to simultaneously construct networks with and without symmetrical flow and select backup facilities for every demand possibly affected by a disruption
(iii) An improved linear model with the use of indicator constraints and
(iv) Efficient Particle Swarm Optimization (PSO) based meta-heuristics for large instances.
The remainder of this paper is organised as follow. In Sect. 2, we describe the problem and present non-linear, linear and improved linear formulations. The details of the three proposed PSO based algorithms are presented in Sect. 3 followed by computational results and discussions in Sect. 5. Concluding remarks are presented in Sect. 6.".
In section 2 he explains (in detail, with almost 100 equations) that a Reliable Single Allocation p-Hub Location problem (RSApHL-I) involves a Mixed Integer Quadratic Programing (MIQP) formulation, and includes $n^2+4n^3+n^4$ variables. His approach is similar to the technique proposed by Chaovalitwongse et al. (2004) to linearize multi-quadratic 0–1 programming problems. He demonstrates that the problems could be solved more effectively and with significantly less computational efforts when all Big-M formulations are replaced with indicator constraints.
He concludes section 2 by saying:
"Due to the limitation in solving the above model formulation to optimality for realistically sized problem instances with CPLEX, one way forward is to design an efficient metaheuristic. In this study, we developed three metaheuristics based on the well-known evolutionary algorithm of Particle Swarm Optimisation (PSO). These algorithms are discussed in section 3.".
"A large portion of the published research works on PSO is dedicated to continuous optimization problems whereas there is a lack of research dealing with its discrete counterpart namely combinatorial optimization problems (Liu et al. 2008). To the best of our knowledge, our proposed algorithms is the first PSO-based optimisation technique to be developed for solving this interesting though challenging class of hub location problem.
In this paper two variations of a basic PSO and a hybrid PSO algorithm that incorporate crossover and memory are proposed. The two versions which we refer to as PSO-v1 and PSO-v2, differ in terms of solutions representation and the mechanism by which a continuous particle position is transformed into one or more discrete solution(s). The hybrid PSO algorithm on the other hand, is an extension of PSO-v2 which incorporates a short-term memory and a crossover operator, both of which are in our knowledge, new ingredients that are successfully embedded for the first time into the PSO search methodology.".
In section 6 the paper concludes:
"In this study, a mixed integer quadratic programming (MIQP) formulation is proposed for the reliable single allocation hub location problem with heterogeneous probability of hub failure. The model assigned a backup facility to every demand point in the network for a fast and a low-cost recovery after hub failure. The objective function of the model minimises the weighted transportation cost in regular situation and the expected cost resulted from hub failures.
The proposed formulation is developed in such way to maintain the size of the problem after linearization. To improve the computational efficiency of the resulting MILP model, all constraints initially formulated using Big Ms are substituted with logical constraints known as Indicator Constraints. Using the proposed formulation, we solved all small instances and some medium size problem instances to optimality and provided lower and/or upper bounds for other problems. To tackle large instances, three novel algorithms based on PSO namely PSO-v1, PSO-v2 and H-PSO are proposed. In PSO-v1 and PSO-v2, we investigated the effect of two different solution representations and the mechanisms through which a continuous particle is transferred to a discrete solution. It was found that PSO-v2 that maps a continuous particle to a single discrete solution performs better than PSO-v1. The performance of the PSO-v2 is then further improved by introducing two new interesting features namely a crossover operator and a memory to keep track of good solutions that are explored during the search. To our knowledge, this is the first time the hybridisation of PSO with those attributes is examined.".
O'Kelly's paper: "Network Hub Location" (2009) is a quick read with a basic overview of the subject.
Campbell and O’Kelly provide an excellent survey on the literature of hub location research in their paper: "Twenty-Five Years of Hub Location Research", published in Transportation Science 46(2):153-169 (January 2013).
Wikipedia: Facility Location Problem, Distributional cost-effectiveness analysis, Allocative efficiency, and Geometric median.
LocalSolver 9.0 documentation: P-median (Code examples for LSP, Python, C++, C#, and Java)
IBM Forum - OPL using CPLEX Optimizer