Enabling Hybrid Modeling in Commercial MBS Software: A Force-Level Approach

1. Introduction

Accurate modeling of dynamic systems is a central task in engineering. Traditional approaches are typically classified as white-box, black-box, or gray-box modeling, each with inherent advantages and limitations. [1,2]

White-box models are based entirely on physical principles, such as Newton’s laws or balance equations, and require detailed knowledge of system parameters and interactions. They offer transparency, interpretability, and extrapolation capabilities, but often suffer from modeling inaccuracies due to simplifications or unknown effects, such as friction or wear. [1,2]

Black-box models, by contrast, rely purely on data. Machine learning techniques are used to learn input–output relationships without explicit knowledge of the underlying physics. These models are flexible and can capture complex phenomena, but they require large datasets, lack physical interpretability, and often perform poorly outside the range of the training data. [1,2]

Gray-box models blend both approaches by embedding data-driven components into physical structures, or vice versa. This hybridization improves generalization while capturing unmodeled dynamics. Examples include physics-informed neural networks (PINNs), which incorporate physical laws into the training loss, and surrogate models that approximate detailed simulations to reduce computational costs and enable an efficient analysis of macroscopic behavior. [1,2]

Hybrid modeling, a specific form of gray-box modeling, emphasizes a close interplay between physical and data-driven components. Neither model dominates; instead, both work in tandem to bridge the gap between idealized models and real-world behavior. [1,3]

Commercial MBS software such as MSC Adams or Simpack already offer interfaces for co-simulation, e.g. with MATLAB/Simulink. These are often used for control engineering or complex subsystems, but in principle also allow the integration of data-driven models for model correction. It has been shown that hybrid approaches can be successfully used for modeling specific components such as rubber-metal bushings [4]. This can be easily implemented as a discrepancy of force using a co-simulation. In this article, however, the focus is on cases in which modeling differences cannot be clearly assigned to individual elements, but result from general model simplifications or unknown effects. The aim is to systematically compensate for such discrepancies at the force level and thus to be able to model complex dynamic systems more realistically.

2. Hybrid Modeling

For a dynamic system, measurements $(t_k,\overline{\underline{x}}\left(t_k\right))$ are given. Furthermore, a white-box model in form of an equation of motion

$$\dot{\tilde{\underline{x}}}=\underline{f}(t,\underline{x})$$

(1)

is given, which does not fully capture the actual behavior of the system and therefore only provides an approximation $\tilde{\underline{x}}$. This may be due to unmodeled phenomena or imprecise parameterization. As a result, a discrepancy

$$\Delta \dot{\underline{x}}=\dot{\overline{\underline{x}}}-\dot{\tilde{\underline{x}}}$$

(2)

between measurement $\dot{\overline{\underline{x}}}$ and model $\dot{\tilde{\underline{x}}}$ arises.

In [3] a continuous hybrid modeling approach was proposed, in which the discrepancy (2) is modeled using a black-box method. This correction term is then added to the white-box model, resulting in the hybrid formulation

$${\dot{\underline{x}}}^{*}=\underline{f}(t,\underline{x})+\Delta \dot{\underline{x}}(t,\underline{x}).$$

(3)

For multibody systems (MBS), the equation of motion is generally given by the second-order differential equation

$$\underline{M}(\underline{q},t)\ddot{\underline{q}}+{\underline{k}}^{\mathrm{c}}(\underline{q},\dot{q},t)={\underline{k}}^{\mathrm{e}}(\underline{q},\dot{\underline{q}},t),$$

(4)

where $\underline{M}$ is the mass matrix, ${\underline{k}}^{\mathrm{c}}$ contains the generalized centrifugal and Coriolis forces, and ${\underline{k}}^{\mathrm{e}}$ contains the generalized applied forces. [5]

Then the introduced hybrid model can be implemented in the form

$$\begin{array}{ccccc}\underbrace{\left[\begin{array}{c}\dot{\underline{q}}\\ \ddot{\underline{q}}\end{array}\right]}& =& \underbrace{\left[\begin{array}{c}\dot{\underline{q}}\\ -{\underline{M}}^{-1}{\underline{k}}^{\mathrm{c}}(\underline{q},\dot{\underline{q}})+{\underline{M}}^{-1}{\underline{k}}^{\mathrm{e}}(\underline{q},\dot{\underline{q}})\end{array}\right]}& +& \underbrace{\left[\begin{array}{c}\underline{0}\\ \Delta \ddot{\underline{q}}\end{array}\right]},\\ \dot{\underline{x}}& =& \underline{f}(t,\underline{x})& +& \Delta \dot{\underline{x}}(t,\underline{x}).\end{array}$$

(5)

However, in many practical cases – particularly when using commercial simulation software such as MSC Adams or Simpack – direct modification of the equations of motion is not possible. [5,6]

2.1. Force Variation

To address this limitation, the black-box correction can be applied at the force level instead of at the acceleration level. The governing equation is reformulated as

$$\underline{M}\ddot{\underline{q}}+{\underline{k}}^{\mathrm{c}}(\underline{q},\dot{\underline{q}})+\Delta \underline{F}(\underline{q},\dot{\underline{q}})={\underline{k}}^{\mathrm{e}}(q,\dot{q}),$$

(6)

from which the force discrepancy can be computed as

$$\Delta \underline{F}(\underline{q},\dot{\underline{q}})=\underline{M}(\underline{q},t)\ddot{\underline{q}}+{\underline{k}}^{\mathrm{c}}(\underline{q},\dot{\underline{q}},t)-{\underline{k}}^{\mathrm{e}}(\underline{q},\dot{\underline{q}},t).$$

(7)

This formulation allows hybrid correction to be introduced using a force element, without the need to alter the underlying integration scheme. It is therefore well suited for implementation in commercial simulation tools like MSC Adams (Section 2.3.1) or Simpack (Section 2.3.2).

2.2. Application Example: Single-Mass Oscillator

A single-mass oscillator with external excitation, as shown in Figure 1a), is used as an application example. The mass m is connected to the environment by a non-linear spring and a non-linear damper and is excited with a harmonic force $F\left(t\right)=\widehat{F}sin\left(\omega t\right)$. In addition, a frictional force $F_{\mathsf{\mu }}$ acts between the mass and the ground. The displacement of the mass is denoted by q. The free-body diagram can be seen in Figure 1b). The complete equation of motion is given by

$$m\ddot{q}+F_c\left(q\right)+F_d\left(q\right)+F_{\mathsf{\mu }}\left(\dot{q}\right)=F\left(t\right),$$

(8)

with

$$\begin{array}{cc}\hfill F_c\left(q\right)& =c{q}^3,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill F_d\left(\dot{q}\right)& =d\left|\dot{q}\right|\dot{q},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill F_{\mathsf{\mu }}\left(\dot{q}\right)& =F_{\mathrm{N}}\mathsf{\mu }\mathrm{sgn}\left(\dot{q}\right)=mg\mathsf{\mu }\mathrm{sgn}\left(\dot{q}\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill F\left(t\right)& =\widehat{F}sin\left(2\mathsf{\pi }ft\right).\hfill \end{array}$$

(12)

The white-box model chosen as the basis for hybrid modeling neglects the dissipative effects, so that the reduced equation of motion is given by

$$m\ddot{q}+F_c\left(q\right)=F\left(t\right).$$

(13)

The force discrepancy can then be calculated via

$$\Delta \overline{F}=-m\ddot{\overline{q}}-F_c\left(\overline{q}\right)+F\left(t\right)$$

(14)

and serves as target for the black-box model. The input consists of the displacement q and the velocity $\dot{q}$.

2.3. Implementation

Synthetic measurement data is generated using MATLAB. Therefore, the equation of motion (8) is solved numerically using the Dormand-Prince method from the Runge-Kutta family. Additionally, a measuring noise of $1.5\%$ of the maximum values is added. The results can be seen in Figure 2 in black. The numerical solution of the white-box model (13) is shown in gray.

Equation (14) is used to calculate the discrepancies $\Delta F$, which serve as target. As input, the state variables q and $\dot{q}$ are used. A feedforward, fully connected neural network for regression with the MATLAB default setting is used as the black box model. Randomly, $70\%$ of the data is used for training and the remaining $30\%$ for testing. The root mean square error (RMSE)

$$\mathrm{RMSE}(\Delta \tilde{F},\Delta \overline{F})=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(\Delta \tilde{F}\left(t_i\right)-\Delta \overline{F}\left(t_i\right))}^2}$$

(15)

between the training data $\Delta \overline{F}$ and the prediction $\Delta \tilde{F}$ is $0.90$ for the training data set, and $0.81$ for the test data set.

The results can be seen in Figure 3. The trained model is saved in a *.mat file for later use with MSC Adams and Simpack.

2.3.1. MSC Adams

The white-box model (13) is built in MSC Adams. The elements used and the setup can be seen in Figure 4a)-b). It consists of a point mass and three SForce elements: the nonlinear spring force F_spring, the excitation F, and the discrepancy force F_Delta. Three state variables are defined: the displacement state_q and the velocity state_dq of the point mass, as well as the force Delta_F, which is the value of the force F_Delta. To calculate Delta_F, the black-box model is integrated by a co-simulation with MATLAB/Simulink. For this purpose, the model is exported using the controls plugin. As input, the state variable Delta_F is selected, and the state variables state_q and state_dq are selected as outputs. [7]

In MATLAB/Simulink, the exported S-function model "MSC Software" can be coupled with the black-box model using the "RegressionNeuralNetwork Predict" block. The model structure in Simulink is shown in Figure 4c).

2.3.2. Simpack

As for the application in MSC Adams, the white-box model (13) is now implemented in Simpack 2024x. The elements used and the setup can be seen in Figure 5a)-b) It consists of a point mass and three force elements. The point mass is connected to the initial system via a joint Element, allowing only motion in one direction. The force elements are a non-linear spring force FSpring connecting the point mass to the origin, an excitation force FExcitation that provides the harmonic excitation, and a controllable input Fdelta. As a model output, there are two stat variables defined, the displacement Yq and the velocity Ydq of the point mass. The force variable UIdeltaF is defined as the model input. For the calculation of the discrepancy force deltaF the data-driven model is integrated via a co-simulation with MATLAB/Simulink. To start a co-simulation between Simpack and MATLAB/Simulink, a "simat" block is added to Simulink to generate the interface. In the "simat" settings, the absolute path of the Simpack model is entered. Similar to the hybrid model using MSC Adams, the "simat" block is coupled with the "RegressionNeuralNetwork Predict" block. The model structure in Simulink is shown in Figure 5c). The command server for co-simulation must also be started in Simpack. Then the co simulation of the two programs can be started from MATLAB/Simulink.

3. Results

The results using the different software are shown in Figure 6, a quantitative evaluation can be found in Table 1. The RMSE is used for the displacement q and for the velocity $\dot{q}$.The synthetic measurements are shown in red. For comparison, the results of the continuous method (MATLAB) are also shown in black. A neural network with the default settings was trained as the data-driven model. Since the targets of the black-box models for $\Delta F$ and $\Delta \ddot{q}$ have different units, the normalized root mean square error (NRMSE) [8]

$$\mathrm{NRMSE}(\Delta \tilde{F},\Delta \overline{F})={\displaystyle \frac{\mathrm{RMSE}\left((\Delta \tilde{F},\Delta \overline{F})\right)}{max(\Delta \overline{F})-min(\Delta \overline{F})}}$$

(16)

for $\Delta F$ and analogous for $\Delta \ddot{q}$ is given in Table 1.

The results are close to the measurements, only at the beginning there is a slight phase shift in which the hybrid model lags behind the measurements. This is probably due to the transient response of the system and its underrepresentation in the data.

The results of the hybrid modeling MSC Adams (green) and Simpack (magenta) are almost congruent with each other and with the measurements. The differences can be attributed to the use of different solvers.

4. Discussion

Overall, the results of the hybrid modeling using both the original continuous method $\Delta \dot{x}$ from [3] and the variation $\Delta F$ are very good and also very similar. Since no optimization of the hyperparameters took place, the results presented here do not allow any statement to be made as to which of the methods provides more accurate results.

5. Conclusions

In this work, a variation of the continuous method for hybrid modeling based on force-level discrepancies was presented. This enables the correction of modeling inaccuracies using data-driven elements without modifying the system’s differential equations or interfering with the numerical solver, thus allowing seamless integration into commercial MBS software such as MSC Adams or Simpack.

While previous hybrid modeling efforts have focused on clearly identifiable components - such as rubber-metal bushings, the approach presented here addresses more general model inaccuracies that cannot be assigned to a specific force element. This makes it suitable for complex multibody systems where simplifications, parameter uncertainties, or unknown effects affect the model as a whole.

The study demonstrates the general feasibility of this approach, though it has not yet been applied to real measurements. As system complexity increases, challenges will arise, particularly in localizing and interpreting force discrepancies. In such cases, expert knowledge remains crucial.

Nevertheless, the promising results with synthetic data, combined with compatibility across established simulation tools, provide a solid foundation for future work. In particular, the ability to embed hybrid models component-wise offers new potential for realistic simulation and efficient modeling of complex mechanical systems.