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 Horn
            and James B Rawlings and Moritz Diehl},
  Title = {{CasADi} -- {A} software framework for nonlinear optimization
           and 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) # equivalent

C_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.............:        2

Total number of variables............................:        3
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  6.2500000e+01 0.00e+00 9.00e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.8457621e+01 1.48e-02 4.10e+01  -1.0 4.10e-01   2.0 1.00e+00 1.00e+00f  1
   2  7.8031530e+00 3.84e-03 8.76e+00  -1.0 2.63e-01   1.5 1.00e+00 1.00e+00f  1
   3  7.1678278e+00 9.42e-05 1.04e+00  -1.0 9.32e-02   1.0 1.00e+00 1.00e+00f  1
   4  6.7419924e+00 6.18e-03 9.95e-01  -1.0 2.69e-01   0.6 1.00e+00 1.00e+00f  1
   5  5.4363330e+00 7.03e-02 1.04e+00  -1.7 8.40e-01   0.1 1.00e+00 1.00e+00f  1
   6  1.2144815e+00 1.52e+00 1.32e+00  -1.7 3.21e+00  -0.4 1.00e+00 1.00e+00f  1
   7  1.0214718e+00 2.51e-01 1.17e+01  -1.7 1.33e+00   0.9 1.00e+00 1.00e+00h  1
   8  3.1864085e-01 1.04e-03 7.53e-01  -1.7 3.58e-01    -  1.00e+00 1.00e+00f  1
   9  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  ls
  10  0.0000000e+00 0.00e+00 0.00e+00  -1.7 3.19e-01    -  1.00e+00 1.00e+00h  1

Number 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+00


Number 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.007

EXIT: Optimal Solution Found.
           S  :   t_proc      (avg)   t_wall      (avg)    n_eval
       nlp_f  |  59.00us (  5.36us)  14.26us (  1.30us)        11
       nlp_g  | 132.00us ( 12.00us)  29.65us (  2.70us)        11
  nlp_grad_f  |  80.00us (  6.67us)  19.74us (  1.64us)        12
  nlp_hess_l  |  60.00us (  6.00us)  13.88us (  1.39us)        10
   nlp_jac_g  |  67.00us (  5.58us)  16.22us (  1.35us)        12
       total  |  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]

4.7 for 循环等价

在 CasADi 中使用表达式图建模时，通常会使用宿主语言（C++/Python/Matlab）的 for 循环结构。

图的大小将随循环的大小线性增长 n 的线性增长，表达式图的构建时间和使用该表达式的函数的初始化时间也将随之增长。

我们提供了一些特殊的结构来改善这种复杂性。

4.7.1 映射

假设您想在矩阵 $X\in {\mathbb{R}}^{n\times N}$ 的所有列上重复计算函数 $f:{\mathbb{R}}^n\to {\mathbb{R}}^m$，并将所有结果汇总到结果矩阵 $Y\in {\mathbb{R}}^{m\times N}$ 中。

N = 4
X = MX.sym("X",1,N)

print(f)

ys = []
for i in range(N):
  ys.append(f(X[:,i]))

Y = hcat(ys)
F = Function('F',[X],[Y])
print(F)
f:(i0)->(o0) SXFunction
F:(i0[1x4])->(o0[1x4]) MXFunction

总函数 $F:{\mathbb{R}}^{n\times N}\to {\mathbb{R}}^{m\times N}$ 可以通过映射构造得到：

F = f.map(N)
print(F)
map4_f:(i0[1x4])->(o0[1x4]) Map

在对 F 进行评估时，可以指示 CasADi 进行并行处理。在下面的示例中，我们为映射任务分配了 2 个线程。

F = f.map(N,"thread",2)
print(F)
threadmap2_map2_f:(i0[1x4])->(o0[1x4]) ThreadMap

映射操作支持具有多个输入/输出（也可以是矩阵）的原始函数 f。聚合总是横向进行的。

映射操作的图形大小和初始化时间不变。

4.7.2. 折叠 Fold

如果每个 for 循环的迭代都依赖于前一次迭代的结果，则适用折叠结构。在下文中，x 变量充当累加器，其初始化为 $x_0\in {\mathbb{R}}^n$

x = x0
for i in range(N):
  x = f(x)

F = Function('F',[x0],[x])
print(F)
F:(i0)->(o0) MXFunction

对于给定函数 $f:{\mathbb{R}}^n\to {\mathbb{R}}^n$，结果函数 $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ 可以通过折叠构造得到：

F = f.fold(N)
print(F)
fold_f:(i0)->(o0) MXFunction

如果需要将中间累加器值作为输出 (${\mathbb{R}}^n\to {\mathbb{R}}^{n\times N}$)，请使用 mapaccum 代替 fold。

fold/mapaccum 运算支持带有多个输入/输出（也可以是矩阵）的原始函数 f。第一个输入和输出用于累加，其余输入在迭代过程中逐列读取。

fold/mapaccum 操作的图形大小和初始化时间与 n 成对数关系。

4.7.3. 条件运算
通过构建条件函数实例，可以在 CasADi 表达式图中对表达式进行条件评估。该函数需要一些现有的函数实例 $f_1,f_2,f_n$ 以及一个 "默认 "函数 $f_{def}$。所有这些函数必须具有相同的输入和输出特征，即具有相同维数的输入和输出：

x = SX.sym("x")
f0 = Function("f0",[x],[sin(x)])
f1 = Function("f1",[x],[cos(x)])
f2 = Function("f2",[x],[tan(x)])
f_cond = Function.conditional('f_cond', [f0, f1], f2)
print(f_cond)
f_cond:(i0,i1)->(o0) Switch
x = SX.sym('x');
f0 = Function('f0',{x},{sin(x)});
f1 = Function('f1',{x},{cos(x)});
f2 = Function('f2',{x},{tan(x)});
f_cond = Function.conditional('f_cond', {f0, f1}, f2);
disp(f_cond);
f_cond:(i0,i1)->(o0) Switch

结果是一个新的函数实例，具有相同的输入/输出特征，但多了一个与索引相对应的输入。对它的计算相当于$f_{\mathrm{cond}}(c,x_1,x_2,\dots ,x_m)=\{\begin{array}{ll}{f}_0(x_1,x_2,\dots ,x_m)& \mathrm{if}c=0,\\ f_1(x_1,x_2,\dots ,x_m)& \mathrm{if}c=1,\\ f_2(x_1,x_2,\dots ,x_m)& \mathrm{otherwise}\end{array}$
上述函数可以是缺失的（即空指针 Function()），在这种情况下，所有输出都将被求值为 NaN。请注意，计算是 "短路 "的，即只对相关函数进行计算。这也适用于任何导数计算。

一种常见的特殊情况是，除了默认情况外只有一种情况。这等同于 if-then-else 语句，可以用速记法来书写：

x = SX.sym("x")
f_true = Function("f_true",[x],[cos(x)])
f_false = Function("f_false",[x],[sin(x)])
f_cond = Function.if_else('f_cond', f_true, f_false)
print(f_cond)
f_cond:(i0,i1)->(o0) Switch
x = SX.sym('x');
f_true = Function('f_true',{x},{cos(x)});
f_false = Function('f_false',{x},{sin(x)});
f_cond = Function.if_else('f_cond', f_true, f_false);
disp(f_cond);
f_cond:(i0,i1)->(o0) Switch

请注意，如果在基于梯度的优化算法中使用此类条件表达式，可能会导致非平滑表达式无法收敛。

脚注

[1] 对于有自由结束时间的问题，可以通过引入一个额外参数来缩放时间，并用 t 代替从 0 到 1 的无量纲时间变量

五、 生成 C 代码

CasADi 中函数对象的数值计算通常在虚拟机中进行，虚拟机是 CasADi 符号框架的一部分。但 CasADi 也支持为大量函数对象子集生成独立的 C 代码。

C 代码生成之所以有趣，有以下几个原因：

加快评估时间。根据经验，使用代码优化标志编译的自动生成代码的数值评估速度比在 CasADi 虚拟机中执行的相同代码快 4 到 10 倍。

允许在未安装 CasADi 的系统（如嵌入式系统）上编译代码。编译生成的代码只需要一个 C 编译器。

调试和剖析功能。生成的代码本质上是虚拟机评估的镜像，如果某个操作速度较慢，使用 gprof 等剖析工具分析生成的代码时很可能会显示出来。通过观察代码，还可以发现哪些操作可能以次优方式完成。如果代码很长，需要很长时间才能编译，则表明某些函数需要拆分成更小的嵌套函数。

5.1. 生成代码的语法

生成的 C 代码可以像调用函数实例的生成成员函数一样简单。

x = MX.sym('x',2)
y = MX.sym('y')
f = Function('f',[x,y],\
      [x,sin(y)*x],\
      ['x','y'],['r','q'])
f.generate('gen.c')
print(open('gen.c','r').read())
/* This file was automatically generated by CasADi 3.6.4.
 *  It consists of: 
 *   1) content generated by CasADi runtime: not copyrighted
 *   2) template code copied from CasADi source: permissively licensed (MIT-0)
 *   3) user code: owned by the user
 *
 */
#ifdef __cplusplus
extern "C" {
#endif

/* How to prefix internal symbols */
#ifdef CASADI_CODEGEN_PREFIX
  #define CASADI_NAMESPACE_CONCAT(NS, ID) _CASADI_NAMESPACE_CONCAT(NS, ID)
  #define _CASADI_NAMESPACE_CONCAT(NS, ID) NS ## ID
  #define CASADI_PREFIX(ID) CASADI_NAMESPACE_CONCAT(CODEGEN_PREFIX, ID)
#else
  #define CASADI_PREFIX(ID) gen_ ## ID
#endif

#include <math.h>

#ifndef casadi_real
#define casadi_real double
#endif

#ifndef casadi_int
#define casadi_int long long int
#endif

/* Add prefix to internal symbols */
#define casadi_copy CASADI_PREFIX(copy)
#define casadi_f0 CASADI_PREFIX(f0)
#define casadi_s0 CASADI_PREFIX(s0)
#define casadi_s1 CASADI_PREFIX(s1)

/* Symbol visibility in DLLs */
#ifndef CASADI_SYMBOL_EXPORT
  #if defined(_WIN32) || defined(__WIN32__) || defined(__CYGWIN__)
    #if defined(STATIC_LINKED)
      #define CASADI_SYMBOL_EXPORT
    #else
      #define CASADI_SYMBOL_EXPORT __declspec(dllexport)
    #endif
  #elif defined(__GNUC__) && defined(GCC_HASCLASSVISIBILITY)
    #define CASADI_SYMBOL_EXPORT __attribute__ ((visibility ("default")))
  #else
    #define CASADI_SYMBOL_EXPORT
  #endif
#endif

void casadi_copy(const casadi_real* x, casadi_int n, casadi_real* y) {
  casadi_int i;
  if (y) {
    if (x) {
      for (i=0; i<n; ++i) *y++ = *x++;
    } else {
      for (i=0; i<n; ++i) *y++ = 0.;
    }
  }
}

static const casadi_int casadi_s0[6] = {2, 1, 0, 2, 0, 1};
static const casadi_int casadi_s1[5] = {1, 1, 0, 1, 0};

/* f:(x[2],y)->(r[2],q[2]) */
static int casadi_f0(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, int mem) {
  casadi_int i;
  casadi_real *rr;
  const casadi_real *cs;
  casadi_real *w0=w+0, w1;
  /* #0: @0 = input[0][0] */
  casadi_copy(arg[0], 2, w0);
  /* #1: output[0][0] = @0 */
  casadi_copy(w0, 2, res[0]);
  /* #2: @1 = input[1][0] */
  w1 = arg[1] ? arg[1][0] : 0;
  /* #3: @1 = sin(@1) */
  w1 = sin( w1 );
  /* #4: @0 = (@1*@0) */
  for (i=0, rr=w0, cs=w0; i<2; ++i) (*rr++)  = (w1*(*cs++));
  /* #5: output[1][0] = @0 */
  casadi_copy(w0, 2, res[1]);
  return 0;
}

CASADI_SYMBOL_EXPORT int f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, int mem){
  return casadi_f0(arg, res, iw, w, mem);
}

CASADI_SYMBOL_EXPORT int f_alloc_mem(void) {
  return 0;
}

CASADI_SYMBOL_EXPORT int f_init_mem(int mem) {
  return 0;
}

CASADI_SYMBOL_EXPORT void f_free_mem(int mem) {
}

CASADI_SYMBOL_EXPORT int f_checkout(void) {
  return 0;
}

CASADI_SYMBOL_EXPORT void f_release(int mem) {
}

CASADI_SYMBOL_EXPORT void f_incref(void) {
}

CASADI_SYMBOL_EXPORT void f_decref(void) {
}

CASADI_SYMBOL_EXPORT casadi_int f_n_in(void) { return 2;}

CASADI_SYMBOL_EXPORT casadi_int f_n_out(void) { return 2;}

CASADI_SYMBOL_EXPORT casadi_real f_default_in(casadi_int i) {
  switch (i) {
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const char* f_name_in(casadi_int i) {
  switch (i) {
    case 0: return "x";
    case 1: return "y";
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const char* f_name_out(casadi_int i) {
  switch (i) {
    case 0: return "r";
    case 1: return "q";
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const casadi_int* f_sparsity_in(casadi_int i) {
  switch (i) {
    case 0: return casadi_s0;
    case 1: return casadi_s1;
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const casadi_int* f_sparsity_out(casadi_int i) {
  switch (i) {
    case 0: return casadi_s0;
    case 1: return casadi_s0;
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT int f_work(casadi_int *sz_arg, casadi_int* sz_res, casadi_int *sz_iw, casadi_int *sz_w) {
  if (sz_arg) *sz_arg = 4;
  if (sz_res) *sz_res = 3;
  if (sz_iw) *sz_iw = 0;
  if (sz_w) *sz_w = 3;
  return 0;
}


#ifdef __cplusplus
} /* extern "C" */
#endif

这将创建一个 C 语言文件 gen.c，其中包含函数 f 及其所有依赖函数和所需的辅助函数。我们将在第 5.1 节再次讨论如何使用该文件，生成代码的结构将在下文第 5.3 节描述。

通过使用 CasADi 的 CodeGenerator 类，可以生成包含多个 CasADi 函数的 C 文件：

f = Function('f',[x],[sin(x)])
g = Function('g',[x],[cos(x)])
C = CodeGenerator('gen.c')
C.add(f)
C.add(g)
C.generate()

generate 函数和 CodeGenerator 构造函数都将可选的选项字典作为参数，以便自定义代码生成。两个有用的选项是：main（生成主入口）和 mex（生成 mexFunction 入口）：

f = Function('f',[x],[sin(x)])
opts = dict(main=True, \
            mex=True)
f.generate('gen.c',opts)

这样就可以分别从命令行和 MATLAB 中执行函数，详见下文第 5.2 节。

如果计划在某些 C/C++ 应用程序中直接针对生成的代码进行部件链接，一个有用的选项是 with_header，它可以控制创建一个头文件，其中包含外部链接函数的声明，即生成代码的 API，将在下文第 5.3 节中介绍。

以下是代码生成器类的可用选项列表：

5.2 使用生成的代码

生成的 C 代码可以多种方式使用：

可以将代码编译成动态链接库（DLL），然后使用 CasADi 的外部函数创建 Function 实例。用户也可以选择依靠 CasADi 即时进行编译。

生成的代码可以编译成 MEX 函数，并从 MATLAB 中执行。

生成的代码可以通过命令行执行。

用户可以静态或动态地将生成的代码部件链接到自己的 C/C++ 应用程序，访问生成代码的 C API。

代码可以编译成动态链接库，然后用户可以使用 Linux/OS X 上的 dlopen 或 Windows 上的 LoadLibrary 手动访问 C API。

下文将对此进行详细说明。

5.2.1. CasADi 的外部函数

外部命令允许用户从动态链接库中创建一个 Function 实例，其入口点由第 5.3 节所述的 C API 描述。由于自动生成的文件是独立的
[1]，在 Linux/OSX 上，编译就像下发命令一样简单：

gcc -fPIC -shared gen.c -o gen.so

或使用 MATLAB 的 system 命令或 Python 的 os.system 命令。或者，也可以使用 MATLAB 的 system 命令或 Python 的 os.system 命令。假设 gen.c 已按上一节所述创建，那么我们就可以创建一个函数 f，如下所示：

f = external('f', './gen.so')
print(f(3.14))
0.00159265

我们还可以使用 CasADi 的导入器类（Importer class），依靠 CasADi 即时进行编译。这是一个插件类，在撰写本文时有两个支持的插件，即 “clang”（调用 LLVM/Clang 编译器框架）和 “shell”（通过命令行调用系统编译器）：

C = Importer('gen.c','shell')
f = external('f',C)
print(f(3.14))
0.00159265

我们将在第 6.3 节再次讨论外部函数。

5.2.2. 从 MATLAB 调用生成的代码

执行生成代码的另一种方法是将代码编译成 MATLAB MEX 函数并从 MATLAB 中调用。前提是在代码生成过程中将 mex 选项设置为 “true”（参见第 5.1 节）。生成的 MEX 函数将函数名作为第一个参数，然后是函数输入：

%mex gen.c -largeArrayDims  % Matlab
mex gen.c -DMATLAB_MEX_FILE % Octave

disp(gen('f', 3.14))
1.5927e-03

请注意，默认情况下执行结果总是 MATLAB 稀疏矩阵。可以设置编译器标志 -DCASASI_MEX_ALWAYS_DENSE 和 -DCASASI_MEX_ALLOW_DENSE 来影响这种行为。

5.2.3. 从命令行调用生成的代码

另一个选择是从 Linux/OSX 命令行执行生成的代码。如果在代码生成过程中将主选项设置为 “true”（参见第 5.1 节），则可以这样做。这对使用 gprof 等工具对生成的代码进行剖析非常有用。

执行生成的代码时，函数名称将作为命令行参数传递。所有输入的非零条目需要通过标准输入传递，函数将通过标准输出返回所有输出的非零条目：

# Command line
echo 3.14 3.14 > gen_in.txt
gcc gen.c -o gen
./gen f < gen_in.txt > gen_out.txt
cat gen_out.txt # returns 0.00159265 0.00159265

5.2.4. 与 C/C++ 应用程序生成的代码部件链接¶。

生成的代码可以直接与 C/C++ 应用程序的部件链接。如果在代码生成过程中将 with_header 选项设置为 “true”，则会生成一个头文件，其中包含文件中所有暴露入口点的声明。使用该头文件需要了解 CasADi 的代码生成 API，详见下文第 5.3 节。未暴露的符号会以文件特定的前缀作为前缀，允许应用程序针对多个生成函数进行部件链接，而不会有符号冲突的风险。

5.2.5. 从 C/C++ 应用程序中动态加载生成的代码¶。

上述方法的一个变种是将生成的代码编译到共享库中，但直接访问暴露的符号，而不是依赖 CasADi 的外部函数。这也需要了解生成代码的结构。

在 CasADi 示例库中，codegen_usage.cpp 演示了如何做到这一点。

5.3. 生成代码的应用程序接口

生成代码的 API 由许多具有外部链接的函数组成。除了实际执行外，还有内存管理以及输入和输出元信息的函数。下文将介绍这些函数。下面，假设我们要访问的函数名称是 fname。要了解这些函数在代码中的实际样子以及何时被调用，请参考 codegen_usage.cpp 示例。

5.3.1. 引用计数

void fname_incref(void);
void fname_decref(void);

生成的函数可能需要在首次调用前读入某些数据或初始化某些数据结构。由 CasADi 表达式生成的函数通常不需要这样做，但当生成的代码包含对外部函数的调用时，可能需要这样做。同样，内存在使用后也可能需要重新分配。

为了跟踪所有权，生成的代码包含两个用于增加和减少引用计数器的函数。它们分别被命名为 fname_incref 和 fname_decref。这些函数没有输入参数，返回 void。

通常，第一次调用 fname_incref 时可能会进行一些初始化，随后的调用只会增加一些内部计数器。另一方面，fname_decref 会减少内部计数器，当计数器为零时，将进行取消分配（如果有的话）。

5.3.2. 输入和输出的数量

casadi_int fname_n_in(void);
casadi_int fname_n_out(void);

调用 fname_n_in 和 fname_n_out 函数可分别获得函数的输入和输出数。这些函数不需要输入，并返回输入或输出的数量（casadi_int 是 long long int 的别名）。

5.3.3. 输入和输出的名称

const char* fname_name_in(casadi_int ind);
const char* fname_name_out(casadi_int ind);

函数 fname_name_in 和 fname_name_out 返回特定输入或输出的名称。它们获取输入或输出的索引（从索引 0 开始），并返回一个带有空端 C 字符串名称的 const char*。如果函数执行失败，将返回一个空指针。

5.3.4. 输入和输出的稀疏性模式

const casadi_int* fname_sparsity_in(casadi_int ind);
const casadi_int* fname_sparsity_out(casadi_int ind);

调用 fname_sparsity_in 和 fname_sparsity_out 可以分别获得给定输入或输出的稀疏性模式。这些函数获取输入或输出索引，并返回一个指向常整数字段（const casadi_int*）的指针。这是 CasADi 使用的压缩列存储（CCS）格式的紧凑表示，参见第 3.5 节。指向的整数字段结构如下：

前两项分别是行数和列数。在下文中称为 nrow 和 ncol。

如果第三个条目为 1，则表示该模式为密集模式，其余条目将被丢弃。

如果第三个条目为 0，则该条目加上随后的 ncol 条目构成每列的非零偏移量，即下文中的 colind。例如 i 将由 colind[i] 到 colind[i+1] 之间的非零指数组成。最后一项 colind[ncol] 将等于非零点数 nnz。

最后，如果稀疏性模式不密集，即 nnz ≠ nrow * ncol，那么最后 nnz 条目将包含行索引。

如果失败，这些函数将返回一个空指针。

5.3.5. 内存对象

函数可能包含一些可变内存，例如，用于缓存最新因式分解或跟踪评估统计。当多个函数需要调用同一个函数而不发生冲突时，每个函数都需要使用不同的内存对象。这对于在共享内存架构上进行并行计算尤为重要，在这种情况下，每个线程都应访问不同的内存对象。

void* fname_alloc_mem(void);

分配一个内存对象，并将其传递给数值计算。

int fname_init_mem(void* mem);

(重新）初始化内存对象。成功后返回 0；

int fname_free_mem(void* mem);

释放内存对象。成功后返回 0；

5.3.6. 工作向量

int fname_work(casadi_int* sz_arg, casadi_int* sz_res, casadi_int* sz_iw, casadi_int* sz_w);

为了在占用较少内存的情况下高效地执行计算，用户需要传递四个工作数组。函数 fname_work 返回这些数组的长度，这些数组的条目类型分别为 const double*、double*、casadi_int 和 double。

函数失败时的返回值为非零。

5.3.7. 数值评估

int fname(const double** arg, double** res,
          casadi_int* iw, double* w, void* mem);

最后，函数 fname 执行实际的计算。它的输入参数包括四个工作矢量和一个使用 fname_alloc_mem 创建的内存对象（如果没有，则为 NULL）。工作矢量的长度必须至少等于 fname_work 命令提供的长度，内存对象的索引必须严格小于 fname_n_mem 返回的值。

函数输入的非零点由 arg 工作矢量的第一个条目指向，并且在求值过程中保持不变。同样，输出的非零值由 res 工作向量的首项指向，并且也保持不变（即指针不变，实际值不变）。

函数失败时的返回值为非零。

脚注

[1]
当生成的函数代码包含对外部函数的调用时是个例外。

六、用户自定义函数对象¶

在某些情况下，使用 CasADi 符号重写用户函数是不可能或不切实际的。为了解决这个问题，CasADi 提供了许多方法来嵌入调用用 CasADi 语言（C++、MATLAB 或 Python）或 C 语言定义的 "黑盒 "函数。使用 CasADi 符号定义的函数几乎总是更高效，尤其是涉及导数计算时，因为通常可以利用更多的结构。

根据具体情况，用户可以通过多种不同方式实现自定义函数对象，下文将详细介绍：

子类化 FunctionInternal 第 6.1 节

子类化回调 第 6.2 节

导入外部函数 第 6.3 节

即时编译 C 语言字符串 第 6.4 节

用查找表替换函数调用 第 6.5 节

6.1. 函数内部子类

第 4 章中介绍的所有函数对象在 CasADi 中都是作为继承于 FunctionInternal 抽象基类的 C++ 类实现的。原则上，熟悉 C++ 编程的用户可以实现一个继承自 FunctionInternal 的类，并重载该类的虚拟方法。这样做的最佳参考资料是 C++ API 文档，选择 "switch to internal "来公开内部 API。

由于 FunctionInternal 并不属于稳定、公开的 API，我们建议一般情况下不要这样做，除非计划贡献 CasADi 的源代码。

6.2 子类化 Callback¶

Callback 类为 FunctionInternal 提供了公共 API，从该类继承与直接从 FunctionInternal 继承效果相同。得益于跨语言多态性，无论是 Python、MATLAB/Octave 还是 C++，都可以实现 Callback 的公开方法。

派生类由以下部分组成：

一个构造函数或一个替代构造函数的静态函数

多个虚拟函数，均为可选函数，可通过重载获得所需的行为。其中包括使用 get_n_in 和 get_n_out 计算输入和输出的数量，使用 get_name_in 和 get_name_out 计算它们的名称，以及使用 get_sparsity_in 和 get_sparsity_out 计算它们的稀疏性模式。

一个可选的 init 函数，在对象构造过程中调用。

数值计算函数

可选的导数函数。你可以选择使用有限差分（enable_fd），提供完整的雅各布因子（has_jacobian、get_jacobian），或者选择提供正向/反向敏感度（has_forward、get_forward、has_reverse、get_reverse）。

有关函数的完整列表，请参阅 Callback 的 C++ API 文档。另请参阅 callback.py 示例。

下文将介绍不同语言的用法。

6.2.1. Python* 语言

class MyCallback(Callback):
  def __init__(self, name, d, opts={}):
    Callback.__init__(self)
    self.d = d
    self.construct(name, opts)

  # Number of inputs and outputs
  def get_n_in(self): return 1
  def get_n_out(self): return 1

  # Initialize the object
  def init(self):
     print('initializing object')

  # Evaluate numerically
  def eval(self, arg):
    x = arg[0]
    f = sin(self.d*x)
    return [f]

实现应包括一个构造函数，该构造函数应以调用基类构造函数 Callback.init(self) 开始，以调用 self.construct(name, opts) 初始化对象构造结束。

该函数可与任何 CasADi 内置函数一样使用，但有一个重要的注意事项，即当嵌入图形时，类的所有权不会在所有引用之间共享。因此，当计算中仍然需要 Python 类时，用户务必不要让它离开作用域。

# Use the function
f = MyCallback('f', 0.5)
print(f(2))

x = MX.sym("x")
print(f(x))
initializing object
0.841471
f(x){0}

6.2.2. MATLAB¶

在 MATLAB 中，自定义函数类可在 MyCallback.m 文件中定义如下：

classdef MyCallback < casadi.Callback
  properties
    d
  end
  methods
    function self = MyCallback(name, d)
      self@casadi.Callback();
      self.d = d;
      construct(self, name);
    end

    % Number of inputs and outputs
    function v=get_n_in(self)
      v=1;
    end
    function v=get_n_out(self)
      v=1;
    end

    % Initialize the object
    function init(self)
      disp('initializing object')
    end

    % Evaluate numerically
    function arg = eval(self, arg)
      x = arg{1};
      f = sin(self.d * x);
      arg = {f};
    end
  end
end

该函数可以像任何内置的 CasADi 函数一样使用，但对于 Python 来说，类的所有权不会在所有引用之间共享。因此，用户必须避免在类实例仍在使用时将其删除，例如将其设置为持久化。

% Use the function
f = MyCallback('f', 0.5);
res = f(2);
disp(res)

x = MX.sym('x');
disp(f(x))

6.2.3. C++¶

在 C++ 中，语法如下：

#include "casadi/casadi.hpp"
using namespace casadi;
class MyCallback : public Callback {
  // Data members
  double d;
public:
  // Constructor
  MyCallback(const std::string& name, double d,
             const Dict& opts=Dict()) : d(d) {
    construct(name, opts);
  }

  // Destructor
  ~MyCallback() override {}

  // Number of inputs and outputs
  casadi_int get_n_in() override { return 1;}
  casadi_int get_n_out() override { return 1;}

  // Initialize the object
  void init() override() {
    std::cout << "initializing object" << std::endl;
  }

  // Evaluate numerically
  std::vector<DM> eval(const std::vector<DM>& arg) const override {
    DM x = arg.at(0);
    DM f = sin(d*x);
    return {f};
  }
};

以这种方式创建的类可以像其他函数实例一样使用，但有一个重要区别，即用户负责管理该类的内存。

int main() {
  MyCallback f("f", 0.5);
  std::vector<DM> arg = {2};
  std::vector<DM> res = f(arg);
  std::cout << res << std::endl;
  return 0;
}

6.3. 导入外部函数

CasADi 外部函数的基本用法已在第 5.2 节的自动生成代码中进行了演示。同样的函数也可用于导入用户定义的函数，只要该函数也使用第 5.3 节所述的 C API。

下文将对此进行进一步阐述。

6.3.1. 默认函数

通常无需定义第 5.3 节中定义的所有函数。如果没有定义 fname_incref 和 fname_decref，则假定不需要内存管理。如果没有提供输入和输出的名称，它们将被赋予默认名称。一般情况下，稀疏性模式默认为标量模式，除非函数对应于另一个函数的导数（见下文），在这种情况下，稀疏性模式默认为密集模式且维度正确。

此外，如果没有实现 fname_work，则假定不需要工作矢量。

6.3.2. 作为注释的元信息

如果您依赖于 CasADi 的即时编译器，您可以在 C 代码中以注释的形式提供元信息，而不是实现实际的回调函数。

元信息的结构如下：

/*CASADIMETA
:fname_N_IN 1
:fname_N_OUT 2
:fname_NAME_IN[0] x
:fname_NAME_OUT[0] r
:fname_NAME_OUT[1] s
:fname_SPARSITY_IN[0] 2 1 0 2
*/

6.3.3. 导数

外部函数可以通过提供导数计算函数来实现微分。在导数计算过程中，CasADi 将在同一文件/共享库中寻找符合特定命名规则的符号。例如，您可以通过执行名为 jac_fname 的函数，为名为 fname 的函数指定相对于所有输入的所有输出的雅各比。同样，您可以通过提供名为 fwd1_fname 的函数来指定计算一个正向导数的函数，其中 1 可以用 2、4、8、16、32 或 64 代替，以便一次计算多个正向导数。对于反向模式方向导数，可将 fwd 替换为 adj。

这是一项试验性功能。

6.4. 即时编译 C 语言字符串

在上一节中，我们介绍了如何指定一个包含数值计算函数和元信息的 C 语言文件。如前所述，CasADi 与 Clang 的接口可以对该文件进行即时编译。用户只需将源代码指定为一个 C 语言字符串即可。

body =\
'r[0] = x[0];'+\
'while (r[0]<s[0]) {'+\
' r[0] *= r[0];'+\
'}'

f = Function.jit('f',body,\
      ['x','s'],['r'],\
      {"compiler":"shell"})
print(f)
f:(x,s)->(r) JitFunction

Function.jit 的这四个参数是必须的： 函数名称、字符串形式的 C 代码以及输入和输出名称。在 C 代码中，输入/输出名称对应于 casadi_real_t 类型的数组，其中包含函数输入和输出的非零元素。默认情况下，所有输入和输出都是标量（即 1-by-1 且密集）。要指定不同的稀疏性模式，请提供两个额外的函数参数，其中包含稀疏性模式的向量/列表：

sp = Sparsity.scalar()
f = Function.jit('f',body,\
     ['x','s'], ['r'],\
     [sp,sp], [sp],\
     {"compiler":"shell"})
print(f)
f:(x,s)->(r) JitFunction

两种变体都接受以选项字典形式出现的第 5 个（或第 7 个）可选参数。

6.5. 使用查找表

可以使用 CasADi 的插值函数创建查找表。不同的插值方案以插件形式实现，类似于 nlpsol 或 integrator 对象。除了标识符名称和插件外，插值函数还需要一组网格点和相应的数值。

内插函数调用的结果是一个可微分的 CasADi 函数对象，可以通过调用 MX 参数嵌入到 CasADi 计算图形中。此外，此类图形完全支持 C 代码生成。

目前，插值法有两个插件： 线性 "和 “bspline”。它们的作用类似于 MATLAB/Octave 的插值，方法设置为 "线性 "或 “样条”–对应于多线性插值和非结边界条件下的样条插值（默认为立方）。

在使用 bspline 的情况下，将在构建时寻找符合所提供数据的系数。或者，您也可以使用更低级的 Function.bspline 来自行提供系数。默认情况下，bspline 的每个维度的阶数都是 3。您可以通过阶数选项偏离默认值。

我们将介绍一维和二维版本的插值法语法，但实际上该语法可用于任意维数。

6.5.1. 一维查找表

与 SciPy 中的相应方法相比，CasADi/Python 中的一维样条曲线拟合方法如下：

import casadi as ca
import numpy as np
xgrid = np.linspace(1,6,6)
V = [-1,-1,-2,-3,0,2]
lut = ca.interpolant('LUT','bspline',[xgrid],V)
print(lut(2.5))
# Using SciPy
import scipy.interpolate as ip
interp = ip.InterpolatedUnivariateSpline(xgrid, V)
print(interp(2.5))
-1.35
-1.3500000000000005

在 MATLAB/Octave 中，相应的代码为

xgrid = 1:6;
V = [-1 -1 -2 -3 0 2];
lut = casadi.interpolant('LUT','bspline',{xgrid},V);
lut(2.5)
% Using MATLAB/Octave builtin
interp1(xgrid,V,2.5,'spline')
ans =

-1.35

ans = -1.3500

请特别注意，内插法的网格和数值参数必须是数值参数。

6.5.2. 二维查找表

在 Python 中，我们可以得到以下二维查找表，也可以与 SciPy 进行比较，以供参考：

xgrid = np.linspace(-5,5,11)
ygrid = np.linspace(-4,4,9)
X,Y = np.meshgrid(xgrid,ygrid,indexing='ij')
R = np.sqrt(5*X**2 + Y**2)+ 1
data = np.sin(R)/R
data_flat = data.ravel(order='F')
lut = ca.interpolant('name','bspline',[xgrid,ygrid],data_flat)
print(lut([0.5,1]))
# Using Scipy

interp = ip.RectBivariateSpline(xgrid, ygrid, data)
print(interp.ev(0.5,1))
0.245507
0.24550661674668917

或在 MATLAB/Octave 中与内置函数进行比较：

xgrid = -5:1:5;
ygrid = -4:1:4;
[X,Y] = ndgrid(xgrid, ygrid);
R = sqrt(5*X.^2 + Y.^2)+ 1;
V = sin(R)./R;
lut = interpolant('LUT','bspline',{xgrid, ygrid},V(:));
lut([0.5 1])
% Using Matlab builtin
interpn(X,Y,V,0.5,1,'spline')
ans =

0.245507

ans = 0.2455

特别要注意的是，values 参数必须扁平化为一维数组。

6.6. 使用有限差分进行衍生计算

CasADi 3.3 为所有函数对象引入了有限差分计算支持，特别是包括第 6.2 节、第 6.3 节或第 6.4 节中定义的外部函数（对于第 6.5 节中的查找表，可使用解析导数）。

除 Function.jit 外，有限差分导数默认为禁用，要启用有限差分导数，必须将 "enable_fd "选项设置为 True/true：

f = external('f', './gen.so',\
   dict(enable_fd=True))

e = jacobian(f(x),x)
D = Function('d',[x],[e])
print(D(0))
1

参见第 5.1 节。

如果没有分析导数，"enable_fd "选项可以让 CasADi 使用有限差分。要强制 CasADi 使用有限差分，可以将 “enable_fd”、"enable_reverse "和 "enable_jacobian "设置为 “假”/“假”，分别对应于 CasADi 所使用的三种分析导数信息。

默认方法是中心差法，步长由函数的舍入误差和截断误差估算值决定。您可以通过将 "fd_method "选项设置为 “forward”（对应一阶正向差分）、“backward”（对应一阶反向差分）和 “smoothing”（二阶精确非连续性避免方案）来更改方法，适用于需要计算域边缘导数的情况。通过设置 "fd_options "选项，可以获得有限差分的其他算法选项。

七、 DaeBuilder 类

CasADi 中的 DaeBuilder 类是一个辅助类，旨在方便复杂动力系统的建模，以便日后与最优控制算法配合使用。该类的抽象层级低于物理建模语言，如 Modelica（参见第 7.2 节），但仍高于直接使用 CasADi 符号表达式的层级。特别是，它可以用来连接以功能模拟接口（FMI）格式提供的物理模型，也可以用作构建特定领域建模环境的构件。

DaeBuilder 原则上有三种不同的使用方法：
可以使用该类逐步构建微分代数方程（DAE）结构系统，然后将其用于连接 CasADi 中的仿真或优化。该类还支持以 FMI 格式导出模型，以便在其他工具中使用。

还可以使用该类导入 FMI 格式的现有模型。从 CasADi 3.6 开始，我们支持导入标准 FMU，其中的模型方程只能通过调用 DLL 来获得。FMU 导入支持已针对具有挑战性的 FMU 进行了测试，但截至本文撰写时仍在开发中。

我们可以使用该类以 XML 符号格式导入现有模型。JModelica.org 工具支持这种格式，OpenModelica 也试验性地支持这种格式。由于生成这种格式的模型的工具有限，因此这种格式没有得到积极的开发。

7.1. 构建 DaeBuilder 实例

下面是一个简单的 DAE，它对应于一个受二次空气摩擦项和重力作用的可控火箭，火箭在耗尽燃料后质量会逐渐减小：
$$\begin{array}{c}h=v,h(0)=0\\ \dot{v}=(u-a{v}^2)/m-g,v(0)=0\\ \dot{m}=-b{u}^2,m(0)=1\end{array}$$
其中，三种状态分别对应高度、速度和质量。
是火箭的推力，（a、 b ) 为参数。

要构建该问题的 DAE 公式，请从一个空的 DaeBuilder 实例开始，按如下步骤逐步添加输入和输出表达式。

dae = DaeBuilder('rocket')
# Add input expressions
a = dae.add_p('a')
b = dae.add_p('b')
u = dae.add_u('u')
h = dae.add_x('h')
v = dae.add_x('v')
m = dae.add_x('m')
# Add output expressions
hdot = v
vdot = (u-a*v**2)/m-g
mdot = -b*u**2
dae.set_ode('h', hdot)
dae.set_ode('v', vdot)
dae.set_ode('m', mdot)
# Specify initial conditions
dae.set_start('h', 0)
dae.set_start('v', 0)
dae.set_start('m', 1)
# Add meta information
dae.set_unit('h','m')
dae.set_unit('v','m/s')
dae.set_unit('m','kg')

其他输入和输出表达式也可以类似方式添加。有关函数的完整列表，请参见 DaeBuilder 的 C++ API 文档。

7.2. 从 Modelica 符号导入 OCP*

注：Modelon 不再提供 JModelica.org。但其闭源后继代码 OCT 保留了 CasADi 的互操作性。关于如何使用 OCT 生成 CasADi 表达式的详细信息，请参阅 OCT 的用户指南。以下文字指的是传统的 Modelica 互操作性支持。

7.2.1. 从 XML 文件的传统符号导入

除了直接在 CasADi 中建模（如上所述），还有一种方法是使用高级物理建模语言（如 Modelica）来指定模型。为此，CasADi 提供了与开源 JModelica.org 编译器的互操作性，该编译器是专门为优化控制而编写的。从 JModelica.org 导入模型有两种不同的方式：使用 JModelica.org 的 CasadiInterface 或通过 DaeBuilder 的 parse_fmi 命令。

我们推荐使用前一种方法，因为它正在得到积极维护，有关如何提取 CasADi 表达式的详细信息，请参阅 JModelica.org 的用户指南。

下面，我们将概述使用 parse_fmi 的传统方法。

7.2.2. 从传统方法导入 modelDescription.xml 文件¶。

要了解如何使用 Modelica 导入，请查看 CasADi 示例集中的 thermodynamics_example.py。

假设 Modelica/Optimica 模型 ModelicaClass.ModelicaModel 定义在 file1.mo 和 file2.mop 文件中，Python 编译命令为

from pymodelica import compile_jmu
jmu_name=compile_jmu('ModelicaClass.ModelicaModel', \
  ['file1.mo','file2.mop'],'auto','ipopt',\
  {'generate_xml_equations':True, 'generate_fmi_me_xml':False})

这将生成一个 jmu 文件，实质上是一个压缩文件，其中包含 modelDescription.xml 文件。该 XML 文件包含最优控制问题的符号表示，可在标准 XML 编辑器中查看。

from zipfile import ZipFile
sfile = ZipFile(jmu_name','r')
mfile = sfile.extract('modelDescription.xml','.')

一旦有了 modelDescription.xml 文件，就可以使用 parse_fmi 命令将其导入：

ocp = DaeBuilder()
ocp.parse_fmi('modelDescription.xml')

7.3. 函数工厂

一旦一个 DaeBuilder 被提出，并可能被重新表述为一个令人满意的形式，我们就可以生成与 sec-daebuilder_io 中列出的输入和输出表达式相对应的 CasADi 函数。例如，要为第 7.1 节中的火箭模型的 ODE 右侧创建一个函数，只需提供所创建函数的显示名称、输入表达式列表和输出表达式列表：

f = dae.create('f',\
     ['x','u','p'],['ode'])

利用命名规则，我们还可以创建雅各布函数，例如 "ode "输出相对于 "x "的雅各布函数：

f = dae.create('f',\
     ['x','u','p'],\
     ['jac_ode_x'])

首先使用 add_lc 创建输出表达式的命名线性组合，然后请求其 Hessian，即可提取具有二阶信息的函数：

dae.add_lc('gamma',['ode'])
hes = dae.create('hes',\
  ['x','u','p','lam_ode'],\
  ['hes_gamma_x_x'])

也可以简单地从 DaeBuilder 实例中提取符号表达式，然后手动创建 CasADi 函数。例如，dae.x 包含与 "x "相对应的所有表达式，dae.ode 包含与 "ode "相对应的表达式，等等。

脚注

[1]
FMI 开发小组。Functional Mock-up Interface for Model Exchange and Co-Simulation. https://www.fmi-standard.org/, July 2014. 规范，FMI 2.0。第 3.1 节，第 71-72 页。

八、使用 CasADi 进行优化控制

CasADi 可用于使用多种方法解决最优控制问题（OCP），包括直接（又称离散化-优化）和间接（又称优化-离散化）方法、一次求解（如搭配）方法以及需要嵌入 ODE 或 DAE 初值求解器的射击方法。作为用户，通常需要编写自己的 OCP 求解器，而 CasADi 的目标是通过提供功能强大的高级构建模块，尽可能简化这一过程。由于是自己编写求解器（而不是调用现有的 "黑盒 "求解器），因此对如何求解 OCP 的基本了解是必不可少的。Biegler [1] 或 Betts 最近出版的教科书对数值最优控制做了很好的、自成一体的介绍。
[1] 或 Betts
[2] 或 Moritz Diehl 的数值最优控制讲义。

8.1. 一个简单的测试问题

为了说明某些方法，我们将考虑以下测试问题，即驱动一个范德尔波尔振荡器到原点，同时试图最小化二次费用：
$\begin{array}{rl}& \text{minimize:}x(\cdot )\in {\mathbb{R}}^2,u(\cdot )\in \mathbb{R}{\int }_{t=0}^T\left(x_0^2+x_1^2+u^2\right)dt\\ & \text{subject to:}\\ & \{\begin{array}{rl}{\dot{x}}_0& =\left(1-x_1^2\right)x_0-x_1+u\\ {\dot{x}}_1& =x_0\\ -1.0\le u\le 1.0,x_1\ge -0.25& \end{array}\mathrm{for}0\le t\le T\\ & x_0(0)=0,x_1(0)=1,\end{array}$

与 T = 10 .

在 CasADi 的示例集
[3] 中，您可以找到使用各种不同方法求解最优控制问题的代码。

下面，我们将讨论三种最重要的方法，即直接单射、直接多射和直接搭配。

8.1.1. 直接单射法

在直接单次射击法中，控制轨迹使用某种片断平滑近似值（通常是片断常数）进行参数化。

使用控制的显式表达，我们就可以从优化问题中消除整个状态轨迹，最终只对离散化的控制进行 NLP。

在 CasADi 的示例集中，您可以找到分别用于 Python 和 MATLAB/Octave 的 direct_single_shooting.py 和 direct_single_shooting.m 代码。这些代码实现了直接单次射击法，并使用 IPOPT 进行求解，依靠 CasADi 计算导数。为了从连续时间动力学获得离散时间动力学，使用 CasADi 符号实现了一个简单的固定步 Runge-Kutta 4 (RK4) 积分器。类似这样的简单积分器代码在优化控制中通常很有用，但必须注意准确解决初值问题。

该代码还展示了如何用一种更先进的积分器（即 SUNDIALS 软件包中的 CVODES 积分器）取代 RK4 方案，后者实现了一种可变步长、可变阶次的后向微分公式 (BDF) 方案。这种高级积分器适用于较大系统、具有刚性动力学的系统、DAEs 以及检查较简单方案的一致性。

8.1.2. 直接多重射击

在 CasADi 示例集中的 direct_multiple_shooting.py 和 direct_multiple_shooting.m 代码实现了直接多重射击法。这种方法与直接单次射击法非常相似，但将某些射击节点的状态作为 NLP 中的决策变量，并包含相等约束以确保轨迹的连续性。

由于将问题 "提升 "到一个更高的维度通常会提高收敛性，因此直接多重射击法通常优于直接单次射击法。用户还可以使用已知的状态轨迹猜测进行初始化。

缺点是所求解的 NLP 会变得更大，不过这通常会被更稀疏的事实所弥补。

8.1.3. 直接定位

最后，direct_collocation.py 和 direct_collocation.m 代码实现了直接搭配法。在这种情况下，整个状态轨迹的参数化（如片断低阶多项式）将作为 NLP 的决策变量。这样就完全不需要离散时间动力学公式了。

直接配准的 NLP 比直接多重射击的 NLP 更大，但也更稀疏。

脚注

[1]
Lorenz T. Biegler, Nonlinear Programming： Concepts, Algorithms, and Applications to Chemical Processes , SIAM 2010

[2]
John T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming , SIAM 2001

[3]
您可以在 CasADi 的下载区以名为 examples_pack.zip 的压缩包形式获取此合集。

九、Opti 栈

Opti 栈是一组 CasADi 辅助类，提供了数学 NLP 符号之间的紧密对应关系，例如
$\begin{array}{rlr}& \text{minimize}& \\ & & (y-x^2{)}^2\\ & x,y& \\ & \text{subject to}& x^2+y^2=1\end{array}$
和计算机代码：