CasADi - 最优控制开源 Python/MATLAB 库2-阿里云开发者社区

CasADi - 最优控制开源 Python/MATLAB 库1：https://developer.aliyun.com/article/1583976

3.8 线性代数

CasADi 支持数量有限的线性代数运算，例如线性方程组的求解：

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

(A\b)

3.9 微积分 - 算法微分

正向和反向模式方向导数的计算成本与计算 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

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 非线性求根问题

请考虑下面的方程组：

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：

对于常微分方程的求解器来说，第二个方程和代数变量 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 为例：

使用 "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

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

0.1724

4.4.2 灵敏度分析

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

4.5 非线性规划

注释

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

与 CasADi 一起发布或连接到 CasADi 的 NLP 求解器可以求解以下形式的参数化 NLP：

4.5.1 创建 NLP 求解器

NLP 解算器使用 CasADi 的 nlpsol 函数创建。不同的求解器和接口作为插件实现。请考虑以下形式的所谓罗森布洛克（Rosenbrock）问题：

使用 "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 高级接口

要通过高级界面解决这个问题，我们只需用 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：

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

CasADi - 最优控制开源 Python/MATLAB 库3：https://developer.aliyun.com/article/1584026

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