CVXPY约束公式

Question

我试图用CVXPY从凸优化的附加练习中解决Stephen提出的等周问题(7.14)。问题的表述是：

约束准则如下：

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

约束有for循环，当我试图根据CVXPY文档描述问题时，我得到了以下错误

无效语法

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

如何在CVXPY约束中使用循环？

Answer 1

您需要用矩阵向量等式和不等式来表示约束，这些约束遵循磷酸氢钙协议，适用于cvxpy。为了详细说明，我可以在这个问题中看到三种约束：(我将假设基于0的索引，以便为其余的答案提供编程方便。)

将y作为N+1维优化变量。

1. 不动点等式约束：这些约束基本上将y向量的一些索引设置为一组给定的值。请注意，零边界条件y[0] == 0和y[N] == 0也属于这种情况。
2. 周长约束：这将使用逐次差异计算。最后，我们将1的平方根的和加上差分的平方设置为小于L。这部分可能是按照cvxpy协议编写的最复杂的部分。
3. 曲率约束：这也涉及类似于上面类似的连续差分的计算，但是编写起来要容易得多，因为您将看到，这只是一个矩阵向量乘法类型的约束，就像第一个约束一样。

现在让我们用代码编写约束。

必要进口：

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

1.平等制约因素：

它们基本上从y中选择一些索引，并将它们设置为给定的值。这可以按以下方式实施：

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是形状(F.shape[0] + 2) x (N+1)的二进制矩阵，它在每一行中恰好有一列设置为1，给出了(N+1)维向量y的索引，而e包含该索引y的值。

2.周界限制：

为此，我们需要对用于y[i+1] - y[i]的i = 0, 1, 2, . . . , N-1表单进行连续的差异。请注意，我们可以类似地构造一个向量，将此N逐次差异作为其元素。通过矢量化计算，我们可以很容易地对这个向量执行平方根和其他操作。这里我们构造了一个N x (N+1)矩阵M，当乘以y时，它会给出N的差异。

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

矩阵M将是循环矩阵。我只需将其转置，并移除最后一行以获得所需的矩阵。您可以看到这个帖子来了解如何创建这样一个矩阵。M看起来是这样的：

 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.曲率限制：

和上一次的矩阵计算完全一样。只需重复并沿行滑动[1, -2, 1]！

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

我除以矩阵本身中的h**2。

示例：

我从这本书本身的网站上拿了这个例子。数据也是可用的这里。

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)

问题的制定和解决：

我把它留给你们去理解，我是如何组合矩阵的，特别是周长的约束。这并不困难，但可能需要您进行一些练习，这取决于您对矢量化的适应程度。磷酸氢钙页面是一个非常好的起点。

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')

对于上面的例子，我以0.1594237500556726的形式得到了答案，曲线如下所示：

​

为了验证正确性，我用了很少的其他人为的测试用例来检查这个解决方案。然而，可能有其他更有效的解决方案，形成这个问题的不同，甚至可能有一些意想不到的或尴尬的错误！请让我知道，如果有什么错误，或您发现任何难以理解的答案。

Answer 2

试着分裂成两半：

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) ]