CasADi - 最优控制开源 Python/MATLAB 库

系列文章目录

---

---

前言

---

一、介绍

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 模块已成功通过 Octave（4.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 storage，CCS）–来存储矩阵。下面，我们将介绍这一框架的最基本类别。

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((10+sq(x)))

您也可以在不使用任何符号基元的情况下创建常量 SX 实例：

• B1 = SX.zeros(4,5)： 一个全为零的 4 乘 5 的密集空矩阵

• B2 = SX(4,5)： 一个全部为零的稀疏 4×5 空矩阵

• B4 = SX.eye(4)： 对角线上有 1 的稀疏 4×4 矩阵

B1: @1=0, 
[[@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: @1=1, 
[[@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 矩阵。

• (仅限 Python）SX(list)： 创建列向量 (n-1例如 SX([1,2,3,4])（注意 Python 列表与 MATLAB/Octave 水平连接之间的区别，两者都使用方括号语法）。

• (仅限 Python）SX(列表的列表)： 用列表中的元素创建密集矩阵，例如 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 *。

遗憾的是，Octave（4.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) # equivalent

SX 的更多使用示例可在 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)

@1=3, 
[[((@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 不同的是，这些基本运算并不局限于标量一元运算或二元运算 ($\mathbb{R}\to \mathbb{R}$ 或 $\mathbb{R}\times \mathbb{R}\to \mathbb{R}$). 相反，用于形成 MX 表达式的基本运算可以是一般的多稀疏矩阵值输入、多稀疏矩阵值输出函数：${\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}.$

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_0+x_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]]
@1=1, 
[[@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, y@y)

[[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))

@1=1, 
[[@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 个较小的表达式，除了要拆分的表达式外，还需要提供一个长度为 n+1 的偏移向量。偏移向量的第一个元素必须是 0，最后一个元素必须是列数。其余元素必须按不递减的顺序排列。拆分操作的输出 i 将包含偏移量[i]≤c<偏移量[i+1]的列 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>:=\mathrm{tr}(AB)={\sum }_{i,j}{A}_{i,j}{B}_{i,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:{\mathbb{R}}^N\to {\mathbb{R}}^M;$
$$y=f(x),$$

前向模式方向导数可用于计算雅可比时间矢量乘积：
$$\hat{y}=\frac{\partial f}{\partial x}\hat{x}.$$
同样，反向模式方向导数也可用于计算雅可比转置时间矢量乘积：
$$\bar{x}={\left(\frac{\partial f}{\partial x}\right)}^{\mathrm{T}}\bar{y}.$$
正向和反向模式方向导数的计算成本与计算 f(x) 成正比，与 x 的维数无关。

CasADi 还能高效生成完整、稀疏的雅可比。其算法非常复杂，但主要包括以下步骤：

• 自动检测雅可比的稀疏性模式

• 使用图着色技术找到构建完整雅可比方程所需的一些正向和/或方向导数

• 用数字或符号计算方向导数

• 组合完整的雅可比

黑塞矩阵（Hessian Matrix）的计算方法是，首先计算梯度，然后执行与上述相同的步骤，利用对称性计算梯度的雅可比矩阵。

3.9.1 语法

Jacobian 的表达式是通过语法得到的：

A = SX.sym('A',3,2)
x = SX.sym('x',2)
print(jacobian(A@x,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:{\mathbb{R}}^2\times \mathbb{R}\to {\mathbb{R}}^2\times {\mathbb{R}}^2,(x,y)\mapsto (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 中为 dict，MATLAB 中为 struct，C++ 中为 std::mapstd::string,MatrixType）：

res = f(x=1.1, y=3.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 非线性求根问题

请考虑下面的方程组：
$$\begin{array}{rl}{g}_0(z,x_1,x_2,\dots ,x_n)& =0\\ g_1(z,x_1,x_2,\dots ,x_n)& =y_1\\ g_2(z,x_1,x_2,\dots ,x_n)& =y_2\\ ⋮& \\ g_m(z,x_1,x_2,\dots ,x_n)& =y_m,\end{array}$$

其中第一个方程通过隐函数定理唯一定义了 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(x+z)
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：
$$\begin{array}{ll}\dot{x}=f_{\mathrm{ode}}(t,x,z,p,u),& x(0)=x_0\\ 0=f_{\mathrm{alg}}(t,x,z,p,u)& \\ \dot{q}=f_{\mathrm{quad}}(t,x,z,p,u),& q(0)=0\end{array}$$

对于常微分方程的求解器来说，第二个方程和代数变量 z 必须不存在。

CasADi 中的积分器是一个函数，它接受初始时间 x0 的状态、一组参数 p 和控制 u 以及代数变量（仅适用于 DAE）z0 的猜测，并返回若干输出时间的状态向量 xf、代数变量 zf 和正交状态 qf。控制向量 u 被假定为片断常数，其网格离散度与输出网格相同。

免费提供的 SUNDIALS 套件（与 CasADi 一起发布）包含两个常用积分器 CVodes 和 IDAS，分别用于 ODE 和 DAE。这些积分器支持前向和旁向灵敏度分析，通过 CasADi 的 Sundials 接口使用时，CasADi 将自动计算雅各布信息，这是 CVodes 和 IDAS 使用的反向微分公式 (BDF) 所需要的。此外，还将自动计算前向和邻接灵敏度方程。

4.4.1 创建积分器

积分器使用 CasADi 的积分器功能创建。不同的积分器方案和接口以插件的形式实现，基本上是在运行时加载的共享库。

以 DAE 为例：
$$\begin{array}{c}\dot{x}=z+p,\\ 0=z\cos (z)-x\end{array}$$

使用 "idas "插件的积分器可以使用以下语法创建：

x = SX.sym('x'); z = SX.sym('z'); p = SX.sym('p')
dae = {'x':x, 'z':z, 'p':p, 'ode':z+p, '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[]) IdasInterface

x = SX.sym('x'); z = SX.sym('z'); p = SX.sym('p');
dae = struct('x',x,'z',z,'p',p,'ode',z+p,'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_f=1$ 的积分器，即单一输出时间。我们可以使用初始条件 $x(0)=0$、参数 $p=0.1$ 和初始时间代数变量的猜测值对函数对象进行评估 $z(0)=0$ 如下所示

r = F(x0=0, z0=0, p=0.1)
print(r['xf'])

0.1724

请注意，时间跨度始终是固定的。可以通过在构造函数的 DAE 后添加两个参数，将其改为默认值 $t_0=0$ 和 $t_f=1$。$t_f$ 可以是一个单独的值，也可以是一个值向量。它可以包括初始时间。

4.4.2 灵敏度分析

从使用角度来看，积分器的行为就像本章前面通过表达式创建的函数对象一样。您可以使用函数类中的成员函数生成与方向导数（正向或反向模式）或完整雅各布因子相对应的新函数对象。然后对这些函数对象进行数值评估，以获取灵敏度信息。文档中的示例 “sensitivity_analysis”（可在 CasADi 的 Python、MATLAB 和 C++ 示例集中找到）演示了如何使用 CasADi 计算简单 DAE 的一阶和二阶导数信息（正向过正向、正向过相邻、相邻过相邻）。

4.5 非线性规划

注释

本节假定读者已熟悉上述大部分内容。第 9 章中还有一个更高级别的界面。该界面可以单独学习。

与 CasADi 一起发布或连接到 CasADi 的 NLP 求解器可以求解以下形式的参数化 NLP：
$$\begin{array}{rlrl}minimize:& x& f(x,p)& \\ \mathrm{subject}to:& x_{\mathrm{lb}}\le & x& \le x_{\mathrm{ub}}\\ \mathrm{subject}to:& g_{\mathrm{lb}}\le & g(x,p)& \le g_{\mathrm{ub}}\end{array}$$

其中，$x\in {\mathbb{R}}^{nx}$ 是决策变量，而 $p\in {\mathbb{R}}^{np}$ 是一个已知参数向量。

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,\lambda )=f(x)+{\lambda }^{\mathrm{T}}g(x)$ 与 x 有关）。

4.5.1 创建 NLP 求解器

NLP 解算器使用 CasADi 的 nlpsol 函数创建。不同的求解器和接口作为插件实现。请考虑以下形式的所谓罗森布洛克（Rosenbrock）问题：
$$\begin{array}{r}\mathrm{minimize}:x^2+100z^2\\ x,y,z\\ subject\; to:z+(1-x{)}^2-y=0\end{array}$$

使用 "ipopt "插件，可以用以下语法创建该问题的求解器：

x = SX.sym('x'); y = SX.sym('y'); z = SX.sym('z')
nlp = {'x':vertcat(x,y,z), 'f':x**2+100*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[]) IpoptInterface

x = SX.sym('x'); y = SX.sym('y'); z = SX.sym('z');
nlp = struct('x',[x;y;z], 'f',x^2+100*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],\lbg=0, ubg=0)
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.2500000e+01 0.00e+00 9.00e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   01  1.8457621e+01 1.48e-02 4.10e+01  -1.0 4.10e-01   2.0 1.00e+00 1.00e+00f  12  7.8031530e+00 3.84e-03 8.76e+00  -1.0 2.63e-01   1.5 1.00e+00 1.00e+00f  13  7.1678278e+00 9.42e-05 1.04e+00  -1.0 9.32e-02   1.0 1.00e+00 1.00e+00f  14  6.7419924e+00 6.18e-03 9.95e-01  -1.0 2.69e-01   0.6 1.00e+00 1.00e+00f  15  5.4363330e+00 7.03e-02 1.04e+00  -1.7 8.40e-01   0.1 1.00e+00 1.00e+00f  16  1.2144815e+00 1.52e+00 1.32e+00  -1.7 3.21e+00  -0.4 1.00e+00 1.00e+00f  17  1.0214718e+00 2.51e-01 1.17e+01  -1.7 1.33e+00   0.9 1.00e+00 1.00e+00h  18  3.1864085e-01 1.04e-03 7.53e-01  -1.7 3.58e-01    -  1.00e+00 1.00e+00f  19  3.7092062e-66 3.19e-01 2.57e-32  -1.7 5.64e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls10  0.0000000e+00 0.00e+00 0.00e+00  -1.7 3.19e-01    -  1.00e+00 1.00e+00h  1Number of Iterations....: 10(scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00Number 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)$ 必须是 x 的凸二次函数，约束函数 $g(x,p)$ 必须与 x . 如果函数分别不是二次函数和线性函数，则在当前线性化点求解，该点由 x 的 "初始猜测 "给出 x .

如果目标函数不是凸函数，求解器可能找不到解，也可能找不到解，或者解不是唯一的。

为说明语法，我们考虑下面的凸 QP：

$$\begin{array}{rl}minimize:& x^2+y^2\\ x,y& \\ subject\; to:& x+y-10\ge 0\end{array}$$

要通过高级界面解决这个问题，我们只需用 qpsol 代替 nlpsol，并使用 QP 求解器插件，如 CasADi 分布式 qpOASES：

x = SX.sym('x'); y = SX.sym('y')
qp = {'x':vertcat(x,y), 'f':x**2+y**2, 'g':x+y-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[]) MXFunction

x = SX.sym('x'); y = SX.sym('y');
qp = struct('x',[x;y], 'f',x^2+y^2, 'g',x+y-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(lbg=0)
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.000000e+00   |    QP SOLVED     |     0   |     1   
x_opt:  [5, 5]

4.6.2 低级接口

而底层接口则解决以下形式的 QPs：

$$\begin{array}{rlrl}minimize:& x& \frac{1}{2}{x}^{T}Hx+g^{T}x& \\ \mathrm{subject}to:& x_{\mathrm{lb}}\le & x& \le x_{\mathrm{ub}}\\ \mathrm{subject}to:& a_{\mathrm{lb}}\le & Ax& \le a_{\mathrm{ub}}\end{array}$$

以这种形式编码问题 (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]) QpoasesInterface

qp = 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(h=H, g=g, \a=A, lba=lba)
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.000000e+00   |    QP SOLVED     |     0   |     1   
x_opt:  [5, 5]