CVXPY 约束公式
CVXPY constraints formulation
我正在尝试使用 CVXPY 解决来自 Stephen Boyd Additional Exercises for Convex Optimization 的 Isoperimetric 问题 (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 文档来表述问题时，出现以下错误
Invalid syntax
cp.abs(( y[i+2] - 2 * y[i+1] + y[i]) / h**2) <= C for i in np.arange(1,199) ]
^
如何在 CVXPY 约束中使用循环？
尝试一分为二：
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) ]
您需要根据 cvxpy 的 DCP 协议，根据矩阵向量等式和不等式来表达约束。
详细来说，我可以在这个问题中看到三种约束：（为了编程方便，我将假设基于 0 的索引用于其余答案。）
取y作为N+1维优化变量
- 不动点等式约束：这些基本上将
y 向量的一些索引设置为一组给定值。请注意，零边界条件 y[0] == 0 和 y[N] == 0 也属于此。
- 周界约束：这是用逐次差分来计算的。最后我们设置类似 sum of the square roots of 1 plus the squares of differences小于
L。这部分可能是遵循 cvxpy 协议编写的最复杂的部分。
- Curvature constraints：这也涉及类似于上述连续差异的计算，但是这更容易编写，因为您会看到这只是一个矩阵向量乘法类型像第一个一样的约束。
现在让我们在代码中编写约束。
必要的进口：
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。周边限制：
为此，对于 i = 0, 1, 2, . . . , N-1，我们需要 y[i+1] - y[i] 形式的连续差异。请注意，我们可以类似地构造一个具有此 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 将是一个 circulant matrix. I simply transposed it and removed the last row to get the desired matrix. You can see this 以了解如何创建这样的矩阵。 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。
示例：
我从本书本身的站点上获取了这个示例。数据也可用 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)
问题的表述和解决方案：
我留给您了解我是如何组装矩阵来制定约束的，尤其是周长约束。这并不困难，但可能需要您进行一些练习，具体取决于您对矢量化的熟悉程度。 DCP 页面是一个很好的起点。
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，曲线看起来是这样的：
我已经用其他一些人为的测试用例检查了这个解决方案以验证正确性。然而，可能有其他更有效的解决方案以不同方式表述这个问题，或者这里甚至可能存在一些意想不到或令人尴尬的错误！如果有错误或您发现答案中有任何难以理解的地方，请随时告诉我。
我正在尝试使用 CVXPY 解决来自 Stephen Boyd Additional Exercises for Convex Optimization 的 Isoperimetric 问题 (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 文档来表述问题时，出现以下错误
Invalid syntax
cp.abs(( y[i+2] - 2 * y[i+1] + y[i]) / h**2) <= C for i in np.arange(1,199) ]
^
如何在 CVXPY 约束中使用循环？
尝试一分为二：
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) ]
您需要根据 cvxpy 的 DCP 协议，根据矩阵向量等式和不等式来表达约束。
详细来说，我可以在这个问题中看到三种约束：（为了编程方便，我将假设基于 0 的索引用于其余答案。）
取y作为N+1维优化变量
- 不动点等式约束：这些基本上将
y向量的一些索引设置为一组给定值。请注意，零边界条件y[0] == 0和y[N] == 0也属于此。 - 周界约束：这是用逐次差分来计算的。最后我们设置类似 sum of the square roots of 1 plus the squares of differences小于
L。这部分可能是遵循cvxpy协议编写的最复杂的部分。 - Curvature constraints：这也涉及类似于上述连续差异的计算，但是这更容易编写，因为您会看到这只是一个矩阵向量乘法类型像第一个一样的约束。
现在让我们在代码中编写约束。
必要的进口：
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。周边限制：
为此，对于 i = 0, 1, 2, . . . , N-1，我们需要 y[i+1] - y[i] 形式的连续差异。请注意，我们可以类似地构造一个具有此 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 将是一个 circulant matrix. I simply transposed it and removed the last row to get the desired matrix. You can see this 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。
示例：
我从本书本身的站点上获取了这个示例。数据也可用 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)
问题的表述和解决方案：
我留给您了解我是如何组装矩阵来制定约束的，尤其是周长约束。这并不困难，但可能需要您进行一些练习，具体取决于您对矢量化的熟悉程度。 DCP 页面是一个很好的起点。
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，曲线看起来是这样的：
我已经用其他一些人为的测试用例检查了这个解决方案以验证正确性。然而，可能有其他更有效的解决方案以不同方式表述这个问题，或者这里甚至可能存在一些意想不到或令人尴尬的错误！如果有错误或您发现答案中有任何难以理解的地方，请随时告诉我。