您需要根据cvxpy 的DCP 协议以矩阵向量等式和不等式的形式表达约束。
详细地说，我可以在这个问题中看到三种约束：（我将假设基于 0 的索引以方便编程以供其余答案使用。）
以y为N+1维度优化变量。
-
定点等式约束：这些基本上将
y 向量的一些索引设置为一组给定值。请注意，y[0] == 0 和 y[N] == 0 的零边界条件也属于此范围。
-
周界约束：这将使用连续的差异来计算。最后，我们将 1 的 平方根 加上差的 平方 的 sum 设置为小于
L。这部分可能是按照cvxpy 协议编写的最复杂的部分。
-
曲率约束：这也涉及类似于上面的连续差分的计算，但这更容易编写，因为您会看到这只是一个矩阵向量乘法类型的约束，就像第一个一样。李>
现在让我们在代码中编写约束。
必要的进口：
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 将是一个circulant 矩阵。我只是将其转置并删除了最后一行以获得所需的矩阵。您可以查看此post 以了解如何创建一个这样的矩阵。 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，曲线如下所示：
我已经用其他几个人为的测试用例检查了这个解决方案以验证正确性。但是，可能还有其他更有效的解决方案以不同的方式来表述这个问题，或者这里甚至可能出现一些意想不到或令人尴尬的错误！如果有任何错误或您发现答案中有任何难以理解的地方，请随时告诉我。