一文搞懂什么是 序列凸优化 SCP

概念介绍

序列凸优化(Sequential Convex Programming) SCP
一种用于解决非线性优化问题的迭代方法。在这一方法中，原始的非线性问题被近似为一系列更易处理的凸优化子问题，每个子问题都相对原始问题有一个凸近似。SCP的基本思想是通过迭代求解这些凸子问题来逼近原始非线性问题的最优解。
简单来说： 很多实际问题的目标函数和约束条件都是非线性的，而非线性函数往往带来非凸性，若将非线性函数在某个点进行近似，则可将问题转换为线性约束和线性目标函数的凸优化问题，从而更容易求解。
但是由于线性近似存在较大误差，所以只能在小范围内进行近似，所以需要迭代求解，在上一次迭代的求解结果处进行线性化，再进行本轮优化，如此往复迭代进行，即可获得最终优化结果。

SCP基本求解步骤

SCP的基本步骤
1.初始点选择： 选择一个初始点，从这个点开始迭代。
2.线性化或凸近似： 在当前迭代点处，将非线性目标函数和约束条件进行线性化或者采用其他凸近似方法，从而构建一个凸优化子问题。
3.求解凸子问题： 利用现有的凸优化算法（如内点法、梯度下降等）求解该凸子问题。
4.更新迭代点： 将子问题的解作为新的迭代点。
5.收敛性检查： 检查迭代是否达到预设的收敛标准（例如，目标函数值的变化小于某个阈值），如果满足则停止迭代，否则返回第2步。

SCP实践——非线性约束最小二乘凸优化求解

设定优化问题和约束条件：
这里目标函数是凸函数 无需进行变换
$$\begin{array}{rl}\text{minimize}& x_1^2+x_2^2\\ \text{subject to}& x_1^2-x_2\le 0\\ & x_1+x_2^2\ge 1\end{array}$$
可见目标函数是凸函数，无需进行变换。
图中蓝色区域和紫色区域的交集便是约束条件，显然，这是一个非凸约束，无法直接求解。
注意横坐标为x1 纵坐标为x2

首先假设初始位置为(0.5,0.5)
在(0.5,0.5)位置对约束条件进行线性化得到线性约束：
$$\begin{array}{rl}\text{minimize}& x_1^2+x_2^2\\ \text{subject to}& x_1-x_2\le -0.25\\ & x_1+x_2\ge 1.25\end{array}$$

图中蓝色区域和紫色区域为线性化之后的约束，曲线为原始约束边界。
对于$x_1=0.5$和$x_2=0.5$的情况和原始非线性约束大致一致。
所以在此时蓝色和紫色的交集中作为本轮迭代的约束条件即可进行优化求解。

CVXPY 问题求解

定义目标函数(目标问题)：

objective = cp.Minimize(cp.sum_squares(x))

定义初始估计点：

x_k = np.array([0.5, 0.5])  # Starting point

开始迭代：

for iteration in range(max_iters):
    linearized_constraint_1 = 2 * x_k[0] * (x[0] - x_k[0]) - x[1] <= -x_k[0]**2
    linearized_constraint_2 = x[0] + 2 * x_k[1] * (x[1] - x_k[1]) >= 1 - x_k[1]**2

    problem = cp.Problem(objective, [linearized_constraint_1, linearized_constraint_2])
    problem.solve()
    
    x_k = x.value
    
    current_obj_value = problem.value
    if abs(prev_obj_value - current_obj_value) < tolerance:
        print(f"Converged after {iteration + 1} iterations.")
        break

    prev_obj_value = current_obj_value
else:
    print("Reached maximum number of iterations without convergence.")

简要说明：
对于每轮迭代，首先计算当前估计值(来自上一次迭代结果)的位置的线性化函数，如果目标函数存在需要线性化的步骤，则在此处进行。
之后根据这个线性化约束构建Problem并进行求解。
求解之后更新最优点，以便下一次线性化在此处展开。
最后是对于终止迭代的条件判断，通过计算本次结果与上一次结果的向量的模来表示收敛程度，若差值基本为0，表示已经收敛，停止迭代。

求解结果

Converged after 5 iterations.
Optimal value of x: [0.50154162 0.70601918]
Optimal objective value: 0.7500070860957275

将求解结果进行绘制：

可见，该凸优化方法找到了最优解。

完整python程序

import cvxpy as cp
import numpy as np

# Define the problem parameters
max_iters = 20  # Maximum number of iterations
tolerance = 1e-4  # Convergence tolerance

# Initial guess for the variables
x_k = np.array([0.5, 0.5])  # Starting point

# Create variables
x = cp.Variable(2)

# Objective function
objective = cp.Minimize(cp.sum_squares(x))

# Store previous objective value
prev_obj_value = float('inf')

for iteration in range(max_iters):
    # Linearize the non-convex constraint x1^2 - x2 <= 0 around x_k
    linearized_constraint_1 = 2 * x_k[0] * (x[0] - x_k[0]) - x[1] <= -x_k[0]**2

    # Linearize the non-convex constraint x1 + x2^2 >= 1 around x_k
    linearized_constraint_2 = x[0] + 2 * x_k[1] * (x[1] - x_k[1]) >= 1 - x_k[1]**2

    # Define the problem with the linearized constraints
    problem = cp.Problem(objective, [linearized_constraint_1, linearized_constraint_2])

    # Solve the convex problem
    problem.solve()

    # Update the iterate
    x_k = x.value

    # Check for convergence
    current_obj_value = problem.value
    if abs(prev_obj_value - current_obj_value) < tolerance:
        print(f"Converged after {iteration + 1} iterations.")
        break

    prev_obj_value = current_obj_value
else:
    print("Reached maximum number of iterations without convergence.")

# Output the results
print("Optimal value of x:", x_k)
print("Optimal objective value:", current_obj_value)

关于线性化