留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

不等式约束PEIV模型的最优性条件及SQP算法

谢建 龙四春 周璀

谢建, 龙四春, 周璀. 不等式约束PEIV模型的最优性条件及SQP算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
引用本文: 谢建, 龙四春, 周璀. 不等式约束PEIV模型的最优性条件及SQP算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
XIE Jian, LONG Sichun, ZHOU Cui. Optimality Conditions of Inequality Constrained Partial EIV Model and the SQP Algorithm[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
Citation: XIE Jian, LONG Sichun, ZHOU Cui. Optimality Conditions of Inequality Constrained Partial EIV Model and the SQP Algorithm[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297

不等式约束PEIV模型的最优性条件及SQP算法

doi: 10.13203/j.whugis20180297
基金项目: 

国家自然科学基金 41704007

国家自然科学基金 41474014

国家自然科学基金 41604012

湖南科技大学博士启动基金 E51673

详细信息
    作者简介:

    谢建,博士,主要从事测量平差及数据处理研究。hsiejian841006@163.com

  • 中图分类号: P207

Optimality Conditions of Inequality Constrained Partial EIV Model and the SQP Algorithm

Funds: 

The National Natural Science Foundation of China 41704007

The National Natural Science Foundation of China 41474014

The National Natural Science Foundation of China 41604012

the Doctoral Scientific Foundation of Hunan University of Science and Technology E51673

More Information
    Author Bio:

    XIE Jian, PhD, specializes in surveying adjustment and data processing. E-mail:hsiejian841006@163.com

计量
  • 文章访问数:  815
  • HTML全文浏览量:  173
  • PDF下载量:  31
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-13
  • 刊出日期:  2020-07-30

不等式约束PEIV模型的最优性条件及SQP算法

doi: 10.13203/j.whugis20180297
    基金项目:

    国家自然科学基金 41704007

    国家自然科学基金 41474014

    国家自然科学基金 41604012

    湖南科技大学博士启动基金 E51673

    作者简介:

    谢建,博士,主要从事测量平差及数据处理研究。hsiejian841006@163.com

  • 中图分类号: P207

摘要: 基于约束非线性规划理论的最优性条件,推导了不等式约束PEIV(partial errors-in-variables)模型在加权最小二乘准则下取得最优解的一阶必要条件和二阶充分条件,以此作为算法设计的依据和检核解最优性的标准。根据序列二次规划算法,将非线性目标函数和约束方程在近似值处用泰勒级数展开,转换为二次规划子问题,采用积极约束算法同时估计模型参数和系数阵元素。数值模拟算例和线性回归的结果表明,新算法可行有效,具有良好的计算效率。

English Abstract

谢建, 龙四春, 周璀. 不等式约束PEIV模型的最优性条件及SQP算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
引用本文: 谢建, 龙四春, 周璀. 不等式约束PEIV模型的最优性条件及SQP算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
XIE Jian, LONG Sichun, ZHOU Cui. Optimality Conditions of Inequality Constrained Partial EIV Model and the SQP Algorithm[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
Citation: XIE Jian, LONG Sichun, ZHOU Cui. Optimality Conditions of Inequality Constrained Partial EIV Model and the SQP Algorithm[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
  • 在大地测量坐标转换[1-2]、曲线曲面拟合[3]、变形监测[4-5]和摄影测量[6]等领域,变量含误差(errors-in-variables,EIV)模型及相应的加权整体最小二乘理论(weighted total least squares, WTLS)由于能考虑设计阵的随机特性,成为数十年来的研究热点。EIV模型的常用算法有奇异值分解类算法、拉格朗日乘子迭代法、线性近似方法[7-10]。Xu等[11]将EIV模型扩展成部分EIV(partial EIV, PEIV)模型,它能反映系数矩阵的结构,当随机量较少时,待估参数也随之减少,推导了PEIV模型的算法和近似方差。Shi等[12]给出了PEIV模型解的经典最小二乘形式,并分析了算法的复杂度。将测量实际中的可靠不等式约束纳入PEIV模型,就形成了不等式约束的PEIV(inequality constrained PEIV, ICPEIV)模型及不等式约束整体最小二乘(inequality constrained WTLS, ICWTLS)理论[13-14]。De Moor[15]给出了不等式约束TLS问题的KKT(Karush-Kuhn-Tucker)条件,转化为广义线性补(linear complementary programming, LCP)问题, 采用组合方法求解,但仅考虑了非负约束与有界约束情形。Zhang等[16]采用穷举搜索法,计算3所有等式约束EIV模型的解,选单位权方差最小的解作为最终解,其计算量随约束数的增加呈指数级增长。以上研究均未考虑不等权和相关观测。Fang[17]提出了一般权阵下的ICWTLS问题, 推导了最优解的一阶必要条件和二阶充分条件,采用有效集(active set, AS)法和序列二次规划算法(sequential quadratic programming, SQP)求ICWTLS解。Zeng等[14]和De Moor[15]首先采用惩罚函数法求解ICPEIV模型,然后研究了不等式约束同时包含模型参数和系数阵元素的情况,采用线性近似方法,将ICWTLS问题化为不等式约束最小二乘(inequality constrained least squares, ICLS)模型,采用LCP方法求解,并分析了解的近似精度。

    目前针对ICPEIV模型的研究,尚未给出其获得最优解的条件,已有线性化方法的计算效率有待提高。针对上述问题,本文研究了ICPEIV模型获得最优解的条件,提出了ICPEIV模型的SQP算法,并与线性化方法的计算效率对比分析,验证了算法的有效性和可行性。

    • ICPEIV模型可以表示成[13]

      $$ {\mathit{\boldsymbol{y}} = \left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) + {\mathit{\boldsymbol{e}}_y}} $$ (1a)
      $$ {\mathit{\boldsymbol{a}} = \mathit{\boldsymbol{\bar a}} + {\mathit{\boldsymbol{e}}_a}} $$ (1b)
      $$ {\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} \ge 0} $$ (1c)

      式中,$ \mathit{\boldsymbol{y}}$和${{\mathit{\boldsymbol{e}}_y}} $分别为$n $维观测向量和误差; $\mathit{\boldsymbol{\beta }} $为$ m$维模型参数; $\mathit{\boldsymbol{a}} $是系数矩阵中随机元素所构成的$t$维列向量; $ {\mathit{\boldsymbol{\bar a}}}$和$ {{\mathit{\boldsymbol{e}}_a}}$分别是相应的真值和误差向量; 将系数矩阵中的非随机元素保留,将随机元素以0值替代后,组成的常数矩阵向量化,得到$ nm$维常数向量$\mathit{\boldsymbol{h}} $; $\mathit{\boldsymbol{B}} $是$ n m \times t$维常数矩阵,其形式由系数矩阵中非随机元素的个数及元素间的相关性决定; $ {{\mathit{\boldsymbol{I}}_n}}$是$ n$维单位矩阵; $ \mathit{\boldsymbol{G}}$为$ s \times m$维约束矩阵; $ \mathit{\boldsymbol{d}}$为常数向量。相应的误差随机模型为:

      $$ \begin{array}{c} \boldsymbol{D}\left(\boldsymbol{e}_{y}\right)=\sigma_{0}^{2} \boldsymbol{W}^{-1}, \boldsymbol{D}\left(\boldsymbol{e}_{a}\right)=\boldsymbol{\sigma}_{0}^{2} \boldsymbol{\omega}^{-1} \end{array} $$ (1d)
      $$ \operatorname{cov}\left(\boldsymbol{e}_{y}, \boldsymbol{e}_{a}\right)=0 $$ (1e)

      式中,$\boldsymbol{\sigma}_{0}^{2} $为单位权方差; $\boldsymbol{W} $为$ n$维观测值权矩阵; $ \boldsymbol{\omega} $为$t $维随机误差权矩阵。在加权最小二乘准则下,可以得到如下的约束极小化目标函数:

      $$ \left\{ {\begin{array}{*{20}{c}} {\min :\mathit{\Phi} (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }}) = {{(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}})}^{\rm{T}}}\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}}) + }\\ {{{\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right)}^{\rm{T}}} \times }\\ {\mathit{\boldsymbol{W}}\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right)}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\quad \mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} \ge 0} \end{array}} \right. $$ (2)
    • 对于约束非线性极小化问题式(2),为了获得最优解的充要条件,一般是在不等式约束中引入含$ s$个非负元素的基变量$\mathit{\boldsymbol{\eta }} $,使下式成立[18]

      $$ \mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} - \mathit{\boldsymbol{\eta }} = 0 $$ (3)

      联合式(2)和式(3),构造Lagrange目标函数为:

      $$ \begin{array}{l} \psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }}) = {(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}})^{\rm{T}}}\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}}) + \\ {\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right)^{\rm{T}}} \times \\ \mathit{\boldsymbol{W}}\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right) - \\ {\mathit{\boldsymbol{\lambda }}^{\rm{T}}}(\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} - \mathit{\boldsymbol{\eta }}) \end{array} $$ (4)

      式中,$\mathit{\boldsymbol{\lambda }} $为$s $维Lagrange乘子向量。目标函数式(4)分别对未知向量$ {\mathit{\boldsymbol{\bar a}}}$、$ \mathit{\boldsymbol{\beta }}$、$ \mathit{\boldsymbol{\lambda }}$求偏导数,并令其值为0,可以得到:

      $$ {\frac{{\partial \psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}}} = 2\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{{\hat {\bar{ a}}}}} - \mathit{\boldsymbol{a}}) + 2{\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat \beta }} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{W}} \times }\\ {\left( {\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B{\hat {\bar{ a}}}}}) - \mathit{\boldsymbol{y}}} \right) = 0} $$ (5a)
      $$ {\frac{{\partial \psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}}} = 2\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{h}}_1^{\rm{T}} + {{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_1^{\rm{T}}}\\ {\mathit{\boldsymbol{h}}_2^{\rm{T}} + {{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_2^{\rm{T}}}\\ \vdots \\ {\mathit{\boldsymbol{h}}_m^{\rm{T}} + {{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_m^{\rm{T}}} \end{array}} \right]\mathit{\boldsymbol{W}} \times }\\ \left\{ {\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B{\hat {\bar{ a}}}}}) - \mathit{\boldsymbol{y}}} \right\} - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} = 0 $$ (5b)
      $$ \frac{{\partial \mathit{\boldsymbol{\psi }}(\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\lambda }}}} = \mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} + \mathit{\boldsymbol{\eta }} = 0 $$ (5c)

      式中, $\mathit{\boldsymbol{h}} = {\left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{h}}_1^{\rm{T}}}&{\mathit{\boldsymbol{h}}_2^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{h}}_m^{\rm{T}}} \end{array}} \right]^{\rm{T}}}, \mathit{\boldsymbol{B}} = {\left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{B}}_1^{\rm{T}}}&{\mathit{\boldsymbol{B}}_2^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{B}}_m^{\rm{T}}} \end{array}} \right]^{\rm{T}}} $, 具体表达式见文献[11]。根据最优化理论,ICPEIV模型在最优解处的Kuhn-Tucker (KT)必要条件为[18]

      $$ {\left( {\mathit{\boldsymbol{\omega }} + \mathit{\boldsymbol{S}}_\beta ^{\rm{T}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{S}}_\beta }} \right)\mathit{\boldsymbol{\hat a}} - \mathit{\boldsymbol{\omega a}} + }\\ {\mathit{\boldsymbol{S}}_\beta ^{\rm{T}}\mathit{\boldsymbol{W}}\left( {\sum\limits_{i = 1}^m {{\mathit{\boldsymbol{h}}_i}} {{\mathit{\boldsymbol{\hat \beta }}}_i}} \right) + \mathit{\boldsymbol{S}}_\beta ^{\rm{T}}\mathit{\boldsymbol{W}}y = 0} $$ (6a)
      $$ {\left( {{\mathit{\boldsymbol{N}}_h} + {\mathit{\boldsymbol{N}}_B} + {\mathit{\boldsymbol{N}}_{Bh}} + {\mathit{\boldsymbol{N}}_{hB}}} \right)\mathit{\boldsymbol{\hat \beta }} - }\\ {{\mathit{\boldsymbol{u}}_h} - {\mathit{\boldsymbol{u}}_B} - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} = 0} $$ (6b)
      $$ {\mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{d}}_i} \ge 0, i = 1, 2 \cdots s} $$ (6c)
      $$ {{{\hat \lambda }_i}\left( {\mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{g}} - {d_i}} \right) = 0, i = 1, 2 \cdots s} $$ (6d)
      $$ {{{\hat \lambda }_i} \ge 0, i = 1, 2 \cdots s} $$ (6e)

      式中,$ {\mathit{\boldsymbol{S}}_\beta }、{\mathit{\boldsymbol{N}}_h}、{\mathit{\boldsymbol{N}}_B}、{\mathit{\boldsymbol{N}}_{Bh}}、{\mathit{\boldsymbol{N}}_{hB}}、{\mathit{\boldsymbol{u}}_h}、{\mathit{\boldsymbol{u}}_B}$的具体表达式见文献[11]。式(6a)和式(6b)分别由式(5a)和式(5b)推导得到,给出了最优解和Lagrange乘子的关系。式(6c)是约束的分量形式,令$\mathit{\boldsymbol{G}} = {\left( {{\mathit{\boldsymbol{g}}_1}\quad {\mathit{\boldsymbol{g}}_2} \cdots {\mathit{\boldsymbol{g}}_s}} \right)^{\rm{T}}}, \mathit{\boldsymbol{g}}_i^{\rm{T}} = \left( {{\mathit{\boldsymbol{g}}_{i1}}\quad {\mathit{\boldsymbol{g}}_{i2}} \cdots {\mathit{\boldsymbol{g}}_{im}}} \right) $表示$\mathit{\boldsymbol{G}} $的第$i $行,相应的常数项为$ \mathit{\boldsymbol{d}} = {\left( {{d_1}\quad d2 \cdots {d_s}} \right)^{\rm{T}}}$。式(6d)称为线性互补条件,当$ \mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{d}}_i} > 0$时,$ {\lambda _i} = 0$,是无效约束; 当$\mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{d}}_i} = 0 $时,$ {\lambda _i} \ge 0$,是有效约束[18]。式(6e)表示Lagrange乘子均为正。KT条件是式(2)取得最优解的必要条件。当获得式(2)的解以后,若需要判断其是否为最优解,还需要用到如下的二阶充分条件。

    • 首先定义目标函数对未知向量的二阶Hessian矩阵为:

      $$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{H}} = \frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{{\partial ^2}{{\left( {{{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}, {\mathit{\boldsymbol{\beta }}^{\rm{T}}}} \right)}^{\rm{T}}}}} = }\\ {\left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}\partial {{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}}}}&{\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}\partial {\mathit{\boldsymbol{\beta }}^{\rm{T}}}}}}\\ {\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}\partial {{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}}}}&{\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}\partial {\mathit{\boldsymbol{\beta }}^{\rm{T}}}}}} \end{array}} \right]} \end{array} $$ (7)

      Hessian矩阵式(7)中的子矩阵可由式(5a)和式(5b)求二阶偏导数获得,分别为:

      $$ \begin{array}{*{20}{c}} {\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}\partial {{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}}} = 2\mathit{\boldsymbol{\omega }} + }\\ {2{\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat \beta }} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{W}}\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{B}}} \end{array} $$ (8a)
      $$ \begin{array}{*{20}{l}} {\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}\partial {\mathit{\boldsymbol{\beta }}^{\rm{T}}}}} = 2{\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\beta }} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{W}} \times }\\ {\left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{h}}_1} + {\mathit{\boldsymbol{B}}_1}\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}&{{\mathit{\boldsymbol{h}}_2} + {\mathit{\boldsymbol{B}}_2}\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}& \cdots &{\left. {{\mathit{\boldsymbol{h}}_m} + {\mathit{\boldsymbol{B}}_m}\mathit{\boldsymbol{{\hat {\bar{ a}}}}}} \right] + }\\ {2\left[ {\mathit{\boldsymbol{B}}_1^{\rm{T}}\mathit{\boldsymbol{W{\hat {\bar{ A}}}\hat \beta }}} \right.}&{\mathit{\boldsymbol{B}}_2^{\rm{T}}\mathit{\boldsymbol{W{\hat {\bar{ A}}}\hat \beta }}}& \cdots &{\left. {\mathit{\boldsymbol{B}}_m^{\rm{T}}\mathit{\boldsymbol{W{\hat {\bar{ A}}}\hat \beta }}} \right] - }\\ {2\left[ {\mathit{\boldsymbol{B}}_1^{\rm{T}}Wy} \right.}&{\mathit{\boldsymbol{B}}_2^{\rm{T}}Wy}& \cdots &{\mathit{\boldsymbol{B}}_m^{\rm{T}}Wy} \end{array}} \right]} \end{array} $$ (8b)
      $$ \frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}\partial {{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}}} = {\left( {\frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}\partial {\mathit{\boldsymbol{\beta }}^{\rm{T}}}}}} \right)^{\rm{T}}} $$ (8c)
      $$ \begin{array}{l} \frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}\partial {\mathit{\boldsymbol{\beta }}^{\rm{T}}}}} = 2\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \mathit{\boldsymbol{h}}_1^{\rm{T}} + {{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_1^{\rm{T}}\\ \mathit{\boldsymbol{h}}_2^{\rm{T}} + {{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_2^{\rm{T}} \end{array}\\ \vdots \\ {\mathit{\boldsymbol{h}}_m^{\rm{T}} + {{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_m^{\rm{T}}} \end{array}} \right]\mathit{\boldsymbol{W}} \times \\ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{h}}_1} + {\mathit{\boldsymbol{B}}_1}\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}&{{\mathit{\boldsymbol{h}}_2} + {\mathit{\boldsymbol{B}}_2}\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}& \cdots &{{\mathit{\boldsymbol{h}}_m} + {\mathit{\boldsymbol{B}}_m}\mathit{\boldsymbol{{\hat {\bar{ a}}}}}} \end{array}} \right] \end{array} $$ (8d)

      式中,$\mathit{\boldsymbol{{\hat {\bar{ A}}}\hat \beta }} = \left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B{\hat {\bar{ a}}}}}), \mathit{\boldsymbol{{\hat {\bar{ A}}}}} = {\mathop{\rm Ivec}\nolimits} (\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B{\hat {\bar{ a}}}}}) $,${\rm Ivec} $是将$nm $维向量转化为$ n \times m$维矩阵的操作算子。将上述式(8a)至式(8d)代入式(7),就得到了目标函数的Hessian矩阵。由文献[18]可得ICPEIV模型的二阶充分条件为:若式(2)的目标函数和不等式约束二次连续可微,最优解满足一阶必要条件式(6),且有:

      $$ {\mathit{\boldsymbol{\mu }}^{\rm{T}}}\mathit{\boldsymbol{H\mu }} > 0 $$ (9)

      则可行点$ \mathit{\boldsymbol{\hat u}} = {\left( {{{\mathit{\boldsymbol{{\hat {\bar{ a}}}}}}^{\rm{T}}}, {{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}}} \right)^{\rm{T}}}$是ICPEIV模型的局部最优解。式中, $ \mathit{\boldsymbol{\mu }}$是可行有效方向,由下式定义:

      $$ \mathit{\boldsymbol{D}} = \left\{ {\mathit{\boldsymbol{\mu }}|\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\mu }} \ne 0}\\ {{\mathit{\boldsymbol{G}}_{\vartheta i}}\mathit{\boldsymbol{\mu }} = 0, i \in \mathit{\boldsymbol{I}}\;\;{\rm{ and}}\;\;{\rm{ }}{\mathit{\boldsymbol{\lambda }}_i} > 0}\\ {{\mathit{\boldsymbol{G}}_{\vartheta i}}\mathit{\boldsymbol{\mu }} \ge 0, i \in \mathit{\boldsymbol{I}}\;\;{\rm{ and}}\;\;{\rm{ }}{\mathit{\boldsymbol{\lambda }}_i} = 0} \end{array}} \right\} $$ (10)

      式中, $ \mathit{\boldsymbol{I}}$表示有效约束指标集; ${{\mathit{\boldsymbol{\lambda }}_i}} $表示对应的Lagrange乘子。令${\mathit{\boldsymbol{G}}^\prime } = {0_{s \times t}}, \quad {\mathit{\boldsymbol{d}}^\prime } = {0_{s \times 1}}, \quad \mathit{\boldsymbol{q}} = {\left( {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{d}}^{\prime {\rm{T}}}}}&{{\mathit{\boldsymbol{d}}^{\rm{T}}}} \end{array}} \right)^{\rm{T}}}$, 则式(1c)等价于:

      $$ \left(\mathit{\boldsymbol{G}}^{\prime} \quad \mathit{\boldsymbol{G}}\right) \mathit{\boldsymbol{u}}-\mathit{\boldsymbol{q}} \geqslant 0 $$ (11)

      则式(10)中$ {{\mathit{\boldsymbol{G}}_{\vartheta i}}}$表示矩阵${\mathit{\boldsymbol{G}}_\vartheta } = \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{G}}^\prime }}&\mathit{\boldsymbol{G}} \end{array}} \right) $中只保留有效约束对应的行所组成的矩阵。首先计算参数的最优值${\mathit{\boldsymbol{\hat u}}} $,将其代入式(7)计算Hessian矩阵; 然后确定有效约束,求式(10)的通解,将其代入式(9)判断解的最优性。

    • 对于约束非线性规划问题:

      $$ \left\{ {\begin{array}{*{20}{l}} {\min f(\mathit{\boldsymbol{u}})}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}{\mathit{\boldsymbol{G}}_\vartheta }\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{q}} \ge 0} \end{array}} \right. $$ (12)

      给定初始值$ {\mathit{\boldsymbol{u}}^{\left( j \right)}}$,计算函数$ {f(\mathit{\boldsymbol{u}})}$在$ {\mathit{\boldsymbol{u}}^{\left( j \right)}}$处的梯度$\mathit{\boldsymbol{g}}\left( {{\mathit{\boldsymbol{u}}^{\left( j \right)}}} \right) $和Hessian矩阵$\mathit{\boldsymbol{H}}\left( {{\mathit{\boldsymbol{u}}^{\left( j \right)}}} \right) $。将目标函数和约束条件在$ {\mathit{\boldsymbol{u}}^{\left( j \right)}}$处用泰勒级数展开保留二次项,则式(12)转化为QP问题[18]

      $$ \left\{ {\begin{array}{*{20}{c}} {\min f\left( {{\mathit{\boldsymbol{u}}^{(j)}} + {\rm{d}}\mathit{\boldsymbol{u}}} \right) = \frac{1}{2}{\rm{d}}{\mathit{\boldsymbol{u}}^{\rm{T}}}H\left( {{\mathit{\boldsymbol{u}}^{(j)}}} \right){\rm{d}}\mathit{\boldsymbol{u}} + }\\ {{g^{\rm{T}}}\left( {{\mathit{\boldsymbol{u}}^{(j)}}} \right){\rm{d}}\mathit{\boldsymbol{u}} + f\left( {{\mathit{\boldsymbol{u}}^{(j)}}} \right)}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\quad {\mathit{\boldsymbol{G}}_\vartheta }{\mathit{\boldsymbol{u}}^{(j)}} - \mathit{\boldsymbol{q}} + {\mathit{\boldsymbol{G}}_\vartheta }{\rm{d}}\mathit{\boldsymbol{u}} \ge 0} \end{array}} \right. $$ (13)

      采用AS法求解式(13),得到$ {{\rm{d}}\mathit{\boldsymbol{u}}}$的最优值,令新的近似值$\boldsymbol{u}^{(j+1)}=\boldsymbol{u}^{(j)}+\mathrm{d} \boldsymbol{u}(j=1, 2 \cdots N) $,形成新的QP子问题。其计算步骤为:

      1) 令$\boldsymbol{u}^{(0)}=\left(\overline{\boldsymbol{a}}^{(0) \mathrm{T}} \quad \boldsymbol{\beta}^{(0) \mathrm{T}}\right)^{\mathrm{T}} $,计算$ {\mathit{\boldsymbol{g}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{u}}^{(0)}}} \right)$和${\mathit{\boldsymbol{H}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{u}}^{(0)}}} \right) $,在可行集内选取初值${\rm{d}}{\mathit{\boldsymbol{u}}_0} $,利用AS法求解QP子问题;

      2) 若$\parallel {\rm{d}}\mathit{\boldsymbol{u}}\parallel \le \varepsilon $($\varepsilon $是任意给定的小的正数),则停止迭代; 若$ \parallel {\rm{d}}\mathit{\boldsymbol{u}}\parallel > \varepsilon $,令$ \boldsymbol{u}^{(j)}=\boldsymbol{u}^{(j-1)}+\mathrm{d} \boldsymbol{u}$,转步骤1)。

      对于ICPEIV模型,目标函数$f\left( \mathit{\boldsymbol{u}} \right) = f{\left( {{{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}\quad {\mathit{\boldsymbol{\beta }}^{\rm{T}}}} \right)^{\rm{T}}} $由式(2)确定,而梯度函数$ g\left( \mathit{\boldsymbol{u}} \right) = {\left( {\frac{{\partial {\mathit{\Phi} ^{\rm{T}}}(\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}}}\frac{{\partial {\mathit{\Phi} ^{\rm{T}}}(\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}}}} \right)^{\rm{T}}}$和Hessian矩阵$ \boldsymbol{H}\left( \mathit{\boldsymbol{u}} \right)$由式(14)确定,由$ \mathit{\Phi} (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})$和$\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }}) $的关系易得:

      $$ \frac{{\partial \mathit{\Phi} (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}}} = \frac{{\partial \psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}}} $$ (14a)
      $$ \frac{{\partial \mathit{\Phi} (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}}} = \frac{{\partial \psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}}} + {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} $$ (14b)
      $$ \frac{{{\partial ^2}\mathit{\Phi} (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{{\partial ^2}{{\left( {{{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}, {\mathit{\boldsymbol{\beta }}^{\rm{T}}}} \right)}^{\rm{T}}}}} = \frac{{{\partial ^2}\psi (\mathit{\boldsymbol{\bar a}}, \mathit{\boldsymbol{\beta }})}}{{{\partial ^2}{{\left( {{{\mathit{\boldsymbol{\bar a}}}^{\rm{T}}}, {\mathit{\boldsymbol{\beta }}^{\rm{T}}}} \right)}^{\rm{T}}}}} = \mathit{\boldsymbol{H}} $$ (14c)

      采用SQP算法求解问题式(2)时,$ {\mathit{\boldsymbol{\bar a}}}$的初始值等于系数矩阵元素,而模型参数$\mathit{\boldsymbol{\beta }} $的初始值由相应的ICLS解确定。根据文献[18]的结论,假设$ {\mathit{\boldsymbol{u}}^*} = {\left( {{\mathit{\boldsymbol{\beta }}^{*{\rm{T}}}}\quad {{\mathit{\boldsymbol{\bar a}}}^{*{\rm{T}}}}} \right)^{\rm{T}}}$是模型式(1)的局部最优解,且满足一阶KT条件式(6)。当$ \lambda _i^* = 0$时,$ \mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{\beta }} - {d_i} > 0$这一严格补充条件成立,且满足二阶充分条件式(9)。若$\left( {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}^{(k)}}}&{{\mathit{\boldsymbol{\lambda }}^{(k)}}} \end{array}} \right) $充分接近最优解$\left( {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}^*}}&{{\mathit{\boldsymbol{\lambda }}^*}} \end{array}} \right) $,则存在子问题式(13)的局部最优解,使它和非线性规划问题式(2)有相同的有效集。SQP算法是一种牛顿型算法,具有二阶收敛速率。算法的不确定度主要由AS法求解QP子问题的不确定性来决定。若QP问题式(13)中的Hessian矩阵对称正定,且每步迭代中的有效约束线性无关,那么AS法在有限步内能得到问题式(13)的唯一全局最优解。由于AS法的每一迭代步是解等式约束LS问题,当Hessian矩阵病态时会引起解的不稳定性。

    • 算例采用文献[14]提出的线性近似方法和§2提出的SQP法计算,并进行比较。前者是将式(1)在待求参数$ {\mathit{\boldsymbol{\bar a}}}$和$\mathit{\boldsymbol{\beta }} $的近似值处线性化,转换为ICLS平差模型后迭代计算。与文献[14]不同的是,本文在解近似ICLS模型时采用AS法取代文献[14]中的LCP方法。这是为了与SQP方法中解QP子问题采用AS法保持一致,从而便于分析两种算法的计算效率。

    • 对于ICEIV模型:

      $$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{y}} = \left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_A}} \right)\mathit{\boldsymbol{\beta }} + {\mathit{\boldsymbol{e}}_y}}\\ {\mathit{\boldsymbol{G\beta }} \ge \mathit{\boldsymbol{d}}} \end{array}} \right. $$ (15)

      已知真值$ \mathit{\boldsymbol{\tilde A}} = {\left[ {\begin{array}{*{20}{c}} {0.25}&{0.25}&{0.5}&1\\ 1&1&1&1 \end{array}} \right]^{\rm{T}}}, \mathit{\boldsymbol{\tilde \beta }} = {\left[ {\begin{array}{*{20}{l}} {0.6}&{0.4} \end{array}} \right]^{\rm{T}}} $,计算观测值的真值为$ \mathit{\boldsymbol{\tilde y}} = {\left[ {\begin{array}{*{20}{l}} {0.55}&{0.55}&{0.7}&1 \end{array}} \right]^{\rm{T}}}$。约束矩阵$ \mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} 1&0\\ { - 1}&{ - 1} \end{array}} \right]$,约束项$ \mathit{\boldsymbol{d}} = {\left[ {\begin{array}{*{20}{l}} 0&{ - 1} \end{array}} \right]^{\rm{T}}}$。分别对系数阵和观测项加入均值为0、方差为0.2的随机误差。首先采用AS法求解相应的ICLS解作为初值,然后将式(15)转换成模型式(1),分别采用线性化方法和SQP方法求解,迭代终止标准均取$ \varepsilon=1.0 \times 10^{-8}$。参数估值、迭代次数和真值与估值的偏差列于表 1(其中$ {{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{ICLS }}}}$表示ICLS解,$\mathit{\boldsymbol{\hat \beta }}_{{\rm{ICWTLS }}}^1 $和$ \mathit{\boldsymbol{\hat \beta }}_{{\rm{ICWTLS }}}^2$分别表示采用线性化方法和SQP方法的ICWTLS解)。

      表 1  不等式约束PEIV模型的真值模拟实验结果

      Table 1.  Simulation Results of ICPEIV Model with True Values

      统计项 ${{\hat \beta }_1} $ ${{\hat \beta }_2} $ 外迭代次数 平均内迭代次数 $\parallel \hat \beta - \tilde \beta {\parallel ^2} $
      $ {\tilde \beta }$ 0.600 000 0.400 000 0
      $ {{\hat \beta }_{{\rm{ICLS }}}}$ 0.577 568 0.422 432 ${4.19 \times {{10}^{ - 3}}} $
      $\hat \beta _{{\rm{ICWTLS }}}^1 $ 0.590 989 0.409 011 13 2 $\begin{array}{l} 1.62 \times 10^{-4} \end{array} $
      $\hat \beta _{{\rm{ICWTLS }}}^2 $ 0.590 990 0.409 010 4 2 $ {1.62 \times {{10}^{ - 4}}}$

      从真值和估计的差值范数来看,ICWTLS解考虑了系数矩阵和右端观测项的误差,且不等式约束信息正确,从而得到了更接近真值的解。ICLS解未顾及系数阵元素误差,估计精度低于ICWTLS解。采用两种方法得到的ICWTLS解估值接近一致,说明SQP方法可行。两种方法的内迭代都是相同维数的QP子问题,外迭代都是代入新的参数估值作为初值,内外迭代的复杂度相近。比较两种方法的内外迭代次数,可见SQP方法具有更好的计算效率。

    • 对于ICEIV模型式(15),观测方程的系数矩阵、右端向量、不等式约束矩阵及约束项见文献[17]中的算例5.1。采用与算例1相同的计算方案,结果见表 2

      表 2  不等式约束PEIV模型的模拟实验结果(文献[17]中算例)

      Table 2.  Simulation Results of ICPEIV Model(Example is from Reference[17])

      统计项 $ {{\hat \beta }_1}$ $ {{\hat \beta }_2}$ ${{\hat \beta }_3} $ $ {{\hat \beta }_4}$ 外迭代次数 平均内迭代次数
      $ {{\hat \beta }_{{\rm{ICLS }}}}$ -0.100 000 -0.100 000 0.215 228 0.350 152
      ${{\hat \beta }_{{\rm{WTLS }}}} $ 0.188 761 -0.716 733 0.560 517 0.210 638
      $ \hat \beta _{{\rm{ICWTLS }}}^1$ -0.100 000 -0.100 000 0.168 581 0.399 765 23 2
      $ \hat \beta _{{\rm{ICWTLS }}}^2$ -0.100 000 -0.100 000 0.168 579 0.399 766 5 2

      表 2中${{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{WTLS }}}} $表示WTLS解,可以看到它不满足约束关系。两种方法得到的ICWTLS解结果一致,SQP方法的外迭代次数明显小于线性化方法。算例中第2、4、6个约束为有效约束,相应的Lagrange乘子为$ \mathit{\boldsymbol{\lambda }} = {\left( {\begin{array}{*{20}{l}} {0.4817}&{0.1424}&{0.5049} \end{array}} \right)^{\rm{T}}}$,满足一阶必要条件。根据式(7)计算Hessian矩阵$ \boldsymbol{H}\left( \mathit{\boldsymbol{u}} \right)$,然后进行特征值分解,存在一个负特征值为$ \sigma_{\min }=-0.0217$。因此,$ \boldsymbol{H}\left( \mathit{\boldsymbol{u}} \right)$不是正定矩阵。根据式(10),组成如下线性方程组:

      $$ \left[\begin{array}{ccc} 0 & 0 & 0 \\ \vdots & \vdots & \vdots \\ 0 & 0 & 0 \\ 0.1987 & -1 & 0 \\ 0.1988 & 0 & -1 \\ 0.4450 & 0 & 0 \\ 0.4186 & 0 & 0 \end{array}\right]^{\mathrm{T}}\mathit{\boldsymbol{\mu }} = 0 $$ (16)

      相应的通解为:

      $$ \mathit{\boldsymbol{\mu }} = {\left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\mu }}_a^{\rm{T}}}&{\mathit{\boldsymbol{\mu }}_\beta ^{\rm{T}}} \end{array}} \right]^{\rm{T}}} $$ (17)

      其中, $ {\mathit{\boldsymbol{\mu }}_a} = {\left[ {\begin{array}{*{20}{c}} {{\mu _1}}&{{\mu _2}}& \cdots &{{\mu _{20}}} \end{array}} \right]^{\rm{T}}}, {\mathit{\boldsymbol{\mu }}_\beta } = {\left[ {\begin{array}{*{20}{c}} 0&0&{\begin{array}{*{20}{c}} {{\mu _{21}}}&{ - 1.0691}&{{\mu _{21}}} \end{array}} \end{array}} \right]^{\rm{T}}}$。计算${\mathit{\boldsymbol{\mu }}^{\rm{T}}}\mathit{\boldsymbol{H\mu }} $,采用无约束最优化算法求解$ s(\mathit{\boldsymbol{\mu }}) = {\mathit{\boldsymbol{\mu }}^{\rm{T}}}\mathit{\boldsymbol{H\mu }}$的极值,当$\mathit{\boldsymbol{\mu }} = 0 $时, 取极小值$ {s_{\min }}(\mu ) = 0$。因此,当$ \mathit{\boldsymbol{\mu }} \ne 0$时,$ s(\mathit{\boldsymbol{\mu }}) > 0$恒成立。因此,ICWTLS解满足二阶充分条件,是严格局部最优解。

    • 直线拟合在大地测量中很常见,它是根据一系列含有误差的离散点,在平面上拟合出最佳的直线。选取文献[14]中算例5.2的直线拟合模型,两种算法计算的结果如表 3所示。

      表 3  直线拟合的不等式约束整体最小二乘解

      Table 3.  ICWTLS Solution of Line Fitting

      统计项 ${{\hat \beta }_1} $ $ {{\hat \beta }_2}$ 外迭代次数 平均内迭代次数
      $ {{\hat \beta }_{{\rm{ICLS }}}}$ 1.950 000 0.450 000
      $ {{\hat \beta }_{{\rm{WTLS }}}}$ 2.202 232 0.498 898
      $ \hat \beta _{{\rm{ICWTLS }}}^1$ 2.025 044 0.500 088 5 3
      $ \hat \beta _{{\rm{ICWTLS }}}^2$ 2.025 044 0.500 088 4 2

      表 3中可以看到,线性化方法的外迭代次数和内迭代次数都大于SQP方法。算例中最后一个约束为有效约束,Lagrange乘子$\lambda=1.3669 $。对$\mathit{\boldsymbol{H}}\left( \mathit{\boldsymbol{u}} \right) $进行特征值分解,最小特征值$ \sigma_{\min }=1.9418$,是正定矩阵,对任意向量$ \mathit{\boldsymbol{\mu }}$都有$ s(\mathit{\boldsymbol{\mu }}) = {\mathit{\boldsymbol{\mu }}^{\rm{T}}}\mathit{\boldsymbol{H\mu }}{\rm{ > 0}}$恒成立。因此,ICWTLS解是严格局部最优解。

      从上述3个算例看出,采用SQP算法求解ICWTLS解可行有效,具有较好的计算效率。它能同时得到模型参数的最优值和设计矩阵元素的改正值。SQP算法在求解QP子问题时,能自动满足一阶必要KKT条件。不同于无约束规划问题取得最优解的充分条件是要求目标函数的二阶Hessian矩阵正定,约束最优化问题的二阶充分条件与约束函数本身有关。在§3.2中,当目标函数的Hessian矩阵不正定时,ICWTLS解仍是局部最优解。在§3.3中,当Hessian矩阵正定时,不需要解线性方程组(17),就可以直接判定ICWTLS解局部最优。

    • PEIV模型能反映系数矩阵的结构特性,处理系数阵中随机元素和非随机元素共存的情况,只需求解模型参数和随机系数阵元素的估值,较之EIV模型有更广泛的适应性和更好的计算效率。合理的不等式约束能提高参数估计的可靠性和精度。本文以ICPEIV模型为基础,研究了ICWTLS解的最优性条件和算法。首先构建了ICWTLS问题的Lagrange目标函数,根据约束非线性规划理论推导了ICWTLS解取得最优值的一阶必要条件和二阶充分条件。然后,提出了一种求解ICWTLS解的SQP算法。最后,采用数值算例和平面直线拟合算例进行计算,证明该算法是一种有效可行的算法,具有良好的计算效率。SQP算法在迭代计算中自动满足一阶必要条件。当最优解处的Hessian矩阵正定时,二阶充分条件恒成立; 当Hessian矩阵不是正定矩阵时,若二次型式(9)成立,则ICWTLS解仍是局部最优解。

参考文献 (18)

目录

    /

    返回文章
    返回