导函数使用说明

---

介绍模型中导函数的使用说明。

概述

在模型翻译分析阶段，Modelica工具采用符号变换技术实现将原始陈述式方程系统约减为可供求解器算法解算的数学模型。其中，符号求导（表达式求导）是影响模型数学解析式生成的重要因素，轻则造成后续求解精度差、计算效率低，重则直接导致模型翻译失败。由于符号求导是Modelica工具的内置操作，建模者对其的感知往往很微弱，既不清楚其应用场景，也不了解其影响域，因此难以在建模时作出针对性优化。通常情况下，Modelica工具能够很好的处理各类常规数学表达式的符号求导，然而当求导表达式为自定义函数或外部函数时，如果缺乏足够的信息，Modelica工具是无法执行这类表达式求导的。例如对于一个Modelica自定义函数，Modelica工具需要知道该函数的连续可导性，即是否关于某个输入参数可导以及几阶可导；对于一个外部函数，其函数实现对Modelica工具是未知的，更是无法对其执行符号求导。针对上述情况，Modelica规范定义了导函数相关内容，通过在模型自定义函数中添加导函数相关注解可帮助Modelica工具从模型中获取必要的求导信息。本文详细解读这部分导函数内容，帮助建模者更好的掌握导函数的使用技巧。

符号求导

符号求导是指对数学表达式中的符号进行求导操作，而不是对具体数值进行求导。例如一个代数表达式 关于变量 求导，结果为一个新的代数表达式 。符号求导是Modelica模型翻译过程中重要的符号变换技术手段，主要用在 DAE 方程降指标以及推导非线性方程组的解析雅克比。然而，模型中一些复杂或不可导的表达式可能导致符号求导失败，如下例所示：

model Case16_1

  function f1
    input Real u;
    output Real y;
  algorithm
    y := exp(u);
    y := y + sin(u);
  end f1;

  Real x;
  Real y;
equation
  x = ceil(time);		   // 情形一：ceil函数的导数不存在
  //x = f1(time + 1);	 // 情形二：自定义函数f1未显式提供导函数且未提供smoothOrder注解
  der(x) = y + 1;

end Case16_1;

如果在 DAE 降指标过程中符号求导失败，如上例所示，则模型分析失败并将错误信息打印至输出面板。

如果在推导非线性方程的解析雅克比时符号求导失败，则模型分析统计信息中将出现数值雅克比，即非线性方程组求解过程中将通过数值差分的方式计算雅克比矩阵，相较于解析雅克比，后者在计算精度以及求解效率等方面的表现都更差。如下例所示：

model Case16_2

  function f1	"y = u^2 + 2*u + 1"
    input Real u;
    output Real y;
  algorithm
    y := u * u;
    y := y + 2 * u + 1;
  end f1;

  Real x(start = 0.5);
  Real y(start = 0.5);
equation
  0 = f1(x);				      // 解析雅克比推导失败，采用数值雅克比
  y ^ 2 + 2 * y + 1 = 0;	// 解析雅克比推导成功
end Case16_2;

符号求导失败产生数值雅克比：

符号求导失败的几种常见情况

1）表达式中包含不可导的内置操作符

包括部分数学函数：div、mod、rem、ceil、floor；

包括部分返回实型值的特殊操作符：delay、spatialDistribution 。

2）表达式为自定义函数或外部函数

如果表达式为自定义函数或外部函数，符号求导需要推导该函数关于其输入变量的导函数，若用户建模时未显式提供导函数且未提供smoothOrder注解，则符号求导会失败，如上述Case16_1情形二。本文后续重点说明如何构建自定义函数的导函数，以及如何通过smoothOrder注解使Modelica工具自动构造导函数。

3）表达式在算法区或 when 方程中

若待求导的表达式位于算法区或为when方程，符号求导会失败。对于算法区的符号求导，可尝试将算法区内容放置于方程区，如下例所示，或将其中内容封装成自定义函数，并构造其导函数，详见后续章节内容。

model Case16_3
// DAE降指标失败
  Real x;
  Real y;
algorithm
  x :=  y - time;			// 将算法区内容置于方程区，或将其封装成自定义函数并构造其导函数
equation
  // x = y - time;			
  der(x) + der(y) = 0.1;
end Case16_3;

DAE 降指标失败：

model Case16_4
// 非线性方程组解析雅克比推导失败
  Real x;
  Real y;
algorithm
  x :=  y - time;			// 将算法区内容置于方程区，或将其封装成自定义函数并构造其导函数
equation
  // x = y - time;			
  x ^ 2 + y = 0.5 * time;
end Case16_4;

非线性方程组的部分内容位于算法区：

针对上述常见失败中的第二种情况，即当表达式为自定义函数或外部函数，用户在建模时可通过两种方式辅助实现函数的符号求导：

1）在自定义函数中添加smoothOrder注解

​ Modelica工具执行自定义函数求导时必须清楚该函数是否可导，而smoothOrder的作用就是标记该函数几阶可导。

2）声明导函数

​ 一种更直接的方式是建模者声明出导函数，并将其绑定给原函数，Modelica工具就可以直接使用无需额外符号求导。

smoothOrder

在声明Modelica自定义函数（非外部函数）时为其添加smoothOrder注解，则在模型翻译过程中Modelica工具将自动推导并构造出该自定义函数的导函数。注意，自定义函数中的复杂或不可导表达式会导致自动推导失败。

smoothOrder注解仅在自定义函数中有效，其使用方式为smoothOrder=n，其中n表示该自定义函数关于其所有实型输入变量均n阶连续可导。smoothOrder注解具有传递性，即由Modelica工具自动构造的导函数中仍具有smoothOrder注解，但其可导阶数变成n-1，高阶导函数以此类推。DAE 降指标时可能要求自定义函数的高阶导函数，因此给定smoothOrder阶数应符合实际，避免Modelica工具：

1）推导不存在的高阶导函数；

2）高阶导函数存在，但因smoothOrder=0而停止推导。

如下例所示，函数 f1 关于输入变量的任意阶可导，本例 DAE 降指标过程中仅需 f1 的 2 阶导函数，故 f1 中添加smoothOrder=2后成功翻译；若给定smoothOrder=1，则会导致模型翻译失败。

model Case16_5
  // Case16_1添加smoothOrder注解
  function f1
    input Real u;
    output Real y;
  algorithm
    y := exp(u);
    y := y + sin(u);
    annotation(smoothOrder=2);	// 添加smoothOrder注解，表明f1函数2阶连续可导
  end f1;
  
  Real x;
  Real y;
  Real z;
equation
  x = f1(time + 1);	           // DAE降指标时自动推导f1的各阶导函数
  der(x) = y + 1;
  der(y) = z;		               // 要得到der(y)，需获取f1的2阶导函数

end Case16_5;

降指标过程中执行的方程求导信息可通过勾选仿真设置-模型翻译-输出指标约减时的微分方程信息打印至输出面板。

导数方程信息：

​

下例演示自定义函数添加smoothOrder注解后成功推导非线性方程的解析雅克比：

model Case16_6
  // Case16_2添加smoothOrder注解
  function f1	"y = u^2 + 2*u + 1"
    input Real u;
    output Real y;
  algorithm
    y := u * u;
    y := y + 2 * u + 1;
    annotation(smoothOrder=1);    
  end f1;

  Real x(start = 0.5);
  Real y(start = 0.5);
equation
  0 = f1(x);				        // 解析雅克比推导成功
  y ^ 2 + 2 * y + 1 = 0;	  // 解析雅克比推导成功
end Case16_6;

Modelica工具自动构造导函数的过程中需要遍历原函数体，并递归地执行表达式求导，由于外部函数的具体实现对于Modelica工具是不可见的，因此 smoothOrder 不适用于外部函数。如果原函数体中调用了其它自定义函数，这些被调函数也必须提供导函数或添加smoothOrder注解，否则表达式求导将失败。

model Case16_7
  // Example005添加被调函数
  function f1
    input Real u;
    output Real y;
  algorithm
    y := exp(u);
    y := f2(y,u);				        // 调用自定义函数
    annotation(smoothOrder=2);	 // 添加smoothOrder注解，表明f1函数2阶连续可导
  end f1;
  
  function f2
  	input Real u1;
  	input Real u2;
  	output Real y;
  algorithm
  	y := u1;
  	y := y + sin(u2);
  	annotation(smoothOrder=2);	 // smoothOrder<2也会导致表达式求导失败
  end f2;
  
  Real x;
  Real y;
  Real z;
equation
  x = f1(time + 1);	            // DAE降指标时自动推导f1的各阶导函数
  der(x) = y + 1;
  der(y) = z;		                // 要得到der(y)，需获取f1的2阶导函数

end Case16_7;

​

导函数

基于smoothOrder自动构造导函数的方式存在一定的局限性，比如不能应用于外部函数、要求自定义函数对所有实型输入变量均满足smoothOrder阶连续可导。相比于Modelica工具自动构造导函数，一种更为直接的方式是由用户定义导函数，从而避免前者的局限性。为此，Modelica规范提供derivative注解用于自定义函数绑定其它自定义函数为其导函数。本节主要介绍两方面内容：

1）如何使用derivative注解绑定导函数；

2）如何声明导函数。

derivative

Modelica规范定义的derivative注解详见 MLS3.6 12.7 节。其语法结构如下图所示：

derivative 语法结构：

derivative注解包含两部分，分别为：

1）通过derivative=name绑定导函数；

2）通过derivative-constraints属性提供导函数限制。

• 绑定导函数

定义自定义函数时可将其它自定义函数声明为其导函数，具体做法是通过添加注解derivative=name，其中name为导函数的函数名，注意，每个自定义函数只能绑定一个导函数，且导函数名name必须是按导函数规则正确声明的自定义函数，导函数声明详见后续章节。如下例所示，通过在自定义函数 f1 中添加注解annotation(derivative=f2)，表示导数关系 。

model Case16_8

  function f1	"y = exp(u) + sin(u)"
    input Real u;
    output Real y;
  algorithm
    y := exp(u);
    y := y + sin(u)
    annotation(derivative=f2);	// 绑定导函数f2，即f1'=f2
  end f1;
  
  function f2	"dy = exp(u)*du + cos(u)*du"
  	input Real u;				       // 原始输入变量
  	input Real du;				      // 输入变量的导数
  	output Real dy;
  algorithm
  	dy := exp(u) * du;
    dy := dy + cos(u) * du;
  end f2;
  
  Real x;
  Real y;
equation
  x = f1(time + 1);	            // der(x)=der(f1(time+1))——>der(x)=f2(time+1,1)
  der(x) = y + 1;

end Case16_8;

• derivative-constraint

derivative-constraint提供了三种属性，用于限制导函数。

1）order

order表示当前绑定的自定义函数是几阶导函数，其用法为oder=n且 n>0，order默认值为 1，下例中 为一阶导函数， 为二阶导函数。导函数的阶数影响其函数声明的输入形参，详见导函数声明章节。

model Case16_9

  function f1	"y = exp(u) + sin(u)"
    input Real u;
    output Real y;
  algorithm
    y := exp(u);
    y := y + sin(u)
    annotation(derivative=f2);	         // 绑定导函数f2，即f1'=f2，默认order=1，因此f2为一阶导函数
  end f1;

  function f2	"dy = exp(u)*du + cos(u)*du"
  // 一阶导函数
  	input Real u;				                // 原始输入变量
  	input Real du;				               // 输入变量的导数
  	output Real dy;
  algorithm
  	dy := exp(u) * du;
    dy := dy + cos(u) * du;
    annotation(derivative(order=2)=f3);	// 绑定导函数f3，即f2'=f3，order=2，因此f3为二阶导函数
  end f2;
  
  function f3	"ddy = exp(u)*du^2 + exp(u)*ddu - sin(u)*du^2 + cos(u)*ddu"
  // 二阶导函数
  	input Real u;				                // 原始输入变量
  	input Real du;				               // 输入变量的一阶导数
  	input Real ddu;				              // 输入变量的二阶导数
  	output Real ddy;
  algorithm
  	ddy := exp(u) * du ^ 2 + exp(u) * ddu;
    ddy := ddy -sin(u) * du ^ 2 + cos(u) * ddu;
  end f3;

  Real x;
  Real y;
  Real z;
equation
  x = f1(time + 1);	 
  der(x) = y + 1;
  der(y) = z;
end Case16_9;

2）zeroDerivative

zeroDerivative 属性用于指定自定义函数调用时输入实参始终与自变量相互独立的输入参数，其用法为zeroDerivative ={，zeroDerivative=}， 表示自定义函数声明的输入参数名，如果有多个该类型输入参数则以逗号分隔。给定多元函数 ，多元函数求导遵循链式法则，其数学表达如下：

若输入参数中存在与自变量 相互独立的输入变量 ，则 ，导函数中表达式项 ，因此在定义导函数时可以忽略对该输入参数求导，同时导函数的参数列表中忽略该输入参数的导数，如下例所示。zeroDerivative属性主要用于忽略原函数输入参数中那些实型参量、非时变的记录变量（具有实型成员变量）等。

model Case16_10

  function f1	"y = p * exp(u) + sin(u)"
    input Real u;
    input Real p;				                             // 输入参数p表示外部传入的一个参量，因此与自变量t相互独立
    output Real y;
  algorithm
    y := p * exp(u);
    y := y + sin(u)
    annotation(derivative(zeroDerivative = p) = f2);	// 绑定导函数f2，指明f2的输入形参中可忽略dp
  end f1;

  function f2	"dy = p * exp(u)*du + cos(u)*du"
  	input Real u;				                             // 原始输入变量
  	input Real p;				                             // 原始输入变量，实际为参量，dp/dt=0，故忽略输入变量dp
  	input Real du;				                            // 输入变量的导数
  	output Real dy;
  algorithm
  	dy := p * exp(u) * du;		                        // 忽略dp/dt相关表达式
    dy := dy + cos(u) * du;
  end f2;

  parameter Real p = 1;
  Real x;
  Real y;
equation
  x = f1(time + 1, p);	 
  der(x) = y + 1;
end Case16_10;

​ 3）noDerivative

​ noDerivative 属性用于一种复合函数的求导优化。假设一个多元函数 被调用时，其中某个实参是关于其它参数的函数式，即形成复合函数 ，根据求导链式法则可知，建模时需要同时提供函数 和 的导函数。通过数学推导发现，可以通过定义一个函数等价为复合函数的导函数：

​ 即复合函数的导函数 ，这种形式可以看成是 ，其中 或不存在。因此，在声明自定义函数 时，添加导函数注解可使用noDerivative=y或zeroDerivative=y忽略导函数中的 dy 参数。注意， 并不是真的为 0 或不存在，而是通过数学变换隐含在了函数 中。还需注意的一点是，这种方式是针对复合函数 的优化，因此自定义函数 声明时绑定的导函数实际是复合函数的导函数，若模型中还存在非复合情况 的函数调用，且模型翻译时对该表达式进行求导，显然函数 绑定的导函数存在冲突，故该属性使用时需特别标注使用限制。

model Case16_11

  function f	"z = sin(x)^2 + y^2"
    input Real x;
    input Real y;
    output Real z;
  algorithm
    z := sin(x) * sin(x);
    z := z + y * y;
    annotation(derivative(noDerivative = y) = h);	 // 绑定导函数h，指明h的输入形参中可忽略dy
  end f;

  function g	"y = cos(x)+1"
    input Real x;
    output Real y;
  algorithm
    y := cos(x);
    y := y + 1;
  end f;	

  function h	"dz = (2*sin(x)*cos(x) - 2*y*sin(x))*dx"
  	input Real x;				                           // 原始输入变量
  	input Real y;				                           // 原始输入变量
  	input Real dx;				                          // 输入变量的导数
  	output Real dz;
  algorithm
  	dz := 2 * sin(x) * cos(x) - 2 * y * sin(x);		
    dz := dz * dx;
  end h;

  Real x;
  Real y;
  Real z;
equation
  x = time + 1;
  y = f(x, g(x));	                                  // 复合函数求导，无需提供函数g的导函数
  der(y) = z + 1;
end Case16_11;

声明导函数

原函数的输入输出参数可以包含任意简单类型（Real、Boolean、Integer、String、enumeration）或记录类型。求导操作仅对连续实型变量才有意义，这要求：

1）原函数的输入参数至少包含一个与自变量相关的实型变量；

2）原函数的输出参数至少包含一个实型变量。给定多元函数 ，根据多元函数链式法则求导可得：

由于输入变量是在函数调用时由外界传入，因此输入变量关于自变量的导数值对导函数而言是未知的。在声明导函数时，所有输入变量关于自变量的导数值也必须作为输入参数提供给导函数，各阶导函数的函数调用形式如下所示：

基于上述推导的函数调用形式，根据导函数的阶次不同，声明导函数时其输入与输出参数列表应按照如下规则：

声明一阶导函数：

1）首先声明原函数所有输入参数（如上述 ），然后按顺序声明上述Real类型输入参数的导数（如 ， 为布尔类型故无需添加导数变量）。注意，原函数与导函数中相同输入参数的名字、类型以及声明序必须保持一致；

2）按照原函数中输出参数的顺序声明Real类型输出参数的导数（如 ）。

声明 n>1 阶导函数：

1）在 n-1 阶导函数的输入参数基础上声明原函数实型输入参数的 n 阶导数（如 n=2 时增加声明 ），同样应保持声明序；

2）声明输出参数与 n-1 阶导函数相似，差异在于该输出参数实际表示的是 n 阶导数。

导函数的函数体则根据链式法则递归地对各表达式求导，如下例所示：

model Case16_12

  function f0	"原函数"
    input Real u;		           // 实型输入参数
    input Integer n;	          // 非实型输入参数
    output Real y;		          // 实型输出参数
  algorithm
    y := 0;
    for i in 1 : n loop
      y := y + u;
    end for;
    y := y + sin(u);
    annotation(derivative=f1);	// 绑定导函数f1
  end f0;

  function f1	"一阶导函数"
    input Real u;		           // 原函数实型输入参数
    input Integer n;	          // 原函数非实型输入参数
    input Real du;		          // 实型输入参数u的一阶导数
    output Real dy;	           // 实型输出参数y的一阶导数
  algorithm
    dy := 0;
    for i in 1 : n loop
      dy := dy + du;
    end for;
    dy := dy + cos(u) * du;
    annotation(derivative(order=2)=f2);
  end f1;
  
  function f2	"二阶导函数"
    input Real u;		           // 原函数实型输入参数
    input Integer n;	          // 原函数非实型输入参数
    input Real du;		          // 实型输入参数u的一阶导数
    input Real ddu;		         // 实型输入参数u的二阶导数
    output Real ddy;            // 实型输出参数y的二阶导数
  algorithm
    ddy := 0;
    for i in 1 : n loop
      ddy := ddy + ddu;
    end for;
    ddy := ddy - sin(u) * du ^ 2 + cos(u) * ddu;
  end f2;
  
  Real x;
  Real y;
  Real z;

equation
  x = f0(time + 1,5);
  der(x) = y + 1;
  der(y) = z;
end Case16_12;

假设原函数的输入参数中包含记录类型，且该记录中包含实型和非实型变量。根据导函数声明规范，导函数的输入参数中还必须包含该记录中实型变量的导数，这种情况下可声明一个新记录类型，该类型中仅包含原记录中的实型变量导数，如下例中记录类型ThermodynamicState与ThermoDynamicState_der ：

记录类型输入参数：

如下给出外部函数的导函数实现示例，外部函数的导函数可以是Modelica自定义函数，也可以是另一外部函数。

model Case16_13
  // 外部函数的导函数声明为自定义函数
  function f	"外部函数"
    input Real u1;
    input Real u2;
    output Real y;
  external "C"
    y = ExtFunc(u1,u2);
    annotation(derivative=f_der,	// 绑定导函数f_der
    Include="
    #ifdef __cplusplus
    extern \"C\" {
    #endif

    double ExtFunc(double u1, double u2)
    {
      return u1 * u2;
    }

    #ifdef __cplusplus
    }
    #endif
    ");
    
  end f;

  function f_der	"自定义导函数"
    input Real u1;
    input Real u2;
    input Real du1;
    input Real du2;
    output Real dy;
  algorithm
    dy := 0;
    dy := u2 * du1 + u1 * du2;
  end f_der;

  Real x,y;
  equation
  x = f(time, 1);
  der(x) = y;

end Case16_13;

model Case16_14
  // 外部函数的导函数声明为另一外部函数
  function f	"外部函数"
    input Real u1;
    input Real u2;
    output Real y;
  external "C"
    y = ExtFunc(u1,u2);
    annotation(derivative=f_der,	// 绑定导函数f_der
    Include="
    #ifdef __cplusplus
    extern \"C\" {
    #endif

    double ExtFunc(double u1, double u2)
    {
      return u1 * u2;
    }
    
    #ifdef __cplusplus
    }
    #endif
    ");

  end f;

  function f_der	"外部函数"
    input Real u1;
    input Real u2;
    input Real du1;
    input Real du2;
    output Real dy;
  external "C"
    dy = ExtFunc_der(u1,u2,du1,du2);
    annotation(    Include="
    #ifdef __cplusplus
    extern \"C\" {
    #endif

    double ExtFunc_der(double u1, double u2, double du1, double du2)
    {
      return u2 * du1 + u1 * du2;
    }
    
    #ifdef __cplusplus
    }
    #endif
    
    ");

  end f_der;

  Real x,y;
  equation
  x = f(time, 1);
  der(x) = y;

end Case16_14;

总结

本文对Modelica模型翻译过程中的符号求导进行了简要介绍，并进一步对Modelica规范中的导函数内容进行了详细解读。通过阅读本文，可进一步加深对下述内容的理解：

1）符号求导的应用场景

主要用在 DAE 方程降指标以及推导非线性方程组的解析雅克比。DAE 降指标过程中符号求导失败将直接导致模型翻译失败，输出面板的提示信息有助于建模者修改原模型。推导非线性方程组的解析雅克比失败后输出面板的模型统计信息中将出现数值雅克比。

2）符号求导失败的几种常见情况

3）函数表达式求导的前提条件是Modelica工具得确定该函数是可导的

4）Modelica自定义函数声明时添加smoothOrder注解表明该函数几阶可导，Modelica工具基于此信息自动构造其导函数

5）注解smoothOrder=n，n 的值要切合实际，避免Modelica工具做无谓的符号求导

6）smoothOrder的局限性

smoothOrder不适用于外部函数，并且要求该函数关于其所有实型输入参数均连续可导。

7）建模者可以自行实现导函数并通过derivative注解绑定给原函数

8）derivative-constraint属性的含义与使用方法

9）如何声明导函数，其参数列表应满足的特定规则