Optimal design of adaptive robust control for a planar two-DOF redundantly actuated parallel robot

An adaptive robust control combined with a multi-objective parameter optimization method for a parallel robot under unknown uncertainty is proposed. In the active joint space, the model of the parallel robot can be obtained by combining the closed-chain constraint force and the open-chain system’s dynamic equation. Based on the Udwadia–Kalaba theory, the closed-chain constraint force imposed by the end effector can be calculated if no uncertainty. If there is uncertainty, we propose an adaptive robust control, based on a set of feasible design parameters, which guarantees deterministic performances, including uniform boundedness and uniform ultimate boundedness. To seek the optimal choice of the design parameters, among the feasible pool. We examine the system performance, which includes the transient performance, the steady state performance, and the control cost. By a fuzzy-theoretic D-operation, a performance index is formulated. The problem of choosing the optimal control parameters is equivalent to the problem of finding the minimum value of the performance index. An illustrative example demonstrates the superiority of the adaptive robust control. Mass uncertainty and external disturbance torque are introduced to test the effectiveness of the proposed control strategy.


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 highprecision 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. Constraintfollowing 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 fourfold. Firstly, the planar two-DOF parallel robot is 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 non-holonomic constraints) in an analytical form without any auxiliary variables. This is very different from the earlier work such as Lagrangian approach in [5].
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. Compared with [25], this paper verifies the ability of the proposed control method to deal with uncertain disturbances such as mass uncertainty and external disturbance torque. 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 gains design problem could be solved by finding the global minimum of the system performance index. Compared with previous similar studies, our system performance index could deal with multi-objective optimization problems without introducing more complex theories.
Fourthly, weighting factors are cited to reflect different actual situations. A larger weight ratio means that we pay more attention to the steady-state performance of the system, whereas a smaller weight ratio represents that we pay more attention to the control cost of the system. Therefore, the proposed system performance index function can meet different actual needs.
2 Dynamic modeling of two-DOF parallel robot

Constraints of mechanical system
The dynamic equation of a general mechanical system without constraints is given by [26]: here, t ∈ R is the independent variable, ϑ ∈ 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(·) and F(·) are considered continuous. In practice, unconstrained system is often affected by constraint force. Suppose there are m(m ≤ n) constraints in the following form: Equation (2) is the constraint in the first-order form (Pfaffian representation). The functions A li (·) and c l (·) are both C 1 . These constraints maybe in holonomic form and/or nonholonomic form. 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 first-order 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 A(ϑ, t) ∈ R m×n , b(θ, ϑ, t) ∈ R m×1 , rank(A) ≥ 1.

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

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 b, 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θ act as [27] where B(ϑ, t) = A(ϑ, t)M 1/2 (ϑ, t) and the symbol "+" denotes the Moore-Penrose generalized inverse. By combining Eq. (1) and Eq. (8), we have Comparing Eq. (7) and Eq. (9), we can easily get 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.

Dynamic modeling based on the active joint space
Consider a two-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.
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 subsystems (as shown in 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 [27] 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 satisfy the corresponding condition of Eq. (1). For simplicity, we use M i denotes M i (ϑ i , t) and C i denotes C i (θ i , ϑ i , t). 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 Eq. (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 closedchain 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 Here x Ei and y Ei are the coordinates of the i − th subsystem's end effector. By differentiating the constraint Eq. (16) twice, we can convert the zero-order form constraints into the second-order form constraints as where Based on the Udwadia-Kalaba equation Eq. (10), the constraint force in the analytical form is where . So far, we have obtained the dynamic equation of the open-chain system and the kinematic constraints 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 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. (19) 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. (19) into the active joint space, which only contains the input of active joints. Equation (19) can be rewritten as where By kinematic analysis [25], we havė The detailed expression of N is By differentiating Eq. (21), the acceleration relationship between the active joints and the passive joints is Combining Eq. (20), Eq. (21) and Eq. (23), we have where E means the identity matrix. Let P = E N , then multiplying both sides of Eq. (24) by P T , we get Since the passive joint has no input torque, (25) 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 , δ) denotes Coriolis/Centrifugal force matrix, F(θ a , ϑ a , δ) means the constraints imposed by the end effector. The functions M(·), C(·) and F(·) 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 constants here, E means the identity matrix. Assumption 2 was generally considered as a fact (as, for example, in the popular textbook "Control of Robot Manipulators in Joint Space" [28]), rather than an assumption. This means Assumption 2 covers a rather broad range of mechanical manipulators. It was not until in Reference [29] that we uncovered some exceptions, hence we now state this as an assumption. In addition, in Reference [29], we know that if there are only revolute joints in the system without slide joints or prismatic joints, the inertia matrix M should be uniformly ultimately upper bounded.
Assumption 3 (i) Consider the initial state vector ϑ a 0 which consists of ϑ a 0 i , 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 (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 Σ i ∈ R, which is characterized by a membership function (iii) Furthermore, the function δ(·) is Lebesgue measurable.
By the fuzzy description, the system with uncertainty is called the fuzzy dynamical system [22]. 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ψ(t) ∈ R n (maybe unknown) and a known function Π(θ a , ϑ a ,ψ) : R n × R n × (0, ∞) n → R + , the n and n represent the corresponding dimensions, such that (ii) The function Π(θ a , ϑ a ,ψ) is C 2 (means twotimes continuously differentiable) and concave. For any parameterψ, it is a non-decreasing function, such that for any 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.
Assumption 5 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 Ξ i ∈ R, which is characterized by a membership function U Ξ i : Ξ i → [0, 1]. That means In order to find a meaningful reference value about the fuzzy numberψ in the fuzzy set W to accurately describe the function, we introduced the D-operation. The definition is as follow Definition 1 [22] There is a function ζ(ψ) : Ξ → R withψ as a parameter, the D-operation of ζ(ψ) is defined as In a sense, the D-operation D[ψ] takes an average value of ζ(ψ) over U Ξ (ψ). In the special case that ζ(ψ) =ψ, this is reduced to the well-known centerof-gravity defuzzification method.

Adaptive robust control design
Although the bound of uncertainty is unknown, the existence of such bound is doubtless. Then we assume that the bound is governed a function named Π(θ a , ϑ a , ψ). The structure of this function is known, but one of the parametersψ is unknown. We estimate the value of ψ throughψ so that the value of Π(θ a , ϑ a , ·) will be infinitely close to the actual value. The accurate bound of uncertainty will be obtained. The estimated valuẽ ψ(t) is used in the controller design, and it satisfies a given adaptive law. That iṡ where k > 0 is constant,

Remark 7
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. The matrix S is the scale coefficient, which determines the response speed of system to improve the regulation accuracy. For S, the appropriate value is generally selected from small to large.
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). Let Then we have the following theorem.
Proof Select a Lyapunov function candidate for the dynamical system as following However the expression of function V is complex and difficult for reading. To make the proving procedures more concise, we define (θ e + Sϑ e ) T M(θ e + Sϑ e ) + ϑ T e β Sϑ e , (40) 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 with Based on Eq. (42), we know λ min (Λ i ) = λ M min , so Moreover, where λ 1 = min{λ M min , k/2}. So, we obtain the possible lower bound of V . Next, we analyze the upper bound. where Based on Eq. (46) 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 = k/2, the information of upper bound is obtained.
To illustrate that V is decreasing, we analyze the derivative of it. The derivative of V with respect to t is given bẏ To make the proving procedures more concise, we define all parts of V related to ϑ e andθ e as V ϑ e , and all parts of V related to ψ as V ψ . TheV ϑ e andV ψ are the derivatives of V ϑ e and V ψ with respect to t. We analyze the two parts of V separately. First, Recalling Eq. (33) and Eq. (38), we geṫ V ϑ e ≤ (θ e + Sϑ e ) T Π(θ a , ϑ a ,ψ) Let λ 3 = min{S 2 ii , 1}, Eq. (52) can be rewritten aṡ The details of the other partV ψ are as followṡ Considering the adaptive law in Eq. (37) and Assumption 5, we havė , ϑ a ,ψ) · θ e + Sϑ e − (ψ −ψ) T kψ.
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 socalled 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. (58) and Eq. (65), we havė The analysis of this differential inequality is explained in detail in Chen [30]. For further analysis of this differential inequality, the following is needed.

Definition 2 [22] Consider a scalar function ω(s(t), t)
with the scalars s(t) and t in some open connected set G, 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. (67) with any t ∈ [t 0 , t 1 ].

Theorem 2 [22] Supose that ϕ(t) is continuous on an open connected set G and there exists a unique solution of the differential equalitẏ
If s(t) and ϕ(t) are the solutions of Eq. (67) and Eq. (68) respectively and the initial condition satisfies Theorem 3 [22] For the differential inequality Eq. (67) and the differential equation Eq. (68), for all (x 1 , t), (x 2 , t) ∈ G, 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. (67) also satisfies the following inequality Exploring the solutions of Eq. (67) 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 nonunique 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 aṡ with ϕ(t 0 ) = V (t 0 ). Then, by solving the differential Eq. (71), we have Therefore, combining Eq. (70) and Eq. (72), we get It should be noted that the initial condition does not necessarily to be t 0 = 0. Because V (t) ≤ ϕ(t) is satisfied at t 0 = 0, so for any t ≥ t 0 , the inequality V (t) ≤ ϕ(t) 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 needed. As mentioned before, the value ofψ is characterized by a membership function. Then 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). Now we define the system performance index function named J as 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. The weighting of the J 1 is normalized to be unity. 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 secondorder partial derivative of J , we further have It should be declared that for any α and β, we have Based on Eq. (80), we know Δ 1 > 0, Δ 3 > 0. So for any solution α l and β l satisfying Eq. (81), J (α l , β l , V 0 ) is the local minimum.
Remark 9 Since for any α and β, Δ 1 > 0 and Δ 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. (81) 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. Figure 3 shows the main design procedure of Fig. 3 The main design procedure of the proposed control scheme 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. (19). 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 .
Step 7: Calculate parameters ρ 0 ∼ ρ 7 , then solve Eq. (81) to get the optimal value α opt and β opt .  (π t), 0.1 cos(π t), 0.1 cos(π t)] T , the actual output of the system is shown in Fig. 4. Under input disturbance τ 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.

Coupling relationship and singularity analysis
We also explore the coupling relationship of the three branches and the singular points in the system. Through the dynamic coupling analysis method in [31], since the inertia matrix M is not a diagonal dominant matrix, there will be coupling torque T ci in the active branch chain i(i = 1, 2, 3) such as The dynamic coupling characteristics of the active branches chain are described according to the inertia matrix. A dynamic coupling index (DC I ), which is used to describe the dynamic coupling degree of each active branch, is given by The DC I of the active branches is defined between [0, 1]. If DC I is closer to 1, it means that the dynamic coupling degree of the corresponding branch is greater. Otherwise, the dynamic coupling degree is smaller. Corresponding to the motion of Fig. 4, the dynamic coupling index of each active branch are all around 0.25. That is, the coupling degree of each branch of the parallel robot is relatively not very great under the selected structure size.
Next, we explore the singular points in the system. From Gosselin's singularity analysis [32] of parallel robots, we know that the relationship between the input and output coordinates could be written as: Differentiating Eq. (84) with respect to t leads to the relationship between the input and output speeds: Simplify Eq. (85) into matrix form as follows: As stated above, singularities occur in configurations where either H or Q becomes singular. Thus, for general closed-loop kinematic chains, a distinction can be made between three kinds of singularities that have different physical interpretations.
(ii) The second kind of singularity occurs when det(H T H ) = 0, then we have the following : (iii) The third kind of singularity occurs when det(Q) = 0 and det(H T H ) = 0 become simultaneously singular. That is We define the range of the actual robot's motion as the workspace. The workspace of the parallel robot is the green area shown in Fig. 5. Through singularity analysis, we know that only the first singularity can occur in the real workspace since the dimensions of the parallel manipulator have been determined. This singularity occurs when the passive link and the active link are on the same line. When the first singularity occurs, the end-effector is at the boundary of the workspace.
As a result, there will be no singularity in the system as long as the end-effector does not move to the boundary.

Verification of proposed control
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 By the membership function Eq. (94) and Doperation in Definition 2, we can obtain the optimal gain α opt and β opt by Eq. (81) 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.

Remark 10
As the ratio μ/ν has a great influence on the optimal solution, we prefer to use the different weight ratio μ/ν directly to 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. How to select this parameter needs to be determined according to different actual needs.
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 system performance index function J with respect to α, β the uncertainty caused by uncertain mass changes and a disturbanceτ to describe the uncertainty caused by disturbance torque. All the simulation results are as follows. From Fig. 6, it can be concluded that the solution of ∂ J ∂α = ∂ J ∂β = 0 is unique. In other words, the optimal gain is also unique. Figure 7 shows the relationship between the system performance index function J and (α, β). Combining with Fig. 7 and Eq. (80), 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.
The proposed method depicts superior performance in the presence of external disturbance, initial position error, and mass change. Comparing with traditional PD control, Fig. 8 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 disturbance torque is added when t = 4s and the detailed results are shown from Figs. 9, 10, 11, 12 and 13. It can be seen that the proposed method is more robust in dealing with external disturbance. Finally, we choose Δm = 0.25 sin(π t) to describe the uncertainty caused by mass variation. 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. From Figs. 9, 10, 11, 12 and 13, 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. By the way, when the ratios μ/ν are the same and the values μ and ν are different, the optimal gains is different. Generally speaking, with the increase of the μ and ν, the corresponding optimal gains is smaller. When the ratios are the same and the values are definitely different, we also make corresponding comparisons which is shown in Fig. 14. The corresponding optimal gains can make the systems outputs track the ideal trajectory, and the tracking accuracy is relatively high. So we often pay more attention to the ratio in the control design.
On the other hand, when the system is affected by external disturbance, the proposed control method can react quickly and maintain stable tracking perfor-  mance. On the contrary, the classical PD control is not strong enough in dealing with external disturbance, and cannot quickly converge to a stable state. From the perspective of tracking error, regardless of the disturbance, the proposed control method has better tracking accuracy. In terms of the ability to deal with 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.
From Figs. 15, 16, 17, 18 and 19, the most obvious difference is that the accumulative errors of the proposed control method are much smaller than the classic PD control. Regardless of the 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 dis-  turbance. However, the classic PD control has no such regularity, so the proposed adaptive robust control is more stable. Generally speaking, the proposed adap- Fig. 19 Accumulative error of ϑ a3 under different control tive 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 two-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 openchain 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.
It is novel and effective to use Udwadia-Kalaba theory to solve the modeling problem of parallel robot. However, to solve the problem of closed-chain constraint force by the Udwadia-Kalaba theory needs a premise, that is, the initial conditions must satisfy the constraint as much as possible. Therefore, the validity of the closed-chain constraint force must be ensured before the controller is designed. In addition, there are many aspects that need to be considered in the proposed control method. In the future, the joint friction and the dynamic model of the servomotor will be considered into the dynamics modeling of parallel robot. The optimal gain design will consider the interrelationship between the gains, rather than just solving the minimum value of the system performance index function.