7.1  引言

    数学模型是观测对象各影响因素相互关系的定量描述。在获得实验数据并做了整理之后，就要建立数学模型。这一工作在科学研究中有着十分重要的意义。

    人们选用的模型函数可以是经验的，可以是半经验的，也可以是理论的。模型函数选定之后，需要对其中的参数进行估值并确定该估值的可靠程度。对于线性模型，待求参数可用线性最小二乘法求得，即用前一章中介绍的确定线性回归方程的方法。对于非线性模型，通常是通过线性化处理而化为线性模型，用线性最小二乘法求出新的参数，从而再还原为原参数。这种方法在处理经验模型时，简便易行，具有一定的实用价值。但要注意到，这样做是使变换后的新变量的残差平方和(即剩余平方和)最小，这并不能保证做到使原变量的残差平方和也达最小值。因此，得到的参数估计值就不一定是最佳的估计值。可见在求理论模型的参数时，这种线性化的方法尚有其不足之处。此外，还有些数学模型无法线性化，所以用线性化的方法是行不通的。为此，需要一种对非线性模型通用的(不管是经验模型还是理论模型，不管这个模型能否线性化)，能够得到参数最佳估计值的参数估计方法。

    在工程中，特别是在化学工程中的数学模型大多是非线性、多变量的。设为变量x1，x2，…，xp，的函数，含有个参数b1，b2，…，bm，则非线性模型的一般形式可表示为:

              f（x1，x2，…，xp ；b1，b2，…，bm ）                 (7.1)

或写为                     (7.2)

式中x为p维自变量向量，b为m维参数向量。

    设给出n组观测数据

              x1 ，x2 ，… ，xn

              y1 ，y2 ，… ，yn

    我们的目的是由此给出模型式(7.2)中的参数b的最佳估计值。可以证明，这个最佳估计值就是最小二乘估计值。

    按最小二乘法原理，b应使Q值为最小，即

或写成                      (7.3)

现在的问题是根据已知的数学模型和实验数据,求出使残差平方和最小,即目标函数式(7.3)取极小值时的模型参数向量b。这显然是一个最优化的数学问题,可以采用逐次逼近法求解。这种处理方法实质上是逐次线性化法或某种模式的搜索法。在下面各节中将介绍几个适用方法。

7.2  高斯——牛顿法

    高斯——牛顿法是非线性最小二乘法的最基本方法。其基本思想是把非线性模型函数在一局部范围内进行泰勒一阶展开，作为原函数的线性近似，代入目标函数则成为线性最小二乘法，因为它有解析解，所以可求出它的精确的极小值，以此点作为下一次线性近似的出发点。这样反复采用同样方法逐次逼近真正的极小点，从而得到最佳的估计参数向量b。所以此法实质是逐次线性化法。

    参数估计的目标函数为式(7.3)。求解参数向量b的思路是这样的，先给定初值b(0)，然后一次次修正

                     (7.4)

式中角标k代表迭代回次，Δ(k)代表第k次的修正量，每次修正须保证Q值下降，这样一步步逼近目标函数的极小值，最后得到b的解。

    为确定Δ，对非线性函数，在初值点b(0)处进行泰勒一阶展开得:

                             (7.5)

或简记为                 (7.6)

代入式(7.3)，得                     (7.7)

当b(0)给定后，和都是自变量x的函数，根据实验点均可计算求得。因此目标函数化为对未知的Δ的线性函数，并可用线性最小二乘法求得Q(Δ)的极小点。

    由极值条件，对式(7.7)求导数，并令其为零，得

                (k=1，2，…，m)   

由此得

             (k=1，2，…，m)   (7.8)

令

                                 (7.9)

                          (7.10)

                      (7.11)

代入式(7.8)则有 

                       (7.12)

令

                           (7.13)

                             (7.14)

则有线性方程组

                              (7.15)

解之便得Δ，从而代入式(7.4)，得到修正的b。经过反复迭代直至的值很小，满足误差要求为止，这时的b即为所求。

    以上计算过程归纳如下：

    (1)人为地给定初值b(0)；    

    (2)求偏导值，构成矩阵；

    (3)计算并构成向量；

    (4)计算并构成矩阵和向量；

    (5)求解线性方程组得到；

    (6)按式(7.4)修正；

    (7)若满足误差要求则结束，否则返回步骤(2)。

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

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

7.3  单纯形法及其Excel程序

高斯——牛顿法或麦夸特法都需要一阶导数矩阵，当函数关系复杂时，就会给运算带来很大的困难。这时可以采用不需要利用一阶导数矩阵的单纯形法。

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

    单纯形法的思路是先算出若干点处的目标函数值，然后进行比较，从它们之间的大小情况来判断函数的变化趋势，按照一定的模式确定搜索方向。这个模式就是利用单纯形进行反射、扩张、压缩和收缩等方法进行搜索，不断形成新的单纯形，直到单纯形缩得很小时，便得到最优值。

    那么，什么是单纯形呢？所谓单纯形就是一定的空间中的最简单的图形。如二维空间(平面上)，单纯形为三角形，三维空间为四面体，m维空间为由m +1个顶点而构成的最简单的图形。如果m +1个顶点的距离都相等，则称为正规的单纯形，简称正单纯形。

    为了讨论的方便，把式(7.3)表示为

                            (7.22)

    下面给出单纯形法的计算步骤及公式。

1．初始单纯形的生成

    给定初始点

                                          (7.23)

其余几个点为

                              (7.24)

式中

               (7.25)

式中h为单纯形边长，一般取0.5≤h≤15。对于二维情况，有

这样形成的初始单纯形为正单纯形。

2．反射

    计算单纯形各顶点的目标函数值

                      (i=0，2，…，m )                    (7.26)

比较这m+1个点的目标函数值的大小，先找出Q的最大、最小及次最大点的值，分别记为、、即

                       (i=0，1，…，m )      (7.27)

                        (i=0，1，…，m )      (7.28)

                    (i=0，1，…，m；i≠H )    (7.29)

计算除外m个点的重心，即反射中心，即

               (7.30)

求出的反射点，即

               (7.31)

式中α为反射系数，一般为对称反射，取α=1。

3．扩张

    计算反射点的目标函数值

                        (7.32)

若<,则进行扩张。令

                                            (7.33)

式中>1，称为扩张系数，一般取=1.2～2.0。

计算扩张点的目标函数值

                        (7.34)

若<,则扩张成功。令

                (7.35)

否则扩张失败。令

                (7.36)

    总之，总可得一新点，用代替构成新的单纯形返回第二步，再次搜索。

4．压缩和收缩

    若≥，则反射失败，进行压缩。令

                                            (7.37)

式中0<<1和≠0.5，称为压缩系数，一般取=0.25或0.75，若=0.5，则=，造成降维。计算

                        (7.38)

若<,则压缩成功。用代替构成新的单纯形返回第二步，再次搜索。

    若≥，则反射压缩都失败，则要进行收缩。收缩的公式为

                           i=0，1，…，m        (7.39)

式中上角标k表示计算的次数。以代替构成新的单纯形返回第二步再次搜索。

5．收敛要求

    继续上述过程，直至

                  (7.40)

或

                  (7.41)

为止。式中、为充分小的正数。

    单纯形法需要的总搜索步数较多，但每一步计算中除计算一次目标函数值外，只是一些简单的代数运算。若使用计算机，运算时间就能大大缩短。单纯形法对初值的要求不是很严，即使所选用的初值远离真值也可以收敛。这一点是单纯形法的突出优点。关于初始单纯形边长的选取，需根据具体问题由经验给出。

    在使用单纯形法进行参数估值时应注意下面几个问题：

    第一，估计参数开始时，单纯形的边长应与估计的参数具有相近的数量级，保证参数估计有较大的空间。

    第二，当一次估计找不到最佳参数值时，可把当前估计的参数值作为初值，接着再估计。单纯形的边长将随着参数接近最佳值而缩小。

    第三，如果有多个参数要同时估计，并且它们的数量级不同，最好对参数进行数量级转换，使估计的参数的数量级接近。例如估计下式

中的参数和，若已知的数量级为，而的数量级为，估计前把上式改写成

即令，这样的数量级也变成为，使和具有相同的数量级，可同时估计出和，然后再求出。这样，单纯形的边长和参数、可同步接近最佳值。

7.4  非线性模型参数估值的Excel“规划求解”法

前面介绍的几种方法，当用到直接在工作表输入公式的方法编程时，如果对Excel软件了解不深，则是有一定难度的。而直接用Excel的工具——“规划求解”就简单多了。本小节就介绍这种方法。

7.5.1  基本原理

在7.1节已经说过，非线性模型参数估值要解决的问题是：根据已知的数学模型和实验数据，求出目标函数式(7.3)最小时的模型参数向量b的值。解决这个问题的方法是采用逐次逼近的最优化方法求解。而用Excel工具——“规划求解”，可以直接解决这一问题。

Excel的工具——“规划求解”的功能是，返回可变单元格的一组数据和目标单元格的值。利用这一功能解决非线性模型参数估值时，首先在Excel工作表中给定b的初值、目标函数的Excel公式，然后使用工具——“规划求解”，设置好其中的参数，点击“求解”，可以非常快速地求出满足目标函数最小时的模型参数向量b的值。

7.5.2  Excel工具——“规划求解”简介

1．安装“规划求解”

打开Excel工作簿后，在工具菜单中如果没有“规划求解”，此时可参照1.11节，安装“加载宏”。安装完毕后，即可使用“规划求解”工具了。此时操作：工具→“规划求解”，将显示“规划求解参数”对话框，再点击选项，又显示“规划求解选项”对话框，下面将介绍这两个对话框。

2．“规划求解”参数对话框

在这里设置“规划求解”的所有参数。

(1)设置目标单元格

在此指定目标单元格(放置目标函数所在单元格的引用)。经求解后，这个单元格将获得某一特定数值、最大值或最小值。这个单元格必须包含公式。

(2)等于

在此指定是否需要对目标单元格求取最大值、最小值或某一指定数值。如果需要指定数值，请在右侧“编辑”框中键入。对于非线性模型参数估值问题，在此指定最小值或指定“值为”选项并在右侧“编辑”框中键入“0”。

(3)可变单元格

在此指定可变单元格(单个或多个单元格)。求解时对其中的数值不断进行调整，直到满足约束条件，并且在“设置目标单元格”编辑框中指定的单元格达到目标值。可变单元格必须直接或间接与目标单元格相联系。

(4)约束

在此列出了当前的所有约束条件。

(5)添加

显示“添加约束”对话框。

(6)更改

显示“改变约束”对话框。

(7)删除

删除选定的约束条件。

(8)求解

对定义好的问题进行求解。

(9)关闭

关闭对话框，不进行“规划求解”。但保留通过“选项”、“添加”、“更改”或“删除”按钮所做的修改。

(10)选项

显示“规划求解选项”对话框。在其中装入或保存“规划求解”模型，并对求解运算的高级属性进行设定。

(11)全部重设

清除“规划求解”中的当前设置，将所有的设置恢复为初始值。

3．“规划求解”选项对话框

在本对话框中，可以设定“规划求解”过程的一些高级属性、装入或保存“规划求解”定义，以及为线性或非线性规划设置参数。每一选项都有默认设置，可以满足大多数情况下的要求。

(1)最长运算时间

在此设定求解过程的时间。可输入的最大值为 32767(秒)，默认值为100(秒)，可以满足大多数小型“规划求解”要求。

(2)迭代次数

在此设定求解过程中迭代运算的最多次数，求解过程的时间。可输入的最大值为 32767，默认值为100次，可满足大多数小型“规划求解”要求。

(3)精度

在此输入用于控制求解精度的数字，以确定约束条件单元格中的数值是否满足目标值或上下限。精度值必须为小数(0 到 1 之间)，输入的数值越小，精度越高。例如 0.0001 比 0.01 精度高。

(4)允许误差

在此输入满足整数约束条件的目标单元格求解结果与最佳结果间的允许百分偏差。这个选项只应用于具有整数约束条件的问题。设置的允许误差值越大，求解过程就越快。

(5)收敛度

在此输入收敛度数值，当最近5次迭代时，目标单元格中数值的变化小于“收敛度”编辑框中设置的数值时，“规划求解”停止运行。收敛度只应用于非线性规划问题，并且必须由一个在 0(零) 和 1 之间的小数表示。设置的数值越小，收敛度就越高。例如，0.0001 表示比 0.01 更小的相对差别。收敛度越小，“规划求解”得到结果所需的时间就越长。

(6)采用线性模型

当模型中的所有关系都是线性的，并且希望解决线性优化问题时，选中此复选框可加速求解进程。

(7)显示迭代结果

如果选中此复选框，每进行一次迭代后都将中断“规划求解”，并显示当前的迭代结果。

(8)自动按比例缩放

当输入和输出值数量差别很大时，可以使用此功能。

(9)假定非负

对于在“添加约束”对话框的“约束值”编辑框中没有设置下限的可变单元格，假定其下限为 0(零)。

(10)估计

指定在每个一维搜索中用来得到基本变量初始估值的逼近方案。

正切函数，使用正切向量线性外推。

二次方程，用二次函数外推法，提高非线性规划问题的计算精度。

(10)导数

指定用于估计目标函数和约束函数偏导数的差分方案。

向前差分，用于大多数约束条件数值变化相对缓慢的问题。

中心差分，用于约束条件变化迅速，特别是接近限定值的问题。虽然此选项要求更多的计算,但在“规划求解”不能返回有效解时也许会有帮助。通常选定此选项,可以提高解的精度。

(11)搜索

指定每次的迭代算法，以确定搜索方向。

牛顿法，用准牛顿法迭代需要的内存比共轭法多，但所需的迭代次数少。

共轭法，比牛顿法需要的内存少，但要达到指定精度需要较多次的迭代运算。当问题较大或内存有限，或单步迭代进程缓慢时，用此选项。

7.5.3  程序33——“规划求解”法

打开一个新工作簿，以“程序33”为文件名存盘，

命名：

选中E7单元格，然后命名。

工具→规划求解，出现“规划求解”对话框，接着按下面的步骤操作：

(1)在“设置目标单元格”的编辑框中输入$E$4；

(2)在“等于”的右侧选中“最小值”；

(3)在“可变单元格”的编辑框中输入$D$3:$L$3，然后点击“选项”，出现“规划求解选项”对话框；

(4)在“导数”选项中选中“中心差分”；

(5)在“搜索”选项中选中“共轭法”后点击“确定”，这时又回到“规划求解选项”对话框；

(6)点击“求解”，很快显示“规划求解结果”对话框，点击“确定”即可。存盘结束。

显示结果：

见图7-5。

说明：

本程序在用于其他不多于5个参数的非线性模型参数估值时，先在数据区输入数据，接着按实际模型修改命名YG，然后操作：工具→规划求解，出现“规划求解”对话框，点击“求解”，出现“规划求解结果”对话框，点击“确定”即可。对于多于5个参数的非线性模型参数估值问题，可参照本程序稍加修改即可。

习题七

    1  已知数学模型和观测数据如下：

            1.00   1.00   1.00   2.72   2.72   2.72

            0.00   1.00   2.00   0.00   1.00   2.00

             0.00   3.00   12.00  2.00   5.00  14.00

要求≤0.01，初值。

     2  用下表中的实验数据，试用单纯形法对方程

中的参数和进行估值。取， =0.01。