第7章 非线性模型参数估值

第7章 非线性模型参数估值7.1 引言数学模型是观测对象各影响因素相互关系的定量描述。

在获得实验数据并做了整理之后，就要建立数学模型。

这一工作在科学研究中有着十分重要的意义。

人们选用的模型函数可以是经验的，可以是半经验的，也可以是理论的。

模型函数选定之后，需要对其中的参数进行估值并确定该估值的可靠程度。

对于线性模型，待求参数可用线性最小二乘法求得，即用前一章中介绍的确定线性回归方程的方法。

对于非线性模型，通常是通过线性化处理而化为线性模型，用线性最小二乘法求出新的参数，从而再还原为原参数。

这种方法在处理经验模型时，简便易行，具有一定的实用价值。

但要注意到，这样做是使变换后的新变量y '的残差平方和(即剩余平方和)最小，这并不能保证做到使原变量y 的残差平方和也达最小值。

因此，得到的参数估计值就不一定是最佳的估计值。

可见在求理论模型的参数时，这种线性化的方法尚有其不足之处。

此外，还有些数学模型无法线性化，所以用线性化的方法是行不通的。

为此，需要一种对非线性模型通用的(不管是经验模型还是理论模型，不管这个模型能否线性化)，能够得到参数最佳估计值的参数估计方法。

在工程中，特别是在化学工程中的数学模型大多是非线性、多变量的。

设y ˆ为变量x 1，x 2，…，x p ，的函数，含有m 个参数b 1，b 2，…，b m ，则非线性模型的一般形式可表示为:=y ˆf （x 1，x 2，…，x p ；b 1，b 2，…，b m ） (7.1) 或写为 ),(ˆb x f y= (7.2) 式中x 为p 维自变量向量，b 为m 维参数向量。

设给出n 组观测数据x 1 ，x 2 ，… ，x n y 1 ，y 2 ，… ，y n我们的目的是由此给出模型式(7.2)中的参数b 的最佳估计值。

可以证明，这个最佳估计值就是最小二乘估计值。

按最小二乘法原理，b 应使Q 值为最小，即∑==-=ni i i yy Q 12min )ˆ( 或写成 ∑==-=ni i i f y Q 12min )],([b x (7.3)现在的问题是根据已知的数学模型和实验数据,求出使残差平方和最小,即目标函数式(7.3)取极小值时的模型参数向量b 。

这显然是一个最优化的数学问题,可以采用逐次逼近法求解。

这种处理方法实质上是逐次线性化法或某种模式的搜索法。

在下面各节中将介绍几个适用方法。

7.2 高斯——牛顿法高斯——牛顿法是非线性最小二乘法的最基本方法。

其基本思想是把非线性模型函数在一局部范围内进行泰勒一阶展开，作为原函数的线性近似，代入目标函数则成为线性最小二乘法，因为它有解析解，所以可求出它的精确的极小值，以此点作为下一次线性近似的出发点。

这样反复采用同样方法逐次逼近真正的极小点，从而得到最佳的估计参数向量b 。

所以此法实质是逐次线性化法。

参数估计的目标函数为式(7.3)。

求解参数向量b 的思路是这样的，先给定初值b (0)，然后一次次修正)()()1(k k k Δb b +=+ (7.4) 式中角标k 代表迭代回次，Δ(k )代表第k 次的修正量，每次修正须保证Q 值下降，这样一步步逼近目标函数的极小值min Q ，最后得到b 的解。

为确定Δ，对非线性函数，在初值点b (0)处进行泰勒一阶展开得:∑=∂∂+≈+mj j j i i i b f f f 1)0()0()0(Δ),(),(b x Δb x (7.5)或简记为 ∑=∆∂∂+≈mj j ji ii b f f f 1)0()0( (7.6)代入式(7.3)，得 ∑∑==∆∂∂--=ni mj j ji ii b f f y Q 121)0()0(][ (7.7) 当b (0)给定后，)0(i f 和]/[)0(j i b f ∂∂都是自变量x 的函数，根据实验点均可计算求得。

因此目标函数化为对未知的Δ的线性函数，并可用线性最小二乘法求得Q (Δ)的极小点。

由极值条件0)/(=∂∂ΔQ ，对式(7.7)求导数，并令其为零，得∑∑===∂∂∆∂∂---n i m j k i j j ii i b f b f f y 11)0()0()0(0][2 (k=1，2，…，m)由此得∑∑∑===∂∂-=∆∂∂∂∂ni mj ni kii i j k i j i b f f y b f b f 111)0()0()0()0()( (k=1，2，…，m) (7.8) 令⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂⋯∂∂∂∂⋯⋯⋯⋯∂∂⋯∂∂∂∂∂∂⋯∂∂∂∂=m n n n m m b f b f b f b f b f b f b f b f b f )0(2)0(1)0()0(22)0(21)0(2)0(12)0(11)0(10A (7.9) )0(i i i f y r -= (7.10) ],,,[21n T r r r ⋯=r (7.11)代入式(7.8)则有r A ΔA A TT 000= (7.12) 令00A A A T = (7.13) r A D T 0= (7.14)则有线性方程组D A Δ= (7.15)解之便得Δ，从而代入式(7.4)，得到修正的b 。

经过反复迭代直至Δ的值很小，满足误差要求为止，这时的b 即为所求。

以上计算过程归纳如下： (1)人为地给定初值b (0)；(2)求偏导值j i b f ∂∂/)0(，构成矩阵0A ； (3)计算并构成向量r ；(4)计算并构成矩阵A 和向量D ； (5)求解线性方程组得到Δ； (6)按式(7.4)修正b ；(7)若Δ满足误差要求则结束，否则返回步骤(2)。

在求解过程中之所以需要反复地迭代和修正，是因为泰勒级数展开式(7.6)只是近似的式子，故得到的b 也是近似的。

高斯——牛顿法对初值的要求是比较严格的，若初值选取不当，很可能得不到结果，即“发散”。

7.3 单纯形法及其Excel 程序高斯——牛顿法或麦夸特法都需要一阶导数矩阵0A ，当函数关系复杂时，就会给运算带来很大的困难。

这时可以采用不需要利用一阶导数矩阵的单纯形法。

我们的目的是寻求一组参数b 使目标函数Q 值为最小，寻求到的b 值称为最优估计值，简称最优值。

单纯形法的思路是先算出若干点处的目标函数值，然后进行比较，从它们之间的大小情况来判断函数的变化趋势，按照一定的模式确定搜索方向。

这个模式就是利用单纯形进行反射、扩张、压缩和收缩等方法进行搜索，不断形成新的单纯形，直到单纯形缩得很小时，便得到最优值。

那么，什么是单纯形呢？所谓单纯形就是一定的空间中的最简单的图形。

如二维空间(平面上)，单纯形为三角形，三维空间为四面体，m 维空间为由m +1个顶点而构成的最简单的图形。

如果m +1个顶点的距离都相等，则称为正规的单纯形，简称正单纯形。

为了讨论的方便，把式(7.3)表示为)(b Q Q = (7.22) 下面给出单纯形法的计算步骤及公式。

1．初始单纯形的生成给定初始点T m b b b ],,,[210⋯=b (7.23)其余几个点为⎪⎪⎭⎪⎪⎬⎫+⋯++=⋯⋯⋯+⋯++=+⋯++=T m m T m T m p b q b q b q b p b q b q b q b p b ],,,[],,,[],,,[21212211b b b (7.24)式中⎪⎪⎭⎪⎪⎬⎫-+=-++=h m m q h m m m p 211211 (7.25)式中h 为单纯形边长，一般取0.5≤h ≤15。

对于二维情况，有h q h p m 2588.0,9659.0,2=== 这样形成的初始单纯形为正单纯形。

2．反射计算单纯形各顶点的目标函数值)(i i Q Q b = (i=0，2，…，m ) (7.26) 比较这m +1个点的目标函数值的大小，先找出Q 的最大、最小及次最大点的b 值，分别记为H b 、L b 、G b 即)(max )(i H H Q Q Q b b == (i =0，1，…，m ) (7.27) )(min )(i L L Q Q Q b b == (i =0，1，…，m ) (7.28) )(max )(i G G Q Q Q b b == (i =0，1，…，m ；i ≠H ) (7.29) 计算除H b 外m 个点的重心，即反射中心C b ，即∑=-=mi H i C m 0)(1b b b (7.30)求出H b 的反射点R b ，即)(H C C R b b b b -+=α(7.31) 式中α为反射系数，一般为对称反射，取α=1。

3．扩张计算反射点的目标函数值)(R R Q Q b = (7.32)若R Q <G Q ,则进行扩张。

令)(H R H E b b b b -+=μ (7.33) 式中μ>1，称为扩张系数，一般取μ=1.2～2.0。

计算扩张点的目标函数值)(E E Q Q b = (7.34) 若E Q <R Q ,则扩张成功。

令E S E S Q Q ==,b b (7.35) 否则扩张失败。

令R S R S Q Q ==,b b (7.36)总之，总可得一新点S b ，用S b 代替H b 构成新的单纯形返回第二步，再次搜索。

4．压缩和收缩若R Q ≥G Q ，则反射失败，进行压缩。

令)(H R H S b b b b -+=λ (7.37) 式中0<λ<1和λ≠0.5，λ称为压缩系数，一般取λ=0.25或0.75，若λ=0.5，则S b =C b ，造成降维。

计算)(S S Q Q b = (7.38)若S Q <G Q ,则压缩成功。

用S b 代替H b 构成新的单纯形返回第二步，再次搜索。

若S Q ≥G Q ，则反射压缩都失败，则要进行收缩。

收缩的公式为2/)()()()1(k L k i k i b b b +=+ i =0，1，…，m (7.39)式中上角标k 表示计算i b 的次数。

以)1(+k i b 代替i b 构成新的单纯形返回第二步再次搜索。

5．收敛要求 继续上述过程，直至121)(ε<+-∑=m Q Q mL ii (7.40)或22)(ε<-∑=mQ Q mii (7.41)为止。

式中1ε、2ε为充分小的正数。

单纯形法需要的总搜索步数较多，但每一步计算中除计算一次目标函数值外，只是一些简单的代数运算。

若使用计算机，运算时间就能大大缩短。

单纯形法对初值的要求不是很严，即使所选用的初值远离真值也可以收敛。

这一点是单纯形法的突出优点。

关于初始单纯形边长的选取，需根据具体问题由经验给出。

在使用单纯形法进行参数估值时应注意下面几个问题：第一，估计参数开始时，单纯形的边长应与估计的参数具有相近的数量级，保证参数估计有较大的空间。