Second-Order Cone Programming Algorithm - MATLAB & Simulink (original) (raw)

Definition of Second-Order Cone Programming

A second-order cone programming problem has the form

subject to the constraints

f, x, b,beq, lb, and ub are vectors, and A and Aeq are matrices. For each_i_, the matrix_A_soc(i), the vectors _b_soc(i) and_d_soc(i), and the scalar γ(i) are in a second-order cone constraint that you create using secondordercone.

In other words, the problem has a linear objective function and linear constraints, as well as a set of second-order cone constraints of the form ‖Asoc(i)⋅x−bsoc(i)‖≤dsocT(i)⋅x−γ(i).

coneprog Algorithm

The coneprog solver uses the algorithm described in Andersen, Roos, and Terlaky [1]. This method is an interior-point algorithm similar to the Interior-Point linprog Algorithm.

Standard Form

The algorithm starts by placing the problem in standard form. The algorithm adds nonnegative slack variables so that the problem has the form

subject to the constraints

The solver expands the sizes of the linear coefficient vector_f_ and linear constraint matrix A to account for the slack variables.

The region K is the cross product of Lorentz cones Equation 1 and the nonnegative orthant. To convert each convex cone

to a Lorentz cone Equation 1, create a column vector of variables _t_1,_t_2, …,t n+1:

Here, the number of variables n for each cone_i_ is the number of rows in_A_soc(i). By its definition, the variable vector t satisfies the inequality

Equation 1 is the definition of a Lorentz cone in (n+1) variables. The variables t appear in the problem in place of the variables_x_ in the convex region K.

Internally, the algorithm also uses a rotated Lorentz cone in the reformulation of cone constraints, but this topic does not address that case. For details, see Andersen, Roos, and Terlaky [1].

When adding slack variables, the algorithm negates variables, as needed, and adds appropriate constants so that:

Dual Problem

The dual cone is

The dual problem is

such that

for some

A dual optimal solution is a point (y,s) that satisfies the dual constraints and maximizes the dual objective.

Homogeneous Self-Dual Formulation

To handle potentially infeasible or unbounded problems, the algorithm adds two more variables τ and κ and formulates the problem as homogeneous (equal to zero) and self-dual.

Ax−bτ=0ATy+s−fτ=0−fTx+bTy−κ=0 (2)

along with the constraints

Here, K¯ is the cone K adjoined with the nonnegative real line, which is the space for (x;τ). Similarly K¯* is the cone K* adjoined with the nonnegative real line, which is the space for (s;κ). In this formulation, the following lemma shows that τ is the scaling for feasible solutions, and κ is the indicator of an infeasible problem.

Lemma ([1] Lemma 2.1)

Let (x, τ, y,s, κ) be a feasible solution of Equation 2 along with the constraints in Equation 3.

In summary, for feasible problems, the variable τ scales the solution between the original standard form problem and the homogeneous self-dual problem. For infeasible problems, the final iterate (x, y, s,τ, κ) provides a certificate of infeasibility for the original standard form problem.

Start Point

The start point for the iterations is the feasible point:

Central Path

The algorithm attempts to follow the central path, which is the parameterized solution to the following equations for_γ_ decreasing from 1 toward 0.

Ax−bτ=γ(Ax0−bτ0)ATy+s−cτ=γ(ATy0+s0−fτ0)−fTx+bTy−κ=γ(−fTx0+bTy0−κ0)XSe=γμ0eτκ=γμ0. (4)

The central path begins at the start point and ends at an optimal solution to the homogeneous self-dual problem.

Andersen, Roos, and Terlaky [1] show in Lemma 3.1 that the complementarity condition_xTs_ = 0, where_x_ and s are in a product of Lorentz cones L, is equivalent to the condition

for every cone i. Here_Xi_ = mat(xi),xi is the variable associated with the Lorentz cone i,Si = mat(si), and_ei_ is the unit vector [1,0,0,…,0] of the appropriate dimension. This discussion shows that the central path satisfies the complementarity condition at its end point.

Search Direction

To obtain points near the central path as the parameter γ decreases from 1 toward 0, the algorithm uses Newton's method. The variables to find are labeled (x, τ,y, s, κ). Let_dx_ represent the search direction for the x variables, and so on. Then the Newton step solves the following linear system, derived from Equation 4.

The algorithm obtains its next point by taking a step in the_d_ direction.

for some step α∈[0,1].

For both numerical stability and accelerated convergence, the algorithm scales the step according to a suggestion in Nesterov and Todd [8]. Also, the algorithm corrects the step according to a variant of Mehrotra's predictor-corrector [7]. (For further details, see Andersen, Roos, and Terlaky [1].)

Step Solver Variations

The preceding discussion relates to the LinearSolver option with the value 'augmented' specified. The solver has other values that change the step calculation to suit different types of problems.

Iterative Display and Stopping Conditions

At each iteration k, the algorithm computes three relative convergence measures:

You can view these three statistics at the command line by specifying iterative display.

options = optimoptions('coneprog','Display','iter');

All three should approach zero when the problem is feasible and the solver converges. For a feasible problem, the variable_κk_ approaches zero, and the variable τk approaches a positive constant.

One stopping condition is somewhat related to the gap infeasibility. The stopping condition is when the following optimality measure decreases below the optimality tolerance.

| Optimalityk=|fTxk−bTyk|τk+|bTyk|=|fTxk/τk−bTyk/τk|1+|bTyk/τk|. | (5) | | ----------------------------------------------------------------- | ------- |

This statistic measures the precision of the objective value.

The solver also stops and declares the problem to be infeasible under the following conditions. The three relative infeasibility measures are less than_c_ = ConstraintTolerance, and

If bTyk > 0, then the solver declares that the primal problem is infeasible. If fTxk < 0, then the solver declares that the dual problem is infeasible.

The algorithm also stops when

and

In this case, coneprog reports that the problem is numerically unstable (exit flag -10).

The remaining stopping condition occurs when at least one infeasibility measure is greater than ConstraintTolerance and the computed step size is too small. In this case, coneprog reports that the search direction became too small and no further progress could be made (exit flag -7).

References

[1] Andersen, E. D., C. Roos, and T. Terlaky. On implementing a primal-dual interior-point method for conic quadratic optimization. Math. Program., Ser. B 95, pp. 249–277 (2003). https://doi.org/10.1007/s10107-002-0349-3

[2] Andersen, K. D.A modified schur-complement method for handling dense columns in interior-point methods for linear programming. ACM Transactions on Mathematical Software (TOMS), 22(3):348–356, 1996.

[3] Ben-Tal, Aharon, and Arkadi Nemirovski. Convex Optimization in Engineering: Modeling, Analysis, Algorithms. (1998).

[4] Goldfarb, D. and K. Scheinberg. A product-form cholesky factorization method for handling dense columns in interior point methods for linear programming. Mathematical Programming, 99(1):1–34, 2004.

[5] Goldfarb, D. and K. Scheinberg. Product-form cholesky factorization in interior point methods for second-order cone programming. Mathematical Programming, 103(1):153–179, 2005.

[6] Luo, Zhi-Quan, Jos F. Sturm, and Shuzhong Zhang. Duality and Self-Duality for Conic Convex Programming. (1996).

[7] Mehrotra, Sanjay. “On the Implementation of a Primal-Dual Interior Point Method.” SIAM Journal on Optimization 2, no. 4 (November 1992): 575–601. https://doi.org/10.1137/0802028.

[8] Nesterov, Yu. E., and M. J. Todd. “Self-Scaled Barriers and Interior-Point Methods for Convex Programming.”Mathematics of Operations Research 22, no. 1 (February 1997): 1–42. https://doi.org/10.1287/moor.22.1.1.

See Also

coneprog | secondordercone

Topics