app网站开发定制,wordpress编辑器汉,响应式网站模板免费,软件工程师是程序员吗系列文章目录 文章目录 系列文章目录前言一、介绍1.1 CasADi 是什么#xff1f;1.2 帮助与支持1.3 引用 CasADi1.4 阅读本文档 二、获取与安装三、符号框架3.1 符号 SX3.1.1 关于命名空间的说明3.1.2 C 用户注意事项 3.2 DM3.3 符号 MX3.4 SX 和 MX 混合使用3.5 稀疏类3.5.1 获…系列文章目录 文章目录 系列文章目录前言一、介绍1.1 CasADi 是什么1.2 帮助与支持1.3 引用 CasADi1.4 阅读本文档 二、获取与安装三、符号框架3.1 符号 SX3.1.1 关于命名空间的说明3.1.2 C 用户注意事项 3.2 DM3.3 符号 MX3.4 SX 和 MX 混合使用3.5 稀疏类3.5.1 获取并设置矩阵中的元素 3.6 运算操作3.7 属性查询3.8 线性代数3.9 微积分 - 算法微分3.9.1 语法 四、函数对象4.1 调用函数对象4.2 将 MX 转换为 SX4.3 非线性求根问题4.4 初值问题和灵敏度分析4.4.1 创建积分器4.4.2 灵敏度分析 4.5 非线性规划4.5.1 创建 NLP 求解器 4.6 二次规划4.6.1 高级接口4.6.2 低级接口 前言 一、介绍
CasADi 是一款开源软件工具用于数值优化特别是最优控制即涉及微分方程的优化。该项目由 Joel Andersson 和 Joris Gillis 在鲁汶工程大学工程优化中心 (OPTEC) 在读博士生在 Moritz Diehl 的指导下发起。
本文档旨在简要介绍 CasADi。阅读之后您应该能够在 CasADi 的符号框架中制定和处理表达式使用算法微分高效生成导数信息为常微分方程ODE或微分代数方程DAE系统设置、求解和执行正向及辅助敏感性分析以及制定和求解非线性程序NLP问题和最优控制问题OCP。
CasADi 可用于 C、Python 和 MATLAB/Octave性能几乎没有差别。一般来说Python API 的文档最好比 MATLAB API 稍为稳定。C API 也很稳定但对于 CasADi 入门来说并不理想因为文档有限而且缺乏 MATLAB 和 Python 等解释型语言的交互性。MATLAB 模块已成功通过 Octave4.0.2 或更高版本测试。
1.1 CasADi 是什么
CasADi 最初是一个算法微分AD工具使用的语法借鉴了计算机代数系统CAS这也是其名称的由来。虽然算法微分仍是该工具的核心功能之一但其范围已大大扩展增加了对 ODE/DAE 集成和敏感性分析、非线性编程的支持以及与其他数值工具的接口。从目前的形式来看CasADi 是一款基于梯度的数值优化通用工具主要侧重于最优控制而 CasADi 只是一个名称没有任何特殊含义。
需要指出的是CasADi 并不是一个传统的 AD 工具它几乎不需要任何修改就能从现有的用户代码中计算出导数信息。如果您有一个用 C、Python 或 MATLAB/Octave 编写的现有模型您需要准备好用 CasADi 语法重新实现该模型。
其次CasADi 不是计算机代数系统。虽然符号内核确实包含了越来越多的符号表达式操作工具但与合适的 CAS 工具相比这些功能非常有限。
最后CasADi 并不是一个 “最优控制问题求解器”它不允许用户输入 OCP然后再给出解决方案。相反CasADi 试图为用户提供一套 “构件”只需少量编程工作就能高效地实现通用或专用的 OCP 求解器。
1.2 帮助与支持
如果您发现了一些简单的错误或缺少一些您认为我们比较容易添加的功能最简单的方法就是写信到位于 http://forum.casadi.org/ 的论坛。我们会定期检查论坛并尝试尽快回复。我们对这种支持的唯一期望是当您在科学工作中使用 CasADi 时请引用我们的内容参见第 1.3 节。
如果您需要更多帮助我们随时欢迎您与我们进行学术或工业合作。学术合作通常以共同撰写同行评审论文的形式进行而产业合作则包括协商签订咨询合同。如果您对此感兴趣请直接与我们联系。
1.3 引用 CasADi
如果您在发表的科学著作中使用了 CasADi请注明出处
Article{Andersson2018,Author {Joel A E Andersson and Joris Gillis and Greg Hornand James B Rawlings and Moritz Diehl},Title {{CasADi} -- {A} software framework for nonlinear optimizationand optimal control},Journal {Mathematical Programming Computation},Year {2018},
}1.4 阅读本文档
本文档的目的是让读者熟悉 CasADi 的语法并为构建数值优化和动态优化软件提供易于使用的构件。我们的解释主要是程序代码驱动的几乎不提供数学背景知识。我们假定读者已经对优化理论、微分方程初值问题的求解以及相关编程语言C、Python 或 MATLAB/Octave有一定的了解。
我们将在本指南中并列使用 Python 和 MATLAB/Octave 语法并指出 Python 界面更稳定、文档更完善。除非另有说明否则 MATLAB/Octave 语法也适用于 Octave。我们会尽量指出语法不同的情况。为了方便在两种编程语言之间切换我们还在第 10 章列出了主要区别。
二、获取与安装
CasADi 是一款开源工具可在 LGPL 许可下使用LGPL 是一种允许免版税使用的许可允许在商业闭源应用程序中使用该工具。LGPL 的主要限制是如果您决定修改 CasADi 的源代码而不仅仅是在应用程序中使用该工具那么这些修改CasADi 的 “衍生作品”也必须在 LGPL 下发布。
CasADi 的源代码托管在 GitHub 上其核心部分由独立的 C 代码编写只依赖 C 标准库。它与 Python 和 MATLAB/Octave 的前端功能齐全使用 SWIG 工具自动生成。这些前端不太可能导致明显的效率损失。CasADi 可在 Linux、OS X 和 Windows 上使用。
如需最新的安装说明请访问 CasADi 的安装部分http://install.casadi.org/。
pip install casadi三、符号框架
CasADi 的核心是一个自足的符号框架允许用户使用受 MATLAB 启发的 一切皆矩阵 语法构建符号表达式即矢量被视为 n-by-1 矩阵标量被视为 1-by-1 矩阵。所有矩阵都是稀疏的并使用通用稀疏sparse格式–压缩列存储compressed column storageCCS–来存储矩阵。下面我们将介绍这一框架的最基本类别。
3.1 符号 SX
SX 数据类型用于表示矩阵其元素由一系列一元和二元运算组成的符号表达式构成。要了解其实际运行情况请启动一个交互式 Python shell例如在 Linux 终端或 Spyder 等集成开发环境中输入 ipython或启动 MATLAB 或 Octave 的图形用户界面。假设 CasADi 已正确安装则可以按如下方式将符号导入工作区
from casadi import *现在使用语法创建一个变量 x
x MX.sym(x)这将创建一个 1-by-1 矩阵即一个包含名为 x 的符号基元的标量。多个变量可以有相同的名称但仍然是不同的。标识符就是返回值。你也可以通过向 SX.sym 提供额外参数来创建向量或矩阵值的符号变量
y SX.sym(y,5)
Z SX.sym(Z,4,2)[y_0, y_1, y_2, y_3, y_4]
[[Z_0, Z_4], [Z_1, Z_5], [Z_2, Z_6], [Z_3, Z_7]]分别创建一个 5×1 矩阵即向量和一个 4×2 矩阵的符号基元。
SX.sym 是一个返回 SX 实例的静态函数。在声明变量后表达式就可以直观地形成了
f x**2 10
f sqrt(f)
print(f)sqrt((10sq(x)))您也可以在不使用任何符号基元的情况下创建常量 SX 实例 B1 SX.zeros(4,5) 一个全为零的 4 乘 5 的密集空矩阵 B2 SX(4,5) 一个全部为零的稀疏 4×5 空矩阵 B4 SX.eye(4) 对角线上有 1 的稀疏 4×4 矩阵
B1: 10,
[[1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1]]
B2:
[[00, 00, 00, 00, 00], [00, 00, 00, 00, 00], [00, 00, 00, 00, 00], [00, 00, 00, 00, 00]]
B4: 11,
[[1, 00, 00, 00], [00, 1, 00, 00], [00, 00, 1, 00], [00, 00, 00, 1]]请注意带有结构零的稀疏矩阵与带有实际零的稠密矩阵之间的区别。打印带有结构零的表达式时这些零将表示为 00以区别于实际零 0。
以下列表总结了构建新 SX 表达式的最常用方法 SX.sym(name,n,m) 创建一个 n-m 的符号基元 SX.zeros(n,m) 创建一个 n-m 全为零的密集矩阵 SX(n,m) 创建一个 n-m 结构零的稀疏矩阵 SX.ones(n,m) 创建一个 n-m 全为 1 的密集矩阵 SX.eye(n) 创建一个 n-n 对角矩阵对角线上为 1其他地方为结构零。 SX(scalar_type) 创建一个标量1 乘 1 矩阵其值由参数给出。此方法可以显式使用例如 SX(9)也可以隐式使用例如 9 * SX.nes(2,2)。 SX(matrix_type) 以 NumPy 或 SciPy 矩阵在 Python 中或密集或稀疏矩阵在 MATLAB/Octave 中的形式创建矩阵。例如在 MATLAB/Octave 中SX([1,2,3,4]) 表示行向量SX([1;2;3;4]) 表示列向量SX([1,2;3,4]) 表示 2×2 矩阵。这种方法可以显式或隐式使用。 repmat(v,n,m) 重复表达式 v n 纵向重复 m repmat(SX(3),2,1) 将创建一个所有元素为 3 的 2 乘 1 矩阵。 (仅限 PythonSX(list) 创建列向量 (n-1例如 SX([1,2,3,4])注意 Python 列表与 MATLAB/Octave 水平连接之间的区别两者都使用方括号语法。 (仅限 PythonSX(列表的列表) 用列表中的元素创建密集矩阵例如 SX([[1,2],[3,4]])或行向量1-by n 矩阵。
3.1.1 关于命名空间的说明
在 MATLAB 中如果省略了 import casadi.* 命令您仍然可以使用 CasADi方法是在所有符号前加上软件包名称例如用 casadi.SX 代替 SX前提是路径中包含 casadi 软件包。出于排版原因我们在下文中不会这样做但请注意在用户代码中这样做通常更为可取。在 Python 中这种用法相当于发布 import casadi 而不是 from casadi import *。
遗憾的是Octave4.0.3 版没有实现 MATLAB 的 import 命令。为了解决这个问题我们提供了一个简单的函数 import.m可以放在 Octave 的路径中从而实现本指南中使用的紧凑语法。
3.1.2 C 用户注意事项
在 C 中所有公共符号都定义在 casadi 命名空间中并要求包含 casadi/casadi.hpp 头文件。上述命令相当于
#include casadi/casadi.hpp
using namespace casadi;
int main() {SX x SX::sym(x);SX y SX::sym(y,5);SX Z SX::sym(Z,4,2)SX f pow(x,2) 10;f sqrt(f);std::cout f: f std::endl;return 0;
}3.2 DM
DM 与 SX 非常相似但不同之处在于非零元素是数值而不是符号表达式。语法也相同但 SX.sym 等函数没有对应的语法。
DM 主要用于在 CasADi 中存储矩阵以及作为函数的输入和输出。它不用于计算密集型计算。为此请使用 MATLAB 中的内置密集或稀疏数据类型、Python 中的 NumPy 或 SciPy 矩阵或基于表达式模板的库如 C 中的 eigen、ublas 或 MTL。类型之间的转换通常很简单
C DM(2,3)C_dense C.full()
from numpy import array
C_dense array(C) # equivalentC_sparse C.sparse()
from scipy.sparse import csc_matrix
C_sparse csc_matrix(C) # equivalentSX 的更多使用示例可在 http://install.casadi.org/ 的示例包中找到。有关该类及其他特定函数的文档请在 http://docs.casadi.org/ 上查找 “C API”并搜索有关 casadi::Matrix 的信息。
3.3 符号 MX
让我们用上面的 SX 来执行一个简单的操作
x SX.sym(x,2,2)
y SX.sym(y)
f 3*x y
print(f)
print(f.shape)13,
[[((1*x_0)y), ((1*x_2)y)], [((1*x_1)y), ((1*x_3)y)]]
(2, 2)可以看到该操作的输出是一个 2×2 矩阵。请注意乘法和加法是按元素进行的并且为结果矩阵的每个条目创建了新的表达式SX 类型。
现在我们将介绍第二种更通用的矩阵表达式类型 MX。与 SX 类型一样MX 类型允许建立由一系列基本运算组成的表达式。但与 SX 不同的是这些基本运算并不局限于标量一元运算或二元运算 ( R → R \mathbb{R}\to\mathbb{R} R→R 或 R × R → R \mathbb{R}\times\mathbb{R}\to\mathbb{R} R×R→R). 相反用于形成 MX 表达式的基本运算可以是一般的多稀疏矩阵值输入、多稀疏矩阵值输出函数 R n 1 × m 1 × . . ⋅ min R n N × m N → R p 1 × q 1 × . . ⋅ × R p M × q M . \mathbb{R}^{n_{1}\times m_{1}}\times..\cdot\min{\mathbb{R}^{n_{N}\times m_{N}}}\to\mathbb{R}^{p_{1}\times q_{1}}\times..\cdot\times\mathbb{R}^{p_{M}\times q_{M}}. Rn1×m1×..⋅minRnN×mN→Rp1×q1×..⋅×RpM×qM.
MX 的语法与 SX 相同
x MX.sym(x,2,2)
y MX.sym(y)
f 3*x y
print(f)
print(f.shape)((3*x)y)
(2, 2)请注意使用 MX 符号结果只包含两个运算一个乘法运算和一个加法运算而 SX 符号则包含八个运算结果矩阵中每个元素两个运算。因此在处理具有许多元素的天然向量或矩阵值的运算时MX 更经济。正如我们将在第 4 章看到的MX 的通用性也更强因为我们允许调用无法用基本运算展开的任意函数。
MX 支持获取和设置元素使用的语法与 SX 相同但实现方式却截然不同。例如测试打印 2×2 符号变量左上角的元素
x MX.sym(x,2,2)
print(x[0,0])x[0]输出应理解为等于 x 的第一个即 C 中的索引 0结构非零元素的表达式这与上述 SX 案例中的 x_0 不同后者是矩阵第一个索引 0位置的符号基元的名称。
在尝试设置元素时也会出现类似的结果
x MX.sym(x,2)
A MX(2,2)
A[0,0] x[0]
A[1,1] x[0]x[1]
print(A:, A)A: (project((zeros(2x2,1nz)[0] x[0]))[1] (x[0]x[1]))对输出结果的解释是从一个全零稀疏矩阵开始一个元素被分配到 x_0。然后将其投影到不同稀疏度的矩阵中再将另一个元素赋值给 x_0x_1。
刚刚所见的元素访问和赋值都是可用于构造表达式的操作示例。其他操作示例包括矩阵乘法、转置、连接、调整大小、重塑和函数调用。
3.4 SX 和 MX 混合使用
在同一个表达式图形中不能将 SX 对象与 MX 对象相乘也不能执行任何其他操作将两者混合。不过您可以在 MX 图形中调用由 SX 表达式定义的函数。第 4 章将对此进行演示。混合使用 SX 和 MX 通常是个好主意因为由 SX 表达式定义的函数每次操作的开销要低得多因此对于自然写成标量操作序列的操作来说速度要快得多。因此SX 表达式旨在用于低级运算例如第 4.4 节中的 DAE 右侧而 MX 表达式则起到粘合剂的作用可用于制定 NLP 的约束函数其中可能包含对 ODE/DAE 积分器的调用也可能因为太大而无法扩展为一个大表达式。
3.5 稀疏类
如上所述CasADi 中的矩阵使用压缩列存储CCS格式存储。这是稀疏矩阵的标准格式可以高效地进行线性代数运算如元素运算、矩阵乘法和转置。在 CCS 格式中稀疏性模式使用维数行数和列数和两个向量进行解码。第一个向量包含每列第一个结构非零元素的索引第二个向量包含每个非零元素的行索引。有关 CCS 格式的更多详情请参阅 Netlib 上的 “线性系统求解模板”。请注意CasADi 对稀疏矩阵和密集矩阵都使用 CCS 格式。
CasADi 中的稀疏性模式以稀疏性类实例的形式存储稀疏性类是按引用计数的这意味着多个矩阵可以共享相同的稀疏性模式包括 MX 表达式图以及 SX 和 DM 实例。稀疏性类也是缓存的这意味着可以避免创建相同稀疏性模式的多个实例。
以下列表总结了构建新稀疏性模式的最常用方法
Sparsity.dense(n,m) 创建一个密集的 n-m 密度模式
稀疏性n,m 创建稀疏 n-m 稀疏模式
Sparsity.diag(n) 创建对角线 n-n 的对角线
Sparsity.upper(n) 创建一个上三角 n-n 稀疏性模式
Sparsity.lower(n) 创建一个下三角 n-n 稀疏性模式
稀疏性类可用于创建非标准矩阵例如
print(SX.sym(x,Sparsity.lower(3)))[[x_0, 00, 00], [x_1, x_3, 00], [x_2, x_4, x_5]]3.5.1 获取并设置矩阵中的元素
要获取或设置 CasADi 矩阵类型SX、MX 和 DM中的一个元素或一组元素我们在 Python 中使用方括号在 C 和 MATLAB 中使用圆括号。按照这些语言的惯例索引在 C 和 Python 中从 0 开始而在 MATLAB 中则从 1 开始。在 Python 和 C 中我们允许使用负索引来指定从末尾开始计算的索引。在 MATLAB 中使用 end 关键字可以从末尾开始计算索引。
索引可以使用一个索引或两个索引。使用两个索引时可以引用特定行或行集和特定列或列集。使用一个索引时您可以从左上角开始逐列到右下角引用一个元素或一组元素。所有元素都会被计算在内无论它们在结构上是否为零。
M SX([[3,7],[4,5]])
print(M[0,:])
M[0,:] 1
print(M)M SX([[3,7],[4,5]])
print(M[0,:])
M[0,:] 1
print(M)
[[3, 7]]
11,
[[1, 1], [4, 5]]与 Python 的 NumPy 不同CasADi 分片不是左侧数据的视图相反分片访问会复制数据。因此矩阵 m 在下面的示例中完全没有改变
M SX([[3,7],[4,5]])
M[0,:][0,0] 1
print(M)[[3, 7], [4, 5]]下文将详细介绍矩阵元素的获取和设置。讨论适用于 CasADi 的所有矩阵类型。
通过提供行列对或其扁平化索引从矩阵左上角开始按列排列来获取或设置单个元素
M diag(SX([3,4,5,6]))
print(M)M diag(SX([3,4,5,6]))
print(M)
[[3, 00, 00, 00], [00, 4, 00, 00], [00, 00, 5, 00], [00, 00, 00, 6]]print(M[0,0])
print(M[1,0])
print(M[-1,-1])3
00
6分片访问意味着一次设置多个元素。这比一次设置一个元素要有效得多。您可以通过提供一个start , stop , step三元组来获取或设置切片。在 Python 和 MATLAB 中CasADi 使用标准语法
print(M[:,1])
print(M[1:,1:4:2])[00, 4, 00, 00][[4, 00], [00, 00], [00, 6]]在 C 中可以使用 CasADi 的 Slice 辅助类。在上面的例子中这分别意味着 M(Slice(),1) 和 M(Slice(1,-1),Slice(1,4,2)) 。
列表访问与切片访问类似但效率可能低于切片访问
M SX([[3,7,8,9],[4,5,6,1]])
print(M)
print(M[0,[0,3]], M[[5,-6]])[[3, 7, 8, 9], [4, 5, 6, 1]]
[[3, 9]] [6, 7]3.6 运算操作
CasADi 支持大多数标准算术运算如加法、乘法、幂、三角函数等
x SX.sym(x)
y SX.sym(y,2,2)
print(sin(y)-x)[[(sin(y_0)-x), (sin(y_2)-x)], [(sin(y_1)-x), (sin(y_3)-x)]]在 C 和 Python 中但不是在 MATLAB 中标准乘法运算使用 被保留给元素乘法运算在 MATLAB 中为 .。对于矩阵乘法使用 A B 或 (mtimes(A,B) in Python 3.4)
print(y*y, yy)[[sq(y_0), sq(y_2)], [sq(y_1), sq(y_3)]]
[[(sq(y_0)(y_2*y_1)), ((y_0*y_2)(y_2*y_3))], [((y_1*y_0)(y_3*y_1)), ((y_1*y_2)sq(y_3))]]按照 MATLAB 的习惯当参数之一是标量时使用 * 和 .* 的乘法运算是等价的。
转置在 Python 中使用语法 A.T在 C 中使用语法 A.T()在 MATLAB 中使用语法 A’ 或 A.
print(y)
print(y.T)[[y_0, y_2], [y_1, y_3]][[y_0, y_1], [y_2, y_3]]重塑是指改变行数和列数但保留元素的数量和非零点的相对位置。这是一种计算成本很低的操作使用的语法是
x SX.eye(4)
print(reshape(x,2,8))11,
[[1, 00, 00, 00, 00, 1, 00, 00], [00, 00, 1, 00, 00, 00, 00, 1]]连接意味着水平或垂直堆叠矩阵。由于 CasADi 采用列为主的元素存储方式因此水平堆叠矩阵的效率最高。实际上是列向量即由单列组成的矩阵也可以高效地垂直堆叠。在 Python 和 C 中可以使用函数 vertcat 和 horzcat输入参数数量可变进行垂直和水平连接在 MATLAB 中可以使用方括号进行连接
x SX.sym(x,5)
y SX.sym(y,5)
print(vertcat(x,y))[x_0, x_1, x_2, x_3, x_4, y_0, y_1, y_2, y_3, y_4]print(horzcat(x,y))print(horzcat(x,y))
[[x_0, y_0], [x_1, y_1], [x_2, y_2], [x_3, y_3], [x_4, y_4]]这些函数还有一些变种它们的输入是列表Python或单元数组Matlab
L [x,y]
print(hcat(L))[[x_0, y_0], [x_1, y_1], [x_2, y_2], [x_3, y_3], [x_4, y_4]]水平分割和垂直分割是上述水平连接和垂直连接的逆运算。要将一个表达式横向拆分为 n 个较小的表达式除了要拆分的表达式外还需要提供一个长度为 n1 的偏移向量。偏移向量的第一个元素必须是 0最后一个元素必须是列数。其余元素必须按不递减的顺序排列。拆分操作的输出 i 将包含偏移量[i]≤c偏移量[i1]的列 c。下面是语法演示
x SX.sym(x,5,2)
w horzsplit(x,[0,1,2])
print(w[0], w[1])[x_0, x_1, x_2, x_3, x_4] [x_5, x_6, x_7, x_8, x_9]顶点分割操作的原理与此类似但偏移向量指的是行
w vertsplit(x,[0,3,5])
print(w[0], w[1])[[x_0, x_5], [x_1, x_6], [x_2, x_7]]
[[x_3, x_8], [x_4, x_9]]请注意对于上述垂直分割始终可以使用切片元素访问来代替水平和垂直分割
w [x[0:3,:], x[3:5,:]]
print(w[0], w[1])[[x_0, x_5], [x_1, x_6], [x_2, x_7]]
[[x_3, x_8], [x_4, x_9]]对于 SX 图形这种替代方法完全等效但对于 MX 图形当需要使用所有分割表达式时使用 horzsplit/vertsplit 的效率要高得多。
内积定义为 A , B : tr ( A B ) ∑ i , j A i , j B i , j \lt A,B\gt :\operatorname{tr}(A B)\sum_{i,j}\,A_{i,j}\,B_{i,j} A,B:tr(AB)∑i,jAi,jBi,j创建方法如下
x SX.sym(x,2,2)
print(dot(x,x))(((sq(x_0)sq(x_1))sq(x_2))sq(x_3))上述许多操作也是为稀疏性类第 3.5 节定义的如 vertcat、horzsplit、转置、加法返回两个稀疏性模式的联合和乘法返回两个稀疏性模式的交集。
3.7 属性查询
您可以通过调用适当的成员函数来检查矩阵或稀疏性模式是否具有特定属性例如
y SX.sym(y,10,1)
print(y.shape)(10, 1)请注意在 MATLAB 中obj.myfcn(arg) 和 myfcn(obj, arg) 都是调用成员函数 myfcn 的有效方法。从风格的角度来看后一种可能更可取。
矩阵 A 的一些常用属性如下
最后几个查询是允许假负值返回的查询示例。A.is_constant() 为真的矩阵保证为常数但如果 A.is_constant() 为假则不能保证为非常数。我们建议您在首次使用某个函数之前先查看该函数的 API 文档。
3.8 线性代数
CasADi 支持数量有限的线性代数运算例如线性方程组的求解
A MX.sym(A,3,3)
b MX.sym(b,3)
print(solve(A,b))(A\b)3.9 微积分 - 算法微分
CasADi 最核心的功能是算法或自动微分AD。对于一个函数 f : R N → R M ; f:\mathbb{R}^{N}\to\mathbb{R}^{M}; f:RN→RM; y f ( x ) , yf(x), yf(x),
前向模式方向导数可用于计算雅可比时间矢量乘积 y ^ ∂ f ∂ x x ^ . {\hat{y}}{\frac{\partial f}{\partial x}}{\hat{x}}. y^∂x∂fx^. 同样反向模式方向导数也可用于计算雅可比转置时间矢量乘积 x ˉ ( ∂ f ∂ x ) T y ˉ . {\bar{x}}\left({\frac{\partial f}{\partial x}}\right)^{\mathrm{T}}{\bar{y}}. xˉ(∂x∂f)Tyˉ. 正向和反向模式方向导数的计算成本与计算 f(x) 成正比与 x 的维数无关。
CasADi 还能高效生成完整、稀疏的雅可比。其算法非常复杂但主要包括以下步骤 自动检测雅可比的稀疏性模式 使用图着色技术找到构建完整雅可比方程所需的一些正向和/或方向导数 用数字或符号计算方向导数 组合完整的雅可比
黑塞矩阵Hessian Matrix的计算方法是首先计算梯度然后执行与上述相同的步骤利用对称性计算梯度的雅可比矩阵。
3.9.1 语法
Jacobian 的表达式是通过语法得到的
A SX.sym(A,3,2)
x SX.sym(x,2)
print(jacobian(Ax,x))四、函数对象
CasADi 允许用户创建函数对象在 C 术语中通常称为函数。这包括由符号表达式定义的函数、ODE/DAE 积分器、QP 求解器、NLP 求解器等。
函数对象通常使用以下语法创建
f functionname(name, arguments, ..., [options])名称主要是一个显示名称会显示在错误信息或生成的 C 代码注释中。随后是一组参数参数取决于类。最后用户可以传递一个选项结构用于自定义类的行为。选项结构在 Python 中是一个字典类型在 MATLAB 中是一个结构在 C 中是 CasADi 的 Dict 类型。
一个函数可以通过传递一个输入表达式列表和一个输出表达式列表来构建
x SX.sym(x,2)
y SX.sym(y)
f Function(f,[x,y],\[x,sin(y)*x])
print(f)f:(i0[2],i1)-(o0[2],o1[2]) SXFunction它定义了一个函数 f : R 2 × R → R 2 × R 2 , ( x , y ) ↦ ( x , sin ( y ) x ) . f:\mathbb{R}^{2}\times\mathbb{R}\to\mathbb{R}^{2}\times\mathbb{R}^{2},\quad(x,y)\mapsto(x,\sin(y)x). f:R2×R→R2×R2,(x,y)↦(x,sin(y)x). 请注意CasADi 中的所有函数对象包括上述对象都是多矩阵值输入、多矩阵值输出。
MX 表达式图的工作方式相同
x MX.sym(x,2)
y MX.sym(y)
f Function(f,[x,y],\[x,sin(y)*x])
print(f)f:(i0[2],i1)-(o0[2],o1[2]) MXFunction在使用类似表达式创建函数时建议将输入和输出命名如下
f Function(f,[x,y],\[x,sin(y)*x],\[x,y],[r,q])
print(f)f:(x[2],y)-(r[2],q[2]) MXFunction命名输入和输出是首选原因有很多 无需记住参数的数量或顺序 可以不设置不存在的输入或输出 语法更易读不易出错。
对于不是直接从表达式创建的函数实例稍后会遇到输入和输出会自动命名。
4.1 调用函数对象
MX 表达式可能包含对函数派生函数的调用。调用函数对象既可以进行数值计算也可以通过传递符号参数将对函数对象的调用嵌入到表达式图中参见第 4.4 节。
要调用一个函数对象必须按正确的顺序传递参数
r0, q0 f(1.1,3.3)
print(r0:,r0)
print(q0:,q0)r0: [1.1, 1.1]
q0: [-0.17352, -0.17352]或以下参数及其名称这将产生一个字典Python 中为 dictMATLAB 中为 structC 中为 std::mapstd::string,MatrixType
res f(x1.1, y3.3)
print(res:, res)res: {q: DM([-0.17352, -0.17352]), r: DM([1.1, 1.1])}调用函数对象时评估参数的维数但不一定是稀疏性模式必须与函数输入的维数相匹配但有两个例外 行向量可以代替列向量反之亦然。 无论输入维度如何标量参数总是可以传递的。这意味着将输入矩阵的所有元素都设置为该值。
当函数对象的输入数量较大或不断变化时除上述语法外另一种语法是使用调用函数该函数接收 Python list / MATLAB 单元数组或者 Python dict / MATLAB struct。返回值将具有相同的类型
arg [1.1,3.3]
res f.call(arg)
print(res:, res)
arg {x:1.1,y:3.3}
res f.call(arg)
print(res:, res)res: [DM([1.1, 1.1]), DM([-0.17352, -0.17352])]
res: {q: DM([-0.17352, -0.17352]), r: DM([1.1, 1.1])}4.2 将 MX 转换为 SX
如果 MX 图形定义的函数对象只包含内置运算如加法、平方根、矩阵乘法等元素运算以及对 SX 函数的调用则可以使用语法将其转换为纯粹由 SX 图形定义的函数
sx_function mx_function.expand()这可能会大大加快计算速度但也可能造成额外的内存开销。
4.3 非线性求根问题
请考虑下面的方程组 g 0 ( z , x 1 , x 2 , … , x n ) 0 g 1 ( z , x 1 , x 2 , … , x n ) y 1 g 2 ( z , x 1 , x 2 , … , x n ) y 2 ⋮ g m ( z , x 1 , x 2 , … , x n ) y m , \begin{array}{r l}{g_{0}(z,x_{1},x_{2},\ldots,x_{n})}{{}0}\\ {g_{1}(z,x_{1},x_{2},\ldots,x_{n})}{{}y_{1}}\\ {g_{2}(z,x_{1},x_{2},\ldots,x_{n})}{{}y_{2}}\\ {\vdots}\\{{} g_{m}(z,x_{1},x_{2},\ldots,x_{n})}{{}y_{m},}\end{array} g0(z,x1,x2,…,xn)g1(z,x1,x2,…,xn)g2(z,x1,x2,…,xn)⋮gm(z,x1,x2,…,xn)0y1y2ym,
其中第一个方程通过隐函数定理唯一定义了 z 是 x1、ldots、xn 的函数其余方程定义了辅助输出 y1、ldots、ym。
给定一个用于评估 g0、ldots、gm 的函数 g我们可以使用 CasADi 自动生成函数 G : { z g u e s s , x 1 , x 2 , ⋅ ⋅ , x n } → { z , y 1 , y 2 , ⋅ ⋅ , y m } . { G}:\{z_{\mathrm{guess}},x_{1},x_{2},\cdot\cdot,x_{n}\}\to\{z,y_{1},y_{2},\cdot\cdot,y_{m}\}. G:{zguess,x1,x2,⋅⋅,xn}→{z,y1,y2,⋅⋅,ym}. 该函数包括对 z 的猜测以处理解非唯一的情况。其语法假设 n m 1 为简单起见为
z SX.sym(x,nz)
x SX.sym(x,nx)
g0 sin(xz)
g1 cos(x-z)
g Function(g,[z,x],[g0,g1])
G rootfinder(G,newton,g)
print(G)G:(i0,i1)-(o0,o1) Newton其中寻根函数需要显示名称、求解器插件名称此处为简单的全步牛顿法和残差函数。
CasADi 中的寻根对象是微分对象导数可以精确计算到任意阶。
参见
寻根器的 API
4.4 初值问题和灵敏度分析
CasADi 可用于解决 ODE 或 DAE 中的初值问题。所使用的问题表述是带有二次函数的半显式 DAE x ˙ f o d e ( t , x , z , p , u ) , x ( 0 ) x 0 0 f a l g ( t , x , z , p , u ) q ˙ f q u a d ( t , x , z , p , u ) , q ( 0 ) 0 \begin{array}{l l}{{\dot{x}f_{\mathrm{ode}}(t,x,z,p,u),}}{{x(0)x_{0}}}\\ {{0f_{\mathrm{alg}}(t,x,z,p,u)}}{{}}\\ {{\dot{q}f_{\mathrm{quad}}(t,x,z,p,u),}}{{q(0)0}}\end{array} x˙fode(t,x,z,p,u),0falg(t,x,z,p,u)q˙fquad(t,x,z,p,u),x(0)x0q(0)0
对于常微分方程的求解器来说第二个方程和代数变量 z 必须不存在。
CasADi 中的积分器是一个函数它接受初始时间 x0 的状态、一组参数 p 和控制 u 以及代数变量仅适用于 DAEz0 的猜测并返回若干输出时间的状态向量 xf、代数变量 zf 和正交状态 qf。控制向量 u 被假定为片断常数其网格离散度与输出网格相同。
免费提供的 SUNDIALS 套件与 CasADi 一起发布包含两个常用积分器 CVodes 和 IDAS分别用于 ODE 和 DAE。这些积分器支持前向和旁向灵敏度分析通过 CasADi 的 Sundials 接口使用时CasADi 将自动计算雅各布信息这是 CVodes 和 IDAS 使用的反向微分公式 (BDF) 所需要的。此外还将自动计算前向和邻接灵敏度方程。
4.4.1 创建积分器
积分器使用 CasADi 的积分器功能创建。不同的积分器方案和接口以插件的形式实现基本上是在运行时加载的共享库。
以 DAE 为例 x ˙ z p , 0 z cos ( z ) − x \begin{array}{c}{{\dot{x} zp,}}\\ {{0 z\,\cos(z)-x}}\end{array} x˙zp,0zcos(z)−x
使用 idas 插件的积分器可以使用以下语法创建
x SX.sym(x); z SX.sym(z); p SX.sym(p)
dae {x:x, z:z, p:p, ode:zp, alg:z*cos(z)-x}
F integrator(F, idas, dae)
print(F)F:(x0,z0,p,u[0],adj_xf[],adj_zf[],adj_qf[])-(xf,zf,qf[0],adj_x0[],adj_z0[],adj_p[],adj_u[]) IdasInterfacex SX.sym(x); z SX.sym(z); p SX.sym(p);
dae struct(x,x,z,z,p,p,ode,zp,alg,z*cos(z)-x);
F integrator(F, idas, dae);
disp(F)F:(x0,z0,p,u[0],adj_xf[],adj_zf[],adj_qf[])-(xf,zf,qf[0],adj_x0[],adj_z0[],adj_p[],adj_u[]) IdasInterface这将导致从 t 0 0 t_{0}0 t00 到 t f 1 t_{f}1 tf1 的积分器即单一输出时间。我们可以使用初始条件 x ( 0 ) 0 x(0)0 x(0)0、参数 p 0.1 p0.1 p0.1 和初始时间代数变量的猜测值对函数对象进行评估 z ( 0 ) 0 z(0)0 z(0)0 如下所示
r F(x00, z00, p0.1)
print(r[xf])0.1724请注意时间跨度始终是固定的。可以通过在构造函数的 DAE 后添加两个参数将其改为默认值 t 0 0 t_{0}0 t00 和 t f 1 t_{f}1 tf1。 t f t_{f} tf 可以是一个单独的值也可以是一个值向量。它可以包括初始时间。
4.4.2 灵敏度分析
从使用角度来看积分器的行为就像本章前面通过表达式创建的函数对象一样。您可以使用函数类中的成员函数生成与方向导数正向或反向模式或完整雅各布因子相对应的新函数对象。然后对这些函数对象进行数值评估以获取灵敏度信息。文档中的示例 “sensitivity_analysis”可在 CasADi 的 Python、MATLAB 和 C 示例集中找到演示了如何使用 CasADi 计算简单 DAE 的一阶和二阶导数信息正向过正向、正向过相邻、相邻过相邻。
4.5 非线性规划 注释 本节假定读者已熟悉上述大部分内容。第 9 章中还有一个更高级别的界面。该界面可以单独学习。 与 CasADi 一起发布或连接到 CasADi 的 NLP 求解器可以求解以下形式的参数化 NLP minimize: x f ( x , p ) subject to: x l b ≤ x ≤ x u b subject to: g l b ≤ g ( x , p ) ≤ g u b \begin{array}{r l r l}{\operatorname*{minimize:}}x{{f(x,p)}}\\ {\operatorname*{subject\;to:}}{{x_{\mathrm{lb}}\leq}}{{x}}{{\leq x_{\mathrm{ub}}}}\\ {\operatorname*{subject\;to:}}{{g_{\mathrm{lb}}\leq}}{{g(x,p)}}{{\leq g_{\mathrm{ub}}}}\end{array} minimize:subjectto:subjectto:xxlb≤glb≤f(x,p)xg(x,p)≤xub≤gub
其中 x ∈ R n x x\in\mathbb{R}^{nx} x∈Rnx 是决策变量而 p ∈ R n p p\in\mathbb{R}^{np} p∈Rnp 是一个已知参数向量。
CasADi 中的 NLP 求解器是一个函数它接受参数值 §、边界 (lbx、ubx、lbg、ubg) 和对原始二元解 (x0、lam_x0、lam_g0) 的猜测并返回最优解。与积分器对象不同NLP 求解器函数目前在 CasADi 中不是可微分函数。
与 CasADi 接口的 NLP 求解器有多个。最流行的是 IPOPT它是一种开源的原始二元内点法已包含在 CasADi 安装中。其他需要安装第三方软件的 NLP 求解器包括 SNOPT、WORHP 和 KNITRO。无论使用哪种 NLP 求解器界面都会自动生成求解 NLP 所需的信息这些信息可能取决于求解器和选项。通常情况下NLP 求解器需要一个函数来给出约束函数的雅各布矩阵和拉格朗日函数的赫塞斯矩阵 L ( x , λ ) f ( x ) λ T g ( x ) L(x,\lambda)f(x)\lambda^{\mathrm{T}}\,g(x) L(x,λ)f(x)λTg(x) 与 x 有关。
4.5.1 创建 NLP 求解器
NLP 解算器使用 CasADi 的 nlpsol 函数创建。不同的求解器和接口作为插件实现。请考虑以下形式的所谓罗森布洛克Rosenbrock问题 minimize : x 2 100 z 2 x , y , z subject to: z ( 1 − x ) 2 − y 0 \begin{array}{r l}{\operatorname*{minimize}:~~~~~~~x^{2}100z^{2}}\\ {x,y,z}\\ {\operatorname*{subject \ to:}~~~z(1-x)^{2}-y0}\end{array} minimize: x2100z2x,y,zsubject to: z(1−x)2−y0
使用 ipopt 插件可以用以下语法创建该问题的求解器
x SX.sym(x); y SX.sym(y); z SX.sym(z)
nlp {x:vertcat(x,y,z), f:x**2100*z**2, g:z(1-x)**2-y}
S nlpsol(S, ipopt, nlp)
print(S)S:(x0[3],p[],lbx[3],ubx[3],lbg,ubg,lam_x0[3],lam_g0)-(x[3],f,g,lam_x[3],lam_g,lam_p[]) IpoptInterfacex SX.sym(x); y SX.sym(y); z SX.sym(z);
nlp struct(x,[x;y;z], f,x^2100*z^2, g,z(1-x)^2-y);
S nlpsol(S, ipopt, nlp);
disp(S)S:(x0[3],p[],lbx[3],ubx[3],lbg,ubg,lam_x0[3],lam_g0)-(x[3],f,g,lam_x[3],lam_g,lam_p[]) IpoptInterface创建求解器后我们可以使用 [2.5,3.0,0.75] 作为初始猜测通过评估函数 S 来求解 NLP
r S(x0[2.5,3.0,0.75],\lbg0, ubg0)
x_opt r[x]
print(x_opt: , x_opt)This is Ipopt version 3.14.11, running with linear solver MUMPS 5.4.1.Number of nonzeros in equality constraint Jacobian...: 3
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 2Total number of variables............................: 3variables with only lower bounds: 0variables with lower and upper bounds: 0variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0inequality constraints with only lower bounds: 0inequality constraints with lower and upper bounds: 0inequality constraints with only upper bounds: 0iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls0 6.2500000e01 0.00e00 9.00e01 -1.0 0.00e00 - 0.00e00 0.00e00 01 1.8457621e01 1.48e-02 4.10e01 -1.0 4.10e-01 2.0 1.00e00 1.00e00f 12 7.8031530e00 3.84e-03 8.76e00 -1.0 2.63e-01 1.5 1.00e00 1.00e00f 13 7.1678278e00 9.42e-05 1.04e00 -1.0 9.32e-02 1.0 1.00e00 1.00e00f 14 6.7419924e00 6.18e-03 9.95e-01 -1.0 2.69e-01 0.6 1.00e00 1.00e00f 15 5.4363330e00 7.03e-02 1.04e00 -1.7 8.40e-01 0.1 1.00e00 1.00e00f 16 1.2144815e00 1.52e00 1.32e00 -1.7 3.21e00 -0.4 1.00e00 1.00e00f 17 1.0214718e00 2.51e-01 1.17e01 -1.7 1.33e00 0.9 1.00e00 1.00e00h 18 3.1864085e-01 1.04e-03 7.53e-01 -1.7 3.58e-01 - 1.00e00 1.00e00f 19 3.7092062e-66 3.19e-01 2.57e-32 -1.7 5.64e-01 - 1.00e00 1.00e00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls10 0.0000000e00 0.00e00 0.00e00 -1.7 3.19e-01 - 1.00e00 1.00e00h 1Number of Iterations....: 10(scaled) (unscaled)
Objective...............: 0.0000000000000000e00 0.0000000000000000e00
Dual infeasibility......: 0.0000000000000000e00 0.0000000000000000e00
Constraint violation....: 0.0000000000000000e00 0.0000000000000000e00
Variable bound violation: 0.0000000000000000e00 0.0000000000000000e00
Complementarity.........: 0.0000000000000000e00 0.0000000000000000e00
Overall NLP error.......: 0.0000000000000000e00 0.0000000000000000e00Number of objective function evaluations 11
Number of objective gradient evaluations 11
Number of equality constraint evaluations 11
Number of inequality constraint evaluations 0
Number of equality constraint Jacobian evaluations 11
Number of inequality constraint Jacobian evaluations 0
Number of Lagrangian Hessian evaluations 10
Total seconds in IPOPT 0.007EXIT: Optimal Solution Found.S : t_proc (avg) t_wall (avg) n_evalnlp_f | 59.00us ( 5.36us) 14.26us ( 1.30us) 11nlp_g | 132.00us ( 12.00us) 29.65us ( 2.70us) 11nlp_grad_f | 80.00us ( 6.67us) 19.74us ( 1.64us) 12nlp_hess_l | 60.00us ( 6.00us) 13.88us ( 1.39us) 10nlp_jac_g | 67.00us ( 5.58us) 16.22us ( 1.35us) 12total | 31.09ms ( 31.09ms) 7.77ms ( 7.77ms) 1
x_opt: [0, 1, 0]4.6 二次规划
CasADi 提供求解二次方程程序 (QP) 的接口。支持的求解器有开源求解器 qpOASES与 CasADi 一起发布、OOQP、OSQP 和 PROXQP以及商用求解器 CPLEX 和 GUROBI。
在 CasADi 中求解 QP 有两种不同的方法即使用高级接口和低级接口。下文将对这两种方法进行介绍。
4.6.1 高级接口
二次规划的高级接口与非线性规划的接口相似即预期问题的形式为 (4.5.1)限制条件是目标函数 f ( x , p ) f(x,p) f(x,p) 必须是 x 的凸二次函数约束函数 g ( x , p ) g(x,p) g(x,p) 必须与 x . 如果函数分别不是二次函数和线性函数则在当前线性化点求解该点由 x 的 初始猜测 给出 x .
如果目标函数不是凸函数求解器可能找不到解也可能找不到解或者解不是唯一的。
为说明语法我们考虑下面的凸 QP minimize: x 2 y 2 x , y subject to: x y − 10 ≥ 0 \begin{array}{r l}{\operatorname*{minimize:}}{{}x^{2}y^{2}}\\ {x,y} \\ {\operatorname*{subject \ to:}}{{}xy-10\geq0}\end{array} minimize:x,ysubject to:x2y2xy−10≥0
要通过高级界面解决这个问题我们只需用 qpsol 代替 nlpsol并使用 QP 求解器插件如 CasADi 分布式 qpOASES
x SX.sym(x); y SX.sym(y)
qp {x:vertcat(x,y), f:x**2y**2, g:xy-10}
S qpsol(S, qpoases, qp)
print(S)S:(x0[2],p[],lbx[2],ubx[2],lbg,ubg,lam_x0[2],lam_g0)-(x[2],f,g,lam_x[2],lam_g,lam_p[]) MXFunctionx SX.sym(x); y SX.sym(y);
qp struct(x,[x;y], f,x^2y^2, g,xy-10);
S qpsol(S, qpoases, qp);
disp(S)S:(x0[2],p[],lbx[2],ubx[2],lbg,ubg,lam_x0[2],lam_g0)-(x[2],f,g,lam_x[2],lam_g,lam_p[]) MXFunction创建的求解器对象 S 将与使用 nlpsol 创建的求解器对象具有相同的输入和输出签名。由于解是唯一的因此提供初始猜测并不那么重要
r S(lbg0)
x_opt r[x]
print(x_opt: , x_opt)#################### qpOASES -- QP NO. 1 #####################Iter | StepLength | Info | nFX | nAC ---------------------------------------------------------------- 0 | 1.866661e-07 | ADD CON 0 | 1 | 1 1 | 8.333622e-10 | REM BND 1 | 0 | 1 2 | 1.000000e00 | QP SOLVED | 0 | 1
x_opt: [5, 5]4.6.2 低级接口
而底层接口则解决以下形式的 QPs minimize: x 1 2 x T H x g T x subject to: x l b ≤ x ≤ x u b subject to: a l b ≤ A x ≤ a u b \begin{array}{r l r l}{\operatorname*{minimize:}}x{\textstyle{\frac{1}{2}}x^{T}\,H\,xg^{T}\,x}\\ {\operatorname*{subject\;to:}}{{x_{\mathrm{lb}}\leq}}{{x}}{{\leq x_{\mathrm{ub}}}}\\ {\operatorname*{subject\;to:}}{{a_{\mathrm{lb}}\leq}}{Ax}{{\leq a_{\mathrm{ub}}}}\end{array} minimize:subjectto:subjectto:xxlb≤alb≤21xTHxgTxxAx≤xub≤aub
以这种形式编码问题 (4.6.1)省略了无穷大的界限非常简单
H 2*DM.eye(2)
A DM.ones(1,2)
g DM.zeros(2)
lba 10.H 2*DM.eye(2);
A DM.ones(1,2);
g DM.zeros(2);
lba 10;在创建求解器实例时我们不再传递 QP 的符号表达式而是传递矩阵的稀疏性模式 H 和 A . 由于我们使用了 CasADi 的 DM 类型因此只需查询稀疏性模式即可
qp {}
qp[h] H.sparsity()
qp[a] A.sparsity()
S conic(S,qpoases,qp)
print(S)S:(h[2x2,2nz],g[2],a[1x2],lba,uba,lbx[2],ubx[2],x0[2],lam_x0[2],lam_a0,q[],p[])-(x[2],cost,lam_a,lam_x[2]) QpoasesInterfaceqp struct;
qp.h H.sparsity();
qp.a A.sparsity();
S conic(S,qpoases,qp);
disp(S)S:(h[2x2,2nz],g[2],a[1x2],lba,uba,lbx[2],ubx[2],x0[2],lam_x0[2],lam_a0,q[],p[])-(x[2],cost,lam_a,lam_x[2]) QpoasesInterface与高层接口相比返回的函数实例将具有不同的输入/输出签名其中包括矩阵 H 和 A :
r S(hH, gg, \aA, lbalba)
x_opt r[x]
print(x_opt: , x_opt)#################### qpOASES -- QP NO. 1 #####################Iter | StepLength | Info | nFX | nAC ---------------------------------------------------------------- 0 | 1.866661e-07 | ADD CON 0 | 1 | 1 1 | 8.333622e-10 | REM BND 1 | 0 | 1 2 | 1.000000e00 | QP SOLVED | 0 | 1
x_opt: [5, 5]
