CVXPY constraints formulation

Question

I am trying to solve the Isoperimetric problem (7.14) from Additional Exercises for Convex Optimization by Stephen Boyd using CVXPY. The problem formulation is:

The code of constraints is given below:

constraints = [ y[1] == 0,
                y[N] == 0,
                y[F] == yfixed[F],
                cp.abs(( y[i+2] - 2 * y[i+1] + y[i]) / h**2) <= C for i in np.arange(1,199) ] #not using the first constraint here

The constraints have for loops, and when I tried to formulate the problem according to the CVXPY documentation, I got the following error

Invalid syntax

cp.abs(( y[i+2] - 2 * y[i+1] + y[i]) / h**2) <= C for i in np.arange(1,199) ]
                                                  ^

How to use the loops in CVXPY constraints?

swag2198 · Accepted Answer

You need to express the constraints in terms of matrix-vector equalities and inequalities which follow the DCP protocol for cvxpy. To elaborate, I can see three kinds of constraints in this problem: (I am going to assume 0-based indexing for programming convenience for the rest of the answer.)

Take y as the N+1 dimensional optimization variable.

Fixed point equality constraints: These basically set some indices of the y vector to a set of given values. Note that the zero boundary conditions y[0] == 0 and y[N] == 0 also fall under this.
Perimeter constraint: This is to be computed using successive differences. And finally we set something like the sum of the square roots of 1 plus the squares of the differences to be less than L. This part is probably the most involved one to write following the cvxpy protocols.
Curvature constraints: This also involves a calculation similar to successive differences like above but this is much easier to write as you'll see this is simply a matrix-vector multiplication type constraint like the first one.

Now let's write the constraints in code.

Necessary imports:

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
from scipy.linalg import circulant

1. Equality constraints:

These basically pick some indices from y and set those to given values. This can be implemented as follows:

def equality_constraints(N, F, vals):
    '''
    Sets some indices (F) in the y vector to given values. Also sets 
    y[0] == 0 and y[N] == 0.
    Returns E (a matrix) and e (a vector) such that E @ y == e.
    '''
    E = np.zeros((F.shape[0]+2, N+1)) # '+2' for selecting y[0] and y[N] also
    e = np.zeros(F.shape[0]+2)
    
    E[0, 0] = 1
    E[-1, -1] = 1
    if F.shape[0]:
        E[1:-1:1][np.arange(F.shape[0]), F] = 1
        e[1:-1:1] = vals
    return E, e

E is a binary matrix of shape (F.shape[0] + 2) x (N+1) and it has exactly one column set to 1 in each row, giving an index for the (N+1) dimensional vector y and e contains the value for that index of y.

2. Perimeter constraint:

For this, we need successive differences of the form y[i+1] - y[i] for i = 0, 1, 2, . . . , N-1. Note that we can similarly construct a vector having this N successive differences as its elements. We can perform the square root and other operations on this vector easily using vectorized computations. Here we are constructing an N x (N+1) matrix M, which when multiplied by y will give the N differences.

def length_matrix(N):
    '''
    Returns L with [-1, 1, 0, . . . , 0] as first row and sliding it
    to the right to get the following rows.
    '''
    val = np.array([-1, 1])
    offsets = np.array([0, 1])
    col0 = np.zeros(N+1)
    col0[offsets] = val
    
    M = circulant(col0).T[:-(len(val) - 1)]
    return M

The matrix M will be a circulant matrix. I simply transposed it and removed the last row to get the desired matrix. You can see this post to know how to create one such matrix. M looks like this:

 array([[-1.,  1.,  0., ...,  0.,  0.,  0.],
        [ 0., -1.,  1., ...,  0.,  0.,  0.],
        [ 0.,  0., -1., ...,  0.,  0.,  0.],
        ...,
        [ 0.,  0.,  0., ...,  1.,  0.,  0.],
        [ 0.,  0.,  0., ..., -1.,  1.,  0.],
        [ 0.,  0.,  0., ...,  0., -1.,  1.]])

3. Curvature constraints:

Exactly same matrix calculation like the last one. Just repeat and slide [1, -2, 1] along the rows!

def curvature_constraints(N, C, h):
    '''
    Returns D and C_vec to be used as D @ y <= C and D @ y >= -C 
    as curvature constraints.
    '''
    val = np.array([1, -2, 1])
    offsets = np.array([0, 1, 2])
    col0 = np.zeros(N+1)
    col0[offsets] = val
    
    D = circulant(col0).T[:-(len(val) - 1)]
    D /= h**2
    C_vec = np.ones(D.shape[0]) * C
    return D, C_vec

I am dividing by h**2 in the matrix itself.

Example:

I have taken this example from the site of the book itself. The data is also available here.

L = 1.5
a = 1
C = 15
N = 200
h = a/N

F = np.array([20,40,140,180]) # fixed points
vals = np.array([0.1, 0.15, 0.15, 0.2])

# Declare an array for plotting purposes
yfixed = np.zeros(N+1)
yfixed[F] = vals

x = np.linspace(0, a, N+1)

Problem Formulation and Solution:

I am leaving it for you to understand how I have assembled the matrices in formulating the constraints, especially the perimeter one. This is not difficult, but might require you some practice depending on how comfortable you are with vectorization. The DCP page is a very good place to start.

y = cp.Variable(N+1)

E, e = equality_constraints(N, F, vals)
M = length_matrix(N)
D, d = curvature_constraints(N, C, h)

constraints = [
    E @ y == e,
    h * cp.sum(cp.norm(cp.vstack([(M @ y)/h, np.ones(N)]), p = 2, axis = 0)) <= L,
    D @ y <= d,
    D @ y >= -d
]

objective_function = h * cp.sum(y)
objective = cp.Maximize(objective_function)

problem = cp.Problem(objective, constraints)
problem.solve()

plt.plot(0, 0, 'ko')
plt.plot(a, 0, 'ko')
for i in F:
    plt.plot(x[i], yfixed[i], 'bo')

plt.plot(x, y.value) # y.value gives the value of the cp Variable
plt.savefig('curve.png')

I got the answer as 0.1594237500556726 for the above example and the curve looks like this:

I have checked this solution with few other contrived test cases to verify correctness. However, there might be other more efficient solutions formulating this problem differently or there might even be some unexpected or embarrassing errors here! Feel free to let me know in case there is some error or you find anything difficult to understand in the answer.