系统识别与模型参数估计优化

自然界中的水资源系统是一种十分复杂的系统，定量描述水资源系统的数学模型往往也比较复杂，甚至面临着这样或那样的问题。

不管采用什么方法来建立水资源系统模型，都将面临着一个共同问题，即自然界的水循环规律还不能被完全认识，获得的水资源知识还存在这样或那样的未知部分（包括结构或参数）。它们需要通过系统观测的数据信息加以判明或予以估计，其处理适当与否不仅直接影响到所建模型的质量，而且涉及到水资源系统特征的正确分析、水资源系统模型的实际应用及其评价等。为了解决这一问题，可以引用系统识别理论，即用系统分析的原理和方法建立可以定量描述水资源系统的数学模型，研究和处理验前缺乏或部分缺乏的系统信息，以及受各种随机性、模糊性等不确定性因素干扰的水资源系统。

1.系统识别的基本概念

条件：

目标：

要求：

式中：ε为数学模型与原型间的误差向量；E为模拟的目标函数；‖·‖为量度误差“距离”的某种范数；D0为数学模型的可行域。

应进一步说明的是：

（1）系统识别有三个基本成分，即模型类φ=｛·｝、观测资料信息D0=｛·｝和系统识别准则它们是描述识别问题的重要部分。只有当三者被确定之后，在数学上才进一步归结为模型或参数的最优化、最优估计等问题。

（2）模型类φ=｛·｝是由结构加参数或函数组成的，故识别内容可分为模型的结构识别与参数识别两部分。只有当结构在验前完全设定的情况下，问题才简化为参数识别。

（3）识别过程有三个环节，即确定基本成分（φ，D0，J）、完成数学求解和检验及应用。

2.系统识别的一般程序

（1）首先依据研究的问题（水文分析、预报，水库调度、水资源评价、规划、管理等），由水文学、水资源学原理和系统方法确定模型类。

（2）其次，选择适当的系统识别准则，利用可以观测到的原型信息（即水文、气象资料）去估计模型中未知部分，并做出必要的检验。

（3）最后应用于实际或反馈作再次改进。

3.系统识别的误差准则

原型观测值Y和模型计算值YM之间的误差通常用泛函数或范数表示：

E=E（Y，YM）

系统识别的误差准则通常是使误差泛函数（或范数）E取最小值。比如，最小二乘准则是

式中：ε（t）=Y（t）-YM（t）；Y、YM、ε均为定义在区间［0，T］上的函数。

4.系统识别最优估计方法

系统识别大多数涉及到数学上的模型最优化问题。关于求解模型最优化问题的方法有多种，其中常用的有最小二乘法、最速下降法等。

5.模型参数估计

如前所述，当模型结构完全确定之后，系统识别问题就转化为参数识别问题。即在输入、输出多组观测数据的基础上，运用系统识别最优估计方法，确定模型参数的工作，就称为模型参数估计。它的含义是：给定m组输入、输出，寻找1组参数，使数学模型的输出和真实物理系统的输出之间的误差最小。

在系统数学模型建立过程中，模型的参数估计是非常重要的一个环节。模型的可靠性与精确度如何，直接与参数估计的正确与否密切相关。模型参数估计的方法也有多种，如最小二乘法、步长加速法、惩罚函数法等。这里仅就常用的几种方法作简单介绍。

（1）矛盾线性方程组的最小二乘法：

在现行的各种水资源系统模型中，有一类较为特殊的模型“相对待定量为‘显式’的线性方程组”。如下式模型方程就是一个典型的线性方程模型：

Y=a1x1+a2x2+…+anxn

有些模型可以通过一定的转换，变成线性模型。比如，下式非线性模型方程：

Y=a0f1a1（x）f2

a2（x）…fnan（x）

通过取自然对数，就转化成线性方程：

lnY=lna0+a1lnf1（x）+a2lnf2（x）+…+anlnfn（x）

对于一般的线性模型：

yj=a1jx1+a2jx2+…+anjxn（j=1，2，…，m）

式中：xi为n个待定参变量（未知）；aij为m组输入观测信息（已知）；yj为m组输出观测信息（已知）。

线性方程组用矩阵表示：

其中

模型参数识别问题常归结为已知A和B求X，且m＞n。所以，它不符合由n个方程求n个未知数的定解条件。从m个已知量yi决定n个未知量x1，x2，…，xn来看，只要已知值的个数m大于惟一决定未知量所必需的方程个数，得到的是一个矛盾方程组。就这样的矛盾方程组，对任何n维向量X，总有

Am×nXn×1-Bm×1≠0（当m≠n）

令其差为残差ε（X），即(www.zuozong.com)

AX-B=ε（X）

式中：ε（X）=（ε1，ε2，…，εm）T。

一般说，有无穷多种ε（X）的取法都可以找到相应的X，使上式成为严格的等式。当然，我们总希望找到使ε在某种意义下达到最小的那个解X，这实际是处理矛盾方程组的最优化问题。如果我们寻找使残差ε1，ε2，…，εm的平方和最小的解X*，这类问题就是线性最小二乘问题。

设E（X）为目标函数，当寻求某一向量X*，使下式成立

E（X）=‖AX-B‖2=min

或

‖ε‖2=min

就称X*为原矛盾方程组的最小二乘解。

对于矛盾方程组AX=B，如果系数矩阵非奇异（矩阵行列式不为零），可求得惟一的最小二乘解：

X*=（ATA）-1（ATB）

（2）直接搜索的步长加速法：

这里介绍的步长加速法，是一种追寻目标函数脊线特征加速前进的直接搜索方法。该方法的特点是：在搜寻极小点的过程中采用的是试验方法，依据周围测试的若干点的情况，判断下山最快的路径。这类方法对目标函数的导数性质并不提要求，只要目标函数可以计算就行了。这种方法多半应用于二维或三维目标函数的极值问题，如

E（X，Y）=min

式中：X、Y为自变量（模型参数）。

步长加速法的基本思路是：依据一定的步长（ΔX，ΔY）作搜索移动。如果任意一组独立变量（X，Y）在已做过的移动试验中，是使目标函数值朝最优方向改进［E（X，Y）值减小］，那么这样的移动就还有继续进行的价值，我们将这种朝最优方向的改进看作为沿目标函数的谷线下降。如果在搜寻中认为谷线是直的，则沿谷线方向加大步长前进将是有效的。因此，步长加速法在搜寻中有两种类型的移动。一种是探测性的移动，目的是探求有利的方向；另一种是模式性的移动，它可以说是循前一步所探测的有利方向所做的一种加速。二者交替进行就构成步长加速法的搜索过程。在实际使用时，由于目标函数响应曲面的分布并不知道，一般在开始时仅按一定步长ΔX或ΔY试测，随着重复的见效逐步增大步长；如果遇到不良的试验，则缩短步长或退回原处；如果谷线的方向改变了，则重新按新的方向开始搜寻。当靠近E（X，Y）的极值点时，步长可以变得较小；以便获得改进的方向，搜索到规定的精度要求极值点为止。

将探测性移动的起点称为参考点，用R0，R1，…表示。探测性移动之后的终点称为基点，用B0，B1，…表示。两相邻的基点可构成一个矢量，沿该矢量的方向加速得到的便是一个新的参考点。模式性移动是以基点为起点，参考点移动后的那点为终点。两种移动交替进行，构成步长加速前进过程。一般的流程可简单地看做下列程序：

B0=R0→B1→R1→…

简单介绍计算步骤如下：

第一步：先确定搜索步长ΔX，ΔY及起始点（X0，Y0）。计算起始点的目标函数E（X0，Y0），并赋值

第五步：重复以上第二、第三、第四步。如果某一参考点经过探索所获得的点并不比原基点目标值小，则要返回到原来基点继续探索，重复以上各步。

第六步：若回到原基点后再经探索仍失败，那么该基点就是在所给定的步长ΔX，ΔY条件下寻找的最优解（Xm，Ym）。

在实际计算中，如果认为精度不够或希望求得更为精确的极值点，可以缩小ΔX，ΔY步长，以（Xm，Ym）为起始点，从第一步起继续搜寻下去，直至得到满意的解。

步长加速法在性质上属于无约束条件的测试寻优算法。它适用于模型参数之间无明显交互作用条件下的参数识别问题。

（3）有约束最优化问题的惩罚函数法：

在处理有约束最优化问题时，有一类常用的有效方法就是惩罚函数法（Penalty Function Method）。该方法的基本思路是：将一个有约束的最优化问题

转化为无约束的最优化问题来求解。采用的手段是在原目标函数E（X）上加一个或多个约束函数，并去掉这些约束条件，使原来有约束的极值问题变成一个无约束的极值问题：

称上式中函数minW（X，μ，γ）为惩罚函数，简称罚函数。H［hi（X）］和G［gj（X）］分别称为等式约束和不等式约束的惩罚项；μi和γj是惩罚项所连带的罚因子。

对于等式约束，典型的选择形式是当hi（X）→0时，有H［hi（X）］→0，一般取

为对于不等式约束，也有多种选择。其中，一种方法是取

式中：δj为单位步长函数，即定义为

如果按照上面对约束条件的处理形式，罚函数就变成

从上式中我们能分析到，罚函数的功能是对不满足约束条件的解起到某种制约作用，因为当解X不满足约束条件时，罚函数中的惩罚项非零，而为一个很大的正数，这时，W函数值可能要比E（X）数值大得多。显然，欲使W函数值取极小值，只有当惩罚项的数值很小时才有可能。这就迫使原来可能远离可行域D0的解变得越来越靠近D0。这种作用可以理解为是对不满足约束的解X实施的某种惩罚。相反，假如解满足约束条件，则惩罚项都为零，这时，W=E（X），表明函数E（X）不受任何惩罚。

按照这种方式，如果随着罚因子的分量逐渐加重，序列解X（k）收敛，则只需解无约束的优化问题，就能得到有约束的优化问题的解。一般计算步骤如下：

第一步：给定初始的近似点X（0）。一般取无约束函数E（X）的极值点。

第二步：判别X（0）是否满足可行域D0的约束条件。如果全部满足，则无约束函数E（X）的初始极值点X（0）就是所求的极值点。否则，如果不能全部满足，取μ（1）＞0，γ（1）＞0，形成惩罚函数。

第三步：利用无约束最优化方法求解转化后的罚函数式对应μ（1）、γ（1）的极小值X（1）。

第四步：判别X（1）是否满足可行域D0的约束条件。如果全部满足，则极值点X（1）就是所求的极值点，计算停止。如果不能全部满足，取μ（2）＞0，γ（2）＞0，形成惩罚函数。再重复第三步，得到新的极小值X（2）。

第五步：重复第四步，直至得到满意结果为止。