Construction of control barrier function and C 2 reference trajectory for constrained attitude maneuvers
CConstruction of control barrier function and C reference trajectory forconstrained attitude maneuvers Xiao Tan, and Dimos V. Dimarogonas
Abstract — Constrained attitude maneuvers have numerousapplications in robotics and aerospace. In our previous work,a general framework to this problem was proposed withresolution completeness guarantee. However, a smooth refer-ence trajectory and a low-level safety-critical controller werelacking. In this work, we propose a novel construction of a C continuous reference trajectory based on B´ezier curves on SO (3) that evolves within predetermined cells and eliminatesprevious stop-and-go behavior. Moreover, we propose a novelzeroing control barrier function on SO (3) that provides asafety certificate over a set of overlapping cells on SO (3) whileavoiding nonsmooth analysis. The safety certificate is given asa linear constraint on the control input and implemented inreal-time. A remedy is proposed to handle the states where thecoefficient of the control input in the linear constraint vanishes.Numerical simulations are given to verify the advantages of theproposed method. I. I
NTRODUCTION
The study of the attitude (orientation) control problemarised from early space and aerial applications and becameprevalent in autonomous robotic systems. A recent trend inthis field is to study this problem using Lie group theory[ ? ], [1]–[3], motivated by the fact that there exists noattitude parametrization other than the special orthogonalgroup SO (3) that both globally and uniquely represents therotational space and avoids singularities and the unwindingphenomenon. Many safe-critical applications, such as spacetelescopes observing some celestial regions while avoidingbright stars [4], and the anisotropy sensitive imaging andcommunication payloads on UAVs and AUVs, motivatefurther study of the attitude planning and control problemin the presence of orientation constraints (i.e., unfeasiblerotational regions).There exist two main approaches for the constrained at-titude maneuver problem: the potential-function [5]–[7] andthe planning based methods [8]–[10]. By delicately designinga potential function, the feedback controller utilizes thenegative gradient to guide the rotational movement. Gen-erally speaking, potential-function based methods are easyto implement, but the state trajectory may get stuck at localminima (where the gradient vanishes) and requires convexityof the safe regions. On the contrary, planning-based methodstry to find a feasible trajectory leading to the target state andthen a tracking controller is utilized. This approach, however, *This work was supported by the H2020 ERC Starting Grant BUCOPH-SYS, the SSF COIN project, the Swedish Research Council (VR) and theKnut och Alice Wallenberg Foundation.Xiao Tan and Dimos V. Dimarogonas are with the Division of Decisionand Control Systems, KTH Royal Institute of Technology, Stockholm,Sweden. [email protected], [email protected] mainly suffers from the lack of completeness guarantees, i.e.,derive a solution if it exists, and safety guarantees, i.e., acertificate that the actual trajectory will not derivate fromthe reference and remain in the safe region.In our previous work [11], a hierarchy framework wasproposed for the constrained attitude maneuver problemconsisting of 1) discretizing the rotation group SO (3) intofinite overlapping cells, 2) planning over the cells, and 3)reference trajectory generation and tracking control. Viewingthe sampling step as the resolution level, we guarantee afeasible path is to be found in finite time when one exists atthat resolution. However, the reference trajectory in [11] isconstructed by the concatenation of geodesic paths and has toreach zero velocities at end points for each sub-maneuvers.This is a potential drawback as it requires the vehicle tostop and go from configuration to configuration. Moreover,no safety guarantee is developed for the low-level trackingcontroller.In this work, we follow the framework in [11] and furtherconstruct a C reference trajectory and develop a safetycertificate by designing zeroing control barrier functionson SO (3) . The C reference trajectory is generated by aB´ezier curve on SO (3) . By choosing the controlling pointscarefully, we show that the constructed curve is of C continuity, connects the initial and target orientations, andevolves within the set of given cells.This paper has two additional contributions: 1) Notingthat the safety region is a union of a set of overlappingcells, we formulate a smooth control barrier function andthus avoid the nonsmooth analysis as the case in [12]. Theproposed formulation is at the cost of shrinking the safetyregion and this conservativeness can be explicitly adjusted bya user-defined parameter; 2) Since the Lie derivative of thecontrol barrier function candidate vanishes at certain states,existing high-order control barrier function design methods[13], [14] are not directly applicable. To address this issue,we introduce a remedy with a proof to render the constrainton the control input feasible for all states in the safety region.All results are illustrated though relevant simulations.II. P RELIMINARIES AND PROBLEM FORMULATION
The set of real, non-negative real, and positive integernumbers are denoted as R , R ≥ , N , respectively. R n denotesthe n -dimensional Euclidean space. The 2-norm of a vector x ∈ R n is (cid:107) x (cid:107) := √ x (cid:62) x . I is the -dimensional identitymatrix. The Frobenius norm of A is defined as (cid:107) A (cid:107) F = tr ( A (cid:62) A ) , where tr ( · ) denotes the trace of a matrix. The Liederivatives of a function h ( x ) for the system ˙ x = f ( x ) + a r X i v : . [ ee ss . S Y ] A ug ( x ) u are denoted by L f h := ∂h∂x f ( x ) and L g h := ∂h∂x g ( x ) ,respectively. A continuous function α : ( − b, a ) → ( −∞ , ∞ ) is said to belong to extended class K for some a, b > if itis strictly increasing and α (0) = 0 [15].Any rotation matrix is an element of the Special Or-thogonal group SO (3) := { R ∈ R × | R (cid:62) R = RR (cid:62) = I, det ( R ) = 1 } which, when associated with the matrixmultiplication operation, forms a Lie group. The associatedLie algebra, denoted by so (3) , consists of the set of allskew-symmetric × matrices, i.e., so (3) := { Ω ∈ R × :Ω (cid:62) = − Ω } . The Lie bracket for so (3) is given as [ V, W ] =
V W − W V , for any
V, W ∈ so (3) . The map [( · )] × : R → so (3) and its inverse map ∨ : so (3) → R are explicitlydefined as x = (cid:16) x x x (cid:17) [( · )] × −−− (cid:42)(cid:41) −−− ( · ) ∨ [ x ] × = (cid:16) − x x x − x − x x (cid:17) . TheLie algebra so (3) allows to represent rotation matrices on SO (3) via the matrix exponential exp( · ) . For [ x ] × ∈ so (3) , exp([ x ] × ) = I + sin( (cid:107) x (cid:107) ) (cid:107) x (cid:107) [ x ] × + − cos( (cid:107) x (cid:107) ) (cid:107) x (cid:107) [ x ] × when x (cid:54) = 0 ,and exp([ x ] × ) = I otherwise [16]. For all rotation matrices R with tr ( R ) (cid:54) = − , the exponential map admits an inverselogarithmic map given by log( R ) = θ ( R )2 sin( θ ( R )) ( R − R (cid:62) ) when R (cid:54) = I , and log( R ) = 0 otherwise, where θ ( R ) :=arccos (( tr ( R ) − / is the rotation angle associated to R [16].A path F ( · ) in A connecting R ∈ SO (3) and R ∈ SO (3) , where A is a subset of SO (3) , is defined as acontinuous function F : [ a, b ] → A with F ( a ) = R and F ( b ) = R . If there exists such a path F ( · ) , we say ( R , R ) is connected . If any two points in A are connected, then wecall the set A path-connected . Given any R , R ∈ SO (3) with tr ( R (cid:62) R ) (cid:54) = − , the geodesic path between R and R is R ( τ ) = R exp( τ log( R (cid:62) R )) , ≤ τ ≤ . The angular distance between R , R is given by d ( R , R ) := (cid:107) log ( R R (cid:62) ) (cid:107) .In [11], we proposed a general framework for the con-strained attitude maneuver problem consisting of SO (3) space partitioning, planning, and reference trajectory gener-ation. We briefly recall the results here. Let the samplingset U := { R , R , . . . , R i , . . . , R n } be a finite set with n elements in SO (3) and let N (cid:48) := { , , . . . , n } be anindex set. For each i ∈ N (cid:48) , define the cell region S i asthe open ball centered at R i with a radius θ ∈ (0 , π/ ,i.e., S i := { R ∈ SO (3) : d ( R, R i ) < θ } , ∀ i ∈ N (cid:48) . Theneighborhood set N i of R i is defined as N i := { R ∈ U : d ( R, R i ) < θ, R (cid:54) = R i } , ∀ i ∈ N (cid:48) . Cells S i , S j are adjacent if S i ∩ S j (cid:54) = ∅ .By choosing U and θ such that the conditions in [11,Theorem 1] are satisfied, we guarantee that for an arbitrarycell S i ∈ N (cid:48) , it has adjacent cells; the center points ofadjacent cells are strictly outside of S i ; the union of thecells covers the whole SO (3) space. Mathematically,i. For all i ∈ N (cid:48) , N i (cid:54) = ∅ ;ii. For all i, j ∈ N (cid:48) , i (cid:54) = j , we have R j / ∈ S i ;iii. For all R i ∈ U , and all R j ∈ N i , θ < d ( R i , R j ) < θ ;iv. ∪ i ∈N (cid:48) S i = SO (3) . Lemma 1 ([11]) . For any cell S i , i ∈ N (cid:48) and two arbitrarypoints R i , R i ∈ S i , the geodesic path between R i and R i is within S i , i.e., for any R i ∈ S i , R i ∈ S i , R ( τ ) = R i exp( τ log( R (cid:62) i R i )) ∈ S i , ≤ τ ≤ , holds. Lemma 2 ([11]) . The geodesic path between any twoneighboring sampling rotations R i and R j is within S i ∪ S j ,i.e., R ( τ ) = R i exp( τ log( R (cid:62) i R j )) ∈ S i ∪ S j , ≤ τ ≤ . We approximate a generic safe attitude zone on SO (3) by a set of cells { S i } , i ∈ N , N ⊂ N (cid:48) and a graphsearch algorithm is utilized that gives out a sequence of cellswhenever feasible at the given resolution level. Without lossof generality, by re-labeling the cells, we assume that theinitial orientation R ∈ S , the target orientation R f ∈ S m , S i and S i +1 are adjacent cells for i ∈ { , · · · , m − } .Based on Lemmas 1 ,
2, a center-to-center attitude maneuverwas then designed, as illustrated in blue dash line in Fig. 1.Though the reference trajectory is guaranteed to be withinthe feasible region, it is not favorable in practice as it requiresthe rigid-body to reach zero velocities at the end points foreach sub-maneuvers.
A. Problem formulation
The attitude dynamics of a rigid body are given by (cid:26) ˙ R = R [ ω ] × ,J ˙ ω + [ ω ] × Jω = u, (1)where the attitude R ∈ SO (3) , ω ∈ R is the angularvelocity in the body-fixed frame, J is the constant and knowninertia matrix and u ∈ R is the input torque. Given a set ofcells { S i } , i ∈ N , N ⊂ N (cid:48) , we call a trajectory γ : t (cid:55)→ R ( t ) is safe if the trajectory always evolves within ∪ i ∈N S i .The control scheme consists of two parts: reference gen-eration and trajectory tracking. Problem 1. ( reference generation) Given a set of cells { S i } , i ∈ N , N ⊂ N (cid:48) such that ∪ i ∈N S i is path-connected.For any given R , R f ∈ ∪ i ∈N S i , find a C curve γ : R ≥ → ∪ i ∈N S i such that dγ/dt (0) = dγ/dt ( T ) =0 , D γ/dt (0) = D γ/dt ( T ) = 0 , γ (0) = R , γ ( t ) = R f , ∀ t ≥ T . Problem 2. ( trajectory tracking) Given a C curve γ : R ≥ → ∪ i ∈N S i , design a control law u for the system (1)such that R ( t ) ∈ ∪ i ∈N S i for t ≥ and lim t →∞ R ( t ) = γ ( t ) .In the following, we will solve Problem 1 and Problem 2in Section III and Section IV, respectively.III. B ´ EZIER CURVE CONSTRUCTION OVER CELLS
In this section, we construct a reference trajectory based onB´ezier curve on SO (3) that solves Problem 1. B´ezier curveis chosen here because De Casteljau Algorithm, which gen-erates B´ezier curves, is in essence a geometric construction, For a curve γ : R → SO (3) , D γ/dt represents the geometricacceleration instead of the second-order total derivatives, following theterminology in [16]. SO (3) manifold, while othersplines are not defined / easy to compute on SO (3) . A. De Casteljau Algorithm on SO (3) We briefly recall De Casteljau Algorithm from [17] asfollows. Taking the geodesics on SO (3) as the analog ofstraight lines, De Casteljau Algorithm connects two pointsin SO (3) via an iterative linear interpolation process. Let n + 1 ordered points of SO (3) be { x , x , · · · , x n } . Thesequence of curves is defined recursively on SO (3) as x ki ( τ ) = x k − i − ( τ ) exp( τ log([ x k − i − ( τ )] (cid:62) x k − i ( τ ))) ,k = 0 , , · · · , n, i = k, k + 1 , · · · , n, (2)where x i ( τ ) = x i . The B´ezier curve is then given by x nn ( τ ) = x n − n − ( τ ) exp( τ log([ x n − n − ( τ )] (cid:62) x n − n ( τ ))) . (3) Lemma 3 ( [18]) . Let n + 1 ordered points of SO (3) be { x , x , · · · , x n } . The corresponding B´ezier curve generatedfrom (3) satisfies the following boundary conditions: x nn (0) = x , x nn (1) = x n ,ddτ x nn ( τ ) | τ =0 = nx V , ddτ x nn ( τ ) | τ =1 = nx n V n − ,D dτ x nn ( τ ) | τ =0 = n ( n − x Υ − ( V − V ) ,D dτ x nn ( τ ) | τ =1 = n ( n − x n Υ − ( V n − − V n − ) , (4) where V i = log( x (cid:62) i x i +1 ) ∈ so (3) , i = 0 , , · · · , n − , Υ − and Υ − are respectively the inverses of theoperators Υ ( W ) = (cid:82) exp( u ad V ) W du, Υ ( W ) = (cid:82) exp( − u ad V n − ) W du.
For any W ∈ so (3) , the operator Υ : so (3) → so (3) isgiven explicitly by the power series (cid:82) exp( u ad V ) W du = (cid:82) W + u [ V , W ] + u [ V , [ V , W ]] + · · · du . The operator Υ ( W ) is given in a similar way. It can be easily verified thatboth operators Υ , Υ are linear transformations on so (3) ,i.e., Υ i ( aW ) = a Υ i ( W ) , Υ i ( W + V ) = Υ i ( W ) + Υ i ( V ) for W, V ∈ so (3) , a ∈ R , i = 1 , . In [18], it is shown that theinverse operator Υ − i always exists for i = 1 , .Lemma 3 introduces the analytical expression of thevelocity and geometric acceleration at the boundary pointthat will facilitate our construction of the reference trajectorywith C continuity. Remark 1.
Note that the De Casteljau algorithm in (2) is notwell-defined for arbitrary points x , x , · · · , x n on SO (3) when tr ([ x k − i − ( τ )] (cid:62) x k − i ( τ )) = − occurs. B. B´ezier curve construction in one cell
In this subsection, we demonstrate the procedure to designthe controlling points and the properties of the constructedB´ezier curve.Given a cell S i , i ∈ N with center point x and twoarbitrary points x , x ∈ S i , the curve c x ,x ,x : [0 , → SO (3) is generated as follows: first add controlling points x , x as the midpoints of x , x , and x , x , respectively; then applying De Casteljau algorithm with n = 4 . Theconstruction is given in Algorithm 1: Algorithm 1
B´ezier curve construction in one cell.
Input: start point x , cell center x , end point x Output: curve c x ,x ,x x ← x exp(1 / x (cid:62) x )) x ← x exp(1 / x (cid:62) x )) calculate a sequence of curves recursively as in (2) giventhe ordered points { x , x , x , x , x } with n = 4 return c x ,x ,x ← x Noticing that V = log( x (cid:62) x ) = 1 / x (cid:62) x ) , V =log( x (cid:62) x ) = 1 / x (cid:62) x ) , we have V = V . Similarly, V = V . From Lemma 3, we can easily check that c x ,x ,x (0) = x , c x ,x ,x (1) = x . The velocities at theboundary point are ddτ c x ,x ,x (0) = 2 x log( x (cid:62) x ) , ddτ c x ,x ,x (1) = 2 x log( x (cid:62) x ) , (5)and the geometric accelerations are given by D dτ c x ,x ,x (0) = 12 x Υ − ( V − V ) = 0 , D dτ c x ,x ,x (1) = 12 x n Υ − ( V − V ) = 0 , (6)noticing that V − V = V − V = 0 and Υ − , Υ − beinglinear transformations.In addition to these explicitly expressed velocities andgeometric accelerations at the endpoints, we have anothernice property of the constructed curve c x ,x ,x . Proposition 1.
Given arbitrary n + 1 ordered points { x , x , x , · · · , x n } such that x i ∈ S , i = 0 , , · · · , n ,where S is a ball region in SO (3) with radius θ ∈ (0 , π/ .The B´ezier curve x nn ( τ ) generated from (3) always exists andevolves in S , i.e., x nn ( τ ) ∈ S, ≤ τ ≤ . Proof.
This can be shown by induction on k . For k = 0 ,as x i ( τ ) = x i , we have x i ( τ ) ∈ S for i = 0 , , · · · , n andd ( x i − ( τ ) , x i ( τ )) < π , which means log([ x i − ( τ )] (cid:62) x i ( τ )) is well-defined and from Lemma 1, x i ( τ ) ∈ S . For any k ∈{ , , · · · , n } , assume that x k − i ( τ ) ∈ S for i = k − , · · · , n ,then x ki ( τ ) = x k − i − ( τ ) exp( τ log([ x k − i − ( τ )] (cid:62) x k − i ( τ ))) iswell-defined as d ( x k − i − ( τ ) , x k − i ( τ )) < π . As x ki ( τ ) lies inthe geodesic path between x k − i − ( τ ) and x k − i ( τ ) , Lemma1 dictates that x ki ( τ ) ∈ S . Thus, by induction, we obtain x nn ( τ ) ∈ S, ≤ τ ≤ . A straightforward conclusion is that the curve c x ,x ,x constructed from Algorithm 1 is well-defined and evolveswithin the cell. C. B´ezier curve construction in a set of cells
Now we apply Algorithm 1 to generate a curve evolvingamong a set of cells. In the following, we use the notation c x ,x ,x ( τ ) : [0 , → SO (3) to denote the curve generatedfrom Algorithm 1 given the three points x , x , x .3 roposition 2. Assume that R , R f are the initial and targetorientations, respectively, R ∈ S , R f ∈ S m , and assumethere exists a sequence of cells S S · · · S m such that S i S i +1 are adjacent for i = 1 , , · · · , m − . Then, a curve c : R ⊃ [0 , m ] → SO (3) defined as c ( τ ) = c R ,R ,R , ( τ ) , τ ∈ [0 , ,c R i − ,i ,R i ,R i,i +1 ( τ − i + 1) , τ ∈ [ i − , i ) ,i ∈ { , , · · · , m − } ,c R m − ,m ,R m ,R f ( τ − m + 1) , τ ∈ [ m − , m ] , (7) where R i is the center of cell S i , R i,i +1 := R i exp(1 / R (cid:62) i R i +1 )) , has the following properties: i. c (0) = R , c ( m ) = R f ; ii. c ( τ ) is a C curve; iii. c ( τ ) ∈ ∪ mi =1 S i for τ ∈ [0 , m ] .Proof. Property i can be straightforwardly verifiedby substituting τ = 1 , , · · · , m − and the factthat c x ,x ,x (0) = x , c x ,x ,x (1) = x . As c ( τ ) is a continuous and piecewise smooth curve, weneed to check the left/right velocity/acceleration at τ = 1 , , · · · , m − . By differentiating (7) andusing (5), we obtain, for i ∈ { , , , · · · , m − } , dc ( τ ) dτ | τ = i − = 2 R i,i +1 log( R (cid:62) i R i,i +1 ) , dc ( τ ) dτ | τ = i + =2 R i,i +1 log( R (cid:62) i,i +1 R i +1 ) . Note that since R i,i +1 = R i exp(1 / R (cid:62) i R i +1 )) , we get log( R (cid:62) i R i,i +1 ) = 1 / R (cid:62) i R i +1 ) , log( R (cid:62) i,i +1 R i +1 ) =log(exp( − / R (cid:62) i R i +1 )) R (cid:62) i R i +1 ) =1 / R (cid:62) i R i +1 ) . As the left and right derivative coincideat τ = 1 , , · · · , m − , c ( τ ) is at least a C curve. From(6), the geometric acceleration at τ = 1 , , · · · , m − satisfies D c ( τ ) dτ | τ = i − = D c ( τ ) dτ | τ = i + = 0 . Thus, c ( τ ) is a C curve.Property iii can be verified piecewise. Rewrite the curvesegments in (7) in the form of c R i ,R i ,R i ( τ ) , τ ∈ [0 , , i ∈{ , , · · · , m } . Proposition 1 implies c R i ,R i ,R i ( τ ) ∈ S i .Thus, the concatenation of the curve segments is containedin the union of the cells, which completes the proof. D. Time re-parameterization
In order to obtain a reference trajectory that solves Prob-lem 1, let τ be a smooth function of time, i.e., τ : R ≥ → [0 , m ] that rescales the trajectory c : [0 , m ] → SO (3) to thetime domain γ := c ◦ τ : R ≥ → SO (3) .Numerous smooth transition functions are known. Herewe adopt one from [19] that fits our needs. s ( x ) = x ∈ ( −∞ , , ρ ( x ) ρ ( x )+ ρ (1 − x ) x ∈ [0 , , x ∈ [1 , ∞ ) (8)with ρ ( x ) := (1 /x ) e − /x . Theorem 1.
Given a sequence of cells S , S , · · · , S m suchthat R ∈ S , R f ∈ S m , S i S i +1 are adjacent. Choose c : [0 , m ] → SO (3) defined in (7) and τ ( t ) := ms ( t/T ) with s ( · ) in (8) . The curve γ := c ◦ τ : R ≥ → SO (3) is continuously differentiable, and satisfies γ (0) = R , γ ( T ) = R f , dγ/dt (0) = dγ/dt ( T ) = 0 , D γ/dt (0) = D γ/dt ( T ) = 0 , γ ( t ) ∈ ∪ i ∈N S i for t ≥ .Proof. Since c is C continuous and s ( · ) is smooth, γ isalso C continuous. From Proposition 2 and the proper-ties of s ( · ) , it can be checked that γ (0) = R , γ ( T ) = R f , D γ/dt (0) = D γ/dt ( T ) = 0 , γ ( t ) ∈ ∪ i ∈N S i for t ≥ . Noticing the fact that s ( x ) is a smooth function with d i s/dx i (0) = d i s/dx i (1) = 0 for any integer i (proven in[19]), we conclude that dγ/dt (0) = dγ/dt ( T ) = 0 . Remark 2.
Although in this work we set the initial and ter-minal velocities to be zero, the presented method can be di-rectly extended to solve interpolation problems with non-zerovelocity boundary conditions by manipulating the controllingpoints in the cells S , S m and the time-reparametrizationfunction s ( · ) . This is a straightforward extension and detailsare omitted here.We demonstrate in Fig. 1 the constructed reference curve c ◦ τ (red line) and the curve from [11] (blue dash line) forcomparison. The data is given in the simulation section. It isseen that a smoother attitude maneuver is obtained comparedto that of [11]. Figure 2 further shows that the trajectory in[11] needs to reach zero velocities at intermediate points,which is avoided in the new construction. The maximalangular velocity magnitude has also decreased compared tothat of [11].IV. C ONTROL BARRIER FUNCTION DESIGN
In this section, we present the procedure to construct azeroing control barrier function that guarantees the actualattitude trajectory evolves within ∪ i ∈N S i .We start the barrier function design from one cell. For anarbitrary cell S i , define a function r i : SO (3) → R r i ( R ) = (cid:15) − (cid:107) R i − R (cid:107) F / , (9)where constant (cid:15) := 4 sin ( θ/ , R i , θ are the cell centerand radius of cell S i , respectively. It is easy to show that r i ( R ) > if and only if R ∈ S i , in view of the fact that (cid:107) v − w (cid:107) F = 2 √ d ( v, w ) / holds for v, w ∈ SO (3) .If we need to constrain the trajectory in cell S i , r i ( R ) is anatural candidate as a zeroing control barrier function as itindicates how far the state is from the cell boundary. Notethat there are many alternatives r i ( R ) to (9), for example, r i ( R ) = θ − d ( R, R i ) . The reason we choose r i ( R ) as in(9) is merely to simplify the expression of its derivatives, asshown later.To ensure the actual attitude trajectory evolves within ∪ i ∈N S i , we need for every time instant t ≥ , there existsat least one i ∈ N that R ( t ) ∈ S i , i.e., max i ∈N ( r i ( R ( t ))) > , for t ≥ . (10)This max operation would lead to nonsmooth analysis anda complex formulation [12]. In the following, we will showhow to circumvent the nonsmooth analysis.4 a) Trajectory of x -axis. (b) Trajectory of y -axis. (c) Trajectory of z -axis. Fig. 1: Comparison of the trajectories of body-fixed axes: c ◦ τ in red and the one from [11] in blue. A ngu l a r v e l o c i t y m agn i t ude (r ad / s ) Fig. 2: Time histories of the reference angular velocitymagnitude of the trajectory c ◦ τ in red and the one from[11] in blue.Define h ( R ) = (cid:88) i ∈N s ( r i ( R ) /(cid:15) ) − δ, (11)where δ > is a user-defined constant, and s ( · ) is givenin (8). The associated constrained set is thus C Rh = { R ∈ SO (3) : h ( R ) ≥ } . Since ∪ i ∈N S i = { R : h ( R ) > − δ } , it is straightforward that C Rh ⊂ ∪ i ∈N S i , and theconstant δ represents the safety margin. The conservativenessis illustrated in Fig. 3 in the planar case. For any given C curve c ◦ τ : R ≥ → ∪ i ∈N S i , we can find a safety margin(i.e., δ ) such that the curve c ◦ τ evolves within C Rh . In thefollowing, we thus assume c ◦ τ ( t ) ∈ C Rh for t ≥ .We embed the attitude dynamics in (1) in a higher dimen-sional Euclidean space as ˙ x := f ( x ) + gu, (12)where x = ( r , r , · · · , r , ω , ω , ω ) ∈ R , f ( x ) = (cid:0) r ω − r ω ; r ω − r ω ; r ω − r ω ; r ω − r ω ; r ω − r ω ; r ω − r ω ; r ω − r ω ; r ω − r ω ; r ω − r ω ; J − ( − [ ω ] × Jω ) (cid:1) ∈ R , g = (cid:18) × J − (cid:19) . . . . . . . . . . . . . . . Fig. 3: Illustration of C h in the planar case with differentconservativeness δ (cid:48) s.This is equivalent to (1) by rewriting the attitude dynamics ina vectorized manner. Note that for r , r , · · · , r , there ex-ist implicit equality constraints since they are elements of arotation matrix. We denote the corresponding − dimensionalsubmanifold C T SO (3) := { x ∈ R : (cid:16) x x x x x x x x x (cid:17) ∈ SO (3) } .Moreover, if x (0) ∈ C T SO (3) , then given any control signal u which is Lipschitz continuous in x , the solution of thedynamical system (12) satisfies x ( t ) ∈ C T SO (3) for t ≥ .This fact can be easily obtained considering that (12) and (1)are equivalent. h in (11) is thus a function of the system states x , in particular, of the states ( x , x , · · · , x ) . The associateconstrained set is C h := { x ∈ C T SO (3) : h ( x ) ≥ } .For all x ∈ C h , we obtain L g h = 0 , and L g L f h ( x ) mayvanish at some points in C h (see Appendix for derivations).Here we note that the higher-order control barrier functiondesign developed in [13], [14] is not directly applicableas a result of this issue. More specifically, to render theset C h forward invariant, the existing methods enforce itby requiring ˙ h ( x ) + α ( h ( x )) ≥ for all x ∈ C h . Re-denote h ( x ) := ˙ h ( x ) + α ( h ( x )) = L f h + α ( h ( x )) . Inorder to show h ( x ) ≥ , it is again enforced by a linearconstraint L g L f h ( x ) u + L f h ( x ) + β ( h ( x )) ≥ on u with α, β extended class K functions. However, this inequalityconstraint may not be feasible when L g L f h ( x ) vanishes.A key observation regarding the singularity set D = { x ∈ C h : L g L f h ( x ) = 0 } is given below:5 roposition 3. Let D = { x ∈ C h : L g L f h ( x ) = 0 } . Then,there exists a constant ξ > such that inf x ∈D h ( x ) ≥ ξ .Proof. See Appendix.In the following a remedy is derived to handle the singu-larity set D = { x ∈ C h : L g L f h ( x ) = 0 } effectively. Denotethe associated set C h,ξ = { x ∈ C T SO (3) : h ( x ) ≥ ξ } .Let χ ( · ) be a twice differentiable function satisfying thefollowing properties χ (0) = 0 ,χ ( a ) = 1 , for a ≥ , dχdι ( a ) > , for a < . (13)Then we smoothly truncate h ( x ) by the upper bound ξ , i.e., b ( x ) = χ ( h ( x ) /ξ ) (14)It is easy to verify that C h = C b := { x ∈ C T SO (3) : b ( x ) ≥ } . Thus, in the following we show the forward invarianceof the set C b instead.We adopt the procedure of the higher-order control barrierfunction design as in [13], [14]. The idea is briefly presentedhere: from Brezis version of Nagumo’s theorem (see [20,Theorem 4.7] and [21, Theorem 4] for a detailed account),the forward invariance of the set C b is guaranteed by showingthat on the boundary of C b , the system states are directed intothe interior or along the boundary of the constrained set. Thiscondition is enforced as ˙ b ( x ) ≥ − α ( b ( x )) for all x ∈ C b ,where α is a continuously differentiable, extended class K function. Note that L g b = 0 , then b ( x ) := ˙ b ( x )+ α ( b ( x )) = L f b + α ( b ( x )) is still a function of the state x . To guaranteethe forward invariance of the set C b := { x ∈ C T SO (3) : b ( x ) ≥ } , using Brezis version of Nagumo’s theoremagain, we have the new condition ˙ b ( x ) ≥ − β ( b ( x )) forall x ∈ C b ∩ C b , where β is a continuously differentiable,extended class K function. Thus, the condition we willenforce in real-time is given as L g b ( x ) u + L f b ( x ) + β ( b ( x )) ≥ (15)for all x ∈ C b ∩ C b . The feasibility result is as follows: Proposition 4.
The inequality condition on u ∈ R L g b ( x ) u + L f b ( x ) + β ( b ( x )) ≥ (16) is feasible for all x ∈ C b ∩ C b .Proof. We examine the condition in two cases. If x ∈ C h,ξ ∩ C b , we have L g b = 0 , L f b = 0 , thus inequality in (16)is equivalent to β ( b ( x )) ≥ , which is trivially satisfied. If x ∈ ( C h \ C h,ξ ) ∩ C b , L g b (cid:54) = 0 , thus we can always finda u that satisfies (16).Suppose a nominal bounded control input u nom ( x ) , Lip-schitz continuous in x , has been designed for the attitudedynamics and the closed-loop solution tracks the constructedreference trajectory γ . We modify the control input online to account for the safety constraints. Concretely, the controlleris given by the quadratic program below: u ( x ) = arg min u ∈ R (cid:107) u − u nom (cid:107) s.t. L g b ( x ) u + L f b ( x ) + β ( b ( x )) ≥ . (17)This formulation reflects that the safety constraint has prior-ity over the tracking mission. Theorem 2.
For the attitude control system in (1) , thecontroller (17) renders the set C b ∩ C b forward invariant.Proof. The feasibility of the linear inequality constraint on u is guaranteed by Proposition 2 for every x ∈ C b ∩ C b .The solution to the quadratic program (17) has a closed-formsolution, given by the KKT condition [22], as u ( x ) = u nom + µL (cid:62) g b ( x ) (18)with µ = , if L g b u nom + β ( b ) + L f b ≥ , − L g b u nom − β ( b ) − L f b (cid:107) L g b (cid:107) , otherwise.This is derived from considering whether the constraintin (17) is activated or not and thus omitted here. Viewing L g b in (26), the property of χ ( · ) and Proposition 3, weobtain L g b = 0 if and only if when x ∈ C h,ξ ∩ C b , and,in the meanwhile, L g b u nom + β ( b ) + L f b ≥ is triviallysatisfied for x ∈ C h,ξ ∩ C b . Thus the solution is well-definedfor every x ∈ C b ∩ C b .The solution in (18) can be viewed as u ( x ) = v ( x ) + v ( v ( x )) v ( x ) (19)with v ( x ) = u nom ( x ) , v ( s ) = (cid:110) , if s ≥ s, if s< , v ( x ) = L g b u nom + β ( b ) + L f b , v ( x ) = − L (cid:62) g b (cid:107) L g b (cid:107) . For x ∈ ( C b \ C h,ξ ) ∩ C b , L g b ( x ) (cid:54) = 0 , we obtain v , v , v , v are locally Lipschitz continuous and thus u ( x ) is locallyLipschitz continuous in ( C b \ C h,ξ ) ∩ C b . Furthermore, for x ∈ C h,ξ ∩ C b , we have u ( x ) = u nom ( x ) and thus u ( x ) islocally Lipschitz continuous in C h,ξ ∩ C b .Now we show that the control input u ( x ) is continuousat the boundary between ( C b \ C h,ξ ) ∩ C b and C h,ξ ∩ C b .Assume a Cauchy sequence of points { x i } i =1 , , , ··· ⊂ ( C h \ C h,ξ ) ∩ C b , lim i →∞ x i = x ∈ ∂C h,ξ ∩ C b . Viewing L f b ( x ) = 0 , b ( x ) = L f b ( x ) + α ( b ( x )) , we get b ( x ) = α ( b ( x )) . Examining closer, it further derives b ( x ) = α (1) > . We then obtain lim i →∞ u ( x i ) = u ( x ) , viewing the closed-form solution in (18) and the factsthat u nom ( x i ) is bounded by definition, lim i →∞ L g b ( x i ) =0 , lim i →∞ L f b ( x i ) = 0 , lim i →∞ β ( b ( x i )) = β ( b ( x )) > . Thus u ( x ) is Lipschitz continuous in x for all C b ∩ C b ,which guarantees the existence and uniqueness of the systemsolution.For all x ∈ C b ∩ C b , we have ∂b∂x ( f ( x ) + g ( x ) u ) = b ≥ ∂b ∂x ( f ( x ) + g ( x ) u ) ≥ (20)6hus, the vector field lies in the tangent cone of set C b ∩ C b for all x ∈ C b ∩ C b . Applying Brezis version of Nagumo’sTheorem [21, Theorem 4] and noticing the locally Lipschitzvector field, we obtain the set C b ∩ C b is thus forwardinvariant. Remark 3.
Compared to the nonsmooth barrier function de-sign in [12], we formulate a smooth control barrier functionand thus avoid the nonsmooth analysis. This formulationcomes at the cost of conservativeness in terms of the setdifference between ∪ i ∈N S i and C Rh . Note that the conser-vativeness can be explicitly adjusted by choosing a proper δ .Based on this smooth control barrier function, we restore thesolvable optimization problem in (17). This set difference canbe viewed as a safety margin in many robotic applications. Remark 4.
Although [23] has studied the application ofbarrier functions in constrained attitude control problem, theproposed framework in this paper is generally more advan-tageous as 1)zeroing instead of reciprocal barrier function isused, which is well-defined even outside of the safety setand is guaranteed to be robust to model perturbations [15];2)here we deal with safety regions of arbitrary shape and thefeasibility to the online optimization is guaranteed.
Remark 5.
In Theorem 2 we guarantee the forward in-variance of the set C b ∩ C b instead of C h . This does notcause conservativeness. For any h ( x (0)) > , or equivalently, b ( x (0)) > , there always exists an extended class K func-tion α ( · ) such that b ( x (0)) = L f b ( x (0)) + α ( b ( x (0))) > .In this way, C b is constructed such that x (0) ∈ C b ∩ C b .V. S IMULATIONS
In this section, we demonstrate the favorable propertiesof the constructed reference trajectory and the designedzeroing control barrier function. The scenario is given asfollows: the inertia matrix of the rigid body is given by J = (cid:104) . . − . .
06 5 . . − .
03 0 .
01 0 . (cid:105) kg · m . The target attitude is set as R f = I , and the center points of the sampling cells are givenas R = exp(15 ◦ / ◦ × π [ e ] × ) , R = exp(30 ◦ / ◦ × π [ e ] × ) R , R = exp(30 ◦ / ◦ × π [0 , . , . × ) R ,and the initial attitude R = exp(10 ◦ / ◦ × π [ e ] × ) R .The radius of the cells is set as θ = 0 . rad (20 ◦ ) . Thesettling time is T = 40 s . Based on these data, we show theconstructed reference trajectory in red in Fig. 1.In what follows, we use the saturated controller from [2]as the nominal controller in (17): u nom = J ˜ R (cid:62) ˙ ω r + [ ˜ R (cid:62) ω r ] × J ˜ R (cid:62) ω r − k ( ˜ R − ˜ R (cid:62) ) ∨ − k tanh(˜ ω ) , (21)where ˜ R = R (cid:62) r R , ˜ ω ( t ) = ω − ˜ R (cid:62) ω r , R r , ω r are thereference orientation and reference angular velocity obtainedfrom the constructed trajectory γ , respectively, k , k > are tuning gains and tanh( · ) is the element-wise hyperbolictangent function. It is shown in [2] that the control law(21) achieves asymptotic convergence of the attitude trackingerror from almost all initial conditions. In the simulations, we augment the control sig-nal in (21) with an additive signal u add = 0 . ∗ (cid:0) sin(2 π t − ) , sin( π t − ) , − sin( π t − ) (cid:1) for the time inter-val t ∈ [20 , . This control signal simulates, for example,a human input to the system that could lead to a deviationfrom the reference trajectory and may even drive the statesout of the safe cells. The controller parameters are set as k = 0 . , k = 0 . . The parameters in the control barrierfunction are chosen as δ = 0 . , ξ = 0 . , α ( x ) = β ( x ) = x, χ ( x ) = (cid:40) ( x − + 1 , if x ≤ , , if x > . The simulation results are shown in the following. Fig. 4shows the trajectories in three cases: 1) no additive signal isapplied and the control barrier function exists (in blue); 2)additive signal is applied and control barrier function doesnot exist (in dark red); 3) additive signal is applied andcontrol barrier function exists (in yellow). It is shown thatwithout the additive control signal, the system trajectory issimilar to the reference trajectory in Fig. 1. However, whenthe additive signal exists and only the controller in (21)is applied, the state deviates from the previous trajectoryand runs out of the safety cells. Once the control barrierformulation is applied, the resulting trajectory remains in thesafety set. This is further supported by the time history of b ( x ) in Fig. 5. VI. C ONCLUSION
In this paper, we construct a C reference trajectoryon SO (3) and develop a safety certificate utilizing thecontrol barrier function formulation for constrained attitudecontrol problems, following the framework of our previouswork in [11]. To construct the reference trajectory, we firstdesign the controlling points for B´ezier curve generationon SO (3) , which is then time re-parametrized to satisfyboundary conditions. The reference trajectory is shown tobe C continuous, connect the initial and target orientations,and evolve within the predefined safe regions. Moreover, asmooth control barrier function is designed over a set ofoverlapping cells to circumvent the non-smooth analysis inprevious works. The safety certificate is given as a linearconstraint on the control input. This paper also provides aremedy to handle the states when the singularity of the linearconstraint occurs. A PPENDIX
We collect in this appendix all the results supporting thederivations and claims of the main part of the paper.From simplicity, denote the auxiliary variables q =( x , x , · · · , x ) , ω = ( ω , ω , ω ) and the elements in R i as (cid:32) x i x i x i x i x i x i x i x i x i (cid:33) , and let [ A ] i,j be the ( i, j ) th element of matrix A . Thus the state variable is rewritten as x = ( q, ω ) . As q isthe stacked vector of the rotation matrix R , we use q and R interchangeably to denote the attitude state in the following.We obtain that, for r i ( q ) defined in (9),7 a) Trajectory of x -axis. (b) Trajectory of y -axis. (c) Trajectory of z -axis. Fig. 4: Comparison of the trajectories of body-fixed axes in three cases.
Time t (sec) -0.500.511.52 b ( t ) safety barrier, no extra input extra input, no safety barrier extra input and safety barrier Fig. 5: Time histories of the control barrier function b ( x ) inthe three cases. (cid:20) ∂r i ∂x (cid:21) j = (cid:26) − x j + x ij , j = 1 , , · · · , , , j = 10 , , . (22)From h ( q ) given in (11), we further have ∂h∂x = ( ∂h∂q , ∂h∂ω ) with ∂h∂q = (cid:80) i ∈N ∂s ( r i ( q ) /(cid:15) ) ∂q := (cid:15) (cid:80) i ∈N dsd η i ∂ ( r i ( q )) ∂q , ∂h∂ω = (cid:80) i ∈N ∂s ( r i ( q ) /(cid:15) ) ∂ω = 0 , where η i ( q ) := r i ( q ) /(cid:15) for brevity.With f in (12), we further obtain L f h = ∂h∂x · f = 1 (cid:15) (cid:88) i ∈N dsd η i ∂ ( r i ( x )) ∂x · f = 1 (cid:15) (cid:88) i ∈N dsd η i ω (cid:62) e i − e i e i − e i e i − e i := 1 (cid:15) (cid:88) i ∈N dsd η i ω (cid:62) e i ( q ) , (23)where e ij,k ( q ) = [ R (cid:62) R i ] j,k for j, k = 1 , , .Similarly, we have L g h = ∂h∂x · g . Noticing that g in (12)and ∂h∂ω = 0 , we obtain L g h = ∂h∂q · × + ∂h∂ω · J − = 0 . Moreover, we can calculate that L g L f h = ∂L f h∂x · g = ∂L f h∂q · × + ∂L f h∂ω · J − = (cid:15) ∂ (cid:80) i ∈N dsd η i ω (cid:62) e i ( q ) ∂ω · J − . Note that η i ( q ) only relies on q , and thus L g L f h = 1 (cid:15) (cid:88) i ∈N dsd η i ∂ω (cid:62) e i ( q ) ∂ω · J − = 1 (cid:15) (cid:88) i ∈N dsd η i e i ( q ) (cid:62) J − (24) Proof of Proposition 3.
Noting that J is positive definite, x ∈ D := { x ∈ C h : L g L f h = 0 } implies that (cid:80) i ∈N dsd η i e i ( q ) (cid:62) = 0 . Further noticing that the state caneither be in one cell or in the intersection of two cells, weanalyze these two cases separately.1) If there exists a cell S j , j ∈ N such that r j ( x ) > , and r k ( x ) = 0 for all k ∈ N , k (cid:54) = j , then (cid:80) i ∈N dsd η i e i ( q ) = 0 ⇒ dsd η j e j ( q ) = 0 ⇒ e j ( q ) = 0 viewing the property of dsd η i . Considering the definitionof e i ( q ) in (23), we obtain R (cid:62) R i = I , i.e., R = R i ,which obviously lies inside C h .2) If there exist two cells S j , S k , j, k ∈ N such that R ∈ S j ∩ S k , i.e., r j > , r k > , then the condition (cid:80) i ∈N dsd η i e i ( q ) = 0 is equivalent to dsd η j e j ( q ) + dsd η k e k ( q ) = 0 . (25)As x lies in the intersection of two cells, we denote exp([ v j ] × ) := R (cid:62) R j , exp([ v k ] × ) := R (cid:62) R k for some v j , v k ∈ R . We can verify that e j ( q ) is parallel to v j ,and e k ( q ) is parallel to v k , respectively. Thus, in orderto fulfill condition (25), it suffices that v j (cid:107) v k , whichmeans that R lies on the geodesic path between R j and R k . Notice again that R ∈ S j ∩ S K , from Lemma2, we obtain that the system state x lies in the interiorof C h .Thus, the singular point set D is composed of all thecenter points of the sampling cells and certain points onthe geodesic path between R j , R k , where S j and S k areadjacent, j, k ∈ N . It can be checked that there exists a ξ > such that h ( R i ) = s (1) > ξ > and h ( R τ ) > ξ > for R τ = R i exp( τ log( R (cid:62) i R j )) , ≤ τ ≤ , and thus weobtain inf x ∈D h ( x ) ≥ ξ.
8n the following, we derive the explicit expressions for L f b, L f b , L g b that are used in the analysis and simula-tions in previous sections. From (14), we obtain L f b = dχdι ( h ( x ) ξ ) ∂h/ξ∂x · f = ξ dχdι L f h and L g b = ξ dχdι L g h = 0 , where ι := h ( x ) /ξ for brevity.As b ( x ) = L f b + α ( b ( x )) , we obtain L f b = ξ (cid:16) ξ d χdι ( L f h ) + dχdι L f h (cid:17) + dαdb ( b ( x )) L f b with L f h = (cid:15) (cid:80) i ∈N (cid:16) d sd η i ∂ η i ∂x ω (cid:62) e i ( q ) + dsd η i ∂ω (cid:62) e i ( q ) ∂x (cid:17) · f = (cid:15) (cid:80) i ∈N (cid:16) (cid:15) d sd η i ( ω (cid:62) e i ( q )) + dsd η i ∂ω (cid:62) e i ( q ) ∂x · f (cid:17) and L g b = L g L f b + dαdb ( b ( x )) L g b = L g L f b = 1 ξ (cid:18) ξ d χdι L g hL f h + dχdι L g L f h (cid:19) = 1 ξ dχdι L g L f h (26)with L g L f h given in (24).A CKNOWLEDGMENT
The authors thank Prof. Fatima Silva Leite for the inspiringdiscussion on the De Casteljau algorithm on SO (3) .R EFERENCES[1] D. E. Koditschek, “Application of a new Lyapunov function toglobal adaptive attitude tracking,” in
Proceedings of the 27th IEEEConference on Decision and Control, , 1988, pp. 63–68.[2] T. Lee, “Robust adaptive attitude tracking on SO (3) with an appli-cation to a quadrotor uav,” IEEE Transactions on Control SystemsTechnology , vol. 21, no. 5, pp. 1924–1930, 2012.[3] S. Berkane, A. Abdessameud, and A. Tayebi, “Hybrid global expo-nential stabilization on SO (3) ,” Automatica , vol. 81, pp. 279–285,2017.[4] H. B. Hablani, “Attitude commands avoiding bright objects and main-taining communication with ground station,”
Journal of Guidance,Control, and Dynamics , vol. 22, no. 6, pp. 759–767, 1999.[5] U. Lee and M. Mesbahi, “Feedback control for spacecraft reorientationunder attitude constraints via convex potentials,”
IEEE Transactionson Aerospace and Electronic Systems , vol. 50, no. 4, pp. 2578–2592,2014.[6] S. Kulumani and T. Lee, “Constrained geometric attitude control on SO (3) ,” International Journal of Control, Automation and Systems ,vol. 15, no. 6, pp. 2796–2809, 2017.[7] Q. Hu, B. Chi, and M. R. Akella, “Anti-unwinding attitude control ofspacecraft with forbidden pointing constraints,”
Journal of Guidance,Control, and Dynamics , pp. 1–13, 2018.[8] E. Frazzoli, M. Dahleh, E. Feron, and R. Kornfeld, “A randomizedattitude slew planning algorithm for autonomous spacecraft,” in
AIAAGuidance, Navigation, and Control Conference , 2001.[9] H. C. Kjellberg and E. G. Lightsey, “Discretized quaternion con-strained attitude pathfinding,”
Journal of Guidance, Control, andDynamics , vol. 38, no. 11, pp. 713–718, 2015.[10] J. D. Biggs and L. Colley, “Geometric attitude motion planningfor spacecraft with pointing and actuator constraints,”
Journal ofGuidance, Control, and Dynamics , vol. 39, no. 7, pp. 1672–1677,2016.[11] X. Tan, S. Berkane, and D. V. Dimarogonas, “Constrained attitudemaneuvers on SO(3) : Rotation space sampling, planning and low-level control,”
Automatica , vol. 112, p. 108659, 2020.[12] P. Glotfelter, J. Cort´es, and M. Egerstedt, “Nonsmooth barrier func-tions with applications to multi-robot systems,”
IEEE control systemsletters , vol. 1, no. 2, pp. 310–315, 2017.[13] W. S. Cortez and D. V. Dimarogonas, “Correct-by-design controlbarrier functions for Euler-Lagrange systems with input constraints,”in , 2020.[14] W. Xiao and C. Belta, “Control barrier functions for systems with highrelative degree,” in . IEEE, 2019, pp. 27–34. [15] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrierfunction based quadratic programs for safety critical systems,”
IEEETransactions on Automatic Control , vol. 62, no. 8, pp. 3861–3876,2016.[16] F. Bullo and A. D. Lewis,
Geometric control of mechanical systems:modeling, analysis, and design for simple mechanical control systems .Springer Science & Business Media, 2004, vol. 49.[17] F. Park and B. Ravani, “B´ezier curves in Riemannian manifolds and liegroups with kinematics applications,”
Journal of mechanical designs ,vol. 117, no. 1, pp. 36–40, 1995.[18] P. Crouch, G. Kun, and F. S. Leite, “The De Casteljau algorithm onLie groups and spheres,”
Journal of Dynamical and Control Systems ,vol. 5, no. 3, pp. 397–429, 1999.[19] S. R. Lindemann and S. M. LaValle, “Simple and efficient algorithmsfor computing smooth, collision-free feedback laws over given celldecompositions,”
The International Journal of Robotics Research ,vol. 28, no. 5, pp. 600–621, 2009.[20] F. Blanchini and S. Miani,
Set-Theoretic Methods in Control , ser.Systems & Control : Foundations & Applications. Birkh¨auser, 2015.[21] R. Redheffer, “The theorems of Bony and Brezis on flow-invariantsets,”
The American Mathematical Monthly , vol. 79, no. 7, pp. 740–747, 1972.[22] S. Boyd, S. P. Boyd, and L. Vandenberghe,
Convex optimization .Cambridge university press, 2004.[23] G. Wu and K. Sreenath, “Safety-critical and constrained geometriccontrol synthesis using control Lyapunov and control barrier func-tions for systems evolving on manifolds,” in . IEEE, 2015, pp. 2038–2044.. IEEE, 2015, pp. 2038–2044.