486 lines
1.0 MiB
Plaintext
486 lines
1.0 MiB
Plaintext
|
{
|
|||
|
|
"cells": [
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "c667fb2b",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"# SAC论文复现学习\n",
|
|||
|
|
"找了一篇Frank L. Lewis经典的Actor–Critic-Based Optimal Tracking for Partially Unknown Nonlinear Discrete-Time Systems的论文来复现学习，主要是为了理解SAC算法在连续动作空间中的应用，该论文考虑了带约束情况下的最优跟踪问题，使用了Actor-Critic方法来解决。"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 43,
|
|||
|
|
"id": "ba88673c",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"import torch\n",
|
|||
|
|
"import torch.nn as nn\n",
|
|||
|
|
"import numpy as np\n",
|
|||
|
|
"import math\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 设置随机种子，保证每次运行的结果可复现\n",
|
|||
|
|
"torch.manual_seed(42)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 超参数与系统参数设置\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"gamma = 0.3 # 折扣因子\n",
|
|||
|
|
"u_bar = 0.4 # 控制饱和上限 (物理极限)\n",
|
|||
|
|
"R = 1.0 # 控制输入惩罚权重\n",
|
|||
|
|
"lr_c = 0.05 # Critic 网络的学习率\n",
|
|||
|
|
"lr_a = 0.05 # Actor 网络的学习率\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 状态惩罚矩阵 Q1\n",
|
|||
|
|
"# 论文中 Q=20I，由于增广状态 X = [e, r]，且只有误差 e 需要惩罚，所以后两项为0\n",
|
|||
|
|
"Q1 = torch.diag(torch.tensor([20.0, 20.0, 0.0, 0.0], dtype=torch.float32))"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "1125d88b",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 系统动态方程与参考轨迹发生器\n",
|
|||
|
|
"对应的是论文里 Section V 的仿真算例。我们要把离散非线性系统状态方程（公式 70）和正弦命令发生器（公式 71）实现出来 。另外，由于我们后面在 Actor 更新时需要用到输入增益矩阵 $G(X)$，所以这里单独写一个函数提取它。"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 44,
|
|||
|
|
"id": "807b6976",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 系统动态与参考轨迹定义\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"def system_step(x, u):\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" 非线性离散时间系统漂移与输入动态\n",
|
|||
|
|
" x: 当前状态向量 [x1, x2]\n",
|
|||
|
|
" u: 控制输入标量\n",
|
|||
|
|
" 返回: 下一时刻的状态 [x1_next, x2_next]\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" x1_next = -0.8 * x[1]\n",
|
|||
|
|
" x2_next = -0.45 * x[0] - torch.sin(x[1]) + 0.2 * x[1] * u\n",
|
|||
|
|
" return torch.stack([x1_next, x2_next])\n",
|
|||
|
|
"\n",
|
|||
|
|
"def ref_step(r):\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" 命令发生器动态 (产生振荡参考轨迹)\n",
|
|||
|
|
" r: 当前参考状态向量 [r1, r2]\n",
|
|||
|
|
" 返回: 下一时刻的参考状态\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" r1_next = -r[0]\n",
|
|||
|
|
" r2_next = -r[1]\n",
|
|||
|
|
" return torch.stack([r1_next, r2_next])\n",
|
|||
|
|
"\n",
|
|||
|
|
"def get_G_matrix(x):\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" 计算输入增益矩阵 G(X)\n",
|
|||
|
|
" 对应于被控系统中控制量 u 前面的状态系数矩阵\n",
|
|||
|
|
" 用于 Actor 目标控制量的解析计算\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" return torch.tensor([[0.0], [0.2 * x[1]], [0.0], [0.0]], dtype=torch.float32)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "effcf549",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 构建 Actor 和 Critic 神经网络\n",
|
|||
|
|
"论文中采用的是单隐藏层神经网络结构，隐藏层有 5 个神经元，激活函数为 $\\tanh(.)$ 。需要特别注意的是，Actor 网络的输出层我们直接加上了 $\\bar{u} \\tanh(.)$ 的限幅，以此来严格保证输出控制量永远不会超出你设定的物理极限 。"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 45,
|
|||
|
|
"id": "ebeed72d",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 神经网络定义与初始化\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"class CriticNet(nn.Module):\n",
|
|||
|
|
" def __init__(self):\n",
|
|||
|
|
" super(CriticNet, self).__init__()\n",
|
|||
|
|
" # Critic网络: 4个输入 -> 5个隐藏层节点(tanh) -> 1个标量输出(线性)\n",
|
|||
|
|
" self.fc1 = nn.Linear(4, 5)\n",
|
|||
|
|
" self.fc2 = nn.Linear(5, 1)\n",
|
|||
|
|
" # 将权重初始化在 -1 到 1 之间\n",
|
|||
|
|
" nn.init.uniform_(self.fc1.weight, -1.0, 1.0)\n",
|
|||
|
|
" nn.init.uniform_(self.fc2.weight, -1.0, 1.0)\n",
|
|||
|
|
" \n",
|
|||
|
|
" def forward(self, X):\n",
|
|||
|
|
" h = torch.tanh(self.fc1(X))\n",
|
|||
|
|
" V = self.fc2(h)\n",
|
|||
|
|
" return V\n",
|
|||
|
|
"\n",
|
|||
|
|
"class ActorNet(nn.Module):\n",
|
|||
|
|
" def __init__(self):\n",
|
|||
|
|
" super(ActorNet, self).__init__()\n",
|
|||
|
|
" # Actor网络: 4个输入 -> 5个隐藏层节点(tanh) -> 1个控制量输出\n",
|
|||
|
|
" self.fc1 = nn.Linear(4, 5)\n",
|
|||
|
|
" self.fc2 = nn.Linear(5, 1)\n",
|
|||
|
|
" nn.init.uniform_(self.fc1.weight, -1.0, 1.0)\n",
|
|||
|
|
" nn.init.uniform_(self.fc2.weight, -1.0, 1.0)\n",
|
|||
|
|
" \n",
|
|||
|
|
" def forward(self, X):\n",
|
|||
|
|
" h = torch.tanh(self.fc1(X))\n",
|
|||
|
|
" # 核心：输出层包含物理饱和约束 u_bar * tanh(.)\n",
|
|||
|
|
" u = u_bar * torch.tanh(self.fc2(h))\n",
|
|||
|
|
" return u\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 实例化网络与对应的随机梯度下降优化器\n",
|
|||
|
|
"critic = CriticNet()\n",
|
|||
|
|
"actor = ActorNet()\n",
|
|||
|
|
"optimizer_c = torch.optim.SGD(critic.parameters(), lr=lr_c)\n",
|
|||
|
|
"optimizer_a = torch.optim.SGD(actor.parameters(), lr=lr_a)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "f177c2b5",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 定义非二次型性能代价函数\n",
|
|||
|
|
"这部分是实现带输入约束控制的核心。通过引入非二次型积分代价 $W(u)$，使得控制算法在优化的同时兼顾输入饱和 。我们用 PyTorch 提供的基础数学算子来实现公式 (20) 的积分解析解。"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 46,
|
|||
|
|
"id": "19039547",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 辅助函数：非二次型积分代价计算\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"def non_quadratic_cost(u):\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" 计算包含输入约束的非二次型积分代价 W(u)\n",
|
|||
|
|
" 防止控制量在优化过程中无限制增长\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" # 稍微截断控制量，防止在反三角函数中出现数值溢出 (除以0或虚数)\n",
|
|||
|
|
" u_clamped = torch.clamp(u, -u_bar + 1e-5, u_bar - 1e-5)\n",
|
|||
|
|
" ratio = u_clamped / u_bar\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 积分的解析解公式\n",
|
|||
|
|
" integral = (ratio * torch.atanh(ratio)) + 0.5 * torch.log(1 - ratio**2)\n",
|
|||
|
|
" return 2 * (u_bar**2) * R * integral"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "840adc63",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 在线 Actor-Critic 同步学习主循环\n",
|
|||
|
|
"这是整个算法的引擎。我们在这里严格贯彻了刚才讨论的时间平移逻辑，利用滞后一步的状态数据来进行贝尔曼误差的计算和权重的反向传播。为了保证系统被充分探索，我们在输出动作上叠加了一个微小的探测噪声。"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 47,
|
|||
|
|
"id": "12444dfc",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"当前步数 10: Critic 损失 = 664.4067, Actor 损失 = 0.0018\n",
|
|||
|
|
"当前步数 20: Critic 损失 = 51.7882, Actor 损失 = 0.0008\n",
|
|||
|
|
"当前步数 30: Critic 损失 = 1.2308, Actor 损失 = 0.0004\n",
|
|||
|
|
"当前步数 40: Critic 损失 = 0.0579, Actor 损失 = 0.0001\n",
|
|||
|
|
"当前步数 50: Critic 损失 = 0.7083, Actor 损失 = 0.0001\n",
|
|||
|
|
"当前步数 60: Critic 损失 = 0.7751, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 70: Critic 损失 = 0.4341, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 80: Critic 损失 = 0.0691, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 90: Critic 损失 = 0.0727, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 100: Critic 损失 = 0.4509, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 110: Critic 损失 = 0.7882, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 120: Critic 损失 = 0.6399, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 130: Critic 损失 = 0.2181, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 140: Critic 损失 = 0.0008, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 150: Critic 损失 = 0.2530, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 160: Critic 损失 = 0.6776, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 170: Critic 损失 = 0.7690, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 180: Critic 损失 = 0.4190, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 190: Critic 损失 = 0.0433, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 200: Critic 损失 = 0.0860, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 210: Critic 损失 = 0.4970, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 220: Critic 损失 = 0.7905, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 230: Critic 损失 = 0.6150, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 240: Critic 损失 = 0.1838, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 250: Critic 损失 = 0.0039, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 260: Critic 损失 = 0.2894, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 270: Critic 损失 = 0.7027, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 280: Critic 损失 = 0.7510, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 290: Critic 损失 = 0.3837, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 300: Critic 损失 = 0.0299, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 310: Critic 损失 = 0.1111, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 320: Critic 损失 = 0.5306, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 330: Critic 损失 = 0.7897, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 340: Critic 损失 = 0.5876, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 350: Critic 损失 = 0.1564, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 360: Critic 损失 = 0.0110, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 370: Critic 损失 = 0.3220, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 380: Critic 损失 = 0.7205, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 390: Critic 损失 = 0.7392, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 400: Critic 损失 = 0.3486, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 410: Critic 损失 = 0.0171, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 420: Critic 损失 = 0.1345, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 430: Critic 损失 = 0.5621, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 440: Critic 损失 = 0.7964, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 450: Critic 损失 = 0.5534, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 460: Critic 损失 = 0.1282, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 470: Critic 损失 = 0.0198, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 480: Critic 损失 = 0.3579, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 490: Critic 损失 = 0.7433, Actor 损失 = 0.0000\n",
|
|||
|
|
"当前步数 500: Critic 损失 = 0.7149, Actor 损失 = 0.0000\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 核心在线学习主循环\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"total_steps = 500\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 初始化历史状态序列 (给一个初始偏差，观察跟踪效果)\n",
|
|||
|
|
"x_k_minus_1 = torch.tensor([1.0, -1.0], dtype=torch.float32)\n",
|
|||
|
|
"r_k_minus_1 = torch.tensor([0.5, -0.5], dtype=torch.float32)\n",
|
|||
|
|
"\n",
|
|||
|
|
"for k in range(1, total_steps + 1):\n",
|
|||
|
|
" # --- 第1部分：前向交互与数据采集 ---\n",
|
|||
|
|
" # 计算 k-1 时刻的误差和增广状态 X(k-1)\n",
|
|||
|
|
" e_k_minus_1 = x_k_minus_1 - r_k_minus_1\n",
|
|||
|
|
" X_k_minus_1 = torch.cat([e_k_minus_1, r_k_minus_1])\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 网络生成策略，并加入探测噪声以满足持续激励(PE)条件\n",
|
|||
|
|
" u_k_minus_1 = actor(X_k_minus_1)\n",
|
|||
|
|
" noise = torch.tensor([0.05 * math.sin(k)], dtype=torch.float32)\n",
|
|||
|
|
" u_applied = torch.clamp(u_k_minus_1 + noise, -u_bar, u_bar) \n",
|
|||
|
|
" \n",
|
|||
|
|
" # 物理系统演化：得到当前真正的 k 时刻状态\n",
|
|||
|
|
" x_k = system_step(x_k_minus_1, u_applied[0])\n",
|
|||
|
|
" r_k = ref_step(r_k_minus_1)\n",
|
|||
|
|
" e_k = x_k - r_k\n",
|
|||
|
|
" X_k = torch.cat([e_k, r_k])\n",
|
|||
|
|
" \n",
|
|||
|
|
" # --- 第2部分：Critic 网络时移更新 ---\n",
|
|||
|
|
" optimizer_c.zero_grad()\n",
|
|||
|
|
" \n",
|
|||
|
|
" V_k_minus_1 = critic(X_k_minus_1)\n",
|
|||
|
|
" # 截断 k 时刻价值的梯度，把它当作目标常数\n",
|
|||
|
|
" V_k = critic(X_k).detach() \n",
|
|||
|
|
" \n",
|
|||
|
|
" # 计算历史这一步的真实单步代价\n",
|
|||
|
|
" state_cost = torch.matmul(torch.matmul(X_k_minus_1.T, Q1), X_k_minus_1)\n",
|
|||
|
|
" ctrl_cost = non_quadratic_cost(u_applied[0])\n",
|
|||
|
|
" U_k_minus_1 = state_cost + ctrl_cost\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 必须将 td_target 彻底切断联系，不让它把梯度传导给 Actor\n",
|
|||
|
|
" td_target = (U_k_minus_1 + gamma * V_k).detach()\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 计算 TD 误差并反向传播更新 Critic\n",
|
|||
|
|
" td_target = U_k_minus_1 + gamma * V_k\n",
|
|||
|
|
" loss_c = 0.5 * (td_target - V_k_minus_1)**2\n",
|
|||
|
|
" loss_c.backward()\n",
|
|||
|
|
" optimizer_c.step()\n",
|
|||
|
|
" \n",
|
|||
|
|
" # --- 第3部分：Actor 网络时移更新 ---\n",
|
|||
|
|
" optimizer_a.zero_grad()\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 重新前向传播以构建计算图\n",
|
|||
|
|
" u_pred = actor(X_k_minus_1)\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 切断 X_k 与过去系统的联系，仅把它当做一个自变量来求偏导\n",
|
|||
|
|
" X_k_for_grad = X_k.clone().detach().requires_grad_(True)\n",
|
|||
|
|
" V_k_grad = critic(X_k_for_grad)\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 使用 torch.autograd.grad 纯净地求导，避免对 Critic 权重二次污染\n",
|
|||
|
|
" dV_dX = torch.autograd.grad(outputs=V_k_grad, inputs=X_k_for_grad)[0]\n",
|
|||
|
|
" \n",
|
|||
|
|
" \n",
|
|||
|
|
" # 计算目标控制量 (结合输入增益矩阵 G)\n",
|
|||
|
|
" G_matrix = get_G_matrix(x_k_minus_1)\n",
|
|||
|
|
" inner_term = - (gamma / 2.0) * (1.0 / (u_bar * R)) * torch.matmul(G_matrix.T, dV_dX)\n",
|
|||
|
|
" u_target = u_bar * torch.tanh(inner_term).detach()\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 最小化 Actor 输出与解析目标控制量的误差\n",
|
|||
|
|
" loss_a = 0.5 * (u_target - u_pred)**2\n",
|
|||
|
|
" loss_a.backward()\n",
|
|||
|
|
" optimizer_a.step()\n",
|
|||
|
|
" \n",
|
|||
|
|
" # --- 第4部分：时间轴步进 ---\n",
|
|||
|
|
" x_k_minus_1 = x_k.detach()\n",
|
|||
|
|
" r_k_minus_1 = r_k.detach()\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 每 10 步打印一次学习进度\n",
|
|||
|
|
" if k % 10 == 0:\n",
|
|||
|
|
" print(f\"当前步数 {k}: Critic 损失 = {loss_c.item():.4f}, Actor 损失 = {loss_a.item():.4f}\")"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 48,
|
|||
|
|
"id": "f8fec043",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAAHfCAYAAAAY1Jq/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsXXd8VMXafs6WJJveOymE3kORHrqACCKiSJFiQezdDy4qgugFUaoFBSmCAoLUgNKRLr23dNJ7NnX7fH8cZzLn7G4AFbnqee7P3w275+zMnDPzzlufEQghBAoUKFCgQIECBQoUKFCg4G8H1b3ugAIFChQoUKBAgQIFChQo+H1QDDoFChQoUKBAgQIFChQo+JtCMegUKFCgQIECBQoUKFCg4G8KxaBToECBAgUKFChQoECBgr8pFINOgQIFChQoUKBAgQIFCv6mUAw6BQoUKFCgQIECBQoUKPibQjHoFChQoECBAgUKFChQoOBvCsWgU6BAgQIFChQoUKBAgYK/KRSDToECBQoUKFCgQIECBQr+plAMOgUKFCj4iyEIAnr27Hlb165YsQKCIGDFihV3tU9/R6xZswbx8fHw8vKCIAh49dVX73WXFChwiDtZ8woUKFBwp1AMOgUKFPwjIAjCHf2nGEh/DDExMZLnqVKp4Ovriy5duuDzzz+HxWK5q+0fO3YMo0ePRkVFBZ577jlMmzYNAwYMuKttKvjrUFZWhjlz5mD06NFo1qwZNBoNBEHAnj177nXX/idw9epVTJs2DQ899BCioqLYOrzb606BAgX/m9Dc6w4oUKBAwZ+BadOm2X02f/586PV6vPLKK/D19ZV816ZNm7+mY38QDz/8MDp16oSwsLB73RWHoM/WarUiLS0NP/74I44dO4a9e/di48aNd63d7du3gxCCb7/9Fl26dLlr7Si4N0hPT8fbb78NAIiMjERgYCDy8/Pvca/+d7Bz507MmDEDarUaDRs2hJubGwwGw73ulgIFCu4RFINOgQIF/wi8//77dp+tWLECer0er776KmJiYv7yPv0Z8PHxgY+Pz73uhlPIn+2UKVPQoUMHbNq0Cb/88gt69OhxV9rNyckBAISHh9+V31dwbxEdHY09e/YgPj4e/v7+GD9+PFauXHmvu/U/g4EDB6Jz585o1aoVdDodYmJikJGRca+7pUCBgnsEJeVSgQIF/zr07NkTgiDAZDJhxowZaNy4MVxdXTF+/HgAgF6vx5w5c9C7d29ERkbCxcUFQUFBGDJkCI4dO+b0d69du4Ynn3wSMTExcHV1RXBwMLp3744vv/zytvo1Z84cqFQqdO3aFSUlJQCc19DFxMQgJiYGVVVVeOuttxAVFQVXV1c0aNAAs2fPBiHE7vcJIViwYAGaNWsGNzc3RERE4MUXX4Rer2e/90fRvHlzVit04sQJ9nlJSQmmTJmCpk2bQqfTwcfHB3369MGuXbvsfoMf888//4yePXvCx8eHfSYIApYvXw4AiI2NZelm6enp7DdOnz6NRx55BMHBwXB1dUV0dDSef/555Obm2rU3fvx4CIKA1NRULFq0iCnJdBz0+7S0NHz22Wfs+cXExOCjjz5iz3r9+vW477774OHhgeDgYLz44ouoqamxa2/z5s0YM2YMGjVqBA8PD3h4eKBdu3ZYuHAhbDab0/6lp6fjq6++QsuWLeHm5oaQkBBMnDgRer3e4bvIysrCyy+/jIYNG0Kn08Hf3x/33XcfPvjgA4fXvvjii6hfvz5cXV0REBCAIUOG4OTJkw5/+05QWlrK1sTp06cl39lsNvTq1QuCIGDVqlXscz8/P/Tp0wf+/v5/uP3fs55pzVtRUREmTpyIsLAwuLq6onnz5mzuyWEymfDBBx8gLi4Orq6uiI2NxTvvvAOj0XhH/U1PT4cgCBg/fjxu3LiBESNGIDg4GCqVCgcOHAAANG7cGB07doROp7uj31agQME/E0qEToECBf9aPPLIIzh58iQGDhyIoUOHIjg4GIBYnzJ16lQkJCRg0KBB8PPzw82bN7F161b89NNP2LZtm1291vbt2/Hoo4/CaDRiwIABGDlyJMrKynD+/Hl8/PHHeO6555z2w2az4dVXX8WiRYswbNgwfPfdd3Bzc7tl/81mM/r374+cnBwMHDgQGo0GmzdvxuTJk2EwGOzSUF944QV8+eWXCA8Px8SJE+Hi4oKtW7fixIkTMJvN0Gq1v+Mp2oMaOIIgAAAyMjLQs2dPpKeno3v37hgwYACqqqqQmJiIAQMG4KuvvsIzzzxj9zsbNmzAzz//jIEDB2LSpEnIyMhAmzZtMG3aNGzevBnnz5+XpNPS/09MTMQjjzwCQgiGDx+O6OhonD59Gl9++SW2bNmCw4cPIzY21q69V155BYcOHcKgQYPwwAMPQK1WS75/8803ceDAAQwePBj3338/tm7diqlTp8JkMsHf3x+TJ0/G0KFD0b17d+zevRuff/45rFarnUE/efJkqFQqdOzYEREREdDr9di3bx9eeeUVnDx5UmLY8Hj77bexc+dO1v7+/fuxZMkSJCcnY9++fZJrT506hf79+6OkpAQJCQkYNmwYqqurceXKFbz//vt499132bVnzpzB/fffj5KSEvTv3x/Dhg1DUVERNm/ejG7dumHTpk144IEHJL9P360jx4Ecfn5+WLNmDRISEjBixAicPXsWXl5eAIDp06fjwIEDGD9+PJ544olb/tbvwe9Zz4BYx9e1a1e4uLhg+PDhMBqNWL9+PZ588kmoVCqMGzeOXUsIwWOPPYYtW7YgLi4OL774IkwmE5YtW4aLFy/+rn6npKSgY8eOaNSoEUaPHo2amhp4e3v/7uegQIGCfzCIAgUKFPxDER0dTQCQtLQ0yec9evQgAEjLli1JYWGh3X1lZWUOP8/MzCRhYWGkSZMmks8LCwuJt7c30Wq15MCBAw7v4wGA9OjRgxBCSE1NDRk2bBgBQF588UVitVol1y5fvpwAIMuXL3c4toEDB5Lq6mr2eX5+PvHx8SE+Pj7EZDKxzw8ePEgAkEaNGpHS0lL2udFoJN27dycASHR0tF3fncHZs7106RLR6XQEADl48CAhRHzegiCQNWvWSK4tLS0lrVu3Jm5ubiQvL89uzIIgkJ9++slh++PGjXPYfkVFBfH39ycqlYq1TzFr1iwCgPTr18/hb4WHh5PU1FSnbUVHR5OsrCxJ/wMCAoi7uzsJDAwkV65cYd8ZDAbStGlT4uLiQvLz8yW/l5ycbNeG1WolY8eOJQDI8ePHHbZfr149kpGRwT43m83s3f3666/sc6PRSGJiYggA8t1339m1xc9Hs9lM4uLiiKurq93czc7OJuHh4SQ0NJQYDAbJdwDInaoQs2fPJgDI448/TgghZN++fUSlUpGmTZuSqqqqOu+lz2D37t131CYhd76eCakd31NPPUUsFgv7/PLly0StVpOmTZtKrv/uu+8IANKpUydSU1PDPi8uLib169eXrPlbIS0tjbU/ZcqU27qHrkez2Xxb1ytQoOCfBcWgU6BAwT8WtzLoNm/efMe/+dJLLxEAEsX6k08+IQDIyy+/fFu/QZW74uJi0rVrVyIIApk9e7bDa29l0CUlJdndQw2Dixcvss+eeuopAoCsXLnS7vrDhw//boPulVdeIdOmTSPvvPMOGT16NDPmHn74YUIIIefOnSMAyPDhwx3+zubNmwkA8vnnn9uNeejQoU7bd2bQrV69mgAgI0eOtLvHbDYzQ4d/f/S35s+fX2dbS5cutftuwoQJBAB599137b57//33CQCHRr4jnD59mgAg06dPd9j+kiVL7O5ZtmwZAUAWLVrEPtuwYQMBQIYMGXLLNunzf/PNNx1+P3/+fAKAbN++XfL51atXydWrV29nWAw2m40MGDCAACAfffQRCQsLIzqdjly4cOGW9/4Rg64uOFrPhIhr1N3dnej1ert7EhISCABSUVHBPuvbty8BQPbt22d3PZ3Pd2rQhYSE2BnSzqAYdAoU/LuhpFwqUKDgX4v77rvP6XdHjhzBggULcOzYMRQUFMBkMkm+z87ORlRUFADg+PHjAESigttFfn4
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1000x500 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAAHfCAYAAAAY1Jq/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXWcVtXWx7/nqeluZoYZultCOkVE7MJAMTCvYlxfBb1iXr3ea+fVKzYqNoqKgiggUlLSndMdzzx19vvHfmKfZ2YURLDO7/OBOc8+sdbavfZea21NCCEwYcKECRMmTJgwYcKECRN/OFh+awZMmDBhwoQJEyZMmDBhwsQvg6nQmTBhwoQJEyZMmDBhwsQfFKZCZ8KECRMmTJgwYcKECRN/UJgKnQkTJkyYMGHChAkTJkz8QWEqdCZMmDBhwoQJEyZMmDDxB4Wp0JkwYcKECRMmTJgwYcLEHxSmQmfChAkTJkyYMGHChAkTf1CYCp0JEyZMmDBhwoQJEyZM/EFhKnQmTJgwYcKECRMmTJgw8QeFqdCZMGHCxDGGpmkMHz78kJ59+eWX0TSNl19++ajy9EfErFmz6NWrF3FxcWiaxtSpU39rlkyYaBKH0+ZNmDBh4nBhKnQmTJj4U0DTtMP6ZypIR4b8/HxDflosFhITExk4cCBPP/00Xq/3qNJfunQpF1xwATU1NVx99dXcddddnHjiiUeVpoljhzVr1jBjxgwGDRpEVlYWDoeD7OxsJk6cyA8//PBbs/ebQgjB559/zt/+9jd69uxJUlISkZGRdOjQgalTp1JUVPRbs2jChIljDNtvzYAJEyZM/Bq46667GqU99thjVFVVccMNN5CYmGi417Nnz2PD2BHi9NNPZ8CAAWRlZf3WrDSJQN76fD527drFe++9x9KlS5k/fz7vv//+UaP76aefIoTg1VdfZeDAgUeNjonfBldddRXLli2jT58+nHHGGcTGxrJmzRreeust3n33Xd5++23OOOOM35rN3wQul4tx48bhcDgYOnQoo0ePxufzsWDBAh5//HHeeustFi1aRLt27X5rVk2YMHGMoAkhxG/NhAkTJkwcDeTn57Nnzx527dpFfn7+b81OEJqmMWzYMBYuXPhbs/KL0Vzebtiwgb59++J0Olm4cCHDhg07KvQvvfRSZs6c+bsrWxO/Dp588knGjRtH27ZtDelvvPEGF154ISkpKRw8eBCHw/EbcXh4+DXbvMfj4V//+hfXXHMNSUlJwXRd17nmmmt4/vnnOfnkk5kzZ84R0zJhwsQfA6bJpQkTJv5yGD58OJqm4Xa7ueeee+jQoQMRERFccsklAFRVVfHwww8zcuRIcnJycDgcpKWlccopp7B06dJmv7t582YuvfRS8vPziYiIID09nSFDhvDss88eEl8PP/wwFouFQYMGUV5eDjTvQ5efn09+fj51dXX8/e9/p2XLlkRERNC2bVseeughmlqrE0Lw+OOP07lzZyIjI8nOzua6666jqqoq+L0jRZcuXYK+QsuXLw+ml5eXc/vtt9OpUyeioqJISEhg1KhRzJs3r9E3VJk///xzhg8fTkJCQjBN0zRmzpwJQKtWrYJmn7t37w5+Y9WqVZx55pmkp6cTERFBXl4e11xzDQUFBY3oXXLJJWiaxs6dO3nyySfp3r07UVFRQTkC93ft2sVTTz0VzL/8/HweeOCBYF7Pnj2bfv36ERMTQ3p6Otdddx1Op7MRvQ8//JALL7yQ9u3bExMTQ0xMDH369OGJJ55A1/Vm+du9ezfPP/883bp1IzIykoyMDKZMmUJVVVWTZbF//36uv/562rVrR1RUFMnJyfTr14977723yWevu+46WrduTUREBCkpKZxyyimsWLGiyW8fDioqKoJtYtWqVYZ7uq4zYsQINE3jtddeC6b/7W9/a6TMAVxwwQW0a9eOsrIy1q9ff0j0f0l7Dvi8lZaWMmXKFLKysoiIiKBLly7BuhcOt9vNvffeS5s2bYiIiKBVq1bccccduFyuQ+IzgN27d6NpGpdccglbt27l3HPPJT09HYvFwsKFC7Hb7UyfPt2gzAFYLBb+8Y9/APyhF4tMmDBx+DBNLk2YMPGXxZlnnsmKFSsYN24cp512Gunp6QBs2rSJ6dOnM3ToUMaPH09SUhJ79+7l448/5rPPPmPOnDmN/LU+/fRTzj77bFwuFyeeeCITJ06ksrKStWvX8q9//Yurr766WT50XWfq1Kk8+eSTnHHGGbzxxhtERkb+LP8ej4exY8dy8OBBxo0bh81m48MPP+S2226joaGhkRnqtddey7PPPkuLFi2YMmUKDoeDjz/+mOXLl+PxeLDb7b8gFxsjoOBomgbAnj17GD58OLt372bIkCGceOKJ1NXV8cknn3DiiSfy/PPPc8UVVzT6zrvvvsvnn3/OuHHjuOqqq9izZw89e/bkrrvu4sMPP2Tt2rUGc9rA308++YQzzzwTIQRnnXUWeXl5rFq1imeffZaPPvqIxYsX06pVq0b0brjhBhYtWsT48eM56aSTsFqthvu33HILCxcuZMKECZxwwgl8/PHHTJ8+HbfbTXJyMrfddhunnXYaQ4YM4csvv+Tpp5/G5/M1Uuhvu+02LBYL/fv3Jzs7m6qqKhYsWMANN9zAihUrDIqNiltvvZUvvvgiSP/rr7/mhRdeYPv27SxYsMDw7MqVKxk7dizl5eUMHTqUM844g/r6ejZu3MiMGTO48847g8/+8MMPnHDCCZSXlzN27FjOOOMMSktL+fDDDxk8eDAffPABJ510kuH7gbI9FCOfpKQkZs2axdChQzn33HNZvXo1cXFxANx9990sXLiQSy65hIsuuuhnvwUE66nNdmhTmF/SngEqKysZNGgQDoeDs846C5fLxezZs7n00kuxWCxcfPHFwWeFEJxzzjl89NFHtGnThuuuuw63281LL710yIpnOHbs2EH//v1p3749F1xwAU6nk/j4+J9853DzxoQJE38SCBMmTJj4kyIvL08AYteuXYb0YcOGCUB069ZNlJSUNHqvsrKyyfR9+/aJrKws0bFjR0N6SUmJiI+PF3a7XSxcuLDJ91QAYtiwYUIIIZxOpzjjjDMEIK677jrh8/kMz86cOVMAYubMmU3KNm7cOFFfXx9MLyoqEgkJCSIhIUG43e5g+rfffisA0b59e1FRURFMd7lcYsiQIQIQeXl5jXhvDs3l7Y8//iiioqIEIL799lshhMxvTdPErFmzDM9WVFSIHj16iMjISFFYWNhIZk3TxGeffdYk/YsvvrhJ+jU1NSI5OVlYLJYg/QAefPBBAYgxY8Y0+a0WLVqInTt3NksrLy9P7N+/38B/SkqKiI6OFqmpqWLjxo3Bew0NDaJTp07C4XCIoqIiw/e2b9/eiIbP5xOTJk0SgPj++++bpJ+bmyv27NkTTPd4PMGyW7ZsWTDd5XKJ/Px8AYg33nijES21Pno8HtGmTRsRERHRqO4eOHBAtGjRQmRmZoqGhgbDPUAc7hTioYceEoA477zzhBBCLFiwQFgsFtGpUydRV1d3SN9YunSpAER2drbwer2H9M7htmchQvJddtllBjobNmwQVqtVdOrUyfD8G2+8IQAxYMAA4XQ6g+llZWWidevWhjb/c9i1a1eQ/u23335I7wQQqOOBPDZhwsRfA6ZCZ8KEiT8tfk6h+/DDDw/7m3/7298EYJhY//vf/xaAuP766w/pG4HJXVlZmRg0aJDQNE089NBDTT77cwrdtm3bGr0TUAzWr18fTLvssssEIF555ZVGzy9evPgXK3Q33HCDuOuuu8Qdd9whLrjggqAyd/rppwshhFizZo0AxFlnndXkdz788EMBiKeffrqRzKeddlqz9JtT6F5//XUBiIkTJzZ6x+PxBBUdtfwC33rsscd+ktaLL77Y6N7kyZMFIO68885G92bMmCGAJpX8prBq1SoBiLvvvrtJ+i+88EKjd1566SUBiCeffDKY9u677wpAnHLKKT9LM5D/t9xyS5P3H3vsMQGITz/91JC+adMmsWnTpkMRKwhd18WJJ54oAPHAAw+IrKwsERUVJdatW3dI75e
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1000x500 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA2QAAAHfCAYAAAA/ehtXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAlU9JREFUeJztnXd4FNX6x7/bsumBNBIghA7Si/SOIiBXRRQrSrkKXhHEdtVrQ8WCFa4oyEUREctP7IoFqRaUjtJbgBASIAnpbcv5/bHsZrZmN7s7u9n9fp4nj2dnzpzPmT0v4757Zs4qhBAChBBCCCGEEEJkRxnoDhBCCCGEEEJIuMKEjBBCCCGEEEICBBMyQgghhBBCCAkQTMgIIYQQQgghJEAwISOEEEIIIYSQAMGEjBBCCCGEEEICBBMyQgghhBBCCAkQTMgIIYQQQgghJEAwISOEEEIIIYSQAMGEjBBCGhgKhQLDhw93q+57770HhUKB9957z699IsSXzJ07FwqFAhs3bgx0VwghxO8wISOEEJiSHE/+mOB4R8uWLfke+5i1a9figQcewGWXXYakpCQoFAoMHjw40N0KSnJycvDGG29g7NixaNmyJbRaLZKSkjBq1Ch8/vnnge4eISTMUAe6A4QQEgw89dRTdtsWLFiA4uJi3HvvvWjUqJHVvh49esjTMS+59tpr0b9/f6Snpwe6Kw5x9N6aaSjvcbDw5ptv4quvvkJkZCTatm2LwsLCQHcpaHnjjTcwf/58tGrVCiNGjEBaWhpOnjyJzz//HD///DPuu+8+vPbaa4HuJiEkTGBCRgghMN0iZct7772H4uJizJkzBy1btpS9T74gISEBCQkJge6GUxryextsPPzww3juuefQsWNHZGdno1WrVoHuUtDSt29fbNy4EcOGDbPafuDAAfTv3x+vv/46br31VvTu3TtAPSSEhBO8ZZEQQjxk+PDhUCgUqKmpwTPPPIMOHTpAq9ViypQpAIDi4mK8/PLLGDlyJJo3b46IiAikpKTg6quvxpYtW5y2e/DgQUybNs1yC1VqaiqGDBmCxYsXu9Wvl19+GUqlEoMGDbLMjjh7hqxly5Zo2bIlysvL8dBDD6FFixbQarVo27Yt5s+fDyGEXftCCCxcuBCdOnVCZGQkmjVrhnvuuQfFxcWW9vyF9Dx++OEHDB8+HAkJCVAoFG7tB0zj8uijj6JDhw6IjIxE48aNMXr0aPz88892vo0bN0KhUGDu3LnYunUrxo0bh8TERCgUCpw4cQIA8Ndff+Hmm2+2jFdKSgp69eqFOXPmQKfT+eS8P/roI4wYMQKNGjVCZGQkLrnkEsybNw/V1dV2dQcMGIDOnTtDpVJ57d2wYQOmT5+OTp06IT4+HlFRUejSpQuefvppVFVV2dWXPvO1evVq9O3bF9HR0UhMTMRNN92EnJwch54dO3ZgzJgxiIuLQ3x8PC6//HKX/0ZcUVFRgRdeeAE9evRATEwMYmNjMWDAAHz00Ud2dSdMmGCXjAHAJZdcghtvvBEA+PwaIUQ2OENGCCH15LrrrsO2bdswduxYjB8/HqmpqQBM37I/9thjGDp0KMaNG4fGjRvj1KlT+Prrr/H999/jm2++wZgxY6za+u677zBx4kRUV1djzJgxuPnmm1FUVIQ9e/bgpZdewr/+9S+n/TAajZgzZw7eeOMNTJgwAatWrUJkZGSd/dfpdBg9ejTOnDmDsWPHQq1W48svv8QjjzyCqqoqu9s4Z86cicWLF6Np06aYPn06IiIi8PXXX2Pr1q3Q6XTQaDT1eBc9Y/Xq1fjhhx8wduxY3HXXXTh58qRb+4uKijBo0CDs378fffr0wZw5c5Cfn4//+7//wxVXXIHFixdjxowZdr4tW7bghRdewODBgzFt2jTk5+cjIiICf/31F/r16weFQoGrr74arVq1QklJCY4ePYq33noL8+bNs3o/zImho0TXGdOmTcPy5cvRvHlzXHfddWjUqBH++OMPPPHEE1i3bh3Wrl0Ltdo//xufP38+Dh48iIEDB2LcuHGoqqrCb7/9hrlz52Ljxo34+eefHSZ+b731Fr7++mtcffXVGDZsGP7880988skn2LNnD3bv3g2tVmup+/vvv+Pyyy9HTU0NJkyYgLZt22L37t0YPnw4Ro4c6VF/i4qKMHLkSOzatQu9evXCtGnTYDQa8eOPP+KWW27Bvn37MG/ePLfaMo+bv95bQgixQxBCCHFIZmamACCysrKstg8bNkwAEF27dhXnz5+3O66oqMjh9uzsbJGeni46duxotf38+fMiPj5eaDQasXHjRofHSQEghg0bJoQQorKyUkyYMEEAEPfcc48wGAxWdZcvXy4AiOXLlzs8t7Fjx4qKigrL9rNnz4qEhASRkJAgampqLNs3b94sAIj27duLCxcuWLZXV1eLIUOGCAAiMzPTru/OMPvvvfde8dRTTzn8q6ystDsPhUIhvv/+e7v26to/ffp0AUBMnz5dGI1Gy/bDhw+L+Ph4ERERYTXOGzZsEAAEALFkyRK79u6//34BQHz55Zd2+woLC+3GwdyWu5jP59prr7UaHyGEeOqppwQAsWDBAqfHZ2VlCQBi0KBBbjulHDt2zOp9MvP4448LAOLjjz922Ke4uDjx119/We27+eabBQDxySefWLYZjUbRoUMHh+/hggULLO/Xhg0b3Orv5MmTBQAxf/58q+2VlZVi9OjRQqFQiF27dtXZTnFxsWjSpIlQKBRi//79brkJIcRbmJARQogT6krIHH0Yr4tZs2YJAOLkyZOWba+88ooAIGbPnu1WG+aErKCgQAwaNEgoFAq7D6Jm6krIjhw5YnfM7bffLgCIv//+27Ltn//8pwAgVqxYYVf/119/rXdC5upPmviZz2P8+PEuz9PR/urqahEdHS1iY2NFQUGB3X5zkvH0009btpkTsh49ejj0mROyH3/80a3zPXDggDhw4IBbdYUQokePHkKtVlu9B2b0er1ISkoSffr0cXq8twmZMwoKCgQAMXXqVKvt5oTssccesztm/fr1AoB44IEHLNvMMTN06FC7+nq9XrRp08bthCw/P1+oVCpx6aWXOty/e/duAUA89NBDLtsxGo1i4sSJAoC4++676/QSQoiv4Hw8IYTUk759+zrd99tvv2HhwoXYsmULzp07h5qaGqv9OTk5aNGiBQDgjz/+AACMHTvWbffZs2cxaNAgHD9+HB988AFuueUWj/ufkJCAtm3b2m3PyMgAAFy4cMGybdeuXQDgcBn1/v371/v2rqysLI+ePXP1njvbf+jQIVRUVGDQoEFITEy02z9y5EjMmzfPco7u+G688UYsXLgQ48ePx/XXX4/LL78cgwYNQps2bRzW79ixo8t+S6moqMCePXuQnJyMBQsWOKyj1Wpx4MABt9v0lPLycixcuBBffPEFDh8+jNLSUqvbLZ09E3bppZfabXMUTzt37gQAh89xqVQqDB48GMeOHXOrr9u2bYPBYLA882eL+Xm+ut6vBx54AJ9++imGDBnCFRYJIbLChIwQQupJWlqaw+1ffPEFrr/+ekRGRmLUqFFo06YNYmJioFQqsXHjRmzatMlqUYaioiIAQLNmzdx25+XloaSkBM2bN6/3b005W27enFwZDAbLtuLiYgBAkyZN7OqrVCokJSXVqw+e4uw9d7Xf3HdnS/+bt5vHwR1f37598csvv+C5557D6tWrsXLlSgBAhw4d8NRTT+Hmm2922U9XXLhwAUIInD9/Hk8//XS926kvOp0OI0eOxNatW9GlSxfceOONSElJsTxb9fTTTztcVARwHFOexhNQ9zhLKSgoAGBKzLZt2+a0XllZmdN9//73v/H6669j6NCh+O6776yedSOEEH/DhIwQQuqJdAU/KU888QQiIiKwfft2XHLJJVb7ZsyYgU2bNlltM3+IzcnJQdeuXd1yd+/eHXfccQemTJmCoUOHYv369WjdurXnJ+Em8fHxAEwzc7Yeg8GAgoICjxLK+uLsPXe137zsf15ensNjcnNzreq56xswYAC+/fZbVFdXY8eOHfjhhx/wxhtv4JZbbkFKSgouv/x
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1000x500 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"import matplotlib.pyplot as plt\n",
|
|||
|
|
"\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 算法验证与结果可视化\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 切换为评估模式 (在此阶段不需要计算梯度)\n",
|
|||
|
|
"actor.eval()\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 验证阶段的步数\n",
|
|||
|
|
"eval_steps = 500\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 重新初始化状态，用于纯净的闭环测试\n",
|
|||
|
|
"x_eval = torch.tensor([1.5, -1.5], dtype=torch.float32)\n",
|
|||
|
|
"r_eval = torch.tensor([1.4, -1.4], dtype=torch.float32)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 用于存放画图数据的列表\n",
|
|||
|
|
"x1_history, x2_history = [], []\n",
|
|||
|
|
"r1_history, r2_history = [], []\n",
|
|||
|
|
"e1_history, e2_history = [], []\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 纯净闭环控制循环 (无梯度更新，无探测噪声)\n",
|
|||
|
|
"with torch.no_grad():\n",
|
|||
|
|
" for k in range(eval_steps):\n",
|
|||
|
|
" # 记录当前时刻的数据\n",
|
|||
|
|
" x1_history.append(x_eval[0].item())\n",
|
|||
|
|
" x2_history.append(x_eval[1].item())\n",
|
|||
|
|
" r1_history.append(r_eval[0].item())\n",
|
|||
|
|
" r2_history.append(r_eval[1].item())\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 记录跟踪误差 e = x - r\n",
|
|||
|
|
" e_eval = x_eval - r_eval\n",
|
|||
|
|
" e1_history.append(e_eval[0].item())\n",
|
|||
|
|
" e2_history.append(e_eval[1].item())\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 拼接增广状态\n",
|
|||
|
|
" X_eval = torch.cat([e_eval, r_eval])\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 使用训练好的 Actor 网络计算控制律 (注意：这里去掉了 noise)\n",
|
|||
|
|
" u_eval = actor(X_eval)\n",
|
|||
|
|
" # 加上物理限幅保护\n",
|
|||
|
|
" u_applied = torch.clamp(u_eval, -u_bar, u_bar)\n",
|
|||
|
|
"\n",
|
|||
|
|
" # 系统状态演化\n",
|
|||
|
|
" x_eval = system_step(x_eval, u_applied[0])\n",
|
|||
|
|
" r_eval = ref_step(r_eval)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"# 使用 Matplotlib 绘制跟踪效果图\n",
|
|||
|
|
"# ==========================================\n",
|
|||
|
|
"time_steps = range(eval_steps)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 设置全局字体大小，方便放入大论文\n",
|
|||
|
|
"plt.rcParams.update({'font.size': 12})\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 绘制状态 x1 与参考轨迹 r1 的对比图 (对应论文 Fig. 6)\n",
|
|||
|
|
"plt.figure(figsize=(10, 5))\n",
|
|||
|
|
"plt.plot(time_steps, x1_history, label='x1 (System State)', color='red', linewidth=1.5)\n",
|
|||
|
|
"plt.plot(time_steps, r1_history, label='r1 (Reference)', color='green', linestyle='--', linewidth=1.5)\n",
|
|||
|
|
"plt.title('Tracking Performance: x1 and r1')\n",
|
|||
|
|
"plt.xlabel('Time step')\n",
|
|||
|
|
"plt.ylabel('x1 and r1')\n",
|
|||
|
|
"plt.legend(loc='upper right')\n",
|
|||
|
|
"plt.grid(True)\n",
|
|||
|
|
"plt.xlim([0, eval_steps])\n",
|
|||
|
|
"plt.show()\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 绘制状态 x2 与参考轨迹 r2 的对比图 (对应论文 Fig. 7)\n",
|
|||
|
|
"plt.figure(figsize=(10, 5))\n",
|
|||
|
|
"plt.plot(time_steps, x2_history, label='x2 (System State)', color='red', linewidth=1.5)\n",
|
|||
|
|
"plt.plot(time_steps, r2_history, label='r2 (Reference)', color='green', linestyle='--', linewidth=1.5)\n",
|
|||
|
|
"plt.title('Tracking Performance: x2 and r2')\n",
|
|||
|
|
"plt.xlabel('Time step')\n",
|
|||
|
|
"plt.ylabel('x2 and r2')\n",
|
|||
|
|
"plt.legend(loc='upper right')\n",
|
|||
|
|
"plt.grid(True)\n",
|
|||
|
|
"plt.xlim([0, eval_steps])\n",
|
|||
|
|
"plt.show()\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 新增：绘制误差曲线 e1 和 e2\n",
|
|||
|
|
"plt.figure(figsize=(10, 5))\n",
|
|||
|
|
"plt.plot(time_steps, e1_history, label='e1 = x1 - r1', color='blue', linewidth=1.5)\n",
|
|||
|
|
"plt.plot(time_steps, e2_history, label='e2 = x2 - r2', color='orange', linewidth=1.5)\n",
|
|||
|
|
"plt.axhline(0.0, color='black', linestyle='--', linewidth=1.0)\n",
|
|||
|
|
"plt.title('Tracking Errors: e1 and e2')\n",
|
|||
|
|
"plt.xlabel('Time step')\n",
|
|||
|
|
"plt.ylabel('Tracking error')\n",
|
|||
|
|
"plt.legend(loc='upper right')\n",
|
|||
|
|
"plt.grid(True)\n",
|
|||
|
|
"plt.xlim([0, eval_steps])\n",
|
|||
|
|
"plt.show()"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"metadata": {
|
|||
|
|
"kernelspec": {
|
|||
|
|
"display_name": "hybrid_engine",
|
|||
|
|
"language": "python",
|
|||
|
|
"name": "python3"
|
|||
|
|
},
|
|||
|
|
"language_info": {
|
|||
|
|
"codemirror_mode": {
|
|||
|
|
"name": "ipython",
|
|||
|
|
"version": 3
|
|||
|
|
},
|
|||
|
|
"file_extension": ".py",
|
|||
|
|
"mimetype": "text/x-python",
|
|||
|
|
"name": "python",
|
|||
|
|
"nbconvert_exporter": "python",
|
|||
|
|
"pygments_lexer": "ipython3",
|
|||
|
|
"version": "3.10.20"
|
|||
|
|
}
|
|||
|
|
},
|
|||
|
|
"nbformat": 4,
|
|||
|
|
"nbformat_minor": 5
|
|||
|
|
}
|