参数约束

2024-07-08

参数约束（通用7篇）

参数约束篇1

1 概述

传统CAD系统无法很好的支持功能设计、概念设计等操作, 为了弥补传统CAD系统的不足, 更好地满足概念化设计的要求、提高设计效率, 在20世纪70年代, CAD系统首次引入了人工智能的思想, 其主要特征和标志为参数化和变量化技术。

参数化技术与变量化技术的核心是:当设计者指定了草图的尺寸和拓扑关系后, 系统能够自动生成相应的结果图形。这一求解过程称为几何约束求解。几何约束求解对CAD系统的性能起着至关重要的作用, 并可应用于其他众多的工程领域。

几何约束求解的数学思想是:将几何约束转化为代数方程组, 然后对方程组进行求解, 得到问题的解。几何约束涉及大量几何元素, 在转化为代数方程组时会产生大型非线性方程组。目前还没有求解大型非线性方程组的稳定方法, 为了提高求解效率、减小求解复杂度, 一个基本思想是将大型系统分解为若干个较小规模的子系统, 分别求解, 然后将各个求解结果组装为整个系统的最终解。几何约束的求解方法融合了很多领域的知识, 如几何、图论、矩阵、组合优化、人工智能和离散数学等, 根据不同的求解方式, 可将几何约束求解方法分为:基于代数的求解方法、基于几何约束图的求解方法和基于人工智能的规则推理法。

2 约束求解方法

2.1 基于代数的求解方法

基于代数的求解方法是将几何约束问题转化为非线性方程组来进行求解。这种方法利用变量和方程表达几何约束系统, 具有很强的通用性和直观性。代数求解法中常用的方法是数值法和符号法。

数值法是直接将几何约束转化为非线性方程组, 然后采用整体数值迭代求解。Newton-Raphson迭代法是使用最多的一种数值计算方法, 通过迭代过程Xk+1=Xk-J (Xk) -1F (Xk) 逐渐逼近方程组的解, 这里J (Xk) 是F (Xk) 在点Xk的Jacobi矩阵, F= (f1, f2, …, fm) T:是方程组向量, X= (x1, x2, …, xm) T是变元向量。这种方法的求解速度快, 使用范围广, 但只能求得一个解。为了求得问题的全部解, 可以将几何约束满足问题转化为优化问题进行求解。Joan利用遗传算法来求解几何约束问题, 但是这种方法只能处理规模较小的问题。欧阳应秀提出将混沌方法嵌入BFGS算法的混合求解方法, 能求解欠约束和过约束情况。

用符号法求解几何约束时, 首先将几何约束映射为代数方程组, 然后利用消元法化简代数方程组, 并采用数值回代的方法进行求解。符号法主要可以分为Grobner基法、吴方法以及结式法。高小山将吴方法应用于几何约束求解。这种方法可以判断完备约束、欠约束和过约束的情况, 但是时间复杂度和空间复杂度都较大。

2.2 基于几何约束图的求解方法

基于几何约束图的方法是对几何约束图进行分析, 将问题分解为可分别求解的子问题, 然后将各个子问题的解合并为最终问题的解。根据分析策略的不同, 可将基于几何约束图的方法分为构造法和自由度分析法。

在构造法中, 根据几何约束图的不同分解形式, 可以将构造法分为递归分割法和递归装配法。递归分割法是通过一定的规则, 将几何约束系统分割为若干子系统, 再继续递归分割各个子系统。Owen首次提出递归分割法, 该方法以图论中的关节点、双连通和三连通等概念为基础, 与深度优先搜索和广度优先搜索算法相结合, 自顶向下分割几何约束图。Joan提出了Deficit的概念, 在保持约束图的Deficit不变的条件下, 分解二叉树s-tree。张桂芳将三维完备几何约束图分解为k-连通子图, 如果其中一个子图是完备约束的, 则求解该子图后, 可以将这些子图的解合并为原问题的解。递归装配法以迭代的方式将刚性体组合成更大的刚性体。如果约束系统中包含完备约束和欠约束子系统, 则递归装配法返回一个完备约束组件的最大集合。Podgorelec通过变换点和直线之间的角度和距离约束, 将复杂约束转化为简单约束;将复杂几何元素映射为辅助线和辅助点, 来处理复杂几何元素。

Kramer以刚体运动学为基础提出了几何约束系统的自由度分析法, 给出了链路搜索和环路搜索算法, 并提出:“几何约束满足问题的核心是实现约束的最大分解”。

彭小波分析了几何约束系统的奇异性, 提出了一种基于有向约束图的欠约束系统的求解的方法, 通过分解几何约束来减小约束求解的规模。蒋丹东引入了点簇的概念, 综合了基于图和基于规则的约束求解方法, 运用图论的原理建立了基于点簇的图形约束模型。李彦涛引入形状自由度的概念, 将剪枝和凝聚相结合, 实现欠约束和完备约束系统的分解, 并使用解析法和数值法对约束进行求解。在欠约束系统分解方面, Lee将非尺规可构造构型分离出来并利用数值计算法进行求解, 利用分类规则找出欠约束子图, 并得到求解序列。Joan通过添加几何约束, 将欠约束转化为可求解的完备几何约束问题。高小山以D-Tree的分解方法为基础, 通过自动添加几何约束, 将欠约束问题转化为可求解的完备几何约束问题, 使分解后的几何约束问题具有最小的规模。

2.3 基于人工智能的规则推理法

基于人工智能的规则推理法利用谓词和一系列重写规则来构造几何元素, 从而将几何约束问题转化为可构造的形式。

Aldefeld采用基于符号推理和操作的专家系统来建立规则体系, 用一阶谓词描述几何约束关系, 然后利用推理机对知识库进行规则匹配, 从而构造出整个图形。

日本东京大学Kimura的研究小组提出了基于规则的构造方法, 通过两种方法得到约束:一种是由系统模型自动生成, 另一种是由用户交互定义。这种方法考虑了系统的多解情况, 但没有考虑约束的一致性检验问题。Verroust提出了面向规则的平行四边形变换方法, 该方法能够解决多边形的所有完备约束问题, 但应用范围有限。高曙明引入已知元素和已知约束的概念, 利用普通算法实现几何推理。Dufourd采用延拓法对几何约束系统进行自动符号推理, 并提出以树的形式表达解空间, 通过剪枝法处理多解性问题。Schreck利用距离比率约束来替换距离约束, 将几何约束问题转变为相似群中的不变形式, 并进行构造求解。Podgorelec给出了分析几何约束是否可解的规则, 并提出了简化几何约束和冗余校验的方法。Fabre将符号法和数值法相结合来简化几何约束, 但该方法仍然以牛顿迭代法为基础, 求解的稳定性和效率有待进一步提高。

3 结论

本文分析介绍了国内外现有的几何约束求解方法。几何约束求解是CAD造型系统的关键技术, 对CAD系统的性能起着至关重要的作用。在造型系统中, 合理有效的几何约束求解机制将会有效地提高造型效率。因此, 研究几何约束求解方法具有重要的现实意义。

摘要：本文介绍了参数化造型中的几何约束求解方法的发展状况, 从基于代数的求解方法、基于几何约束图的求解方法和基于人工智能的规则推理法等方面分析了目前造型系统中采用的几何约束求解方法。

关键词：约束求解,几何约束图,人工智能

参考文献

[1]高雪瑶.语义特征造型的与历程无关技术的研究[D].哈尔滨:哈尔滨理工大学, 2009.[1]高雪瑶.语义特征造型的与历程无关技术的研究[D].哈尔滨:哈尔滨理工大学, 2009.

参数约束篇2

基于非线性约束条件研究了一类不能求解析解的微分方程系统模型的`参数估计方法.通过非线性系统模型和观测模型的线性化,用线性模型理论推导了非线性约束条件下的参数估计方法,设计了全局最优解的高斯牛顿迭代算法.某卫星轨道确定算例表明:该参数估计法的卫星轨道确定精度较传统最小二乘法提高约50%.?

作者：刘靖潘晓刚 LIU Jing PAN Xiao-gang 作者单位：刘靖,LIU Jing(湖南涉外经济学院,湖南,长沙,410205)

潘晓刚,PAN Xiao-gang(国防科学技术大学,理学院,湖南,长沙,410073)

参数约束篇3

全局优化问题在科学计算、工程技术、经济管理等领域得到越来越广泛的关注和应用, 近些年来, 人们相继提出一些求解全局优化问题的算法, 填充函数法属于其中的确定型算法, 最早由Ge[1]提出, 用来求解无约束多极值函数的全局极小点。这种算法的关键是构造一个称为填充函数的辅助函数, 其基本思想是借助辅助函数, 从目标函数的当前局部极小点找到另一个目标函数值更小的局部极小点, 直至找到全局最小点为止。

Ge之后, 不少学者相继提出了不少性质良好的填充函数, 但基本上都是针对无约束优化问题的, 对于求解约束优化问题的填充函数虽有讨论, 如文献[2,3,4,5,6,7], 但都有各自不同程度的缺陷。本文在无强制性条件下给出了一类性质良好的求解带一般约束优化问题的单参数填充函数, 并讨论了其填充性质。

1 基本概念

考虑带约束全局优化问题 (P) : min f (x) , s.t.x∈S={x∈Rn|gi (x) ≤0, i∈Γ}, 其中f (x) , gi (x) , i∈Γ:Rn→R是连续可微函数, Γ={1, 2, …, m}是下标集。对上述问题假设如下。

假设1 问题 (P) 的不同局部极小点的个数可以是无限的, 但不同局部极小值个数是有限的。

假设2 用L (p) 表示问题 (P) 的局部极小点集合, 若x*是问题 (P) 的局部极小点, 则Lx*={x∈L (p) |f (x) =f (x*) }是一个有界闭集, 且问题 (P) 存在全局极小点。

注:在很多全局优化的实际问题中, f (x) 的强制性条件不一定满足, 因此这里对问题 (P) 的假设中没有考虑到目标函数f (x) 的强制性条件。下面针对问题 (P) 给出有约束优化问题的填充函数定义。

定义1[4] 函数p (x, x*, a) 称为f (x) 在局部极小点x*处的填充函数, 如果它满足:

(1) 在Rn空间中, x*是p (x, x*, a) 的严格局部极大点;

(2) 对x∈S1∩S且x≠x*或x∈Rn＼S有∇p (x, x*, a) ≠0, 这里S1={x∈Rn|f (x) ≥f (x*) };

(3) 如果x*不是全局极小点, 那么p (x, x*, a) 一定在S2={x∈S|f (x) <f (x*) }上有局部极小点。

以上定义的填充函数的意义为:对于有约束最优化问题, 首先要考虑的就是求出来的点是否可行。由条件 (2) 知, 要找的点一定在可行域内;其次, 对于全局最优化问题的填充函数法, 关心的是比当前局部极小点更好的那些点, 由条件 (2) 知, 要找的点一定不会再比当前局部极小点差的水平集上达到。若x*是局部极小点但不是全局极小点, 由定义中的条件 (1) 和条件 (3) , 则可以从x*的邻域中的任意一点出发, 用求无约束最优化问题的极小化方法极小化填充函数, 总能找到填充函数的局部极小点x*0, 由条件 (3) 知, f (x*0) <f (x*) , 再由条件 (2) , 对x∈S1∩S且x≠x*或x∈Rn＼S有∇p (x, x*, a) ≠0, 得知x*0∈S。

2 填充函数及其性质

针对问题 (P) , 设x*是当前局部极小点, 构造单参数填充函数如下。

1) p (x, x*, a) =-φ (1+‖x-x*‖3) +amin [0, max (f (x) -f (x*) , gi (x) , i∈Γ) ]。

其中a>0是参数, 函数φ (t) 满足:1) φ (0) =0, 2) φ′ (t) >0, 这样的函数有 $1 - e^{- t}, \arctan t, \frac{t}{1 + t}$ 等。

当参数a>0充分大时, 下面的几个定理表明p (x, x*, a) 是满足定义1的一类填充函数。

定理1 对任意a>0, x*是p (x, x*, a) 的严格局部极大点。

证明因为x*是f (x) 的局部极小点, 则存在它的一个邻域N (x*, δ) (δ>0) , 使得∀x∈N (x*, δ) ∩S, 有f (x) ≥f (x*) , gi (x) ≤0, i∈Γ。

下面分两种情况来证明对∀x∈N (x*, δ) , p (x, x*, a) <p (x*, x*, a) 成立。

(1) 当x∈N (x*, δ) ∩S, x≠x*时, 由于f (x) ≥f (x*) , 则有min [0, max (f (x) -f (x*) , gi (x) , i∈Γ) ]=0于是p (x, x*, a) =-φ (1+‖x-x*‖3) <-φ (1) =p (x*, x*, a) 。

(2) 当x∈N (x*, δ) ∩ (Rn＼S) 时, 则至少存在一个指标i0∈Γ, 使得gi0 (x) >0, 故min [0, max (f (x) -f (x*) , gi (x) , i∈Γ) ]=0, 于是同 (1) 有结论成立。

综上, p (x, x*, a) <p (x*, x*, a) 对于∀x∈N (x*, δ) 都成立。因此, x*是p (x, x*, a) 的严格局部极大点。

定理2 若x*是f (x) 的局部极小点, 则p (x, x*, a) 在 (S∩S1) ＼{x*}或Rn＼S上有∇p (x, x*, a) ≠0成立。

证明易得对∀x∈S∩S1或Rn＼S都有min [0, max (f (x) -f (x*) , gi (x) , i∈Γ) ]=0, 此时, p (x, x*, a) =-φ (1+‖x-x*‖3) , 显然∇p (x, x*, a) ≠0。

定理3 若x*是f (x) 的局部极小点, 但不是f (x) 的全局极小点, 且cl (int S) =cl (S) , 则当a>0充分大时, 一定存在x*0∈S2, 使得x*0是p (x, x*, a) 的局部极小点。

证明因为x*是f (x) 的局部极小点但非全局极小点, 则存在f (x) 的另一个局部极小点x*1, 使得f (x*1) <f (x*) , gi (x*1) ≤0, i∈Γ。

由于f (x) , gi (x) , i∈Γ是连续函数, 且cl (int S) =cl (S) , 则一定存在x*2∈Rn, 使得f (x*2) <f (x*) , gi (x*2) <0, 故有

p (x*2, x*, a) =-φ (1+‖x*2-x*‖3) +amax {f (x*2) -f (x*) , gi (x*2) , i∈Γ}。

显然Lx*1⊂S, x*2∈int S∩S2。

另一方面, 对∀x∈∂S, 至少存在一个指标i1∈Γ, 使得gi1 (x) =0, 于是当x∈∂S时有min [0, max (f (x) -f (x*) , gi (x) , i∈Γ]=0, 从而p (x, x*, a) =-φ (1+‖x-x*‖3) 。

所以, 当a>0充分大时, 对∀x∈∂S, 有p (x*2, x*, a) <p (x, x*, a) 。显然, S＼∂S是开集, 于是当a>0充分大时存在一点x*0∈S＼∂S使得

$\begin{array}{l} \min_{x \in S} p (x ‚ x^{*} ‚ a) = \min_{x \in S \partial S} p (x ‚ x^{*} ‚ a) = \\ p (x_{0}^{*} ‚ x^{*} ‚ a) \leq p (x_{2}^{*} ‚ x^{*} ‚ a) ‚ \end{array}$

故x*0∈int S, 且f (x*0) <f (x*) 。

定理4 若x*是f (x) 的全局极小点, 则对∀x∈S, x≠x*都有p (x, x*, a) <0。

证明由于x*是f (x) 的全局极小点, 则对所有的x∈S都有f (x) ≥f (x*) 成立。因此, 由定理1, 对∀x∈S, x≠x*有p (x, x*, a) <0成立。

定理5 任给x1, x2∈Rn, 若f (x1) ≥f (x*) , f (x2) ≥f (x*) , 则‖x2-x*‖>‖x1-x*‖当且仅当p (x2, x*, a) <p (x1, x*, a) 。

证明若f (x1) ≥f (x*) , f (x2) ≥f (x*) , 则min [0, max (f (xj) -f (x*) , gi (x) , i=∈Γ) ]=0, j=1, 2, 此时有p (xj, x*, a) =-φ (1+‖xj-x*‖3) , 显然有p (x2, x*, a) <p (x1, x*, a) , 反之亦然。

定理6 如果x1, x2∈Rn并且满足f (x1) ≥f (x*) ≥f (x2) 和‖x2-x*‖>‖x1-x*‖, 则p (x2, x*, a) <p (x1, x*, a) 。

证明由条件, p (x2, x*, a) =-φ (1+‖x2-x*‖3) +amax [f (x2) -f (x*) , gi (x2) , i∈Γ]。

且p (x1, x*, a) =-φ (1+‖x1-x*‖3) , 显然结论是成立的。

注:定理4-定理6说明目标函数当前的极小点x*被淘汰, 而新的填充函数的极小点将会被找到, 并且目标函数在这个点的函数值不比当前极小点处的函数值大。这也正是本文的填充函数所具有的良好性质之一。

3 算法和数值实验

3.1 算法

根据上述理论, 参考文献[7]给出下述算法。

(0) 初始步:选取初始点xk∈S;选取A>0作为可接受的终止参数;令a=1, k=1;

(1) 以xk为初始点, 应用局部下降算法求得问题 (P) 的一个局部极小点, 记作x*k;

(2) 选取初始点列{x $_{k + 1}^{i}$ :i∈Γ}, 使得对于某个δk>0有x $_{k + 1}^{i}$ ∈S＼N (x*k, δk) ;

(3) 令i=1;

(4) 若i≤m, 令x=x $_{k + 1}^{i}$ , 转步 (5) ;否则, 转步 (7) ;

(5) 若f (x) <f (x*k) , 且x∈S, 则以x为初始点, 应用已有局部下降算法求问题 (P) 的局部极小点x*k+1, 使得f (x*k+1) <f (x*k) , 令x*k=x*k+1, k=k+1, 转步 (2) ;否则, 转步 (6) ;

(6) 沿方向D1=-∇p (x, x*k, a) 或者 $D_{2} = - \frac{\nabla f (x)}{∥ \nabla f (x) ∥} - \frac{\nabla p (x, x_{k}^{*}, a)}{∥ \nabla p (x, x_{k}^{*}, a) ∥}$ , 找到一个新的点x, 使得p (x, x*, a) 能下降到一定程度。若在极小化过程中, x超出S的边界, 则令i=i+1, 转步 (4) , 否则, 转步 (5) ;

(7) 令a=2a。若a≤A, 转步 (3) ;否则, 算法终止, 视x*k为问题 (P) 的全局极小点。

注:下降方向D1、D2的合理性及优越性在文献[7]中有详细的说明, 不再赘述。

3.2 实验

通过两个算例的数值实验验证算法的可行性与有效性。算例都是在Matlab7.1.0运行环境中实现的。

算例1

$\min f (x) = 4 x_{1}^{2} - 4 x_{2}^{2} + 4 x_{2}^{4} - 2.1 x_{1}^{4} + \frac{1}{3} x_{1}^{6} - x_{1} x_{2}$

$s . t . x \in {(x_{1} ‚ x_{2}) | - 3 \leq x_{1} ‚ x_{2} \leq 3}$ 全局极小点x*= (0.090 1, 0.712 2) T, 全局极小值f (x*) =-1.133 1。

算例2

min f (x) = (x $_{1}^{6}$ -16x $_{1}^{5}$ +86x $_{1}^{4}$ -176x $_{1}^{3}$ +105x $_{1}^{2}$ ) /100+x $_{2}^{2}$ (x $_{2}^{2}$ -6x2+8) (x $_{2}^{2}$ -14x2+48) /100。

$s . t . x \in {(x_{1} ‚ x_{2}) | x_{1} + x_{2} \leq 10 ‚ x_{1} ‚ x_{2} \geq 0 ‚ 0 \leq x_{1} ‚ x_{2} \leq 8}$ 。

全局极小点x*= (2.174 0, 7.341 9) T, 全局极小值f (x*) =-9.123 1。

用k表示迭代次数, xk1表示第k次迭代初始点, x*k表示第k次迭代的局部最优解。迭代过程见表1。

参考文献

[1] Ge R P, Qin Y F. A class of filled functions for finding a global minimizer of a function of several variables.Journal of Optimization Theory and Applications, 1987;54 (2) :241—252

[2]Yang Yongjian, Shang Youlin.A filled function method for uncon-strained global optimization.Shanghai:Science College, Shanghai Univevs:ty, 2006;503—504

[3] Zhang L S, Kong Chi, Duan N G, et al.A new filled function method for global optimization.Global Optim, 2004;28:17—43

[4] Wang Weixiang, Shang Youlin, Zhang Liansheng. A filled function method with one parameter for box constrained global optimization.Applied Mathematics and Computation, 2007;194:54—66

[5] Wang Changyu, Li Duan. Unified theory of augmented Lagrangian methods for constrain-ed global optimization.J Glob Optim, 2009;44:433—458

[6] Wang Weixiang, Shang Youlin, Zhang Liansheng. A filled function method with one parameter for constrained global optim-ization.Journal of Engineering Mathemati-cs, 2008;25 (5) :795—803

参数约束篇4

不同的函数模型,对应着不同的平差方法,反之亦然,平差方法的不同,其对应的函数模型也不尽相同[1~3],本文首先简单介绍了条件平差、间接平差和附有约束条件的间接平差这3种经典平差方法的函数模型,推导了这3种平差模型的参数估计公式。尤其对于间接平差模型而言,由于其误差方程列立性强、精度评定便利,便于计算机编程等优势,因而被大多数人所采用,但是如果约束方程中的某些未知参数,不出现在误差方程中,这时求得的法方程是秩亏的,其逆阵不存在,因此附有约束条件的间接平差的传统模型已不能适用于这种情况的计算。

为此,本文提出了附有约束条件的间接平差扩展模型,给出了该平差模型的函数公式,用两种方法详细推导了该模型的参数估计及其精度评定公式,并用实际观测数据进行了验证,最后得出一些有意义的结论。

1 附有约束条件的间接平差扩展模型

条件平差:利用观测值之间的r个几何条件建立的条件方程为函数模型的平差方法。其线性化后的函数模型为:

式(1)中,B为r×m阶系数矩阵,rank(B)=r,r为条件方程的个数,V为m维观测值的改正数向量,WB为r维闭合差向量。

按最小二乘准则原理,可求得改正数向量和法方程为:

式中,P为观测量权阵;Q为观测量协因数阵;K为r维联系数向量。

间接平差:n个观测值分别表示成t个相互独立参数的函数而建立的函数模型,间接平差的误差方程式为:

式(4)中,V为n维观测量改正数;A为n×t阶未知参数的系数阵,rank(A)=t;为t维独立未知参数向量;l为n维观测值常数向量。利用最小二乘原理求得参数向量的解为:

附有约束条件的间接平差:如果误差方程式(4)的未知参数满足如下s个线性约束条件:

式(6)中,CΤ为s×1阶行满秩阵;wx为s维约束方程常数向量。则在最小二乘准则下,由式(4)和式(6)求得的法方程为:

式中,N=ATPA。

解式(7)得到未知参数向量的解为:

附有约束条件的间接平差的扩展模型:设在附有约束条件的间接平差中未知参数的个数为u,其中误差方程中有t个独立的未知参数向量x1,而约束方程中出现的u-t个未知参数x2不包含在误差方程中,则式(4)可以表示为:

其中,法方程矩阵为秩亏矩阵,其逆阵不存在,因此附有约束条件的间接平差的传统模型已不能适用于这种情况的计算。

为了克服法方程的不可逆,在附加约束条件时对式(9)的函数模型可以改写为:

式中,A为n×t阶系数矩阵;为t维未知参数向量;C1Τ为s×t阶未知参数向量系数矩阵;C2Τ为s×(u-t)阶未知参数向量系数矩阵;为u-t维未知参数向量。

在式(10)中方程的个数为n+s个,未知参数的个数为u个。当n+s>u时,可利用最小二乘原理求解未知参数及其协方差阵。

对式(10)中未知参数估计可以采取以下两种方法。

解法一:对式(10)分两步解。设P为观测量权阵,∑为观测量的方差协方差矩阵,按照间接平差公式求出式(10)中第一式参数和方差协方差阵为:

引入式(10)中第二式的条件后,设的改正数为V1,则式(10)中第二式的方程可写为:

设V'=C1ΤV1,则式(13)变为:

根据的协方差阵式(12),求得式(14)的协方差阵为:

按间接平差求得式(14)中为:

将式(11)代入式(16)中得:

将式(19)代入式(13)得:

按条件平差公式(2)求得式(21)的改正数V1为:

从而求得改正后的为:

解法二:采用附有约束条件的间接平差求解:

根据求条件极值的理论,组成函数:

对K、求偏导数并令其为零,则有:

解式(26)求得K为:

式(27)乘以C2同时减去式(26)中的第二式,可得:

将式(28)代入式(27)中得:

将式(29)代入式(26)的第一式后得:

从式(23)、(30)和式(19)、(28)可知,采用两种方法推导的该平差模型的参数估计公式是完全相同的。

同时根据式(19)和式(23)可得的方差协方差矩阵为:

以上采用两种方法详细推导了附有约束条件的间接平差扩展模型的参数估计及其精度评定公式。

2 实例计算与分析

图1为控制网观测示意图,中间的圆形物体为观测物体,T_A、T_B、T_C、T_D点为布设在该观测物体周围的GPS控制点,GPS控制点采用Trimble双频GNSS接收机观测两个时段,每个时段至少观测12 h。以该物体周围的四个控制点为基准,当该观测物体旋转过程中,使用不少于3台TCA2003全站仪同时观测该物体上的标志点,获取该标志点在地面网中的常规观测量(水平方向、垂直角和边长)。其中水平方向和垂直角的观测精度为0.7″,测距精度为1mm+1ppm。同时这些标志点的运动轨迹满足一空间平面方程和圆球方程,如图2所示。

根据地面网观测量可以列出其在空间直角坐标系中的误差方程[9,10],误差方程中包含了水平方向定向角、大气折光和未知坐标参数。

根据标志点的坐标满足平面方程和球面方程的条件建立约束方程,在该约束方程中,除了包括地面网误差方程中的未知坐标参数外,还包括其他的未知参数,如平面方程的法向量系数、圆球的中心坐标和半径,很显然采用传统的附有约束条件的间接平差模型无法解算,因此采用本文提出的附有约束条件的间接平差扩展模型可以获得所有未知参数估计值及其精度。

四个控制点的GPS观测数据可以采用GPS数据处理软件解算,解算后获取4个控制点在空间直角坐标系中的坐标。以这4个控制点为基准,获取标志点在三维空间直角坐标系中的坐标。具体如下:

地面网中n个标志点的水平方向、垂直角和边长观测量的误差方程的列举具体可参考文献[9,10]。这3类观测量的误差方程分别可以表示为:

式中,Lij为水平方向值,VLij为测站i至测点j水平方向改正数,为测点j在空间直角坐标系中的改正数;(Nj,Ej,Uj)为测点j在以测站i为原点的站心坐标系中坐标,Nj为北方向坐标,Ej为东方向坐标,Uj为高程值,(Nj0,Ej0,Uj0)为其近似值;zi为测站i的定向角,为测站i定向角参数改正数。

式中,βij为测站i至测点j的垂直角,为其改正数;Ra为地球平均曲率半径;K0为大气折光近似值,为其改正数。

式中,Sij为测站i至测点j的边长观测值,为其改正数。

式(33)、(34)和(35)以矩阵的形式可以表示为:

式中,观测量改正数V由水平方向、垂直角和边长观测量的改正数组成,A是未知参数的系数矩阵;为t维独立未知参数向量改正数,,l为误差方程的常数项。

根据图2可知,每个观测标志点坐标同时满足一个平面方程和一个球面方程。其平面和圆球方程可以写为:

线性化式(37)并以矩阵的形式表示为:

其中:

式中,;w为约束方程的常数项;(a,b,c)为平面方程法向量系数,(a0,b0,c0)为其近似值;(Xc,Yc,Zc)为平面圆中心坐标,(Xc0,Yc0,Zc0)为其近似值;r为平面圆半径,r0为其近似值。

联立式(36)和式(38),可得:

综上可以看出,式(42)和式(10)的形式完全相同,因此可以按照本文介绍的附有限制条件的间接平差扩展模型求得未知参数的估计值及其精度。在解算过程中,水平方向和垂直角与边长之间的权估计采用赫尔默特方差分量估计法。具体的成果如表1~表3所示。

表1给出了标志点的坐标及其中误差;表2给出了平面方程中参数的估值及其中误差;表3给出了圆球方程中参数的估值及其中误差。

以上给出了附有约束条件的间接平差扩展模型的应用,利用该模型的解算出平差方程和圆球方程的参数值及其中误差。

3 结论

本文首先给出了条件平差、间接平差和附有约束条件的间接平差模型的函数模型和平差公式,在此基础上,提出了附有约束条件的间接平差的扩展模型,采用两种方法详细推导了该模型的计算公式,最后采用实例数据验证了该模型的应用,并得出以下结论:

(1)采用两种方式推导的附有约束条件的间接平差扩展模型的参数估计及其精度评定公式完全一致。

(2)在附有约束条件的间接平差中,约束方程的未知参数都会出现在误差方程中。而在附有约束条件的间接平差的扩展模型中,部分未知参数不会出现在误差方程中,因此该模型适用于求解出现在约束方程而不出现在误差方程中的未知参数。

摘要：在条件平差、间接平差、附有约束条件的间接平差模型的基础之上,对附有约束条件的间接平差函数模型进行了改进,提出了附有约束条件的间接平差扩展模型,以便解算存在于约束方程而不出现在误差方程中的未知参数,并用两种方法详细推导了该模型的参数估计及其精度评定公式,并证实了两种推导方法的一致性,最后用实际数据验证了该模型的有效性。

关键词：条件平差,间接平差,约束条件,扩展模型,参数估计,精度评定

参考文献

[1]王穗辉.误差处理与测量平差[M].上海:同济大学出版社,2009:1~166.Wang Suihui.Error theory and surveying adjustment[M].Shanghai:Tongji University Press,2009:1~166.(in Chinese)

[2]武汉大学测绘学院测量平差学科组.误差理论与测量平差基础(第二版)[M].武汉:武汉大学出版社,2009:1~178.Surveying Adjustment Subject Group of School of Geodesy and Geomatics,Wuhan University.Error theory and foundation of surveying adjustment(The second edition)[M].Wuhan:Wuhan University Press,2009:1~178.(in Chinese)

[3]黄维彬.近代测量平差理论应用[M].北京:解放军出版社,1992:1~43.Huang Weibin.Modern surveying adjustment theory and its application[M].Beijing:The People’s Liberation Army Press,1992:1~43.(in Chinese)

[4]崔希璋,於宗俦,陶本藻等.广义测量平差(第二版)[M].北京:测绘出版社,1992:46~84.Cui Xizhang,Yu Zongtao,Tao Benzao et al.Generalzied surveying adjustment(the second edition)[M].Beijing:Surveying and Mapping Press,1992:46~84.(in Chinese)

[5]陶本藻.现代平差模型及其应用[J].南京信息工程大学学报(自然科学版),2009,1(1):27~31.Tao Benzao.Modern Adjustment Model and its Application[J].Journal of Nanjing University of Information Science and Technology(Natural Science Edition),2009,1(1):27~31.(in Chinese)

[6]王新洲.广义平差的概括模型[J].武汉测绘科技大学学报,2000,25(3):257~260.Wang Xinzhou.The comprehensive model of generalized surveying adjustment[J].Journal of Wuhan Technical University of Surveying and Mapping,2000,25(3):257~260.(in Chinese)

[7]陶本藻.测量数据处理的统计理论与方法[M].北京:测绘出版社,2007:134~190.Tao Benzao.Stastistical theory and method of survey data procession[M].Beijing:Surveying and Mapping Press,2007:134~190.(in Chinese)

[8]沈云中,陶本藻.实用测量数据处理方法(第二版)[M].北京:测绘出版社,2012:1~20.Shen Yunzhong,Tao Benzao.Practical methods for surveying data processing(The second edition)[M].Beijing:Surveying and Mapping Press,2012:1~20.(in Chinese)

[9]沈云中,陈廷武.上海天文台并址站的空间归心测量[J].同济大学学报(自然科学版),2006,34(2):217~222.Shen Yunzhong,Chen Tingwu.Determination of space coordinate differences of co-location sites in Shanghai observatory[J].Journal of Tongji University(Nature Science),2006,34(2):217~222.(in Chinese)

参数约束篇5

全局最优化在许多领域有广泛的应用, 如计算机科学, 经济管理, 资源管理, 工程设计, 生物工程等.自70年代以来有关全局最优化的新的理论及计算方法层出不穷.人们已提出的有效全局优化方法可以分成两类:确定性方法, 如打洞函数法[3,4,5,10], 填充函数法[1,2]等.不确定方法, 又称随机类方法, 如模拟退火算法[11], 遗传算法[12]等.因此, 研究全局最优化方法, 既具有十分重要的理论意义, 又具有广泛的直接应用前景.

填充函数法最早是由西安交通大学的葛仁傅教授在文章[1]中首先提出的, 它的基本思想是:在目标函数f (x) 的当前的局部极小点x*1处构造填充函数P (x) , 如果P (x) 的一个极小点 $\bar{x}$ 在比x*1所在盆谷更低的盆谷中, 则以 $\bar{x}$ 为初始点极小化f (x) 可得到f (x) 的比x*1更低的极小点x*2, 再用x*2代替x*1可找到更低的极小值点.重复以上过程, 在一定条件下结束.最终可以找到f (x) 的全局极小点x*g.

文献[7,8,9]针对两个参数不易调节的问题利用文献[1]的定义加以修改, 构造出了单参数的填充函数.文中是在以上文献的基础上给出的一种有效的单参数填充函数.

本文的结构如下:第2节是预备知识, 给出填充函数的定义以及一些假设条件;第3节给出一个单参数的填充函数, 通过证明它的性质说明所给的函数是一个填充函数;第4节给出填充函数的算法, 并用数值试验结果来说明算法的有效性和可行性;第5节给出结论.

2 预备知识

我们考虑如下无约束最优化问题:

${\begin{cases} \min f (x), \\ s . t x \in R^{n} . \end{cases} (2.1)$

首先我们做如下假设:

假设1f (x) 在Rn上连续可微, f (x) 的局部极小点个数可以有无限个, 但其局部极小值个数只有有限个.

假设2f (x) 是一个强制函数, 即当x→+∞时, 有f (x) →+∞.

显然, 由假设2可知, 存在这样一个强紧集Ω⊂Rn, 它的内部包含f (x) 的所有极小点.为方便起见, 设Ω被一些常数ci, di, i=1, …, n所确定, 特别, 不妨令Ω={ (x1, …, xn) |ci≤xi≤di, i=1, …, n}, 这里ci, di, i=1, …, n是常数.所以, 原极小化问题等价于如下的极小化问题:

${\begin{cases} \min f (x), \\ s . t x \in Ω . \end{cases} (2.2)$

接下来我们介绍几个概念.

定义2.1 一个连通区域B*称为f (x) 在孤立局部极小点x*的盆谷, 是指x*∈B*, 而且从B*内任一点出发, f (x) 的最速下降轨迹一定趋向于x*, 但从B*外的任一点出发, f (x) 的最速下降轨迹一定不趋向于x*.类似地, 称B*为f (x) 的孤立局部极大点x*处的山丘, 若B*为-f (x) 在x*处的盆.

如果f (x) 的两个局部极小点x*1和x*2处的函数值满足f (x*1) ≤f (x*2) , 称x*1处的盆比x*2处的盆低;否则称x*1处的盆比x*2处的盆高.显然有这样的结论:如果B*是x*的盆谷, 那么对∀x∈B*且x≠x*, 有f (x) >f (x*) .

定义2.2 设x*是f (x) 的一个局部极小点, x*处的简单盆S*是一个含于B*内的一个连通区域, 对∀x∈S*且x≠x*, 不等式 (x-x*) T∇f (x) >0成立, 类似地可以定义简单山丘.

定义2.3F (x, q) 称为极小化问题 (2.1) 的对应于局部极小点x*处一个填充函数, 如果F (x, q) 满足如下性质:

(1) x*是F (x, q) 的一个严格局部极大点;

(2) ∇F (x, q) ≠0, 当x∈S1时, S1={x|f (x) ≥f (x*) , x∈Ω＼x*};

(3) 如果x*不是f (x) 的一个全局极小点, 且S2={x|f (x) <f (x*) , x∈Ω}≠Ø, 那么存在一个点 $\bar{x}$ ∈S2是F (x, q) 的一个局部极小点.

3 填充函数的构造及其性质

我们构造在f (x) 当前局部极小点x*处的填充函数的形式如下:

F (x, q) =-ln (f (x) -f (x*) +1) -q (x-x*) 2. (3.1)

这里, x*是f (x) 当前局部极小点, q>0且足够大.

对任意的x∈Ω, 令d (x) =x-x*;L1=max∇f (x) , 其中x∈Ω;D=minf (x) -f (x*) , 其中x∉S1, S1是x*包含的简单盆.

下面我们将证明F (x, q) 满足定义2.3.

定理3.1 假设x*是函数f (x) 的一个局部极小点, 则x*一定是F (x, q) 的一个严格局部极大点.

证明因为x*是函数f (x) 的一个局部极小点, 所以存在x*的一个领域N (x*, δ) (δ>0) 使得对任意的x∈N (x*, δ) 都有f (x) ≥f (x*) .故, 对任意的x∈N (x*, δ) /{x*}都有

F (x, q) =-ln (f (x) -f (x*) +1) -q (x-x*) 2<0=F (x*, q) . (3.1)

所以, x*是函数F (x, q) 的一个严格局部极大点.

定理3.2 对任意的d∈Ω, 且f (x) >f (x*) , 若dT∇f (x) ≥0, dT (x-x*) >0, 或dT∇f (x) >0, dT (x-x*) ≥0, 那么, d都是F (x, q) 在x*处的一个下降方向.特别沿x-x*方向, F (x, q) 是下降的.

证明由 (3.1) 可知:

$\begin{array}{l} d^{Τ} \nabla F (x, q) \\ = - d^{Τ} (\frac{\nabla f (x)}{f (x) - f (x^{*})} + 2 q | x - x^{*} |) \\ = - (\frac{d^{Τ} \nabla f (x)}{f (x) - f (x^{*})} + 2 q d^{Τ} (x - x^{*})) . (3.2) \end{array}$

由定理2的条件可得:dT∇F (x, q) <0, 故d是F (x, q) 在x*处的下降方向.特别地, 当x∈S1时, F (x, q) 沿方向x-x*是下降的, 则S1变成F (x, q) 的一个山丘.

定理3.3 若f (x) >f (x*) , 且dT∇f (x) <0, dT (x-x*) >0时, 有

$q > \frac{- d^{Τ} \nabla f (x)}{2 d^{Τ} (x - x^{*}) (f (x) - f (x^{*}))} = a (x), (3.3)$

则d是F (x, q) 在x*处的下降方向, 特别在x-x*方向成立.

证明只需要证明dT∇F (x, q) <0即可.由 (3.2) 可知, 要使dT∇F (x, q) <0, 只要 $\frac{d^{Τ} \nabla f (x)}{f (x) - f (x^{*})} + 2 q d^{Τ} (x - x^{*}) > 0$ 即可.在题设条件下, 只要 (3.3) 式成立, 就有dT∇F (x, q) <0, 故d是F (x, q) 在x*处的下降方向, 特别在x-x*方向成立.即当 $q > \frac{L_{1}}{D}$ 时, d是F (x, q) 在x*处的下降方向.

定理3.2和定理3.3说明, 在f (x) 的一个极小点x*1的盆S2中, 只要S2比S1高, 即f (x) >f (x*) , 则至少沿方向x-x*, x∈S2, F (x, q) 总是下降的, 所以F (x, q) 在S2中不可能有任何极小点或鞍点, 即∇F (x, q) ≠0.

定理3.4 若f (x) >f (x*) , 且dT∇f (x) <0, dT (x-x*) >0时, 有q<a (x) , 则d是F (x, q) 在x*处的一个上升方向.

证明仿定理3.3的证明可得, 在给定的条件下有dT∇F (x, q) >0.

定理3.5 当 $\frac{L_{1}}{D} < q < a (x) (3.4)$

时, F (x, q) 在比S1低的盆中必有极小点或鞍点.

证明由定理3.3和定理3.4, 可得定理3.5.此外, 在比x*的盆S1低的盆S2中, 当f (x) -f (x*) →0+时, a (x) →+∞, 这时 (3.3) 式右侧趋于正无穷大.这个性质保证, 不管q多大, (3.4) 式都成立, 同时也保证了全局极小点不会丢失.

4 算法和数值试验

求解问题 (2.1) 全局最优解的新的单参数填充函数算法[6]如下:

步0.选取M>0作为q的终止值.选取方向ei, i=1, …, k0和整数k0≥2n, 这里n是变量的个数.选取一个初始点x $_{1}^{0}$ ∈Ω.令k∶=1.

步1.从初始点x $_{k}^{0}$ 出发, 用局部极小化方法得到目标函数f (x) 的一个局部极小点, 记为:x*k.取一个初始参数q∶=q0, 令i=1, δ>0 (δ可适当选取) .

步2.令F (x, q) =-ln (f (x) -f (x*) +1) -q (x-x*) 2. (4.1)

步3.令 $\bar{x}_{k}^{*} = x_{k}^{*} + δ e_{i}$ .如果 $f (\bar{x}_{k}^{*}) < f (x_{k}^{*})$ , 则令 $x_{k + 1}^{0} ∶ = \bar{x}_{k}^{*}, k ∶ = k + 1$ 转步1.

步4.以 $\bar{x}_{k}^{*}$ 为初始点, 用局部极小化方法解极小化问题 (4.1) , 令 $\bar{x}$ q, x*k为所得极小点, 如果 $\bar{x}$ q, x*k∈int Ω, 则令 $x_{k + 1}^{0} ∶ = \bar{x}_{q, x_{k}^{*}}, k ∶ = k + 1$ 转步1;如果 $\bar{x}$ q, x*k∈∂Ω, 则转步5.

步5.如果q<M, 则令i∶=1, 增加q, 转步4.否则, q:=q0, 转步6.

步6.如果i<k0, 令i∶=i+1, 转步3.否则, 停止, x*k已经是极小化问题的一个全局极小点.

我们验算如下算例:

以下算例的局部极小值点都是在Matlab 7.0的工具箱fmincon, Windows XP, Celeron (R) CPU 2.80 GHZ上得到的.每个算例的数值结果都分别用表格给出, 在运算或绘制的表格中我们用到如下记号:

ei, i=1, …, n:其第i个元素为1, 其它元素为0.

k:表示极小化问题 (2.1) 的局部极小化过程的次数.

q:表示用于寻找第k+1个局部极小值点的参数.

x $_{k}^{0}$ :表示原极小化问题 (2.1) 的第k次极小化过程的初始点.

x*k:表示原极小化问题 (2.1) 的第k个极小值点.

f (x*k) :表示原极小化问题 (2.1) 的第k个极小值点处的函数值.

time:表示算法停止时所占用CPU的时间.

数值试验结果如下:

算例1 Goldstein-Price问题

我们取Ω={-3≤xi≤3|i=1, 2}.以 (1, 1) 为初始点, 用上述算法得到全局极小值点为x*= (0.0000, -1.0000) ;对应的全局极小值为f (x*) =3.0000.

算例2 Three-Hump Camel-back问题

f (x) =2x $_{1}^{2}$ -1.05x $_{1}^{4}$ +x $_{1}^{6}$ /6-x1x2+x $_{2}^{2}$ .

我们取Ω={-3≤xi≤3|i=1, 2}.以 (2, 2) 为初始点, 用上述算法得到全局极小值点为x*=1.0e-005 (-0.4384, 0.1846) ;对应的全局极小值为f (x*) =4.9931e-011.

算例3 Six-Hump Camel-back问题

f (x) =4x $_{1}^{2}$ -2.1x $_{1}^{4}$ +x $_{1}^{6}$ /3-x1x2-4x $_{2}^{2}$ +4x $_{2}^{4}$ .

我们取Ω={-3≤xi≤3|i=1, 2}.以 (0, 0) 为初始点, 用上述算法得到全局极小值点为x*= (0.0898, 0.7127) ;对应的全局极小值为f (x*) =-1.0316.

5 结论

本文给出了一个单参数的填充函数, 这个填充函数是满足定义2.3的, 并给出了相应的算法, 并且数值试验表明了这个算法的可行性和有效性.

参考文献

[1]R.P.Ge.A filled function method for finding aglobal minimizer of a function of several varia-bles[J].Mathematical Programming 46 (1990) 191-204.

[2]R.P.Ge, Y.F.Qin.The global convexizedfilled functions for globally optimization[J].Applied Mathematics and Computation 35 (1990) 131-158.

[3]Levy, A.V.and Montalvo, A..The tunnelingalgorithm for the global minimization of func-tions[J].SLAM J.Sci.&Stat.Comput., 1985, 6 (1) , 15-29.

[4]Yao, T..Dynamic Tunneling Algorithm forGlobal Optimization[J].IEEE Transactions onSystems, Man, and Cybernetics, 1989, 19 (5) , 1222-1230.

[5]Barben, J., Protopopescu, V.and Reister, D..TRUST:A deterministic algorithm for globaloptimization[J].Science 276 (1997) , 1094-1097.

[6]吴至友.全局优化的几种确定性方法[D].上海大学, (2003) 70-91.

[7]王忠, 王永军.用于全局优化的一种有效的单参填充函数[J].Journal of Inner MongoliaNormal University (Natural Science Edition) 35 (2006) 308-312.

[8]You-lin Shang, Ding-guo Pu, Ai-ping Jiang.Finding global minimizer with one-parameterfilled function on unconstrained global optimi-zation[J].Applied Mathematics and Computa-tion 191 (2007) 176-182.

[9]Xian Liu.Two new classes of filled functions[J].Applied Mathematics and Computation149 (2004) 577-588.

[10]Cetin, B.C., Barben, J.and Burdick, J.W..Terminal Repeller Unconstrained SubenergyTunneling (TRUST) for Fast Global Optimi-zation[J].Journal of Optimization Theory andApplications, 1993, 77 (1) , 97-126.

[11]VAPNIK VLADIMIR N.The nature of sta-tistical learning theory[M].New York:Springer-Verlag, 1995.

参数约束篇6

现在的大多数商用CAM软件没有考虑加工切削力和被加工表面的形状特征等物理特性。利用这些CAM软件生成的刀具轨迹在加工时会出现切削力突变的现象。切削力的变化会引起变形误差、热变形误差及振动误差。这些误差都直接影响工件成形误差的大小[1]。许多学者为了提高工件加工质量和加工效率, 针对切削参数优化进行了大量研究。翟玉山等[2]以材料去除率为约束来优化二维铣削的进给率;Wang[3]将材料去除率作为控制指标来改变数控程序中的进给率;刘长清等[4]在对加工过程进行仿真并预测切削力的基础上, 利用粒子群优化算法来优化数控程序中的进给率和主轴转速, 达到对加工过程优化和控制的目的;Lim等[5,6]提出了使用加工模拟程序帮助NC编程人员规划走刀路径和选择进给率的思想;Kim等[7]利用仿真获取刀具的有效直径, 通过改变主轴转速来达到以恒速度切削的目的。

笔者在前人研究成果的基础上, 首先建立了球头铣刀主切削力的力学模型, 然后在切削力基本恒定的约束下, 利用线性经验公式对数控加工中的NC代码进行切削进给率的优化, 以达到减小甚至消除加工中切削力突变的目的, 从而减小或者消除切削力突变给机床带来的损坏, 减缓刀具的磨损, 防止刀具折断、扎伤工件, 提高加工件的质量。仿真结果表明, 采用优化以后的NC代码能显著缩短加工时间, 提高加工效率。

1 球头铣刀力学模型

文献[8,9,10]建立了多个切削力力学模型。为简化问题, 本文采用文献[10]中的主切削力模型, 该模型比较简单, 计算所得主切削力的准确程度依赖于单位切削力的经验值, 可以满足本文研究的要求。主切削力作用在铣刀外圆的切线上, 它消耗了机床动力的大部分功率[11]。

1.1 球头铣刀主切削力模型

主切削力模型如图1所示。当apo≤R时, 如图1a所示, 切削力公式为[10]

$\begin{array}{l} F = \frac{a_{e o} f_{z} Ζ p_{2}}{2 π} [\arcsin (1 - \frac{h}{R}) - \\ \arcsin (1 - \frac{a_{p o}}{R})] + \frac{f_{z} Ζ p_{1} h}{π} (1) \end{array}$

$h = R - \sqrt{R^{2} - (a_{e o} / 2)^{2}}$

式中, F为铣削主切削力;aeo为铣削行距;fz为铣刀每齿进给量;Z为铣刀齿数;apo为铣削深度;R为铣刀球面半径;p1、p2分别为不同切削层的单位切削力;h为球头铣刀铣削行距间的铣削残留面积高度。

当apo>R时, 如图1b所示, 切削力公式为[10]

$F = \frac{f_{z} Ζ}{π} [\frac{a_{e o} p_{3}}{2} (\frac{a_{p o}}{R} - 1) + \frac{a_{e o} p_{2}}{2} \arcsin (1 - \frac{h}{R}) + p_{1} h] (2)$

式中, p3为 (apo-R, 0) 切削层之间的单位切削力。

单位切削力p4与切削层厚度的关系为

p4=p5/huav (3)

其中, huav为切削层公称厚度;u为指数, 表示hav对切削层单位面积切削力的影响程度。当切削层公称厚度和公称宽度各为1mm时, p5为切削层单位切削力[10]。

(1) 当R-apo≤z≤0时, hav=aeofz/[Rarccos (1-aeo/R) ], 将hav代入式 (3) 得单位切削力p4:

p4=p3=p5[Rarccos (1-aeo/R) / (aeofz) ]u (4)

(2) 当R-h≤z≤R时, hav=2fz/π, 将hav代入式 (3) 得单位切削力p4:

p4=p1=p5[π/ (2fz) ]u (5)

(3) 当0<z<R-h时, 切削层厚度hav随z的变化而变化。

但由于指数u很小, 因此hav随z变化不大, 则以球头铣刀参与切削的切削层厚度hav的平均值havp来计算单位切削力, p4=p5/huavp。

当apo≤R时, havp=∫ $_{R - a_{p o}}^{R - h}$ hav/ (apo-h) dz, 所以

p4=p2=p5/huavp=p5/[∫R-hR-apohav/ (apo-h) dz]u (6)

当apo>R时, havp=∫ $_{0}^{R - h}$ hav/ (R-h) dz, 所以

p4=p2=p5/huavp=p5/[∫ $_{0}^{R - h}$ hav/ (R-h) dz]u (7)

$h_{a v} = a_{e o} f_{z} / [\sqrt{R^{2} - z^{2}} \arccos (1 - a_{e o} / \sqrt{R^{2} - z^{2}})] (8)$

1.2 刀轴不垂直于被加工曲面时切削力的计算

在用球头铣刀进行铣削的时候, 如果刀轴不垂直于被加工曲面, 如图2所示, 则参与切削的切削刃起点不是刀尖的顶点, 需要对1.1节中的铣削力F进行修正, 设F′=F-ΔF。F′为修正后的主切削力, ΔF为从球头铣刀刀尖开始到参与切削的切削刃起点这一段的切削力。

当hh≤h时

ΔF1=∫ $_{R - h_{h}}^{R}$ (fzZp1/π) dz=fzZp1hh/π (9)

当hh>h时

$\begin{array}{l} Δ F_{2} = \int_{R - h}^{R} \frac{f_{z} Ζ p_{1}}{π} d z + \int_{R - h_{h}}^{R - h} \frac{a_{e o} f_{z} Ζ p_{2}}{2 π \sqrt{R^{2} - z^{2}}} d z = \\ \frac{f_{z} Ζ p_{1} h}{π} + \frac{a_{e o} f_{z} Ζ p_{2}}{2 π} [\arcsin (1 - \frac{h}{R}) - \arcsin (1 - \frac{h_{h}}{R})] (10) \end{array}$

所以修正以后的切削力公式如下:

当apo≤R且hh≤h时

$\begin{array}{l} F^{'} = F - Δ F_{1} = \frac{a_{e o} f_{z} Ζ p_{2}}{2 π} [\arcsin (1 - \frac{h}{R}) - \\ \arcsin (1 - \frac{a_{p o}}{R})] + \frac{f_{z} Ζ p_{1} (h - h_{h})}{π} (11) \end{array}$

当apo≤R且hh>h时

$\begin{array}{l} F^{'} = F - Δ F_{2} = \frac{a_{e o} f_{z} Ζ p_{2}}{2 π} [\arcsin (1 - \frac{h_{h}}{R}) - \\ \arcsin (1 - \frac{a_{p o}}{R})] (12) \end{array}$

当apo>R且hh≤h时

$\begin{array}{l} F^{'} = F - Δ F_{1} = \frac{a_{e o} f_{z} Ζ p_{3}}{2 π} (\frac{a_{p o}}{R} - 1) + \\ \frac{a_{e o} f_{z} Ζ p_{2}}{2 π} \arcsin (1 - \frac{h}{R}) + \frac{f_{z} Ζ p_{1} (h - h_{h})}{π} (13) \end{array}$

当apo>R且hh>h时

$\begin{array}{l} F^{'} = F - Δ F_{1} = \frac{a_{e o} f_{z} Ζ p_{3}}{2 π} (\frac{a_{p o}}{R} - 1) + \\ \frac{a_{e o} f_{z} Ζ p_{2}}{2 π} \arcsin (1 - \frac{h_{h}}{R}) (14) \end{array}$

2 切削力优化公式

文献[12]给出了关于切削进给率和期望切削力之间的一个线性关系式:

$f_{\lim, i} = (F_{\lim} - F_{1, i}) \frac{f_{2} - f_{1}}{F_{2, i} - F_{1, i}} + f_{1} (15)$

i=1, 2, 3, …, C

式中, C为刀位轨迹中刀位点的序号;f1为给定的进给率, f2=2f1 (为了获得线性关系) ;F1, i为第i个刀位点进给率为f1时的切削力;F2, i为第i个刀位点进给率为f2时的切削力;flim, i为经过优化得到的第i个刀位点的进给率;Flim为设定的期望切削力。

该关系式可以使得切削力保持在期望切削力附近。

3 数控加工刀轴控制方式

Erdim等[12]分析了球头铣刀的进给运动, 见图3。从图3中可以看出, 球头铣刀在进行下行斜面铣削加工、上行斜面铣削加工以及平面铣削加工时, 切削厚度是不同的, 因此, 切削力也是不同的。在自由曲面的铣削加工中, 下行斜面铣削加工、上行斜面铣削加工和平面铣削加工状态会交替出现, 从而造成切削力的不断变化, 恶化了切削条件。

在加工自由曲面的过程中, 如果条件允许, 应使刀轴矢量始终垂直于被加工曲面, 刀轴矢量与进给速度矢量夹角为90°。这样, 切削厚度波动范围缩小, 减小了铣削过程中切削力的突变, 使得切削平稳。但是需要注意的是, 当刀轴垂直于工件表面时, 球头铣刀刀尖的线速度为零, 与加工件处于摩擦状态, 这样会增加刀具在切削过程中的发热, 影响加工质量, 所以在优化参数时应考虑将刀轴适当倾斜一个角度, 以改善发热的不良状态。

随着刀具轴线相对于刀具铣削方向倾斜的角度β的增大, 切削力变化的趋势是减小的;但是当倾角β达到15°之后, 随着倾角的增大, 切削力的减小趋势不再明显[13]。

4 切削力基本恒定约束下切削进给率的优化仿真

4.1 加工曲面的选取

选取加工曲面时需要考虑如下问题:①考虑机床刀轴转速的变化, 曲面的曲率变化率不能太大;②曲面的波峰波谷不可以太接近, 以免造成干涉。

本文选取的材料为镁合金, 其参数为p5=25MPa, u=0.19[14]。选取的加工曲面的截面曲线为椭圆, 长半轴长为100mm, 短半轴长为5mm。截选了宽为160mm的上半部分来进行加工仿真。为便于研究, 所选取的整个曲面的曲率较小。

设定的切削参数如下:进给率f=3mm/s, 铣削行距aeo=3mm。选用半径为5mm的球头铣刀, 铣刀齿数Z=3, 选用五轴数控机床, 采用的刀轴控制方式为Pattern Surface。

4.2 仿真结果

设定优化后切削力的期望值为400N, 根据球头铣刀主切削力模型和Master CAM生成的NC代码计算切削力F1, i, 然后将进给率f加倍, 其他切削参数保持不变, 根据式 (1) 和式 (2) 计算切削力F2, i, 并根据式 (15) 来计算优化后的进给率flim, i, 最后将计算出的进给率插入到原来的NC代码中进行优化。优化前后的切削力变化曲线如图4所示。切削力的比较见表1。

从表1中的数据可以看出, 优化前, 当进给率为f时切削力的变化幅值为230.231N, 当进给率为2f时, 切削力的变化幅值为403.644N, 而优化后的切削力的变化幅值仅为6.936N, 是优化前进给率为3mm/s时变化幅值的3%, 是进给率为6mm/s时的变化幅值的1.7%。

图5所示为优化前后的进给率变化曲线, 优化前的进给率恒定, 为3mm/s, 比较进给率可知, 优化后加工效率明显提高。

从图4中可以看出, 在恒定的进给率下, 随着刀具的移动, 切削力先减小、后增大, 造成切削力波动。利用式 (15) 对进给率进行优化后, 在优化前主切削力较小的地方采用较大的进给率, 从而保证了切削力基本维持恒定, 如图4中优化后曲线所示。

本文就该优化方法与普通三轴自由曲面加工和三轴Highfeed优化加工方式的效率作了比较。切削加工参数如下:p5=275MPa, u=0.19, 进给率f=3mm/s, 铣削行距aeo=3mm, 采用半径R=5mm的球头铣刀, 铣刀齿数Z=3。三种方式的加工时间见表2, 加工时的切削力如图6所示。

从图6a可以看出, 在切削力基本恒定的约束下, 对NC代码进行进给率优化后, 切削力基本维持不变, 从图6b中可以看出使用Master CAM的Highfeed模块对程序进行优化, 优化后的切削力变化较大, 且在加工中出现2个波峰。

高速切削加工的加工时间比普通三轴加工方式缩短了84%, 但是切削力有较大波动;而本文提出的优化加工方式, 不仅在加工时间方面比普通三轴加工方式缩短了42%, 而且切削力在加工过程中基本保持不变。与普通的三轴加工方式相比, 本文建立的优化加工方式在时间与切削力方面都有很大的优势。

切削力基本恒定约束下自由曲面优化加工考虑了切削加工中实际的物理过程, 可以减小机床因为力的波动而产生的振动以及超负荷, 可以延长刀具的寿命且避免刀具折损, 优于Highfeed加工方式。

5 结束语

针对铣削加工中常用的球头铣刀, 建立了球头铣刀的主切削力模型, 该模型考虑了切削厚度对单位切削力的影响, 并对当参与切削的切削刃的起始点不是球头铣刀球头顶点时切削力的计算进行了修正, 探讨了在切削力基本恒定约束下切削进给率的优化。在切削力基本恒定约束下, 切削进给率的优化减小了加工中切削力的波动, 减小了加工中的设备损坏率和刀具的磨损, 延长了它们的使用寿命, 节约了成本, 同时优化以后显著缩短了加工时间, 提高了生产效率。

摘要：针对球头铣刀的铣削特点, 建立了铣削的力学模型。利用经验公式, 在保持切削力基本恒定的约束下, 对加工自由曲面的NC代码进行了进给率的优化, 从而减小切削力的波动, 提高数控加工件的质量和加工效率, 延长刀具的使用寿命并避免刀具折损等现象。建立了在切削力基本恒定的约束下的加工优化原型系统, 并用实例进行了验证。

参数约束篇7

在我国实施经济发展方式转变的战略过程中, 积极推进产业与资源环境协调发展是减轻生态环境压力, 有效利用资源的必然选择和有效路径。然而, 产业与资源环境之间存在各种胁迫与约束, 是国内外学者研究的重点领域。国外如福雷斯特[1]提出产业环境的概念, 倡导政府在实施产业结构转换或调整过程中充分考虑产业发展与环境保护的相互协调性。赫尔曼·戴利[2]倡导在产业结构调整中尽量减少资源消耗量大的产业的规模, 尽量扩大无污染或污染小的产业。他认为人类具有克服环境资源稀缺的能力。米勒和布莱尔[3]利用投入产出方法分析能源使用和环境问题, 而且对行业活动中的能源投入和污染物排放进行定量分析, 为产业结构调整提供了政策支持。国内如张景云[4]认为沿江地区工业结构不合理是资源供应紧张和环境质量恶化等问题产生的根源。要突破资源环境约束, 只有协调好环境保护与产业发展之间的矛盾, 才能实现可持续发展。此外, 王力[5]、张少兵[6]对产业结构优化的环境, 卓锰钢[7]对产业结构优化的资源, 杨艳琳[8]对产业转型的资源环境进行研究。姚聪丽[9]在理论上分析了资源环境与工业化发展的关系, 在此基础上建立了资源环境下工业化发展的模型。张燕[10]基于产业地位划分法探讨了中国产业升级路径中的资源环境。

这些研究主要是从定性的角度来分析产业与资源环境的关系, 如何定量地揭示绿洲产业与资源环境间的关系, 特别是量化资源环境对产业发展的约束, 研究更少。绿洲是在干旱气候条件、水文地貌、人类活动诸因素综合作用下形成的特殊景观, 是干旱区的精华和人类活动的载体。新疆绿洲占干旱区绿洲总面积的72% , 因此, 以新疆作为研究对象具有一定的代表性。本文采用变参数状态空间模型, 以新疆为例, 对产业发展受资源环境的约束进行了实证研究, 对推进干旱区绿洲产业可持续发展至关重要。

2 可变参数模型的状态空间模型构建

传统变量之间的关系一般采用固定参数模型来描述, 只能反映变量间的静态关系, 难以反映变量间由于时间变化以及各种各样的外界冲击导致的变量间关系的改变程度, 因此, 需要考虑采用可变参数模型。根据状态空间分析的思想, 状态空间分析法是现代控制理论的基础, 是对系统内部以及系统之间输入输出关系的描述, 同时适用于时变系统、非线性系统和随机控制系统, 基于量测方程和状态方程, 以数学模型对系统的一种完全描述。产业结构受资源环境的约束也是随时间发生变化的, 因此, 为测度产业结构受资源环境的约束程度, 利用状态空间模型 ( Harry[11], Hamilton[12]) 构造的可变参数模型如下:

方程 ( 1) ( 2) ( 3) 就是状态空间模型。状态空间模型是动态模型的一般形式, 由量测方程和状态方程构成。量测方程 ( 1) 中, xt是具有随机系数的解释变量的集合, 随机系数向量 βt是状态向量, 称为可变参数。向量 μt, εt是量测方程和状态方程的扰动项, 假定为相互独立的, 且服从均值为0, σ2方差为和协方差矩阵为Q的正态分布。Rt是系统矩阵, 能随时间改变, 但是都是可以预先确定。βt是不可观测变量, 必须利用可观测变量yt和xt来估计。状态向量方程一般采用卡尔曼滤波 ( Kalman Filtering) 的迭代算法进行最优值估计, 当新的观测值得到, 可以利用Kalman滤波连续地修正状态向量的估计值, 体现出状态空间模型估计参数的时变性。

根据变参数状态空间模型量测方程和状态方程, 得到预测方程:

一旦得到新的预测值yt , 就能够修正 βt的估计at / t - 1, 更新方程是:

给出一步向前状态条件均值, 我们还可以得到的一步向前最小均方误差估计:

一步向前预测误差可以通过下面的公式得到:

预测误差的方差被定义为:

Kalman滤波提供了状态向量的最优估计。当所有的T个观测值都已处理, Kalman滤波基于信息集YT, 产生当前状态向量和下一时间期间状态向量的最优估计。这个估计包含了产生未来状态向量和未来观测值的最优预测所需的所有信息。

3 新疆产业发展的资源环境约束实证检验

3. 1 变量选取与数据说明

样本研究期间确定为1995 - 2009 年。数据来自《新疆统计年鉴》、《新疆年鉴》相关年份。产业发展很大程度体现为产业结构的不断优化, 因此, 本文选择产业结构优化率 ( ISR) 反映产业发展水平, 以第二产业产值与第三产业产值的和占GDP的比重表示, 用以表征因变量, 工业废气排放 ( FQ) 、工业废水排放 ( FS) 、工业固体废物排放量 ( GF) 、一产能源效率 ( YNX) 、二产能源效率 ( ENX) 、三产能源效率 ( SNX) 、农业水耗 ( NHS) 、非农业水耗 ( FNHS) 反映资源环境影响情况, 是随时间变化的自变量。其中产业能源效率为产业产值与产业能源消费之比。为了保证数据的可比性和容易得到平稳序列, 同时削弱可能的异方差, 对数据取自然对数处理。

3. 2 平稳及协整检验

要建立产业发展资源环境约束的可变参数状态空间模型, 要求变量是平稳的且变量间存在协整关系。首先, 采用扩展的Dickey - Fuller ( ADF) 的检验方法对变量进行单位根检验, 利用AIC与SC准则确定变量的滞后阶数。检验结果如表1 所示。

注: △代表一阶差分; 括号内前两个字符表示检验的类型 ( c: 含常数项, t: 含趋势项, ) , 第3 个字符表示滞后的阶数

从表1 可以看出, 对变量的对数取一阶差分后为平稳序列, 均在5% 的水平上显著, 通过ADF检验, 即各变量均为I ( 1) 阶单整序列, 为了研究平稳时间序列间的长期关系, 需要确定它们之间是否具有协整关系。我们选取适合多变量协整检验的Jo-hansen协整法对变量进行协整分析, 协整检验结果如表2。

注: * 表示5% 显著水平下的临界值, 以迹检验结果为准

从表2 可以看出, 产业发展水平ISR与二产能耗ENH的统计量19. 23 大于5% 显著水平下的临界值18. 39, 拒绝没有协整方程的原假设, 5. 87 大于5% 显著水平下的临界值3. 84, 拒绝至多有1 个协整方程的原假设, 认为至少存在2 个协整方程, 因此, 产业发展水平ISR与二产能耗ENH存在2 个协整关系。同理, 产业发展水平ISR与其他变量存在1 个协整关系, 故产业发展水平ISR与各变量间存在长期的动态均衡关系, 可进行变参数状态空间模型的估计。

3. 3 时变参数模型估计结果

根据可变参数状态空间模型, 代入新疆产业发展与资源环境相关数据估计结果如下:

式 ( 12) 表示新疆产业发展可变参数状态空间模型的量测方程, 下方括号中的数字表示Z统计量, βit分别代表资源环境变量对产业发展的动态约束关系。

表3 为资源环境各解释变量时变系数变动统计描述。资源环境对产业发展的动态约束关系如图1、2、3 所示。

图1 反映了生态环境对产业发展的约束路径。结合表3 和图1, 环境对于新疆产业发展的约束自1995 年以来, 工业废气排放对产业发展的约束一直最强, 工业固废排放对新疆产业发展的约束居其次, 工业废水排放对新疆产业发展的约束较弱, 波动幅度较大。从产业发展约束的各生态因子来看, 工业废气排放对产业发展的约束从1995 - 1997 年迅速增加, 于1997 年达到极大值0. 1445, 1998 - 2009 年间对产业结构的约束较稳定, 波动性最小, 标准差仅为0. 0209; 工业固废排放对产业发展的约束1996年增大, 1997 年约束减小, 为负值。1998 年以后, 工业固废对产业发展的约束较为稳定; 工业废水排放对产业发展的约束较弱, 自1997 年以后都为负值, 并于2003 年达到最小值0. 3611, 总体来看工业废水排放对产业发展的约束波动性最大, 标准差达到0. 1347, 且有逐渐增大的趋势。

从图2 可以看出, “九五”以来, 三产能源效率对产业发展的约束波动较大。1999 年以前, 一产能源效率对新疆产业发展的约束最大, 第三产业能源效率约束程度其次, 二产能源效率约束较弱;2000 年以后, 二产能源效率约束程度迅速增加, 一产能源效率约束程度逐渐减弱。1999 年以前, 三次产业能源效率对产业发展的约束程度较为接近, 1999 年以后, 一产能源效率对产业发展的约束程度变动较平稳, 且有逐渐减弱的趋势; 二产能源效率对产业发展的约束程度自1999 年迅速增强, 2002 年达到极大值后约束程度减弱, 2005 年后保持较高程度的约束水平, 标准差较大; 三产能源效率对产业发展的约束程度1999 年以前为正值, 2000 年后对产业结构的约束程度减弱, 2000 - 2007 年均为负值, 2001 年达到最低约束水平, 2008 年以来约束程度又有增大的趋势。

由图3 可知, 农业耗水与非农业耗水对新疆产业发展约束程度的变化较稳定, 且非农业耗水对产业发展的约束程度较小。结合表3 分析可知, 非农业耗水对产业发展的约束程度除1996 年外, 均为正值, 意味着产业结构的优化与非农业耗水成同向变动的关系。由图3 可知, 非农业耗水对产业发展的约束程度较小, 且保持较平稳的趋势; 而农业耗水与产业发展成反向变动的关系, 农业耗水约束系数较大, 即对产业发展的约束程度较大。产业结构的优化按照国际经验表现为二、三产业所占比重的提高, 新疆农业耗水占三产用水总量的90% 以上, 意味着产业结构越优化受到水资源的约束就越小。

4 结论与建议

4. 1 结论

本文采用可变参数状态空间模型, 利用卡尔曼滤波方法估计并检验新疆产业发展的资源环境约束的变参数关系, 分析结论如下:

( 1) 通过生态环境对产业发展约束的数据分析表明, 重点选取的环境因子中工业废气排放对产业发展的约束最强, 工业废水排放对新疆产业发展的约束较弱。从工业废气排放的行业来看, 电力蒸汽热水生产和供应业、非金属矿物制造业、黑色金属冶炼及压延加工业、石油加工及炼焦业、采选业是新疆向大气排放污染物的重点行业, 以上5 个行业占全行业工业废气污染物排放量的70% 以上, 其中电力蒸汽热水生产和供应业的排污量在全行业中一直位居第1。说明新疆自 “九五”以来产业结构以资源型产业为主, 从全疆排放工业废水的主要行业来看, “八五”末排放量最大的是纺织业、造纸及纸制品业, 到 “十一五”末工业废水排放量最大的是化学纤维制造业、化学原料及化学制品、食品制造业。由于产业结构的调整, 化工业得到迅速发展, 减少了废水的排放。全疆工业固体废物产生量最大的行业是采选业, 其次是电力蒸汽热水的生产和供应业、黑色金属冶炼及压延加工业, 利用率最高的是黑色金属冶炼加工业。从约束力度来看, 一直保持较平稳的趋势, 说明新疆产业结构的调整初显成效, 有很大的优化空间。

( 2) 通过资源对产业发展约束的数据分析表明, 从产业结构调整的能源约束来看, 2000 年以后, 二产能源效率约束程度在逐渐增加, 一产能源效率约束程度逐渐减弱, 三产能源效率约束程度较小。正是由于新疆新型工业化的推进, 优势资源转换战略的实施, 以石油、天然气、煤炭为主的产业结构占据了主导地位, 能源技术效率低下, 能源的约束力度较大。从产业结构调整的水资源约束来看, 农业耗水与非农业耗水对产业发展的约束程度的变化较稳定, 且农业耗水对产业发展的约束程度较大。

4. 2 建议

产业是资源的转换器, 并决定了对生态环境的胁迫, 因此, 资源环境对产业发展的约束, 最终决定于产业的发展水平, 不断合理化、高度化的产业结构是突破资源环境瓶颈的重要路径。对于新疆来讲, 首先, 加快推进新疆高效环保的现代农牧业发展, 继续推进新疆农牧特色产品生产向农牧产品流通和精深加工转变, 同时, 加快特色农牧业和绿色、有机农牧产品基地建设, 发展现代农牧业, 推进新疆农牧业由大变强的优化升级。其次, 新疆以资源为主的产业具有明显的超前发展优势, 应运用信息技术和高新技术, 提高原油和天然气开采、加工与石化后续产业的结合程度, 延长加工工艺流程和产业链, 提高附加值, 按照资源节约, 保护环境的要求, 推进节能减排, 促进新疆传统重化工业向低消耗、轻污染、高素质产业发展。新疆第三产业应重点发展金融、信息、物流、法律、咨询等生产性服务业, 并促进服务业与制造业的渗透与融合。利用新疆的地缘优势, 大力发展现代物流业, 加快构建高效、快捷、通畅、具有较强竞争力的现代物流服务体系, 打造西部物流基地或区域中心, 实现产业结构高度化。另外, 新疆应加快培育和发展战略性新兴产业, 从实际出发, 有计划地培育和发展风电产业、煤炭资源清洁利用产业, 节能环保产业、生物医药, 还有符合新型工业化要求的环保、安全生产等新兴产业, 打造环保产业体系, 逐步形成新的经济增长点。

参考文献

[1]福雷斯特.世界动态学[C].中国环境管理、经济和法学学会, 1984

[2]戴利.超越增长:可持续发展的经济学[M].上海:上海译文出版社, 2001

[3]MILLER R E, BLAIR P D.Input-output Analysis:Foundations and Extensions[J].Englewood Cliffs, 1985, 200-227

[4]张景云.江苏省沿江地区工业可持续发展对策研究——基于资源环境约束条件的分析[D].东南大学, 2006

[5]王力.环境影响下江苏省产业结构优化研究[D].江南大学, 2008

[6]张少兵.环境影响下区域产业结构优化升级研究:以长三角为例[D].华中农业大学, 2008

[7]卓锰钢.基于节能减排的浙江产业结构优化研究[D].浙江工业大学, 2008

[8]杨艳琳.中国中部地区资源环境影响与产业转型研究[J].学习与探索, 2010 (3) :154-157

[9]姚聪莉.资源环境影响下的中国新型工业化道路研究[D].西北大学, 2009

[10]张燕, 路文杰.中国产业升级路径中的资源环境影响——基于产业地位划分法的研究[J].河北科技大学学报:社会科学版, 2010, 10 (3) :11-15

[11]HARRY.Foresting Structural Time Series Models and the Kalman Filter[M].Cambridge University Press, 1999