## 2 Large-Scale MDO by Using PODA and MMRT
### Analysis Graph of Full-Vehicle Stability Properties and Variation Patterns
At the working point, \(q=q^{k}\) and \(\dot{q}=\dot{q}^{k}\), the first order difference formula of Equation (10) is given as follows  
\[\Bigg{[}\frac{\partial\mathbf{F}}{\partial q}\bigg{|}_{q^{k}\dot{q}^{k}}+\frac{1} {h\beta_{0}}\frac{\partial\mathbf{F}}{\partial\dot{q}}\bigg{|}_{q^{k}\dot{q}^{k}} \Bigg{]}\Delta q=-\mathbf{F}(q,\dot{q},t) \tag{11}\]  
where \(\beta_{0}\) is the scalar constant, that is related to the integrator order; \(\Delta q\) is the corrected difference direction; \(h\) is the integration step size; \(-\mathbf{F}\) is the residual part of Equation (10), which represents the imbalance degree of the force system. The left matrix of Equation (11) is the Jacobean matrix on the required quasi-equilibrium state, i.e., the Jacobean matrix of first-index DAEs for Equation (1)  
\[\mathbf{\Gamma}= \begin{bmatrix}\frac{\mathbf{M}}{h\beta_{0}}-\mathbf{C}_{u}&\mathbf{M}_{q} \dot{u}+\mathbf{C}_{qq}^{T}\lambda-\mathbf{Q}_{q}&\mathbf{C}_{q}^{T}\\ \mathbf{I}&-\frac{1}{h\beta_{0}}&\mathbf{0}\\ \mathbf{0}&\mathbf{C}_{q}&\mathbf{0}\end{bmatrix}  \tag{12}\]  
where \(\mathbf{I}\) is the unit matrix; \(u,q,qq\) are subscripts, i.e., the first or second-order partial derivatives of (non-) holonomic constraints \(\mathbf{C}\) to \(u/q\). Accordingly, the MBS eigenvalue solution \(Eig(\mathbf{\Gamma})\) can be obtained by orthogonal decomposition. Some analysis software can provide systemic modal analysis tools, i.e., root-locus graphs of closed-loop nonlinear systems. The disciplinary surrogate models used herein to approach the (near-) linear relationship of wheel-rail contact, i.e., the geometric linearity between wheel-rail contact, are represented by mono-curvature contact between wheel and rail arcs. The HSR practices must grasp the dialectical relationship between (near-) linearity and nonlinearity of wheel-rail contact. The (near-) linear surrogate model of wheel-rail contact is composed of such parameters as \(\lambda_{e}\), which is the equivalent parameter given by the RRD (Rolling Radius Difference) curve of the wheelset, similar to the AC voltage of 220/380 V. The nominal equivalent conicity \(\lambda_{eN}\) refers to the equivalent conicity with a wheelset hunting displacement amplitude of 3 mm under the specific wheel-rail profile matching. For rail 60E1/wheel S1002, with the back-to-back distance of 1360 mm and bottom inclination of 1/20, \(\lambda_{eN}\approx 0.01\) and the track window becomes narrower. While considering rail CN60 KG/wheel XP55 with the back-to-back distance of 1353 mm and bottom inclination of 1/40, \(\lambda_{eN}\approx 0.06\) and track window becomes wider. For the worn wheels and rails, the actual equivalent conicity \(\lambda\)e will be affected by two aspects due to the variation of the RRD curve: firstly, when the negative slope at the zero crossing point causes the local conformal contact at one side of rail gauge corner, which in turn transforms into frictional heat loss and eventually evolves into the local RCF (Rolling Contact Fatigue) failure or plastic flow; secondly, when the discontinuity at the zero crossing point occurs, the local conformal contact is constituted on the top of the raillhead, which will occasionally force the wheelset to form the small-amplitude hunting oscillation and gradually constitute the sufficiently strong external excitation source. Compared with the conventional root locus graph, the analysis graph of full-vehicle stability properties and variation patterns has three important features: i.e., closed-loop pole, stability margin (or critical damping), and convected motion relationship. The convected motion relationship refers to the convected motion formed between relevant modes caused by the complex constraint singularity in the generalized space, which constitutes the corresponding kinetic energy exchange or transformation. Generally, the nonlinear systems no longer meet the following requirements, i.e., modal orthogonality and mode shape symmetry or anti-symmetry.