Accepted: June 13, 2024 | Published Online: June 15, 2024

Citation: Baccouch M, Sun J (2024) Proportional-Integral-Derivative Controller for an n-Link Revolute Robot Manipulator. Int J Robot Eng 7:038

Copyright: © 2024 Baccouch M, et al. This This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

This paper aims to model and guide the motion of an n -link revolute robot arm manipulator by implementing an efficient proportional-integral-derivative (PID) control strategy. The n -link robot arm manipulator we examine shares similarities with a pendulum system comprising n arms and n masses. This system features a stationary point around which n distinct links, each carrying a mass at its end, are interconnected. We first derive the equations of motion using Lagrangian formulation. We then implement a PID controller for the robot manipulator, enabling it to attain designated target positions as required. These equations are represented by a second-order system of nonlinear ordinary differential equations. Due to the absence of closed-form solutions for the equations of motion, we employ the classical fourth-order Runge-Kutta method to approximate the solution of the initial-value problem. Due to the inherent nonlinear nature, achieving precise control over the motion of the n -link robot manipulator at user-defined positions poses a challenging endeavor. In light of this, our primary emphasis centers on controlling the robot manipulator to attain the desired position using the PID controller. Numerous computer simulations are conducted to validate the controller’s performance. Notably, a PID controller is presented as part of the simulations, illustrating how we can achieve balance for the n links with n = 3 − 6 on a moving robot at various angles.

AMS subject classifications (2020): 65L05, 65L06, 70E60, 70B15, 70Q05, 93C85.

Keywords

n -Link robotic arm, Robot Manipulator, Dynamic equations, Lagrangian, PID controller

Introduction

The purpose of this paper is to design a nonlinear proportional-integral-derivative (PID) controller for trajectory tracking of a manipulator robot. Robotic arms play an essential role in various sectors, including manufacturing, transportation, and healthcare. They contribute significantly to enhancing human livelihoods, enabling their endeavors and innovations, ultimately leading to improved quality of life. An overview of commercially available robotic manipulators, sensors, and controllers can be found in [1]. The field of robotic manipulator control is a well-established and promising domain that encompasses research, development, and manufacturing. Industrial robots primarily function as positioning and handling devices. Hence, a valuable robot possesses the capability to control its motion along with the interactive forces and torques occurring between the robot and its surroundings.

Robot Manipulators are composed of links connected by joints into a kinematic chain. These joints predominantly consist of rotary components capable of achieving a diverse range of motions, thereby imparting flexibility and maneuvering capabilities to the robot. The n -link robot arm manipulator closely resembles an n -pendulum system with n arms and n masses, where there exists a fixed point about which n different links with a mass at the end of each are connected upon. Each adjacent mass and link combination forms its own simple pendulum such that the system contains n total simple pendulums, each conjoined with one another to form an n -pendulum system. When released from an arbitrary position without control, each pendulum is able to freely oscillate on the xy -plane. While the system itself is simple dynamically, the behavior exhibited by it is complex and nonlinear. In order to understand this behavior and define the motion of the system, we design a nonlinear PID controller for trajectory tracking of a manipulator robot.

Achieving precision and speed in the movements of a robot manipulator is critical in its efficiency. In this paper, we propose a proportional-integral-derivative (PID) control to achieve this. To define clearly designated parameters and to track the position of the system within those parameters we use a PID controller. The PID controller has three parts: a proportional (P) controller, an integral (I) controller, and a derivative (D) controller. This allows for the optimization of speed, stability, and precision. Subsequently, PID controllers are known to be reliable and efficient, considered a building block in industrial technology and the standard control in industrial robotics as well as research [2].

Other control systems have been examined as well for the n -link manipulator. An overview of various control theory concepts that are used in the control of robots can be found in [1]. In [3], the authors presented a numerical solution of the second-order robot arm control problem using the Runge-Kutta-Butcher algorithm. Research has also formulated using sets of stabilizing decentralized PID controls with Kharitonov’s theorem for n -link manipulators to track positions and was successfully simulated on a two-link manipulator [4]. Recently, in [5], we designed a robust, fast, and practical PID controller for the classical double pendulum system. We first derived the equations of motion for the two-link robot manipulator using the Lagrangian approach. We used the classical fourth-order Runge-Kutta method to approximate the solution of the nonlinear system of second-order ordinary differential equations. We focused mainly on control of the robot manipulator to get the desired position using the PID control approach. In this paper, we extend the approach to an n -link robot manipulator for any n . More recently, an Adaptive Observer Based Neural Control for flexible-joint manipulators has been researched, utilizing the backstepping method with Lyapunov stability theory [6]. A Sliding Mode Control employed with the super-twisting method has also been successfully applied to robot manipulators, demonstrating fast error convergence and angular velocity approximation for each link [7]. Others have used the Lyapunov-based Control Scheme for linear n -link manipulators mounted on mobile sliders to navigate links [8].

In this paper, we propose a PID control for an n -link robot manipulator to guide its motion from an initial position to a target position on a plane in Figure 1. We first derive the total potential and kinetic energies for the n -link system. Then, we explicitly find the so-called Lagrangian. The Euler-Lagrangian equation [9] is then used to define the torques for each link. Our proposed PID controller uses the classical fourth order Runge-Kutta method [10] to approximate the angles and torques of the links relative to the x -axis from a system of nonlinear differential equations and initial conditions. The numerical approximations are performed by Java code. The system of differential equations our method relies on is based on the equations of motion for the n -link robot manipulator. These equations of motion are derived through Lagrangian dynamics, specifically from the Euler-Lagrange equation that is based on the total kinetic and potential energies of the n -link robot manipulator system. Through this control derivation, we are able to obtain the most precise and efficient guidance of the manipulator’s motion to its designated target positions. This is demonstrated through subsequent simulations of the control that we perform.

The paper is organized as follows. Section 2 presents the derivation of the equations of motion for the n -link robotic manipulator using the Lagrangian approach. Furthermore, the procedure for obtaining numerical solutions using the classical fourth-order Runge-Kutta method is detailed. The design of the PID controller for the n -pendulum system is expounded upon in Section 3 . Section 4 provides a series of numerical examples intended to evaluate the performance of the PID controller. Concluding remarks are offered in Section 5 .

Dynamics Equations

The robot manipulator is a device composed of multiple conjoined arms, joints, and links integrated with a control mechanism designed to manipulate objects and perform physical tasks without manual human intervention. Robot manipulators are consequently vital to the world of industry and everyday life in assembly, automation, manufacturing, logistics, and more.

A typical application involving an industrial manipulator with n revolute links is shown in Example 6.1. This robot can be viewed as a pendulum made of n masses m ₁ , m ₂ , . . . , m _n and n rods of lengths ℓ ₁ , ℓ ₂ , . . . , ℓ _n .

For our model problem, we consider a Cartesian coordinate system with origin placed at the base of the robot.

The Equations of Motion

We may understand the n -link robot manipulator as an n -link pendulum system on a two-dimensional plane. This system consists of n weightless bars, having lengths of ℓ ₁ , ℓ ₂ , . . . , ℓ _n . There exists a mass at the end of each bar such that there are n masses in total connected by the bars. Denote these n masses as m ₁ , m ₂ , . . . , m _n . Each bar may rotate freely, where the first bar rotates about the origin and each subsequent bar rotates about the endpoint of its preceding bar, producing n degrees of freedom. The angles for each of these degrees of freedom produced by the bar relative to the x-axis may be denoted as θ ₁ , θ ₂ , . . . , θ _n . Such a system exhibits chaotic motion that may be studied through Lagrangian dynamics. In order to derive the equations of motion through Lagrangian dynamics, we must first solve for the kinetic energy and the potential energy of the system.

Let ( x _i , y _i ) denote the position of the mass m _i for each i = 1,2,…, n . The position ( x _i , y _i ) of m _i at time t is given by the equations:

$$x_i={\displaystyle \sum _{j=1}^i{\mathcal{l}}_j\cos {\theta }_j,}{y}_i={\displaystyle \sum _{j=1}^i{\mathcal{l}}_j\sin {\theta }_j,}i=1,2,\mathrm{...},n.\text{ (4}\text{.1)}$$

The velocity of mass m _i , denoted as v _i , is given by

$$v_i=\sqrt{{\dot{x}}_i^2+{\dot{y}}_i^2},i=1,2,\mathrm{...},n,\text{ (4}\text{.2)}$$

where ${\dot{x}}_i$ and ${\dot{y}}_i$ are, respectively, the derivative of x _i and y _i with respect to t . Since

$${\dot{x}}_i=-{\displaystyle \sum _{j=1}^i{\mathcal{l}}_j{\dot{\theta }}_j\sin {\theta }_j,}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}{\dot{y}}_i={\displaystyle \sum _{j=1}^i{\mathcal{l}}_j{\dot{\theta }}_j\cos {\theta }_j,}i=1,2,\mathrm{...},n,$$

we have

$$v_i=\sqrt{{\left(-{\displaystyle \sum _{j=1}^i{\mathcal{l}}_j{\dot{\theta }}_j\sin {\theta }_j}\right)}^2+{\left({\displaystyle \sum _{j=1}^i{\mathcal{l}}_j{\dot{\theta }}_j\cos {\theta }_j}\right)}^2},i=1,2,\mathrm{...},n.\text{ (4}\text{.3)}$$

Consequently, we can express the total kinetic energy of the system, KE, as

$$\begin{array}{l}KE={\displaystyle \sum _{i=1}^n\frac{1}{2}{m}_i{v}_i^2}\text{ (4}\text{.4)}\\ ={\displaystyle \sum _{i=1}^n\frac{1}{2}{m}_i\left({\left(-{\displaystyle \sum _{j=1}^i{\mathcal{l}}_j{\dot{\theta }}_j\sin {\theta }_j}\right)}^2+{\left({\displaystyle \sum _{j=1}^i{\mathcal{l}}_j{\dot{\theta }}_j\cos {\theta }_j}\right)}^2\right)}\\ ={\displaystyle \sum _{i=1}^n\frac{1}{2}{m}_i}\left({\displaystyle \sum _{j=1}^i{\displaystyle \sum _{k=1}^i{\mathcal{l}}_j{\mathcal{l}}_k{\dot{\theta }}_j{\dot{\theta }}_k\sin {\theta }_j\sin {\theta }_k}+}{\displaystyle \sum _{j=1}^i{\displaystyle \sum _{k=1}^i{\mathcal{l}}_j{\mathcal{l}}_k{\dot{\theta }}_j{\dot{\theta }}_k\cos {\theta }_j\cos {\theta }_k}}\right)\\ =\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{\displaystyle \sum _{k=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_k{\dot{\theta }}_j{\dot{\theta }}_k\left(\sin {\theta }_j\sin {\theta }_k+\cos {\theta }_j\cos {\theta }_k\right)}}}\\ =\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{\displaystyle \sum _{k=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_k{\dot{\theta }}_j{\dot{\theta }}_k\cos \left({\theta }_j-{\theta }_k\right)}}.}\end{array}$$

We may also find the total potential energy of the system by summing the gravitational potential energies of each mass m _i to get

$$PE={\displaystyle \sum _{i=1}^n{m}_{i}g{y}_i={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{m}_{i}g{\mathcal{l}}_j\sin {\theta }_j.}}}\text{ (4}\text{.5)}$$

Next, we define the Lagrangian function $\mathcal{L}$ , also called the Lagrangian, which is quantity that characterizes the state of a physical system. The Lagrangian function is just the kinetic energy (energy of motion) minus the potential energy (energy of position) i.e. ,

$$\mathcal{L}=KE-PE=\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{\displaystyle \sum _{k=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_k{\dot{\theta }}_j{\dot{\theta }}_k\cos ({\theta }_j-{\theta }_k)-{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{m}_{i}g{\mathcal{l}}_j\sin {\theta }_j.}}}}}\text{ (4}\text{.6)}$$

With these definitions, the Euler-Lagrange equations (or Lagrange’s equations) of motion for a conservative system are given by

$$\mathcal{L}=KE-PE=\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{\displaystyle \sum _{k=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_k{\dot{\theta }}_j{\dot{\theta }}_k\cos ({\theta }_j-{\theta }_k)-{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^i{m}_{i}g{\mathcal{l}}_j\sin {\theta }_j.}}}}}\text{ (4}\text{.6)}$$

where τ _q is the torque applied to the q ^th link. The derivations for the Lagrangian and the Euler-Lagrange equation may be found in [11].

Using (4.6), we have, for any q = 1,2,…, n ,

$$\begin{array}{l}\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_q}={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\dot{\theta }}_j\cos ({\theta }_q-{\theta }_j)}}\\ \frac{\partial \mathcal{L}}{\partial {\theta }_q}=-{\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\dot{\theta }}_j{\dot{\theta }}_q\sin ({\theta }_q-{\theta }_j)}}-{\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q}\\ \frac{d}{dt}\left[\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_q}\right]={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q\left[{\ddot{\theta }}_j\cos ({\theta }_q-{\theta }_j)-{\dot{\theta }}_j({\dot{\theta }}_q-{\dot{\theta }}_j)\sin ({\theta }_q-{\theta }_j)\right].}}\end{array}$$

Substituting these relations into (4.7), we get the equation of motion

$$\begin{array}{l}{\tau }_q=\frac{d}{dt}\left[\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_q}\right]-\frac{\partial \mathcal{L}}{\partial {\theta }_q}\\ ={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q\left[{\ddot{\theta }}_j\cos ({\theta }_q-{\theta }_j)-{\dot{\theta }}_j({\dot{\theta }}_q-{\dot{\theta }}_j)\sin ({\theta }_q-{\theta }_j)\right]}}\\ +{\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\dot{\theta }}_j{\dot{\theta }}_q\sin ({\theta }_q-{\theta }_j)+{\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q}}}\\ ={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q\left[{\ddot{\theta }}_j\cos ({\theta }_q-{\theta }_j)-{\dot{\theta }}_j({\dot{\theta }}_q-{\dot{\theta }}_j)\sin ({\theta }_q-{\theta }_j)+{\dot{\theta }}_j{\dot{\theta }}_q\sin ({\theta }_q-{\theta }_j)\right]}}\\ +{\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q}\\ ={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q\left[{\ddot{\theta }}_j\cos ({\theta }_q-{\theta }_j)+{\dot{\theta }}_j^2\sin ({\theta }_q-{\theta }_j)\right]+}}{\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q},\end{array}$$

for q = 1,2,…,n. Thus, we get the following n nonlinear second-order system of ordinary differential equations

$${\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\ddot{\theta }}_j\cos ({\theta }_q-{\theta }_j)+}}{\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\dot{\theta }}_j^2\sin ({\theta }_q-{\theta }_j)+}}{\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q={\tau }_q.}\text{ (4}\text{.8)}$$

Define

$$\theta =\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\\ ⋮\\ {\theta }_n\end{array}\right],\dot{\theta }=\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\\ ⋮\\ {\dot{\theta }}_n\end{array}\right],\ddot{\theta }=\left[\begin{array}{c}{\ddot{\theta }}_1\\ {\ddot{\theta }}_2\\ ⋮\\ {\ddot{\theta }}_n\end{array}\right],F=\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\\ ⋮\\ {\tau }_n\end{array}\right],c(\theta ,\dot{\theta })=\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_n\end{array}\right],G(\theta )=\left[\begin{array}{c}{g}_1\\ g_2\\ ⋮\\ g_n\end{array}\right],$$

$$M(\theta )=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ a_{21}& a_{22}& a_{23}& \cdots & a_{2n}\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ a_{n1}& a_{n2}& a_{n3}& \cdots & a_{nn}\end{array}\right],$$

where

$$\begin{array}{l}{a}_{qk}={\displaystyle \sum _{i=\min (k,q)}^n{m}_i{\mathcal{l}}_k{\mathcal{l}}_q\cos ({\theta }_q-{\theta }_k),}q,k=1,2,\mathrm{...},n,\\ c_q={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\dot{\theta }}_j^2\sin ({\theta }_q-{\theta }_j),}}q=1,2,\mathrm{...},n,\\ g_q={\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q,q=1,2,\mathrm{...},n.}\end{array}$$

Then a more compact formulation of (4.8) is given by

$$M(\theta )\ddot{\theta }+c(\theta ,\dot{\theta })+G(\theta )=F,\text{ (4}\text{.9)}$$

where M is the n × n mass (or inertia) matrix, θ is the n × 1 trajectory vector (the generalized coordinates), $\dot{\theta }$ and $\ddot{\theta }$ are the first and second derivative of θ , c( θ, $\dot{\theta }$ ) is the n × 1 Centrifugal/Coriolis force containing velocity-dependent torques, G ( θ ) is the n × 1 gravitational torque, and F contains the input joint torques.

Since the kinetic energy is positive, vanishing only when the generalized velocity equals zero, the inertia matrix M ( θ ) is also positive definite. Thus, the matrix M is invertible. Multiplying (4.9) from the left by M ^{− 1} , we obtain the normal form of the dynamics equations

$$\ddot{\theta }=f(t,\theta ,\dot{\theta }),\text{ (4}\text{.10a)}$$

where

$$f(t,\theta ,\dot{\theta })=M^{-1}(\theta )F-M^{-1}(\theta )c(\theta ,\dot{\theta })-M^{-1}(\theta )G(\theta ).\text{ (4}\text{.10b)}$$

A Numerical Scheme

In this subsection, we use the fourth-order Runge-Kutta (RK4) method to approximate the solution of the following initial-value problem (IVP)

$$\ddot{\theta }=f(t,\theta ,\dot{\theta }),t\in \left[0,T\right],\text{ (4}\text{.11a)}$$

$$\theta (0)=\alpha ,\dot{\theta }(0)=\beta ,\text{ (4}\text{.11b)}$$

where the initial conditions $\alpha ={\left[{\alpha }_1,{\alpha }_2,\mathrm{...},{\alpha }_n\right]}^t$ and $\beta ={\left[{\beta }_1,{\beta }_2,\mathrm{...},{\beta }_n\right]}^t$ are given.

To apply the RK4 method, we first convert the second-order system of equations into a first-order system of equations. To do this, we introduce the new variables

$$\begin{array}{l}{u}_i={\theta }_i,i=1,2,\mathrm{...},n,\\ u_{n+i}=\dot{\theta }{}_i,i=1,2,\mathrm{...},n.\end{array}$$

Since ${\dot{u}}_i={\dot{\theta }}_i=u_{n+1}$ and ${\dot{u}}_{n+i}={\ddot{\theta }}_i=f_i$ , i = 1,2,…, n , we get the first-order system of equations

$$\dot{u}=s(t,u),t\in \left[0,T\right],\text{ (4}\text{.12a)}$$

$$u(0)=u^{(0)},\text{ (4}\text{.12b)}$$

where

$$\begin{array}{l}u={\left[u_1,u_2,\mathrm{...},u_n,u_{n+1},\mathrm{...},u_{2n}\right]}^t,\\ s={\left[s_1,s_2,\mathrm{...},s_n,s_{n+1},\mathrm{...},s_{2n}\right]}^t\\ ={\left[u_{n+1},u_{n+2},\mathrm{...},u_{2n},f_1(t,u{}_1,u{}_2,\mathrm{...},u{}_n,u{}_{n+1},\mathrm{...},u{}_{2n}),\mathrm{...},f_n(t,u{}_1,u{}_2,\mathrm{...},u{}_n,u{}_{n+1},\mathrm{...},u{}_{2n})\right]}^t,\\ u^{(0)}=\left[{\alpha }^t,{\beta }^t\right]={\left[{\alpha }_1,{\alpha }_2,\mathrm{...},{\alpha }_n,{\beta }_1,{\beta }_2,\mathrm{...},{\beta }_n\right]}^t.\end{array}$$

Next, we use the fourth-order Runge-Kutta method [12] to approximate the solution u to the IVP (4.12). We first discretize the computational interval [0, T ] by non-overlapping intervals I _i = ( t _i-1 , t _i ), i = 1,2,…, m such that 0 = t ₀ < t ₁ < ∙ ∙ ∙ < t _m = T . For simplicity we assume that the partition is uniform. Let us define the length of I _i as $h_i=x_i-x_{i-1}=\frac{T}{m}$ . Then for each i = 1,2,…, m we compute k ₁ , k ₂ , k ₃ , k ₄ using

$$\begin{array}{l}{k}_1=s(t{}_i,u^{(i)}),\\ k_2=s(t{}_i+\frac{h}{2},u^{(i)}+\frac{h}{2}{k}_1),\\ k_3=s(t{}_i+\frac{h}{2},u^{(i)}+\frac{h}{2}{k}_2),\\ k_4=s(t{}_i+\frac{h}{2},u^{(i)}+h{k}_3),\end{array}$$

And finally the RK4 approximation of u(ⁱ⁺¹) at t_i+1 is given by

$$u^{(i+1)}=u^{(i)}+\frac{h}{6}(k_1+2k_2+2k_3+k_4),i=\mathrm{1,2,....},m,$$

For some arbitrary step-size h > 0. The initial condition is given by

$$u^{(0)}={\left[u_1,u_2,\mathrm{...},u_n,u_{n+1},\mathrm{...},u_{2n}\right]}^t={\left[{\theta }_1(0),{\theta }_2(0),\mathrm{...},{\theta }_n(0),{\dot{\theta }}_1(0),\mathrm{...},{\dot{\theta }}_n(0)\right]}^t.$$

PID Controller

For the robot manipulator to reach target positions, we require the stabilization of its motion and the ability to fix it at a particular position. In order to achieve this, we must implement a control method. Our control method utilizes a computed torque technique based on our equations of motion to hold and move each link of the n -link robot manipulator at an angle θ _i with respect to the x-axis, such that it may obtain a specified position.

We utilize the proportional-integral-derivative (PID) control in this paper. The PID control is one of the most effective and popular control systems. There are three parts to the control system: the Proportional control (P), the Integral control (I), and the Derivative control (D). The PID control works through utilizing a continuously calculated error value V _E = V _SP −V _PV . where V_sp is a desired setpoint and V_pv is measured process variable. The Proportional, Integral, and Derivative controls may now be expressed with constants K _P , K _I , K _D , as K _P V _E $K_I{\displaystyle \int V_{E}dt}$, and $K_D\frac{d{V}_E}{dt}$ , respectively.

In this paper, we define the setpoint as θ _f and the process variable as θ ( t ) such that our error function e ( t ) is e ( t ) = θ _f − θ ( t ). The three components are each responsible for their own purposes in the overall control. The Proportional control is the main drive in the system and acts proportional to the error value V _E . The Integral control is meant to account for residual error produced by the Proportional control over time, integrating it to eliminate it. The Derivative control estimates future error values based on its rate of change, reducing overshoot and ringing.

To formulate our PID control approach, we consider the dynamics of n -link robot manipulator described by the nonlinear equation

$$M(\theta )\ddot{\theta }+c(\theta ,\dot{\theta })+G(\theta )=\tau ,\text{ (5}\text{.1)}$$

where M( θ ) = ( m _qk ) _{1≤ q,k ≤ n} is an n × n symmetric positive definite inertia matrix, $c(\theta ,\dot{\theta })={[c_1,c_2,\mathrm{...},c_n]}^t$ is an n × 1 Coriolis and centrifugal vector, G ( θ ) = [ g ₁ , g ₂ ,…, g _n ] ^t is an n × 1 gravity vector of the manipulator, θ = [ θ ₁ , θ ₂ ,…, θ _n ] ^t is the n × 1 vector representing joint angular positions, and τ = [τ ₁ , τ ₂ , …, τ _n ] ^t is the n × 1 vector of applied joint torques. The entries of M ( θ ), $c(\theta ,\dot{\theta })$ , and G ( θ ) are, respectively,

$$\begin{array}{l}{m}_{qk}={\displaystyle \sum _{i=\min (k,q)}^n{m}_i{\mathcal{l}}_k{\mathcal{l}}_q\cos ({\theta }_q-{\theta }_k),q,k=1,2,\mathrm{...},n,}\\ c_q={\displaystyle \sum _{i=q}^n{\displaystyle \sum _{j=1}^i{m}_i{\mathcal{l}}_j{\mathcal{l}}_q{\dot{\theta }}_j^2\sin ({\theta }_q-{\theta }_j),}}q=1,2,\mathrm{...},n,\\ g_q={\displaystyle \sum _{i=q}^n{m}_{i}g{\mathcal{l}}_q\cos {\theta }_q,}q=1,2,\mathrm{...},n.\end{array}$$

Solving for $\ddot{\theta }$, we obtain

$$\ddot{\theta }=-M^{-1}(\theta )\left[c(\theta ,\dot{\theta })+G(\theta )\right]+\widehat{\tau },\text{ (5}\text{.2)}$$

where $\widehat{\tau }=M^{-1}(\theta )\tau ={\left[{\widehat{\tau }}_1,{\widehat{\tau }}_2,\mathrm{...},{\widehat{\tau }}_n\right]}^t$ . We define the position tracking error of the q ^th joint as

$$e_q(t)={\theta }_{qf}-{\theta }_q(t),\text{ (5}\text{.3)}$$

For all q = 1,2,…, n where ${\theta }_{q_f}$ denotes the q ^th joint’s desired constant position and θ _q is the actual joint position.

A PID controller for the q ^th joint is described by

$${\widehat{\tau }}_q(t)=K_{P_q}{e}_q(t)+K_{D_q}{\dot{e}}_q(t)+K_{I_q}{\displaystyle {\int }_0^t{e}_q(s)ds,}\text{ (5}\text{.4)}$$

where $K_{P_q}$ , $K_{D_q}$ and $K_{I_q}$ are, respectively, the proportional, derivative, and integral gains of the q ^th joint controller. The control problem is to provide a complete solution to the constant gain stabilizing control parameters $K_{P_q}$ , $K_{D_q}$ and $K_{I_q}$ such that the position error e _q ( t ) reduces to zero with time, i.e. , $\underset{t\to \infty }{\lim }{e}_q(t)=0$.

The input to the system (3.2) with the PID controller is $\widehat{\tau }={\left[{\widehat{\tau }}_1,{\widehat{\tau }}_2,\mathrm{...},{\widehat{\tau }}_n\right]}^t$, where for each q = 1,2,...,n,

$$\begin{array}{l}{\widehat{\tau }}_q(t)=K_{P_q}{e}_q(t)+K_{D_q}{\dot{e}}_q(t)+K_{I_q}{\displaystyle {\int }_0^t{e}_q(s)ds}\text{ (5}\text{.5)}\\ =K_{P_q}({\theta }_{q_f}-{\theta }_q)-K_{D_q}{\dot{\theta }}_q+K_{I_q}{\displaystyle {\int }_0^t({\theta }_{q_f}-{\theta }_q)ds,}\\ \end{array}$$

where we used the fact that ${\dot{e}}_q=-{\dot{\theta }}_q$ since ${\theta }_{q_f}$ is a constant denoting the desired target angle.

To implement the PID controller, we introduce new variables x ₁ , x ₂ ,…, x _n , as follows

$$x_q(t)={\displaystyle {\int }_0^t{e}_q(s)ds=}{\displaystyle {\int }_0^t({\theta }_{q_f}-{\theta }_q)ds,}q=1,2,\mathrm{...},n,$$

So that

$${\dot{x}}_q(t)=e_q(t)={\theta }_{q_f}-{\theta }_q,q=1,2,\mathrm{...},n.$$

Thus, ${\widehat{\tau }}_q(t)$ can be written as

$${\widehat{\tau }}_q(t)=K_{P_q}({\theta }_{q_f}-{\theta }_q)-K_{D_q}{\dot{\theta }}_q+K_{I_q{x}_q}.\text{ (5}\text{.6)}$$

With this notation, we rewrite (5.2) as

$$\ddot{\theta }=z(t,x,\theta ,\dot{\theta }),\text{ (5}\text{.7)}$$

Where $z(t,x,\theta ,\dot{\theta })=-M^{-1}(\theta )\left[c(\theta ,\dot{\theta })+G(\theta )\right]+\widehat{\tau }.$ We note that

$$\widehat{\tau }=\left[\begin{array}{c}{\widehat{\tau }}_1\\ {\widehat{\tau }}_2\\ ⋮\\ {\widehat{\tau }}_n\end{array}\right]=\left[\begin{array}{c}{K}_{P_1}({\theta }_{1f}-{\theta }_1)-K_{D_1}{\dot{\theta }}_1+K_{I_1}x{}_1\\ K_{P_2}({\theta }_{2f}-{\theta }_2)-K_{D_2}{\dot{\theta }}_2+K_{I_2}x{}_2\\ ⋮\\ K_{P_n}({\theta }_{nf}-{\theta }_{2n})-K_{D_n}{\dot{\theta }}_n+K_{I_n}x{}_n\end{array}\right]=K_P({\theta }_f-\theta )-K_D\dot{\theta }+K_{I}x,$$

Where

$$K_s=\left[\begin{array}{ccccc}{K}_{s_1}& 0& 0& \dots & 0\\ 0& K_{s_2}& 0& \dots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& 0& 0& \dots & K_{s_n}\end{array}\right],{\theta }_f=\left[\begin{array}{c}{\theta }_{1f}\\ {\theta }_{2f}\\ ⋮\\ {\theta }_{nf}\end{array}\right],x=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right],$$

where S stands for either P or D or I .

To numerically implement the PID controller, we convert the second-order system of differential equations (5.7) into a first-order system of differential equations by introducing the following variables:

$$u_i=x_i,u_{n+i}={\theta }_i,u_{2n+i}={\dot{\theta }}_i,i=1,2,\mathrm{...},n.$$

Differentiating them with respect to t , we get

$$\begin{array}{l}{\dot{u}}_i={\dot{x}}_i={\theta }_{if}-{\theta }_i={\theta }_{if}-u_{n+i},\\ {\dot{u}}_{n+i}={\dot{\theta }}_i=u_{2n+i},\\ {\dot{u}}_{2n+i}={\ddot{\theta }}_i=z_i(t,u_1,u_2,\mathrm{...},u_{3n}),i=1,2,\mathrm{...},n,\end{array}$$

where z _i are the entries of z i.e., we used the notation $z(t,x,\theta ,\dot{\theta })={\left[z_1,z_2,\mathrm{...},z_n\right]}^t$ .

Define

$$u={\left[u_1,\mathrm{...},u_n,u_{n+1},\mathrm{...},u_{2n},u_{2n+1},\mathrm{...},u_{3n}\right]}^t$$.

Then we obtain

$$\dot{u}=g(t,u),u(0)=u^{(0)},\text{ (5}\text{.8)}$$

where

$$g={\left[{\theta }_{1f}-u_{n+1},\mathrm{...},{\theta }_{nf}-u_{2n},u_{2n+1},\mathrm{...},u_{3n},z_1(t,u_1,u_2,\mathrm{...},u_{3n}),\mathrm{...},z_n(t,u_1,u_2,\mathrm{...},u_{3n})\right]}^t.$$

Finally, the fourth-order Runge-Kutta method can be used to approximate u and thus we deduce approximations to the angles θ _i , the errors e _i , and the torques τ _i . The fourth-order Runge-Kutta method consists of generating a sequence { u ^{( i )} } as follows: Given the initial condition u (0) = u ⁽⁰⁾ , where

$$u^{(0)}={\left[x_1(0),x_2(0),\mathrm{...},x{}_n(0),{\theta }_1(0),{\theta }_2(0),\mathrm{...},{\theta }_n(0),{\dot{\theta }}_1(0),{\dot{\theta }}_2(0),\mathrm{...},{\dot{\theta }}_n(0)\right]}^t,$$

we generate u ^{( i +1)} for i 0,1,…,m using

$$u^{(i+1)}=u^{(i)}+\frac{h}{6}(k_1+2k_2+2k_3+k_4),$$

where

$$\begin{array}{l}{k}_1=g(t_i,u^{(i)}),\\ k_2=g(t_i+\frac{h}{2},u^{(i)}+\frac{h}{2}{k}_1),\\ k_3=g(t_i+\frac{h}{2},u^{(i)}+\frac{h}{2}{k}_2),\\ k_4=g(t_i+h,u^{(i)}+h{k}_3).\end{array}$$

Here, $h=\frac{T}{m}$ is the step size and t _i = ih for I = 0,1,…,m.

Numerical Experiments

In this section, we present several numerical results to show the performance of the PID controller. We consider n -link robot manipulators with n = 3,4,5,6. The case n = 2 is considered in [5].

Example 6.1. (Simulation results for the three-link robot manipulator): In this example, we consider a 3-link robot manipulator. We take the following parameters: m ₁ = m ₂ = m ₃ = 1, l ₁ = 2, l ₂ = l ₃ = 1. We choose the final position as ${\theta }_{1f}=\frac{\pi }{2},{\theta }_{2f}=0,{\theta }_{3f}=0.$ The initial conditions are taken as follows

$$\left[\begin{array}{c}{\theta }_1(0)\\ {\theta }_2(0)\\ {\theta }_3(0)\end{array}\right]=\left[\begin{array}{c}\frac{\pi }{2}\\ \frac{\pi }{2}\\ \frac{\pi }{2}\end{array}\right],\left[\begin{array}{c}{\dot{\theta }}_1(0)\\ {\dot{\theta }}_2(0)\\ {\dot{\theta }}_3(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right].$$

We choose the following PID parameters

$$K_{P_1}=K_{P_2}=K_{P_3}=30,K_{D_1}=15,K_{D_2}=K_{D_3}=10,K_{I_1}=K_{I_2}=K_{I_3}=20.$$

The initial and final positions of the robot manipulator are shown in Figure 2, Figure 3, Figure 4, and Figure 5 demonstrate θ _i ( t ), e _i ( t ) and τ _i ( t ) over time interval [0,30] for i = 1,2,3.

Example 6.2. (Simulation results for the three-link robot manipulator): In this example, we repeat the previous example (Example 6.1) with the following values: m ₁ = m ₂ = 1, m ₃ = 3, l ₁ = 2, l ₂ = 4, l ₃ = 1. The final position is $({\theta }_{1f},{\theta }_{2f},{\theta }_{3f})=\left(\frac{\pi }{2},\frac{\pi }{2},0\right).$ The initial conditions are

$$\left[\begin{array}{c}{\theta }_1(0)\\ {\theta }_2(0)\\ {\theta }_3(0)\end{array}\right]=\left[\begin{array}{c}\frac{\pi }{2}\\ 0\\ \frac{\pi }{2}\end{array}\right],\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\\ {\dot{\theta }}_3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right].$$

The PID parameters are

$$K_{P_1}=K_{P_2}=K_{P_3}=30,K_{D_1}=15,K_{D_2}=K_{D_3}=10,K_{I_1}=K_{I_2}=K_{I_3}=10.$$

The initial and final positions of the robot manipulator can be observed in Figure 6. To further illustrate the dynamics, Figure 7, Figure 8, and Figure 9 provide insight into the variations of θ _i ( t ), e _i ( t ) and τ _i ( t ) for i = 1,2,3, spanning the time interval [0, 30] (initial 30 seconds).

Example 6.3. (Simulation results for the three-link robot manipulator): We repeat the previous example with the following values: m ₁ = 1, m ₂ = m ₃ = 2, l ₁ = l ₂ = l ₃ = 1. The final position is $({\theta }_{1f},{\theta }_{2f},{\theta }_{3f})=\left(\frac{\pi }{3},0,\frac{\pi }{4}\right).$ The initial conditions are

$$\left[\begin{array}{c}{\theta }_1(0)\\ {\theta }_2(0)\\ {\theta }_3(0)\end{array}\right]=\left[\begin{array}{c}\frac{\pi }{6}\\ 0\\ \frac{\pi }{5}\end{array}\right],\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\\ {\dot{\theta }}_3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right].$$

The PID parameters are

$$K_{P_1}=K_{P_2}=K_{P_3}=20,K_{D_1}=15,K_{D_2}=K_{D_3}=10,K_{I_1}=K_{I_2}=K_{I_3}=20.$$

Figure 10 displays the initial and final positions of the robot manipulator. Figure 11, Figure 12 and Figure 13 illustrate the temporal variations of θ _i ( t ), e _i ( t ) and τ _i ( t ) for i = 1,2,3, over the time interval [0, 30] (during the first 30 seconds).

Example 6.4. (Simulation results for the four-link robot manipulator). In this experiment, we consider a 4-link robot manipulator. We utilize the following values: m ₁ = 4, m ₂ = 3, m ₃ = 2, m ₄ = 1, l ₁ = 4, l ₂ = 3, l ₃ = 2, l ₄ = 1. The final position is $\left({\theta }_{1f},{\theta }_{2f},{\theta }_{3f},{\theta }_{4f}\right)=\left(\frac{\pi }{2},0,0,\frac{3\pi }{2}\right)$ .

The initial conditions are taken as follows

$$\left[\begin{array}{c}{\theta }_1(0)\\ {\theta }_2(0)\\ {\theta }_3(0)\\ {\theta }_4(0)\end{array}\right]=\left[\begin{array}{c}\frac{\pi }{2}\\ \frac{\pi }{2}\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\\ {\dot{\theta }}_3\\ {\dot{\theta }}_4\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\\ x_4(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right].$$

Our PID parameters are

$$K_{P_1}=K_{P_2}=K_{P_3}=K_{P_4}=30,K_{D_1}=15,K_{D_2}=K_{D_3}=K_{D_4}=10,K_{I_1}=K_{I_2}=20,K_{I_3}=K_{I_4}=10.$$

Figure 14 presents the initial and final positions of the robot manipulator. Additionally, Figure 15, Figure 16, Figure 17 and Figure 18 depict the temporal variations of θ _i ( t ), e _i ( t ) and τ _i ( t ) for i = 1,2,3,4, within the time interval [0, 30] (during the initial 30 seconds).

Example 6.5. (Simulation results for the five-link robot manipulator). Here, we consider a 5-link robot manipulator. We take m ₁ = 25, m ₂ = 20, m ₃ = 15, m ₄ = 10, m ₅ = 5, l ₁ = 5, l ₂ = 4, l ₃ = 3, l ₄ = 2, and l_5=1. The final position is $({\theta }_{1f},{\theta }_{2f},{\theta }_{3f},{\theta }_{4f},{\theta }_{5f})=\left(\frac{\pi }{4},0,\frac{\pi }{4},0,\frac{\pi }{4}\right).$ The initial conditions are

$$\left[\begin{array}{c}{\theta }_1(0)\\ {\theta }_2(0)\\ {\theta }_3(0)\\ {\theta }_4(0)\\ {\theta }_5(0)\end{array}\right]=\left[\begin{array}{c}\frac{3\pi }{2}\\ \frac{3\pi }{2}\\ \frac{3\pi }{2}\\ \frac{3\pi }{2}\\ \frac{3\pi }{2}\end{array}\right],\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\\ {\dot{\theta }}_3\\ {\dot{\theta }}_4\\ {\dot{\theta }}_5\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\\ x_4(0)\\ x_5(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right].$$

For this experiment, we choose the following PID parameters $K_{P_1}=K_{P_2}=K_{P_3}=K_{P_4}=K_{P_5}=30$, $K_{D_1}=15,K_{D_2}=K_{D_3}=10,K_{D_4}=K_{D_5}=5$, $K_{I_1}=K_{I_2}=20,K_{I_3}=K_{I_4}=10,\text{and}{K}_{I_5}=5$ . The initial and final positions of the robot manipulator are shown in Figure 19, Figure 20, Figure 21, Figure 22, Figure 23, and Figure 24 illustrate the variations of θ _i ( t ), e _i ( t ) and τ _i ( t ) for i = 1,2,…,5, across the time interval [0, 30] (during the first 30 seconds).

Example 6.6. (Simulation results for the six-link robot manipulator): In this final experiment, we consider a 6-link robot manipulator. For this test we utilize the following values: m ₁ = 6, m ₂ = 5, m ₃ = 4, m ₄ = 3, m ₅ = 2, m ₆ = 1, l ₁ = 6, l ₂ = 5, l ₃ = 4, l ₄ = 3, l ₅ = 2, and l ₆ = 1. The target position is

$$({\theta }_{1f},{\theta }_{2f},{\theta }_{3f},{\theta }_{4f},{\theta }_{5f},{\theta }_{6f})=\left(\frac{\pi }{2},\frac{\pi }{3},\frac{\pi }{4},\frac{\pi }{5},\frac{\pi }{6},\frac{\pi }{7}\right).$$

The initial conditions are

$$\left[\begin{array}{c}{\theta }_1(0)\\ {\theta }_2(0)\\ {\theta }_3(0)\\ {\theta }_4(0)\\ {\theta }_5(0)\\ {\theta }_6(0)\end{array}\right]=\left[\begin{array}{c}\frac{\pi }{2}\\ \frac{\pi }{2}\\ \frac{\pi }{2}\\ \frac{\pi }{2}\\ \frac{\pi }{2}\\ \frac{\pi }{2}\end{array}\right],\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\\ {\dot{\theta }}_3\\ {\dot{\theta }}_4\\ {\dot{\theta }}_5\\ {\dot{\theta }}_6\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\\ x_4(0)\\ x_5(0)\\ x_6(0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right].$$

The PID parameters are $K_{P_1}=K_{P_2}=K_{P_3}=K_{P_4}=K_{P_5}=K_{P_6}=30$, $K_{D_1}=15,K_{D_2}=K_{D_3}=K_{D_4}=K_{D_5}=K_{D_6}=10$ , $K_{I_1}=K_{I_2}=20$ , $K_{I_3}=K_{I_4}=K_{I_5}=K_{I_6}=10$ . The initial and final positions of the robot manipulator are resented in Figure 25. Figure 26, Figure 27, Figure 28, Figure 29, Figure 30, and Figure 31 illustrate the behaviours of θ _i ( t ), e _i ( t ) and τ _i ( t ) for i = 1,2,…,6 over the time interval [0,30] (initial 30 seconds).

Conclusion

In this paper, we introduced a Proportional-Integral-Derivative (PID) control approach for an n-link robot manipulator, aimed at guiding its motion to achieve target positions from given initial conditions. We first derived the system’s equations of motion using the Lagrangian formulation. Subsequently, the PID is implemented, and the resulting second-order system of ordinary differential equations is approximated using the classical fourth-order Runge-Kutta method for systems. The efficiency of the control strategy is demonstrated through a series of numerical experiments, involving the simulation of the manipulator’s motion from an initial position to a desired target position. This paper focuses on the n-link robot manipulator in two-dimensional space. The PID controller introduced in this paper can be extended to accommodate robots with prismatic joints. In the future, we are planning to extend our approach to the n -link robot manipulator in three-dimensional space, where we will consider other typical joint types, including revolute (the attached links rotate about a common axis), prismatic (the attached links translate about a common axis) and spherical (the attached links rotate about a point).

Conflict of Interest

The authors declare that they have no conflicts of interest.

Ethical Statement

The authors declare that the material is their original work, which has not been previously published and is not currently being considered elsewhere.

Data Availability Statement

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

---