# Embedded Optimal Control of Robot Manipulators with Passive Joints

Authors: Alberto Olivares; Ernesto Staffetti
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper studies the optimal control problem for planar underactuated robot manipulators with two revolute joints and brakes at the unactuated joints in the presence of gravity. The presence of a brake at an unactuated joint gives rise to two operating modes for that joint: free and braked. As a consequence, there exist two operating modes for a robot manipulator with one unactuated joint and four operating modes for a manipulator with two unactuated joints. Since these systems can change dynamics, they can be regarded as switched dynamical systems. The optimal control problem for these systems is solved using the so-called embedding approach. The main advantages of this technique are that assumptions about the number of switches are not necessary, integer or binary variables do not have to be introduced to model switching decisions between modes, and the optimal switching times between modes are not unknowns of the optimal control problem. As a consequence, the resulting problem is a classical continuous optimal control problem. In this way, a general method for the solution of optimal control problems for switched dynamical systems is obtained. It is shown in this paper that it can deal with nonintegrable differential constraints.


---

## Body

## 1. Introduction

In underactuated manipulators the number of available control inputs is strictly less than the number of the degrees of freedom of the robot. However, the control problem for different underactuated manipulators may have different levels of difficulty.This paper studies the optimal control problem for planar underactuated robot manipulators with two revolute joints that move in a vertical plane and, therefore, are subject to the gravity force. The presence of two revolute joints in the mechanical system will be denoted byRR. Both possible models of planar underactuatedRRrobot manipulators under gravity are considered, namely, the model in which only the shoulder joint is actuated, denoted byRR¯, which is usually called Pendubot, and the model in which only the elbow joint is actuated, denoted byR¯R, which is usually called Acrobot. The mechanical system obtained by removing both motors from theRRrobot manipulator under gravity is also considered, which is actually a double pendulum and will be denoted byR¯R¯.Underactuated robots are mechanical systems with second order nonholonomic constraints, because the dynamic equation of the unactuated part of the mechanical system is a second order differential constraint which, in general, is nonintegrable. In the absence of gravity it is not possible to integrate even partially this second order differential constraint in the dynamic model of theRR¯robot manipulator. However, in the presence of this second order nonholonomic constraint the system is controllable. On the contrary, in the dynamic model of theR¯Rrobot manipulator without gravity this differential constraint is completely integrable. It can be converted into a holonomic constraint that makes the system not controllable. As a consequence, the trajectory planning problem has a solution only for particular initial and final states. In the presence of gravity, the dynamic equation of the unactuated part of the mechanical system is a second order nonintegrable differential constraint for both the Pendubot and the Acrobot models and, therefore, both systems are controllable.A comprehensive review of the control properties of underactuated robot manipulators can be found in [1], [2, Chapter 42].In [3] the chaotic features of the double pendulum are studied by means of numerical methods. It is shown that, depending on the initial conditions, this system has a periodic, quasiperiodic, and chaotic behaviour, in which the motion of the system is unpredictable and small differences in initial conditions induce very different trajectories.For both the Pendubot and the Acrobot models, the control objective is, in general, to drive the manipulator away from the stable equilibrium state and steer it at an unstable equilibrium state. Since both theRR¯and theR¯Runderactuated robot manipulators are linearly controllable, it is not difficult to control them around an equilibrium point. However, due to the gravitational drift, the region of the state space where the robot can be kept in equilibrium is reduced and consists of two disjoint manifolds. Moving between these two manifolds requires appropriate swing-up maneuvers, whose synthesis has been tackled by several control [4] and optimal control [5] techniques.The underactuated robot manipulators studied in this paper are supposed to be equipped with brakes at the unactuated joints. The Pendubot and the Acrobot with brakes at the unactuated joints will be denoted byRR¯bandR¯bR, respectively, whereas theR¯R¯mechanical system with brakes at both unactuated joints will be denoted byR¯bR¯b. This mechanical system can be regarded as a brake-actuated mechanical system and will be referred to as Brakebot. Thus, since each passive joint has two operating modes, unbraked and braked, there are two operating modes for the Acrobot and the Pendubot and four operating modes for the Brakebot.Since the Brakebot is a passive system it can only move from higher to lower energy states and, therefore, the control objective is, in general, to drive the manipulator away from an unstable equilibrium state and steer it at a stable equilibrium state with appropriate swing-down maneuvers eliminating periodic, quasiperiodic, or chaotic behaviour. For the synthesis of these maneuvers, optimal control techniques can be employed. The cost functional to be optimized could be, for instance, the duration of the maneuver, the number of switches of the dynamical system between unbraked and braked modes, or the total duration of the evolution of the system in braked mode.The problem studied in this paper can be stated as follows: given an initial state and a final state, find the sequences of modes, the corresponding trajectory, and available control inputs that satisfy the dynamic equation of the robot manipulator and steer the system between the initial and the final states optimizing a cost functional during the motion. Since the optimal sequence of modes has to be determined, this problem is actually an optimal control problem of a switched dynamical system.Switched systems are particular hybrid systems [6] whose continuous state does not exhibit jumps at the switching instants. Solving the optimal control problem for a switched system implies finding the optimal sequence of switching instants, the optimal sequence of discrete modes, and the optimal value for the continuous control input.The optimal control problem of switched dynamical systems can be formulated using integer or binary variables to model the choice of the optimal mode sequence. Optimal control problems involving binary or integer variables are called Mixed-Integer Optimal Control (Mioc) problems [7]. In [8], recent progress in theoretical bounds, reformulations, and algorithms for this class of problems have been reviewed and a comprehensive algorithm, based on the solution of a sequence of purely continuous problems and simulations, has been proposed. A hybrid optimal control technique based on Mioc has been presented in [9] and applied to several mechanical systems including a Pendubot with a brake at the unactuated joint.Hybrid control techniques have been applied to both the Pendubot and the Acrobat models in [10] and [11], respectively.In [12], the optimal control problem of switched systems has been solved via embedding into a continuous optimal control problem. It has been shown that, for quite a general class of optimal control problems for switched systems, the computational complexity of the Embedded Optimal Control (Eoc) method is no greater than that of smooth optimal control problems. With this technique the switched system is embedded into a larger set of systems and the optimal control problem is formulated in the latter. The main advantages of this method are that assumptions about the number of switches are not necessary, integer or binary variables do not have to be introduced to model switching decisions between modes, and the optimal switching times between modes are not unknowns of the optimal control problem. As a consequence, the resulting problem is a classical continuous optimal control problem. This greatly reduces the computation time to find a solution with respect to the Mioc technique. In [13] the Eoc approach introduced in [12] has been extended to incorporate hybrid behaviour due to autonomous switches and to systems with an arbitrary number of modes. In [14] the Eoc method has been compared with other optimal control methods for switched systems over several dynamic models with autonomous and controlled switches. A comprehensive treatment of all the theoretical results obtained with the Eoc approach to the optimal control problem of switched mechanical systems can be found in [15, Chapter 31].In this paper, the optimal control of the Pendubot, Acrobot, and Brakebot has been studied using the Eoc approach. The obtained results have been compared with those obtained using a Multiphase Mixed-Integer Optimal Control (Mmioc) method. Due to the nature of the problem, only controlled switches have been considered. The results have shown that the Eoc approach is effective in solving optimal control problems for these underactuated mechanical systems, which are highly nonlinear and have nonholonomic constraints.This paper is organized as follows. In Section2 the dynamic models of the Pendubot, Acrobot, and Brakebot are described and in Section 3 their control properties are discussed. The optimal control problem for these dynamic systems is stated in Section 4 and in Section 5 the Eoc approach is introduced. In Section 6 the results of the application of this method to several instances of optimal control problem are reported. Finally, Section 7 contains the conclusions.


## 2. Dynamic Model of Underactuated Manipulators

The general dynamic model of a robot manipulator is described by the following second order differential equation:(1)Bθθ¨+Cθ,θ˙θ˙+h+eθ=u-Gθw,where the first term of this equation,B(θ)θ¨, represents the inertial forces due to acceleration at the joints and the second term,C(θ,θ˙)θ˙, represents the Coriolis and centrifugal forces. The third term represents both the effects of friction at the joints and the effects of friction induced by electromagnetic friction brakes at the unactuated joints. The terme(θ)represents the potential forces such as elasticity and gravity. The matrixG(θ)on the right-hand side maps the external forces/torqueswto forces/torques at the joints. Finally,urepresent the forces/torques at the joints, which are the control variables of the system. Explicit time dependence is omitted in this section to simplify the notation.The links are supposed to be rigid, as well as the transmission elements. No external forces are supposed to be acting on the mechanical system. Under these hypotheses, the dynamic model of the robotic system reduces to(2)Bθθ¨+Cθ,θ˙θ˙+h+eθ=u.Only the effects of friction induced by the electromagnetic friction brakes at the unactuated joints will be taken into account. This means thathrepresents the braking torques induced by electromagnetic friction brakes at the unactuated joints.A planar verticalRRmanipulator is composed of two homogeneous links and two revolute joints moving in a vertical plane{x,y}, as shown in Figure 1, whereliis the length of linki,riis the distance between jointiand the mass center of linki,miis the mass of linki, andIziis the barycentric inertia with respect to a vertical axis,z, of linki, fori=1,2. In this case the two matricesB(θ)andC(θ,θ˙)and the vectore(θ)have the form(3)Bθ=α+2βcos⁡θ2δ+βcos⁡θ2δ+βcos⁡θ2δ,(4)Cθ,θ˙=-βsin⁡θ2θ˙2-βsin⁡θ2θ˙1+θ˙2βsin⁡θ2θ˙10,(5)eθ=m1r1+m2l1gcos⁡θ1+m2r2gcos⁡θ1+θ2m2r2gcos⁡θ1+θ2, whereθ=(θ1,θ2)Tis the vector of configuration variables,θ1being the angular position of link 1 with respect to the x-axis of the reference frame{x,y}andθ2the angular position of link 2 with respect to link 1 as illustrated in Figure 1. The vectorθ˙=(θ˙1,θ˙2)Tis the vector of angular velocities andθ¨=(θ¨1,θ¨2)Tis the vector of angular accelerations. The vectoru=(u1,u2)Tis the vector of control inputs, in whichu1is the torque applied by the actuator at Joint1, andu2is the torque applied by the actuator at Joint2. Constantg=9.80665 [m/s2] is the gravity acceleration. The parametersα,β, andδin (3) and (4) have the following expressionsα=Iz1+Iz2+m1r12+m2(l12+r22),β=m2l1r2, andδ=Iz2+m2r22.Figure 1
A planarRRrobot manipulator that moves in a vertical plane. In the Pendubot model Joint 2 is not actuated whereas in the Acrobot model Joint 1 is not actuated. Both underactuated planarRRrobot manipulator models are supposed to be equipped with brakes at the unactuated joints. In the Brakebot model both joints are not actuated and equipped with brakes.The dynamic model (2) can be rewritten in the following more explicit form:(6)B11θB12θB21θB22θθ¨1θ¨2+C11θ,θ˙C12θ,θ˙C21θ,θ˙C22θ,θ˙θ˙1θ˙2+h1h2+e1θe2θ=u1u2.A robot manipulator is said to be underactuated when the number of actuators is less than the degree of freedom of the mechanical system. Thus, the dynamic model (6) can be rewritten for RR¯robot manipulator with a brake at the unactuated joint as follows:(7)B11θB12θB21θB22θθ¨1θ¨2+C11θ,θ˙C12θ,θ˙C21θ,θ˙C22θ,θ˙θ˙1θ˙2+0h2+e1θe2θ=u10,and for R¯Rrobot manipulator with a brake at the unactuated joint as follows:(8)B11θB12θB21θB22θθ¨1θ¨2+C11θ,θ˙C12θ,θ˙C21θ,θ˙C22θ,θ˙θ˙1θ˙2+h10+e1θe2θ=0u2.In passive joints equipped with brakes, the braking torquehi, i∈1,2, can be described as a friction-velocity map with load dependency and lockup behavior [16]:(9)hi=dθ˙i+c+qfcl·sgn⁡θ˙iifθ˙i≥ɛ,heiifθ˙i<ɛ,hei<hs+qfcl,hs+qfcl·sgn⁡heiotherwise.Heredis the viscous friction coefficient,cis the load-independent Coulomb friction torque,qis the friction load dependency, andheiis the external torque at jointi.fclis the brake clamp force, andhsis the load-independent static friction torque, which are supposed to have the same values at the two brakes. The parameterɛdefines a small zero velocity bound in accordance with Karnopp’s model for zero velocity detection. A continuously differentiable version of this friction model [17] has been used for the experiments.Equations (7) and (8) withh2=0and withh1=0represent the unbraked operational modes of the Pendubot and the Acrobot, respectively. In this case, the second equation in (7) and the first equation in (8) describe the dynamics of the unactuated part of the mechanical system, which are second order differential constraints without input variables. Activating brakes at the passive joints of the Pendubot and the Acrobot gives rise to the braked operational modes of these two mechanical systems which thus have two operational modes, namely, (1) unactuated joint free and (2) unactuated joint braked.The Brakebot has the following dynamic model:(10)B11θB12θB21θB22θθ¨1θ¨2+C11θ,θ˙C12θ,θ˙C21θ,θ˙C22θ,θ˙θ˙1θ˙2+h1h2+e1θe2θ=00.Since all the joints of the Brakebot are unactuated, it has four operational modes depending on which joints are braked and which are free, namely, (1) both joints free, (2) the first joint free and the second joint braked, (3) the first joint braked and the second joint free, and (4) both joints braked. The corresponding dynamic equations can be obtained from (10) by setting one or both braking torques to zero.


## 3. Control Properties of UnderactuatedRR Robot Manipulators

In this section, the main control properties of planar underactuatedRRrobot manipulators in the presence of gravity will be described. For a more general description of the control properties of underactuated robot manipulators see [1].Optimal control approaches to trajectory planning assume that there exists a control input that steers the system between two given states. Thus, controllability is the most important aspect to be checked before studying the optimal control problem of a dynamic system. If in the trajectory planning problem the duration of the motionTis not assigned, the existence of a finite-time solution for any state(θF,θ˙F)in a neighborhood of(θI,θ˙I)is equivalent for the robotic system to the property of local controllability at(θI,θ˙I). If local controllability holds at any state, then the system is controllable and the trajectory planning problem is solvable for any pair of initial and final states.For underactuatedRRrobot manipulators, controllability is related to integrability of the second order nonholonomic constraints. The second order differential constraint in (7) and (8) may either be partially integrable to a first order differential equation or completely integrable to a holonomic equation. Necessary and sufficient integrability conditions are given in [18]. If the second order differential constraint in (7) or (8) is not partially integrable, it is possible to steer the corresponding system between equilibrium points. This occurs for planar underactuatedRR¯robot manipulators in the absence of gravity and for both theRR¯andR¯Rin the presence of gravity. If the second order differential constraint in (7) or (8) is completely integrable to a holonomic constraint, the motion of the corresponding mechanical system is restricted to a1-dimensional submanifold of the configuration space which depends on the initial configuration. This occurs for the planar underactuatedR¯Rrobot manipulators without the effects of gravity [18]. For this robot model, the trajectory planning problem has solution only for particular initial and final states.Thus, when the second order differential constraint in (7) or (8) is not partially or completely integrable, the corresponding mechanical system is controllable. However, several aspects of controllability can be studied which characterize planar underactuatedRRrobot manipulators.The equilibrium conditions for theRR¯robot manipulator are (11)e1θ1,θ2=gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2=u1,e2θ1,θ2=gm2r2cos⁡θ1+θ2=0, for theR¯Rrobot manipulator are (12)e1θ1,θ2=gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2=0,e2θ1,θ2=gm2r2cos⁡θ1+θ2=u2, and, finally, for theR¯R¯robot manipulator are (13)e1θ1,θ2=gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2=0,e2θ1,θ2=gm2r2cos⁡θ1+θ2=0. Not all the solutions of these equations correspond to stable equilibrium configurations. Foru1=0andu2=0these dynamical systems have three unstable configurations; namely,(θ1,θ2)=(π/2,0), which will be called Up Configuration,(θ1,θ2)=(-π/2,π)and(θ1,θ2)=(π/2,-π), which will be called Mid Configurations, and one stable equilibrium configuration,(θ1,θ2)=(-π/2,0), which will be called Down Configuration.Both planar underactuatedRR¯andR¯Rrobot manipulators in the presence of gravity are linearly controllable [1]. Linear controllability means that the linear approximation of the system around an equilibrium point is controllable. Both planar underactuatedRR¯andR¯Rrobot manipulators in the presence of gravity are neither small-time locally controllable nor small-time locally configuration controllable [19, 20]. Finally, both planar underactuatedRR¯andR¯Rrobot manipulators in the presence of gravity are not kinematically controllable [19–21]. The lack of kinematic controllability implies that in general a configuration is not reachable by means of a sequence of kinematic motions, that is, feasible paths in the configuration space which may be followed with an arbitrary timing law.


## 4. The Optimal Control Problem

In this section, following [13], the optimal control problem for switched dynamical systems will be described.In switched dynamical systems, the vector fields that describe the evolution of the system undergo discontinuous jumps depending on the state and the input. These switches are called autonomous. Discontinuous jumps in the evolution of switched dynamical systems can also be induced by a control input. For this reason these switches are called controlled.The evolution of a switched dynamical system subject to autonomous and controlled switches can be described by the following four quantities: (1) the discrete state or modeξ(t)∈Dξ={1,2,…,dξ}, (2) the continuous statex(t)∈Rn, (3) the switching control inputv(t)∈Dv={1,2,…,dv}, which can be regarded as a discrete control input, and (4) the continuous control inputu(t)∈Rm. In this work, we only consider memoryless systems for which the autonomous switches depend on the continuous statex(t)and the inputu(t), and not on the current discrete stateξ(t). The evolution of the discrete state of the memoryless system is described by a piecewise continuous functionη:Rn×Rm→Dξ, which selects the discrete stateξof the system based on the continuous statex(t)and the continuous control inputu(t); that is,(14)ξ+t=ηx,u.Letf(i,j), i∈Dξ, j∈Dv, be a collection ofC1vector fieldsf(i,j):Mi→Rn, whereMi⊆Rn×Rm, i∈Dξ, is the set of pairs(x,u)corresponding to the discrete statei∈Dξ; that is,Mi={(x,u)∈Rn×Rm∣η(x,u)=i}. The evolution of the continuous statex(t)is described by(15)x˙=fηxt,ut,vtxt,ut,xt0=x0.This means that, at eacht≥t0and for each discrete stateξ(t)∈Dξ, the switching control inputv(t)∈Dvselects the vector field that governs the evolution of the continuous state. Note that in (15) the switching control inputv(t)does not affect the autonomous switches. Systems described by (15) are called systems with decoupled switches. The continuous control inputu(t)is assumed to belong to a convex and compact setΩ⊆Rm, andu(t)andv(t)are supposed to be measurable functions.Since (14) completely determines the discrete stateξ(t)based onx(t)andu(t), a piecewiseC1vector fieldfj(x(t),u(t))=f(η(x(t),u(t)),j)(x(t),u(t))can be defined for eachj∈Dvand (15) can be rewritten as follows:(16)x˙=fvtxt,ut,xt0=x0.We are interested in computing optimal control laws for systems described by (16). This entails finding (1) the optimal sequence of switching control inputs, that is, the optimal sequence of values of the functionv(t), (2) the optimal switching instants,t1,t2,…,tn, including their number,n, and (3) the optimal continuous control over each interval[ti,ti+1), i=0,1,2,…,n.Mmioc techniques can be employed to find the optimal solution of this problem. However, the computational complexity of the resulting problem limits the application of these techniques in practice. It has been shown in [12] that for quite a general class of optimal control problems of switched dynamical systems, a technique called embedding can be applied, whose computational complexity is no greater than that of continuous optimal control problems. This technique will be described in Section 5.Consider the system described by (16), wherev(t)andu(t)are control variables of the optimal control problem. The cost functional of the optimal control problem, defined over the interval[t0,tf], can be expressed in the following form:(17)JCt0,x0,u,v=gt0,x0,tf,xf+∫t0tfFηxt,ut,vt0xt,utdt,where(t0,x(t0))∈T0×B0,(tf,x(tf))∈Tf×Bf, andT0,B0,Tf, andBfare the boundary sets fort0,x(t0),tf,x(tf).gis a real-valued function defined in the neighbourhood B=T0×B0×Tf×Bf, which is assumed to be a compact set, andF(i,j)0:Mi→R, i∈Dξ, j∈Dv, areC1functions.By (14), which determines the discrete stateξ(t)based onx(t)andu(t), we can define a new piecewiseC1functionFj0(x(t),u(t))=F(η(x(t),u(t)),j)0(x(t),u(t))and rewrite the cost functional as follows:(18)JCt0,x0,u,v=gt0,x0,tf,xf+∫t0tfFvt0xt,utdt.Thus, the optimal control problem for the above defined class of switched systems, called Switched Optimal Control Problem (Socp), can be stated as follows: (19)minu,v⁡JCt0,x0,u,v subject to (1) (16), (2) boundary conditions(t0,x(t0),tf,x(tf))∈B, and (3)  u(t)∈Ω⊆Rm, andv(t)∈Dvfort∈[t0,tf].


## 5. Embedding Approach

The approach described in [12] to solve optimal control problems of switched dynamical systems consists in embedding system (16) into a larger set of systems. Sincev(t)∈Dv={1,2,…,dv}definedvnew variablesvi∈[0,1], i∈Dv, such that(20)∑i=1dvvit=1.Letuibe the control input for each vector fieldfi,i∈Dvin (15). Define a new system(21)x˙t=∑i=1dvvitfixt,uit,xt0=x0and the associated cost functional(22)JEt0,x0,u,v=gt0,x0,tf,xf+∫t0tf∑i=1dvvitFi0xt,uitdt.The Socp is converted into the following Embedded Optimal Control Problem (Eocp): minimize the functional (22) over all functionsvianduisubject to the following constraints: (1) (21), (2) boundary conditions(t0,x(t0),tf,x(tf))∈B, (3) for eacht∈[t0,tf]and eachi∈Dv,vi(t)∈[0,1]andui(t)∈Ω, and (4) for eacht∈[t0,tf],∑i=1dv‍vi(t)=1.This is an optimal control problem without integer or binary variables. Therefore, classical necessary and sufficient conditions of optimal control theory can be applied to solve it. However, to guarantee the existence of a solution some additional hypotheses are needed. We assume that the functionsfi(x(t),ui(t))are affine in the continuous control variable, that is,fi(x(t),ui(t))=Ai(x)+Bi(x)ui, and that the functionsFi0(x(t),ui(t))are convex inuifor each(t,x).It has been shown that, once a solution of the Eocp has been obtained, either the solution is of the switched type, that is, exactly one of thevi’s being 1and all the others being0, or suboptimal trajectories of the Socp can be constructed that can approach the value of the cost for the Eocp arbitrarily closely and satisfy the boundary conditions withinϵ, with arbitraryϵ>0. A thorough discussion about the relationship between Socp’s and Eocp’s solutions can be found in [12]. The numerical method employed to transcribe the Eocp into Nlp problem is described in [22].


## 6. Implementation and Results

Several numerical experiments have been carried out for theRR¯b,R¯bR, andR¯bR¯bplanar underactuated robot manipulators in the presence of gravity. The corresponding dynamical systems are described by vector fields that are affine in the continuous control variables, thus fulfilling the sufficient conditions for the existence of a solution stated in the previous section.The results obtained using the Eoc approach have been compared to those obtained using a Mmioc technique [7]. In the Mmioc approach, to find the optimal sequence of modes an iterative process has been used. First, the maximum number of switches between modes has been estimated. This number determines the maximum number of phases that are considered. Then, an optimal control problem has been formulated, in which the unknown switching times are part of the state, and binary variables have been introduced to represent the choice of the active mode in each phase. Finally, a third-degree Gauss-Lobatto direct collocation method has been used to tackle the dynamic equations. The resulting Minlp problem has been solved using a nonlinear programming based branch-and-bound algorithm specifically tailored to the problem. Details about the implementation of this method can be found in [23]. As pointed out in [8], special care has to be taken to treat the cases in which phase lengths diminish during the optimization procedure. Therefore, the phases of zero duration have been eliminated.


### 6.1. Pendubot with Brake

In this section, the optimal control problem for a Pendubot with a brake at the unactuated joint is studied. In this robot model, the second joint is not actuated; thusu=(u1,0)T. In this case, it is not possible to integrate partially or completely the nonholonomic constraint because the conditions for partial integrability [18] are not fulfilled. Hence, the system is controllable. The numerical results of the application of the Eoc method for optimal control to two different boundary value problems for this system will be described. In both cases the functional to be minimized is the energy consumption during the motion.The Pendubot has two operational modes, depending on whether the unactuated joint is braked or not, namely,(i)
Mode 1: Joint 2 free,(ii)
Mode 2: Joint 2 braked. From (7), the dynamic model of a Pendubot with a brake at the unactuated joint is (23)α+2βcos⁡θ2θ¨1+δ+βcos⁡θ2θ¨2-βsin⁡θ22θ˙1θ˙2+θ˙22+gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2=u1,δ+βcos⁡θ2θ¨1+δθ¨2+βsin⁡θ2θ˙12+gm2r2cos⁡θ1+θ2+h2=0. In these equationsu1is the first component of the input vectoru=(u1,u2)T. When the Pendubot is operating in Mode 1, that is, with Joint 2 free,h2=0. Otherwise, it takes the form (9). To convert this second order differential model into the first order model (16) the following change of variables is introduced:(24)x1=θ1,x2=θ2,x3=θ˙1,x4=θ˙2,with the following additional differential relations:(25)x˙1=x3,x˙2=x4,among the components of state variablex=(x1,x2,x3,x4)T. The result is (26)α+2βcos⁡x2x˙3+δ+βcos⁡x2x˙4-βsin⁡x22x3x4+x42+gl1m2+m1r1cos⁡x1+gm2r2cos⁡x1+x2=u1,δ+βcos⁡x2x˙3+δx˙4+βsin⁡x2x32+gm2r2cos⁡x1+x2+h2=0. Since the Pendubot has two dynamical modes,Dv={1,2},2vector fields,f1=(f1,1,f1,2,f1,3,f1,4)Tandf2=(f2,1,f2,2,f2,3,f2,4)T, must be considered. The dynamic models of the Pendubot operating in Mode 1 and Mode 2 are the following: (27)x˙1=x3=f1,1,x˙2=x4=f1,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2gm2r2cos⁡x1+x2+βx32sin⁡x2+δu1,1+βx3+x42sin⁡x2α-δδ-β2cos⁡x22=f1,3,x˙4=u1,1δ+βcos⁡x2-gl1m2+m1r1cos⁡x1δ+βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+βx42cos⁡x2sin⁡x2+βx3x3+x4sin⁡2x2·δ-α+δ+β2cos⁡x22-1=f1,4,x˙1=x3=f2,1,x˙2=x4=f2,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δh2+u2,1+βx3+x42sin⁡x2α-δδ-β2cos⁡x22=f2,3,x˙4=-h2+gm2r2cos⁡x1+x2+βx32sin⁡x2δ+δ+βcos⁡x2-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δh2+u2,1+βx3+x42sin⁡x2·δδ-α+δ+β2cos⁡x22-1=f2,4.In these equations,u1,1is the first component of vectoru1=(u1,1,u1,2)T, the input vector for vector fieldf1, andu2,1is the first component of vectoru2=(u2,1,u2,2)T, the input vector for vector fieldf2. Moreover, two new variables,vi∈[0,1], i∈{1,2}, such thatv1+v2=1, must be introduced. Thus, (21) takes the form (28)x˙1=v1f1,1+v2f2,1,x˙2=v1f1,2+v2f2,2,x˙3=v1f1,3+v2f2,3,x˙4=v1f1,4+v2f2,4, and the cost functional (22) becomes (29)JE=∫t0tfv1u1,12+v2u2,12dt.This objective functional fulfills the sufficient conditions for the existence of a solution stated in Section 5.Inequality constraints on the state variables have been also considered. In particular, variableθ2has been bounded to take into account possible range limitations of the mechanical system. Thus, the inequality constraints-π[rad]≤θ2≤π[rad]have been introduced in the model. Inequality constraints on the available control variable have also been considered. In particular, the constraints-25[Nm]≤u1≤25[Nm]have been introduced in the model to take into account physical limitations of the actuators.The following parameters of the dynamic model have been used:Iz1=Iz2=1[kgm2],l1=l2=1[m], andm1=m2=1[kg].


#### 6.1.1. Swing to the Up Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(30)θ1tI=-π2rad,θ1tF=π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.6 [s], respectively.The energy consumption, the resolution time, the number of phases, the sequence of modes, and the corresponding switching times between them obtained with both the Eoc and Mmioc methods are reported in Table 1. Figure 2 shows the sequence of configurations of the robot at timest=k3.6/64, k=0,1,…,64, corresponding to the optimal solution of the boundary value problem obtained with a discretization of the interval[tI,tF]into128subintervals. Configurations that correspond to the braked mode of the mechanical systems are represented in red. It is easy to see that the braked mode is used in the solution just before swinging and when the system is reaching its final state. Figure 3 shows the corresponding control and state variables obtained with the Eoc method (in red) and the Mmioc approach (in black). In the Eoc solution, six phases can be distinguished, although two of them have a very short duration. The optimal sequence of modes obtained has been(0,1,0,1,0,1)which corresponds to braking the unactuated joint in phases2,4, and6and not braking in phases1,3, and5. In the Mmioc solution,4phases can be distinguished. The optimal sequence of modes obtained has been(0,1,0,1)which corresponds to braking the unactuated joint in phases2and4and not braking in phases1and3. It is easy to see from Figure 3(f) that the optimal sequence of modes obtained with the Eoc method coincides with that obtained with the Mmioc approach except for the fifth short phase that does not exist in the Mmioc solution. As a consequence, the continuous solutions almost coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been192.71 [J], whereas the value obtained with the Mmioc approach has been195.532 [J], higher than the previous one. The computational time to find a solution with the Eoc method has been2.84221 [s] whereas it has been26.984307 [s] with the Mmioc approach, higher than the previous one.Table 1
Results for the swing-up maneuver for the Pendubot with brake:tf=3.6 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
192.71
2.84221
(0, 1, 0, 1, 0, 1)
t
1
=
1.01250
t
2
=
2.08125
t
3
=
3.29062
t
4
=
3.54375
t
5
=
3.57187
Mmioc
195.532
26.984307
(0, 1, 0, 1)
t
1
=
0.98867
t
2
=
2.07424
t
3
=
3.29758Figure 2
Sequence of configurations of the optimal solution for the swing-up maneuver for the Pendubot with brake. The configurations in which the Pendubot operates in Mode 1 and Mode 2 are represented in black and red, respectively. The corresponding control and state variables are represented in Figure3.Figure 3
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-up maneuver of the Pendubot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 2. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach. The dashed vertical lines represent switches between phases.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


#### 6.1.2. Swing to the Down Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(31)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.6 [s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 2, Figures 4 and 5 with the same interpretation as in Section 6.1.1. It is easy to see from Figure 5(f) that the optimal sequence of modes obtained with the Eoc method coincides with that obtained with the Mmioc approach but the time location of the corresponding phases does not. As a consequence, the continuous solutions do not coincide. Nevertheless, the value of the objective functional for the optimal solution obtained with the Eoc method has been21.5397[J]whereas the value obtained with the Mmioc approach has been21.5739[J], very similar to the previous one. The computational time to find a solution with the Eoc method has been1.57952[s], shorter than computational time to find a solution with the Mmioc approach, which has been28.65894[s].Table 2
Results for the swing-down maneuver for the Pendubot with brake:tf=3.6 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
21.5397
1.57952
(0, 1, 0, 1, 0, 1)
t
1
=
1.01250
t
2
=
1.37813
t
3
=
2.25000
t
4
=
2.67188
t
5
=
3.26250
Mmioc
21.5739
28.65894
(0, 1, 0, 1, 0, 1)
t
1
=
1.11980
t
2
=
1.53434
t
3
=
2.41782
t
4
=
2.81834
t
5
=
3.31516Figure 4
Sequence of configurations of the optimal solution for the swing-down maneuver for the Pendubot with brake. The corresponding control and state variables are represented in Figure5.Figure 5
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Pendubot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 4. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


### 6.2. Acrobot with Brake

In this section, the optimal control problem of an Acrobot with a brake at the unactuated joint is studied. In this robot model the first joint is not actuated; thusu=(0,u2)T. In this case, it is not possible to integrate partially or completely the nonholonomic constraint because the conditions for partial integrability [18] are not fulfilled. Hence, the system is controllable. The numerical results of the application of our method for optimal control to two different boundary value problems for this system will be described. In both cases the functional to be minimized is the energy consumption during the motion.The Acrobot has two operational modes, depending on whether the unactuated joint is braked or not, namely,(i)
Mode 1: Joint 1 free,(ii)
Mode 2: Joint 1 braked. From (8), the dynamic model of an Acrobot with a brake at the unactuated joint is (32)α+2βcos⁡θ2θ¨1+δ+βcos⁡θ2θ¨2-βsin⁡θ22θ˙1θ˙2+θ˙22+gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2+h1=0,δ+βcos⁡θ2θ¨1+δθ¨2+βsin⁡θ2θ˙12+gm2r2cos⁡θ1+θ2=u2. In these equationsu2is the second component of the input vectoru=(u1,u2)T. When the Acrobot is operating in Mode 1, that is, with Joint 1 free,h1=0. Otherwise, it takes the form (9). To convert this second order differential model into the first order model (16) the change of variables (24) and the additional differential relations (25) are introduced. The result is (33)α+2βcos⁡x2x˙3+δ+βcos⁡x2x˙4-βsin⁡x22x3x4+x42+gl1m2+m1r1cos⁡x1+gm2r2cos⁡x1+x2+h1=0,δ+βcos⁡x2x˙3+δx˙4+βsin⁡x2x32+gm2r2cos⁡x1+x2=u2. Since the Acrobot has two dynamical modes,Dv={1,2}, two vector fields,f1andf2, must be considered. The dynamic models of the Acrobot operating in Mode 1 and Mode 2 are the following: (34)x˙1=x3=f1,1,x˙2=x4=f1,2,x˙3=-δu1,2-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2-u1,2+gm2r2cos⁡x1+x2+βx32sin⁡x2α-δδ-β2cos⁡x22=f1,3,x˙4=-gl1m2+m1r1cos⁡x1δ+βcos⁡x2-u1,2α+2βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+β2x32+2x3x4+x42cos⁡x2sinx2·δ-α+δ+β2cos⁡x22-1=f1,4,x˙1=x3=f2,1,x˙2=x4=f2,2,x˙3=-δh1+u2,2-δgl1m2+m1r1cos⁡x1+βδx3+x42sinx2+βcos⁡x2-x6+gm2r2cosx1+x2+βx32sinx2α-δδ-β2cos⁡x22=f2,3,x˙4=1δu2,2-gm2r2cos⁡x1+x2-βx32sin⁡x2+δ+βcos⁡x2-δh1+u2,2-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2-u2,2+gm2r2cos⁡x1+x2+βx32sin⁡x2·δ-α+δ+β2cos⁡x22-1=f2,4.In these equations,u1,2is the second component of vectoru1=(u1,1,u1,2)T, the input vector for vector fieldf1, andu2,2is the second component of vectoru2=(u2,1,u2,2)T, the input vector for vector fieldf2. Moreover, two new variables,vi∈[0,1], i∈{1,2}, such thatv1+v2=1, must be introduced. Thus, (21) takes the form (35)x˙1=v1f1,1+v2f2,1,x˙2=v1f1,2+v2f2,2,x˙3=v1f1,3+v2f2,3,x˙4=v1f1,4+v2f2,4, and the cost functional (22) becomes (36)JE=∫t0tfv1u1,22+v2u2,22dt. This objective functional fulfills the sufficient conditions for the existence of a solution stated in Section 5.Variableθ2has been bounded, to take into account possible range limitations of the mechanical system. Thus, the inequality constraints-π[rad]≤θ2≤π[rad]have been introduced in the model together with the constraints-25[Nm]≤u2≤25[Nm]to take into account physical limitations of the actuators. The same parameters of the dynamic model of Section 6.1 have been used.


#### 6.2.1. Swing to the Up Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(37)θ1tI=-π2rad,θ1tF=π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.4[s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 3, Figures 6 and 7 with the same interpretation as in the previous sections. It is easy to see from Figure 7(f) that the optimal sequence of modes obtained with Eoc method does not coincide with that obtained with the Mmioc approach, since the first unbraked phase does not exist in the Mmioc solution. As a consequence, the continuous solutions do not coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been113.968[J]whereas the value obtained with the Mmioc approach has been113.027[J], a bit lower than the previous one. The computational time to find a solution with the Eoc method has been20.2848[s], whereas it has been32.137628[s]with the Mmioc approach, higher than the previous one.Table 3
Results for the swing-up maneuver for the Acrobot with brake:tf=3.4 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
113.027
20.2848
(0, 1, 0, 1, 0, 1)
t
1
=
0.42500
t
2
=
0.61094
t
3
=
1.48750
t
4
=
1.88594
t
5
=
2.94844
Mmioc
113.968
32.137628
(1, 0, 1, 0, 1)
t
1
=
0.48917
t
2
=
1.49831
t
3
=
1.72721
t
4
=
3.11023Figure 6
Sequence of configurations of the optimal solution for the swing-up maneuver for the Acrobot with brake. The corresponding control and state variables are represented in Figure7.Figure 7
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-up maneuver for the Acrobot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 6. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


#### 6.2.2. Swing to the Down Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(38)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.8 [s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 4, Figures 8 and 9 with the same interpretation as the previous sections. It is easy to see from Figure 9(f) that the optimal sequence of modes obtained with Eoc method coincides with that obtained with the Mmioc approach except for the first short phase that does not exist in the Mmioc solution. Nevertheless, the continuous solutions almost coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been5.71916[J]whereas the value obtained with the Mmioc approach has been5.73743[J], slightly higher than the previous one. The computational time to find a solution with the Eoc method has been4.7636[s], whereas it has been19.428886[s]with the Mmioc approach, higher than the previous one.Table 4
Results for the swing-down maneuver for the Acrobot with brake:tf=3.8 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
5.71916
4.7636
(0, 1, 0, 1, 0, 1)
t
1
=
0.02969
t
2
=
0.62344
t
3
=
1.78125
t
4
=
2.79062
t
5
=
3.05781
Mmioc
5.73743
19.428886
(1, 0, 1, 0, 1)
t
1
=
0.72684
t
2
=
1.78373
t
3
=
2.82161
t
4
=
3.06599Figure 8
Sequence of configurations of the optimal solution for the swing-down maneuver for the Acrobot with brake. The corresponding control and state variables are represented in Figure7.Figure 9
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Acrobot with brake. The corresponding sequence of configurations of the robot is represented in Figure 8. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


### 6.3. Brakebot

In this section, the optimal control problem for a Brakebot is studied. In this robot model both joints are not actuated; thusu=(0,0)T. In this case the functional to be minimized is the duration of the motion.The Brakebot has four operational modes depending on which joints are braked and which can rotate freely, namely,(i)
Mode 1: Joint 1 free, Joint 2 free,(ii)
Mode 2: Joint 1 free, Joint 2 braked,(iii)
Mode 3: Joint 1 braked, Joint 2 free,(iv)
Mode 4: Joint 1 braked, Joint 2 braked. From (10), the model of the Brakebot is (39)α+2βcos⁡θ2θ¨1+δ+βcos⁡θ2θ¨2-βsin⁡θ22θ˙1θ˙2+θ˙22+gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2+h1=0,δ+βcos⁡θ2θ¨1+δθ¨2+βsin⁡θ2θ˙12+gm2r2cos⁡θ1+θ2+h2=0. This general dynamic model can be easily particularized to describe the four operating modes of the Brakebot. For example, when the Brakebot is operating in Mode 1, that is, with both Joint 1 and Joint 2 free,h1=h2=0. When the Brakebot is operating in Mode 2, that is, with Joint 1 free and Joint 2 braked,h1=0andh2takes the form (9), and so on. To convert this second order differential model into the first order model (16) the change of variables (24) and the additional differential relations (25) are introduced. The result is (40)α+2βcos⁡x2x˙3+δ+βcos⁡x2x˙4-βsin⁡x22x3x4+x42+gl1m2+m1r1cos⁡x1+gm2r2cos⁡x1+x2+h1=0,δ+βcos⁡x2x˙3+δx˙4+βsin⁡x2x32+gm2r2cos⁡x1+x2+h2=0. Since the Brakebot has four dynamical modes,Dv={1,2,3,4}, four vector fields,fi=(fi,1,fi,2,fi,3,fi,4)T, i=1,…,4, must be considered. The dynamic models of the Brakebot operating in Mode 1, Mode 2, Mode 3, and Mode 4 are the following: (41)x˙1=x3=f1,1,x˙2=x4=f1,2,x˙3=-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2gm2r2cos⁡x1+x2+βx32sin⁡x2α-δδ-β2cos⁡x22=f1,3,x˙4=-gl1m2+m1r1cos⁡x1δ+βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+β2x32+2x3x4+x42cos⁡x2sin⁡x2·δ-α+δ+β2cos⁡x22-1=f1,4,x˙1=x3=f2,1,x˙2=x4=f2,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δh2+βx3+x42sin⁡x2·α-δδ-β2cos⁡x22-1=f2,3,x˙4=-gl1m2+m1r1cos⁡x1δ+βcos⁡x2+h2α+2βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+β2x32+2x3x4+x42cos⁡x2sinx2·δ-α+δ+β2cos⁡x22-1=f2,4,x˙1=x3=f3,1,x˙2=x4=f3,2,x˙3=-δh1-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2gm2r2cosx1+x2+βx32sin⁡x2α-δδ-β2cos⁡x22=f3,3,x˙4=1δ-gm2r2cos⁡x1+x2-βx32sin⁡x2+δ+βcos⁡x2-δh1-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2gm2r2cos⁡x1+x2+βx32sin⁡x2·δ-α+δ+β2cos⁡x22-1=f3,4,x˙1=x3=f4,1,x˙2=x4=f4,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δ-h1+h2+βx3+x42sin⁡x2·α-δδ-β2cos⁡x22-1=f4,3,x˙4=-h2+gm2r2cos⁡x1+x2+βx32sin⁡x2δ+δ+βcos⁡x2-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δ-h1+h2+βx3+x42sin⁡x2·δδ-α+δ+β2cos⁡x22-1=f4,4.Moreover, four new variables,vi∈[0,1], i∈{1,4}, such thatv1+v2+v3+v4=1, must be introduced. Thus, (21) takes the form (42)x˙1=v1f1,1+v2f2,1+v3f3,1+v4f4,1,x˙2=v1f1,2+v2f2,2+v3f3,2+v4f4,2,x˙3=v1f1,3+v2f2,3+v3f3,3+v4f4,3,x˙4=v1f1,4+v2f2,4+v3f3,4+v4f4,4, and the cost functional (22) is (43)JE=tF-tI. This objective functional fulfills the sufficient conditions for the existence of a solution stated in Section 5.Inequality constraints on the state variables have been also considered. In particular, variableθ2has been bounded to take into account possible range limitations of the mechanical system. Thus, the inequality constraints-π[rad]≤θ2≤π[rad]have been introduced in the model.


#### 6.3.1. Swing to the Down Configuration

In this numerical experiment the minimum time problem has been solved with the following initial and final conditions:(44)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial time has beentI=0.The results obtained with both the Eoc and Mmioc methods are reported in Table 5, Figures 10 and 12. Mode 1 has been represented in green, Mode 2 in blue, Mode 3 in orange, and Mode 4 in red. It is easy to see from Figure 11 that the optimal sequence of modes obtained with the Eoc method does not coincide with that obtained with the Mmioc approach. As a consequence, the continuous solutions do not coincide, as well. The value of the objective functional for the optimal solution obtained with the Eoc method has been4.25354[s]whereas the value obtained with the Mmioc approach has been5.00115[J], higher than the previous one. The computational time to find a solution with the Eoc method has been9.29855[s], whereas it has been1270.4665[s]with the Mmioc approach, much higher than the previous one.Table 5
Results for the swing-down maneuver for the Brakebot.
Method
Minimum time [s]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
4.25354
9.29855
(1, 4, 2, 4, 2, 1, 3, 2, 3, 1,  3, 1, 2, 4, 2, 3, 4, 3, 1)
t
1
=
0.03323, t2=0.06646
t
3
=
0.09969, t4=0.23262
t
5
=
0.26585, t6=0.84738
t
7
=
0.89723, t8=1.03016
t
9
=
1.06339, t10=1.54523
t
11
=
1.92739, t12=2.40923
t
13
=
2.49231, t14=2.69170
t
15
=
2.72493, t16=3.88801
t
17
=
4.18708, t18=4.25354
Mmioc
5.00115
1270.4665
(1, 4, 1, 3, 4)
t
1
=
2.26549
t
2
=
3.14850
t
3
=
3.44436
t
4
=
3.80418Figure 10
Sequence of configurations of the optimal solution for the swing-down maneuver for the Brakebot. The configurations in which the Brakebot operates in Mode 1, Mode 2, Mode 3, and Mode 4 are represented in green, blue, orange, and red, respectively. The corresponding control and state variables are represented in Figure7.Figure 11
This figure shows in different colors, namely, Mode 1 in green, Mode 2 in blue, Mode 3 in orange, and Mode 4 in red, the operating modes of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Brakebot. The corresponding sequence of configurations of the robot manipulator is represented in Figure 10. The same figure shows (in black) the operating modes of the optimal solution obtained with a Mmioc approach.
(a)
v
1
(b)
v
2
(c)
v
3
(d)
v
4Figure 12
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Brakebot. The corresponding sequence of configurations of the robot manipulator is represented in Figure 10. The same figure shows (in black) the optimal solution of the same problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2


### 6.4.  Computational Issues

In the Eoc approach the Nlp solver Ipopt [24] has been used (https://projects.coin-or.org/Ipopt), whereas in the Mmioc approach the Minlp solver Bonmin [25] has been used (https://projects.coin-or.org/Bonmin).With both approaches an initial guess of the optimal solution had to be generated. However, with the Mmioc method this generation is much more difficult since, besides the optimal continuous part of the solution, the optimal number of phases, the optimal switching times between them, and the optimal sequence of modes have to be estimated. The initial guess of the optimal switching times has been calculated assuming phases of equal duration. The Mmioc method has shown high sensitivity to the initial guess of the optimal sequence of modes. Additionally, phases of zero length had to be eliminated. Therefore an iterative procedure has been implemented to select the best initial guess and to reduce the number of phases whenever a phase of zero duration appeared in the solution. The numerical experiments have been carried out on a 2.8 GHz Intel Core i7 computer with the Mac OS X Version 10.10.1 operating system and 16 Gb RAM. The Eoc approach has always found solutions with lower or similar performance index than the Mmioc method, and in shorter computational time, especially in the case with four operating modes. Indeed, when the number of modes increases, the Mmioc method becomes slower since it suffers from drawbacks due to the inherent combinatorial nature of the problem, whereas the Eoc approach does not.


## 6.1. Pendubot with Brake

In this section, the optimal control problem for a Pendubot with a brake at the unactuated joint is studied. In this robot model, the second joint is not actuated; thusu=(u1,0)T. In this case, it is not possible to integrate partially or completely the nonholonomic constraint because the conditions for partial integrability [18] are not fulfilled. Hence, the system is controllable. The numerical results of the application of the Eoc method for optimal control to two different boundary value problems for this system will be described. In both cases the functional to be minimized is the energy consumption during the motion.The Pendubot has two operational modes, depending on whether the unactuated joint is braked or not, namely,(i)
Mode 1: Joint 2 free,(ii)
Mode 2: Joint 2 braked. From (7), the dynamic model of a Pendubot with a brake at the unactuated joint is (23)α+2βcos⁡θ2θ¨1+δ+βcos⁡θ2θ¨2-βsin⁡θ22θ˙1θ˙2+θ˙22+gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2=u1,δ+βcos⁡θ2θ¨1+δθ¨2+βsin⁡θ2θ˙12+gm2r2cos⁡θ1+θ2+h2=0. In these equationsu1is the first component of the input vectoru=(u1,u2)T. When the Pendubot is operating in Mode 1, that is, with Joint 2 free,h2=0. Otherwise, it takes the form (9). To convert this second order differential model into the first order model (16) the following change of variables is introduced:(24)x1=θ1,x2=θ2,x3=θ˙1,x4=θ˙2,with the following additional differential relations:(25)x˙1=x3,x˙2=x4,among the components of state variablex=(x1,x2,x3,x4)T. The result is (26)α+2βcos⁡x2x˙3+δ+βcos⁡x2x˙4-βsin⁡x22x3x4+x42+gl1m2+m1r1cos⁡x1+gm2r2cos⁡x1+x2=u1,δ+βcos⁡x2x˙3+δx˙4+βsin⁡x2x32+gm2r2cos⁡x1+x2+h2=0. Since the Pendubot has two dynamical modes,Dv={1,2},2vector fields,f1=(f1,1,f1,2,f1,3,f1,4)Tandf2=(f2,1,f2,2,f2,3,f2,4)T, must be considered. The dynamic models of the Pendubot operating in Mode 1 and Mode 2 are the following: (27)x˙1=x3=f1,1,x˙2=x4=f1,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2gm2r2cos⁡x1+x2+βx32sin⁡x2+δu1,1+βx3+x42sin⁡x2α-δδ-β2cos⁡x22=f1,3,x˙4=u1,1δ+βcos⁡x2-gl1m2+m1r1cos⁡x1δ+βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+βx42cos⁡x2sin⁡x2+βx3x3+x4sin⁡2x2·δ-α+δ+β2cos⁡x22-1=f1,4,x˙1=x3=f2,1,x˙2=x4=f2,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δh2+u2,1+βx3+x42sin⁡x2α-δδ-β2cos⁡x22=f2,3,x˙4=-h2+gm2r2cos⁡x1+x2+βx32sin⁡x2δ+δ+βcos⁡x2-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δh2+u2,1+βx3+x42sin⁡x2·δδ-α+δ+β2cos⁡x22-1=f2,4.In these equations,u1,1is the first component of vectoru1=(u1,1,u1,2)T, the input vector for vector fieldf1, andu2,1is the first component of vectoru2=(u2,1,u2,2)T, the input vector for vector fieldf2. Moreover, two new variables,vi∈[0,1], i∈{1,2}, such thatv1+v2=1, must be introduced. Thus, (21) takes the form (28)x˙1=v1f1,1+v2f2,1,x˙2=v1f1,2+v2f2,2,x˙3=v1f1,3+v2f2,3,x˙4=v1f1,4+v2f2,4, and the cost functional (22) becomes (29)JE=∫t0tfv1u1,12+v2u2,12dt.This objective functional fulfills the sufficient conditions for the existence of a solution stated in Section 5.Inequality constraints on the state variables have been also considered. In particular, variableθ2has been bounded to take into account possible range limitations of the mechanical system. Thus, the inequality constraints-π[rad]≤θ2≤π[rad]have been introduced in the model. Inequality constraints on the available control variable have also been considered. In particular, the constraints-25[Nm]≤u1≤25[Nm]have been introduced in the model to take into account physical limitations of the actuators.The following parameters of the dynamic model have been used:Iz1=Iz2=1[kgm2],l1=l2=1[m], andm1=m2=1[kg].


### 6.1.1. Swing to the Up Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(30)θ1tI=-π2rad,θ1tF=π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.6 [s], respectively.The energy consumption, the resolution time, the number of phases, the sequence of modes, and the corresponding switching times between them obtained with both the Eoc and Mmioc methods are reported in Table 1. Figure 2 shows the sequence of configurations of the robot at timest=k3.6/64, k=0,1,…,64, corresponding to the optimal solution of the boundary value problem obtained with a discretization of the interval[tI,tF]into128subintervals. Configurations that correspond to the braked mode of the mechanical systems are represented in red. It is easy to see that the braked mode is used in the solution just before swinging and when the system is reaching its final state. Figure 3 shows the corresponding control and state variables obtained with the Eoc method (in red) and the Mmioc approach (in black). In the Eoc solution, six phases can be distinguished, although two of them have a very short duration. The optimal sequence of modes obtained has been(0,1,0,1,0,1)which corresponds to braking the unactuated joint in phases2,4, and6and not braking in phases1,3, and5. In the Mmioc solution,4phases can be distinguished. The optimal sequence of modes obtained has been(0,1,0,1)which corresponds to braking the unactuated joint in phases2and4and not braking in phases1and3. It is easy to see from Figure 3(f) that the optimal sequence of modes obtained with the Eoc method coincides with that obtained with the Mmioc approach except for the fifth short phase that does not exist in the Mmioc solution. As a consequence, the continuous solutions almost coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been192.71 [J], whereas the value obtained with the Mmioc approach has been195.532 [J], higher than the previous one. The computational time to find a solution with the Eoc method has been2.84221 [s] whereas it has been26.984307 [s] with the Mmioc approach, higher than the previous one.Table 1
Results for the swing-up maneuver for the Pendubot with brake:tf=3.6 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
192.71
2.84221
(0, 1, 0, 1, 0, 1)
t
1
=
1.01250
t
2
=
2.08125
t
3
=
3.29062
t
4
=
3.54375
t
5
=
3.57187
Mmioc
195.532
26.984307
(0, 1, 0, 1)
t
1
=
0.98867
t
2
=
2.07424
t
3
=
3.29758Figure 2
Sequence of configurations of the optimal solution for the swing-up maneuver for the Pendubot with brake. The configurations in which the Pendubot operates in Mode 1 and Mode 2 are represented in black and red, respectively. The corresponding control and state variables are represented in Figure3.Figure 3
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-up maneuver of the Pendubot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 2. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach. The dashed vertical lines represent switches between phases.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


### 6.1.2. Swing to the Down Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(31)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.6 [s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 2, Figures 4 and 5 with the same interpretation as in Section 6.1.1. It is easy to see from Figure 5(f) that the optimal sequence of modes obtained with the Eoc method coincides with that obtained with the Mmioc approach but the time location of the corresponding phases does not. As a consequence, the continuous solutions do not coincide. Nevertheless, the value of the objective functional for the optimal solution obtained with the Eoc method has been21.5397[J]whereas the value obtained with the Mmioc approach has been21.5739[J], very similar to the previous one. The computational time to find a solution with the Eoc method has been1.57952[s], shorter than computational time to find a solution with the Mmioc approach, which has been28.65894[s].Table 2
Results for the swing-down maneuver for the Pendubot with brake:tf=3.6 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
21.5397
1.57952
(0, 1, 0, 1, 0, 1)
t
1
=
1.01250
t
2
=
1.37813
t
3
=
2.25000
t
4
=
2.67188
t
5
=
3.26250
Mmioc
21.5739
28.65894
(0, 1, 0, 1, 0, 1)
t
1
=
1.11980
t
2
=
1.53434
t
3
=
2.41782
t
4
=
2.81834
t
5
=
3.31516Figure 4
Sequence of configurations of the optimal solution for the swing-down maneuver for the Pendubot with brake. The corresponding control and state variables are represented in Figure5.Figure 5
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Pendubot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 4. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


## 6.1.1. Swing to the Up Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(30)θ1tI=-π2rad,θ1tF=π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.6 [s], respectively.The energy consumption, the resolution time, the number of phases, the sequence of modes, and the corresponding switching times between them obtained with both the Eoc and Mmioc methods are reported in Table 1. Figure 2 shows the sequence of configurations of the robot at timest=k3.6/64, k=0,1,…,64, corresponding to the optimal solution of the boundary value problem obtained with a discretization of the interval[tI,tF]into128subintervals. Configurations that correspond to the braked mode of the mechanical systems are represented in red. It is easy to see that the braked mode is used in the solution just before swinging and when the system is reaching its final state. Figure 3 shows the corresponding control and state variables obtained with the Eoc method (in red) and the Mmioc approach (in black). In the Eoc solution, six phases can be distinguished, although two of them have a very short duration. The optimal sequence of modes obtained has been(0,1,0,1,0,1)which corresponds to braking the unactuated joint in phases2,4, and6and not braking in phases1,3, and5. In the Mmioc solution,4phases can be distinguished. The optimal sequence of modes obtained has been(0,1,0,1)which corresponds to braking the unactuated joint in phases2and4and not braking in phases1and3. It is easy to see from Figure 3(f) that the optimal sequence of modes obtained with the Eoc method coincides with that obtained with the Mmioc approach except for the fifth short phase that does not exist in the Mmioc solution. As a consequence, the continuous solutions almost coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been192.71 [J], whereas the value obtained with the Mmioc approach has been195.532 [J], higher than the previous one. The computational time to find a solution with the Eoc method has been2.84221 [s] whereas it has been26.984307 [s] with the Mmioc approach, higher than the previous one.Table 1
Results for the swing-up maneuver for the Pendubot with brake:tf=3.6 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
192.71
2.84221
(0, 1, 0, 1, 0, 1)
t
1
=
1.01250
t
2
=
2.08125
t
3
=
3.29062
t
4
=
3.54375
t
5
=
3.57187
Mmioc
195.532
26.984307
(0, 1, 0, 1)
t
1
=
0.98867
t
2
=
2.07424
t
3
=
3.29758Figure 2
Sequence of configurations of the optimal solution for the swing-up maneuver for the Pendubot with brake. The configurations in which the Pendubot operates in Mode 1 and Mode 2 are represented in black and red, respectively. The corresponding control and state variables are represented in Figure3.Figure 3
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-up maneuver of the Pendubot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 2. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach. The dashed vertical lines represent switches between phases.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


## 6.1.2. Swing to the Down Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(31)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.6 [s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 2, Figures 4 and 5 with the same interpretation as in Section 6.1.1. It is easy to see from Figure 5(f) that the optimal sequence of modes obtained with the Eoc method coincides with that obtained with the Mmioc approach but the time location of the corresponding phases does not. As a consequence, the continuous solutions do not coincide. Nevertheless, the value of the objective functional for the optimal solution obtained with the Eoc method has been21.5397[J]whereas the value obtained with the Mmioc approach has been21.5739[J], very similar to the previous one. The computational time to find a solution with the Eoc method has been1.57952[s], shorter than computational time to find a solution with the Mmioc approach, which has been28.65894[s].Table 2
Results for the swing-down maneuver for the Pendubot with brake:tf=3.6 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
21.5397
1.57952
(0, 1, 0, 1, 0, 1)
t
1
=
1.01250
t
2
=
1.37813
t
3
=
2.25000
t
4
=
2.67188
t
5
=
3.26250
Mmioc
21.5739
28.65894
(0, 1, 0, 1, 0, 1)
t
1
=
1.11980
t
2
=
1.53434
t
3
=
2.41782
t
4
=
2.81834
t
5
=
3.31516Figure 4
Sequence of configurations of the optimal solution for the swing-down maneuver for the Pendubot with brake. The corresponding control and state variables are represented in Figure5.Figure 5
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Pendubot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 4. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


## 6.2. Acrobot with Brake

In this section, the optimal control problem of an Acrobot with a brake at the unactuated joint is studied. In this robot model the first joint is not actuated; thusu=(0,u2)T. In this case, it is not possible to integrate partially or completely the nonholonomic constraint because the conditions for partial integrability [18] are not fulfilled. Hence, the system is controllable. The numerical results of the application of our method for optimal control to two different boundary value problems for this system will be described. In both cases the functional to be minimized is the energy consumption during the motion.The Acrobot has two operational modes, depending on whether the unactuated joint is braked or not, namely,(i)
Mode 1: Joint 1 free,(ii)
Mode 2: Joint 1 braked. From (8), the dynamic model of an Acrobot with a brake at the unactuated joint is (32)α+2βcos⁡θ2θ¨1+δ+βcos⁡θ2θ¨2-βsin⁡θ22θ˙1θ˙2+θ˙22+gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2+h1=0,δ+βcos⁡θ2θ¨1+δθ¨2+βsin⁡θ2θ˙12+gm2r2cos⁡θ1+θ2=u2. In these equationsu2is the second component of the input vectoru=(u1,u2)T. When the Acrobot is operating in Mode 1, that is, with Joint 1 free,h1=0. Otherwise, it takes the form (9). To convert this second order differential model into the first order model (16) the change of variables (24) and the additional differential relations (25) are introduced. The result is (33)α+2βcos⁡x2x˙3+δ+βcos⁡x2x˙4-βsin⁡x22x3x4+x42+gl1m2+m1r1cos⁡x1+gm2r2cos⁡x1+x2+h1=0,δ+βcos⁡x2x˙3+δx˙4+βsin⁡x2x32+gm2r2cos⁡x1+x2=u2. Since the Acrobot has two dynamical modes,Dv={1,2}, two vector fields,f1andf2, must be considered. The dynamic models of the Acrobot operating in Mode 1 and Mode 2 are the following: (34)x˙1=x3=f1,1,x˙2=x4=f1,2,x˙3=-δu1,2-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2-u1,2+gm2r2cos⁡x1+x2+βx32sin⁡x2α-δδ-β2cos⁡x22=f1,3,x˙4=-gl1m2+m1r1cos⁡x1δ+βcos⁡x2-u1,2α+2βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+β2x32+2x3x4+x42cos⁡x2sinx2·δ-α+δ+β2cos⁡x22-1=f1,4,x˙1=x3=f2,1,x˙2=x4=f2,2,x˙3=-δh1+u2,2-δgl1m2+m1r1cos⁡x1+βδx3+x42sinx2+βcos⁡x2-x6+gm2r2cosx1+x2+βx32sinx2α-δδ-β2cos⁡x22=f2,3,x˙4=1δu2,2-gm2r2cos⁡x1+x2-βx32sin⁡x2+δ+βcos⁡x2-δh1+u2,2-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2-u2,2+gm2r2cos⁡x1+x2+βx32sin⁡x2·δ-α+δ+β2cos⁡x22-1=f2,4.In these equations,u1,2is the second component of vectoru1=(u1,1,u1,2)T, the input vector for vector fieldf1, andu2,2is the second component of vectoru2=(u2,1,u2,2)T, the input vector for vector fieldf2. Moreover, two new variables,vi∈[0,1], i∈{1,2}, such thatv1+v2=1, must be introduced. Thus, (21) takes the form (35)x˙1=v1f1,1+v2f2,1,x˙2=v1f1,2+v2f2,2,x˙3=v1f1,3+v2f2,3,x˙4=v1f1,4+v2f2,4, and the cost functional (22) becomes (36)JE=∫t0tfv1u1,22+v2u2,22dt. This objective functional fulfills the sufficient conditions for the existence of a solution stated in Section 5.Variableθ2has been bounded, to take into account possible range limitations of the mechanical system. Thus, the inequality constraints-π[rad]≤θ2≤π[rad]have been introduced in the model together with the constraints-25[Nm]≤u2≤25[Nm]to take into account physical limitations of the actuators. The same parameters of the dynamic model of Section 6.1 have been used.


### 6.2.1. Swing to the Up Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(37)θ1tI=-π2rad,θ1tF=π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.4[s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 3, Figures 6 and 7 with the same interpretation as in the previous sections. It is easy to see from Figure 7(f) that the optimal sequence of modes obtained with Eoc method does not coincide with that obtained with the Mmioc approach, since the first unbraked phase does not exist in the Mmioc solution. As a consequence, the continuous solutions do not coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been113.968[J]whereas the value obtained with the Mmioc approach has been113.027[J], a bit lower than the previous one. The computational time to find a solution with the Eoc method has been20.2848[s], whereas it has been32.137628[s]with the Mmioc approach, higher than the previous one.Table 3
Results for the swing-up maneuver for the Acrobot with brake:tf=3.4 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
113.027
20.2848
(0, 1, 0, 1, 0, 1)
t
1
=
0.42500
t
2
=
0.61094
t
3
=
1.48750
t
4
=
1.88594
t
5
=
2.94844
Mmioc
113.968
32.137628
(1, 0, 1, 0, 1)
t
1
=
0.48917
t
2
=
1.49831
t
3
=
1.72721
t
4
=
3.11023Figure 6
Sequence of configurations of the optimal solution for the swing-up maneuver for the Acrobot with brake. The corresponding control and state variables are represented in Figure7.Figure 7
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-up maneuver for the Acrobot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 6. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


### 6.2.2. Swing to the Down Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(38)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.8 [s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 4, Figures 8 and 9 with the same interpretation as the previous sections. It is easy to see from Figure 9(f) that the optimal sequence of modes obtained with Eoc method coincides with that obtained with the Mmioc approach except for the first short phase that does not exist in the Mmioc solution. Nevertheless, the continuous solutions almost coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been5.71916[J]whereas the value obtained with the Mmioc approach has been5.73743[J], slightly higher than the previous one. The computational time to find a solution with the Eoc method has been4.7636[s], whereas it has been19.428886[s]with the Mmioc approach, higher than the previous one.Table 4
Results for the swing-down maneuver for the Acrobot with brake:tf=3.8 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
5.71916
4.7636
(0, 1, 0, 1, 0, 1)
t
1
=
0.02969
t
2
=
0.62344
t
3
=
1.78125
t
4
=
2.79062
t
5
=
3.05781
Mmioc
5.73743
19.428886
(1, 0, 1, 0, 1)
t
1
=
0.72684
t
2
=
1.78373
t
3
=
2.82161
t
4
=
3.06599Figure 8
Sequence of configurations of the optimal solution for the swing-down maneuver for the Acrobot with brake. The corresponding control and state variables are represented in Figure7.Figure 9
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Acrobot with brake. The corresponding sequence of configurations of the robot is represented in Figure 8. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


## 6.2.1. Swing to the Up Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(37)θ1tI=-π2rad,θ1tF=π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.4[s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 3, Figures 6 and 7 with the same interpretation as in the previous sections. It is easy to see from Figure 7(f) that the optimal sequence of modes obtained with Eoc method does not coincide with that obtained with the Mmioc approach, since the first unbraked phase does not exist in the Mmioc solution. As a consequence, the continuous solutions do not coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been113.968[J]whereas the value obtained with the Mmioc approach has been113.027[J], a bit lower than the previous one. The computational time to find a solution with the Eoc method has been20.2848[s], whereas it has been32.137628[s]with the Mmioc approach, higher than the previous one.Table 3
Results for the swing-up maneuver for the Acrobot with brake:tf=3.4 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
113.027
20.2848
(0, 1, 0, 1, 0, 1)
t
1
=
0.42500
t
2
=
0.61094
t
3
=
1.48750
t
4
=
1.88594
t
5
=
2.94844
Mmioc
113.968
32.137628
(1, 0, 1, 0, 1)
t
1
=
0.48917
t
2
=
1.49831
t
3
=
1.72721
t
4
=
3.11023Figure 6
Sequence of configurations of the optimal solution for the swing-up maneuver for the Acrobot with brake. The corresponding control and state variables are represented in Figure7.Figure 7
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-up maneuver for the Acrobot with brake. The corresponding sequence of configurations of the robot manipulator is represented in Figure 6. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


## 6.2.2. Swing to the Down Configuration

In this numerical experiment, a boundary value problem has been solved with the following initial and final conditions:(38)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial and final times have beentI=0andtF=3.8 [s], respectively.The results obtained with both the Eoc and Mmioc methods are reported in Table 4, Figures 8 and 9 with the same interpretation as the previous sections. It is easy to see from Figure 9(f) that the optimal sequence of modes obtained with Eoc method coincides with that obtained with the Mmioc approach except for the first short phase that does not exist in the Mmioc solution. Nevertheless, the continuous solutions almost coincide. The value of the objective functional for the optimal solution obtained with the Eoc method has been5.71916[J]whereas the value obtained with the Mmioc approach has been5.73743[J], slightly higher than the previous one. The computational time to find a solution with the Eoc method has been4.7636[s], whereas it has been19.428886[s]with the Mmioc approach, higher than the previous one.Table 4
Results for the swing-down maneuver for the Acrobot with brake:tf=3.8 [s].
Method
Minimum energy [J]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
5.71916
4.7636
(0, 1, 0, 1, 0, 1)
t
1
=
0.02969
t
2
=
0.62344
t
3
=
1.78125
t
4
=
2.79062
t
5
=
3.05781
Mmioc
5.73743
19.428886
(1, 0, 1, 0, 1)
t
1
=
0.72684
t
2
=
1.78373
t
3
=
2.82161
t
4
=
3.06599Figure 8
Sequence of configurations of the optimal solution for the swing-down maneuver for the Acrobot with brake. The corresponding control and state variables are represented in Figure7.Figure 9
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Acrobot with brake. The corresponding sequence of configurations of the robot is represented in Figure 8. The same figure shows (in black) the optimal solution of the same boundary value problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2
(e)
u
1
(f)
v


## 6.3. Brakebot

In this section, the optimal control problem for a Brakebot is studied. In this robot model both joints are not actuated; thusu=(0,0)T. In this case the functional to be minimized is the duration of the motion.The Brakebot has four operational modes depending on which joints are braked and which can rotate freely, namely,(i)
Mode 1: Joint 1 free, Joint 2 free,(ii)
Mode 2: Joint 1 free, Joint 2 braked,(iii)
Mode 3: Joint 1 braked, Joint 2 free,(iv)
Mode 4: Joint 1 braked, Joint 2 braked. From (10), the model of the Brakebot is (39)α+2βcos⁡θ2θ¨1+δ+βcos⁡θ2θ¨2-βsin⁡θ22θ˙1θ˙2+θ˙22+gl1m2+m1r1cos⁡θ1+gm2r2cos⁡θ1+θ2+h1=0,δ+βcos⁡θ2θ¨1+δθ¨2+βsin⁡θ2θ˙12+gm2r2cos⁡θ1+θ2+h2=0. This general dynamic model can be easily particularized to describe the four operating modes of the Brakebot. For example, when the Brakebot is operating in Mode 1, that is, with both Joint 1 and Joint 2 free,h1=h2=0. When the Brakebot is operating in Mode 2, that is, with Joint 1 free and Joint 2 braked,h1=0andh2takes the form (9), and so on. To convert this second order differential model into the first order model (16) the change of variables (24) and the additional differential relations (25) are introduced. The result is (40)α+2βcos⁡x2x˙3+δ+βcos⁡x2x˙4-βsin⁡x22x3x4+x42+gl1m2+m1r1cos⁡x1+gm2r2cos⁡x1+x2+h1=0,δ+βcos⁡x2x˙3+δx˙4+βsin⁡x2x32+gm2r2cos⁡x1+x2+h2=0. Since the Brakebot has four dynamical modes,Dv={1,2,3,4}, four vector fields,fi=(fi,1,fi,2,fi,3,fi,4)T, i=1,…,4, must be considered. The dynamic models of the Brakebot operating in Mode 1, Mode 2, Mode 3, and Mode 4 are the following: (41)x˙1=x3=f1,1,x˙2=x4=f1,2,x˙3=-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2gm2r2cos⁡x1+x2+βx32sin⁡x2α-δδ-β2cos⁡x22=f1,3,x˙4=-gl1m2+m1r1cos⁡x1δ+βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+β2x32+2x3x4+x42cos⁡x2sin⁡x2·δ-α+δ+β2cos⁡x22-1=f1,4,x˙1=x3=f2,1,x˙2=x4=f2,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δh2+βx3+x42sin⁡x2·α-δδ-β2cos⁡x22-1=f2,3,x˙4=-gl1m2+m1r1cos⁡x1δ+βcos⁡x2+h2α+2βcos⁡x2+gm2r2α-δ+βcos⁡x2cos⁡x1+x2+βαx32+δx42x3+x4+β2x32+2x3x4+x42cos⁡x2sinx2·δ-α+δ+β2cos⁡x22-1=f2,4,x˙1=x3=f3,1,x˙2=x4=f3,2,x˙3=-δh1-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2gm2r2cosx1+x2+βx32sin⁡x2α-δδ-β2cos⁡x22=f3,3,x˙4=1δ-gm2r2cos⁡x1+x2-βx32sin⁡x2+δ+βcos⁡x2-δh1-δgl1m2+m1r1cos⁡x1+βδx3+x42sin⁡x2+βcos⁡x2gm2r2cos⁡x1+x2+βx32sin⁡x2·δ-α+δ+β2cos⁡x22-1=f3,4,x˙1=x3=f4,1,x˙2=x4=f4,2,x˙3=-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δ-h1+h2+βx3+x42sin⁡x2·α-δδ-β2cos⁡x22-1=f4,3,x˙4=-h2+gm2r2cos⁡x1+x2+βx32sin⁡x2δ+δ+βcos⁡x2-δgl1m2+m1r1cos⁡x1+βcos⁡x2h2+gm2r2cos⁡x1+x2+βx32sin⁡x2+δ-h1+h2+βx3+x42sin⁡x2·δδ-α+δ+β2cos⁡x22-1=f4,4.Moreover, four new variables,vi∈[0,1], i∈{1,4}, such thatv1+v2+v3+v4=1, must be introduced. Thus, (21) takes the form (42)x˙1=v1f1,1+v2f2,1+v3f3,1+v4f4,1,x˙2=v1f1,2+v2f2,2+v3f3,2+v4f4,2,x˙3=v1f1,3+v2f2,3+v3f3,3+v4f4,3,x˙4=v1f1,4+v2f2,4+v3f3,4+v4f4,4, and the cost functional (22) is (43)JE=tF-tI. This objective functional fulfills the sufficient conditions for the existence of a solution stated in Section 5.Inequality constraints on the state variables have been also considered. In particular, variableθ2has been bounded to take into account possible range limitations of the mechanical system. Thus, the inequality constraints-π[rad]≤θ2≤π[rad]have been introduced in the model.


### 6.3.1. Swing to the Down Configuration

In this numerical experiment the minimum time problem has been solved with the following initial and final conditions:(44)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial time has beentI=0.The results obtained with both the Eoc and Mmioc methods are reported in Table 5, Figures 10 and 12. Mode 1 has been represented in green, Mode 2 in blue, Mode 3 in orange, and Mode 4 in red. It is easy to see from Figure 11 that the optimal sequence of modes obtained with the Eoc method does not coincide with that obtained with the Mmioc approach. As a consequence, the continuous solutions do not coincide, as well. The value of the objective functional for the optimal solution obtained with the Eoc method has been4.25354[s]whereas the value obtained with the Mmioc approach has been5.00115[J], higher than the previous one. The computational time to find a solution with the Eoc method has been9.29855[s], whereas it has been1270.4665[s]with the Mmioc approach, much higher than the previous one.Table 5
Results for the swing-down maneuver for the Brakebot.
Method
Minimum time [s]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
4.25354
9.29855
(1, 4, 2, 4, 2, 1, 3, 2, 3, 1,  3, 1, 2, 4, 2, 3, 4, 3, 1)
t
1
=
0.03323, t2=0.06646
t
3
=
0.09969, t4=0.23262
t
5
=
0.26585, t6=0.84738
t
7
=
0.89723, t8=1.03016
t
9
=
1.06339, t10=1.54523
t
11
=
1.92739, t12=2.40923
t
13
=
2.49231, t14=2.69170
t
15
=
2.72493, t16=3.88801
t
17
=
4.18708, t18=4.25354
Mmioc
5.00115
1270.4665
(1, 4, 1, 3, 4)
t
1
=
2.26549
t
2
=
3.14850
t
3
=
3.44436
t
4
=
3.80418Figure 10
Sequence of configurations of the optimal solution for the swing-down maneuver for the Brakebot. The configurations in which the Brakebot operates in Mode 1, Mode 2, Mode 3, and Mode 4 are represented in green, blue, orange, and red, respectively. The corresponding control and state variables are represented in Figure7.Figure 11
This figure shows in different colors, namely, Mode 1 in green, Mode 2 in blue, Mode 3 in orange, and Mode 4 in red, the operating modes of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Brakebot. The corresponding sequence of configurations of the robot manipulator is represented in Figure 10. The same figure shows (in black) the operating modes of the optimal solution obtained with a Mmioc approach.
(a)
v
1
(b)
v
2
(c)
v
3
(d)
v
4Figure 12
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Brakebot. The corresponding sequence of configurations of the robot manipulator is represented in Figure 10. The same figure shows (in black) the optimal solution of the same problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2


## 6.3.1. Swing to the Down Configuration

In this numerical experiment the minimum time problem has been solved with the following initial and final conditions:(44)θ1tI=π2rad,θ1tF=-π2rad,θ2tI=0rad,θ2tF=0rad,and θ˙1(tI)=θ˙1(tF)=θ˙2(tI)=θ˙2(tF)=0[rad/s]. The initial time has beentI=0.The results obtained with both the Eoc and Mmioc methods are reported in Table 5, Figures 10 and 12. Mode 1 has been represented in green, Mode 2 in blue, Mode 3 in orange, and Mode 4 in red. It is easy to see from Figure 11 that the optimal sequence of modes obtained with the Eoc method does not coincide with that obtained with the Mmioc approach. As a consequence, the continuous solutions do not coincide, as well. The value of the objective functional for the optimal solution obtained with the Eoc method has been4.25354[s]whereas the value obtained with the Mmioc approach has been5.00115[J], higher than the previous one. The computational time to find a solution with the Eoc method has been9.29855[s], whereas it has been1270.4665[s]with the Mmioc approach, much higher than the previous one.Table 5
Results for the swing-down maneuver for the Brakebot.
Method
Minimum time [s]
Resolution time  [s]
Sequence of modes
Switching times [s]
Eoc
4.25354
9.29855
(1, 4, 2, 4, 2, 1, 3, 2, 3, 1,  3, 1, 2, 4, 2, 3, 4, 3, 1)
t
1
=
0.03323, t2=0.06646
t
3
=
0.09969, t4=0.23262
t
5
=
0.26585, t6=0.84738
t
7
=
0.89723, t8=1.03016
t
9
=
1.06339, t10=1.54523
t
11
=
1.92739, t12=2.40923
t
13
=
2.49231, t14=2.69170
t
15
=
2.72493, t16=3.88801
t
17
=
4.18708, t18=4.25354
Mmioc
5.00115
1270.4665
(1, 4, 1, 3, 4)
t
1
=
2.26549
t
2
=
3.14850
t
3
=
3.44436
t
4
=
3.80418Figure 10
Sequence of configurations of the optimal solution for the swing-down maneuver for the Brakebot. The configurations in which the Brakebot operates in Mode 1, Mode 2, Mode 3, and Mode 4 are represented in green, blue, orange, and red, respectively. The corresponding control and state variables are represented in Figure7.Figure 11
This figure shows in different colors, namely, Mode 1 in green, Mode 2 in blue, Mode 3 in orange, and Mode 4 in red, the operating modes of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Brakebot. The corresponding sequence of configurations of the robot manipulator is represented in Figure 10. The same figure shows (in black) the operating modes of the optimal solution obtained with a Mmioc approach.
(a)
v
1
(b)
v
2
(c)
v
3
(d)
v
4Figure 12
This figure shows (in red) the control and state variables of the optimal solution obtained with the Eoc approach for the swing-down maneuver for the Brakebot. The corresponding sequence of configurations of the robot manipulator is represented in Figure 10. The same figure shows (in black) the optimal solution of the same problem obtained with a Mmioc approach.
(a)
θ
1
(b)
θ
2
(c)
θ
˙
1
(d)
θ
˙
2


## 6.4.  Computational Issues

In the Eoc approach the Nlp solver Ipopt [24] has been used (https://projects.coin-or.org/Ipopt), whereas in the Mmioc approach the Minlp solver Bonmin [25] has been used (https://projects.coin-or.org/Bonmin).With both approaches an initial guess of the optimal solution had to be generated. However, with the Mmioc method this generation is much more difficult since, besides the optimal continuous part of the solution, the optimal number of phases, the optimal switching times between them, and the optimal sequence of modes have to be estimated. The initial guess of the optimal switching times has been calculated assuming phases of equal duration. The Mmioc method has shown high sensitivity to the initial guess of the optimal sequence of modes. Additionally, phases of zero length had to be eliminated. Therefore an iterative procedure has been implemented to select the best initial guess and to reduce the number of phases whenever a phase of zero duration appeared in the solution. The numerical experiments have been carried out on a 2.8 GHz Intel Core i7 computer with the Mac OS X Version 10.10.1 operating system and 16 Gb RAM. The Eoc approach has always found solutions with lower or similar performance index than the Mmioc method, and in shorter computational time, especially in the case with four operating modes. Indeed, when the number of modes increases, the Mmioc method becomes slower since it suffers from drawbacks due to the inherent combinatorial nature of the problem, whereas the Eoc approach does not.


## 7. Conclusions

In this paper, the embedding approach has been used to efficiently solve optimal control problems for planar underactuated robot manipulators with two revolute joints and brakes at the unactuated joints, and in the presence of gravity. The control problem of these switched dynamical systems is particularly difficult to solve due to the combinatorial nature of the problem and to the presence of nonholonomic differential constraints. It has been shown that this method does not require using integer variables or computing the optimal switching times. It requires neither guessing the number of modes nor the initial sequence of modes. The results of the numerical experiments have shown the effectiveness of the embedding approach in solving optimal control problems of underactuated mechanical systems.



---

*Source: 348178-2015-09-15.xml*