【2024 Optimal Control 16-745】【Lecture 2】integrators.ipynb功能分析

代码功能分析

---

导入库和项目设置

import Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()

• 功能：激活当前文件夹为 Julia 项目环境，并安装当前项目中缺失的依赖包。

• import Pkg：

  • 导入 Julia 的包管理模块 Pkg，用于管理项目依赖。

• Pkg.activate(@__DIR__)：

  • 激活当前脚本所在目录作为 Julia 项目环境。@__DIR__ 是当前文件的目录路径。

• Pkg.instantiate()：

  • 根据项目的 Project.toml 和 Manifest.toml 文件，安装所有依赖包。如果某些包未安装，这个命令会自动安装它们。

using LinearAlgebra
using PyPlot
using ForwardDiff

功能：加载必要的库：

• using LinearAlgebra：

  • 导入线性代数模块，提供矩阵运算、特征值分解等功能。属于 Julia 标准库的一部分。

• using PyPlot：

  • 导入 PyPlot 模块，用于绘制图形。PyPlot 是 Julia 对 Python Matplotlib 的接口，依赖于 Python 环境。

• using ForwardDiff：

  • 导入 ForwardDiff 模块，提供自动微分功能，常用于数值优化和科学计算。

---

定义摆动动力学函数

function pendulum_dynamics(x)l = 1.0g = 9.81θ = x[1]        # 摆动角度 θθ̇ = x[2]       # 摆动角速度 θ̇θ̈ = -(g / l) * sin(θ)  # 根据摆动力学方程，计算角加速度 θ̈return [θ̇; θ̈]  # 返回角速度和角加速度
end

• 功能：模拟单摆的动力学：
  • 输入 x 是状态向量 [θ, θ̇]，即当前角度和角速度。
  • 计算结果是新的状态向量 [θ̇, θ̈]，即角速度和角加速度。

没有控制项时，动力学方程中的质量项m会相互抵消

通过拉格朗日方程推导，质量$m$在动能和势能中对称出现，最终通过偏导数和求导操作被完全抵消。这反映了单摆在无外力（无控制项）情况下的动力学是一个与质量$m$无关的系统，运动仅由摆长$l$、重力加速度 $g$ 和角度 $\theta$ 决定。

θ 和 θ̇ 是 Julia 中的 Unicode 变量名，可以通过快捷输入实现。在代码中，这些符号只是普通变量名，它们的含义通过上下文和赋值来表达。

---

Forward Euler

定义前向欧拉法求解器

function pendulum_forward_euler(fun, x0, Tf, h)t = Array(range(0, Tf, step=h))  # 定义时间步长数组 tx_hist = zeros(length(x0), length(t))  # 创建状态记录数组，初始化为 0x_hist[:,1] .= x0  # 设置初始状态为 x0for k = 1:(length(t)-1)x_hist[:,k+1] .= x_hist[:,k] + h * fun(x_hist[:,k])  # 使用前向欧拉法计算下一步状态endreturn x_hist, t  # 返回状态历史和时间步长
end

• 功能：
  • 计算单摆运动的数值解，使用前向欧拉法。
  • fun 是状态导数的计算函数（这里是 pendulum_dynamics）。
  • x0 是初始状态 [θ₀, θ̇₀]。
  • Tf 是模拟总时间。
  • h 是时间步长。
  • 返回：
    • x_hist：时间序列中摆的角度和角速度。
    • t：时间序列。

.= 用于广播赋值，使得多个值可以逐元素地赋给目标数组或集合。

前向欧拉法-角度随时间的变化

x0 = [0.1; 0]  # 初始角度为 0.1 弧度，初始角速度为 0
x_hist, t = pendulum_forward_euler(pendulum_dynamics, x0, 5, 0.1)  # 运行模拟，总时间 5 秒，步长 0.1 秒plot(t, x_hist[1,:])  # 绘制角度随时间变化曲线

• 功能：
  1. 设置单摆的初始状态
  2. 使用 pendulum_forward_euler 求解 5 秒内的摆动状态，步长为 0.1 秒。
  3. 绘制角度（x_hist[1,:]）随时间变化的曲线。

可视化图像

生成的图像展示了单摆在模拟时间内的角度随时间的变化，是一个振荡曲线（近似简谐振荡）。

由于采用的是前向欧拉法，随着时间步长的增加，数值解可能会出现一定的误差，尤其是在较长时间的模拟中。

1. 数值方法的稳定性：

   • 前向欧拉法是一种显式的数值积分方法，虽然实现简单，但在处理刚性方程或需要高精度的情况下可能不够稳定。可以考虑使用更高阶的方法，如四阶龙格-库塔法（RK4）。

2. 时间步长的选择：

   • 时间步长 h 影响数值解的精度和稳定性。较小的步长通常提供更高的精度，但会增加计算量。需要根据具体问题权衡选择合适的步长。

• 总模拟时间5s，时间步长0.1s
  -

• 总模拟时间5s，时间步长0.01s

• 

• 总模拟时间5s，时间步长0.01s

• 总模拟时间20s，时间步长0.01s

• 总模拟时间20s，时间步长0.001s

• 总模拟时间100s，时间步长0.001s

计算摆动系统离散化后的线性化矩阵 pendulum_euler_Ad

function pendulum_euler_Ad(x0, h)g = 9.81Ad = [1 h; -g*h*cos(x0[1]) 1]
end

功能解释：

1. function pendulum_euler_Ad(x0, h)

   • 定义了一个名为 pendulum_euler_Ad 的函数，输入参数为：
     • x0：当前状态向量（包括角度和角速度）。
     • h：时间步长。
   • 这个函数计算摆动系统离散化后的线性化矩阵。

2. g = 9.81

   • 定义重力加速度为 9.81 m/s²。

3. Ad = [1 h; -g*h*cos(x0[1]) 1]

   • 定义欧拉法计算的状态转移矩阵：
     • 第一行表示角度更新（线性项）。
     • 第二行表示角速度更新，包含重力和角度余弦的线性化影响。

调用函数并计算特征值

eigvals(pendulum_euler_Ad(0, 0.1))

功能解释：

1. pendulum_euler_Ad(0, 0.1)

   • 计算当初始角度 $x_0=0$，时间步长 $h=0.1$ 时的离散化矩阵。

2. eigvals()

   • 计算该矩阵的特征值，返回的结果为：

     1.0 - 0.31320919526731655im
     1.0 + 0.31320919526731655im
     

   • 表示系统状态转移矩阵的两个复数特征值。

等价于使用使用 ForwardDiff 自动求雅可比矩阵

using ForwardDiffh = 0.1 # 时间步长
Ad = I + h * ForwardDiff.jacobian(pendulum_dynamics, [0.0, 0.0]) # 计算A_d矩阵
eigvals(Ad) # 计算特征值

前向欧拉法-绘制特征值模与时间步长的关系图

eig_norm = zeros(100)
h = LinRange(0, 0.1, 100)
for k = 1:length(eig_norm)eig_norm[k] = max(norm.(eigvals(pendulum_euler_Ad([0;0], h[k]))))
end
plot(h, eig_norm)

功能解释：

1. eig_norm = zeros(100)

   • 初始化一个大小为 100 的数组，用于存储不同时间步长下特征值模的最大值。

2. h = LinRange(0, 0.1, 100)

   • 创建一个从 0 到 0.1 的线性空间，分为 100 等分，表示时间步长的变化范围。

3. for k = 1:length(eig_norm)

   • 开始循环，对每个时间步长计算对应的最大特征值模。

4. pendulum_euler_Ad([0;0], h[k])

   • 计算当前时间步长下的状态转移矩阵。

5. eigvals()

   • 计算当前矩阵的特征值。

6. norm.(eigvals(...))

   • 计算所有特征值的模。

7. max(...)

   • 获取当前时间步长下特征值模的最大值。

8. plot(h, eig_norm)

   • 绘制时间步长 $h$ 和对应最大特征值模的关系图。

图像分析：

• 图像显示特征值模随时间步长 $h$ 的变化：
  • 当 $h\to 0$ 时，特征值模接近 1，系统趋于边际稳定。
  • 随着 $h$ 增大，特征值模逐渐超过 1，系统变得不稳定。

这段代码实现了摆动系统的离散化线性化分析，探索了时间步长对系统稳定性的影响，并通过绘制关系图验证了系统的边界稳定性条件。

---

四阶龙格-库塔方法（RK4）

实现四阶龙格-库塔方法（RK4）fd_pendulum_rk4

function fd_pendulum_rk4(xk, h)f1 = pendulum_dynamics(xk)f2 = pendulum_dynamics(xk + 0.5*h*f1)f3 = pendulum_dynamics(xk + 0.5*h*f2)f4 = pendulum_dynamics(xk + h*f3)return xk + (h/6.0)*(f1 + 2*f2 + 2*f3 + f4)
end

功能：
实现四阶龙格-库塔方法（RK4）的核心步骤：

1. f1 到 f4 是四次不同点计算的斜率估计值。
2. 最终返回加权平均的更新状态，权重为 $1,2,2,1$。

函数 pendulum_rk4

function pendulum_rk4(fun, x0, Tf, h)t = Array(range(0, Tf, step=h))x_hist = zeros(length(x0), length(t))x_hist[:,1] .= x0for k = 1:(length(t)-1)x_hist[:,k+1] .= fd_pendulum_rk4(x_hist[:,k], h)endreturn x_hist, t
end

功能：
用 RK4 模拟摆动系统。

1. x_hist 用于存储状态历史，t 是时间步数组。
2. 初始状态为 x0，逐步调用 fd_pendulum_rk4 更新状态。

RK4法模拟单摆系统角度随时间变化

x0 = [1.0; 0]
x_hist2, t_hist2 = pendulum_rk4(pendulum_dynamics, x0, 10, 0.1)
plot(t_hist2, x_hist2[1,:])

功能：

1. 定义初始状态 x0 = [1.0; 0]（初始角度为 1 弧度，初始角速度为 0）。
2. 运行 RK4 模拟，总时长 10 秒，时间步长 0.1。
3. 绘制角度随时间变化的曲线。
   

使用 ForwardDiff 计算雅可比矩阵

using ForwardDiff
Ad = ForwardDiff.jacobian(x -> fd_pendulum_rk4(x, 0.1), [0; 0])
norm.(eigvals(Ad))

功能：

1. using ForwardDiff:

   • 引入 ForwardDiff 包，用于自动微分计算梯度、雅可比矩阵（Jacobian）和高阶导数。
   • 在这里，它被用来计算一个函数的雅可比矩阵。

2. ForwardDiff.jacobian(x -> fd_pendulum_rk4(x, 0.1), [0; 0]):

   • 雅可比矩阵的计算：

     • x -> fd_pendulum_rk4(x, 0.1) 是一个匿名函数，它将输入 $x$ 映射到通过四阶龙格-库塔方法（fd_pendulum_rk4）计算的状态更新结果。
     • fd_pendulum_rk4 是一个积分器函数，用于模拟系统的动力学行为。它的输入参数：
       • $x$ 是当前状态向量（如 [角度; 角速度]）。
       • 0.1 是时间步长 $h$。
     • 雅可比矩阵是关于状态变量 $x$ 的偏导数矩阵，其元素定义为：
       $$J_{ij}=\frac{\partial f_i}{\partial x_j}$$
       在这里，$f_i$ 表示第 $i$ 个方程的结果，$x_j$ 表示第 $j$ 个输入变量。
       

   • 状态点 [0; 0]:

     • 指定初始状态为 $x=[0;0]$，即摆的角度和角速度都为零。

   • 输出 Ad:

     • 计算得到的雅可比矩阵 $A_d$，表示在 $x=[0;0]$ 处的局部线性化模型。

3. eigvals(Ad):

   • 计算矩阵 $A_d$ 的特征值：
     • 特征值描述了状态变化的动态特性，决定系统的稳定性。
     • 如果特征值的模（绝对值）小于 1，则系统是稳定的。

4. norm.(eigvals(Ad)):

   • 计算矩阵特征值的模（欧几里得范数）。
   • norm. 是广播操作符（dot syntax），对特征值数组逐个计算模。
   • 输出是一个数组，其中每个值表示对应特征值的模。

总结

• 该代码通过自动微分计算四阶龙格-库塔方法在特定状态（$x=[0;0]$）下的雅可比矩阵。
• 通过特征值的模，分析离散时间系统在该状态附近的稳定性。
• 如果所有特征值的模小于 1，系统在该状态下是局部稳定的。

RK4法特征值模随时间步长的变化

eig_norm = zeros(100)
h = LinRange(0, 1, 100)
for k = 1:length(eig_norm)eig_norm[k] = max(norm.(eigvals(ForwardDiff.jacobian(x -> fd_pendulum_rk4(x, h[k]), [0; 0]))))
end
plot(h, eig_norm)

功能：

1. 初始化特征值模数组 eig_norm 和时间步长范围 h。
2. 循环计算不同时间步长下特征值模的最大值。
3. 绘制时间步长和特征值模的关系图，分析稳定性。
   

---

Backward Euler

反向欧拉法 pendulum_backward_euler

function pendulum_backward_euler(fun, x0, Tf, dt)t = Array(range(0, Tf, step=dt))x_hist = zeros(length(x0), length(t))x_hist[:,1] .= x0for k = 1:(length(t)-1)e = 1x_hist[:,k+1] .= x_hist[:,k]while e > 1e-8xn = x_hist[:,k] + dt.*fun(x_hist[:,k+1])e = norm(xn - x_hist[:,k+1])x_hist[:,k+1] .= xnendendreturn x_hist, t
end

该函数 pendulum_backward_euler 使用反向欧拉法Backward Euler Method模拟摆动系统的动力学行为。

工作原理总结

1. 反向欧拉法：

   • 隐式方法，在每个时间步需要迭代求解状态更新。
   • 更加稳定，适合处理刚性系统。

2. 迭代过程：

   • 初始猜测下一步状态值。
   • 利用当前猜测值计算下一步状态的动力学函数。
   • 反复更新，直到误差小于设定的阈值。

函数头

function pendulum_backward_euler(fun, x0, Tf, dt)

• fun: 表示系统的动力学函数，输入当前状态返回状态导数。例如，对于摆动系统，这个函数返回角速度和角加速度。
• x0: 初始状态向量，通常表示 [角度; 角速度]。
• Tf: 总模拟时间。
• dt: 时间步长。

时间数组与状态历史初始化

t = Array(range(0, Tf, step=dt))
x_hist = zeros(length(x0), length(t))
x_hist[:,1] .= x0

1. t:

   • 创建时间数组，从 0 到 Tf，以步长 dt 增加。
   • 每个时间点对应系统的状态更新。

2. x_hist:

   • 初始化一个矩阵，用于存储状态历史。行数为状态变量的数量（如角度和角速度），列数为时间点数量。

3. x_hist[:,1] .= x0:

   • 将第一列（即初始时间点）的状态设置为初始状态 x0。

时间步循环

for k = 1:(length(t)-1)e = 1x_hist[:,k+1] .= x_hist[:,k]

• 遍历所有时间步，k 是当前时间步索引。
• 初始化误差 e = 1（用于收敛判断）。
• 用当前状态 x_hist[:,k] 作为下一状态 x_hist[:,k+1] 的初始猜测值。

反向欧拉迭代

while e > 1e-8xn = x_hist[:,k] + dt.*fun(x_hist[:,k+1])e = norm(xn - x_hist[:,k+1])x_hist[:,k+1] .= xn
end

1. 迭代公式：

   • 反向欧拉法的核心方程：
     $$x_{k+1}=x_k+h\cdot f(x_{k+1})$$
   • 方程是隐式的，因为 $x_{k+1}$ 既在左侧，也在右侧的动力学函数中。

2. 计算新估计值 xn:

   • 使用当前估计值 x_hist[:,k+1] 计算右侧的动力学值 $f(x_{k+1})$。
   • 更新 $x_{k+1}$ 的新估计 xn。

3. 误差计算 e:

   • 计算当前估计值 xn 与之前的值 x_hist[:,k+1] 的欧几里得范数（即误差）。
   • 如果误差小于 $1\times 10^{-8}$，迭代停止，认为 $x_{k+1}$ 已经收敛。

4. 更新 x_hist[:,k+1]:

   • 将新的估计值 xn 赋值给下一步状态 x_hist[:,k+1]。

返回结果

return x_hist, t

• 返回：
  • x_hist: 包含整个模拟过程中每个时间点的状态历史。
  • t: 对应的时间点数组。

反向欧拉法-模拟角度随时间变化

x0 = [1.0; 0]
x_hist3, t_hist3 = pendulum_backward_euler(pendulum_dynamics, x0, 10, 0.001)
plot(t_hist3, x_hist3[1,:])

功能：

1. 定义初始状态 x0 = [1.0; 0]。
2. 用反向欧拉法模拟总时长 10 秒，时间步长为 0.001。
3. 绘制角度随时间变化的曲线，展示系统动力学行为。