插值¶
本节中，我们将通过示例以不同的方法实现一个简单的一维插值法。我们首先会通过完全使用Modelica来实现一维插值法，然后展示一个Modelica结合C语言的替代方案。最后，将讨论这两种方案的优缺点。
Modelica实现¶
函数定义¶
在这个例子中，我们假定插值数据的格式如下：
自变量\(x\) |
因变量\(y\) |
| \(x_1\) | \(y_1\) |
| \(x_2\) | \(y_2\) |
| \(x_3\) | \(y_3\) |
| ... | ... |
| \(x_n\) | \(y_n\) |
这里我们假设\(x_i<x_{i+1}\)。
利用以上数据，对于我们某个感兴趣的自变量\(x\)值，我们的函数应对\(y\)赋一个内插值。这样的函数在Modelica中实现如下：
function InterpolateVector "Interpolate a function defined by a vector"
input Real x "Independent variable";
input Real ybar[:,2] "Interpolation data";
output Real y "Dependent variable";
protected
Integer i;
Integer n = size(ybar,1) "Number of interpolation points";
Real p;
algorithm
assert(x>=ybar[1,1], "Independent variable must be greater than "+String(ybar[1,1]));
assert(x<=ybar[n,1], "Independent variable must be less than "+String(ybar[n,1]));
i := 1;
while x>=ybar[i+1,1] loop
i := i + 1;
end while;
p := (x-ybar[i,1])/(ybar[i+1,1]-ybar[i,1]);
y := p*ybar[i+1,2]+(1-p)*ybar[i,2];
end InterpolateVector;
接下来，我们会一步步地对上述函数功能进行讲解。首先，从变量声明部分开始：
input Real x "Independent variable";
input Real ybar[:,2] "Interpolation data";
output Real y "Dependent variable";
input变量x表示插值函数求解的独立变量值。变量ybar表示插值数据。output变量y表示利用插值函数求得的内插值。函数包含的下一部分内容如下所示：
protected
Integer i;
Integer n = size(ybar,1) "Number of interpolation points";
Real p;
上述函数定义中也包含了各种protected变量的声明。在多项式计算示例中可以看到，函数内部使用了很多有效的中间变量。上述例子中，变量i表示索引变量，变量n表示插值数据点的个数，变量p表示插值函数使用的权重系数。
当完成所有的变量声明后，接下来需要定义的是函数的algorithm区域，如下所示：
algorithm
assert(x>=ybar[1,1], "Independent variable must be greater than "+String(ybar[1,1]));
assert(x<=ybar[n,1], "Independent variable must be less than "+String(ybar[n,1]));
i := 1;
while x>=ybar[i+1,1] loop
i := i + 1;
end while;
p := (x-ybar[i,1])/(ybar[i+1,1]-ybar[i,1]);
y := p*ybar[i+1,2]+(1-p)*ybar[i,2];
上述部分的前两个语句为assert声明，用于验证变量x的值是否在插值区间\([x_1, x_n]\)的范围内。如果变量不在上述区间，语句则会生成相应的错误信息。
函数的其余部分则用于查找满足表达式\(x_i<=x<x_{i+1}\)的i值。当变量i的值确定后，变量x对应的数值可以通过下述插值表达式进行计算：
其中
测试用例¶
现在，我们通过编写一个模型来测试上述函数。作为一个简单的测试用例，我们将对插值函数的返回值进行积分。我们将使用下列数据作为插值函数的基础数据：
| \(x\) | \(y\) |
| 0 | 0 |
| 2 | 0 |
| 4 | 2 |
| 6 | 0 |
| 8 | 0 |
如果绘制这些数据，所使用的插值函数如下图所示：
在下述测试模型中，我们令独立变量x等于变量time。插值函数根据样本数据计算输出变量y的值。然后通过积分y值，计算得到变量z的值。
model IntegrateInterpolatedVector "Exercises the InterpolateVector"
Real x;
Real y;
Real z;
equation
x = time;
y = InterpolateVector(x, [0.0, 0.0; 2.0, 0.0; 4.0, 2.0; 6.0, 0.0; 8.0, 0.0]);
der(z) = y;
annotation (experiment(StopTime=6));
end IntegrateInterpolatedVector;
上述模型的仿真结果如下图所示：
上述插值方法有几个缺点。首先，需要将数据传递到函数所有需求的地方。另外，对于多维插值函数，所需求的数据可能会非常复杂（例如不规则网格）和庞大。因此，以Modelica源代码的方式进行数据存储将会非常不方便，可以采用其他方式进行数据存储，例如使用外部文件。但是，如果使用外部文件进行插值数据的填充，需使用ExternalObject类型。
使用ExternalObject类型¶
ExternalObject类型主要是指定未（或不必）通过Modelica源代码体现的信息。这主要用于指定Modelica源代码以外的数据或状态。这些信息可以是插值函数所用的数据，也可以是其他软件的相关状态。
测试用例¶
在下述例子中，将围绕上面的测试用例进行讲述。在上述应用背景下，可以更清晰的展示如何使用ExternalObject类型。测试用例的Modelica源代码如下所示：
model IntegrateInterpolatedExternalVector
"Exercises the InterpolateExternalVector"
parameter VectorTable vector = VectorTable(ybar=[0.0, 0.0;
2.0, 0.0;
4.0, 2.0;
6.0, 0.0;
8.0, 0.0]);
Real x;
Real y;
Real z;
equation
x = time;
y = InterpolateExternalVector(x, vector);
der(z) = y;
annotation (experiment(StopTime=6));
end IntegrateInterpolatedExternalVector;
上述测试用例与之前测试用例的主要区别在于：插值数据并不是直接传递给插值函数，而是在模型中创建了一个VectorTable类型的变量vector，然后进行数据传递的。对于VectorTable类型，我们可以暂且把它们看作插值数据。在后面的章节中，我们将对其进行详细讨论。除了使用InterpolateExternalVector函数调用插值函数，以及创建vector变量，以替代之前的原始插值数据进行参数传递以外，上述模型的其他部分与之前的测试模型基本是相同的。
对上述测试模型进行仿真，其结果与之前模型的仿真结果是一致的，如下图所示：
ExternalObject类型定义¶
要明白上述模型的实现，我们首先需要了解如何定义VectorTable类型。如前所述，VectorTable是一种ExternalObject类型。在Modelica语言中，它用以代表一种特殊的类型。这通常被称为“隐形”指针。这也就意味着，ExternalObject类型主要表示那些不能直接（通过Modelica）访问的数据。
在本例里，我们通过下列方式定义VectorTable类型：
type VectorTable "A vector table implemented as an ExternalObject"
extends ExternalObject;
function constructor
input Real ybar[:,2];
output VectorTable table;
external "C" table=createVectorTable(ybar, size(ybar,1))
annotation(IncludeDirectory="modelica://ModelicaByExample.Functions.Interpolation/source",
Include="#include \"VectorTable.c\"");
end constructor;
function destructor "Release storage"
input VectorTable table;
external "C" destroyVectorTable(table)
annotation(IncludeDirectory="modelica://ModelicaByExample.Functions.Interpolation/source",
Include="#include \"VectorTable.c\"");
end destructor;
end VectorTable;
注意，VectorTable类型与ExternalObject类型是继承关系。在ExternalObject类型定义时，可以同时包含两个特殊的函数：constructor函数（构造函数）和destructor函数（析构函数）。下面将对其进行详细讲述。
构造函数¶
当创建VectorTable类型时，将调用constructor函数（例如在上述测试用例中，对vector变量的声明）。此时，constructor函数主要用于初始化模型中定义的隐形指针。初始化过程中的所有数据都应以参数的形式传递给constructor函数。在实例化过程中，相同的数据也必须存在（例如vector变量声明时定义的data参数）。
不同于之前的示例，constructor函数的定义有些特殊。它并不包括algorithm区域。通常，algorithm区域是用来计算函数返回值的。相反，constructor函数包含一个external声明，用以表明所定义的函数是通过Modelica以外的语言实现的。在上述示例中，实现的语言是C（在external关键字后以"C"标明）。而且table变量（函数的output，也同时是表示定义的隐形指针）其返回值（变量ybar的值及其数组长度）是通过C语言定义的createVectorTable函数实现的。
紧随createVectorTable函数其后的是annotation。此注解的功能是告诉Modelica编译器在哪里可以找到外部C函数的源代码。
至关重要的一点是，从Modelica编译器的角度来看，VectorTable变量只是createVectorTable函数返回的一个隐形指针。仅仅通过Modelica并不能访问这些指针指向的数据。在下面的例子中，我们可以看到，这些指针可以传递给其他由C语言实现的函数，并且可以以此访问变量VectorTable所代表的数据。
析构函数¶
当ExternalObject类型不再需要时，将调用destructor函数。这使得Modelica编译器可以清除ExternalObject类型运行所占用的内存。模型内ExternalObject类型的实例化通常会一直持续到仿真结束。但是，在函数定义时，如果ExternalObject类型被声明为protected变量，则可以在单个表达式的求值过程中创建和销毁。出于上述原因，我们必须确保由ExternalObject类型分配的任何内存在仿真结束后都被释放。
一般来说，destructor函数也是作为外部函数来实现的。在这种情况下，Modelica在调用destructor函数时会调用C语言定义的destroyVectorTable。上述C函数以VectorTable作为参数。当调用destructor函数时，所有与VectorTable实例相关的内存都将被释放。同样的，函数有着与构造函数相同的注释，以此通知Modelica编译器destoryVectorTable函数源代码的位置。
外部C代码¶
这些外部的C函数通过以下方式实现：
#ifndef _VECTOR_TABLE_C_
#define _VECTOR_TABLE_C_
#include <stdlib.h>
#include "ModelicaUtilities.h"
/*
Here we define the structure associated
with our ExternalObject type 'VectorTable'
*/
typedef struct {
double *x; /* Independent variable values */
double *y; /* Dependent variable values */
size_t npoints; /* Number of points in this data */
size_t lastIndex; /* Cached value of last index */
} VectorTable;
void *
createVectorTable(double *data, size_t np) {
VectorTable *table = (VectorTable*) malloc(sizeof(VectorTable));
if (table) {
/* Allocate memory for data */
table->x = (double*) malloc(sizeof(double)*np);
if (table->x) {
table->y = (double*) malloc(sizeof(double)*np);
if (table->y) {
/* Copy data into our local array */
size_t i;
for(i=0;i<np;i++) {
table->x[i] = data[2*i];
table->y[i] = data[2*i+1];
}
/* Initialize the rest of the table object */
table->npoints = np;
table->lastIndex = 0;
}
else {
free(table->x);
free(table);
table = NULL;
ModelicaError("Memory allocation error\n");
}
}
else {
free(table);
table = NULL;
ModelicaError("Memory allocation error\n");
}
}
else {
ModelicaError("Memory allocation error\n");
}
return table;
}
void
destroyVectorTable(void *object) {
VectorTable *table = (VectorTable *)object;
if (table==NULL) return;
free(table->x);
free(table->y);
free(table);
}
double
interpolateVectorTable(void *object, double x) {
VectorTable *table = (VectorTable *)object;
size_t i = table->lastIndex;
double p;
ModelicaFormatMessage("Request to compute value of y at %g\n", x);
if (x<table->x[0])
ModelicaFormatError("Requested value of x=%g is below the lower bound of %g\n",
x, table->x[0]);
if (x>table->x[table->npoints-1])
ModelicaFormatError("Requested value of x=%g is above the upper bound of %g\n",
x, table->x[table->npoints-1]);
while(i<table->npoints-1&&x>table->x[i+1]) i++;
while(i>0&&x<table->x[i]) i--;
p = (x-table->x[i])/(table->x[i+1]-table->x[i]);
table->lastIndex = i;
return p*table->y[i+1]+(1-p)*table->y[i];
}
#endif
这不是一本关于C语言编程的书。因此我们将不再详细讨论这段代码的语法及其所实现的功能。但是，我们可以对此文件的内容作出如下总结。
首先，名为VectorTable的struct变量内的数据对应Modelica变量VectorTable的内容。这个结构体中不仅仅包括插值数据（以成员变量x、y的形式）。它还包括数据点的个数npoints以及用于缓存最后使用的索引值lastIndex。
其次，可以看到createVectorTable函数对VectorTable的结构进行了分配并对所有的数据进行了初始化。上述实例化过程完成后返回到Modelica中运行。createVectorTable函数定义完成后，紧接着又对destroyVectorTable进行了定义，对createVectorTable函数定义的功能进行了有效的撤销。
最后，对interpolateVectorTable函数进行了定义。此函数是用C语言定义的，将VectorTable结构体及所需插值的独立变量进行参数传递，函数的返回值为利用插值算法计算的插入值。此函数实现的功能与之前定义的InterpolateVector函数所实现的功能几乎完全相同。在 Modelica 运行时，会执行ModelicaFormatError函数，以便报告 C 代码运行过程中出现的错误。对于interpolateVectorTable函数，这些错误函数主要用于实现前面InterpolateVector函数中的断言。在interpolateVectorTable函数中对于变量i的查找确认和Modelica版本基本相同。唯一的不同在于，查找并非每次都从1开始。函数会读取上次调用时找到的i值而开始查找。
插值¶
我们已经了解了如何定义interpolateVectorTable函数。但是目前为止，我们还不知道在什么地方调用了它。前面我们已经提到，除了使用VectorTable对象进行插值数据传递以外，其他的运行机理与InterpolateVector函数相同。在Modelica中调用interpolateVectorTable函数需定义一个简单的Modelica函数，如下所示：
function InterpolateExternalVector
"Interpolate a function defined by a vector using an ExternalObject"
input Real x;
input VectorTable table;
output Real y;
external "C" y = interpolateVectorTable(table, x)
annotation(IncludeDirectory="modelica://ModelicaByExample.Functions.Interpolation/source",
Include="#include \"VectorTable.c\"");
end InterpolateExternalVector;
在前面已经提到了，VectorTable表示隐形指针，而且Modelica代码不能访问VectorTable变量包含的数据。但是，在Modelica中定义的InterpolateExternalVector函数可以调用C语言定义的interpolateVectorTable函数，从而可以获取插值数据，根据插值算法计算相应的插入值。
讨论¶
如前所述，最初的插值方法需要传递大量复杂的插值数据。但是，通过定义VectorTable变量，可以非常方便的用一个变量表示上述插值数据。
关于ExternalObject插值方法，有一点在上述示例中并没有充分的展开讨论。但这点非常重要：初始化数据可以完全独立于Modelica源代码。简便起见，在上述示例中展示的代码中，通过定义数组的方式对VectorTable变量进行了初始化。但是，它也可以简单的通过文件名完成上述初始化工作。createVectorTable函数通过读取相应的文件，用文件中对应的数据对VectorTable结构体进行赋值，完成数据的初始化工作。很多情况下，这种方法不仅可以更方便的管理数据，而且可以利用C语言实现（新的或已有的）更加复杂的算法。
下一章里包含了另一个通过Modelica调用外部C代码的例子。