Vibration control of Stewart platform based on an improved nonlinear proportional-derivative control method

Abstract

This paper proposes an improved switching nonlinear proportional-derivative (S-NPD) control algorithm for a six-degree-of-freedom Stewart platform which has strong nonlinearity and uncertainty characteristics, designed to achieve vibration isolation and suppression over a wide bandwidth range while solving the problem of the difficulty in parameter tuning. The controller is designed based on the idea of variable-structure control, which achieves independent real-time adjustments of the control gains of the six legs and takes the sign of the proportional and differential errors as the switching function to realize the feed-forward correction of the proportional gain. This control input can track the expected value more precisely, thus strengthening the system's ability to cope with high dynamics. To validate the results of the proposed approach, numerical simulations are compared with those obtained from PD and NPD controls. The findings confirmed that the designed S-NPD controller shows commendable stability and robustness, and the vibration attenuation of the upper platform under low- and medium-frequency conditions is substantially improved. In addition, it obtains smoother vibration waveforms and lower control force under random disturbance conditions. Consequently, it is proved that the S-NPD control can obtain better regulation accuracy. Compared with the PD and NPD controls, the S-NPD control achieves up to 64% and 31.4% relative performance improvement in vibration attenuation rate, respectively.

Introduction

The Stewart platform has a high response speed and large range of motion, and the use of the Stewart parallel mechanism for vibration isolation is one of the most crucial applications currently.^1–5 In most studies, the Stewart platform is applied on the ground where the environment is characterized by large mass and mechanical simplicity. However, as the platform application scenarios continue to expand, additional requirements on the platform structural design, such as lightweight construction and reduced volume. Simultaneously, the corresponding control algorithms must be improved to meet the evolving demands of diverse application scenarios.^6–9

The dynamics of the Stewart mechanism include workspace dynamics and joint space dynamics equations in which the workspace dynamics equations are used to describe the dynamics relationship between the driving force of each leg and the platform's positional parameters, while the joint space dynamics equations are used to describe the dynamics relationship between the driving force of each leg and the change in the loudness displacement between the upper and lower parts of each leg. Consequently, control strategies for these dynamic models can correspondingly categorized into two types: workspace-based control and joint-space-based control. Workspace-based control is a multi-input and multi-output approach directly based on the positional attitude deviation of the platform in inertial space. However, joint-space-based control translates the end attitude position information into the joint space of the legs. Owing to the mutual coupling of the Stewart mechanism, the workspace-based positive solution is complicated to solve, and the joint-space-based decentralized control strategy for each actuator is usually adopted.

In current academic research and engineering applications, the proportional-integral-derivative (PID) control algorithm remains a widely adopted strategy for the Stewart platform due to its ease-to-implement characteristics and good control performance. However, with the increasing complexity of the control object structure, its strong nonlinearity, strong coupling and uncertainty content also gradually increase.^8,9 Despite these problems, the performance of PID is severely constrained, and its controller can no longer meet the performance requirements under some high-precision working conditions.

To overcome this shortcoming, based on the classical control theory of PID, scholars have made a series of improvements and research, and many strategies and theories have been put forward one after another. The literatures^10–12 combine fuzzy control, neural networks, and other contents with PID control theory to form an improved PID control that adds functions such as parameter online tuning, and the control performance achieves a favorable improvement effect. Fuzzy PID control requires setting the fuzzy structure first. If the fuzzy structure is not set accurately, it will directly affect the control performance of the system. ¹⁰ The 2-DOF PID enhances the controller's anti-interference ability by adding a set of independent parameters to be adjusted, however at the same time, more undetermined parameters also increase the system's complexity, which reduces its application convenience. ¹³ The nonlinear proportional derivative (NPD) control proposed by Prof. Han Jingqing,^14,15 using the advantage of nonlinear function, to a certain extent, solved the problem that PD control does not apply to complex working conditions, and provided a new method of thinking for the improvement of PID control. In Professor Han's research, it was pointed out that the feedback of error integral in PID control has a certain effect on suppressing constant disturbances. However, when there is no disturbance, integral feedback often deteriorates the dynamic characteristics of the closed-loop. For time-varying disturbances, the suppression ability of integral feedback is not significant enough, and it will increase the complexity of the system and the difficulty of parameter tuning. Therefore, in his subsequent research and improvement (NPD and related research), he omitted integral feedback and only retained PD control. Zheng et al. ¹⁶ derived and analyzed the role of the δ parameter in the NPD control equation, pointed out that a very small value of δ may cause potential jitter, and provided a theoretical basis for the assignment of this parameter. To further improve the performance of NPD control, it is essential to introduce new features to the control algorithm while retaining the advantages of NPD control. The nonlinear function in the error term is well-established and challenging to modify; therefore, the improvement mainly focuses on the control gain term. Two improvement ideas emerge: the first involves the combination of the previously mentioned neural network method. Different neural network control algorithms are distinguished by different neural network functions, among which the more common ones include Sigmoid, ¹⁷ Radical Basis Function (RBF) ¹⁰ and Brain Emotional Learning Based Intelligent Controller (BELBIC). ¹⁸ Neural networks using Sigmoid basis functions have strong fault tolerance and function approximation capabilities, but the energy of Sigmoid functions is infinite, requiring infinite weights, and the network size is difficult to determine. The RBF has a faster learning rate but also requires more radial basis function neurons. Neuroendocrine feedback through simulating hormone regulation, however, there may be slow or overshoot issues while it has self adjustment function. ¹⁹ BELBIC implements neural network patterns by simulating the learning structure of the human brain, which also requires a complex and time-consuming learning process. While this reduces parameterization difficulty and improves the algorithm's applicability, it introduces complexity to the system. Although the neural network reduces the difficulty of parameterization, it also increases the complexity of the system simultaneously. The neural network needs to undergo several training and adjustments to achieve the expected control effect, which is in contradiction with the PD/NPD control. This contradicts the advantages of a simple structure and easy application.

The second idea is based on adaptive gain, allowing real-time adjustment to control gain for further system performance enhancement.^20–23 Compared with the first idea, the adaptive gain has better rapidity, and the control performance of the system can be improved in real-time only by setting the function for the control gain in advance. Armstrong et al. ²⁴ set up two groups of PD gain parameters with different sizes, taking the sign of proportional error and its derivative as the switching condition for the control system of the servo motor and using different parameters under the same and different sign conditions, to realize a simple switching control. However, the simplicity of the structural setup leads to limitations in the robustness of the system. Sliding mode control (SMC) improves the controller's ability to cope with uncertainty by setting the sliding mode surface and switching conditions. Meng et al. ²⁵ proposed a model-based sliding mode controller design that determines the amount of gain of the control in terms of the system's damping and intrinsic frequency. This design reduces the difficulty of parameter tuning and at the same time increases the requirements for the model information and accuracy. Boudjedir et al. ²⁶ proposed a method to combine NPD control with SMC, which is oriented to the application scenario of delta parallel robots, and improved NPD control by taking the proportional and differential error as the switching condition to realize the dynamic change of NPD control parameters. The simulation results show a significant improvement in the effectiveness of the improved control algorithm in tracking a given curve. However, the control law does not take into account the feedback loop, which may cause the controller to be susceptible to conditions such as vibration jittering. Therefore, it is necessary to introduce a feedback link to further strengthen its robustness.

In this paper, the theory of variable structure control is applied to the parameter rectification of NPD control, and the switching nonlinear proportional-derivative (S-NPD) control algorithm is applied to the Stewart platform. Based on the theory of the variable structure, the S-function is proposed to be applied to parameter tuning, and the error variation is estimated by the proposed S-function to enable the advance correction of parameters. In this way, the correction of the control gain becomes more stable and precise. Compared with neural network control, the improved NPD control does not require additional neural network construction, nor does it require a long period of learning. Parameter adjustments are made in real-time during the working process. The subsequent sections of this paper are organized as follows. System modeling, including kinematic and dynamic modeling is introduced in Dynamic modeling section. The improved NPD controller design is introduced in Improved NPD controller design section. Experiments based on a six-degree-of-freedom Stewart platform and results analysis of PD, NPD and improved S-NPD are introduced in Numerical simulation and analysis section, which proves that the improved S-NPD controller can realize a higher vibration attenuation rate of the upper platform and lower input driving force. Finally, Conclusion section presents the conclusion of the research.

Dynamic modeling

Coordinate system definition

In this paper, the following assumptions are made for the models: the lower leg serves as the driving leg with the upper leg as the driven leg; the six legs are identical and orthogonal to each other for easy decoupling; both the upper and lower platforms are treated as rigid bodies, and the load is located on top of the upper platform; ball hinges connect the hinge points between the outriggers and the platform while sliding hinges the upper and lower outriggers at the expansion and contraction points.

The schematic diagram of the parallel platform model, defining the global coordinate system $Oxyz$ is presented in Figure 1. The initial configuration places the centroid of the lower platform at the origin of the global coordinate system. The centroid coordinate systems of the upper and lower at the origin of the global coordinate system. The centroid coordinate system of the lower platform is firmly attached to their respective solid platforms, with origins at points P and B, respectively. Initially, the direction of each axis in the centroid coordinate system of the upper and lower platforms aligns with that of the global coordinate system.

Figure 1.

Mechanical model drawing of parallel platform/structural drawing of parallel platform.

To determine the number of degrees of freedom in a parallel platform, the methodology outlined in the literature ²⁷ is employed. The three rotational degrees of freedom of the upper and lower platforms involve attitude changes, which are usually described using attitude parameters such as Euler angles and quaternions.

According to Leonhard Euler, ²⁸ any rotation can be expressed as a combination of three angles around three axes of rotation in sequence, known as Euler angles. To transform from coordinate system A to coordinate system B using the Euler angle parameter, for system A the rotation angle around the x-axis is $\alpha$ , while the y-axis is $\beta$ and the z-axis is $\gamma$ .

Derivation of the angular velocity of rotation for the B coordinate system yields the following formula for the angular velocity described by quaternions:

$$\mathit{\omega }=2\mathit{G}\dot{\mathit{q}}$$

(1)

where $\mathit{G}$ is a matrix defined by quaternions:

$$\mathit{G}=\left[\begin{array}{cc}\begin{array}{cc}-q_1& q_0\\ -q_2& -q_3\end{array}& \begin{array}{cc}{q}_3& -q_2\\ q_0& q_1\end{array}\\ \begin{array}{cc}-q_3& q_2\end{array}& \begin{array}{cc}-q_1& q_0\end{array}\end{array}\right]$$

(2)

The deviation quaternion ${\mathit{q}}_e$ can be expressed as:

$${\mathit{q}}_e=\left[\begin{array}{cc}\begin{array}{cc}{q}_{d1}& q_{d2}\\ -q_{d2}& q_{d1}\end{array}& \begin{array}{cc}{q}_{d3}& q_{d4}\\ q_{d4}& -q_{d3}\end{array}\\ \begin{array}{cc}-q_{d3}& -q_{d4}\\ -q_{d4}& q_{d3}\end{array}& \begin{array}{cc}{q}_{d1}& q_{d2}\\ -q_{d2}& q_{d1}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{q}_1\\ q_2\end{array}\\ \begin{array}{c}{q}_3\\ q_4\end{array}\end{array}\right]$$

(3)

Kinematic modeling

The main geometric parameters of the platform are described in the figure below. In Figure 2, ${\mathit{r}}_b$ is the position vector of the geometric center B of the lower platform, the position vector ${\mathit{r}}_p$ of the geometric center P of the upper platform ${\mathit{p}}_i$ is the vector from the upper platform from center to leg i with the upper platform constraint position, and ${\mathit{b}}_i$ is the vector from the lower platform from center to leg i with the lower platform constraint position.

Figure 2.

Schematic description of kinematics.

Based on the above kinematic description, the position vectors at the hinge points of the upper and lower platforms can be obtained as:

$${\mathit{r}}_{pi}={\mathit{r}}_p+{\mathit{p}}_i,{\mathit{r}}_{bi}={\mathit{r}}_b+{\mathit{b}}_i$$

(4)

The vector of leg i pointing from the lower platform connection point to the upper platform connection point is

$${\mathit{l}}_i={\mathit{r}}_{pi}-{\mathit{r}}_{bi}={\mathit{r}}_p+{\mathit{p}}_i-{\mathit{r}}_b-{\mathit{b}}_i$$

(5)

In the global coordinate system, the vector of leg i can be further expressed as:

$${\mathit{l}}_i^o={\mathit{r}}_p^o-{\mathit{r}}_b^o+{\mathit{A}}_o^p{\mathit{p}}_i^p-{\mathit{A}}_o^B{\mathit{b}}_i^B$$

(6)

where superscript O represents the coordinate component of the vector under the global coordinate system, superscript P represents the coordinate component of the vector under the upper platform coordinate system, superscript B represents the coordinate component of the vector under the lower platform coordinate system, ${\mathit{A}}_o^B$ denotes the coordinate transformation matrix from the coordinate system O to the coordinate system B, and ${\mathit{A}}_o^P$ denotes the coordinate transformation matrix from the coordinate system O to the coordinate system P.

According to equation (5), the length of the leg i which represents the kinematic inverse solution $l_i$ and the unit vector ${\mathit{e}}_i$ can be obtained as:

$$l_i=|{\mathit{l}}_i|,{\mathit{e}}_i={\mathit{l}}_i/l_i$$

(7)

Derivation of equation (4) gives the velocity and acceleration at the hinge point of the upper and lower platforms as:

$${\dot{\mathit{r}}}_{pi}={\dot{\mathit{r}}}_p+{\mathit{\omega }}_p\times {\mathit{p}}_i,{\dot{\mathit{r}}}_{bi}={\dot{\mathit{r}}}_b+{\mathit{\omega }}_b\times {\mathit{b}}_i$$

(8)

$${\ddot{\mathit{r}}}_{pi}={\ddot{\mathit{r}}}_p+{\mathit{\omega }}_p\times ({\mathit{\omega }}_p\times {\mathit{p}}_i)+{\mathit{\alpha }}_p\times {\mathit{p}}_i,{\ddot{\mathit{r}}}_{bi}={\ddot{\mathit{r}}}_b+{\mathit{\omega }}_b\times ({\mathit{\omega }}_b\times {\mathit{b}}_i)+{\mathit{\alpha }}_b\times {\mathit{b}}_i$$

(9)

where ${\mathit{\omega }}_p$ and ${\mathit{\omega }}_b$ are the rotational angular velocities. ${\mathit{\alpha }}_p$ and ${\mathit{\alpha }}_b$ are the rotational angular accelerations.

Subtracting the two terms of equation (8) and projecting them in the direction of the leg, the rate of change of the leg can be obtained as follows:

$${\dot{l}}_i=({\dot{\mathit{r}}}_{pi}-{\dot{\mathit{r}}}_{bi})\cdot {\mathit{e}}_i=[({\dot{\mathit{r}}}_p+{\mathit{\omega }}_p\times {\mathit{p}}_i)-({\dot{\mathit{r}}}_b+{\mathit{\omega }}_b\times {\mathit{b}}_i)]\cdot {\mathit{e}}_i$$

(10)

Organizing equation (10) leads to:

$${\dot{l}}_i=[{\mathit{e}}_i^T,{({\mathit{p}}_i\times {\mathit{e}}_i)}^T]\left[\begin{array}{c}{\dot{\mathit{r}}}_p\\ {\mathit{\omega }}_p\end{array}\right]-[{\mathit{e}}_i^T,{({\mathit{b}}_i\times {\mathit{e}}_i)}^T]\left[\begin{array}{c}{\dot{\mathit{r}}}_b\\ {\mathit{\omega }}_b\end{array}\right]={\mathit{J}}_{pi}^T\left[\begin{array}{c}{\dot{\mathit{r}}}_p\\ {\mathit{\omega }}_p\end{array}\right]-{\mathit{J}}_{bi}^T\left[\begin{array}{c}{\dot{\mathit{r}}}_b\\ {\mathit{\omega }}_b\end{array}\right]$$

(11)

where ${\mathit{J}}_p={\left[\begin{array}{ccc}{\mathit{e}}_1& \cdots & {\mathit{e}}_6\\ {\tilde{\mathit{p}}}_1{\mathit{e}}_1& \cdots & {\tilde{\mathit{p}}}_6{\mathit{e}}_6\end{array}\right]}_{6\times 6},{\mathit{J}}_b={\left[\begin{array}{ccc}{\mathit{e}}_1& \cdots & {\mathit{e}}_6\\ {\tilde{\mathit{b}}}_1{\mathit{e}}_1& \cdots & {\tilde{\mathit{b}}}_6{\mathit{e}}_6\end{array}\right]}_{6\times 6}$ , ${\mathit{J}}_p$ and ${\mathit{J}}_b$ are the Jacobi matrices for the upper and lower platforms respectively.

In the same way, it can be deduced that:

$${\ddot{l}}_i=\left[{\mathit{e}}_i^T,({\mathit{p}}_i\times {\mathit{e}}_i{)}^T]\left[\begin{array}{c}{\ddot{\mathit{r}}}_p\\ {\mathit{\alpha }}_p\end{array}\right]-\right[{\mathit{e}}_i^T,({\mathit{b}}_i\times {\mathit{e}}_i{)}^T]\left[\begin{array}{c}{\ddot{\mathit{r}}}_b\\ {\mathit{\alpha }}_b\end{array}\right]={\mathit{J}}_{pi}^T\left[\begin{array}{c}{\ddot{\mathit{r}}}_p\\ {\mathit{\alpha }}_p\end{array}\right]-{\mathit{J}}_{bi}^T\left[\begin{array}{c}{\ddot{\mathit{r}}}_b\\ {\mathit{\alpha }}_b\end{array}\right]$$

(12)

Dynamic modeling

The legs are treated as rigid for dynamics analysis, considering that the six-degree-of-freedom parallel platform operates in a space environment; therefore, the effect of gravity is ignored. In addition, friction during the relative motion of the legs has little effect on the overall motion of the parallel platform and is therefore ignored in the modeling process. ²⁸

The equations of rotational balance of the upper and lower legs are:

$$l_{ui}{\mathit{e}}_i\times {\mathit{F}}_{ui}-M_i-(l_U-l_{ui}){\mathit{e}}_i\times {\mathit{F}}_{ci}={\mathit{I}}_{ui}{\mathit{\alpha }}_i+{\mathit{\omega }}_i\times {\mathit{I}}_{ui}{\mathit{\omega }}_i$$

(13)

$$-l_{di}{\mathit{e}}_i\times {\mathit{F}}_{di}+M_i-(l_D-l_{di}){\mathit{e}}_i\times {\mathit{F}}_{ci}={\mathit{I}}_{di}{\mathit{\alpha }}_i+{\mathit{\omega }}_i\times {\mathit{I}}_{di}{\mathit{\omega }}_i$$

(14)

Where ${\mathit{F}}_{ui}$ is the binding force on the upper leg, ${\mathit{F}}_{di}$ is the lower leg and ${\mathit{F}}_{ci}$ is the sliding hinge. $m_{ui}$ and $m_{di}$ are the masses of the upper and lower legs. $l_{ui}$ and $l_{di}$ are the distances from the upper and lower platform hinge points i to centers of gravity of the upper and lower legs i. $l_U$ and $l_D$ are the lengths of the upper and lower legs, $M_i$ the rotational moment at the sliding hinge, $I_{ui}$ , and $I_{di}$ the moment of inertia of the upper and lower legs about the centroid.

Thereby obtain

$$\{\begin{array}{cc}& {\mathit{e}}_i\cdot {\mathit{F}}_{ui}=m_{ui}{\mathit{e}}_i\cdot {\mathit{a}}_{ui}-{\mathit{e}}_i\cdot {\mathit{F}}_{ci}\\ & l_i({\mathit{e}}_i\cdot {\mathit{F}}_{ui}){\mathit{e}}_i-l_i\times {\mathit{F}}_{ui}={\mathit{e}}_i\times ({\mathit{I}}_{ui}+{\mathit{I}}_{di}){\mathit{\alpha }}_i+{\mathit{e}}_i\times {\mathit{\omega }}_i\times ({\mathit{I}}_{ui}+{\mathit{I}}_{di}){\mathit{\omega }}_i\end{array}$$

(15)

Let

$${\mathit{Q}}_{pi}=m_{ui}{\mathit{e}}_i{\mathit{e}}_i^T-{\displaystyle \frac{1}{l_i^2}[m_{ui}{(l_i-l_{ui})}^2+m_{di}{l}_{di}^2]{\tilde{\mathit{e}}}_i{\tilde{\mathit{e}}}_i-{\displaystyle \frac{1}{l_i^2}{\tilde{\mathit{e}}}_i({\mathit{I}}_{ui}+{\mathit{I}}_{di}){\tilde{\mathit{e}}}_i}}$$

(16)

$${\mathit{Q}}_{bi}={\displaystyle \frac{1}{l_i^2}[m_{ui}(l_i-l_{ui})l_{ui}+m_{di}{l}_{di}(l_i-l_{di})]{\tilde{\mathit{e}}}_i{\tilde{\mathit{e}}}_i-{\displaystyle \frac{1}{l_i^2}{\tilde{\mathit{e}}}_i({\mathit{I}}_{ui}+{\mathit{I}}_{di}){\tilde{\mathit{e}}}_i}}$$

(17)

$${\mathit{Q}}_i=-{\displaystyle \frac{1}{l_i}{\tilde{\mathit{e}}}_i({\mathit{I}}_{ui}+{\mathit{I}}_{di}){\mathit{e}}_i}$$

(18)

where ${\mathit{Q}}_i$ represents the inertia term generated during the rotation of the leg, ignoring this term gives the force expression at the upper hinge point as:

$${\mathit{F}}_{ui}={\mathit{Q}}_{pi}({\ddot{\mathit{r}}}_p+{\mathit{\alpha }}_p\times {\mathit{p}}_i)+{\mathit{Q}}_{bi}({\ddot{\mathit{r}}}_b+{\mathit{\alpha }}_b\times {\mathit{b}}_i)$$

(19)

The six-degree-of-freedom parallel platform mainly relies on the sliding hinge at the leg to perform relative motion along the axial direction and realizes the translation or rotation of the upper platform through expansion and contraction. The whole leg force is equated to the main force output from the spring, damping, and actuator, and the gap friction between the legs is ignored. The equivalent model is shown as Figure 3.

Figure 3.

Equivalent model of leg.

The axial force in the leg i can be expressed as follows:

$$F_i=-k_i\cdot \mathrm{\Delta }{l}_i-c_i\cdot \mathrm{\Delta }{\dot{l}}_i+F_{ai}$$

(20)

where $k_i$ is rigidity, $c_i$ is friction, $F_{ai}$ is the axial force generated by the actuator, and $\mathrm{\Delta }{l}_i$ is the relative displacement of the upper and lower legs during expansion and contraction. $\mathrm{\Delta }{\dot{l}}_i$ is the relative speed of the upper and lower legs.

Substituting into the original formula to further simplify, the expression of the axial force of the leg is:

$$\mathit{F}=-k\mathit{E}({\mathit{J}}_p^T{\mathit{q}}_p-{\mathit{J}}_b^T{\mathit{q}}_b)-c\mathit{E}({\mathit{J}}_p^T{\dot{\mathit{q}}}_p-{\mathit{J}}_b^T{\dot{\mathit{q}}}_b)+{\mathit{F}}_{ai}$$

(21)

where $\mathit{E}$ is an identity matrix with dimension of $6\times 6$ . In real space applications, the upper platform model is often omitted and the legs are directly mounted on the load base; thus, the upper platform and load are merged in the derivation of the formula, and the payload mass and inertia are regarded as m and ${\mathit{I}}_m$ , while ${\mathit{r}}_c$ is the position vector of the payload in the global coordinate system.

The equilibrium equation for the parallel platform is:

$$\left[m\mathit{E}+\sum _{i=1}^6{\mathit{Q}}_{pi}m{\tilde{\mathit{r}}}_c^T+\sum _{i=1}^6{\mathit{Q}}_{pi}{\tilde{\mathit{p}}}_i^T\right]{\ddot{\mathit{x}}}_p=\left[\sum _{i=1}^6{\mathit{Q}}_{bi}\sum _{i=1}^6{\mathit{Q}}_{bi}{\tilde{\mathit{b}}}_i^T\right]{\ddot{\mathit{x}}}_b+[{\mathit{e}}_1\cdots {\mathit{e}}_6]\mathit{F}$$

(22)

The equations of rotational balance are:

$$\left[m{\tilde{\mathit{r}}}_c+\sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{pi}{\mathit{I}}_m+m{\tilde{\mathit{r}}}_c{\tilde{\mathit{r}}}_c^T+\sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{pi}{\tilde{\mathit{p}}}_i^T\right]{\ddot{\mathit{x}}}_p=\left[\sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{bi}\sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{bi}{\tilde{\mathit{b}}}_i^T\right]{\ddot{\mathit{x}}}_b+[{\tilde{\mathit{p}}}_1{\mathit{e}}_1\cdots {\tilde{\mathit{p}}}_6{\mathit{e}}_6]\mathit{F}$$

(23)

where $\mathit{F}$ is the force at the sliding hinge of the outrigger. Let

$${\mathit{M}}_p=\left[\begin{array}{c}mE+\sum _{i=1}^6{\mathit{Q}}_{pi}m{\tilde{\mathit{r}}}_c^T+\sum _{i=1}^6{\mathit{Q}}_{pi}{\tilde{\mathit{p}}}_i^T\\ m{\tilde{\mathit{r}}}_c+\sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{pi}{\mathit{I}}_m+m{\tilde{\mathit{r}}}_c{\tilde{\mathit{r}}}_c^T+\sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{pi}{\tilde{\mathit{p}}}_i^T\end{array}\right]$$

(24)

$${\mathit{M}}_b=\left[\begin{array}{cc}\sum _{i=1}^6{\mathit{Q}}_{bi}& \sum _{i=1}^6{\mathit{Q}}_{bi}{\tilde{\mathit{b}}}_i^T\\ \sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{bi}& \sum _{i=1}^6{\tilde{\mathit{p}}}_i{\mathit{Q}}_{bi}{\tilde{\mathit{b}}}_i^T\end{array}\right]$$

(25)

Then the kinetic equation of the platform can be expressed as:

$${\mathit{M}}_p{\ddot{\mathit{x}}}_p={\mathit{M}}_b{\ddot{\mathit{x}}}_b+{\mathit{J}}_p\mathit{F}$$

(26)

Substituting equations (21) into (26) and organizing it yields the workspace-based platform dynamics equation:

$${\mathit{M}}_p{\ddot{\mathit{q}}}_p+{\mathit{C}}_p{\dot{\mathit{q}}}_p+{\mathit{K}}_p{\mathit{q}}_p=\mathit{B}+{\mathit{J}}_p{\mathit{F}}_a$$

(27)

where ${\mathit{C}}_p=c{\mathit{J}}_p{\mathit{J}}_p^T{\dot{\mathit{q}}}_p,{\mathit{K}}_p=k{\mathit{J}}_p{\mathit{J}}_p^T{\dot{\mathit{q}}}_p,\mathit{B}={\mathit{M}}_b{\ddot{\mathit{q}}}_b+c{\mathit{J}}_p{\mathit{J}}_b^T{\dot{\mathit{q}}}_b$ $+k{\mathit{J}}_p{\mathit{J}}_b^T{\mathit{q}}_b$ .

Improved NPD controller design

Nonlinear functions

Compared with the traditional PD control, NPD control offers better rapidity and higher accuracy, capable of compensating for unknown or uncertain disturbances through its nonlinear nature. Confronting strongly coupled, strongly nonlinear, and uncertain systems, NPD control enhances the robustness of the system by introducing a nonlinear function in the error term, leading to a more favorable and stable control performance.

Before proceeding to the controller design, a central concept in NPD control, which is defined as a nonlinear function (fal function), shall be first introduced. When initially proposed by Han Jingqing, several different forms of nonlinear functions were developed to cope with different working conditions.^14,15 Subsequent studies and experiments have revealed that the fal function is the most widely used nonlinear function, defined as follows:^14,15

$$fal(\epsilon ,\alpha ,\delta )=\{\begin{array}{c}{|\epsilon |}^{\alpha }\mathrm{sgn}(e),|\epsilon |>\delta \\ {\displaystyle \frac{\epsilon }{{\delta }^{1-\alpha }},|\epsilon |\le \delta }\end{array}$$

(28)

where $\epsilon$ is the general independent variable, and specifically when using the fal function to describe the NPD control law, the independent variable is proportional error $e_{\mathrm{p}}$ or derivative error $e_{\mathrm{d}}$ . $\alpha$ is an exponential parameter that determines the rate of growth of the function value, and $\delta$ delimits the linear region in the nonlinear feedback.

Based on the different values of α, the nonlinear function is further analyzed. When $\alpha >1$ , the feedback is termed ‘smooth feedback’. During this period, the function has a ‘large gain for large error’, signifying increased sensitivity to changes in larger error values, as shown in Figure 4. Conversely, for $\alpha <1$ , the feedback is termed ‘non-smooth feedback’, characterized by the opposite traits of ‘smooth feedback’. Therefore, in the majority of existing research and applications, $\alpha <1$ is commonly employed.

Figure 4.

Images of the fal function corresponding to different values of $\alpha (\delta =0.001)$ .

For cases where $\alpha <1$ , the function possesses an exceptionally high growth rate in the region proximate to zero bilaterally. This characteristic raises concerns about substantial numerical fluctuations resulting from small sudden changes in error during practical applications. Such fluctuations can affect the control effect and overall system stability. After the introduction of the linear region, the growth rate of the function is then reduced in this part of the region, which can inhibit the possible effects due to the sudden change of the error.

$\delta$ determines the ‘degree of linearity’ of the entire fal function, and when its value is too large, the function will be linearized to a high degree and; thus, lose the advantages of nonlinearity. In practice, δ is usually set to a small value, ranging from 0.001−0.01, and it is also highlighted that a small value of δ may also cause system jitter in the literature. ¹⁶

Improved S-NPD method

Figure 5 shows the overall control scheme of the Stewart platform. Based on the fal function, a generalized expression for NPD control is introduced as^14,15

$$u=K_{\mathrm{p}}fal(e_{\mathrm{p}},{\alpha }_1,{\delta }_1)+K_{\mathrm{d}}fal(e_{\mathrm{d}},{\alpha }_2,{\delta }_2)$$

(29)

Figure 5.

NPD/S-NPD control system block diagram.

Here, the first term corresponds to the proportional (P) component of PD control, the second term corresponds to the differential (D) component of the PD control, $e_{\mathrm{p}}$ represents the error, and $e_{\mathrm{d}}$ is the differential component of the error. $K_{\mathrm{p}}$ and $K_{\mathrm{d}}$ are constant control gain.

In contrast to PD control, where the error term is a simple difference between the actual and desired values, NPD control adjusts this term to be a nonlinear function of the difference modifiable by two parameters: $\alpha$ and $\delta$ . Compared with the linear combinations of PD control, NPD control has more nonlinear combinations, providing the control algorithm with increased flexibility and adaptability.

While generalized NPD theory can initially meet control performance in most cases, specific considerations are essential for applications involving the six-degree-of-freedom Stewart platform, as addressed in this paper. The following points merit attention:

The direct control focus is on the six legs of the platform, each of which has different motion states and error information during coupling. Using identical parameters for controlling all six legs may result in information loss, compromising overall control effectiveness.

For the aforementioned challenge in (1), it is a solution idea to adjust the parameters of the six legs separately. However, this approach significantly escalates the difficulty and complexity of parameter adjustment, posing challenges to the practicality and simplicity of the NPD.

Realising real-time adjustment and online change of parameters is an idea that has been proposed since the research of PID, and the studies in the literatures 15, 16, 20 have proved its feasibility and improvement effect. Therefore, to further improve the performance of NPD control, it is essential to add this function.

To address the above challenge and deficiencies based on the theoretical formulas mentioned in Nonlinear functions section, coupled with the idea of adaptive gain, the parameters before the error term are also set in the form of nonlinear variable. Therefore, the parameters also change with the change of the error, making the nonlinear PD parameter evolve from the original ‘comprehensive optimization’ to the ‘point-by-point optimization’, to achieve the purpose of improving the control effect. This makes the nonlinear PD parameters evolve from the original ‘comprehensive optimization’ to ‘point-by-point optimization’ to achieve the purpose of improving the control effect.

To distinguish the two NPD controls, the constant control gain $K_{\mathrm{p}}$ and $K_{\mathrm{d}}$ are hereinafter turned variable, and a generalized PD parameter equation is further given as follows:

$$\{\begin{array}{c}{K}_{\mathrm{p}}=a_{\mathrm{p}}+b_{\mathrm{p}}f(e_{\mathrm{p}}(t))\\ K_{\mathrm{d}}=a_{\mathrm{d}}+b_{\mathrm{d}}g(e_{\mathrm{d}}(t))\end{array}$$

(30)

where f and g are two nonlinear functions, the specific form needs to be set by different situations so that the PD parameters can be made to vary with the error.

Theorem: the larger the error is, the larger the control force is needed, $K_{\mathrm{p}}$ which is positively correlated with the absolute value of the error; the faster the error changes and the more serious the oscillations are, the higher the stability is needed to reduce the amount of overshooting, $K_{\mathrm{d}}$ is negatively correlated with the speed of the error change.

Assumption: in real engineering applications, the disturbance and control errors suffered by the system when it reaches a steady state are always bounded.

The differential error, which is the derivative function of the proportional error reflects the trend of the latter, providing a basis for feed-forward correction of the parameters, thereby a new function S is defined as follows:

$$S(e)=\{\begin{array}{c}{s}_0,\mathrm{sgn}(e_{\mathrm{p}})=\mathrm{sgn}(e_{\mathrm{d}})\\ s_1,\mathrm{sgn}(e_{\mathrm{p}})\ne \mathrm{sgn}(e_{\mathrm{d}})\end{array}$$

(31)

where $s_0$ and $s_1$ are two constants that determine the range of variation of the proportional gain, which are set to 0.5 and 1.

Based on the proposed S-function formulation, the complete improved PD parameter formulation is given:

$$K_{\mathrm{p}}=\{\begin{array}{c}{K}_{\mathrm{p}0}+K_{\mathrm{p}0}S(e_{\mathrm{p}}){|e_{\mathrm{p}}|}^k\\ K_{\mathrm{pmax}},K_{\mathrm{p}}(e_{\mathrm{p}}(t))\ge K_{\mathrm{pmax}}\end{array}{K}_{\mathrm{d}}=\{\begin{array}{c}{K}_{\mathrm{d}0}-K_{\mathrm{d}1}{|e_{\mathrm{d}}|}^n\\ K_{\mathrm{dmin}},K_{\mathrm{d}}(e_{\mathrm{d}}(t))\le K_{\mathrm{dmin}}\end{array}$$

(32)

Taking $K_{\mathrm{p}0}$ as the base value of the proportional gain, the change of the error $e_{\mathrm{p}}$ is estimated by introducing the S function, and the specific amount of change is calculated by the exponential function to realize the real-time adjustment of the proportional gain.

However, when estimating the differential gain term, calculating the acceleration of the error change is necessary, leading to an increase in the computational and system complexity. Consequently, the S-function is not introduced in the differential term to avoid further complexities.

When the system is subjected to a large transient disturbance, or the sensor signal suffers a large disturbance, the error undergoes abrupt changes. During such instances, the PD parameters are calculated based on Equation (32), likewise change abruptly. To prevent the effects caused by this situation, the saturation zone of the PD parameter changes is set, and the calculated parameter values are limited to the set boundary values when they exceed the upper and lower limits, preventing excessive gain from generating additional jitter, which improves the system's stability.

With the proposed parameter formulas, it becomes apparent that smaller errors correspond to smaller proportional parameters and control gains, while large errors correspond to larger proportional parameters and control gains, thus theoretically realizing more precise control.

Numerical simulation and analysis

Simulation conditions

For the six-degree-of-freedom Stewart platform, the control gain is a six-dimensional vector, which is calculated separately for each leg.

$$\begin{array}{rl}& K_{\mathrm{p}}(e_{\mathrm{p}})=diag\{K_{\mathrm{p}1}(e_{\mathrm{p}1}),K_{\mathrm{p}2}(e_{\mathrm{p}2}),K_{\mathrm{p}3}(e_{\mathrm{p}3}),K_{\mathrm{p}4}(e_{\mathrm{p}4}),\\ & K_{\mathrm{p}5}(e_{\mathrm{p}5}),K_{\mathrm{p}6}(e_{\mathrm{p}6})\}\\ & K_{\mathrm{d}}(e_{\mathrm{d}})=diag\{K_{\mathrm{d}1}(e_{\mathrm{d}1}),K_{\mathrm{d}2}(e_{\mathrm{d}2}),K_{\mathrm{d}3}(e_{\mathrm{d}3}),K_{\mathrm{d}4}(e_{\mathrm{d}4}),\\ & K_{\mathrm{d}5}(e_{\mathrm{d}5}),K_{\mathrm{d}6}(e_{\mathrm{d}6})\}\end{array}$$

(33)

The algorithmic formulation of the final controller can be expressed as follows:

$$u(t)=K_{\mathrm{p}}(e_{\mathrm{p}})fa{l}_1(e_{\mathrm{p}})+K_{\mathrm{d}}(e_{\mathrm{d}})fa{l}_2(e_{\mathrm{d}})$$

(34)

The relevant parameters of the controller are set as follows:

$$\begin{array}{r}{\alpha }_1=0.5,{\alpha }_2=0.75,\delta =0.001,k=0.25,\\ n=0.75,K_{\mathrm{pmax}}=1.5K_{\mathrm{p}0},K_{\mathrm{dmin}}=0.5K_{\mathrm{d}0}\end{array}$$

The relevant parameters of the platform are set in the Table 1

Table 1.

Platform parameter.

Platform parts	Size/m	Mass/kg	Moment of inertia/kg m2
Upper platform	$r_P=0.2$	$m_P=10$	$I_P=[0.1045,0.1045,0.2083]$
Lower platform	$r_B=0.4$	$m_B=10$	$I_B=[0.4170,0.4170,0.8333]$
Upper leg	$l_u=0.3$	$m_u=0.4$	$I_u=[4.33,4.33,0.02]\times {10}^{-3}$
Lower leg	$l_d=0.3$	$m_d=0.6$	$I_d=[11.58,11.58,0.06]\times {10}^{-3}$
Load	$a_l=0.2$	$m_l=20$	$I_l=[0.1333,0.1333,0.1333]$

Results and analysis

The data obtained from the simulation are processed using the absolute amplitude integral formula proposed below:

$$IAA={\int }_0^{t_0}|A(t)|dt$$

(35)

where $t_0$ represents the simulation time and further the attenuation rate of the vibration can be calculated as follows:

$$\eta =1-{\displaystyle \frac{IA{A}_1}{IA{A}_2}}$$

(36)

Ensuring that the system reaches stability within each working cycle or simulation cycle is a prerequisite for various operating conditions experiments. Under this standard, further attention is paid to accuracy indicators, which are reflected in specific indicators, namely the degree of amplitude attenuation. Because the set operating conditions are all time-varying disturbances, including sinusoidal and random superimposed vibrations, rather than simple step vibrations, our existing research mainly focuses on improving accuracy.

Three simulation conditions were set up with the first two being sinusoidal vibration, mainly used to simulate vibration in ideal environments and verify the basic performance of the proposed controller. In literature, ²⁹ sine vibration was also used as input, and they designed a Stewart platform and control system for vibration isolation on ships. The design based on sinusoidal simulation verification also achieved good results. Condition 3 is composed of a series of vibrations with random amplitudes and frequencies within a certain range, which are superimposed to simulate situations that are closer to reality.

Case 1: Simulation of low-frequency disturbances, the vibration is sinusoidal with a frequency of 0.5 Hz and an amplitude of 0.2°. The control gain is set to $K_{\mathrm{p}0}=90,K_{\mathrm{d}0}=10$ . For PD control, $K_{\mathrm{p}}=75,K_{\mathrm{d}}=15$ . The results are shown in Figures 6–9.

Figure 6.

Upper platform vibration image of Case 1.

Figure 7.

Input force image for legs 1 and 2 of Case 1.

Figure 8.

Input force image for legs 3 and 6 of Case 1.

Figure 9.

Input force image for legs 4 and 5 of Case 1.

From Table 2 above it is observed that the improved S-NPD control algorithm on vibration suppression has a better performance, and under Case 1, compared with the original PD and NPD controls, S-NPD has a relative performance improvement of 17.2% and 16.4% respectively.

Table 2.

Summary of Case 1 regulation data.

Control algorithm	Maximum amplitude	Vibration attenuation rate (%)
PD	0.079°	78.19
NPD	0.052°	78.70
S-NPD	0.024°	91.62

The average vibration attenuation rate reflects the steady-state performance of the system throughout the simulation phase. In addition, we also want to focus on the transient performance of the system, so the maximum amplitude is chosen as another evaluation indicator. From the table data, it can be seen that S-NPD has a decrease of 0.055° in maximum amplitude compared to PD. This value itself provides a reference value, as there are certain differences between the specific values in the simulation environment and the real environment. Therefore, we focus on the reduced relative proportion, where 0.055° corresponds to a maximum amplitude attenuation of 69.6%, which is a significant performance improvement.

The driving force of the leg is an auxiliary evaluation indicator. If the required driving force is too high, it indicates that the actuator needs to face significant physical constraints, which may not be achievable in practical applications or significantly increase the complexity of the system. Due to the upper limit of the output driving force of the actuator, excessive driving force required can make the selection of the actuator more stringent. ²⁶

Table 3.

Summary of Case 1 IAE data.

Leg	Control algorithm	IAE
	PD	1.84 × 10−3
Legs 1 and 2	NPD	7.15 × 10−4
	S-NPD	2.16 × 10−4
	PD	2.65 × 10−3
Legs 3 and 6	NPD	7.12 × 10−4
	S-NPD	2.46 × 10−4
	PD	4.48 × 10−3
Legs 4 and 5	NPD	1.54 × 10−3
	S-NPD	5.98 × 10−4

A lower required driving force indicates a looser physical material constraint on the actuator, and at the same time it corresponds to lower energy consumption, which has contribution to practical engineering applications.

An analysis was conducted on the displacement errors of the three types of control corresponding to case 1, the results are shown in Figures 10–12, and the integral absolute error was calculated. Compared with NPD control, S-NPD control achieved a reduction of over 60% in the IAE of leg displacement as shown in Table 3, which is a performance improvement. However, the overall waveform and trend of leg error are similar to the changes in platform attitude. This is because the platform attitude is jointly controlled by six legs, and the vibration attenuation of the platform is our main performance indicator of concern. Therefore, in subsequent working conditions we did not separately analyze the displacement of the legs.

Figure 10.

Displacement error for leg 1 and 2.

Figure 11.

Displacement error for leg 3 and 6.

Figure 12.

Displacement error for leg 4 and 5.

Case 2: Simulation of intermediate frequency disturbances, the vibration is sinusoidal with a frequency of 1 Hz and an amplitude of 0.1°, The control gain is set to $K_{\mathrm{p}0}=90,K_{\mathrm{d}0}=20$ . For PD control, $K_{\mathrm{p}}=75,K_{\mathrm{d}}=15$ . The results are shown in Figures 13–16.

Figure 13.

Upper platform vibration image of Case 2.

Figure 14.

Input force image for legs 1 and 2 of Case 2.

Figure 15.

Input force image for legs 3 and 6 of Case 2.

Figure 16.

Input force image for legs 4 and 5 of Case 2.

Based on the dynamic analysis and simulation verification, the intrinsic frequency of the model is approximated to be 1 Hz. Therefore, the regulation of working Case 2 is the most challenging and the regulation effect is relatively the worst as shown in Table 4. Simultaneously, S-NPD is compared with the original PD control and NPD control, and S-NPD has a relative performance improvement of 64.0% and 31.4% respectively.

Table 4.

Summary of Case 2 regulation data.

Control algorithm	Maximum amplitude	Vibration attenuation rate (%)
PD	0.054°	49.29
NPD	0.041°	61.51
S-NPD	0.021°	80.83

A comprehensive analysis of the vibration pattern of the upper platform shows that the input disturbance is a sinusoidal vibration with favorable periodicity, and the vibration of the upper platform reflects its corresponding attitude. Analyzing the vibration images of the upper platform under the control of the three different algorithms, it is observed that under the low-frequency disturbance (Cases 1 and 2), the overall vibration curves of the PD control are smoother but the peaks of the vibration are higher, in the startup phase, both the NPD control and the S-NPD control correspond to the faster change of the driving force. Thus, resulting in faster amplitude change. However, after reaching the stabilizing stage, the change in the upper platform's attitude under the NPD control tends to be smoother and the stability of the modulated platform is more powerful.

In practical settings, vibrations, especially those of high frequency, often manifest as a composite of various frequencies and amplitudes. To further verify the efficacy of the designed controller in addressing the intricacies of environmental complexity and its ability to suppress vibrations, the following conditions are proposed for continued simulation tests:

Case 3: Simulation of random disturbances, the vibration is superimposed by random sinusoidal vibrations of which the frequency ranges from 5−15 Hz and the amplitude ranges from 0.003°−0.015°. The control gain is set to $K_{\mathrm{p}0}=30,K_{\mathrm{d}0}=5$ . For PD control, $K_{\mathrm{p}}=75,K_{\mathrm{d}}=15$ . The results are shown in Figures 17–20.

Figure 17.

Upper platform vibration image of Case 3.

Figure 18.

Input force image for legs 1 and 2 of Case 3.

Figure 19.

Input force image for legs 3 and 6 of Case 3.

Figure 20.

Input force image for legs 4 and 5 of Case 3.

In simulated real-world conditions featuring random vibration, PD control experiences a further decline in performance, resulting in large vibration peaks. During such instances, the advantages of NPD control become even more pronounced. The improved S-NPD control further optimizes this effect, producing smoother vibration waveforms and enhancing overall regulation stability as shown in Table 5.

Table 5.

Summary of Case 3 regulation data.

Control algorithm	Maximum amplitude	Vibration attenuation rate (%)
PD	0.03782°	66.84
NPD	0.01776°	84.1
S-NPD	0.01633°	92.19

The driving forces of NPD and S-NPD control are similar in the overall change trend, with a close frequency and period. However, the amplitude of the driving force of S-NPD is smaller due to the inclusion of the link of online rectification coefficients, which makes the driving force change in real time according to the regulation needs and improves the control efficiency.

When the vibration frequency is low (Cases 1 and 2), the frequency of change of the driving force of the PD control is low, which makes the driving force lack the capability to make changes in advancing the disturbance, leading to the slower change of the position of the controlled object, and the error generated by the disturbance cannot be regulated on time, restricting the control effect. The composition of the error term of the PD and NPD control is different, and the original error presents a relatively smooth sinusoidal change. In contrast, the NPD control's error term, processed nonlinearly, has a higher change frequency and initial amplitude, providing a more agile response to disturbances. However, once the entire control system stabilizes in a steady state, the driving force of NPD control becomes even lower. In cases where disturbance frequencies vary and represent a complex case of multiple vibrations superimposed on each other (Case 3), the driving forces corresponding to the three controls share similar overall waveforms, while the improved S-NPD control corresponds to the lowest driving force with the peaks, which is in agreement with the conclusions drawn from the simple case tests.

Conclusion

This paper proposes a S-NPD control algorithm designed to facilitate real-time parameter adjustments and improve controller robustness. The proposed S-NPD controller addresses the challenges of information omission and parameter immutability, inheriting the advantages of classical PD control while using the nonlinear advantages of NPD control.

Using the Newton-Euler method, a dynamics model is developed for the Stewart platform, and a series of simulation tests are conducted. The simulations cover cases with low frequencies (0.5 Hz), intrinsic frequencies (1 Hz), and random sinusoidal vibrations (frequencies ranging from 5−15 Hz). Comparative analyses with PD and NPD control reveal that S-NPD control consistently achieves better control effectiveness with lower input driving force in each simulation condition. This becomes particularly advantageous when considering actuator constraints.

Future endeavors can focus on further optimizing switching conditions, and conducting tests that account for flexible structure.