Optimal Design Oriented Adaptive Robust Control for A Planar 2-DOF Redundantly Actuated Parallel Robot

An adaptive robust control combined with a multi-objective parameter optimization method for the parallel robot with unknown uncertainty is proposed. In the active joint space, the accurate dynamic model of the parallel robot can be obtained by combining the closed-chain constraint force and the open-chain system’s dynamic equation. According to the Udwadia-Kalaba theory, the closed-chain constraint force imposed by the end eﬀector can be calculated in a simple way. The proposed adaptive robust control could guarantee deterministic robust performances of the system (the uniform boundedness and uniform ultimate boundedness). To seek suitable weighting factors for the proposed control, a system performance function, which includes the transient performance portion, the steady state performance portion, and the control cost, is introduced. By applying D -operation, the performance function is transformed into a multi-objective function with weighting factors. Mean-while, the problem of choosing the optimal gain is equivalent to the problem of ﬁnding the minimum value of the system performance index function. An illustrative example illustrates the superiority of the proposed modeling method and the proposed adaptive robust control.


Introduction
Compared with the traditional serial robot, the parallel robot is a closed-loop motion structure with multiple kinematic chains, which makes the kinematic and dynamic analysis of the parallel robot tough. But it has superiorities that traditional serial robots do not have, such as high rigidity, low inertia, and so on. Because of the above high performances, it has been widely used in many fields of science and engineering [1][2][3]. In order to realize its superiorities, the dynamic modeling and design of a proper controller are essential. However, the coordinated operation of the end effector by multiple kinematic chains poses a challenge to the design of motion control. Generally, current controls of the parallel robot are mainly composed of two branches, one is the kinematic controller and the other is the dynamic controller.
For the kinematic controller of the parallel robot, there are many classic control methods such as PD controller [4] and nonlinear PD (NPD) controller [5]. PD controller does not need a dynamic solution. Its structure and the calculation are simple. NPD controller is designed to improve system control performance in the nonlinear system. In the high-speed condition, the nonlinear dynamic characteristics of the parallel robots are more prominent, while the control effect of the kinematic controller is often less than satisfactory.
Unlike the kinematic controller, the design of the dynamic controller is based on the actual dynamic equation of the parallel robot. Fichter and Merlet [6] are the first scholars to explore the dynamic modeling of parallel robots and established the dynamic model of the Stewart robot. Some classic approaches, such as the Lagrangian approach [7] and Newton-Euler approach [8], have been widely applied. Others like d'Alembert's principle, Kane equation, and so forth also solved many dynamics problems. Reference [9] established the dynamic model of the redundant actuation parallel robot in generalized space by the application of the Udwadia-Kalaba method. Solving dynamic problems is to provide the design foundation for the controller design.
For dynamic controllers, the system often has a good control effect especially in high-speed and high-precision situations, because it takes the dynamic characteristics of parallel robots into consideration. Reference [10] analyzed the dynamics of the simplified rigid body model of parallel robot. Then, the expression of inertia moment was obtained, and the computed torque method was proposed to improve the augmented PID controller. This method can reduce the difficulty of decoupling the control algorithm. Constraint-following control [11] shows good control ability in the dynamic control strategy. The LQR controller was introduced in [12] which presented a comparison of the time specification performance between the LQR controller and PD controller for a double inverted pendulum system. The sliding mode control [13], due to its variable structure characteristics, is robust to nonlinear systems with uncertain parameters, especially for parallel robots. From the perspective of improving the robustness of the mechanism, reference [14] designed a variable structure sliding mode controller to solve the instability problem.
In a high-speed motion environment, the impact of external uncertainty on system performance cannot be ignored. Hence, a large number of control methods based on dynamic model emerge for the uncertain mechanical system. When the bound of uncertainty has been obtained, many controllers can effectively solve the uncertainty problem such as H ∞ control [15], robust control [16] and so on. Adaptive robust control [17] has good control performance under the unknown uncertainty. The uncertainty of parallel robot was handled by estimating the uncertain parameters online [18][19][20] and the stochastic properties of uncertain terms were explored at the same time [21]. A fuzzy approach is proposed in [22][23][24] to describe the uncertainty by the fuzzy set theoretic. The uncertainty's information tends to be unknown but its boundary is within a fuzzy threshold.
The main contributions of this work are threefold. Firstly, the planar 2-DOF parallel robot is we divided into several subsystems and the open-chain dynamic equation is obtained. In order to quickly and effectively obtain the complete dynamic model, a novel approach to calculate the constraint force imposed by the end effector is introduced. This novel analytical mechanics approach can calculate the constraint force (holonomic constraints and/or nonholonomic constraints) in an analytical form without any auxiliary variables such as Lagrange's multiplier. Secondly, in the high-speed motion environment, the impact of external uncertainty on system performance cannot be ignored. By using a simple function to describe the bound of uncertainty and a given adaptive law to estimate the unknown parameter, an adaptive robust control with weighting factors is proposed. Thirdly, by system performance analysis, the system has the uniform ultimate boundedness property and uniform ultimate boundedness property. However, as the system's control performance is enhanced, the control costs increase accordingly. In order to keep a proper balance between the system performance and the control costs, we need to seek the optimal gains. Then, the optimal gain design problem is described by a system performance index which includes a transient portion, a steady state portion and a control cost part. The optimal gain design problem could be solved by finding the global minimum of the system performance index.

Constraints of mechanical system
The dynamic equation of a general mechanical system without constraints is given by [25]: here, ϑ ∈ R n×1 is an n-dimensional generalized coordinate. In general, it means the joint angles. The velocity vector,θ ∈ R n×1 , is the first derivative of ϑ with respect to time t, andθ ∈ R n×1 represents the acceleration vector. The M (ϑ, t) ∈ R n×n is the inertia matrix which is positive definite and M (ϑ, t)=M T (ϑ, t). The F (θ, ϑ, t) ∈ R n×1 is the given force that includes Coriolis/Centrifugal force, gravitational force, and other control input forces of the unconstrained system. The functions M (ϑ, t) and F (θ, ϑ, t) are considered continuous.
In practice, unconstrained system is often affected by constraint force. Suppose there are m(m < n) constraints in the following form: Eq.(2) is the constraint in the first-order form (Pfaffian representation). The functions A li (·) and c l (·) are both C 1 . These constraints may in holonomic form or nonholonomic form, and may be not integrable in general. The m means the number of constraints attached to the system. The matrix form of those constraints is given by where The expression of constraints likes Eq.(3) is more universal. Generally, the constraint often includes two forms, namely active constraint and passive constraint. The active constraint is the control input which ensures the desired output. The control input of system is the main problem for engineers. The passive constraint is from the structure or the environment which could automatically generate the required constraint force.

Remark 1
The constraints of the first-order form (or zero-order form when A(ϑ, t) = 0) are widely applied to Gibbs and Appell equation, Lagrange's motion equation, Maggi equation, and so on. However, the existing analytical methods are based on some auxiliary variables (such as Lagrange multipliers, projections, or any quasi or auxiliary variables) rather than analytical expressions.
Udwadia and Kalaba introduced second-order constraints (instead of firstorder constraints) into the mechanical system which is believed to be an essential breakthrough in analytical mechanics. By differentiating Eq.(2) with respect to t, we could obtain the second-order form of the constraints as follows: Next, put the second-order derivative of ϑ on the left and move the rest to the right side. Let by combining Eq.(4) and Eq. (5), the matrix form of the second-order constraint is given as where

Assumption 1
The constraint in the second-order form of (6) is consistent, i.e., there exists at least one solutionθ which satisfies the Eq.(6) for the given A(·) and c(·).

Remark 2
The information related to constraints should not be lost for differential operations and it has been illustrated by most control problems such as stabilization, trajectory following control, and so on.
To satisfy the constraint Eq. (6), it is necessary to provide the system Eq.(1) with an active manner constraint force F c (θ, ϑ, t), which satisfies that the virtual work done by the constraint force is zero (the d'Alembert's principle). Therefore, when the mechanical system has additional constraint force F c (θ, ϑ, t), it is called a "constrained mechanical system", and its dynamic equation like Eq.(1) could be rewritten as following: Consider the given A(·) and c(·), there exists at least one solutionθ. All the solutions may be called the possible acceleration. The problem of how to calculate the constraint force is actually to find the only actual acceleration which conforms to the d'Alembert's principle among all possible accelerations. In other words, the corresponding constraint force F c (θ, ϑ, t) ∈ ℜ(A T (ϑ, t)). Udwadia and Kalaba provided the analytical form of the only actual acceleration as [26]q where B(ϑ, t) = A(ϑ, t)M 1 / 2 (ϑ, t). The "+" is the Moore-Penrose generalized inverse. By combining Eq.(1) and Eq.(8), we have (9) Comparing Eq.(7) and Eq.(9), we can easily get Because Udwadia and Kalaba pioneered the method, the equation was named after them, called the Udwadia-Kalaba equation. It could be proved by d'Alembert's principle or extended d'Alembert's principle and Gauss principle. The constraint force ensures the system meets the constraint with analytical expression.
Remark 3 Consider a constrained mechanical system subject to constraints, the constraint force can be obtained by Udwadia-Kalaba equation in analytic expressions, which does not need other auxiliary variables and just based on the same generalized coordinates/velocities of the original system. Consider a 2-DOF redundantly actuated parallel robot in Fig. 1, we will formulate its precise dynamic model first. This parallel robot is composed of three subsystems connected in the same plane. The end effector is installed where the three subsystems are connected to each other. Each subsystem is a two-link mechanism. The active links of the three subsystems, called the active joints or driving joints of the parallel robot, are connected to the base. They are driven by the servo motors installed on the base. The joint between the active link and the passive link is not driven, which is called the passive joint.

Dynamic modeling based on the active joint space
The end effector is installed at the joint of the three passive links. The end effector can move in the X and Y directions and is driven by three active inputs. Therefore, the parallel robot belongs to a redundant drive type parallel robot. The robot can be regarded as a traditional five-link mechanism with a branch chain attached. This redundant branch chain can eliminate the singularity of the five-link mechanism, improve the force transmission performance, and obtain the only forward kinematics solution. At the same time, the redundant drive makes the model of the mechanism more complex. So we introduced a novel Udwadia-Kalaba method to quickly and effectively obtain the dynamic model of the robot. First, the parallel robot is divided into 3  Fig. 2). The m ai and m bi (i=1,2,3) represent the mass of passive link and active link. The C ai and C bi mean the position of the mass center. The r ai and r bi mean the distance between the mass center and the corresponding joint. The l represents the link length. For each simple subsystem, through Lagrangian mechanics, the dynamic equation of the subsystem could be obtained easily. Particularly, since the influence of gravity and friction are ignored in the plane, the dynamic equation of the subsystem is given by [26] where t ∈ R is an independent variable, ϑ i = [ϑ ai , ϑ bi ] T is the generalized coordinates of i-th subsystem. Theθ i andθ i represent the velocity vector and acceleration vector respectively. The τ i = [τ ai ,τ bi ] T is control input torque, denotes Coriolis/Centrifugal force matrix. The functions M i (ϑ i , t) and C i (θ i , ϑ i , t) satisfy the corresponding condition of Eq.(1) and they are given by Here, X i , Y i , Z i are constants related to the physical parameters of the subsystems, i.e.,      where the J ai and J bi represent the moment of inertia relative to the mass center of the corresponding link. Then combine the dynamic models of three two-link mechanisms to get the dynamic equation of the entire open-chain system, we get with The subscript o represents the open-chain system and the superscript indicates the position (the row number and the column number) in the matrix. Consider the open-chain equation (14), the complete dynamic equation cannot be obtained without the closed chain constraint force. The constraint force which puts three subsystems together turns the unconstrained open-chain system into a complete closed-chain system, which is called "constrained mechanical system". By introducing the Udwadia-Kalaba equation, the constraint force can be easily obtained.
Suppose there are some constraints imposed by the end effector in the unconstrained open-chain system. The constraints ensure that the end positions of the three subsystems coincide. The detailed expressions of the constraint equation are By differentiating the constraint Eq.(16) twice, we can convert the zero-order form constraints into the second-order form constraints as where where So far, we have obtained the dynamic equation (Eq. (14)) of the open-chain system and the kinematic constraints (Eq. (19)) of the parallel manipulator. Then, in order to obtain the complete dynamic model, the additional constraint forces must be imposed on the open-chain system. Thus the final dynamic model of the closed-chain system is Remark 4 Now, we have the complete dynamic model of the parallel robot showed in Fig. 1 by Udwadia-Kalaba theory. This novel approach can calculate the constraint force (holonomic constraints and/or non-holonomic constraints) in an analytical form, just by its generalized coordinates. So the constraint force could be expressed by the generalized coordinates of the system itself.
Having the dynamic equation like Eq.(20) is still not enough, because the input torque contains passive joints. However, passive joints do not have input torque. So, another problem is to convert Eq.(20) into the active joint space, which only contains the input of active joints. Eq. (20) can be rewritten as where By kinematic analysis [27], we havė The detailed expression of N is By differentiating Eq. (22), the acceleration relationship between the active joints and the passive joints isθ Combining Eq.(21), Eq. (22) and Eq. (24), we have where E means the identity matrix. Let P = , then multiplying both sides of Eq.(25) by P T , we get Since the passive joint has no input torque, (26) can be simplified to Remark 5 So far, the problem of dynamic modeling is solved. This complete dynamic model is based on the active joint space which provides a practical basis for the design of the following controller.

Mechanical system with uncertainty
Although complete dynamic model based on the active joint space is obtained, there are always many unavoidable uncertain disturbances in reality. The system with uncertainty could be described as where δ ∈ Σ is an uncertain parameter vector. The bounding set Σ is prescribed and compact. The ϑ a = [ϑ a1 , ϑ a2 , ϑ a3 ] T is the active joint angles of the system,θ a andθ a are the corresponding differential of ϑ a over t, τ a = [τ a1 , τ a2 , τ a3 ] T means the active inputs, M (ϑ a , δ) is a positive definite inertia matrix and M (ϑ a , δ)=M T (ϑ a , δ). The C a (θ a , ϑ a , δ) denotes Coriolis/Centrifugal force matrix, F (θ a , ϑ a , δ) means the constraints imposed by the end effector. The functions M (ϑ a , δ), C(θ a , ϑ a , δ) and F (θ a , ϑ a , δ) are continuous. The τ a is the control input.

Assumption 2
The M (ϑ a , δ) is considered to be uniformly positive definite and uniformly ultimately upper bounded. That means there always exists con- here, E means the identity matrix. We state that M (ϑ a , δ) ≤ λ M −max E is not necessarily satisfied and more extensive properties could be found in [23]. But if there are only revolute joints in the system without slide joints or prismatic joints, the inertia matrix M (ϑ a , δ) should be uniformly positive definite and uniformly ultimately upper bounded so that Eq.(29) is strictly satisfied.
Assumption 3 (i) Consider the initial state vector ϑ a0 which consists of ϑ a0i , i = 1, 2, 3, there always exists a fuzzy set named X 0i in a universe of discourse I i ∈ R which is characterized by a membership function U Θi : (ii) Consider the uncertain parameter vector expressed by δ i , i = 1, 2, 3, for any δ i (·), there always exists a fuzzy set named H i in a universe of discourse (iii) Furthermore, the function δ(·) is Lebesgue measurable.
By the fuzzy description, the system with uncertainty is called the fuzzy dynamical system. However, it should be noted that the fuzzy model in our work differs from the Takagi-Sugeno model or other if-then rule-based fuzzy models.
We define the tracking error of the active joint angles as where ϑ ad represents the ideal trajectory of the active joint.

Remark 6
Here, ϑ e (t) and ϑ a (t) are actually equivalent, because the ideal trajectory ϑ ad (t) is known. But for the convenience of expression and the reality of the project, ϑ e (t) needs to be mentioned. Therefore, in the following text, the function represented by ϑ e (t) is actually equivalent to the function represented by ϑ a (t), just a different form of expression.
Assumption 4 (i) There exists a function named g(θ e , ϑ e , δ) which represents the total uncertainty of the whole system. The detailed expression of g(θ e , ϑ e , δ) is According to [23], there exists a corresponding dimensional constant vector (ii) Consider all members of the parameter vectorψ(t), namelyψ(t) i , i = 1, 2, 3, there always exists a fuzzy set named W i in a universe of discourse The function Π(θ a (t), ϑ a (t),ψ(t)) is a suitable upper bound for g(θ e , ϑ e , δ). Such that for all ϑ e (t) It should be noted that the function Π(θ a (t), ϑ a (t),ψ(t)) is not arbitrary. It must meet the following assumptions.
Although the bound of uncertainty is not know, the existence of such bound is doubtless. Then we assume that there is a function named Π(θ a (t), ϑ a (t),ψ(t)). The structure of this function is known, but one of the parametersψ(t) is unknown. We estimate the value ofψ(t) throughψ(t) so that the value of Π(θ a (t), ϑ a (t), ·) will be infinitely close to the actual value. The accurate bound of uncertainty will be obtained. The estimated valueψ(t) is used in the controller design, and it satisfies a given adaptive law. That is where k is a given constant positive scalar, S = diag[S ii ] 3×3 , S ii > 1, i = 1, 2, 3. The scalar k determines the "rate" of adaptation. If k is appropriately large, the adaptive parameters will increase rapidly, so that the control input changes correspondingly in order to compensate for the uncertainty effectively.

Remark 7
Here, g(θ e , ϑ e , δ) represents the total uncertainty of the whole system, in other words, it means the bound of uncertainty. However, it is obvious that the form of g(θ e , ϑ e , δ) is very complicated. How to get the value of g(θ e , ϑ e , δ) is a challenging problem. But using a simple function Π(θ a (t), ϑ a (t), ·) to describe the bound is a reasonable idea.

Adaptive robust control design
We design the adaptive robust control for the system as following: where τ a (t) means the input torque which consists of τ 1 (t) and τ 2 (t). The τ 1 (t) stands for the robust control part whose main contribution is to compensate for the uncertainty. The τ 2 (t) stands for the PD control part, and it ensures that the controller has the ability to track the ideal trajectory. The α and β are positive constants which respectively represent weighting factors of adaptive robust portion and PD portion. Subject to the proposed control input, the system should achieve the following goals. As t approaches to infinity, the tracking error ϑ e (t) and its first-order derivativeθ e (t) approaches to zero. Meanwhile, the parameterψ(t) approaches toψ(t).
In general, Eq.(28) could be rewritten as following for convenience:
Proof Select a Lyapunov function candidate for the dynamical system Eq.(39) as following: To prove V is a legitimate Lyapunov function for the fuzzy dynamical system, we must prove that V is positive definite and decreasing globally. According to Assumption 2, we have here with Based on Eq.(43), we know λ min (Λ i ) = λ M −min , so Moreover, where λ 1 = min{λ M −min , 1 2 k}. So, we obtain the possible lower bound of V . Next, we analyze the upper bound. where Based on Eq.(43) we know λ max (Λ i ) = 2βS ii , So Moreover If k is appropriately large, the adaptive parameters will increase rapidly, so that the control input changes correspondingly in order to compensate for the uncertainty effectively. If we choose k ≥ 2β max(λ S ) and λ 2 = 1 2 k, the information of upper bound is obtained.
To illustrate that V is decreasing, we analyze the derivative of it. The derivative of V is given byV We analyze the two parts of V separately. First, By combining Eq.(33) and Eq.(51), we havė V ϑe = (θ e + Sϑ e ) T (τ a + g) + +2βϑ e T Sθ e .

Π(q(t), q(t),ψ(t)) − Π(q(t), q(t),ψ(t))
) According to Eq.(54) and Eq.(56), the complete expression ofV is as folloẇ According to the inequality −ax 2 − x ≤ 1 4a , a ∈ R, we can geṫ As being stated above, we know k ≥ 2β max(λ S ). Sȯ whereλ = βλ 3 ,ĥ = βλ 3 ψ 2 + 1 4α . Now, we have been proved that V is decreasing if From [24], we can obtain the uniform boundedness performance. That is where In addition, the uniform ultimate boundedness performance is as following: Consider allχ > χ min Remark 8 The value of χ min determines the radius of uniform ultimate boundedness ball. It is clear that χ min will decrease correspondingly as Γ decreases. Meanwhile, Γ is determined by the weighting factors α and β. The Γ will decrease as α and β both increase. So the uniform boundedness performance of the system will be enhanced as α and β increase. However, the bigger α or β does not always mean the better. Because as they increase, the control input of the system also increases, which means greater control cost. Therefore, how to keep a proper balance between the control cost and the control performance is really necessary.

Optimal gain design based on D-operation
The system performance guaranteed by the control input has been analyzed. Through the analysis, it can be seen that the uniform boundedness performance will be enhanced when α and β increase. However, as they increase, the control cost increases accordingly, which becomes another issue. In order to keep a proper balance between the control cost and the system performance, we need to seek the optimal α and β. The so-called balance will be described by a system performance index function, which includes a transient portion, a steady state portion and a control cost portion. The problem of finding the optimal gains can be transformed into the problem of finding the minimum value of the system performance index. First of all, we will further explore the uncertain system. Based on the lower bound and the upper bound of V , we have and hence By combining Eq.(59) and Eq.(66), we havė The analysis of this differential inequality is explained in detail in Chen [28]. For further analysis of this differential inequality, the following is needed.

Definition 1 Consider a scalar function ω(s(t), t) with the scalars s(t) and
t in some open connected set D, we hold a function s(t) with t ∈ [t 0 , t 1 ] is a solution of the differential inequalitẏ if s(t) is continuous and its derivative satisfies Eq.(68) with any t ∈ [t 0 , t 1 ].

If s(t) and φ(t) are the solutions of Eq.(68) and Eq.(69) respectively and the initial condition satisfies s(t
Theorem 3 For the differential inequality Eq.(68) and the differential equation Eq.(69), for all (x 1 , t), (x 2 , t) ∈ D, there always exists a constant L > 0 which ensures the function ω(·) satisfies the Lipschitz condition Then, for all t ∈ [t 0 , t 1 ], any solution s(t) of the differential inequality Eq.(68) also satisfies the following inequality Exploring the solutions of Eq.(68) often involves many difficulties because the solution may be not available or non-unique. On the contrary, studying the upper bound information of the solution may be a feasible way and it can avoid the trouble caused by the non-unique solution. Because the analytic form of upper bound can be solved, it gives an approach to explore the performance index function.
Based on the above theorems and definition, a differential equation is constructed asφ with φ(t 0 ) = V (t 0 ). Then, by solving the differential Eq.(72), we have Therefore, combining Eq.(71) and Eq.(73), we get where V 0 = V e (t 0 )+V ψ (t 0 ), t 0 = 0. It should be noted that the initial condition does not necessarily to be is also satisfied. We assume that t 0 = 0 is just for the control input starts at t = 0. Notice that for each t ∈ [t 0 , ∞), let We can find that for each α, β andψ, η(α, β,ψ, t) approaches to 0 as t → ∞. So η(α, β,ψ, t) reflects the transient performance which depends on the initial state. On the other hand, η ∞ (α, β,ψ, t) represents the steady state performance.
The exact knowledge of the uncertainty may be not available. Therefore, it may be feasible to study η(α, β,ψ, t) and η ∞ (α, β,ψ, t) when analyzing the system performance. Although η(α, β,ψ, t) and η ∞ (α, β,ψ, t) are related to the value ofψ and the knowledge of exact value may not have, the value of ψ is characterized by a known membership function. In order to find a meaningful reference value about the fuzzy numberψ in the fuzzy set to accurately describe the function η(α, β,ψ, t) and η ∞ (α, β,ψ, t), the D-operation is introduced as following.

Definition 2 Consider a fuzzy set
and a function ζ(δ) : ∑ → R with δ as a parameter, the D-operation of ζ(δ) is defined as Remark 9 As mentioned before, the value ofψ is characterized by a membership function. Then the D-operation is introduced to find a meaningful reference value of the fuzzy numberψ in the fuzzy set so that we can accurately describe the function η(α, β,ψ, t) and η ∞ (α, β,ψ, t). A special case of D-operation while ζ(δ) = 1 is equivalent to the well-known center-of-gravity defuzzification method.
Now we define the system performance index function named J as here, (79) The J 1 reflects the transient performance which relates to the initial state, J 2 is related to the steady state performance. The J 3 part will increase as α and β increases. It represents the control cost.The control cost and the steady state performance are the key points to be considered in engineering practice. The parameters µ and ν are the weighting factors between the steady state performance and the control cost. Next, the optimal α and β under different weighting factors µ and ν will be further explored.
Based on the definition of D-operation, we have Under reasonable initial conditions, ρ 0 ∼ ρ 7 are generally normal numbers. The detailed expressions of parameters ρ 0 ∼ ρ 7 are Seeking the optimal parameters α and β can be transformed into a problem of finding the minimum of J. This means To find the minimum value of J, we first calculate the first-order partial derivative of J. For any initial condition From the information of the first-order partial derivative, it is not easy to directly find out when the function takes the minimum value. Therefore, taking the second-order partial derivative of J, we further have It should be declared that for any α and β, we have Based on Eq.(85), we know A > 0, C > 0. So for any solution α l and β l satisfying Eq.(86), J(α l , β l , V 0 ) is the local minimum.
Remark 10 For any α and β, ∂ 2 J ∂α 2 > 0, ∂ 2 J ∂β 2 > 0 are globally satisfied. Therefore, J(α, β, V 0 ) is a global concave function with respect to α and β so that the local minimum of J(α, β, V 0 ) is also the global minimum. Thus, the solution that satisfies Eq.(86) is indeed the optimal solution named α opt and β opt .  Until now, the problems about dynamic modeling, adaptive robust control design and optimal gain design are solved. Fig. 3 shows the main design procedure of proposed control scheme. The entire design process of the proposed control method is summarized as follows: Step 1: Obtain the dynamic equation of the entire open chain system like Eq.(14).
Step 2: Calculate the closed chain constraint force by Udwadia-Kalaba equation.
Step 3: Obtain the complete dynamic model of the parallel robot like Eq. (20). Then convert it into the active joint space.
Step 4: Choose suitable k and bound function Π so that we can formulate the adaptive laws by Eq.(37).
Step 6: Based on the system performance analysis and D-operation, construct the system performance index function J.
All initial values are based on SI units. For any input torque such as τ a = [0.1 cos(πt), 0.1 cos(πt), 0.1 cos(πt)] T , the actual output of the system is shown in Fig. 4 Under random input τ a , Fig. 4 shows the trajectory of the three subsystems' ends on the plane. Due to the constraint force imposed by the end effector, the phenomenon of the end trajectory overlap between the three subsystems shown by the system output is consistent with the actual situation. Therefore, the proposed dynamic modeling method is proved to be effective.
On the other hand, to illustrate that the proposed control has the ability to deal with uncertain interference, we assume that the mass of the links is uncertain and there exists external disturbing force. Let m ai =m ai + ∆m ai , m bi = m bi +∆m bi . Them ai andm bi means the known values, ∆m ai and ∆m bi means unknown values related to uncertainty. We define the uncertainties in m ai and m bi are the same, that is ∆m ai = ∆m bi = ∆m. ∆m is "close to 0.25" with a membership function (sine) as It needs to be stated that the sine function is not the only suitable membership function for the system. Before introducing the fuzzy set to describe the uncertainty of the system, the uncertainty should be studied first. Then the correspondingly suitable membership function for the system could be chosen reasonably.
Next, the bounding function Π is given as here,ψ satisfies the following adaptive laẇ According to step 5 and step 6, we choose k = 12, S = 2, λ 3 = 1. Then, we will calculate the optimal α and β under different weighting factors µ and ν. Based on the fuzzy description, we define the estimated valueψ is also "close to 0.25" with the membership function (sine) as Based on the D-operation, we can get ρ 0 = -0.02298 + 1.24611µ By the membership function Eq.(91) and D-operation in Definition 2, we can obtain the optimal gain α opt and β opt by Eq.(86) under different weighting factors µ and ν. The optimal solutions are summarised in Table 1. From this table, we notice that the optimal gain (α opt , β opt ) will increase as the weighting ratio (µ/ν) increases. Another point that should be noted is that different weight ratio µ/ν represents different needs in actual situations. If the engineers pay more attention to control performance in practice, they can choose a larger weight ratio. On the contrary, if more attention is paid to the control cost, a lower weight ratio should be chosen.
Until now, all parameters have been obtained. Next, we choose the desired trajectory of the output like { x E (t) = b(a cos(ωt) + R cos(ωt) cos(nωt)) y E (t) = b(a sin(ωt) + R sin(ωt) cos(nωt)), here, ω = 1, n = 5, R = 0.5, a = 3, b = 0.06/3.5. Meanwhile, we choose ∆m = 0.25 sin(πt) to describe the uncertainty caused by random mass changes and a randomτ to describe the uncertainty caused by disturbance torque. All the simulation results are as follows.  Fig. 5, it can be concluded that the solution of ∂J ∂α = ∂J ∂β = 0 is unique. In other words, the optimal gain is also unique. Fig. 6 shows the relationship between the system performance index function J and (α, β). Combining with Fig. 6 and Eq.(85), we know that J(α, β, V 0 ) is a global concave function with respect to (α, β). Therefore, the solution satisfying ∂J ∂α = ∂J ∂β = 0 is the unique optimal gain (α opt , β opt ). The minimum value of the system performance index function does not have multiple solutions. In other words, the result of the proposed optimal design is unique. In general, the proposed method has a better control effect through random disturbance, initial position error, and mass change. Comparing with traditional PD control, Fig. 7 shows the actual tracking performance of the proposed method. For the same initial position error, the proposed method has a faster response and smaller steady-state error when tracking the desired trajectory. In addition, a random disturbance torque is given when t = 4s and the detailed results are shown from Fig. 8 to Fig. 12. It can be seen that the proposed method is more robust in dealing with random disturbance. Finally, we choose ∆m = 0.25 sin(πt) to describe the uncertainty caused by random mass changes in the whole process. The tracking performance of the proposed adaptive robust method is more accurate than that of traditional PD control. All the simulation results show that the proposed adaptive robust method has better tracking accuracy, more rapid response, and stronger robustness than traditional PD control.   Fig.8 to Fig.12, the tracking error of the proposed method is smaller than traditional PD control, which indicates better control performance. On the one hand, under the interference of mass changes, the tracking error under the proposed control converges into a neighborhood of zero rapidly in a finite time and then remained within a narrow zone thereafter. This phenomenon proves that the proposed adaptive robust control could guarantee deterministic robust performances of the system (the uniform boundedness and uniform ultimate boundedness). On the contrary, the tracking error under the classic PD control does not converge into a narrow region around zero. The phenomenon of error convergence shows the superiority of the proposed control in tracking performance.
On the other hand, when the system is affected by random disturbance, the proposed control method can react quickly and maintain stable tracking performance. On the contrary, the classical PD control is not strong enough in dealing with random disturbance, and cannot quickly converge to a stable state. From the perspective of tracking error, regardless of the random disturbance, the proposed control method has better tracking accuracy. In terms of the ability to deal with random disturbance, the proposed control method has stronger robustness and stability.
In addition, from the performance index function of the system, we know that a larger weight ratio µ/ν means that in practice we pay more attention to the steady state performance of the system. At the same time, when the system is in stable state, the tracking error of the system has a certain trend, that is, it decreases as the weight ratio µ/ν increases. So we can conclude that the tracking error in stable state will decrease as the weight ratio µ/ν increases when t → ∞. This regularity also proves the validity of the system performance index function.     Fig. 13 to Fig. 17, the most obvious difference is that the accumulative errors of the proposed control method are much smaller than the classic PD control. Regardless of random disturbance, the proposed adaptive robust control is superior to the classical PD control in terms of accumulative errors. This means that the tracking accuracy of the proposed control method is higher and the ability to deal with external uncertainty is stronger in the whole process. Comparing the system with disturbance, the proposed control method has a smaller accumulative error when there is no disturbance. However, the classic PD control has no such regularity, so the proposed adaptive robust control is more stable. Generally speaking, the proposed adaptive robust control has better control performance under different performance indexes.

Conclusion
In this work, a novel dynamic modeling method and a systematic control design approach for a planar 2-DOF redundantly actuated parallel robot are proposed. Based on Udwadia-Kalaba equation, the constraint force imposed by the end effector can be calculated more effectively. This novel analytical mechanics approach can calculate the constraint force (holonomic constraints and/or nonholonomic constraints) in an analytical form without any auxiliary variables such as Lagrange's multiplier. Combining the closed-chain constraint force and the dynamic equation of the open-chain system, we can obtain the complete dynamic model in the active joint space. In order to ensure the performance of the system under external uncertainty, an adaptive robust control is proposed to guarantee the uniform boundedness and the uniform ultimate boundedness. The unknown external uncertainty is described by the fuzzy set theory and the bound of uncertainty is described by an estimated value that satisfies a given adaptive law. The proposed adaptive robust control consists of an adaptive robust control part and a PD control part. In order to seek suitable weighting factors for these two parts, a system performance index function is adopted, which includes the transient portion, the steady state portion, and the control cost. The problem of finding the optimal gain can be transformed into the problem of finding the minimum value of the system performance index function. By analysis, there is a unique solution for the optimal gain. We conclude that the proposed adaptive robust method has better tracking accuracy, more rapid response, and stronger robustness than the traditional control method. The simulation results prove the superiority of the proposed adaptive robust control. Meanwhile, the regularity reflected in the simulation results also proves the validity of the system performance index function.

Conflict of Interest
The authors declare that they have no conflict of interest.