【2024 Optimal Control 16-745】Julia语法

Lecture 2

θ和它的导数符号是通过 Julia 中的变量命名方式实现的

1. 变量 θ 的输入：

   • 在 Julia 中，θ 是一个合法的变量名，就像普通的字母 x 或 y 一样。
   • 要输入 θ，可以使用以下方法：
     • 在 Jupyter Notebook 或 Julia REPL 中直接键入 \theta，然后按 Tab 键，Julia 会自动将其转化为 θ。
     • 这是 Julia 支持的 Unicode 字符的一部分，可以直接作为变量名。

2. 变量 θ̇ 的输入：

   • θ̇ 是带点的变量（表示导数），它也是 Julia 的一个合法变量名。
   • 键入 \theta，后按 Tab 键,键入 \dot，后按 Tab 键

3. 变量 θ̈ 的输入：

   • θ̈ 是加速度（带两个点的变量名）。
   • 键入 \theta，后按 Tab 键,键入 \ddot，后按 Tab 键

• θ 和 θ̇ 是 Julia 中的 Unicode 变量名，可以通过快捷输入实现。
• 在代码中，这些符号只是普通变量名，它们的含义通过上下文和赋值来表达：
  • θ 是角度。
  • θ̇ 是角速度。
  • θ̈ 是角加速度，计算公式是 $-(g/l)\sin (\theta )$。
    

---

. 使得操作符能够逐个元素地应用于整个数组或集合，而不需要显式地使用循环

在 Julia 中，等号前面的 . 是一个 广播（broadcasting） 操作符，它的作用是让一个操作应用于整个数组或容器，而不仅仅是单一的元素。这是 Julia 中非常常见和强大的特性。

x_hist[:,1] .= x0

• x_hist[:,1]：这是对 x_hist 数组的访问，表示提取 x_hist 的第一列（所有行，第一列）。
• .=：这是广播赋值操作符。它表示对 x_hist[:,1] 的每个元素执行赋值操作，等价于逐元素地把 x0 的每个值赋给 x_hist[:,1] 的每个元素。
• x0：这是一个向量，通常是初始状态向量，比如 [0.1; 0]，表示摆的初始角度和角速度。

示例

假设 x0 是 [0.1, 0]，而 x_hist[:,1] 是一个长度为 2 的列向量（例如 [0, 0]），那么：

• 无广播（普通赋值）：
  如果直接使用 =（没有点操作符），会发生维度不匹配的错误，因为你不能将一个列向量直接赋值给数组的某一列。

  x_hist[:,1] = x0  # 这会报错
  

• 使用广播赋值：
  使用 .= 后，x0 的每个元素都会依次赋值给 x_hist[:,1] 的每个元素：

  x_hist[:,1] .= x0  # 正确：x0 的每个元素逐个赋值给 x_hist[:,1]
  

总结

• . 操作符使得赋值操作可以应用到整个数组或容器的每个元素，而不仅仅是简单的元素赋值。
• .= 用于广播赋值，使得多个值可以逐元素地赋给目标数组或集合。

这种广播机制使得在 Julia 中处理数组和矩阵时更加简洁高效。

---

解释匿名函数 x -> fd_pendulum_rk4(x, 0.1)

在 Julia 中，匿名函数是一种无需命名的函数，通常用于临时计算。形式为：

x -> expression

其中，x 是输入参数，expression 是处理该参数的表达式。

代码部分

x -> fd_pendulum_rk4(x, 0.1)

这是一个匿名函数，具体含义如下：

1. 输入参数：

   • x 是匿名函数的输入参数，代表当前状态变量，通常为一个向量。例如，对于摆动系统，x 可能是 $x=[\text{angle};\text{angular velocity}]$（角度和角速度）。

2. 调用 fd_pendulum_rk4:

   • 匿名函数内部调用了 fd_pendulum_rk4 函数，这是一个实现四阶龙格-库塔（RK4）积分方法的函数。
   • fd_pendulum_rk4(x, 0.1) 的两个参数：
     • x: 当前状态（例如 $x=[\text{角度};\text{角速度}]$）。
     • 0.1: 时间步长 $h$，表示模拟的离散时间间隔。

3. 输出结果：

   • 函数返回 RK4 方法计算的状态更新结果，即从状态 $x$ 经一步积分后的新状态 $x_{n+1}$。
   • 该结果是一个新的状态向量，表示在给定时间步长下，系统从状态 $x$ 演化到的下一个状态。

用途

匿名函数 x -> fd_pendulum_rk4(x, 0.1) 的核心作用是：

• 将输入状态 $x$（如摆的当前角度和角速度）映射为通过 RK4 方法计算得到的下一个状态。
• 在调用 ForwardDiff.jacobian 时，匿名函数为 ForwardDiff 提供了所需的输入-输出关系。

示例

假设：

• 初始状态为 $x=[\pi /4;0]$（角度 $\pi /4$，角速度 0）。
• 时间步长 $h=0.1$。

匿名函数的计算流程如下：

# 定义匿名函数
f = x -> fd_pendulum_rk4(x, 0.1)# 调用匿名函数
new_state = f([pi/4; 0])

结果 new_state 是通过 RK4 积分计算出的下一步状态（新的角度和角速度）。

为什么使用匿名函数？

1. 简洁性：

   • 不需要显式定义一个新函数，而是直接将 fd_pendulum_rk4 封装成满足特定需求的函数（固定步长为 0.1）。

2. 灵活性：

   • 匿名函数可以动态封装不同的参数和逻辑。例如，步长可以通过匿名函数灵活指定。

3. 与 ForwardDiff.jacobian 配合：

   • ForwardDiff.jacobian 需要输入一个函数，该函数的输入是状态 $x$，输出是对应的更新结果。匿名函数很好地满足这一要求。

总结

• 匿名函数 x -> fd_pendulum_rk4(x, 0.1) 将状态变量 $x$ 映射为通过 RK4 方法计算得到的下一个状态。
• 它的主要作用是为 ForwardDiff.jacobian 提供输入-输出映射关系，以计算状态更新过程的雅可比矩阵。

---

Lecture 3

∂r \ r 表示 矩阵左除

在 Julia 中，表达式 ∂r \ r 表示 矩阵左除，也就是 求解线性方程组 的一种简洁方式。

线性方程组的求解

线性方程组的一般形式是：

$$A\cdot x=b$$

其中：

• $A$ 是系数矩阵（这里对应 ∂r）。
• $x$ 是未知量向量（这里对应 $\Delta x$）。
• $b$ 是右侧的已知向量（这里对应 r）。

A \ x 的含义是 求解 $x$ 的值，即：

$$x=A^{-1}\cdot b$$

• 这相当于将 $A$ 的逆矩阵 $A^{-1}$ 左乘到 $b$ 上，求解 $x$ 的值。

但是，显式求逆（即 $A^{-1}$）的计算代价很高，且可能会引入数值不稳定性。因此，∂r \ r 使用了一种数值高效的方式解决这个问题。

Julia 的 \ 运算符

• 在 Julia 中，A \ b 是求解线性方程组 $A\cdot x=b$ 的符号，表示“将矩阵 $A$ 左除向量 $b$”。
• 实现时，Julia 使用优化的数值线性代数方法（如 LU 分解、QR 分解或 Cholesky 分解）来高效求解，而不是直接计算矩阵的逆。

数值计算的优势

• 高效性: 求解 $A\cdot x=b$ 的方法通常比显式逆矩阵的计算更高效。
• 数值稳定性: 显式计算逆矩阵可能会导致数值不稳定性（尤其当矩阵接近奇异时），而直接求解方程组能够减少误差。
• 灵活性: Julia 的 \ 运算符会自动选择最适合的分解算法（如 LU、QR 或其他方法）来解决问题，适用于稠密矩阵或稀疏矩阵。

---

梯度和雅可比矩阵在 Julia 中的使用规则

• 梯度 (gradient): 用于标量值函数，返回一个列向量。
• 雅可比矩阵 (jacobian): 用于矢量值函数，返回一个矩阵。
• 在实际使用中，必须明确函数的输入和输出维度，误用可能导致报错。

1. 梯度（gradient）

梯度是用于标量值函数的，它返回的是一个列向量。

例子

using ForwardDiff# 定义一个标量函数
f(x) = x[1]^2 + x[2]^2 + x[3]^2# 对 f 求梯度
x = [1.0, 2.0, 3.0]
grad = ForwardDiff.gradient(f, x)
println("梯度: ", grad)

输出

• 解析:
  • 函数 $f(x)=x_1^2+x_2^2+x_3^2$。
  • 梯度为 $\nabla f=[\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2},\frac{\partial f}{\partial x_3}{]}^T=[2x_1,2x_2,2x_3{]}^T$。
  • 结果为 [2.0, 4.0, 6.0]。

---

2. 雅可比矩阵（Jacobian）

雅可比矩阵用于矢量值函数，返回的是一个矩阵。

例子

using ForwardDiff# 定义一个矢量值函数
g(x) = [x[1]^2, x[1]*x[2], x[2]^2]# 对 g 求雅可比矩阵
x = [1.0, 2.0]
jacobian = ForwardDiff.jacobian(g, x)
println("雅可比矩阵: ")
println(jacobian)

输出

• 解析:

  • 函数 $g(x)=[x_1^2,x_1{x}_2,x_2^2]$。

  • 雅可比矩阵为：

    $J=\left[\begin{array}{cc}\frac{\partial g_1}{\partial x_1}& \frac{\partial g_1}{\partial x_2}\\ \frac{\partial g_2}{\partial x_1}& \frac{\partial g_2}{\partial x_2}\\ \frac{\partial g_3}{\partial x_1}& \frac{\partial g_3}{\partial x_2}\end{array}\right]$=$\left[\begin{array}{cc}2x_1& 0\\ x_2& x_1\\ 0& 2x_2\end{array}\right]$

  • 在 $x=[1.0,2.0]$ 时，结果为：
    $\left[\begin{array}{cc}2.0& 0.0\\ 2.0& 1.0\\ 0.0& 4.0\end{array}\right]$

3. 错误调用的情况

错误调用 gradient 对矢量值函数

如果尝试对矢量值函数调用 gradient 会导致报错，因为梯度只适用于标量值函数。

例子

g(x) = [x[1]^2, x[1]*x[2], x[2]^2]
x = [1.0, 2.0]
grad = ForwardDiff.gradient(g, x)  # 错误

报错信息

• 原因:
  • gradient 只能对标量值函数使用，而这里的 $g(x)$ 是矢量值函数。

错误调用 jacobian 对标量值函数

如果尝试对标量值函数调用 jacobian，理论上应返回梯度的转置，但通常会导致报错。

例子

f(x) = x[1]^2 + x[2]^2
x = [1.0, 2.0]
jacobian = ForwardDiff.jacobian(f, x)  # 错误

报错信息

• 原因:
  • jacobian 期望输入是矢量值函数，而这里的 $f(x)$ 是标量值函数。

---

Lecture 4

Kronecker 积（Kronecker Product）

Kronecker 积的定义

给定两个矩阵 $A\in {\mathbb{R}}^{m\times n}$ 和 $B\in {\mathbb{R}}^{p\times q}$，它们的 Kronecker 积 $A\otimes B$ 是一个大小为 $(mp)\times (nq)$ 的矩阵。具体构造规则如下：

• $A\otimes B$ 将矩阵 $A$ 的每个元素 $a_{ij}$ 替换为该元素与矩阵 $B$ 的乘积 $a_{ij}B$。

数学表达为：
$$A\otimes B=\left[\begin{array}{cccc}{a}_{11}B& a_{12}B& \dots & a_{1n}B\\ a_{21}B& a_{22}B& \dots & a_{2n}B\\ ⋮& ⋮& & ⋮\\ a_{m1}B& a_{m2}B& \dots & a_{mn}B\end{array}\right]$$

Kronecker 积的计算方法

假设：
$$A=\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right],B=\left[\begin{array}{cc}0& 5\\ 6& 7\end{array}\right]$$

计算 $A\otimes B$：
$$A\otimes B=\left[\begin{array}{cc}1\cdot B& 2\cdot B\\ 3\cdot B& 4\cdot B\end{array}\right]=\left[\begin{array}{cc}\left[\begin{array}{cc}0& 5\\ 6& 7\end{array}\right]& \left[\begin{array}{cc}0& 10\\ 12& 14\end{array}\right]\\ \left[\begin{array}{cc}0& 15\\ 18& 21\end{array}\right]& \left[\begin{array}{cc}0& 20\\ 24& 28\end{array}\right]\end{array}\right]$$

最终结果为：
$$A\otimes B=\left[\begin{array}{cccc}0& 5& 0& 10\\ 6& 7& 12& 14\\ 0& 15& 0& 20\\ 18& 21& 24& 28\end{array}\right]$$

Kronecker 积的性质

1. 尺寸：
   如果 $A$ 是 $m\times n$，$B$ 是 $p\times q$，那么 $A\otimes B$ 的大小为 $(mp)\times (nq)$。
2. 分布律：
   $$(A+C)\otimes B=A\otimes B+C\otimes B$$
3. 结合律：
   $$(A\otimes B)\otimes C=A\otimes (B\otimes C)$$
4. 与标量的关系：
   $$(\alpha A)\otimes B=A\otimes (\alpha B)$$

Kronecker 积的应用

1. 生成重复矩阵：

   • 在 Julia 中，kron(ones(m), A) 会生成一个矩阵，其中矩阵 $A$ 的每一行重复 $m$ 次。
   • 类似地，kron(ones(m)', A) 会生成一个矩阵，其中矩阵 $A$ 的每一列重复 $m$ 次。

2. 向量化操作：

   • Kronecker 积常用于将向量化表达与矩阵展开结合。比如，将矩阵 $A$ 的每一项与另一个矩阵 $B$ 关联。

3. 量子计算：

   • Kronecker 积在量子力学中用于描述复合量子系统的状态和操作，比如计算张量积态。

4. 系统理论和信号处理：

   • Kronecker 积用于构造大规模系统矩阵，特别是在多维信号处理中的应用。

在 Julia 中的实现

在 Julia 中，kron 函数用于计算 Kronecker 积，语法为：

C = kron(A, B)

在equality-constraints.ipynb代码中的作用

在 plot_landscape 函数中，Kronecker 积被用来生成网格点：

1. kron(ones(Nsamp), LinRange(-4, 4, Nsamp)')：
   • 生成一个矩阵，其中每一行是从 $-4$ 到 $4$ 的序列，表示 $x$ 坐标。
2. kron(ones(Nsamp)', LinRange(-4, 4, Nsamp))：
   • 生成一个矩阵，其中每一列是从 $-4$ 到 $4$ 的序列，表示 $y$ 坐标。

这样，通过 Kronecker 积可以快速构造二维网格，用于绘制等高线图。