8.2 Problems with Linear Objectives

cpl(c, F[, G, h[, dims[, A, b[, kktsolver]]]])

Solves a convex optimization problem with a linear objective

minimize  cTx
subject to  fk(x) ≤ 0, k = 0,...,m - 1
          Gx ≼ h
          Ax = b.

c is a real single-column dense matrix.

F is a function that evaluates the nonlinear constraint functions. It must handle the following calling sequences.

The linear inequalities are with respect to a cone C defined as a Cartesian product of a nonnegative orthant, a number of second-order cones, and a number of positive semidefinite cones:

C = C0 × C1 × ⋅⋅⋅× CM × CM+1 × ⋅⋅⋅× CM+N

with

          l                                            rk-1                                          {            tk}
C0 = {u ∈ R | uk ≥ 0, k = 1,...,l}, Ck+1 = {(u0,u1) ∈ R×R    | u0 ≥ ∥u1∥2}, k = 0,...,M - 1, Ck+M+1  =  vec(u) | u ∈ S+ , k = 0,...,N - 1.

Here vec(u) denotes a symmetric matrix u stored as a vector in column major order.

The arguments h and b are real single-column dense matrices. G and A are real dense or sparse matrices. The default values for A and b are sparse matrices with zero rows, meaning that there are no equality constraints. The number of rows of G and h is equal to

       M∑-1    N∑- 12
K = l+     rk +    tk.
        k=0     k=0

The columns of G and h are vectors in

                         2         2
Rl × Rr0 × ⋅⋅⋅× RrM -1 × Rt0 × ⋅⋅⋅× RtN-1,

where the last N components represent symmetric matrices stored in column major order. The strictly upper triangular entries of these matrices are not accessed (i.e., the symmetric matrices are stored in the ’L’-type column major order used in the blas and lapack modules).

The argument dims is a dictionary with the dimensions of the cones. It has three fields.

dims[’l’]:
l, the dimension of the nonnegative orthant (a nonnegative integer).
dims[’q’]:
[r0,,rM-1], a list with the dimensions of the second-order cones (positive integers).
dims[’s’]:
[t0,,tN-1], a list with the dimensions of the positive semidefinite cones (nonnegative integers).

The default value of dims is {’l’: h.size[0], ’q’: [], ’s’: []}, i.e., the default assumption is that the linear inequalities are componentwise inequalities.

The role of the optional argument kktsolver is explained in section 8.4.

cpl() returns a dictionary that contains the result and information about the accuracy of the solution. The most important fields have keys ’status’, ’x’, ’snl’, ’sl’, ’y’, ’znl’, ’zl’. The possible values of the ’status’ key are:

’optimal’
In this case the ’x’ entry of the dictionary is the primal optimal solution, the ’snl’ and ’sl’ entries are the corresponding slacks in the nonlinear and linear inequality constraints, and the ’znl’, ’zl’ and ’y’ entries are the optimal values of the dual variables associated with the nonlinear inequalities, the linear inequalities, and the linear equality constraints. These vectors approximately satisfy the Karush-Kuhn-Tucker (KKT) conditions
       T      T    T
c+Df (x) znl+G  zl+A   y = 0,   f(x)+snl = 0, k = 1,...,m,    Gx+sl =  h,   Ax = b,

snl ≽ 0,   sl ≽ 0,  znl ≽ 0,   zl ≽ 0,  sTnlznl + sTl zl = 0.

’unknown’
This indicates that the algorithm terminated before a solution was found, due to numerical difficulties or because the maximum number of iterations was reached. The ’x’, ’snl’, ’sl’, ’y’, ’znl’ and ’zl’ contain the iterates when the algorithm terminated.

The other entries in the output dictionary describe the accuracy of the solution. The entries ’primal objective’, ’dual objective’, ’gap’, and ’relative gap’, give the primal objective cTx, the dual objective calculated as

cT x+ zTnlf(x)+ zTl (Gx - h)+ yT (Ax - b),

the duality gap

 T      T
snlznl + sl zl,

and the relative gap. The relative gap is defined as

      gap                                    gap
- primal objective if primal objective < 0,   dual objective if dual objective > 0,

and None otherwise. The entry with key ’primal infeasibility’ gives the residual in the primal constraints,

----∥(f(x)+-snl,Gx-+-sl --h,Ax---b)∥2----
max{1,∥(f(x0) +1,Gx0 + 1- h,Ax0 - b)∥2}

where x0 is the point returned by F(). The entry with key ’dual infeasibility’ gives the residual

 ∥c+ Df(x)Tznl + GT zl + AT y∥2
max-{1,∥c+-Df-(x-)T1-+-GT1∥-} ≤ ϵfeas.
               0          2

cpl() requires that the problem is strictly primal and dual feasible and that

                  ([ ∑m - 1   2        T                        T ])
rank(A) = p,   rank     k=0 zk∇ fk(x) A    ∇f0(x)  ⋅⋅⋅∇fm -1(x)  G     = n,

for all x and all positive z.

Example: floor planning

This example is the floor planning problem of section 8.8.2 in the book Convex Optimization:

minimize   W + H
subject to Amin,k∕hk - wk ≤ 0, k = 1,...,5
           x1 ≥ 0, x2 ≥ 0, x4 ≥ 0
           x1 + w1 + ρ ≤ x3, x2 + w2 + ρ ≤ x3, x3 + w3 + ρ ≤ x5, x4 + w4 + ρ ≤ x5, x5 + w5 ≤ W
           y2 ≥ 0, y3 ≥ 0, y5 ≥ 0
           y2 + h2 + ρ ≤ y1, y1 + h1 + ρ ≤ y4, y3 + h3 + ρ ≤ y4, y4 + h4 ≤ H, y5 + h5 ≤ H
           hk∕γ ≤ wk ≤ γhk, k = 1,...,5.

This problem has 22 variables

                  5         5          5          5
W,    H,     x ∈ R ,   y ∈ R ,    w ∈ R ,   h ∈ R ,

5 nonlinear inequality constraints, and 26 linear inequality constraints. The code belows defines a function floorplan() that solves the problem by calling cp(), then applies it to 4 instances, and creates a figure.

import pylab  
from cvxopt import solvers, matrix, spmatrix, mul, div  
 
def floorplan(Amin):  
 
    #     minimize    W+H  
    #     subject to  Amink / hk <= wk, k = 1,..., 5  
    #                 x1 >= 0,  x2 >= 0, x4 >= 0  
    #                 x1 + w1 + rho <= x3  
    #                 x2 + w2 + rho <= x3  
    #                 x3 + w3 + rho <= x5  
    #                 x4 + w4 + rho <= x5  
    #                 x5 + w5 <= W  
    #                 y2 >= 0,  y3 >= 0,  y5 >= 0  
    #                 y2 + h2 + rho <= y1  
    #                 y1 + h1 + rho <= y4  
    #                 y3 + h3 + rho <= y4  
    #                 y4 + h4 <= H  
    #                 y5 + h5 <= H  
    #                 hk/gamma <= wk <= gamma*hk,  k = 1, ..., 5  
    #  
    # 22 Variables W, H, x (5), y (5), w (5), h (5).  
    #  
    # W, H:  scalars; bounding box width and height  
    # x, y:  5-vectors; coordinates of bottom left corners of blocks  
    # w, h:  5-vectors; widths and heigths of the 5 blocks  
 
    rho, gamma = 1.0, 5.0   # min spacing, min aspect ratio  
 
    # The objective is to minimize W + H.  There are five nonlinear  
    # constraints  
    #  
    #     -wk + Amink / hk <= 0,  k = 1, ..., 5  
 
    c = matrix(2*[1.0] + 20*[0.0])  
 
    def F(x=None, z=None):  
        if x is None:  return 5, matrix(17*[0.0] + 5*[1.0])  
        if min(x[17:]) <= 0.0:  return None  
        f = -x[12:17] + div(Amin, x[17:])  
        Df = matrix(0.0, (5,22))  
        Df[:,12:17] = spmatrix(-1.0, range(5), range(5))  
        Df[:,17:] = spmatrix(-div(Amin, x[17:]**2), range(5), range(5))  
        if z is None: return f, Df  
        H = spmatrix( 2.0* mul(z, div(Amin, x[17::]**3)), range(17,22), range(17,22) )  
        return f, Df, H  
 
    G = matrix(0.0, (26,22))  
    h = matrix(0.0, (26,1))  
    G[0,2] = -1.0                                       # -x1 <= 0  
    G[1,3] = -1.0                                       # -x2 <= 0  
    G[2,5] = -1.0                                       # -x4 <= 0  
    G[3, [2, 4, 12]], h[3] = [1.0, -1.0, 1.0], -rho     # x1 - x3 + w1 <= -rho  
    G[4, [3, 4, 13]], h[4] = [1.0, -1.0, 1.0], -rho     # x2 - x3 + w2 <= -rho  
    G[5, [4, 6, 14]], h[5] = [1.0, -1.0, 1.0], -rho     # x3 - x5 + w3 <= -rho  
    G[6, [5, 6, 15]], h[6] = [1.0, -1.0, 1.0], -rho     # x4 - x5 + w4 <= -rho  
    G[7, [0, 6, 16]] = -1.0, 1.0, 1.0                   # -W + x5 + w5 <= 0  
    G[8,8] = -1.0                                       # -y2 <= 0  
    G[9,9] = -1.0                                       # -y3 <= 0  
    G[10,11] = -1.0                                     # -y5 <= 0  
    G[11, [7, 8, 18]], h[11] = [-1.0, 1.0, 1.0], -rho   # -y1 + y2 + h2 <= -rho  
    G[12, [7, 10, 17]], h[12] = [1.0, -1.0, 1.0], -rho  #  y1 - y4 + h1 <= -rho  
    G[13, [9, 10, 19]], h[13] = [1.0, -1.0, 1.0], -rho  #  y3 - y4 + h3 <= -rho  
    G[14, [1, 10, 20]] = -1.0, 1.0, 1.0                 # -H + y4 + h4 <= 0  
    G[15, [1, 11, 21]] = -1.0, 1.0, 1.0                 # -H + y5 + h5 <= 0  
    G[16, [12, 17]] = -1.0, 1.0/gamma                   # -w1 + h1/gamma <= 0  
    G[17, [12, 17]] = 1.0, -gamma                       #  w1 - gamma * h1 <= 0  
    G[18, [13, 18]] = -1.0, 1.0/gamma                   # -w2 + h2/gamma <= 0  
    G[19, [13, 18]] = 1.0, -gamma                       #  w2 - gamma * h2 <= 0  
    G[20, [14, 18]] = -1.0, 1.0/gamma                   # -w3 + h3/gamma <= 0  
    G[21, [14, 19]] = 1.0, -gamma                       #  w3 - gamma * h3 <= 0  
    G[22, [15, 19]] = -1.0, 1.0/gamma                   # -w4  + h4/gamma <= 0  
    G[23, [15, 20]] = 1.0, -gamma                       #  w4 - gamma * h4 <= 0  
    G[24, [16, 21]] = -1.0, 1.0/gamma                   # -w5 + h5/gamma <= 0  
    G[25, [16, 21]] = 1.0, -gamma                       #  w5 - gamma * h5 <= 0.0  
 
    # solve and return W, H, x, y, w, h  
    sol = solvers.cpl(c, F, G, h)  
    return  sol[’x’][0], sol[’x’][1], sol[’x’][2:7], sol[’x’][7:12], sol[’x’][12:17], sol[’x’][17:]  
 
pylab.figure(facecolor=’w’)  
pylab.subplot(221)  
Amin = matrix([100., 100., 100., 100., 100.])  
W, H, x, y, w, h =  floorplan(Amin)  
for k in xrange(5):  
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]],  
               [y[k], y[k]+h[k], y[k]+h[k], y[k]], facecolor = ’#D0D0D0’)  
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))  
pylab.axis([-1.0, 26, -1.0, 26])  
pylab.xticks([])  
pylab.yticks([])  
 
pylab.subplot(222)  
Amin = matrix([20., 50., 80., 150., 200.])  
W, H, x, y, w, h =  floorplan(Amin)  
for k in xrange(5):  
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]],  
               [y[k], y[k]+h[k], y[k]+h[k], y[k]], ’facecolor = #D0D0D0’)  
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))  
pylab.axis([-1.0, 26, -1.0, 26])  
pylab.xticks([])  
pylab.yticks([])  
 
pylab.subplot(223)  
Amin = matrix([180., 80., 80., 80., 80.])  
W, H, x, y, w, h =  floorplan(Amin)  
for k in xrange(5):  
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]],  
               [y[k], y[k]+h[k], y[k]+h[k], y[k]], ’facecolor = #D0D0D0’)  
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))  
pylab.axis([-1.0, 26, -1.0, 26])  
pylab.xticks([])  
pylab.yticks([])  
 
pylab.subplot(224)  
Amin = matrix([20., 150., 20., 200., 110.])  
W, H, x, y, w, h =  floorplan(Amin)  
for k in xrange(5):  
    pylab.fill([x[k], x[k], x[k]+w[k], x[k]+w[k]],  
               [y[k], y[k]+h[k], y[k]+h[k], y[k]], ’facecolor = #D0D0D0’)  
    pylab.text(x[k]+.5*w[k], y[k]+.5*h[k], "%d" %(k+1))  
pylab.axis([-1.0, 26, -1.0, 26])  
pylab.xticks([])  
pylab.yticks([])  
 
pylab.show()

PIC