文章目錄

• 1. 常用符號
• • 1.1 SX 符號
  • 1.2 DM 符號
  • 1.3 MX 符號
  • 1.4 SX 和 MX 的混合
  • 1.5 稀疏類
  • • 1.5.1 獲取和設置矩陣中的元素
  • 1.6 算術運算
  • • 1.6.1 乘法
    • 1.6.2 transpose
    • 1.6.3 reshape()
    • 1.6.4 vertcat()/horzcat()
    • 1.6.5 split
    • 1.6.6 Inner product
  • 1.7 查詢屬性
  • 1.8 線性代數
  • 1.9 微積分 - 算法微分
  • • 1.9.1 句法

---

1. 常用符號

1.1 SX 符號

??SX 數據類型用于表示矩陣，其元素由一系列一元和二元運算序列組成的符號表達式組成。要查看它在實踐中的工作原理，請啟動交互式 Python shell（例如，通過從 Linux 終端或在集成開發環境（如 Spyder）中鍵入 ipython）或啟動 MATLAB 或 Octave 的圖形用戶界面。 假設 CasADi 已正確安裝，您可以將該庫導入工作區，如下所示：

from casadi import *

使用以下語法創建一個變量 x：

x = MX.sym("x") print(x)

運行后輸出結果為：

x

上述代碼創建一個 1×1 矩陣，即一個名字叫 x 的符號基元（是個標量）。 這只是顯示名稱，而不是標識符。 多個變量可以具有相同的名稱，但仍然不同。 標識符是返回值。 此外，還可以調用函數 SX.sym （帶參數）來創建向量或矩陣值符號變量：

y = SX.sym('y',5) Z = SX.sym('Z',4,2) print(y) print(z)

運行后輸出結果為：

[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): # 稀疏的4×4矩陣，對角線為1 print(B1) print(B2) print(B4)

運行后輸出結果為：

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，以將它們與實際零 @1 （實際顯式賦值的元素，print后使用@1來表示0，以便區分結構零與實際零之間的區別）區分開來。

除了上述的使用方法外，下表羅列了常用的 SX 表達式：

函數表達式用途

SX.sym(name,n,m)	創建一個 n×m 符號基元
SX.zeros（n，m）	創建一個所有元素全為零的n×m密集矩陣
SX(n,m)	創建一個所有元素為零的 n×m 稀疏矩陣
SX.ones(n,m)	創建一個包含所有 1 的 n×m 密集矩陣
SX.eye(n)	創建一個 n×n 對角矩陣，對角線元素為1，其他地方元素為結構零
SX(scalar_type)	使用參數給定的值創建一個標量（1×1 矩陣），可以顯式使用此方法，例如 SX(9)，或隱式，例如9 * SX.ones(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 矩陣。
SX(list) （僅限 Python）	使用列表中的元素創建列向量（n×1 矩陣），例如SX([1,2,3,4])。（注意 Python 列表和 MATLAB/Octave 水平串聯的區別，兩者都使用方括號語法）
SX(list of list)（僅限 Python）	使用列表中的元素創建一個密集矩陣，例如SX([[1,2],[3,4]]) 或使用 SX([[1,2,3,4]]) 創建一個（1×n 矩陣）行向量。

---

1.2 DM 符號

??DM 與 SX 非常相似，但不同的是非零元素是數值而不是符號表達式。 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/ 的示例包中找到。 有關此類（和其他類）特定功能的文檔，例如 C++ API， 可以在http://docs.casadi.org/上找到。

---

1.3 MX 符號

??下面代碼的輸出是一個 2×2 矩陣，注意看這里的第一句話@1=3，也就意味著第二行的輸出中@1代表3，輸出結果意味著CasADi為矩陣f的每個條目創建了新的表達式（SX 類型）。

x = SX.sym('x',2,2) y = SX.sym('y') f = 3*x + y print(f) print(f.shape)

運行后輸出結果為：

@1=3,
[[((@1x_0)+y), ((@1x_2)+y)],
?[((@1x_1)+y), ((@1x_3)+y)]]
(2, 2)

??除了上面的方式外，我們還可以使用第二種更通用的矩陣表達式類型MX。 MX 類型允許像 SX 一樣構建由一系列基本操作組成的表達式。 但與 SX 不同的是，這些基本運算不限于標量一元或二元運算（$\mathbb{R}\to \mathbb{R}$ 或 $\mathbb{R}\times \mathbb{R}\to \mathbb{R}$）。 取而代之的是，用于形成MX表達式的基本運算可以為通用的多個稀疏矩陣值輸入，多個稀疏矩陣值輸出函數： ${\mathbb{R}}^{n_1\times m_1}\times \dots \times {\mathbb{R}}^{n_N\times m_N}\to {\mathbb{R}}^{p_1\times q_1}\times \dots \times {\mathbb{R}}^{p_M\times q_M}$，例如：

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 等價物有八個（每個元素由兩個構成，f(1,1)由(@1x_0)和y構成）。 因此，當使用具有許多元素的自然向量或矩陣值的操作時，MX 可以更經濟。

??MX支持使用與SX相同的語法來獲取和設置元素，但是實現方式卻大不相同。 例如，測試打印 2×2 符號變量左上角的元素：

x = MX.sym('x',2,2) print(x[0,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。

---

1.4 SX 和 MX 的混合

??不能將SX對象與MX對象相乘，也不能執行任何其他操作以將兩者混合在同一表達式圖中。 但是卻可以在 MX 圖中包括對由 SX 表達式定義的函數的調用。如果需要提升程序的效率，可以混合使用 SX 和 MX，因為由 SX 表達式定義的函數每個操作的開銷要低得多，這使得系統完成一系列標量操作的速度更快。 因此，SX 表達式旨在用于低級操作，而 MX 表達式則用于實現負責的操作，例如實現 NLP 的約束函數。

---

1.5 稀疏類

??CasADi 中的矩陣使用壓縮列存儲 (CCS) 格式進行存儲。 這是稀疏矩陣的標準格式，允許有效執行線性代數運算，例如逐元素運算、矩陣乘法和轉置。 在 CCS 格式中，稀疏模式使用維度（行數和列數）和兩個向量進行解碼。 第一個向量包含每列的第一個結構上非零元素的索引，第二個向量包含每個非零元素的行索引。 有關 CCS 格式的更多詳細信息，請參見例如 Netlib 上線性系統解決方案的模板。 請注意，CasADi 對稀疏矩陣和密集矩陣使用 CCS 格式。

??CasADi 中的稀疏模式存儲為 Sparsity 類的實例，該類是引用計數的，這意味著多個矩陣可以共享相同的稀疏模式，包括 MX 表達式圖以及 SX 和 DM 的實例。 稀疏類也被緩存，這意味著總是避免創建相同稀疏模式的多個實例。

??以下列表總結了構建新稀疏模式的最常用方法：

函數名功能

Sparsity.dense(n,m)	創建一個密集的 n×m 稀疏模式
Sparsity(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]]

1.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)

運行后輸出結果為：

[[3, 7]]
@1=1,
[[@1, @1],
?[4, 5]]

??與 Python 的 NumPy 不同，CasADi 片段不是左側數據的視圖； 相反，片段訪問（Slice access）復制數據。 因此，在以下示例中，矩陣 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)

運行后輸出結果為：

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

??片段訪問（Slice access）意味著一次設置多個元素。 這比一次設置一個元素要高效得多。 可以通過提供 (start , stop , step) 元組來獲取或設置片段。 在 Python 和 MATLAB 中，CasADi 使用標準語法：

print(M[:,1]) print(M[1:,1:4:2]) # 對應(start , stop , step)

運行后輸出結果為：

[00, 4, 00, 00]

[[4, 00],
?[00, 00],
?[00, 6]]

??列表訪問類似于（但效率可能低于）片段訪問：

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]

---

1.6 算術運算

1.6.1 乘法

??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 中的慣例，當任一參數為標量時，使用 * 和 .* 的乘法是等效的。

1.6.2 transpose

??轉置使用 Python 中的語法 A.T、C++ 中的 A.T() 以及 A 或 A 形成。 在 MATLAB 中：

print(y) print(y.T)

運行后輸出結果為：

[[y_0, y_2],
?[y_1, y_3]]

[[y_0, y_1],
?[y_2, y_3]]

1.6.3 reshape()

??重塑（reshape）意味著改變行數和列數，但保留元素數和非零元素的相對位置。 這是一個計算成本非常低的操作，它使用以下語法執行：

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]]

1.6.4 vertcat()/horzcat()

??串聯意味著水平或垂直堆疊矩陣。 由于在 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))

運行后輸出結果為：

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

1.6.5 split

??水平和垂直拆分是上面介紹的水平和垂直串聯的逆操作。 要將表達式水平拆分為 n 個較小的表達式，除了要拆分的表達式之外，還需要提供長度為 n+1 的向量偏移。 偏移向量的第一個元素必須是 0，最后一個元素必須是列數。 其余元素必須按非遞減順序排列。然后，拆分操作的輸出 i 包含具有 $offset[i]\le c<offset[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]

??vertsplit 操作的工作原理類似，但偏移向量指的是行：

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 的效率要高得多。

1.6.6 Inner product

??內積，定義為 $<A,B>:=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))

??許多上述操作也是為稀疏類定義的，經過 vertcat、horzsplit、轉置、加法（返回兩個稀疏模式的并集）和乘法（返回兩個稀疏模式的交集）等操作后返回的矩陣也是稀疏的。

1.7 查詢屬性

??可以通過調用適當的成員函數來檢查矩陣或稀疏模式是否具有特定屬性。 例如查詢矩陣大小：

y = SX.sym('y',10,1) print(y.shape)

運行后輸出結果為：

(10, 1)

??查詢某個矩陣 A 的屬性常用的操作有下面這些：

函數名功能

A.size1()	行數
A.size2()	列數
A.shape	（在 MATLAB 中的“size”）形狀，即對 (nrow,ncol)
A.numel()	元素個數，即 nrow?ncol
A.nnz()	結構上非零元素的數量，如果密集，則等于 A.numel()。
A.sparsity()	檢索對稀疏模式的引用
A.is_dense()	是一個密集矩陣，即沒有結構零
A.is_scalar()	矩陣是否為標量，即維度為 1×1？
A.is_column()	矩陣是一個向量，即維度為 n 比 1？
A.is_square()	矩陣是方陣嗎？
A.is_triu()	是矩陣上三角嗎？
A.is_constant()	矩陣項都是常數嗎？
A.is_integer()	矩陣項都是整數值嗎？

---

1.8 線性代數

??CasADi 支持有限數量的線性代數運算，例如 對于線性方程組的解：

A = MX.sym('A',3,3) b = MX.sym('b',3) print(solve(A,b))

運行后輸出結果為：

(A\b)

1.9 微積分 - 算法微分

??CasADi 最核心的功能是algorithmic（or automatic）微分 (AD)。 對于函數 $f:{\mathbb{R}}^N\to {\mathbb{R}}^M$:
$$y=f(x),$$
前向模式方向導數（forward mode directional derivatives）可用于計算 Jacobian-times-vector 乘積：

$$\hat{y}=\frac{?f}{?x}\hat{x}$$
類似地，反向模式方向導數（reverse mode directional derivatives）可用于計算雅可比轉置時間向量乘積：
$$\hat{x}=(\frac{?f}{?x}{)}^T\hat{y}$$

??正向和反向模式方向導數的計算成本與評估 f(x) 成正比，而與 x 的維數無關。

??CasADi 還能夠有效地生成完整的、稀疏的雅可比行列式。 其算法非常復雜，但基本上包括以下步驟：

• 自動檢測 Jacobian 的稀疏模式
• 使用圖著色技術找到構建完整雅可比行列式所需的一些前向和/或方向導數
• 以數字或符號方式計算方向導數
• 計算雅可比行列式

Hessian 的計算方法是首先計算梯度，然后執行與上述相同的步驟，以與上述相同的方式計算梯度的雅可比矩陣，同時利用對稱性。

1.9.1 句法

??計算雅可比行列式的表達式：

A = SX.sym('A',3,2) x = SX.sym('x',2) print(jacobian(A@x,x))

運行后輸出結果為：

[[A_0, A_3],
?[A_1, A_4],
?[A_2, A_5]]

??當微分表達式為標量時，還可以計算矩陣意義上的梯度：

print(gradient(dot(A,A),A))

運行后輸出結果為：

[[(A_0+A_0), (A_3+A_3)],
?[(A_1+A_1), (A_4+A_4)],
?[(A_2+A_2), (A_5+A_5)]]

請注意，與 jacobian 不同，梯度始終返回一個密集向量。

Hessians 和作為副產品的梯度計算方法如下：

[H,g] = hessian(dot(x,x),x) print('H:', H)

運行后輸出結果為：

H: @1=2,
[[@1, 00],
?[00, @1]]

??對于計算雅可比乘法向量乘積，jtimes 函數——執行前向模式 AD——通常比創建完整雅可比矩陣和執行矩陣向量乘法更有效：

A = DM([[1,3],[4,7],[2,8]]) x = SX.sym('x',2) v = SX.sym('v',2) f = mtimes(A,x) print(jtimes(f,x,v))

運行后輸出結果為：

[(v_0+(3v_1)), ((4v_0)+(7v_1)), ((2v_0)+(8*v_1))]

jtimes 函數可選地計算轉置雅可比次向量乘積，即反向模式 AD：

w = SX.sym('w',3) f = mtimes(A,x) print(jtimes(f,x,w,True))

運行后輸出結果為：

[(((2w_2)+(4w_1))+w_0), (((8w_2)+(7w_1))+(3*w_0))]

總結

以上是生活随笔為你收集整理的CasADi——数据类型详解与基本操作介绍的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。