【2024 Optimal Control 16-745】【Lecture4】equality-constraints.ipynb功能分析

• 代码实现了一个二次优化问题的可视化解法，包括目标函数、约束以及优化路径。
• 提供了两种优化方法：牛顿法和高斯-牛顿法，用于对比效果。
• 利用了自动微分工具 ForwardDiff 来计算约束曲率。

---

环境初始化并导入依赖库

# 激活当前文件夹下的项目环境，并确保依赖安装完成
import Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()

# 导入必要的包
using LinearAlgebra # 提供线性代数工具，例如矩阵运算
using ForwardDiff  # 支持自动微分，用于梯度计算
using PyPlot       # 提供基于Python的Matplotlib绘图库，用于绘制图形

---

定义目标函数和梯度

# 定义二次型的权重矩阵为对角矩阵 Q
Q = Diagonal([0.5; 1])

# 定义目标函数 f(x)，是一个二次型
function f(x)
    return 0.5 * (x - [1; 0])' * Q * (x - [1; 0])
end

# 定义目标函数 f(x) 的梯度 ∇f(x)
function ∇f(x)
    # 梯度公式
    return Q * (x - [1; 0])
end

# 定义目标函数 f(x) 的 Hessian 矩阵 ∇²f(x)
function ∇2f(x)
    # 对于二次型函数，Hessian 是常数矩阵 Q
    return Q
end

• 功能：
  1. Q：一个对角矩阵，定义二次型函数的权重。
  2. f(x)：目标函数，二次型形式。
  3. ∇f(x)：目标函数的梯度，利用矩阵的乘法进行计算。
  4. ∇2f(x)：目标函数的 Hessian 矩阵，对二次型函数来说是常数矩阵 Q。
• 参数：
  • x 是二维向量，表示优化变量。

---

定义约束和其导数

# 定义约束函数 c(x)，这是一个非线性约束
function c(x)
    return x[1]^2 + 2 * x[1] - x[2]
end

# 定义约束函数 c(x) 的雅可比矩阵∂c(x)
function ∂c(x)
    return [2 * x[1] + 2, -1]
end

• 功能：
  1. c(x)：定义约束函数（非线性）。
  2. ∂c(x)：约束函数的雅可比矩阵。
• 参数：
  • x 是二维向量，表示优化变量。
    

---

绘制目标函数和约束的图形

# 定义绘制目标函数和约束的图形
function plot_landscape()
    Nsamp = 20  # 采样点数量
    # 创建网格上的 X 和 Y 点
    Xsamp = kron(ones(Nsamp), LinRange(-4, 4, Nsamp)')
    Ysamp = kron(ones(Nsamp)', LinRange(-4, 4, Nsamp))
    Zsamp = zeros(Nsamp, Nsamp)  # 用于存储 f(x) 的值

    # 遍历每个网格点，计算目标函数 f(x) 的值
    for j = 1:Nsamp
        for k = 1:Nsamp
            Zsamp[j, k] = f([Xsamp[j, k]; Ysamp[j, k]])
        end
    end

    # 绘制目标函数的等高线图
    contour(Xsamp, Ysamp, Zsamp)

    # 绘制约束函数 c(x) 的曲线
    xc = LinRange(-3.2, 1.2, Nsamp)
    plot(xc, xc.^2 + 2.0 .* xc, "y")  
end

# 调用绘图函数，展示目标函数和约束
plot_landscape()

• 功能：
  1. 生成网格数据 Xsamp 和 Ysamp。
  2. 计算每个点上的目标函数值 Zsamp。
  3. 绘制目标函数的等高线图。
  4. 绘制约束函数 c(x) 的曲线。
• 参数：
  • Nsamp 是采样点的数量。
  • LinRange 用于生成从 -4 到 4 的均匀分布数值。

绘制二维等高线图需要网格点

• 等高线图的本质：等高线图需要在二维平面中每个网格点计算函数值，然后通过等值线连接具有相同函数值的点。
• 为此，我们需要在$(x,y)$ 的定义域上定义一个规则的网格：
  1. 每个网格点对应平面上的一个$(x,y)$坐标。
  2. 每个网格点的 $z$ 值由目标函数 $f(x,y)$计算得出。

生成网格的方式

Julia 中的 kron（Kronecker 积）函数被用来构造重复的矩阵。

1. 构造 -4到 4 的序列

LinRange(-4, 4, Nsamp)

• 生成一个长度为 $Nsamp=20$ 的线性序列，表示横轴或纵轴上的采样点。

2. 生成重复行的矩阵

Xsamp = kron(ones(Nsamp), LinRange(-4, 4, Nsamp)')

• 目的：生成一个矩阵，其中每一行都是从 -4到 4 的序列，用于表示二维网格中 $x$ 方向的坐标。
• 分解解释：
  1. ones(Nsamp) 生成一个长度为 $Nsamp$ 的全 1 向量，用于 Kronecker 积操作。
  2. LinRange(-4, 4, Nsamp)' 转置生成的序列，变成 1 行 $Nsamp$列。
  3. kron(ones(Nsamp), LinRange(-4, 4, Nsamp)')：
     • 取 Kronecker 积，相当于把 LinRange(-4, 4, Nsamp)' 的每一行重复 $Nsamp$次，生成一个矩阵，每一行都是相同的 -4到 4 的序列。

3. 生成重复列的矩阵

Ysamp = kron(ones(Nsamp)', LinRange(-4, 4, Nsamp))

• 目的：生成一个矩阵，其中每一列都是从 -4 到 4 的序列，用于表示二维网格中 $y$ 方向的坐标。
• 分解解释：
  1. ones(Nsamp)' 转置生成的全 1 行向量。
  2. LinRange(-4, 4, Nsamp) 是纵轴的序列。
  3. kron(ones(Nsamp)', LinRange(-4, 4, Nsamp))：
     • Kronecker 积相当于把 LinRange(-4, 4, Nsamp) 的每一列重复 $Nsamp$次，生成一个矩阵，每一列都是相同的 -4 到 4 的序列。

结果：二维网格的坐标点

对于 $Nsamp=3$的网格生成，横轴和纵轴的坐标矩阵如下：

横轴 $x$ 坐标矩阵 $Xsamp$：

• LinRange(-4, 4, Nsamp)' 生成了一个行向量：[-4.0, 0.0, 4.0]。
• ones(Nsamp) 生成了一个列向量：[1.0; 1.0; 1.0]（3 行 1 列）。
• kron(ones(Nsamp), LinRange(-4, 4, Nsamp)')：
  • Kronecker 积的规则：将 ones(Nsamp) 的每个元素与行向量 LinRange(-4, 4, Nsamp)' 相乘，并将结果堆叠。
  • 这里，相乘后，每一行都复制了 [-4.0, 0.0, 4.0]，形成了一个 $3\times 3$的矩阵，每一行相同。表示 $x$ 坐标沿着纵轴重复。

纵轴 $y$坐标矩阵 $Ysamp$：

每一列的值相同，表示 $y$ 坐标沿着横轴重复。

网格点组合 $(x,y)$

将 $Xsamp[j,k]$ 和 $Ysamp[j,k]$ 组合起来，得到网格点,对应的二维平面排列为：

(-4.0, -4.0)   (0.0, -4.0)   (4.0, -4.0)
(-4.0,  0.0)   (0.0,  0.0)   (4.0,  0.0)
(-4.0,  4.0)   (0.0,  4.0)   (4.0,  4.0)

---

定义牛顿法步进

# 定义牛顿法的单步更新函数
function newton_step(x0, λ0)
    # 计算目标函数的 Hessian 和约束的二阶导数项
    H = ∇2f(x0) + ForwardDiff.jacobian(x -> ∂c(x)' * λ0, x0)
    C = ∂c(x0)  # 计算约束的导数
    # 形成 KKT 系统，并求解增量 Δz
    Δz = [H C'; C 0] \ [-∇f(x0) - C' * λ0; -c(x0)]
    Δx = Δz[1:2]  # 更新变量 x 的增量
    Δλ = Δz[3]    # 更新拉格朗日乘子 λ 的增量
    return x0 + Δx, λ0 + Δλ  # 返回更新后的 x 和 λ
end

• 功能：
  1. 计算目标函数的 Hessian 和约束的 Jacobian 矩阵。
  2. 形成 KKT 系统并求解增量 $\Delta z=[\Delta x;\Delta \lambda ]$。
  3. 更新变量 $x$ 和拉格朗日乘子 $\lambda$。
• 参数：
  • x0：当前变量的初值。
  • λ0：当前拉格朗日乘子的初值。

使用牛顿法求解

# 初始化优化变量和拉格朗日乘子的初始值
xguess = [-3; 2]  # 初始变量 x 的值
λguess = [0.0]    # 初始拉格朗日乘子 λ 的值

# 绘制初始点和优化路径
plot_landscape()
plot(xguess[1], xguess[2], "rx")  # 标记初始点

# 执行一次牛顿法迭代
xnew, λnew = newton_step(xguess[:, end], λguess[end])
# 更新变量的轨迹
xguess = [xguess xnew]
λguess = [λguess λnew]

# 绘制更新后的路径
plot_landscape()
plot(xguess[1, :], xguess[2, :], "rx")  # 标记新的优化路径

# 计算 Hessian 矩阵和约束的二阶导数，用于调试
H = ∇2f(xguess[:, end]) + ForwardDiff.jacobian(x -> ∂c(x)' * λguess[end], xguess[:, end])

• 功能：
  1. 初始化优化变量 xguess 和拉格朗日乘子 λguess。
  2. 进行一次牛顿法迭代。
  3. 绘制优化路径。
     

---

定义高斯-牛顿法步进

# 定义高斯-牛顿法的单步更新函数
function gauss_newton_step(x0, λ0)
    H = ∇2f(x0)  # 仅使用目标函数的 Hessian
    C = ∂c(x0)  # 计算约束的梯度
    # 形成简化的 KKT 系统，并求解增量 Δz
    Δz = [H C'; C 0] \ [-∇f(x0) - C' * λ0; -c(x0)]
    Δx = Δz[1:2]  # 更新变量 x 的增量
    Δλ = Δz[3]    # 更新拉格朗日乘子 λ 的增量
    return x0 + Δx, λ0 + Δλ  # 返回更新后的 x 和 λ
end

• 功能：
  1. 使用目标函数的 Hessian 矩阵，省略约束的曲率项。
  2. 形成简化的 KKT 系统，求解增量 $\Delta z$。
  3. 更新变量 $x$ 和拉格朗日乘子 $\lambda$。
• 参数：
  • x0：当前变量的初值。
  • λ0：当前拉格朗日乘子的初值。

使用高斯-牛顿法求解

# 重新初始化优化变量和拉格朗日乘子的初始值
xguess = [-3; 2]  # 初始变量 x 的值
λguess = [0.0]    # 初始拉格朗日乘子 λ 的值

# 绘制初始点和优化路径
plot_landscape()
plot(xguess[1], xguess[2], "rx")  # 标记初始点

# 执行一次高斯-牛顿法迭代
xnew, λnew = gauss_newton_step(xguess[:, end], λguess[end])
# 更新变量的轨迹
xguess = [xguess xnew]
λguess = [λguess λnew]

# 绘制更新后的路径
plot_landscape()
plot(xguess[1, :], xguess[2, :], "rx")  # 标记新的优化路径

• 功能：
  1. 初始化优化变量和拉格朗日乘子。
  2. 进行一次高斯-牛顿法迭代。
  3. 绘制优化路径。
• 参数：
  • xguess 和 λguess 是初始值，最终生成优化路径。

原文地址：https://blog.csdn.net/qq_45762996/article/details/144010163

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！