Abstract
Finding optimal trajectory is critical in several applications of robot manipulators. This paper is applied the open-loop optimal control approach for generating the optimal trajectory of the flexible mobile manipulators in point-to-point motion. This method is based on the Pontryagin's minimum principle that by providing a two-point boundary value problem is solved the problem. This problem is known to be complex in particular when combined motion of the base and manipulator, nonholonomic constraint of the base and highly non-linear and complicated dynamic equations as a result of flexible nature of links are taken into account. The study emphasizes on modeling of the complete optimal control problem by remaining all nonlinear state and costate variables as well as control constraints. In this method, designer can compromise between different objectives by considering the proper penalty matrices and it yields to choose the proper trajectory among the various paths. The effectiveness and capability of the proposed approach are demonstrated through simulation studies. Finally, to verify the proposed method, the simulation results obtained from the model are compared with the results of those available in the literature.
1. Introduction
Mobile manipulators due to their extended workspace offer an efficient application for wide areas. However, these systems are usually “power on board” with limited capacity. Hence, using light links can minimize the inertia and gravity effects on actuators. This helps the mobile manipulator to work in energy - efficient manner and pilots us to use of flexible manipulators.
The platform path planning was investigated for a given end-effector path with assuming that the robot is redundant (Chen M.W. & Zalzala A.M.S., 1997). Then, the robot configuration was determined by minimizing the overall actuation torque and maximizing the manipulability criterion. By employing an appropriate mapping, Papadopoulos et al. transformed the three differentials of the nonholonomic constraint into two differentials and designed the path via polynomial functions (Papadopoulos E. G. & Poulakakis J., 2001). Tanner by implementing a potential field technique using point-world topology proposed a methodology to motion planning for multiple mobile manipulators cooperation that is applicable to articulated, non point nonholonomic robots with guaranteed convergence properties (Tanner H.G., 2003). Sun et al. discussed the load distribution and joint trajectory planning of two robots with one flexible link (Sun Q., 1999). The control of two-link flexible manipulators was studied on the basis of quasi-static equation (Matsuno F. & Hatayama M., 1999). Subudhi and Morris presented a dynamic modeling technique for a manipulator with multiple flexible links and flexible joints based on a combined Euler–Lagrange formulation and assumed modes method. Then, they controlled such system by formulating a singularly perturbed model and used it to design a reduced-order controller (Subudhi B. & Morris A.S., 2002). Path planning of a flexible manipulator was formulated as a discrete-time open-loop optimal control problem which its solution was done via discrete dynamic programming (Wilson D.G. et al., 2004). The assumed mode expansion method was used to derive the dynamic equation of fixed base flexible manipulator (Green A. & Sasiadek J.Z., 2004).
In above-mentioned works, only base mobility or link flexibility has been considered, and the synthesis of the mobile base with flexible links has not been studied. However, as explained before, the flexible links manipulator mounted on the mobile platform can exhibit many advantages. Such systems require less material, consume less power, require smaller motors and they reduce the torque capacity; hence, they can improve the efficiency of the system noticeably. The maximum payload of flexible mobile manipulator was determined along the given trajectory using finite element approach. So, the study did not consider finding the optimal path (Korayem M.H. et al., 2008). A computational algorithm was presented to determine the maximum payload. In this study, linearizing the dynamic equations and constraints are done on the basis of Iterative Linear Programming (ILP) approach for flexible mobile manipulators (Korayem M.H. & Ghariblu H., 2004). But, in ILP method, final boundary conditions are not satisfied exactly and obtained results have a significant error in final time, so the end effector can not place in the desired target point. Furthermore, in this method, the linearizing procedure and its convergence to the proper answer is a challenging issue, especially when nonlinear terms are large and fluctuating, e.g. in problems with consideration of flexibility in links or having high-speed motion. As a result, in the previous ILP based works, the link flexibility has not been considered either in the dynamic equations or in the simulation procedure.
Optimal control can be used in both open loop and close-loop strategies. However, because of the off-line nature of the open loop optimal control in spite of the close-loop ones, many difficulties such as system nonlinearities and all types of constraints may be catered for and implemented easily, so it generally used in analyzing nonlinear systems such as trajectory optimization of different types of robots (Mohri et al., 2001; Furuno S. et al., 2003; Sentinella M. R. & Casalino L., 2006). This method is solved by direct and indirect approaches. However, since direct method leads to the approximate solution also this approach is time consuming and quite ineffective due to the large number of parameters involved (Wilson D.G. et al., 2004; Hull D.G., 1997; Park K. J., 2004), indirect methods is a good candidate for the cases where the system has a large number of degrees of freedom or optimization of the various objectives is targeted (Kirk DE., 1970; Bessonnet G. & Chessé S., 2005; Bertolazzi E. et al., 2005; Callies R. & Rentrop P., 2008). This method is often applied to the generation of highly accurate trajectories and in combination with the structural analysis of optimal control problems in robotics.
Korayem et al. used indirect solution of open-loop optimal control approach to find the maximum load carrying capacity of a rigid link mobile manipulators for a given two-end-point task (Korayem M.H. et al., 2009). In spite of the previous studies based on ILP approach, they satisfied boundary conditions exactly and solved the optimization problem explicitly. Also, the complete form of the obtained nonlinear equation was used in their study and unlike the previous studies linearizing the equations was not required. However, they have studied the rigid-link manipulator and they have not considered the flexibility of links.
In the presenting paper, with applying an indirect solution of optimal control problem, path planning of flexible mobile manipulator is studied based on Pontryagin's minimum principle. The paper firstly deals with the modeling of the general flexible mobile manipulators. Although path planning of flexible manipulators is much more complicated task than the rigid ones, the method handles highly non-linear and complicated dynamics of the system. After that, nonholonomic constraints are defined and dynamic equations are obtained in the state space form. Then, the optimal control problem that with employing of Pontryagin-s minimum principle supports the execution of the optimization solution of model is expressed. In comparison with other method the open-loop optimal control method does not require linearizing the equations, differentiating with respect to joint parameters and using of a fixed-order polynomial as the solution form. Moreover, via changing the penalty matrices values, various optimal trajectories with different specifications can be obtained which able the designer to select a suitable path through a set of obtained paths. Then, an application example with two-link mobile manipulator with flexible links is presented and discussed to demonstrate the effective performance of the proposed approach. Also, in order to validate the method, simulations are performed for the high-stiffness flexible link manipulator and compared with the equivalent rigid links done in 0(Korayem M.H. et al., 2009). The comparison shows a good agreement between the results of this study and those reported in the literature. Lastly, the paper is concluded with highlighting the feature properties of the proposed model.
2. Modeling of a general mobile manipulator with multiple flexible links
In order to path planning of a flexible mobile manipulator, a suitable dynamic modeling of the system is needed. Thus, concepts such as base mobility and link flexibility must be explained exactly. The generalized coordinates of the flexible mobile manipulator consist of three parts, the generalized coordinates defining the mobile base motion →q b = (qb1,qb2,…,q bn b ) T , the generalized coordinates of the rigid body motion of links →q r = (q1, q2,…, q k ) T and the generalized coordinates that related to link flexibility →q f = (q11, q12, …, q21,…, q2n f ,…, q21,…, q2n f ,…, qk1,…, q kn f ) T , where n b and n f are base degree of freedom and number of manipulator mode shapes, respectively. In the mobile manipulators, the concept of redundancy expresses when the number of the overall rigid system degrees of freedom (n = n m + n b ) is strictly greater than the end effector degree of freedom in Cartesian space (m). Where, n m is stood for the rigid manipulator degree of freedom. Thus, the mechanical system is redundant if n > m and the order of redundancy is r = n - m. There are different techniques, which can be applied to a robotic system to solve the redundancy resolution. Some of these techniques are based on an optimization criterion such as overall torque minimization, minimum joint motion and so on. Hence, Seraji has used r additional user-defined kinematic constraint equations as a function of the motion variables (Seraji H., 1998). This method results in simple and online coordination of the control of a mobile manipulator during motion. The presenting study follows this method. Hence, some additional suitable kinematic constraint equations to the system dynamics are applied. Results are in simple and on-line coordination of the mobile manipulator during the motion. These constraints undertake the robot movement only in the direction normal to the axis of the driving wheels along with previously specification of the base trajectory during the motion.
The robotic systems with flexible links are continuous dynamical systems characterized by an infinite number of degrees of freedom and are governed by nonlinear coupled, ordinary and partial differential equations. The exact solution of such systems is not feasible practically and the infinite dimensional model imposes severe constraints on the design of controllers as well. However, they are discretized using assumed modes, finite elements or lumped parameter methods.
In the lumped parameter model, which is the simplest one for analysis purpose, the manipulator is modeled as spring and mass system, which does often not yield sufficiently accurate results (Zhu G. et al. 1999; Kim J.S. & Uchiyama M., 2003). Many authors used the finite element method where the elastic deformations are analyzed by assuming a known rigid body motion and later superposing the elastic deformation with the rigid body motion 0(Korayem M.H. et al., 2008). However, in order to solve a large set of differential equations derived by the finite element method, a lot of boundary conditions have to be considered, which are, in most situations, uncertain for flexible manipulators (Yue S. et al., 2001). In assumed mode model formulation, the link flexibility is usually represented by a truncated finite modal series in terms of spatial mode eigen functions and time-varying mode amplitudes. Hence, using this modes of vibration to model robot dynamics establish an appropriate interaction between flexural vibrations and nonlinear dynamics of system. The assumed mode method and the finite element method use either the Lagrangian formulation or the Newton-Euler recursive formulation to model the problem. In the presenting research, all the flexible links are assumed Euler-Bernoulli beams, so the shearing and rotary inertia effects are neglected. In addition, assumed modes method with the generalized Euler–Lagrange formulation is used to derive the dynamic equations of flexible manipulators.
For a general k-link flexible robot, the vibration v i (x i , t) of each link can be obtained through truncated modal expansion, under the planar small deflection assumption of the link as
where k i is the number of modes used to describe the deflection of link i; φ ij (x i ) and e ij (t) are the j th mode shape function and j th modal displacement for the i th link, respectively.
Hence, applying the Lagrangian assumed modes method the dynamic equation of flexible mobile manipulator could be obtained as:
By defining →q = (q b T , q r T , q f T ) T as the vector of generalized coordinate of the system, Eq. (2) can be written in compact form as
where M is the Inertia matrix, C is the vector of Coriolis and centrifugal forces, G describes the gravity effects and U is the generalized force inserted into the actuator.
3. State-space form of dynamic equation
Consider an n DOFs rigid mobile manipulator with generalized coordinates q = [q i ], i = 1,2,…,n and a task described by m task coordinates r j , j = 1,2,…,m with m < n. By applying h holonomic constraints and c non-holonomic constraints to the system, r = h + c redundant DOFs of the system can be directly determined. Therefore, m DOFs of the system is remained to accomplish the desired task. As a result, we can decomposed the generalized coordinate vector as q = [q r q nr ] T , where q r is the redundant generalized coordinate vector determined by applying constraints and q nr is the non-redundant generalized coordinate vector. Also, considering the flexible link manipulators instead of the rigid ones, their related generalized coordinates, q f , are added to the system; therefore, the overall decomposed generalized coordinate vector of system obtain as q = [q r q nrf ] T , where q nrf is the combination vector of q nr and q f .
In this case, the system dynamics can also be decomposed into two parts: one is corresponding to redundant set of variables, q r and the remained nonholonomic set, q nrf . That is,
where by considering the second row in order to path optimization procedure leads to
Using redundancy resolution q r will be obtained as a known vector in terms of the time (t). Therefore A is obtained as a function of time and q nrf and B as a function of time, q r and q̇ nrf .
By defining the state vector as
Eq. (5) can be rewritten in state space form as
where D = M1 and N= M1(C(X1,X2) + G(X1)). Then, optimal control problem is determined the position and velocity variable X 1 (t) and X 2 (t), and the joint torque U(t) which optimize a well-defined performance measure when the model is given in Eq. (7).
4. Formulation of the optimal control problem
Pontryagin's minimum principle provides an excellent tool to calculate optimal trajectories by deriving a two-point boundary value problem. Let the trajectory generation problem be defined as determining a feasible specification of motion, which will cause the robot to move from a given initial state to a given final state. The method presented in this article adapts in a straightforward manner to the generation of such dynamic profiles.
There are known that nonlinear system dynamics stated as Eq. (7) be expressed in the term of states (X), controls (U) and time (t) as
Generating optimal movements can be achieved by minimizing a variety of quantities involving directly or not some dynamic capacities of the mechanical system. A functional is considered as the integral
where the function L may be specified in quite varied manners. There are initial and terminal constraints on the states:
There may also be certain pragmatic constraints (reflecting such concerns as limited actuator power) on the inputs. For example:
According to the minimum principle of Pontryagin (Kirk D.E., 1970), minimization of performance criterion at Eq. (9), is achieved by minimizing the Hamiltonian (H) which is defined as follow:
where Ψ(t) = [Ψ1(t) T Ψ2(t) T ] T is the non-zero costate time vector-function.
Finally, according to the aforementioned principle, stating the costate vector-equation
in addition to the minimality condition for the Hamiltonian as
leads to transform the problem of optimal control into a non-linear multi-point boundary value problem.
Consequently, for a specified payload value, substituting obtained computed control equations from Eqs. (14) and Eq. (11) into Eqs. (13) and (15), leads to sixteen nonlinear ordinary differential equations which with sixteen
boundary conditions given in Eq. (10), constructs a Two Point Boundary Value Problem (TPBVP). Such a problem is solvable with available commands in different software such as MATLAB and MATEMATHICA.
5. Implementation
In this section, implementation is performed for the flexible mobile manipulator consists of a two flexible link manipulator attached at point F(x f , y f ) on the main axis of a differentially driven vehicle as shown in Figure 1. A concentrated payload of mass m p is connected to the second link.

Two-link mobile manipulator with flexible links
In order to define the mobile manipulator with two flexible links the mechanical system generalized coordinates can be chosen as q = [q b q r q f ] = [x f y f θ0 θ1 θ2 e11 e21], where q b = [x f y f θ0] is the base generalized coordinates vector, q r = [θ1 θ2] is the joint angles vector and q f = [e11 e21] is the vector of link modal displacements. Here, one mode shape per each link is considered to show the simplified model of the system with sensible accuracy. By defining r b , r1, r2 and r mp as the position vectors of the base, first link, second link and payload respectively, according to Fig. 1, the expressions of these vectors can be expressed as:
where v iL is the vibration of the end of each link, θ10 = θ0 + θ1 and θ210 = θ0 + θ1 + θ2. v i (x i , t) and v iL (t) can be expressed as follows:
where with considering the simply support mode shape, φ ij and Ψ can be computed as:
where i and j refer to the number of link and mode shape respectively.
Now, using combined Euler-Lagrange formulation and assumed modes method, the set of system dynamic equations is derived as follows:
where because of the motion is in the horizontal plane, the gravity effects, G = (G b , G r , G f ) will not taking into account. The position of the end effector in the world reference frame RFw can be specified with x e and y e as shown in Figure 1. Therefore operational coordinated of the end effector can be chosen as E(x e , y e ) and the end effector degrees of freedom in the Cartesian coordinate system will be m = 2. The rigid system degree of freedom is equal to n = 5, hence the system has redundancy of order r = n - m = 3. Thus, this system needs three additional kinematical constraints for resolving the redundancy problem and guarantees soft and well-organized movement. There is one nonholonomic constraint for the mobile base that undertakes the robot movement only in the direction normal to the axis of the driving wheels:
Hence, the number of kinematical constraints which must be applied to system for redundancy resolution is decreases to two constraints. In this case, with the previously specification of the base trajectory during the motion, the user-specified additional constraints can be considered. The base position coordinates at point F(x f , y f ), gives x f = X1z and y f = X2z. Where, X1z and X2z are functions in terms of time. Now, by considering fifth order polynomial function for the base trajectory along a straight-line path from (x f (0), y f (0)) to (x f (t f ), y f (t f )) during the overall time t f the redundancy solution of the system can be complete. Velocity at start and stop time is considered to be zero. According to the base motion, ẋ f and ẏ f are known, therefore if the base angle at initial time θ0(0) be specified, angular position and angular velocity of the base (θ0(t), θ̇0(t)), can be determined by solving non-holonomic constraint equations. Note that, the initial base angle is considered to be zero, θ0(0) = 0.
Now, after resolving the redundancy of the system and remaining the non-redundant set of dynamic equations (Eq. (5)), by defining the state vectors as
The state space form of the dynamic equations of the system can be written as
where f2(i) can be obtained from Eq. (7). These equations must be completed with boundary conditions that consist in the specified initial and final states as
Other boundary conditions are assumed to be zero 0(Korayem M.H. et al., 2009). A performance index is expressed as some functional evaluated over the trajectory. In the presenting paper, the desired path is optimally designed so that the required movement can be achieved while minimum effort is obtained and velocities are bounded at proper magnitude. Hence, the objective function is considered as:
The main advantage of this cost functional is that the designer can decide on the relative importance among the objective function. For example, in objective function explained in Eq. (24), between the angular velocity and control effort. In this case, by the numerical choice of w i and r i the desire movement can be achieved. This, results to choose the most appropriate path among various optimal paths that satisfy the designing requirements. Now, by considering the costate vector as Ψ = [Ψ1 Ψ2 … Ψ8], the Hamiltonian function is formed as
where ẋ i , i = 1,…, 8 can be substituted from Eq. (22). Using Eq. (13), differentiating the Hamiltonian function with respect to the states, result in costate equations as follows:
The control function in the admissible interval can be computed according to Eq. (14), by differentiating the Hamiltonian function with respect to the torques (τ1, τ2) and setting the derivative equal to zero. The actuators that are used in this study are the permanent magnet D.C. motor with the final bound of control as
where S i = (τ i / ω mi . τ i and ω mi are the stall torque and maximum no-load speed of i th motor respectively.
Typically, applying the Pontryagin's minimum principle results a two-point boundary value problem. Indeed, substituting the computed control equations in the state and costate equations (22) and (26) yields a 4n ordinary differential equation, where the functions X (t) and Ψ(t) must satisfy the 4n boundary conditions.
Numerical libraries offer efficient computing codes for solving non-linear problems of this type. We used the routine BVP4C of the MATLAB® library, which is based on the collocation method. The details of the numerical technique used in MATLAB® to solve the TPBVP are given in (Shampine L. F. et al).
6. Illustrative examples
In this section, simulations are carried out at two cases. In the first case, path planning is performed for the different values of penalty matrices. In the second case, simulation is done for the high stiffness mobile robot and results are compared with the rigid mobile robot (Korayem M.H. et al., 2009). The necessary parameters used in the simulations are summarized in the Table 1.
Simulation Parameters
The load must be carried from an initial point with coordinate (x e = 0.3m, y e = 0.5m) to the final point with coordinate (x e = 2.45m, y e = 0.5m). The optimal trajectory between these two points during the overall time t f =1.9 second is desired. It should be noted that the final load position is not feasible without the base motion. Simultaneously, the mobile base is initially at point (x f = 0.75m, y f = 0.2m, θ f = 0) and moves to final position (x f = 1.6m, y f = 0.5m). L 0 is supposed to be 40 cm. (Korayem M.H. et al., 2009).
In the first case, the payload is considered to be 2 kg. The purpose is to find the optimal path of payload in such a way that the more suitable amount of control value be applied and the angular velocity values of motors be bounded in ±1.5 rad/s.
By considering penalty matrices as: w1 = w2 = 0.1, w3 = w4 = 0 and r1 = r2 = diag (0.01) the first path with appropriate amount of control value are determined, but the angular velocities are greater than |1.5| rad/s. Hence, reducing the maximum velocity magnitude must be considered. Therefore, for decreasing the joint velocities, at cost functional defined in the Eq. (24) the penalty matrix corresponding to them must be increased. A range of values of W = (w, w, 0, 0) used in simulation are given in Table 2. r i remains without changes.
The values of w used in the simulation
End effector trajectories in the XY plane are shown in Figure 2. Figure 3 shows the angular positions of joints with respect to time. This graph shows that by increasing the w, the angular position change to approach approximately to a straight line. Figure 4 shows the angular velocities of the first and second joints. It can be found that by increasing the w, the final values of angular velocity reduce from −2.3 rad/s to −1.5 rad/s. By growing the w, the angular velocities reduce greatly for the first to second cases whereas at the third case a little reduction has been occurred in spite of great increase in w.

End effector trajectory in XY plane

Angular positions of joints

Angular velocities of joints
The computed torqueses are plotted in Figure 5. As it can be seen, increasing the w causes to raise the torques, so that for the last case the torque curves reach to their bounds at the beginning and end of the path. This result is predictable, because, increasing the w decreases the proportion of r i and the result of this is increasing the control values.

Torques of motors
Mode shapes are plotted in Figure 6, shows the flexible nature of the system in case 3. Finally, arm motions with related end effector and base trajectories in this case are plotted in Figure 7.

Mode shapes of links

robot configuration in the Cartesian plane
It can be concluded that, the first path is the optimal path with the least control values, whereas its angular velocity is the largest magnitude and the last path is the optimal path with the smallest amount of velocities, while its control values is the largest amount. Finally, for the considered problem, the optimal path is the third path, which its velocity magnitude is bounded in ±1.5 rad/s interval and the torque values is the lowest.
In contrast, the optimal path with lower effort with greater velocity and the optimal path with greater effort with lower velocity are obtained using this method; accordingly, based on the objective contrast principle, there is not a unique solution that satisfies all the desired objectives simultaneously. Consequently, in this method, designer compromises between different objectives by considering the proper penalty matrices and obtain the finest path that appropriately cope with the designing requirements.
In order to verify the validity of the study, the case proposed in literature is developed with elastic behavior in links (Korayem M.H. et al., 2009). The rational comparison can be made between rigid manipulator and flexible ones, by increasing the stiffness at the flexible link. It leads to reduce the flexible manners of the system and approaches it into rigid links. Hence, the simulation is done by increasing the moments of area of each link in previous case study from 1e-11 to 6e-9. The payload is considered to be 20.2 kg (Korayem M.H. et al., 2009). Figures 8–10 show the comparison between angular positions, angular velocities and torque curves in the rigid and high stiffness case, respectively; as seen from these figures, both simulation outputs are in a good agreement, which validates the method studied in this paper. Finally, mode shapes in this case is shown in Figure 11.

Angular positions of joints

Angular velocities of joints

Torque curve of actuators

Mode shapes of links 1 and 2
7. Conclusion
In this paper, path planning problem of a flexible mobile manipulator in point-to-point motion is formulated based on the open-loop optimal control approach. The first objective of paper is to state the dynamic optimization problem under a quite generalized form for the optimality solution that commonly encountered in the field of path planning of mechanical manipulators. The second objective consists in developing the method for optimizing the applicable case studies, which results. First, concepts such as redundancy and link flexibility are explained and the general flexible links mobile manipulator is modeled. Then, Pontryagin's minimum principle is called for a Hamiltonian state equations and optimal path with general objective function is designed based on TPBVP. Finally, the method is implemented on two case studies and obtained results are compared with those reported in literature. It illustrates the power and capability of the method to overcome the high nonlinearity nature of the optimization problem in spite of using complete form of the obtained nonlinear equations. Highlighting the main contribution of the paper can be presented as:
The proposed approach can be adapted to any general serial manipulator including both non-redundant and redundant systems with link flexibility and base mobility.
In this approach the nonholonomic constraints do not appear in TPBVP directly, unlike the method given in literature (Mohri A. et al. 2001; Furuno S. et al. 2003).
This approach allows completely nonlinear states and control constraints treated without any simplifications.
This paper is the first path generation study that applied at the flexible mobile manipulators with un-predefined trajectory of end-effector.
The obtained results illustrate the power and efficiency of the method to overcome the high nonlinearity nature of the optimization problem, which with other methods, it may be very difficult or impossible.
In this method, boundary conditions are satisfied exactly, while the results obtained by ILP method have a considerable error in final time.
In this method, designer is able to choose the most appropriate path among various optimal paths by considering the proper penalty matrices.
The optimal trajectory and corresponding input control obtained using this method can be used as a reference signal and feed forward command in the closed-loop control of such manipulators.
