A POLYNOMIAL-TIME INTERIOR-POINT METHOD FOR CIRCULAR CONE PROGRAMMING BASED ON KERNEL FUNCTIONS

. We present an interior-point method based on kernel functions for circular cone optimization problems, which has been found useful for describing optimal design problems of optimal grasping manipulation for multi-ﬁngered robots. Since the well-known second order cone is a particular circular cone, we ﬁrst establish an invertible linear mapping between a circular cone and its corresponding second order cone. Then we develop a kernel function based interior-point method to solve circular cone optimization in terms of the corresponding second order cone optimization problem. We derive the complexity bound of the interior-point method and conclude that circular cone optimization is polynomial-time solvable. Finally we illustrate the performance of interior-point method by a real-world quadruped robot example of optimal contact forces taken from the literature [10].


1.
Introduction. Linear conic programming extends linear programming models from linearity to non-linearity and from convexity to non-convexity. Nemirovskii [21] points out that conic programming allows us to reveal rich structure which is usually possessed by a convex program and to exploit this structure in order to process the program efficiently. Up to year of 2000, linear conic optimization mainly focuses on optimization problems over symmetric cones. Güler [16] summarizes the classification of symmetric cones encompassing only five cones out of which three are of interests to optimization research. Renegar [25] points out that these cones are the nonnegative orthant of the Euclidean space R n which leads to linear programming, the Lorentz cone in R n , which leads to second order cone programming and the cone of positive semi-definite matrices that leads to semidefinite programming. These three types of conic programming are first unified into a framework of symmetric conic programming by Nesterov and Todd and studied extensively. One of the powerful tools for solving symmetric conic programming problems is well-known interior-point method (IPM). The reader is referred to [9,21,25] for more details of IPM. Recently an increasing amount attention has been paid to a class of kernel function based interior-point algorithm, which is useful for solving symmetric conic where µ is the static friction coefficient of the substrate. The second example is an optimal grasping manipulation for multi-fingered robots [8,17,18]. The grasping force of the i-th finger is subject to the inequality (1), geometrically, which represents a circular cone, in which the set of all forces is constrained and belong to a circular cone centred about the surface normal (see Figure 2) with half-angle β = tan −1 µ.

Figure 2. Circular cone at a grasp point
Although nonsymmetric cone programming problems are useful to discrete nonconvex quadratic problems and combinatorial NP-hard problems, solving nonsymmetric cone programming efficiently remains challenge. For nonsymmetric cone optimization, there are two difficulties that are encountered by the properties of non-self-dual in the constraint cone. The first one is that verification of an element belongs a copositive cone is an NP-complete problem. The second one is developing interior-point algorithm with polynomial-time computational complexity. Unlike symmetric conic optimization, there are less IPMs for solving nonsymmetric conic optimization. Recently, Nesterov [22] proposed a new theoretical IPM that is based on an extension of the ideas of self-scaled optimization to the non-symmetric conic optimization. He developed a 4n-self-concordant barrier for an n-dimensional pcone, which is an example of nonsymmetric cone. Matsukawa and Yoshise [20] proposed a primal barrier function Phase I algorithm for solving conic optimization problems over the closed convex cone K called doubly nonnegative cone. Skajaa and Ye [26] designed a homogeneous IPM for nonsymmetric convex conic optimization. All of these IPMs are developed based on self-concordant barrier functions for its corresponding cone.
We are motivated by both of the application of circular cone constraints for quadruped robot and multi-fingered robots and the challenge for designing IPMs for a particular nonsymmetric conic optimization problem. In this paper, we present an interior-point method based on kernel functions for circular cone optimization problems, which has been found useful for describing optimal design problems of optimal grasping manipulation for multi-fingered robots. Since the well-known second order cone is a particular circular cone, we first establish an invertible linear mapping between a circular cone and its corresponding second order cone. Then we develop a kernel function based interior-point method to solve circular cone optimization in terms of the corresponding second order cone optimization problem. We derive the complexity bound of the interior-point method and conclude that circular cone optimization is polynomial-time solvable. Finally we illustrate the performance of interior-point method by a real-world quadruped robot example of optimal contact forces taken from the literature [10].
The paper is organized as follows. In Section 2, we establish the geometric and algebraic relationship between the circular cone and the second order cone. Furthermore, based on the spectral factorization of second order cone, we extend it to the circular cone. In Section 3, we consider circular cone programming problems and establish its dual theorem. Then we introduce the optimality condition and central path for circular cone programming problems in Section 4. In Section 5, based on the invertible linear map established between the circular cone and the second order cone, we present an interior-point algorithm for circular cone programming in terms of that of second order cone programming. Section 6 implements our interior-point method by a real-world quadruped robot example of optimal contact forces taken from the literature [10]. Finally, we conclude and give future research in Section 7.
Some notations used throughout the paper are as follows. R n and R n + denote the space of n-dimensional real column vectors and the space of n-dimensional nonnegative real column vectors, respectively. Thus, x = (x 1 , . . . , x n ) T ∈ R n is viewed as a column vector in R n . For any x, y ∈ R n , the Euclidean inner product and norm are denoted x, y = x T y and x = √ x T x. Given a set S, we denote S • by the interior of S. Finally, I n denotes the n × n identity matrix.
2. Second order cone and circular cone. The second order cone in R n , also called the Lorentz cone, is defined to be where x 2:n = (x 2 , ...x n ) T and · 2 denotes the Euclidean norm. The second order cone programming (SOCP) has received considerable attentions from researchers because of its wide range of applications, the particular structure and algebraic properties. IPMs exploit the special structure of SOCP problems for polynomialtime efficience [1]. Given a rotation angle θ ∈ (0, π 2 ), we can define an n-dimensional circular cone that associates with a second order cone in R n .
Obviously, if the rotation angle θ = π 4 , then the circular cone is the second order cone. In other words, the second order cone is a special circular cone. Throughout the paper, we call L n an association cone with C n θ . Like the second order cone, a circular cone is a pointed closed convex cone having hyper-spherical sections orthogonal to its axis of revolution around which the cone is invariant to rotation.
Recall the definition of dual cone of a given cone K ⊆ R n , where ·, · : R n × R n → R denotes an Euclidean inner product on R n . Then, we can verify that We state the definition of symmetric cone as below.
Definition 2.1. Let K ⊆ R n be a proper cone, that is, a pointed closed convex cone with nonempty interior, K is called a symmetric cone, • if K is self-dual, that is, K = K * • K is homogeneous, i.e., for ∀x, y ∈ K • , there exists a invertible mapping G such that G(K) = K and y = G(x).
It is straightforward to verify that the second order cone is a symmetric cone. However, a circular cone is not self-dual, i.e., K = K * unless θ = π 4 . The geometric interpretation of a circular cone, its dual cone and associated second order cone are shown in Figure 3.  2.1. Algebraic representation between second order cone and circular cone. In this subsection, we study the algebraic relationship between a circular cone and its associated second order cone. We first recall the algebraic representation developed by Chen et al [12] between the second order cone L n and the circular cone C n θ . Then we recall the classical spectral factorization for the second order cone L n summarized by Alizadeh in [1]. Finally, we conclude with a spectral factorization for circular cone C n θ . For the detail of algebra of circular cone, please refer to [12,13,29].
Algebraic relationship between the second order cone and the circular cone is established in terms of an invertible linear mapping which plays a crucial rule in designing interior-point algorithm for circular cone programming.
Given a rotation angle θ, let where I n−1 is the n − 1 dimensional unit matrix. It is straightforward to verify that The following theorem states the basic algebraic relationship between the second order cone and the circular cone.  29]). Let L n and C n θ be defined as in (2) and (3), respectively. Then, we have Obviously, we have which yields The above expression (7) is equivalent to By Theorem 2.2, we obtain a new algebraic representation of the dual cone of a circular cone: Obviously, (C n θ ) * is also a circular cone with the rotation angleθ = π 2 − θ. Although the following Theorem 2.3 and Corollary 1 are quite simple, they are useful for designing a primal dual interior-point algorithm for circular cone programming. Theorem 2.3. For any x ∈ C n θ and s ∈ (C n θ ) * , there existsx ∈ L n ands ∈ L n such that Obviously, we have the following corollary.

2.2.
Spectral factorization associated with second order cone and circular cone. In this subsection we recall the spectral factorization of a vector associated with the second order cone L n . Followed by the spectral factorization associated with the circular cone C n θ [29], we specialize the spectral factorizations for circular cone C n θ and its dual cone (C n θ ) * . Using these spectral factorizations, we construct the corresponding barrier functions for C n θ and (C n θ ) * , which are used to design the search direction for a primal dual interior-point algorithm in Section 3.
Alizadeh et al [1] pointed out that the spectral factorization is the special case of the so-called Euclidean Jordan algebra. Given x = (x 1 , x 2:n ) T ∈ R × R n−1 . Then x can be decomposed as where λ 1 , λ 2 and u (1) , u (2) are the spectral values and associated spectral vectors of x given by for i = 1, 2, with w being any vector in R n−1 satisfying w = 1. If x 2:n = 0, the decomposition (9) is unique.
The following theorem is the spectral factorization for a vector x associated with circular cone C n θ .

Theorem 2.4 ([29]
). Given a circular C n θ with θ ∈ (0, π 2 ), for any x ∈ R n , one has where If x 2:n = 0, then w = x 2:n x 2:n and otherwise, then w can be any vector in R n−1 satisfying w = 1.
Note that the above spectral factorization is hold for all x ∈ R n . Throughout the paper, we focus on the C n θ and (C n θ ) * . It is easy to have the next corollary.
3. Circular cone programming and duality theorem. In this paper, we consider the following circular cone programming as follows: Its dual problem becomes where K ⊆ R n and K * ⊆ R n are the Cartesian product of circular cones, i.e., with K j = C nj θj for each j, and n = N j=1 n j . Furthermore, we partition the vectors x, s, c, and the matrix A accordingly as x = (x 1 ; x 2 ; . . . ; x N ), s = (s 1 ; s 2 ; . . . ; s N ) with x j ∈ K j , s j ∈ (K j ) * , c = (c 1 ; c 2 ; . . . ; c N ) with c j ∈ R nj , and A = (A 1 , . . . , A N ) with A j ∈ R m×nj , and b ∈ R m . K j is a circular cone and (K j ) * is its dual cone. Assume that m ≤ n and rank(A) = m. The simplest case for (CCP) is that N = 1, that is the case we consider in Section 5 where we take that K = C n θ . As a matter of fact, it is easy to verify that (14) is the Lagrangian dual (13), using the fact that min{u T x : x ∈ K} = 0 if u ∈ K * and u → −∞, otherwise. We can also check the following weak duality theorem. Proof. Notice that with inequality holds. Remark 1. In (13) and (14), if one replaces C n θ with the nonnegative orthant, (13) and (14) becomes the standard form of linear programming problems. In that case, we have strong duality theorem. However, for circular cone programming, such property no longer holds.
In order to develop an interior-point algorithm, we need the following theorem [5].
Theorem 3.2 (Dual Theorem). Consider a circular cone programming (13) problem and its dual problem (14). The following results hold.
(a) The dual problem is conic, and the dual problem of dual is the primal.   (14) is bounded above and strictly feasible (i.e., A T y + s = c and s ∈ (K * ) • for some (y, s) ), then the primal (14) achieves optimal solution and there is no duality gap. (d) Assume that at least one of the problems (13), (14) is bounded and strictly feasible. Then a primal-dual feasible pair (x, y, s) is a pair of optimal solutions to the respective problems if and only if A useful consequence is the following Corollary 3 (Strong Duality Theorem). Assume that both (13) and (14) are strictly feasible. Then both problems achieve optimality with no duality gap and (d.1), (d.2) are necessary and sufficient for the optimality of a primal-dual feasible pair.

Optimality conditions and central path.
Based on the strong duality theorem of circular cone programming problems, given a rotation angle θ, we consider the simplest case of circular cone programming problems and discuss the optimality conditions as follows.
(CCP) min Its dual problem is Solving the primal dual pair of circular cone programming problems (16) and (17) is equivalent to solve the following optimality conditions We assume that both (CCP) and (CCD) satisfy the interior-point condition, i.e., there exists (x 0 , y 0 , s 0 ) such that Ax 0 = b, x 0 ∈ (C n θ ) • , A T y 0 +s 0 = c, s 0 ∈ ((C n θ ) * ) • . Now we use both the invertible linear mapping and its inverse mapping to transform the feasible sets of (CCP) and (CCD), respectively. x ∈ {x ∈ R n | Ax = b, x ∈ C n θ } and s ∈ {x ∈ R n | A T y + s = c, s ∈ (C n θ ) * }, respectively, letx = A θ x ands = A −1 θ s. Moreover, letÃ = AA −1 θ andc = A −1 θ c. We have the following associated primal dual pair of second order cone programming problems as follows.

YANQIN BAI, PENGFEI MA AND JING ZHANG
The optimality conditions for (20) and (21) are as follows.

5.
Interior-point algorithm for circular cone programming. In this section, we present an interior-point algorithm for solving (13) and (14) and derive the complexity bound.
We first recall an existing interior-point algorithm for SOCP presented by Bai et al in [3] and [4]. Then we develop a corresponding interior-point algorithm for nonsymmetric circular cone programming.
We associate a given kernel function ψ(t) to a real-valued barrier function Ψ(x) on L n , based on the spectral factorization associated with the second order cone L n as follows.
Primal-dual interior-point algorithm for second order cone programming

5.2.
Kernel function based interior-point algorithm for CCP. Now let us present a primal-dual interior-point algorithm for nonsymmetric circular cone programming problems. Given a rotation angle θ, we use the same initial parameters such as τ, , etc, and we use the same barrier function Φ(x, s, µ). The framework of primal-dual interior-point algorithm are as follows.

5.3.
Analysis of the interior-point algorithm for circular cone programming. Because the core of primal-dual interior-point algorithm for CCP-CCD (5.2) is the same as that of primal-dual interior-point algorithm for SOCP (5.1), the analysis of the computational complexity of (5.2) is similar to that of (5.1). Here we briefly recall the analysis of primal-dual interior-point algorithm for SOCP.
Furthermore, the iteration bound for the interior-point algorithm for SOCP states as follows: . Choose a suitable kernel function ψ(t). The iteration bound of the algorithm 5.2 is given by Note that this iteration bound is completely determined by the (eligible) kernel function ψ(t), the parameters η and τ in the algorithm. Also note that when 2N is replaced by n. Therefore, we conclude that the complexity bound for nonsymmetric circular cone programming is as below.
Theorem 5.3. Algorithm 5.2 has the iteration bounds O( √ n log n) log n for largeupdate methods and O( √ n) log n for small-update methods for solving nonsymmetirc circular cone programming (13) and (14), respectively.
In [3], it says that a suitable kernel function is chosen and (20) and (21) are solved by algorithm 5.1, the algorithm has an iteration bound O( √ n log n) log n for for large-update method and O( √ n) log n for small-update method, respectively. Therefore, for nonsymmetirc circular cone programming (13) and (14), it also has same iteration bounds as SOCP (20) and (21).
Remark 2. If we transform circular cone programming (13) and (14) into SOCP (20) and (21), the interior-point algorithms for solving SOCP can be used to solve (20) and (21). 6. Numerical tests. In this section, we use a practical numerical example to illustrate the performance of our interior-point method for solving CCP problems described in Section 5. The numerical testing is carried out on a Lenovo PC (2.6GHz, 4.00GB of RAM) with the use of Matlab R2012a. The parameters are taken as = 10 −5 , τ = 1, η = 0.7. The kernel functions ψ(t) are taken as ψ 1 (t) = t 2 −1 2 − log t, t > 0 and ψ 2 (t) = t 2 −1 2 + t −2 −1 2 , t > 0, respectively. Our example is for quadruped robot appeared in [10]. Example 1. We consider a force optimization problem for quadruped robot in Figure 1. Assume that each leg contacts one point with the ground. The quadruped robot has four contact points, which have positions p i ∈ R 3 (in the global coordinate system), for i = 1, 2, 3, 4. The force applied at a contact point p i is denoted by f i ∈ R 3 , and is given in a local coordinate system. We denote its components as f i = (f ni , f oi , f ti ) T . Thus f ni is the normal (inward) force applied at contact point p i , and (f oi , f ti ) is the tangential force applied at contact point p i . The contact force f i satisfy the following circular cone constraints where µ i = tan θ i is the friction coefficient at contact point p i . The above inequalities are equivalent to f i ∈ C 3 θi , i = 1, 2, 3, 4. The quadruped robot satisfies the force equilibrium condition and torque equilibrium condition. Let Q i be the 3 × 3 matrix that transforms forces in the local coordinate system at p i into the global coordinate system. Thus, the force at contact point p i is Q i f i (in the global coordinate system). The force equilibrium condition is where f ext ∈ R 3 is the total external force that acts on the quadruped robot (in the global coordinate system). The torque applied to the object by the force at contact point p i is given by p i ⊗ Q i f i (in the global coordinate system). The torque equilibrium condition is where τ ext ∈ R 3 is the total external torque that acts on the object (in the global coordinate system). We define the contact force vector f ∈ R 12 as f = (f 1 , f 2 , f 3 , f 4 ), contact matrices A i ∈ R 6×3 as and overall contact matrix A = [A 1 , A 2 , A 3 , A 4 ] ∈ R 6×12 . We collect the external force and torque into a single external wrench ω ext = (−f ext , −τ ext ) ∈ R 6 . Thus, we can write the the equilibrium conditions (28) and (29) as At each contact point i, a direction is c i ∈ R 3 . A set of contact force vectors f i ∈ C 3 θi can be found as a work does as small as possible in c i , while satisfying the torque equilibrium conditions (30). We define the direction vector c ∈ R 12 as c = (c 1 , c 2 , c 3 , c 4 ). The force optimization problem for quadruped robot can be expressed as the following circular cone programming For problem (31), we have the following data: and y = (0, 0, 0, 0, 0, 0).
We choose the rotation angles as θ = π 12 , π 10 , π 8 , π 6 , π 4 , respectively. The coefficients of friction corresponding with rotation angles are tan π 12 , tan π 10 , tan π 8 , tan π 6 , tan π 4 , respectively. The numerical results are summarized in Tables 1 and 2. Pobj, Dobj, Itn and DG denote the primal objective function value, the dual objective function value, the iteration number and the duality gap. Table 1 shows the objective function value, the iteration number and the duality gap. Because the CPU time is almost negligible, we do not report the CPU time. The duality gap implies that desired accuracy is achieved. The simple practical example is shown to justify the efficiency of the proposed algorithm. Observing Table 1, we find that the objective function value gets greater as the angle gets smaller. By (3), the reason for this phenomenon is that the feasible region of the problem (31) gets smaller as the angle gets smaller. On the other hand, the algorithm based on ψ 2 (t) is superior to the algorithm based on ψ 1 (t). Table 2 shows the duality gaps of each iteration based on ψ 2 (t) for θ = π 10 . Observing Table 2, we find that the duality gap between odd iteration and its adjacent odd iteration is reduced with a factor almost 10 after the third iteration.  Further, to verify whether iterations required comply with the complexity result in Section 5, we conduct tests in different errors = 10 −2 , = 10 −8 and = 10 −12 and do not change the other parameters. The numerical results are summarized in Tables 3, 4 and 5. Observing Tables 3, 4 and 5, we find that the iterations required in practice agree with the complexity result in Section 5.
Example 2. Boyd et al in [10] proposed a force optimization problem for multifingered robots, where the grasp force is minimized. Assume that the object is grasped at the M contact points. Then the force optimization problem is formulated as follows: . . , f M ) ∈ R 3M and A ∈ R 6×M . We modify the above example slightly to minimize the work caused by force. For the multi-fingered robots that have M contact points with the object, the optimization problem can be expressed into a circular cone optimization problem    Therefore, it can be solved by our interior-point algorithm efficiently after giving the fixed data of all parameters. 7. Conclusions. In this paper, we have presented an interior-point method based on kernel functions for circular cone optimization problems which can be used to describe optimal design problems of optimal grasping manipulation for multi-fingered robots. Then we have developed a kernel function based interior-point method to solve circular cone optimization in terms of the corresponding second order cone optimization problem. We have concluded that circular cone optimization is polynomial-time solvable. This inspirits us to develop an especial interior-point algorithm for solving general nonsymmetric cone optimization problems, where the nonsymmetric cone constraint may be decomposed as an union of serval circular cones.