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

第7章非线性模型参数估值1第7章非线性模型参数估值7.1引言数学模型是观测对象各影响因素相互关系的定量描述。在获得实验数据并做了整理之后，就要建立数学模型。这一工作在科学研究中有着十分重要的意义。人们选用的模型函数可以是经验的，可以是半经验的，也可以是理论的。模型函数选定之后，需要对其中的参数进行估值并确定该估值的可靠程度。对于线性模型，待求参数可用线性最小二乘法求得，即用前一章中介绍的确定线性回归方程的方法。对于非线性模型，通常是通过线性化处理而化为线性模型，用线性最小二乘法求出新的参数，从而再还原为原参数。这种方法在处理经验模型时，简便易行，具有一定的实用价值。但要注意到，这样做是使变换后的新变量y的残差平方和(即剩余平方和)最小，这并不能保证做到使原变量y的残差平方和也达最小值。因此，得到的参数估计值就不一定是最佳的估计值。可见在求理论模型的参数时，这种线性化的方法尚有其不足之处。此外，还有些数学模型无法线性化，所以用线性化的方法是行不通的。为此，需要一种对非线性模型通用的(不管是经验模型还是理论模型，不管这个模型能否线性化)，能够得到参数最佳估计值的参数估计方法。在工程中，特别是在化学工程中的数学模型大多是非线性、多变量的。设yˆ为变量x1，x2，…，xp，的函数，含有m个参数b1，b2，…，bm，则非线性模型的一般形式可表示为:yˆf（x1，x2，…，xp；b1，b2，…，bm）(7.1)或写为),(ˆbxfy(7.2)式中x为p维自变量向量，b为m维参数向量。设给出n组观测数据x1，x2，…，xny1，y2，…，yn我们的目的是由此给出模型式(7.2)中的参数b的最佳估计值。可以证明，这个最佳估计值就是最小二乘估计值。按最小二乘法原理，b应使Q值为最小，即niiiyyQ12min)ˆ(或写成niiifyQ12min)],([bx(7.3)现在的问题是根据已知的数学模型和实验数据,求出使残差平方和最小,即目标函数式(7.3)取极小值时的模型参数向量b。这显然是一个最优化的数学问题,可以采用逐次逼近法求解。这种处理方法实质上是逐次线性化法或某种模式的搜索法。在下面各节中将介绍几个适用方法。第7章非线性模型参数估值27.2高斯——牛顿法高斯——牛顿法是非线性最小二乘法的最基本方法。其基本思想是把非线性模型函数在一局部范围内进行泰勒一阶展开，作为原函数的线性近似，代入目标函数则成为线性最小二乘法，因为它有解析解，所以可求出它的精确的极小值，以此点作为下一次线性近似的出发点。这样反复采用同样方法逐次逼近真正的极小点，从而得到最佳的估计参数向量b。所以此法实质是逐次线性化法。参数估计的目标函数为式(7.3)。求解参数向量b的思路是这样的，先给定初值b(0)，然后一次次修正)()()1(kkkΔbb(7.4)式中角标k代表迭代回次，Δ(k)代表第k次的修正量，每次修正须保证Q值下降，这样一步步逼近目标函数的极小值minQ，最后得到b的解。为确定Δ，对非线性函数，在初值点b(0)处进行泰勒一阶展开得:mjjjiiibfff1)0()0()0(Δ),(),(bxΔbx(7.5)或简记为mjjjiiibfff1)0()0((7.6)代入式(7.3)，得nimjjjiiibffyQ121)0()0(][(7.7)当b(0)给定后，)0(if和]/[)0(jibf都是自变量x的函数，根据实验点均可计算求得。因此目标函数化为对未知的Δ的线性函数，并可用线性最小二乘法求得Q(Δ)的极小点。由极值条件0)/(ΔQ，对式(7.7)求导数，并令其为零，得nimjkijjiiibfbffy11)0()0()0(0][2(k=1，2，…，m)由此得nimjnikiiijkijibffybfbf111)0()0()0()0()((k=1，2，…，m)(7.8)令mnnnmmbfbfbfbfbfbfbfbfbf)0(2)0(1)0()0(22)0(21)0(2)0(12)0(11)0(10A(7.9))0(iiifyr(7.10)],,,[21nTrrrr(7.11)第7章非线性模型参数估值3代入式(7.8)则有rAΔAATT000(7.12)令00AAAT(7.13)rADT0(7.14)则有线性方程组DAΔ(7.15)解之便得Δ，从而代入式(7.4)，得到修正的b。经过反复迭代直至Δ的值很小，满足误差要求为止，这时的b即为所求。以上计算过程归纳如下：(1)人为地给定初值b(0)；(2)求偏导值jibf/)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)表示为第7章非线性模型参数估值4)(bQQ(7.22)下面给出单纯形法的计算步骤及公式。1．初始单纯形的生成给定初始点Tmbbb],,,[210b(7.23)其余几个点为TmmTmTmpbqbqbqbpbqbqbqbpb],,,[],,,[],,,[21212211bbb(7.24)式中hmmqhmmmp211211(7.25)式中h为单纯形边长，一般取0.5≤h≤15。对于二维情况，有hqhpm2588.0,9659.0,2这样形成的初始单纯形为正单纯形。2．反射计算单纯形各顶点的目标函数值)(iiQQb(i=0，2，…，m)(7.26)比较这m+1个点的目标函数值的大小，先找出Q的最大、最小及次最大点的b值，分别记为Hb、Lb、Gb即)(max)(iHHQQQbb(i=0，1，…，m)(7.27))(min)(iLLQQQbb(i=0，1，…，m)(7.28))(max)(iGGQQQbb(i=0，1，…，m；i≠H)(7.29)计算除Hb外m个点的重心，即反射中心Cb，即miHiCm0)(1bbb(7.30)求出Hb的反射点Rb，即)(HCCRbbbb(7.31)式中α为反射系数，一般为对称反射，取α=1。3．扩张计算反射点的目标函数值)(RRQQb(7.32)第7章非线性模型参数估值5若RQGQ,则进行扩张。令)(HRHEbbbb(7.33)式中1，称为扩张系数，一般取=1.2～2.0。计算扩张点的目标函数值)(EEQQb(7.34)若EQRQ,则扩张成功。令ESESQQ,bb(7.35)否则扩张失败。令RSRSQQ,bb(7.36)总之，总可得一新点Sb，用Sb代替Hb构成新的单纯形返回第二步，再次搜索。4．压缩和收缩若RQ≥GQ，则反射失败，进行压缩。令)(HRHSbbbb(7.37)式中01和≠0.5，称为压缩系数，一般取=0.25或0.75，若=0.5，则Sb=Cb，造成降维。计算)(SSQQb(7.38)若SQGQ,则压缩成功。用Sb代替Hb构成新的单纯形返回第二步，再次搜索。若SQ≥GQ，则反射压缩都失败，则要进行收缩。收缩的公式为2/)()()()1(kLkikibbbi=0，1，…，m(7.39)式中上角标k表示计算ib的次数。以)1(kib代替ib构成新的单纯形返回第二步再次搜索。5．收敛要求继续上述过程，直至121)(mQQmLi0i(7.40)或220)(mQQmii(7.41)为止。式中1、2为充分小的正数。单纯形法需要的总搜索步数较多，但每一步计算中除计算一次目标函数值外，只是一些简单的代数运算。若使用计算机，运算时间就能大大缩短。单纯形法对初值的要求不是很严，即使所选用的初值远离真值也可以收敛。这一点是单第7章非线性模型参数估值6纯形法的突出优点。关于初始单纯形边长的选取，需根据具体问题由经验给出。在使用单纯形法进行参数估值时应注意下面几个问题：第一，估计参数开始时，单纯形的边长应与估计的参数具有相近的数量级，保证参数估计有较大的空间。第二，当一次估计找不到最佳参数值时，可把当前估计的参数值作为初值，接着再估计。单纯形的边长将随着参数接近最佳值而缩小。第三，如果有多个参数要同时估计，并且它们的数量级不同，最好对参数进行数量级转换，使估计的参数的数量级接近。例如估计下式)1(21xbeby中的参数1b和2b，若已知1b的数量级为210，而2b的数量级为110，估计前把上式改写成)1(100021xbeby即令13110bb，这样1b的数量级也变成为110，使1b和2b具有相同的数量级，可同时估计出1b和2b，然后再求出1b。这样，单纯形的边长和参数1b、2b可同步接近最佳值。7.4非线性模型参数估值的Excel“规划求解”法前面介绍的几种方法，当用到直接在工作表输入公式的方法编程时，如果对Excel软件了解不深，则是有一定难度的。而直接用Excel的工具——“规划求解”就简单多了。本小节就介绍这种方法。7.5.1基本原理在7.1节已经说过，非线性模型参数估值要解决的问题是：根据已知的数学模型和实验数据，求出目标函数式(7.3)最小时的模型参数向量b的值。解决这个问题的方法是采用逐次逼近的最优化方法求解。而用Excel工具——“规划求解”，可以直接解决这一问题。Excel的工具——“规划求解”的功能是，返回可变单元格的一组数据和目标单元格的值。利用这一功能解决非线性模型参数估值时，首先在Excel工作表中给定b的初值、目标函数的Excel公式，然后使用工具——“规划求解”，设置好其中的参数，点击“求解”，可以非常快速地求出满足目标函数最小时的模型参数向量b的值。7.5.2Excel工具——“规划求解”简介1．安装“规划求解”打开Excel工作簿后，在工具菜单中如果没有“规划求解”，此时可参照1.11节，安装“加载宏”。安装完毕后，即可使用“规划求解”工具了。此时操作：工具→“规划求解”，将显示“规划求解参数”对话框，再点击选项，又显示“规划求解选项”对话框，下面将介绍这两个对话框。第7章非线性模型参数估值72．“规划求解”参数对话框在这里设置“规划求解”的所有参数。(1)设置目标单元格在此指定目标单元格(放置目标函数所在单元格的引用)。经求解后，这个单元格将获得某一特定数值、最大值或最小值。这个单元格必须包含公式。(2)等于在此指定是否需要对目标单元格求取最大值、最小值或某一指定数值。如果需要指定数值，请在右侧“编辑”框中键入。对于非线性模型参数估值问题，在此指定最小值或指定“值为”选项并