## 5 Solution Approach  
There are several algorithms for solving the SUE problem, such as method of successive averages (MSA) [31], the disaggregated simplicial decomposition method (DSD) [43], and the gradient projection method (GP) [44]. We developed a method based on the MSA algorithm to solve the proposed PSL-SUE problem with train capacity constraints.  
Step 1: Network construction and valid path set.  
The schedule-based network is constructed according to the established timetable of cross-line operation with state and suburban trains. In the case of free flow, the \(k\) shortest paths between each O-D pair are found by using the \(k\) shortest path algorithm and then stored into the valid path set \(K\).  
Step 2: Initialization.  
Initialize all arc flows to 0 and calculate the arcs' impedance under the condition of free flow. The generalized costs of all paths are calculated. Assign the passengers to the paths by using the MNL mode on the zero-flow network, and the initial path flow \(f_{0}\) is obtained and set the number of iteration \(n=0\).  
Step 3: Update arc impedance.  
According to the path flow \(f^{n}\) of the \(n\)th iteration, the new arc flow \(x_{a}^{n}\) is updated by Equation (22). Then, update the impedance of all arcs and generalized cost of all paths based on the generalized travel cost function mentioned above.  
Step 4: Iteration.  
According to the updated generalized cost of paths, the passengers are reassigned to paths by using the PSL model  
\[f_{(i,j)}^{k,n}=q_{(i,j)}\cdot\frac{\exp(-\theta\cdot\mathrm{C}_{(i,j)}^{k,n} +\ln\mathrm{S}_{(i,j)}^{k})}{\sum\limits_{k^{\prime}\in K(i,j)}\exp(-\theta \cdot\mathrm{C}_{(i,j)}^{k^{\prime},n}+\ln\mathrm{S}_{(i,j)}^{k^{\prime}})},k \in K(i,j). \tag{26}\]  
Update the arc flow according to the new path flow, and judge whether the arc flow exceeds the train capacity. If the arc flow exceeds the train capacity, it is necessary to reassign the path flow by reducing it proportionally. If \(x_{a}^{n}>cp,\eta_{a}^{n}<1\), the path flow will be reduced. Otherwise, if \(x_{a}^{n}>cp,\eta_{a}^{n}<1\), the path flow will not change.  
\[cp/x_{a}^{n}=cp/\sum\limits_{(i,j)\in(O,D)}\sum\limits_{k\in K(i,j)}\delta_{( i,j),k}^{q}f_{(i,j)}^{k,n} \tag{27}\]  
\[\eta_{a}^{n}=\min[1,cp/x_{a}^{n}] \tag{28}\]  
\[f_{(i,j)}^{k,n^{\prime}}=f_{(i,j)}^{k,n}\cdot\eta_{a}^{n}. \tag{29}\]  
While the path flow is reduced, the total passenger number of the O-D pair \((i,j)\) will also be reduced, so it is also necessary to enlarge the total number of reduced path flow to the original value  
\[\eta_{(i,j)}^{n}=q_{(i,j)}/\sum\limits_{k\in K(i,j)}f_{(i,j)}^{k,n^{\prime}}. \tag{30}\]  
Thus, at the \(n+1\)th iteration, the new ancillary path flow can be obtained by multiplying the \(n\)th iteration path flow by the reduction factor and amplification factor successively  
\[f_{(i,j)}^{k,n+1}=f_{(i,j)}^{k,n}\cdot\eta_{a}^{n}\cdot\eta_{(i,j)}^{n}. \tag{31}\]  
Step 5: Convergence and output results.  
Stop iterating if the convergence criterion (32) is satisfied. Otherwise, return to Step 3  
\[\sum\limits_{(i,j)\in(O,D)}\sum\limits_{k\in K(i,j)}\frac{\sqrt{\left(f_{(i,j )}^{k,n+1}-f_{(i,j)}^{k,n}\right)^{2}}}{f_{(i,j)}^{k,n}}\leq\epsilon. \tag{32}\]  
The solution approach is shown in Figure 6.