Controller Parameters Tuning Based on Transfer Matrix Method for Multibody Systems

Abstract

Transfer matrix method for multibody systems (MS-TMM) is a rife method to multi-rigid-flexible-body systems dynamics model deduction due to that there are no needs to establish the global dynamics equations of the system. Its basic idea is transferring a state vector between the body input(s) and output(s); this idea is close to the linear theories in control analysis and design. In this paper, three controllers’ parameters tuning techniques for the proposed system model using MS-TMM are utilized; one technique is applied to get the stability regions via the frequency response of MS-TMM derived model. Another technique considers a classical PID controller design through the analysis of step input response of the system, and the last technique can be applied in both time and frequency domains if the model has a known mathematical model. A car suspension system is considered to represent modeling and tuning problems. In-depth study of MS-TMM with control techniques and defining the controllers’ parameters stability regions provide an opportunity to formulate a relationship between MS-TMM and control design for novel control applications due to the powerful strength of MS-TMM dealing with more complex problems of the controlled multibody systems.

1. Introduction

With the increments of complexity of multibody systems and the development of their design and control methods, the need for more elegant formulations of the equations of motion becomes an issue of paramount importance. Many methods and theories for developing the model of the multibody system dynamics and control are presented for such reasons. In transfer matrix method for multibody systems (MS-TMM) there are two cases to deal with control element, one is to express the control force with state of system of previous time such as the delay control, and the second is that the control force is relative to present state of system, such as real-time control systems. Rui et al. have trials to develop new controlled systems depending on deriving the dynamics equations using MS-TMM [1, 2]. Bestle et al. reformulated the car suspension system using MS-TMM as similar as classical control theory's structural diagram [3].

Proportional-Integration-Differentiation (PID) controller is prevalent in industry applications. Reference [4] indicated that more than 90% of feedback control loops are based on PID control and [5] indicated that more than 97% of regulatory controllers utilize the PID algorithm. But the tuning of the controller gains is a problem because many industrial models suffer some burdens such as nonlinearities, higher order, and time delay [6]. The basic idea of MS-TMM is transferring the state vector between the body input(s) and output(s); this idea is close to linear theories in control analysis and design. The control and feedback variables could be considered in the transfer equation of MS-TMM. The tuning of the controller parameters is a great issue to get a stabilized system. However Ziegler-Nichols technique is one of the simplest techniques of PID tuning; there are several disadvantages, such as many trials are needed to find the stability limits and also not all tuning values work efficiently for all situations [7–10]. Shamsuzzoha and Skogestad (2010) developed a new procedure for PI/PID tuning method in closed-loop mode, only one closed-loop step test is required to deduce PI controller values. In this method also it is simple to obtain the PID tuning parameters in one step [11, 12]. This paper is organized as follows: the proposed system is modeled by MS-TMM including the control law in Section 2. The design of PI controller values stability region that achieves a specified gain and phase margins, in addition to tuning the controller parameters using different techniques, is presented in Section 3. Simulation and results of car active suspension system model with the designed controllers are provided in Section 4. Conclusions are offered in Section 5.

2. Problem Formulation

A two-degree-of-freedom car dynamic model is established to simulate the system as shown in Figure 1. It does not only simplify the system design and analysis, but also represents most of full car features [13]. In Figure 1, m₁ is the sprung mass that represents the car body, m₃ is the unsprung mass which stands for the wheel, C₂ is the damping coefficient, and K₂ and K₄ are the springs stiffness. In this section, this model is deduced using MS-TMM. According to MS-TMM, the system is portioned into many elements which are numbered as 1 to 4 from up to down, where body 1 and 3 are lumped masses, element 4 is a linear spring, while the controlled element 2 includes spring and damper in parallel connections. The boundaries at two ends are noted as 0. The positive direction of x axis is also shown in the figure. It should be mentioned that the controlled force f _c serves as an internal force of the system.

Figure 1:

Dynamic model of controlled car.

The transfer direction is stipulated as the positive direction of x axis. The state vector of an arbitrary connection point of the system is defined as ${\mathbf{Z}}_{i,j}={[X,Q_x]}_{i,j}^T$ , where i and j stands for the sequence number of the body element and hinge element, respectively. Thereby, the transfer equation of each element is discussed as follows:

Considering the rise and fall of the road, the boundary end x_4,0 is assumed to be in the form of x_4,0 = Ae ^iΩt , where A is the complex amplitude and Ω is the frequency of the movement. Thus, in a steady-state sense, the whole system experiences the forced vibration in the sine form.

According to MS-TMM for the lumped mass 1 and 3, their transfer equations are

${\mathbf{Z}}_{j,O}={\mathbf{U}}_j{\mathbf{Z}}_{j,I},\hspace{1em}\left(j=\mathrm{1,3}\right),$ (1)

where ${\mathbf{U}}_j=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ m_j{\Omega }^{\mathrm{2}}& \mathrm{1}\end{array}\right]$ and Ω is the external excitation frequency.

For the linear spring 4, the transfer equation is

${\mathbf{Z}}_{\mathrm{4},O}={\mathbf{U}}_{\mathrm{4}}{\mathbf{Z}}_{\mathrm{4},I}$ (2)

with transfer matrix ${\mathbf{U}}_{\mathrm{4}}=\left[\begin{array}{cc}\mathrm{1}& -\mathrm{1}/K_{\mathrm{4}}\\ \mathrm{0}& \mathrm{1}\end{array}\right]$ .

For the controlled element 2, the transfer equation should be

${\mathbf{Z}}_{\mathrm{2},O}={\mathbf{U}}_{\mathrm{2}}{\mathbf{Z}}_{\mathrm{2},I}+{\mathbf{E}}_c{F}_c,$ (3)

where ${\mathbf{U}}_{\mathrm{2}}=\left[\begin{array}{cc}\mathrm{1}& -(\mathrm{1}/(K_{\mathrm{2}}+i{C}_{\mathrm{2}}\Omega ))\\ \mathrm{0}& \mathrm{1}\end{array}\right]$ , ${\mathbf{E}}_c=\left[\begin{array}{c}1/(K_{\mathrm{2}}+i{C}_{\mathrm{2}}\Omega )\\ \mathrm{0}\end{array}\right]$ , and F _c is the control force.

Considering the negative PID feedback control and regarding x_1,0 as the measurement signal, the control force takes the form

$f_c=-\left(k_P{x}_{\mathrm{1,0}}+k_I{\int }_{\mathrm{0}}^t{x}_{\mathrm{1,0}}dt+k_D\frac{d{x}_{\mathrm{1,0}}}{dt}\right).$ (4)

By transformation x = Xe ^iΩt , f _c = F _c e ^iΩt , it yields

$F_c=-\left(k_P+\frac{k_I}{i\Omega }+i\Omega k_D\right)X_{\mathrm{1,0}}=:-C\left(i\Omega \right)X_{\mathrm{1,0}}.$ (5)

Meanwhile there should be

$X_{\mathrm{1,0}}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right]{\left[\begin{array}{c}X\\ Q_x\end{array}\right]}_{\mathrm{1,0}}=:{\mathbf{E}}_m{\mathbf{Z}}_{\mathrm{1,0}}.$ (6)

Regarding the transfer direction and (1)–(3) and (5)–(6), the structure diagram of the controlled system can be drawn in Figure 2.

Figure 2:

The structural diagram of the controlled system [3].

From Figure 2, the overall transfer equation of the controlled system can be written as

$\begin{array}{c}{\mathbf{Z}}_{\mathrm{1,0}}={\mathbf{U}}_{\mathrm{1}}{\mathbf{Z}}_{\mathrm{2,1}}\\ ={\mathbf{U}}_{\mathrm{1}}\left({\mathbf{U}}_{\mathrm{2}}{\mathbf{U}}_{\mathrm{3}}{\mathbf{U}}_{\mathrm{4}}{\mathbf{Z}}_{\mathrm{4,0}}+{\mathbf{E}}_c{F}_c\right)\\ ={\mathbf{U}}_{\mathrm{1}}\left({\mathbf{U}}_{\mathrm{2}}{\mathbf{U}}_{\mathrm{3}}{\mathbf{U}}_{\mathrm{4}}{\mathbf{Z}}_{\mathrm{4,0}}-{\mathbf{E}}_{c}C\left(i\Omega \right){\mathbf{E}}_m{\mathbf{Z}}_{\mathrm{1,0}}\right).\end{array}$ (7)

which finally ends up with

${\mathbf{Z}}_{\mathrm{4,0}}={\left({\mathbf{U}}_{\mathrm{1}}{\mathbf{U}}_{\mathrm{2}}{\mathbf{U}}_{\mathrm{3}}{\mathbf{U}}_{\mathrm{4}}\right)}^{-\mathrm{1}}\left[{\mathbf{I}}_{\mathrm{2}}+C\left(i\Omega \right){\mathbf{U}}_{\mathrm{1}}{\mathbf{E}}_c{\mathbf{E}}_m\right]{\mathbf{Z}}_{\mathrm{1,0}}.$ (8)

This can be rewritten considering the boundary condition of the system ${\mathbf{Z}}_{\mathrm{1,0}}={\left[\begin{array}{c}X\\ \mathrm{0}\end{array}\right]}_{\mathrm{1,0}},{\mathbf{Z}}_{\mathrm{4,0}}={\left[\begin{array}{c}A\\ Q_x\end{array}\right]}_{\mathrm{4,0}}$ as

${\left[\begin{array}{c}A\\ Q_x\end{array}\right]}_{\mathrm{4,0}}=\left[\begin{array}{cc}{U}_{\mathrm{11}}& U_{\mathrm{12}}\\ U_{\mathrm{21}}& U_{\mathrm{22}}\end{array}\right]{\left[\begin{array}{c}X\\ \mathrm{0}\end{array}\right]}_{\mathrm{1,0}}.$ (9)

Thus the frequency response function from the disturbance to the output can be read as

$\frac{X_{\mathrm{1,0}}}{A}=\frac{\mathrm{1}}{U_{\mathrm{11}}}.$ (10)

Setting k _P = k _I = k _D = 0, one can acquire the formula for the uncontrolled system.

From the classical approaches for modeling the transfer function of a passive suspension system using ordinary differential equations, the equations which describe the system can be found in [14], and the transfer equation is

$G\left(s\right)=\frac{K_{\mathrm{4}}{C}_{\mathrm{2}}}{m_{\mathrm{3}}{m}_{\mathrm{1}}}\left(s+\frac{K_{\mathrm{2}}}{C_{\mathrm{2}}}\right)\times (s^{\mathrm{4}}+\left(\frac{C_{\mathrm{2}}}{m_{\mathrm{3}}}+\frac{C_{\mathrm{2}}}{m_{\mathrm{1}}}\right)s^{\mathrm{3}}{+\left(\frac{K_{\mathrm{2}}}{m_{\mathrm{3}}}+\frac{K_{\mathrm{2}}}{m_{\mathrm{1}}}+\frac{K_{\mathrm{4}}}{m_{\mathrm{3}}}\right)s^{\mathrm{2}}+\frac{K_{\mathrm{4}}{C}_{\mathrm{2}}}{m_{\mathrm{3}}{m}_{\mathrm{1}}}s+\frac{K_{\mathrm{4}}{K}_{\mathrm{2}}}{m_{\mathrm{3}}{m}_{\mathrm{1}}})}^{-1}.$ (11)

The described system elements values are given in Table 1 [15].

Table 1:

Elements values of controlled car.

Name	Value	Company
m 3	50	kg
m 1	250	kg
K 2	16	kN/m
K 4	160	kN/m
C 2	1500	N·s/m

For the uncontrolled case, Bode plots for both classical method and MS-TMM results are shown in Figure 3, considering the elements values given in Table 1 and substituting into (10) after eliminating the controller gains, the transfer equation of MS-TMM system as a function of the excitation frequency Ω is

$G\left(j\Omega \right)=\left(K_{\mathrm{4}}{K}_{\mathrm{2}}+\left(j\Omega \right)C_{\mathrm{2}}{K}_{\mathrm{4}}\right)\times (m_{\mathrm{3}}{m}_{\mathrm{1}}{\Omega }^{\mathrm{4}}-j{\Omega }^{\mathrm{3}}{C}_{\mathrm{2}}\left(m_{\mathrm{3}}+m_{\mathrm{1}}\right)-{\Omega }^{\mathrm{2}}(K_{\mathrm{2}}{m}_{\mathrm{3}}+K_{\mathrm{2}}{m}_{\mathrm{1}}+K_{\mathrm{4}}{m}_{\mathrm{3}}){+\left(j\Omega \right)K_{\mathrm{4}}{C}_{\mathrm{2}}+K_{\mathrm{4}}{K}_{\mathrm{2}})}^{-1}.$ (12)

This equation represents the open loop uncontrolled part of the system enclosed in the solid blue line in Figure 2; results from Figure 3 indicate that both techniques are identical.

Figure 3:

Bode plots for MS-TMM and classical method.

3. Controller Parameters Tuning

This section illustrates different approaches for designing both PI and PID controllers and some methods for tuning their parameters. Considering the control scheme shown in Figure 2, at the beginning a PI controller is designed; from (5) the designed PI controller will be at the form of

$C\left(s\right)=k_P+\frac{k_I}{s}=\frac{k_{P}s+k_I}{s}.$ (13)

From the previous literature reviews, instability is the disadvantage of feedback control system. There is a risk that the closed-loop system becomes unstable when using feedback. Thus, analyzing the closed-loop stability of system is essential requirement for feedback control system. The mean obstacle is to indicate the controller parameters or gains stabilizing region [16].

Decomposing (12) numerator and denominator to their even and odd parts, the equation can be rewritten as

$G\left(j\Omega \right)=\frac{N_e\left(-{\Omega }^{\mathrm{2}}\right)+j\Omega N_o\left(-{\Omega }^{\mathrm{2}}\right)}{D_e\left(-{\Omega }^{\mathrm{2}}\right)+j\Omega D_o\left(-{\Omega }^{\mathrm{2}}\right)}.$ (14)

From the control scheme diagram shown in Figure 2, and considering the control law in (13), the closed-loop characteristic polynomial derived in (10) can be separated and solved to zero; then the control parameters can be written as [17]

$k_P=\frac{X\left(\Omega \right)U\left(\Omega \right)-Y\left(\Omega \right)R\left(\Omega \right)}{Q\left(\Omega \right)U\left(\Omega \right)-R\left(\Omega \right)S\left(\Omega \right)},$ (15) $k_I=\frac{Y\left(\Omega \right)Q\left(\Omega \right)-X\left(\Omega \right)S\left(\Omega \right)}{Q\left(\Omega \right)U\left(\Omega \right)-R\left(\Omega \right)S\left(\Omega \right)},$ (16)

where

$\begin{array}{c}Y\left(\Omega \right)=-\Omega D_e\left(-{\Omega }^{\mathrm{2}}\right),\hspace{1em}\hspace{1em}Q\left(\Omega \right)=-{\Omega }^{\mathrm{2}}{N}_{\mathrm{0}}\left(-{\Omega }^{\mathrm{2}}\right),\\ X\left(\Omega \right)={\Omega }^{\mathrm{2}}{D}_{\mathrm{0}}\left(-{\Omega }^{\mathrm{2}}\right),\hspace{1em}\hspace{1em}S\left(\Omega \right)=\Omega N_e\left(-{\Omega }^{\mathrm{2}}\right),\\ U\left(\Omega \right)=\Omega N_o\left(-{\Omega }^{\mathrm{2}}\right),\hspace{1em}\hspace{1em}R\left(\Omega \right)=N_e\left(-{\Omega }^{\mathrm{2}}\right).\end{array}$ (17)

To determine the missing expressions in (15), (16), and (17), the transfer function in (12) can be divided into odd and even portions for both numerator and denominator, and after the comparison with (14), the following variables can be substituted as

$\begin{array}{c}{N}_e\left(-{\Omega }^{\mathrm{2}}\right)=K_{\mathrm{1}}{K}_{\mathrm{3}},\hspace{1em}\hspace{1em}{N}_o\left(-{\Omega }^{\mathrm{2}}\right)=K_{\mathrm{1}}{C}_{\mathrm{3}},\\ D_e\left(-{\Omega }^{\mathrm{2}}\right)=m_{\mathrm{2}}{m}_{\mathrm{4}}{\Omega }^{\mathrm{4}}-(K_{\mathrm{1}}{m}_{\mathrm{4}}+K_{\mathrm{3}}{m}_{\mathrm{2}}+K_{\mathrm{3}}{m}_{\mathrm{4}}){\Omega }^{\mathrm{2}}+K_{\mathrm{1}}{K}_{\mathrm{3}},\\ D_o\left(-{\Omega }^{\mathrm{2}}\right)=-C_{\mathrm{3}}\left(m_{\mathrm{2}}+m_{\mathrm{4}}\right){\Omega }^{\mathrm{2}}+K_{\mathrm{1}}{C}_{\mathrm{3}}.\end{array}$ (18)

For the controller design problems which based on gain and phase margins design, it is inevitable to estimate the stabilizing region of the controller gains. The stability boundary locus in proportional and integral gains plane can be obtained from (15) and (16) after substituting the missing terms from (17) and (18), then verifying the presented model using the table of contents in Table 1.

The region of stability of the proportional gain k _P and the integral gain k _I that cause the closed-loop polynomial to satisfy Hurwitz stability test is shown in Figure 4, where line k _I = 0 divides the controller gains plane into stable and unstable regions.

Figure 4:

Stability region for frequency range Ω → [0, 56.7].

For designing of controller to meet a certain demanded gain and phase margins which are important frequency domain requirements for user specifications, consider a control function:

$G_c\left(j\Omega \right)=M{e}^{-j\phi },$ (19)

where M and φ are gain and phase margins, respectively.

Substituting (17) and (18) into (15) and (16) and considering (19) yields to

$\begin{array}{c}{K}_P=\frac{N_A\cos \left[\phi \right]+N_B\sin \left[\phi \right]}{-M(N_e^{\mathrm{2}}+{\Omega }^{\mathrm{2}}{N}_o^{\mathrm{2}})},\\ K_I=\frac{\Omega N_B\cos \left[\phi \right]-\Omega N_A\sin \left[\phi \right]}{-M(N_e^{\mathrm{2}}+{\Omega }^{\mathrm{2}}{N}_o^{\mathrm{2}})},\end{array}$ (20)

where

$N_A=\left({\Omega }^{\mathrm{2}}{N}_o{D}_o+N_e{D}_e\right),\hspace{1em}{N}_B=\Omega \left(N_o{D}_e-N_e{D}_o\right).$ (21)

Using these formulae makes it easy to satisfy certain conditions and find the stabilization region for the controller parameters. Assuming two cases, one as the phase margin is greater than 45° and gain margin is greater than 1, another as the phase margin is greater than 60° and gain margin is greater than 1.

Figure 5 shows two regions of stabilities for two different conditions. The blue circles locus shows the stability region for a demanded conditions gain margin =1 and phase margin 45°, the maximum frequency for this region is 41.5 rad/sec, the red solid locus shows another design conditions with gain margin =1 and phase margin 60°, and the maximum frequency for this region is 31.5 rad/sec.

Figure 5:

Stability region for two cases.

It is an applicable method to achieve user specified gain and phase margins to define the values of controller parameters with the ability to justify the maximum frequency without the need of sweeping over the parameters and also without using programming for solving inequality problems. This method can be expanded to estimate the stabilizing region of PID controller gains as shown in [17].

From Figure 6, the plotted line in red solid color indicates the output response of the uncontrolled system and also shows the stability region of frequencies shown in red points in the range Ω ∈ [10.2, 15.6] rad/sec, while the blue dashed locus indicates the system response after applying the PI controller which is chosen from the stability region shown in Figure 6 as K _P = 1 and K _I = 2.5; it is obvious that the stability region of frequencies shown in blue squares is expanded in the range Ω ∈ [0, 15.5] rad/sec.

Figure 6:

Bode plot for the controlled and uncontrolled system.

The relationship between the phase margin Φ _M , the overshoot (OS), and the damping ratio ζ can be expressed as:

${\Phi }_M=\mathrm{9}{\mathrm{0}}^{°}-\text{ta}{\text{n}}^{-\mathrm{1}}\frac{\sqrt{-\mathrm{2}{\zeta }^{\mathrm{2}}+\sqrt{\mathrm{1}+\mathrm{4}{\zeta }^{\mathrm{4}}}}}{\mathrm{2}\zeta },$ (22)

where

$\text{OS}\%=e^{-\left(\zeta \pi /\sqrt{\mathrm{1}-{\zeta }^{\mathrm{2}}}\right)}\times \mathrm{100}\%.$ (23)

The relation between Φ _M and ζ is plotted in Figure 7.

Figure 7:

Phase margin versus damping ratio.

So it is noticeable to get different overshoot values due to the changes of the phase margins with the controller design parameters; Figure 9 shows the effect of changing these values as follow.

In Figure 8 the red dash-dot line indicates the uncontrolled system step response; the overshoot is about 46%, and choosing a test point from Figure 3 with values of K _P = 1 and K _I = 25, it is clear that it locates outside the regions of stabilities from the used technique, so it is clear from the figure that the green solid line response overshoot value is about 55%, and unlike the line in blue dot color it is clear that the overshoot value is improved to be 7% when the test point is chosen as K _P = 1 and K _I = 7.5.

Figure 8:

Step response of the system.

Figure 9:

Step response of the system.

Returning back to (10) and substituting jΩ = s, another technique is applied in this section; a developed methodology by Shamsuzzoha and Skogestad requires only one closed-loop step test to obtain PI controller setting. In this method it is also simple to obtain the PID tuning parameters [11, 12].

Considering an ideal PID controller in the form, as shown in (5)

$c\left(s\right)=K_c\left(\mathrm{1}+\frac{\mathrm{1}}{{\tau }_{I}s}+{\tau }_{D}s\right)⟺c\left(s\right)=\left(K_P+\frac{K_I}{s}+K_{D}s\right).$ (24)

K _c is the controller gain and τ _I ,τ _D are the controller integral and derivative time, respectively.

The defined parameters can be calculated using the following basic formulae:

$\begin{array}{c}{K}_c=\frac{\mathrm{2}\tau +\theta }{\mathrm{3}k\theta },\hspace{1em}\hspace{1em}{\tau }_I=\min \left\{\left(r+\frac{\theta }{\mathrm{2}}\right),\mathrm{8}\theta \right\},\\ {\tau }_D=\frac{r\theta }{\mathrm{2}\tau +\theta },\end{array}$ (25)

where k is the gain, τ is the time lag constant, and θ is the time delay. The controller gain is only depending on the overshoot, while the integral and derivative times are function of the system peak time. Applying this approach, the tuning results simulation for the designed problem assuming different values for τ/θ is indicated in Table 2.

Table 2:

Controller tuning parameters.

τ/θ	Controller parameters
	K P	K I	K D
0.1	0.4000	0.6667	0.0333
0.5	0.6667	0.6667	0.1667
2.5	2.0000	0.6667	0.8333
10	7.0000	0.8750	3.3333
20	13.6667	1.7083	6.6667
50	33.6667	4.2083	16.6667
100	67.0000	8.3750	33.3333

An improved definition of (25) parameters were concluded in Shamsuzzoha et al. [12], a closed-loop test with P-controller was applied to adjust overshoot value of 0.3021 for P-only control then assigned the other parameters which leaded to gain K _c = 2.5 with τ _I = 0.2040 and peak time 0.153 for the controller gains.

In spite of Zeigler-Nichols methods’ drawbacks which are mentioned in previous sections, it is still one of the easiest methods for determining the proportional gain K _P , integral time T _i , and derivative time T _d based on the transient response of a given system. Zeigler-Nichols suggested the values of the controller parameters according to the shown formulae in Table 3.

Table 3:

Zeigler-Nichols tuning rules.

Type of controller	K P	T i	T d
P	0.5Kcr	∞	0
PI	0.45Kcr	0.83 Pcr	0
PID	0.6Kcr	0.5 Pcr	0.125 Pcr

For the systems which have a mathematical model, root-locus method can be used to find the critical gain K_cr and the sustained oscillation frequency ω_cr; then get P_cr = 2π/ω_cr. Also the systems’ open loop system (passive suspension) and closed-loop system (Active suspension) can be analyzed in frequency domain to find the bandwidth, gain, and phase margins, and then the design parameters can be found easily.

Considering the previous controlled car problem, the results are found as follows: the critical gain is 5.1, and the sustained oscillation frequency is 56.7 rad/sec. It is the same result shown in Figure 4, where the stability frequency for the case of assuming gain margin 1 and phase margin 0 is Ω → [0, 56.7].

From the previous values, the P, PI, and PID Zeigler-Nichols controller values can be tabulated in Table 4.

Table 4:

Zeigler-Nichols tuning values.

Type of controller	K P	T i	T d
P	2.55	∞	0
PI	2.295	0.091976	0
PID	3.06	0.036568	0.013852

Applying the controller parameters obtained from the different tuning techniques, the system response is plotted in Figure 9.

From Figure 9 and related to the different tuning methods, the systems’ time of the first peak, the peak values, and the overshoot values can be tabulated in Table 5.

Table 5:

System response analysis parameters.

4. Simulations and Results

After defining the controller parameters, it is necessary to test the system performance for some situations. The suspension response for a vehicle coming out of a pothole can be simulated as a step input with amplitude of 0.08 as [18]:

$\text{Amp}=\{\begin{array}{cc}\mathrm{0}⟹t\le \mathrm{0},\hfill & \hfill \\ \mathrm{0.08}⟹t\ge \mathrm{0}.\hfill & \hfill \end{array}$ (26)

The system output displacement is shown in Figure 10.

Figure 10:

The output displacement from the system.

Figure 10 shows that the peak value of the passive suspension system in red solid line is nearly 0.116 m although the input step amplitude is 0.08 m; that is, the overshoot is 45%, and after using the estimated values from the described techniques, the overshoot and response of the system improved although the transient time is increased for some of them but the displacement action is smoothed; Zeigler-Nichols results are plotted in dash-dot line, Shamsuzzoha results are plotted in dash lines, and PI controller is plotted in x marker line.

Another simulation model for random roughness road is applied where the road model is expressed as a differential equation with the formula [19]:

$\ddot{Z}\left(t\right)+avZ\left(t\right)=W\left(t\right),$ (27)

where $\ddot{Z}(t)$ is the road random incentive and the constant “a” depends on road grade; the values of “a” are shown in Table 6.

Table 6:

“a” value under different road grades.

Road grade	A	B	C	D	E
a	0.132	0.1303	0.12	1.1007	0.09

W(t) is the white noise, with zero mean value and amplitude equals 1, and v is the vehicle speed.

Focusing for more details can be shown in Figure 12.

The desired problem is validated using the results of road simulation of grade B and assuming the velocity is 50 km/hr; the results show the effectiveness of the used controllers to reduce the displacement motion as shown in Figure 11 and zoomed in Figure 12, also results illustrated that different control techniques can be applied to MS-TMM, results of the passive suspension response are plotted in blue solid line, Zeigler-Nichols result is plotted in blue dot plot, Shamsuzzohas’ result is plotted in green dash-dot plot, and finally the PI controls’ result is plotted in red dash plot.

Figure 11:

Output displacement of the system.

Figure 12:

Output displacement zoomed in.

5. Conclusion and Future Work

As the increments of complexity of multibody systems and the development of their design and control methods, the need for more elegant formulations of the equations of motion becomes an issue of paramount importance. The results showed identical values for MS-TMM results with the ordinary classical methods. Thus, ongoing researches, MS-TMM strategies, and its related topologies are highly recommended for analyzing and solving the controlled systems due to easy formulations, being systematic to apply, and simple program coding with less computational time. The results also showed that different controller design methods can be impeded with MS-TMM which leads to more availability for other controllers applications; one depends on the frequency response analysis, another considers the step input response of the system, and the last technique can be applied for both time and frequency domains if the model has a known mathematical model. All controllers’ parameters values are applied to a model of controlled car system with brief analysis for the obtained results. The future work will concern with deeply impeding different control techniques with MS-TMM to develop new techniques for the more complex multibody controlled systems considering nonstationary inputs by taking into consideration the system's equations of motion with studying practical systems in frequency domain due to its advantage of indicating inputs and outputs easily such as power spectral density measurements for such kinds of active vibration isolation techniques.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of the paper.