Jack Huang's Blog


  • 首页

  • 标签

  • 归档

  • 搜索

策略梯度方法笔记

发表于 2020-02-10 | 更新于 2021-06-16

使用强化学习实现机器人的连续控制,策略梯度方法是首选。下面即对强化学习中策略梯度方法进行总结。

核心概念

免模型学习(Model-Free) vs 有模型学习(Model-Based)

不同强化学习算法最重要的区分点之一就是智能体是否能完整了解或学习到所在环境的模型。 环境的模型是指一个预测状态转换和奖励的函数。

有模型学习最大的优势在于智能体能够 提前考虑来进行规划,走到每一步的时候,都提前尝试未来可能的选择,然后明确地从这些候选项中进行选择。智能体可以把预先规划的结果提取为学习策略。这其中最著名的例子就是 AlphaZero。这个方法起作用的时候,可以大幅度提升采样效率 —— 相对于那些没有模型的方法。

有模型学习最大的缺点就是智能体往往不能获得环境的真实模型。如果智能体想在一个场景下使用模型,那它必须完全从经验中学习,这会带来很多挑战。最大的挑战就是,智能体探索出来的模型和真实模型之间存在误差,而这种误差会导致智能体在学习到的模型中表现很好,但在真实的环境中表现得不好(甚至很差)。基于模型的学习从根本上讲是非常困难的,即使你愿意花费大量的时间和计算力,最终的结果也可能达不到预期的效果。

使用模型的算法叫做有模型学习,不基于模型的叫做免模型学习。虽然免模型学习放弃了有模型学习在样本效率方面的潜在收益,但是他们往往更加易于实现和调整。

同策略(on-policy) vs 异策略(off-policy)

异策略(off-policy)的代表算法Q-learning,亦称SarasMax,其采样的策略(用于执行,behavior policy) 和更新Q值的策略(用于评估,target policy)不一样,行为策略为贪心策略,而target policy为确定性策略,即选择最Q值最优的action。

同策略(on-policy)的代表算法Sarsa,亦称on-line Q-learning,其采样的策略(用于执行,behavior policy) 和更新Q值的策略(用于评估,target policy)一样,行为策略和目标策略均为贪心策略。Sarsa的每次Q值更新需要知道前一步的状态(state)、前一步的动作(action)、奖赏值(reward)、当前状态(state)、将要执行的动作(action),由此得名Sarsa算法。

在线(online) vs 离线(offline)

在计算机科学中,在线机器学习是一种机器学习的方法,其中数据按顺序可用,并且用于在每个步骤中为将来的数据更新我们的最佳预测器,而不是通过学习生成最佳预测器的批处理学习技术 一次对整个训练数据集。 在线学习是机器学习领域中的一种常用技术,在该领域中,计算无法训练整个数据集是不可行的,因此需要核心算法。 它也用于算法必须动态适应数据中的新模式的情况下,或者在数据本身随时间而变的情况下(例如,股价预测)使用。 在线学习算法可能易于遭受灾难性干扰,这一问题可以通过增量学习方法来解决。

在机器学习中,采用离线学习的系统在初始训练阶段完成后不会改变其对目标函数的近似值。这些系统通常也是渴望学习的示例。

在在线学习中,只有一组可能的元素是已知的,而在离线学习中,学习者则知道这些元素的标识以及它们显示的顺序。

强化学习分类

强化学习简单分类

图1 强化学习简单分类

在机器人学习领域,目前主要有三类有效的免模型的深度强化学习算法:

  • TRPO,PPO
  • DDPG及其拓展(D4PG,TD3等)
  • Soft Q-Learning, Soft Actor-Critic

PPO算法是TRPO(Trust Region Policy Optimization)算法的近似,该算法更能适应大规模的运算,是目前最主流的DRL算法,同时面向离散控制和连续控制,在OpenAI Five上取得了巨大成功。但是PPO是一种on-policy的算法,也就是PPO面临着严重的sample inefficiency,需要巨量的采样才能学习,这对于真实的机器人训练来说,是无法接受的。

DDPG及其拓展则是DeepMind开发的面向连续控制的off policy算法,相对PPO 更sample efficient。DDPG训练的是一种确定性策略deterministic policy,即每一个state下都只考虑最优的一个动作。

Soft Actor-Critic (SAC)是面向Maximum Entropy Reinforcement learning 开发的一种off policy算法,和DDPG相比,Soft Actor-Critic使用的是随机策略stochastic policy,相比确定性策略具有一定的优势(具体后面分析)。Soft Actor-Critic在公开的benchmark中取得了非常好的效果,并且能直接应用到真实机器人上。

符号定义

下表给出强化学习常用符号定义。

符号 含义
$s \in \mathcal{S}$ 状态。
$a \in \mathcal{A}$ 动作。
$r \in \mathcal{R}$ 回报。
$S_{t}, A_{t}, R_{t}$ 一个轨迹中第t个时间步对应的状态、动作以及回报。我可能会偶尔使用$s_t,a_t,r_t$来代替。
$\gamma$ 折扣因子;用于惩罚未来回报中的不确定性;$0<γ≤1$。
$G_{t}$ 累积回报;或者说累积折扣回报;$G_{t}=\sum_{k=0}^{\infty} \gamma^{k} R_{t+k+1}$。
$P\left(s^{\prime}, r\vert s, a\right)$ 在当前状态s下采取动作a后转移到下一个状态 s′ 并得到回报 r 的概率。
$\pi(a\vert s)$ 随机策略(智能体行为逻辑);$\pi_{\theta}( .)$代表由θ参数化的策略。
$μ(s)$ 确定性策略;虽然也可以把确定性策略记为$π(s)$,但是采用一个不同的字母可以让我们更容易分辨一个策略到底是确定性的还是随机的。π或者μ都是强化学习算法要学习的目标。
$V(s)$ 状态-值函数衡量状态s的期望累积回报;$V_{w}( .)$代表由w参数化的状态-值函数。
$V^{\pi}(s)$ 当智能体遵循策略π时状态s的期望累积回报;$V^{\pi}(s)=\mathbb{E}{a \sim \pi}\left[G{t}\vert S_{t}=s\right]$
$Q(s,a)$ 动作-值函数,与状态-值函数类似,但是它衡量在状态s下采取动作a后的期望累积回报;$Q_{w}( .)$代表由w参数化的动作-值函数。
$Q^{\pi}(s, a)$ 与$V^{\pi}(s)$类似,当智能体遵循策略π时,在状态s下采取动作a后的期望累积回报;$Q^{\pi}(s, a)=\mathbb{E}{a \sim \pi}\left[G{t}\vert S_{t}=s,A_{t}=a\right]$
$A(s, a)$ 优势函数,$A(s,a)=Q(s,a)−V(s)$;可以认为优势函数是加强版本的动作-值函数,但是由于它采用状态-值函数作为基准使得它具有更小的方差。

策略梯度方法

强化学习的目标是为智能体找到一个最优的行为策略从而获取最大的回报。策略梯度方法主要特点在于直接对策略进行建模并优化。策略通常被建模为由θ参数化的函数$\pi_{\theta}(a | s)$。回报(目标)函数的值受到该策略的直接影响,因而可以采用很多算法来对θ进行优化来最大化回报(目标)函数。

回报(目标)函数定义如下:
$$ J(\theta)=E_{\tau \sim \pi_{\theta}} [R(\tau)]=\sum_{s \in \mathcal{S}} d^{\pi}(s) V^{\pi}(s)=\sum_{s \in \mathcal{S}} d^{\pi}(s) \sum_{a \in \mathcal{A}} \pi_{\theta}(a | s) Q^{\pi}(s, a) $$

其中$d^{\pi}(s)$代表由$\pi_{\theta}$引出的马尔科夫链的平稳分布(π下的在线策略状态分布)。

使用梯度上升方法,我们可以将参数 $\theta$ 往梯度 $\nabla_{\theta} J(\theta)$ 给出的方向进行改变从而去找到最优的 $\theta$ 使得其对应的策略 $\pi_{\theta}$ 能够给智能体带来最大的期望累积回报。

$$\theta_{k+1} = \theta_k + \alpha \left. \nabla_{\theta} J(\pi_{\theta}) \right|_{\theta_k}.$$

策略性能的梯度 $\nabla_{\theta} J(\pi_{\theta})$ ,通常被称为 策略梯度 ,优化策略的算法通常被称为 策略算法 。

策略梯度定理

$$\begin{aligned} \nabla_\theta J(\theta) &\propto \sum_{s \in \mathcal{S}} d^\pi(s) \sum_{a \in \mathcal{A}} Q^\pi(s, a) \nabla_\theta \pi_\theta(a \vert s) &\ &= \sum_{s \in \mathcal{S}} d^\pi(s) \sum_{a \in \mathcal{A}} \pi_\theta(a \vert s) Q^\pi(s, a) \frac{\nabla_\theta \pi_\theta(a \vert s)}{\pi_\theta(a \vert s)} &\ &= \mathbb{E}_\pi [Q^\pi(s, a) \nabla_\theta \ln \pi_\theta(a \vert s)] & \scriptstyle{\text{; 因为 } (\ln x)’ = 1/x} \end{aligned}$$

$\mathbb{E}{\pi}$代表$\mathbb{E}{s \sim d_{\pi}, a \sim \pi_{\theta}}$,下标表示遵循策略$\pi_{\theta}$(在线策略)时状态以及动作的分布。

深度确定性策略梯度 (DDPG)

DDPG(Lillicrap, et al., 2015)是深度确定性策略梯度(Deep Deterministic Policy Gradient)的缩写,是一个结合了DPG以及DQN的无模型离线演员-评论家算法。DQN(深度Q网络)通过经验回访以及冻结目标网络的方式来稳定Q函数的训练过程。原始的DQN算法只能在离散的动作空间上使用,DDPG算法在学习一个确定性策略的同时通过演员-评论家框架将其扩展到连续的动作空间中。

深度确定性策略梯度算法伪代码

图2 深度确定性策略梯度算法伪代码

近似策略优化PPO

DQN

SAC

参考链接

  1. 策略梯度方法,by Abracadabra.
  2. A (Long) Peek into Reinforcement Learning,by Lilian Weng.
  3. 第三部分:策略优化介绍,by spinningup.
  4. 深度强化学习研究笔记,by jackhuang.
  5. 异策略(Q-learning) v.s. 同策略(Sarsa),by MOMO.
  6. Online_machine_learning,by wikipedia.
  7. 最前沿:深度解读Soft Actor-Critic 算法,by Flood Sung.
  8. 重要性采样(Importance Sampling),by 时雨.
  9. TRPO论文推导,by Ja1r0.
  10. 强化学习进阶 第七讲 TRPO,by 天津包子馅儿.
  11. 强化学习–信赖域系方法:TRPO、PPO,by 秋曾万.
  12. 强化学习(8)——DQN,by 自由而无用.
  13. 理解策略梯度算法,by SIGAI.
  14. SAC论文解读以及简易代码复现,by 已注销.
  15. PPO(Proximal Policy Optimization)近端策略优化算法,by shura_R.
  16. Policy Gradient Algorithms,by lilianweng.
  17. TRPO论文推导,by Ja1r0.
  18. 第九章:连续动作空间的确定性策略,by anesck.

科学研究的标准流程

发表于 2020-02-09

科学研究是每一个科研人必备的技能,那么科学研究应如何入手呢?通常科学研究应遵循如下标准流程:

问题牵引

形式化定义问题

问题等价转换

建立模型

已有数学模型

类似数学模型

重建数学模型

模型求解

收集数据

求解模型

验证模型

模型预测

灵敏度分析

XML解析入门

发表于 2020-01-30

最近在研究编写飞行动力学模型,发现需要使用很多用于插值的数据,这些数据可以是一维向量、二维表格或三维数据。在代码中直接硬编码存储是不合适的,降低程序的灵活性。直接使用文本文档存储也不合适,这些插值数据明显具有结构化的特征。于是想到用XML来存储表示这些数据。下面总结介绍XML解析相关知识。

XML简介

可扩展标记语言(英语:Extensible Markup Language,简称:XML)是一种标记语言。标记指计算机所能理解的信息符号,通过此种标记,计算机之间可以处理包含各种信息的文章等。如何定义这些标记,既可以选择国际通用的标记语言,比如HTML,也可以使用像XML这样由相关人士自由决定的标记语言,这就是语言的可扩展性。XML是从标准通用标记语言(SGML)中简化修改出来的。它主要用到的有可扩展标记语言、可扩展样式语言(XSL)、XBRL和XPath等。

XML结构

XML定义结构、存储信息、传送信息。下例为小张发送给大元的便条,存储为XML。

1
2
3
4
5
6
7
<?xml version="1.0"?>
<小纸条>
<收件人>大元</收件人>
<發件人>小張</發件人>
<主題>問候</主題>
<具體內容>早啊,飯吃了沒? </具體內容>
</小纸条>

每个XML文档都由XML序言开始,在前面的代码中的第一行就是XML序言,。这一行代码会告诉解析器或浏览器这个文件应该按照XML规则进行解析。

但是,根元素到底叫<小纸条>还是<小便条>,则是由文档类型定义(DTD)或XML纲要(XML Schema)定义的。如果DTD规定根元素必须叫<小便条>,那么若写作<小纸条>就不符合要求。这种不符合DTD或XML纲要的要求的XML文档,被称作不合法的XML,反之则是合法的XML。

XML文件的第二行并不一定要包含文档元素;如果有注释或者其他内容,文档元素可以迟些出现。

XML解析器

C++类型XML解析器有:

  • Boost.PropertyTree - A property tree parser/generator that can be used to parse XML/JSON/INI/Info files. [Boost]
  • Expat - An XML parser library written in C. [MIT]
  • Libxml2 - The XML C parser and toolkit of Gnome. [MIT]
  • libxml++ - An XML Parser for C++. [LGPL2]
  • Mini-XML - A small XML parsing library written in ANSI C. [LGPL2 with exceptions]
  • PugiXML - A light-weight, simple and fast XML parser for C++ with XPath support. [MIT]
  • RapidXml - An attempt to create the fastest XML parser possible, while retaining useability, portability and reasonable W3C compatibility. [Boost]
  • TinyXML - A simple, small, minimal, C++ XML parser that can be easily integrating into other programs. [zlib]
  • TinyXML2 - A simple, small, efficient, C++ XML parser that can be easily integrating into other programs. [zlib]
  • TinyXML++ - A completely new interface to TinyXML that uses MANY of the C++ strengths. Templates, exceptions, and much better error * handling. [MIT]
  • Xerces-C++ - A validating XML parser written in a portable subset of C++. [Apache2]

推荐使用TinyXML2。

参考链接

  1. C++解析xml有什么好用的轮子?,by 知乎.
  2. awesome-cpp xml,by fffaraz.
  3. XML与C++对象的相互转化,by Mr_John_Liang.
  4. JSON与XML的区别比较,by SanMaoSpace.

软件开发文档的编写方法

发表于 2020-01-22 | 更新于 2022-07-07

软件开发文档是软件开发过程的输出产物。软件开发过程的不同阶段将产生不同的软件开发文档。例如:软件需求分析阶段将产生软件需求规格说明书,软件概要设计阶段将产生概要设计说明书,软件详细设计阶段将产生详细设计说明书。按照软件工程的原则,软件开发过程输出这些文档的目的是为了保障软件开发的质量,确保软件项目能够按时完成,并保质保量。下面重点介绍各类软件开发文档的编写方法。

软件过程模型

软件过程模型是软件过程的简化表示。典型的软件过程模型有:瀑布模型、增量式开发模型和面向服用的软件工程模型。以瀑布模型为例,其涉及的开发活动如图1所示。

瀑布模型

图1 瀑布模型

各个开发活动对应产出的软件开发文档主要有:

  1. 可行性研究报告
  2. 项目开发计划
  3. 软件需求说明书
  4. 概要设计说明书
  5. 详细设计说明书
  6. 数据库设计说明书
  7. 数据要求说明书
  8. 测试计划
  9. 测试分析报告
  10. 项目开发总结报告
  11. 操作手册
  12. 用户手册
  13. 开发进度月报

软件开发文档

可行性研究报告

可行性研究报告

图2 可行性研究报告

项目开发计划

项目开发计划

图3 项目开发计划

软件需求说明书

软件需求说明书

图4 软件需求说明书

概要设计说明书

概要设计说明书

图5 概要设计说明书

详细设计说明书

详细设计说明书

图6 详细设计说明书

数据库设计说明书

数据库设计说明书

图7 数据库设计说明书

数据要求说明书

数据要求说明书

图8 数据要求说明书

测试计划

测试计划

图9 测试计划

测试分析报告

测试分析报告

图10 测试分析报告

项目开发总结报告

项目开发总结报告

图11 项目开发总结报告

操作手册

操作手册

图12 操作手册

用户手册

用户手册

图13 用户手册

软件开发文档的使用

软件文档分类

软件文档分类

图14 软件文档分类

软件文档读者

软件文档读者

图15 软件文档读者

软件文档使用

软件文档使用

图16 软件文档使用

参考链接

  1. 软件工程,by wikipedia.
  2. 软件需求,概要设计,详细设计(文档)怎么做,做什么?,by 安东尼_Anthony.
  3. 软件工程文档总结,by BONIC.
  4. 国标:计算机软件文档编制规范,by 宋哥.
  5. 软件测试流程,by HenryZ.Tang.

CPlusPlus不常用语法解析

发表于 2020-01-21

近年来C++发展很快,出现了一些新的语法和特性。熟练掌握这些语法和特性,可提高编写C++代码的效率。下面即简要介绍这些C++语法和特性。

const=0

在类声明中,会出现const=0语法,如下所示:

1
2
3
4
5
class Weapon
{
public:
virtual void attack() const = 0;
};

在此处 =0 说明该类成员函数是一个纯虚函数。而将const放在成员函数之后,表示该成员函数禁止修改该类的数据成员(mutable成员除外)。如果您无意中修改了该类的数据成员,编译器会报告一个错误。

参考链接

  1. 关于virtual:c ++:const = 0的方法原型的代码说明,by 码农家园.
  2. C++构造函数和析构函数的调用顺序,by 靖心.

CPlusPlus单元测试框架Catch入门

发表于 2020-01-15 | 更新于 2022-07-11

最近在编写一个飞行力学的类库,随着类数量的增加,代码越来越复杂,质量越来越难以控制,因此引入单元测试,通过自动化测试以保障代码质量,防止因代码修改引入新的Bug。C++已经有一些成熟的代码测试框架,例如:Google Test, Boost.Test, CppUnit, Cute,等等。通过综合分析和比较,最终选择Catch2测试框架。选择该测试框架的原因是其够轻量级,够简单。

Catch2简介

Catch2是轻量级的C++的多范式测试框架。 它也支持Objective-C(也许还有C)。 它主要作为单个头文件分发,尽管某些扩展可能需要其他头文件。

关键特征

  • 快速且非常容易上手。 只需下载catch.hpp,#include它就可以了。
  • 没有外部依赖性。 只要您可以编译C ++ 11并拥有C ++标准库即可。
  • 将测试用例编写为自注册函数(或方法,如果您愿意的话)。
  • 将测试用例划分为几个部分,每个部分都是独立运行的(消除了对夹具的需求)。
  • 使用BDD样式的“时准时限”部分以及传统的单元测试用例。
  • 仅一个核心声明宏可以进行比较。 使用标准C / C ++运算符进行比较-但是完整的表达式已分解,并且记录了lhs和rhs值。
  • 测试使用自由格式的字符串命名-合法标识符中没有其他名称。

核心特征

  • 可以对测试进行标记,以方便地运行临时的测试组。
  • 故障可能(可选)进入Windows和Mac上的调试器。
  • 输出通过模块化报告器对象。 包括基本的文本和XML报告程序。 自定义记者可以轻松添加。
  • 支持JUnit xml输出以与第三方工具(例如CI服务器)集成。
  • 提供了默认的main()函数,但您可以提供自己的控件来进行完全控制(例如,集成到自己的测试运行器GUI中)。
  • 提供了命令行解析器,如果您选择提供自己的main()函数,该解析器仍可以使用。
  • Catch可以测试自己。
  • 替代性断言宏报告失败,但不中止测试用例
  • 浮点公差比较是使用表达性的Approx()语法构建的。
  • 内部和友好的宏是隔离的,因此可以管理名称冲突
  • 匹配器

Catch示例

使用Catch进行单元测试,只需简单掌握TEST_CASE、REQUIRE、SECTION三个宏即可编写绝大部分的测试用例。简单示例如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
TEST_CASE( "vectors can be sized and resized", "[vector]" ) {

std::vector<int> v( 5 );

REQUIRE( v.size() == 5 );
REQUIRE( v.capacity() >= 5 );

SECTION( "resizing bigger changes size and capacity" ) {
v.resize( 10 );

REQUIRE( v.size() == 10 );
REQUIRE( v.capacity() >= 10 );
}
SECTION( "resizing smaller changes size but not capacity" ) {
v.resize( 0 );

REQUIRE( v.size() == 0 );
REQUIRE( v.capacity() >= 5 );
}
SECTION( "reserving bigger changes capacity but not size" ) {
v.reserve( 10 );

REQUIRE( v.size() == 5 );
REQUIRE( v.capacity() >= 10 );
}
SECTION( "reserving smaller does not change size or capacity" ) {
v.reserve( 0 );

REQUIRE( v.size() == 5 );
REQUIRE( v.capacity() >= 5 );
}
}

上述示例中,对于每个SECTION,TEST_CASE都是从头开始执行的,因此,当我们进入每个部分时,我们知道vector的大小为5,容量至少为5。通过REQUIRE宏在顶层确保vector大小和容量的正确性。这是可行的,因为SECTION宏包含一个if语句,该语句回调Catch来查看是否应执行该节。 通过TEST_CASE,每次运行都会执行一个叶子部分。 其他部分将被跳过。 下次执行下一个部分,依此类推,直到没有新的部分为止。

参考链接

  1. Writing Unit Tests with Catch and CMake,by filebox.
  2. Integrating catch2 with CMake and Jenkins,by mariuselvert.
  3. Catch2,by catchorg.
  4. C++单元测试入门,by pezy.
  5. Catch2 - 用于 test 的轻量级库,by Bluemultipl.

四元数与旋转矩阵

发表于 2020-01-15 | 更新于 2021-04-09

四元数是由爱尔兰数学家威廉·卢云·哈密顿在1843年创立出的数学概念。单位四元数(Unit quaternion)可以用于表示三维空间里的旋转。它与常用的另外两种表示方式(三维正交矩阵和欧拉角)是等价的,但是避免了欧拉角表示法中的万向锁问题。比起三维正交矩阵表示,四元数表示能够更方便地给出旋转的转轴与旋转角。

欧拉角

欧拉角(Euler Angles)是一种描述三维旋转的方式,根据欧拉旋转定理,任何一个旋转都可以用三个旋转的参数来表示。但欧拉角的描述方式有很多种,并没有一个统一标准。对于定义一个欧拉角,需要明确的内容包括:

  1. 三个旋转角的组合方式(是xyz还是yzx还是zxy)
  2. 旋转角度的参考坐标系统(旋转是相对于固定的坐标系还是相对于自身的坐标系)
  3. 使用旋转角度是左手系还是右手系
  4. 三个旋转角的记法

旋转角的记法

顺序 飞行器 望远镜 符号 角速度
第一 heading azimuth $θ$ yaw
第二 attitude elevation $ϕ$ pitch
第三 bank tilt $ψ$ roll

旋转矩阵

对于两个三维点 $p_1(x_1, y_1, z_1)$,$p_2(x_2,y_2,z_2)$,由点 $p_1$ 经过旋转矩阵 $R$ 旋转到 $p_2$,则有:

$$R = \left[ \begin{matrix}
r_{11} & r_{12} & r_{13}\
r_{21} & r_{22} & r_{23}\
r_{31} & r_{32} & r_{33}
\end{matrix} \right]$$

$$\left[ \begin{matrix}
x_2 \
y_2 \
z_2
\end{matrix} \right] = R
\left[ \begin{matrix}
x_1 \
y_1 \
z_1
\end{matrix} \right]$$

注:旋转矩阵为正交矩阵$RR^T=E$

  • 绕x轴旋转:

$$R_x(\theta) = \left[ \begin{matrix}
1 & 0 & 0\
0 & cos\theta & -sin\theta\
0 & sin\theta & cos\theta
\end{matrix} \right]$$

  • 绕y轴旋转:

$$R_y(\theta) = \left[ \begin{matrix}
cos\theta & 0 & sin\theta\
0 & 1 & 0\
-sin\theta & 0 & cos\theta
\end{matrix} \right]$$

  • 绕z轴旋转:

$$R_z(\theta) = \left[ \begin{matrix}
cos\theta & -sin\theta & 0\
sin\theta & cos\theta & 0\
0 & 0 & 1
\end{matrix} \right]$$

  • 任意旋转矩阵:

任何一个旋转可以表示为依次绕着三个旋转轴旋三个角度的组合。这三个角度称为欧拉角。

由欧拉角求旋转矩阵

设三个轴$x,y,z$的欧拉角分别为$θ_x,θ_y,θ_z$,正弦值、余弦值分别为$s_x, c_x, s_y, c_y, s_z, c_z$那么旋转矩阵为:

$$R(\theta_x,\theta_y,\theta_z)=R_z(\theta_z)R_y(\theta_y)R_x(\theta_x) = \left[ \begin{matrix}
c_y c_z & c_z s_x s_y - c_x s_z & s_x s_z + c_x c_z s_y\
c_y s_z & c_x c_z + s_x s_y s_z & c_x s_y s_z - c_z s_z\
-s_y & c_y s_x & c_x c_y
\end{matrix} \right]$$

由旋转矩阵求欧拉角

$$R = \left[ \begin{matrix}
r_{11} & r_{12} & r_{13}\
r_{21} & r_{22} & r_{23}\
r_{31} & r_{32} & r_{33}
\end{matrix} \right] = \left[ \begin{matrix}
c_y c_z & c_z s_x s_y - c_x s_z & s_x s_z + c_x c_z s_y\
c_y s_z & c_x c_z + s_x s_y s_z & c_x s_y s_z - c_z s_z\
-s_y & c_y s_x & c_x c_y
\end{matrix} \right]$$

解方程可得:

$$\theta_x = atan2(r_{32}, r_{33})$$
$$\theta_y = atan2(-r_{31}, \sqrt{r_{32}^2+r_{33}^2})$$
$$\theta_z = atan2(r_{21}, r_{11})$$

注意:atan2()为C++中的函数,atan2(y,x) 求的是y/x的反正切,其返回值为[-pi,+pi]之间的一个数。

四元数

三维空间的任意旋转,都可以用绕三维空间的某个轴旋转过某个角度来表示,即所谓的Axis-Angle表示方法。这种表示方法里,Axis可用一个三维向量$(x,y,z)$来表示,$θ$可以用一个角度值来表示,直观来讲,一个四维向量$(θ,x,y,z)$就可以表示出三维空间任意的旋转。注意,这里的三维向量$(x,y,z)$只是用来表示axis的方向朝向,因此更紧凑的表示方式是用一个单位向量来表示方向axis,而用该三维向量的长度来表示角度值$θ$。这样以来,可以用一个三维向量$(θ∗x,θ∗y,θ∗z)$就可以表示出三维空间任意的旋转,前提是其中$(x,y,z)$是单位向量。这就是旋转向量(Rotation Vector)的表示方式,OpenCV里大量使用的就是这种表示方法来表示旋转(见OpenCV相机标定部分的rvec)。

  • 单位向量(x,y,z)旋转θ角度后的四元数:

$$(cos \frac{\theta}{2}, xsin \frac{\theta}{2}, ysin \frac{\theta}{2}, z*sin \frac{\theta}{2})$$

  • 对于三维坐标的旋转,可以通过四元数乘法直接操作,与旋转矩阵操作可以等价,但是表示方式更加紧凑,计算量也可以小一些。

四元数求旋转矩阵

已知四元数:

$$\mathbf{q} = q_0 + q_1 i + q_2 j + q_3 k = [s, \mathbf{v}]$$

利用Rodrigues公式可以由四元数求得旋转矩阵R:

$$R = \left[ \begin{matrix}
1 - 2 q_2^2 - 2 q_3^2 & 2q_1 q_2 - 2q_0 q_3 & 2 q_1 q_3 + 2 q_0 q_2 \
2q_1 q_2 + 2q_0 q_3 & 1 - 2 q_1^2 - 2 q_3^2 & 2 q_2 q_3 - 2 q_0 q_1 \
2 q_1 q_3 - 2 q_0 q_2 & 2 q_2 q_3 + 2 q_0 q_1 & 1 - 2 q_1^2 - 2 q_2^2
\end{matrix} \right ]$$

旋转矩阵求四元数

给出其中一种情况的计算方法:

$$q_0 = \frac{\sqrt{1+r_{11}+r_{22}+r_{33}}}{2}$$
$$q_1 = \frac{r_{32}-r_{23}}{4q_0}$$
$$q_2 = \frac{r_{13}-r_{31}}{4q_0}$$
$$q_3 = \frac{r_{21}-r_{12}}{4q_0}$$

其中要满足 $q_0 \neq 0$,$1+r_{11}+r_{22}+r_{33}>0$,即 $1+tr(R)>0$

四元数进行姿态变换

假设坐标系O1上的点P1(x1, y1, z1), 存在变换矩阵R, 可计算P1点在坐标系O2上的坐标值为P2(x2, y2, z2):
$$P2=R∗P1$$

矩阵R对应的四元数为q, 则使用四元数计算为:
首先三维空间点用一个虚四元数来描述:P1=[0, x1, y1, z1], P2=[0, x2, y2, z2]
则P2和P1将计算关系为:

$$P_2=qP_1q^{-1}$$

参考链接

  1. 四元数,by wikipedia.
  2. 四元数与空间旋转,by wikipedia.
  3. 从旋转矩阵计算欧拉角,by PengChao.
  4. 旋转变换(一)旋转矩阵,by csxiaoshui.
  5. 旋转矩阵、欧拉角、四元数理论及其转换关系,by jason_ql.
  6. 旋转变换(二)欧拉角,by csxiaoshui.
  7. 欧拉角,by wikipedia.
  8. 欧拉角细节/旋转顺序/内旋外旋,by 能儿.
  9. 四元数运算与姿态变换,by Yoyo_wym.
  10. Flight Parameters and Quaternion Math,by mathwork.

JSBSim源代码分析

发表于 2020-01-13

JSBSim是一个开源跨平台的飞行动力学模型(FDM)软件库,用于模拟航空航天飞行器的飞行动力学。 该库已被纳入飞行模拟软件包FlightGear和OpenEaagles。JSBSim可以独立运行,通过命令行参数指定飞行器和初始状态,进行简单情境下的飞行动力学仿真,也可以将JSBSim作为代码库,编程实现飞行器模型加载,设置输入,获得输出。下面将通过分析JSBSim源代码,研究其实现通用飞行动力学模型的方法。

入口分析

下面是JSBSim参考手册中的最简单实例,因JSBSim的不断开发,JSBSim参考手册中该编程实例有点过时,因此进行了少量修改。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <FGFDMExec.h>
#include <sg_path.hxx>

using namespace std;

int main(int argc, char **argv)
{
JSBSim::FGFDMExec FDMExec;
bool result = true;

FDMExec.LoadScript(SGPath::fromUtf8(argv[1]));

while (result) result = FDMExec.Run();
}

从上述代码可知,调用JSBSim的主要方法是利用FGFDMExec类,通过实例化一个FGFDMExec类,就相当于获得了一个运行JSBSim仿真的工具箱,通过这个工具箱就可以调用JSBSim的大部分功能,实现我们要的仿真目标。同时FGFDMExec类通过加载外部飞机的XML脚本,实现飞行动力学模型的通用性。

JSBSim初始化流程

上述JSBSim最简仿真示例中已包含其初始化流程,采用图示如下:

JSBSim初始化流程

图1 JSBSim初始化流程

FGFDMExec初始化

FGFDMExec类在其构造函数中对各个模型进行初始化,具体代码在Allocate函数中:

1
2
3
4
5
6
7
8
9
10
11
12
13
FGFDMExec::FGFDMExec(FGPropertyManager* root, unsigned int* fdmctr)
: Root(root), RandomEngine(new default_random_engine), FDMctr(fdmctr)
{
...
try {
Allocate();
} catch (const string& msg ) {
cout << "Caught error: " << msg << endl;
exit(1);
}

...
}

Allocate函数代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
bool FGFDMExec::Allocate(void)
{
bool result=true;

Models.resize(eNumStandardModels);

// First build the inertial model since some other models are relying on
// the inertial model and the ground callback to build themselves.
// Note that this does not affect the order in which the models will be
// executed later.
Models[eInertial] = new FGInertial(this);

// See the eModels enum specification in the header file. The order of the
// enums specifies the order of execution. The Models[] vector is the primary
// storage array for the list of models.
Models[ePropagate] = new FGPropagate(this);
Models[eInput] = new FGInput(this);
Models[eAtmosphere] = new FGStandardAtmosphere(this);
...

// Assign the Model shortcuts for internal executive use only.
Propagate = (FGPropagate*)Models[ePropagate];
Inertial = (FGInertial*)Models[eInertial];
Atmosphere = (FGAtmosphere*)Models[eAtmosphere];
...

// Initialize planet (environment) constants
LoadPlanetConstants();

// Initialize models
for (unsigned int i = 0; i < Models.size(); i++) {
// The Input/Output models must not be initialized prior to IC loading
if (i == eInput || i == eOutput) continue;

LoadInputs(i);
Models[i]->InitModel();
}

...
return result;
}

Allocate函数代码中需要注意LoadInputs函数,该函数决定各个子模型的初始化顺序,确定各个子模型的输入输出,具体代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
void FGFDMExec::LoadInputs(unsigned int idx)
{
switch(idx) {
case ePropagate:
Propagate->in.vPQRidot = Accelerations->GetPQRidot();
Propagate->in.vUVWidot = Accelerations->GetUVWidot();
Propagate->in.DeltaT = dT;
break;
case eInput:
break;
case eInertial:
Inertial->in.Position = Propagate->GetLocation();
break;
case eAtmosphere:
Atmosphere->in.altitudeASL = Propagate->GetAltitudeASL();
break;
case eWinds:
Winds->in.AltitudeASL = Propagate->GetAltitudeASL();
Winds->in.DistanceAGL = Propagate->GetDistanceAGL();
...
break;
case eAuxiliary:
Auxiliary->in.Pressure = Atmosphere->GetPressure();
Auxiliary->in.Density = Atmosphere->GetDensity();
...
break;
case eSystems:
// Dynamic inputs come into the components that FCS manages through properties
break;
case ePropulsion:
Propulsion->in.Pressure = Atmosphere->GetPressure();
...

break;
case eAerodynamics:
Aerodynamics->in.Alpha = Auxiliary->Getalpha();
...
break;
case eGroundReactions:
// There are no external inputs to this model.
GroundReactions->in.Vground = Auxiliary->GetVground();
...
break;
case eExternalReactions:
// There are no external inputs to this model.
break;
case eBuoyantForces:
BuoyantForces->in.Density = Atmosphere->GetDensity();
...
break;
case eMassBalance:
MassBalance->in.GasInertia = BuoyantForces->GetGasMassInertia();
MassBalance->in.GasMass = BuoyantForces->GetGasMass();
...
break;
case eAircraft:
Aircraft->in.AeroForce = Aerodynamics->GetForces();
Aircraft->in.PropForce = Propulsion->GetForces();
...
break;
case eAccelerations:
Accelerations->in.J = MassBalance->GetJ();
Accelerations->in.Jinv = MassBalance->GetJinv();
...
break;
default:
break;
}
}

FGFDMExec加载飞机配置

FGFDMExec的LoadScript函数在初始化时负责加载飞机配置,用于初始化各个子模型。

1
2
3
4
5
6
7
8
9
10
bool FGFDMExec::LoadScript(const SGPath& script, double deltaT,
const SGPath& initfile)
{
bool result;

Script = new FGScript(this);
result = Script->LoadScript(GetFullPath(script), deltaT, initfile);

return result;
}

FGFDMExec运行

FGFDMExec的Run函数负责飞行动力学模型的计算,其代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
bool FGFDMExec::Run(void)
{
bool success=true;

...

// returns true if success, false if complete
if (Script != 0 && !IntegrationSuspended()) success = Script->RunScript();

for (unsigned int i = 0; i < Models.size(); i++) {
LoadInputs(i);
Models[i]->Run(holding);
}

...

return success;
}

FGFDMExec的Run函数将依次调用各个子模型的Run函数。

参考链接

  1. JSBSim编程实践之入门,by jackhuang.

排序算法总结

发表于 2020-01-05 | 更新于 2020-05-23

排序算法是计算机科学的基石之一,可从时间复杂度、空间复杂度、稳定性、是否原地排序等维度对排序算法进行分类。下面从时间复杂度方面对排序算法进行分类。

O(n^2)算法

冒泡排序

选择排序

插入排序

O(nlogn)算法

希尔排序

快速排序

归并排序

堆排序

O(n)算法

计数排序

桶排序

基数排序

参考链接

  1. 漫画:“排序算法” 大总结,by 小灰.
  2. 分布式哈希表 (DHT) 和 P2P 技术,by luyuhuang.
  3. Gzip 格式和 DEFLATE 压缩算法,by Luyu Huang.

Pandas入门教程

发表于 2020-01-02 | 更新于 2024-11-18

Pandas是一个开源的,BSD许可的库,为Python编程语言提供高性能,易于使用的数据结构和数据分析工具。

Pandas特色

Pandas 适用于处理以下类型的数据:

  • 与 SQL 或 Excel 表类似的,含异构列的表格数据;
  • 有序和无序(非固定频率)的时间序列数据;
  • 带行列标签的矩阵数据,包括同构或异构型数据;
  • 任意其它形式的观测、统计数据集, 数据转入 Pandas 数据结构时不必事先标记。

Pandas数据结构

Pandas 的主要数据结构是 Series(一维数据)与 DataFrame(二维数据),这两种数据结构足以处理金融、统计、社会科学、工程等领域里的大多数典型用例。对于 R 用户,DataFrame 提供了比 R 语言 data.frame 更丰富的功能。Pandas 基于 NumPy 开发,可以与其它第三方科学计算支持库完美集成。

维数 名称 描述
1 Series 带标签的一维同构数组
2 DataFrame 带标签的,大小可变的,二维异构表格

Pandas用法

Pandas用法与Matlab中矩阵操作很类似,熟悉Matlab操作的同学可以很快上手Pandas。

生成对象

生成Series对象:

1
2
3
4
5
6
7
8
9
10
11
In [3]: s = pd.Series([1, 3, 5, np.nan, 6, 8])

In [4]: s
Out[4]:
0 1.0
1 3.0
2 5.0
3 NaN
4 6.0
5 8.0
dtype: float64

生成DataFrame对象:

1
2
3
4
5
6
7
8
9
10
11
In [7]: df = pd.DataFrame(np.random.randn(6, 4), index=dates, columns=list('ABCD'))

In [8]: df
Out[8]:
A B C D
2013-01-01 0.469112 -0.282863 -1.509059 -1.135632
2013-01-02 1.212112 -0.173215 0.119209 -1.044236
2013-01-03 -0.861849 -2.104569 -0.494929 1.071804
2013-01-04 0.721555 -0.706771 -1.039575 0.271860
2013-01-05 -0.424972 0.567020 0.276232 -1.087401
2013-01-06 -0.673690 0.113648 -1.478427 0.524988

列切片

1
df1 = df[df.columns[0:6]];

重命名列名

1
df1.columns['id','name','sex','age','department','work']

过滤行

1
df1 = df1[df1.iloc[:,1]=='YES']

合并表格

pandas 2.0实现数据的合并与拼接,主要有三种方法:

  • join 最简单,主要用于基于索引的横向合并拼接
  • merge 最常用,主要用于基于指定列的横向合并拼接
  • concat最强大,可用于横向和纵向合并拼接
1
2
# 合并两个表,df1 和 df2 表结构一样
df3 = pd.concat([df1,df2])

毫秒解析

1
2
timeDelta = datetime.datetime(2024,11,11) - datetime.datetime(1900,1,1)
df3['datetime'] = pd.to_datetime(df['datetime'], format='%Y-%m-%d %H:%M:%S:%f') + timeDelta

按时间排序

1
df3 = df3.sort_values(by="datetime")

输出CSV表格

1
df3.to_csv(filePath + 'test.csv',index=0)

参考链接

  1. Pandas 中文,by pypandas.
  2. 十分钟入门 Pandas,by pypandas.
  3. Python读取csv文件的三种方式,by 涛声依旧2019.
  4. Python模块化开发组织代码程序示例,by BabyFish13.
  5. Python最佳实践指南!,by Prodesire.
  6. pandas中DataFrame 数据合并,连接(merge,join,concat) ,by Vincent-yuan.
  7. 【已解决】AttributeError: ‘DataFrame‘ object has no attribute ‘append‘,by zhtstar.
  8. pandas DataFrame数据重命名列名的几种方式,by littleRpl.
  9. Pandas 毫秒级时间解析,by 一定波兮.
上一页1…303132…53下一页

Jack Huang

528 日志
68 标签
© 2025 Jack Huang
由 Hexo 强力驱动
|
主题 — NexT.Muse