多项式计算¶

数学背景¶

在能处理任意阶多项式之前，让我们首先思考一下我们如何使用函数计算一条直线上的点。在数学上，我们想定义一个如下所示的函数：

其中\(x\)是自变量，\((x_0,y_0)\)是定义直线的一个点，\((x_1,y_1)\)是定义直线的另一个点。数学上，这种函数可以被定义为如下形式：

为了减少参数数量，让我们假设\(x_0\)和\(y_0\)代表一个单一的点，以矢量\(\vec{p_0}\)表示，\(x_1\)和\(y_1\)代表另一个单一的点，以矢量\(\vec{p_1}\)表示，因此现在函数调用如如下所示：

Modelica描述¶

现在的问题是我们如何将这种数学关系转换为可以在Modelica模型中调用的函数。为了做到这一点，我们必须定义一个新的Modelica函数。

事实证明函数定义与模型定义模型定义非常相似（至少在语法层面）。这里我们展示的是在Modelica中定义的Line函数：

function Line "Compute coordinates along a line"
  input Real x     "Independent variable";
  input Real p0[2] "Coordinates for one point on the line";
  input Real p1[2] "Coordinates for another point on the line";
  output Real y    "Value of y at the specified x";
algorithm
  y := x*(p1[2]-p0[2])/(p1[1]-p0[1])+
       (p1[2]+p0[2]-(p1[1]+p0[1])*(p1[2]-p0[2])/
       (p1[1]-p0[1]))/2.0;
end Line;

函数的所有输入参数均以input限定词作为前缀。函数的输出结果以output限定词作为前缀。函数主体为algorithm区域。返回值（本例中是y）通过algorithm区域计算得到。

所以在本例中，output值y是根据input值x、p0、p1计算得到。注意，此函数中没有return语句。函数将自动返回algorithm区域计算得到的output变量的结果。

有几个在前面几章已经讨论过的事情需要注意。首先，注意函数本身和参数的说明字符串。它们在记录函数和参数的用法方面很有意义。同时需要注意如何使用数组来表达一个二维向量和这些数组在本例中是如何索引的。

Line模型唯一不足之处在于计算y表达式的长度。如果我们可以拆分表达式，则就能改善这种不足。

中间变量¶

为了简化y的表达式，我们需要引入一些中间变量。我们已经看到x，p0，p1是我们可以在函数中使用的变量。我们想引入额外的变量。但是，同时这些变量也不应是参数。换句话说，引入的变量值必须在函数“内部”进行计算。为了实现这个目标，我们创建了一个受保护（protected）变量的集合。这些变量可认为是在函数内部计算得到的。下面的例子展示了使用protected来声明和计算两个内部变量：

function LineWithProtected "The Line function with protected variables"
  input Real x     "Independent variable";
  input Real p0[2] "Coordinates for one point on the line";
  input Real p1[2] "Coordinates for another point on the line";
  output Real y    "Value of y at the specified x";
protected
  Real m = (p1[2]-p0[2])/(p1[1]-p0[1])        "Slope";
  Real b = (p1[2]+p0[2]-m*(p1[1]+p0[1]))/2.0  "Offset";
algorithm
  y := m*x+b;
end LineWithProtected;

此模型引入了两个新变量。其中一个变量m代表了直线的斜率。另一个变量b代表了当条件为x=0时的返回值。通过计算这两个中间变量，使得计算y的表达式变成了非常容易识别的形式：y := m*x+b。

数学定义¶

当然，本节我们的目标是创建一个函数，以此计算任意的多项式。到目前为止，我们已经看到了一些基本的函数。在此基础上让我们继续朝着最终目标前进。我们会构建一个函数，其调用方式如下：

其中，\(x\)依然代表独立变量。\(\vec{c}\)代表系数组成的向量。这样我们的多项式以如下方式计算：

其中，N是传递给函数的系数的个数。在这一点上有两件重要的事情需要注意。首先，\(\vec{c}\)的第一个元素对应于多项式中的最高阶次项。然后，我们通过使用符号并假设\(\vec{c}\)的元素是从1开始编号。以此使函数能更容易地转化为Modelica代码（其数组索引从1开始）。

注意，上文中\(p\)的定义非常容易阅读和理解。但是在处理有限精度的浮点数时，使用递归的方法来计算多项式会更高效和更精确。对于一个\(4\)阶多项式，其表达式为：

这种方法相对高效，因为它仅用到了简单的乘法和加法操作，避免了处理更复杂的求幂运算。同时，这种方法更为精确。因为使用有限精度浮点表示形式时，求幂运算更容易引发舍入误差和截尾误差。

Modelica定义¶

到目前为止，我们已经精确定义了我们希望函数执行的计算方法。剩下的工作只是在Modelica中定义函数。本例中，我们的多项式计算函数在Modelica中描述如下：

function Polynomial "Create a generic polynomial from coefficients"
  input Real x     "Independent variable";
  input Real c[:]  "Polynomial coefficients";
  output Real y    "Computed polynomial value";
protected
  Integer n = size(c,1);
algorithm
  y := c[1];
  for i in 2:n loop
    y := y*x + c[i];
  end for;
end Polynomial;

回顾一下，函数所有的（输入）参数都有input限定词，返回值都有output限定词。与前面的例子一样，我们定义了一个中间变量n，以此提供一种简便的方法去查阅系数向量的长度。本例中，我们也展示了如何使用一个for循环来描述任意阶数的多项式的递归计算。

让我们通过在模型中使用这个函数，以此验证其正确性。使用的Modelica模型如下：

model EvaluationTest1 "Model that evaluates a polynomial"
  Real yf;
  Real yp;
equation
  yf = Polynomial(time, {1, -2, 2});
  yp = time^2-2*time+2;
end EvaluationTest1;

记住，c的第一个元素对应最高阶项。如果我们把多项式的直接计算结果yp和使用我们函数计算的结果yf作对比，我们发现计算结果相等：

微分¶

这个多项式所计算的量完全可能最终被Modelica编译器求微分。下面的例子虽然有些牵强，但是它展示了一个多项式在什么情况下会在一个模型中被微分：

model Differentiation1 "Model that differentiates a function"
  Real yf;
  Real yp;
  Real d_yf;
  Real d_yp;
equation
  yf = Polynomial(time, {1, -2, 2});
  yp = time^2-2*time+2;
  d_yf = der(yf); // How to compute?
  d_yp = der(yp);
end Differentiation1;

本例中我们使用了与上述例子相同的yf和yp方程。yf使用Polynomial计算，yp直接使用多项式计算。我们还额外添加了两个变量，d_yf和d_yp。它们分别代表yf和yp的导数。如果我们尝试编译此模型，编译器很有可能会引发一个关于d_yf方程的错误。其原因是编译器无法计算yf的导数。这是因为与通过简单表达式计算的yp不同，我们在Polynomial函数的背后隐藏了计算yf的细节。通常来说，Modelica工具不会从函数实现推导其导数。即使编译器做到了这点，如何确定任意算法的导数依然很不容易。

因此，接下来的问题就是我们应如何处理这种情况？这不会使得在函数在模型中变得难以使用吗？幸运的是，Modelica给我们提供了一个方法来指定函数导数的计算。这就是通过在函数定义时添加一种叫做annotation的语句来实现的。

在本例中，我们需要derivative标注。因为此标注可以让我们告诉Modelica编译器如何计算函数的导数。为了实现这一点，我们定义了一个新函数，PolynomialWithDerivative用于计算，如下所示：

function PolynomialWithDerivative
  "Create a generic polynomial from coefficients (with derivative information)"
  input Real x     "Independent variable";
  input Real c[:]  "Polynomial coefficients";
  output Real y    "Computed polynomial value";
protected
  Integer n = size(c,1);
algorithm
  y := c[1];
  for i in 2:n loop
    y := y*x + c[i];
  end for;
  annotation(derivative=PolynomialFirstDerivative);
end PolynomialWithDerivative;

注意，这个函数除了高亮显示的行以外，与之前的函数是一致的。换句话说，我们只需要把这行代码添加到我们的函数中：

  annotation(derivative=PolynomialFirstDerivative);

我们通过在函数内添加这行代码以此向Modelica编译器解释如何计算这个函数的导数。这行代码的含义是函数PolynomialFirstDerivative应该用于计算PolynomialWithDerivative的导数。

在讨论函数PolynomialFirstDerivative的执行之前，我们首先要从数学角度明白什么是必需的。回忆我们之前定义的多项式插值函数：

注意\(p\)有两个参数。如果我们希望\(p\)对任意的变量\(z\)求导，我们可以使用链式法则来表示\(p\)对\(z\)的全微分，如下所示：

我们可以从最初\(p\)的定义中导出如下关系。首先，对于\(p\)对\(x\)的偏导数，我们可以得到：

其中\(c'\)定义如下：

其次，对于\(p\)对\(\vec{c}\)的偏导数，我们可以得到：

其中向量\(\vec{d_i}\)是\(NxN\)的单位矩阵的第\(i\)列。

事实表明，出于效率的原因，Modelica编译器提供的\(\frac{\mathrm{d}x}{\mathrm{d}z}\)和\(\frac{\mathrm{d}\vec{c}}{\mathrm{d}z}\)要好于提供函数去计算\(\frac{\partial p}{\partial x}\)和\(\frac{\partial p}{\partial c_i}\)的函数。因此从数学角度来说，Modelica编译器需要的是一个可以用以下参数进行调用的新函数：

因此：

为此，derivative标注必须指定一个与\(df\)参数相同的函数。在我们的例子中，函数PolynomialFirstDerivative将被定义为如下形式：

function PolynomialFirstDerivative
  "First derivative of the function Polynomial"
  input Real x;
  input Real c[:];
  input Real x_der;
  input Real c_der[size(c,1)];
  output Real y_der;
protected
  Integer n = size(c,1);
  Real c_diff[n-1] = {(n-i)*c[i] for i in 1:n-1};
algorithm
  y_der :=PolynomialWithDerivative(x, c_diff)*x_der +
          PolynomialWithDerivative(x, c_der);
end PolynomialFirstDerivative;

请注意我们最初函数的参数如何被复制变成（预期数量的）两倍那么多。第二组参数分别代表了\(\frac{\mathrm{d}x}{\mathrm{d}z}\)和\(\frac{\mathrm{d}\vec{c}}{\mathrm{d}z}\)这些量。注意，这里假设\(z\)是一个标量。因此输入参数的类型是相同的。利用我们关于多项式的偏导数的知识，导数的计算就可通过同一个多项式计算函数来完成。

我们可以通过以下的模型来使用所有的这些函数：

model Differentiation2 "Model that differentiates a function using derivative annotation"
  Real yf;
  Real yp;
  Real d_yf;
  Real d_yp;
equation
  yf = PolynomialWithDerivative(time, {1, -2, 2});
  yp = time^2-2*time+2;
  d_yf = der(yf);
  d_yp = der(yp);
end Differentiation2;

对这个模型进行仿真并且比较结果，我们看到参数yf和yp的对比结果，以及d_yf和d_yp的对比结果，如下：