非线性最小二乘

2024-10-10

非线性最小二乘（共9篇）

非线性最小二乘篇1

1. 分析方法

1.1 线性最小二乘

对于线性回归模型[1]

其中y是因变量向量, x是自变量矩阵, β0, β1是未知参数, ε为随机误差。参数β0, β1一般采用最小二乘法进行估计, 求使达到极小的β0, β1的估计分别为, 则

为所求的回归值, 亦称预测值。

1.2 非线性最小二乘

在实际问题中存在着大量的非线性最小二乘问题。一般地, 非线性观测方程[2]可写为

相应的误差方程为, 残差平方和为

使得上式达到最小的称为得一个非线性最小二乘估计。因为非线性模型比较复杂, 其最小二乘估计往往难以求得, 有些采用数值迭代方法, 往往带有一定的误差。本文讨论用经典的线性最小二乘的方法求非线性最小二乘估计。

2. 案例应用

表1是患者细胞呈现的某种基因的观测值[3], 已知中数据适合模型, a是未知参数。由于a不仅仅是与相乘, 这个问题应该用非线性最小二乘法来解决。

将函数求导, 令各参数对应的偏导数为零, 求出极值点, 即最小二乘估计, 代入模型之后进行预测。绘制图形, 其中*点为数据点, 曲线为回归线, 并比较真实值和模型预测值。

另一种方法是将数据转换为线性最小二乘法问题, 构造新的测量值, 通过模型找到X的函数z。用线性最小二乘法估计变量a, 绘制图形并与非线性最小二乘法结果比较。令

将z1与X进行一元线性回归分析 (过原点) ,

结果为即, 参数a的置信区间为 (1.8656, 2.1928) 。

从图形看用线性最小二乘法所得回归线与真实值近似程度较好。从具体误差数据来看, 也是如此。

3. 结论

本文讨论了将非线性问题转化为线性问题, 通过一个实际案例两个方法进行比较, 说明线性最小二乘在某些情况下优于非线性最小二乘。

参考文献

[1]王松桂, 陈敏, 陈立萍.线性统计模型线性回归与方差分析[M].北京:高等教育出版社, 1999:29-30.

[2]王新洲.非线性模型参数估计理论与应用[M].武汉:武汉大学出版社, 2002.

[3]Frederick R.Adler.Modeling the Dynamics of Life[M].Thomson Brooks/Cole, 2001.

非线性最小二乘篇2

在线性回归模型中普通的最小二乘估计(LSE)许多情形下是不稳健的.本文介绍了一种投影深度函数,深度加权平均和深度加权LSE,这些估计量有符合需要的`稳健性.并讨论了在深度加权LSE情形下线性回归模型的拟合检验问题.

作者：范允征林路 Fan Yunzheng kin Lu 作者单位：范允征,Fan Yunzheng(南通大学理学院,江苏,南通,226007)

林路,kin Lu(山东大学数学与系统科学学院,山东,济南,250100)

非线性最小二乘篇3

总体最小二乘[1]认为平差模型中观测向量和系数矩阵同时含有误差, 并对其进行改正。针对系数矩阵含有误差的平差问题, 在理论上总体最小二乘较最小二乘更为严密。总体最小二乘的经典解法是奇异值分解法 (SVD) , 但由于其计算复杂且不利于编程, 测量学者提出了一些基于平差模型的迭代解法[2~7], 并对其实际应用作了一些研究[8~10]。在测量数据处理中, 线性拟合是总体最小二乘一个最重要的应用, 但总体最小二乘的常规解法 (即奇异值分解法和迭代解法) 都无法顾及线性回归中系数矩阵的常数列而是对其进行了改正, 这对线性回归模型来讲是不合理的。文献[5]以直线必须通过数据点的中心, 才能使其偏差最小为约束条件, 推导了一种求解线性模型的总体最小二乘算法, 虽然该算法提高了线性回归的逼近精度, 但仍然存在偏差。文献[7]通过将数据中心化后再采用奇异值分解法对回归参数进行解算, 则分离了系数矩阵的常数列, 得到的参数估值合理可靠。但其缺陷是计算复杂难以编程实现。鉴于此, 本文在基于平差模型的基础上, 推导了一种求解线性回归参数的总体最小二乘算法, 该算法易于编程实现。算例分析表明其可靠、合理。

1 总体最小二乘常规解法

1.1 SVD解法

线性回归模型为[1,2,6]:

式中:Δ观测量y的真误差, ak为回归参数, 当同时考虑到自变量xk也含有误差, 此时系数矩阵A中的元素含有误差, 则式 (2) 变换为总体最小二乘的平差模型:

采用奇异值分解 (SVD) 法将式 (2) 变换为:

求得参数估值为[1,2,6]:

式 (4) 中, n为非固定列对应参数的个数, σ2n+1为增广矩阵[A L]的最小特征值。其单位权中误差计算公式为σ02=σ2n+1/ (m-n-1) , 其中m为观测个数, n+1为参数的个数。

SVD解是将系数矩阵中所有元素看成是含有误差的, 这对线性回归模型来讲是不合理的。而且SVD解法较为复杂, 不利于测量数据处理的编程实现。

1.2 常规迭代解法

总体最小二乘的迭代解法为[2]:

进行迭代解算时, 首先按最小二乘法计算参数的初值X (0) , 然后根据式 (5) 反复迭代直到‖X (i+1) -X (i) ‖<ε则停止计算输出参数值, 单位权中误差按σ0=V/ (m-n) 计算。迭代解法最大的优点是充分运用到测量平差模型, 其迭代格式简单易于编程。但常规的总体最小二乘迭代解法都是将系数矩阵所有元素看成是含误差的, 这对线性回归模型也是不合理的。

2 线性回归参数的总体最小二乘迭代解法

2.1 算法推导

根据前文所述, 有必要建立基于总体最小二乘的线性回归模型和解法, 线性回归数学模型为:

可将式 (6) 进行等价转换为:

表示为矩阵形式为:

设v=vec (EA) , vec表示向量拉直运算, vec (EA) 即将矩阵从左至右逐列拉直成一列向量。在不考虑权重时, 根据总体最小二乘原理, 相应的误差期望和方差为:

式 (9) 中, 为矩阵的克罗内克积, Im和In分别为m和n阶的单位矩阵, 相应的其平差准则为:

根据式 (8) 与式 (10) 构造拉格朗日目标函数:

式 (11) 中, k为 (n+1) ×1的拉格朗日乘数;, In+1为n+1阶的单位矩阵。根据拉格朗日函数求极值的必要条件, 将式 (11) 分别对v、x求导并令其等于0, 化简整理后可得:

由式 (12) 的第一式可得:

式 (13) 中, vec-1表示vec的逆运算, 即将矩阵重新构造为m× (n+1) 的矩阵。

将式 (12) 的第一式代入到式 (8) 中, 并顾及, 整理后可得:

根据式 (11) 第二式和式 (14) 则可得:

根据式 (15) 即可得到参数^x的表达式:

根据上述推导, 其参数解算步骤如下:

(1) 对回归参数附初值^x (0)

(2) 按式 (17) 计算k值和新的回归参数值

(3) 重复步骤 (2) , 直到‖x (i+1) -x (i) ‖<ε, 则停止迭代

按上述迭代方法即可求得式 (7) 的方程, 相应的线性回归方程的参数解:

2.2 单位权中误差评定

根据线性回归模型式 (7) , 当同时考虑自变量误差即采用总体最小二乘求取回归参数值时, 其单位权中误差的计算公式为:

根据上文的迭代解再结合式 (12) 可得:

则按本文的推导, 其单位权中误差计算公式为:

在计算线性回归单位权中误差时, 按常规的总体最小二乘解法算得的结果比按式 (19) 计算的结果要小, 这与实践是不不相符的。因为, 常规的总体最小二乘解法没有顾及线性回归模型中系数矩阵常数列。而本文推导的解法从理论上比常规总体最小二乘解法更严密, 究其原因, 是因为按常规总体最小二乘解法对线性回归模型的单位权中误差评定有误, 下面则具体进行讨论。

按常规总体最小二乘法将线性回归系数矩阵中的常数列也进行了改正并参与精度评定, 其改正后的模型可表示为:

其单位权中误差评定公式为:

从式 (22) 可以看出, 其将常数列的改正值看成是自变量的改正值, 这显然是不对的。从改正模型来看, 应该将常数列的改正值归入到因变量的改正值中, 则变量的改正值为:

按式 (24) 进行改正值修正后, 再按式 (19) 进行单位权中误差的评定。如此得到的结果才与线性回归模型相符。这也解释了前文提出的疑问, 同时也表明不能单一以单位权中误差来评定一种平差模型和算法的优劣。因为它们的单位权中误差评定方式不一定都相同。

3 算例分析

3.1 算例1

运用文献[6]例5~10的数据, 用三个点 (1, 2) 、 (2, 6) 、 (6, 1) 来拟合一元线性回归方程。运用本文的解法得到的参数估值为6、1。这与文献[5]中将数据中心化后再采用奇异值分解法得到的结果相同, 而按常规总体最小二乘的奇异值分解法和迭代解法得到的参数估值皆为6.74、0.99。这是因为总体最小二乘的常规解法将系数矩阵的常数列也进行了改正, 得到的结果有偏差。而本文的解法则将系数矩阵的常数列分离开不加改正, 计算得到的参数估值较之可靠。

3.2 算例2

为进一步验证本文方法对线性回归模型的统一性和可靠性, 运用MATLAB模拟一平面方程。设平面方程为:z=1.5+x+2y, 在x和y没有误差时求得z的值上加上均值为0, 方差为0.005的随机误差, 组成观测值。然后分别对x和y添加均值为0, 方差为0.03的随机误差, 组成新的观测值, 如表1所示。分别采用最小二乘法 (LS) 、常规总最小二乘法、本文解法解算线性回归方程的参数值, 并按文中所述的单位权中误差评定公式计算其单位权中误差, 结果如表2所示。

从表2中的参数估值结果不难看出, 采用本文的总体最小二乘解法解算的回归参数值与真值最为接近, 而常规的总体最小二乘解法由于将线性回归模型中的系数矩阵中的常数列也进行了改正导致其结果与真实值有偏差, 最小二乘法由于没有考虑到线性回归自变量的误差即系数矩阵元素的误差, 计算的结果与真值偏差最大。这也可以从其单位权中误差中得出, 本文的总体最小二乘解法所得单位权中误差最小, 常规总体最小二乘解法次之, 最小二乘法最小。需要指出的是, 如果按常规总体最小二乘的单位权中误差计算公式得到的单位权中误差为0.0227, 比本文方法得到的值要小, 这就会误以为按常规总体最小二乘解法解算线性回归参数得到的精度比本文要高。其实不然, 这是因为其针对线性回归模型的单位权中误差计算公式有误, 通过本文的纠正得出的值才符合其实际。这也可以通过参数的估值结果与真值的比较中看出。

4 结束语

本文给出了一种线性回归参数估计的总体最小二乘算法, 该算法即能同时考虑线性回归自变量的误差即平差模型系数矩阵元素的误差, 又能顾及线性回归中系数矩阵常数列即对常数项不加改正。这弥补了常规总体最小二乘解法针对线性回归模型解算的不足, 且对基于总体最小二乘的线性回归单位权中误差进行了探讨, 纠正了针对线性回归模型的常规总体最小二乘单位权中误差评定的偏颇。

摘要：根据线性回归模型的特点, 推导了一种解线性回归参数的总体最小二乘算法。并对其单位权中误差的评定进行了探讨, 通过理论推导表明按常规的总体最小二乘解法求得的线性回归单位权中误差与实际不符, 并对其进行了纠正。通过算例分析且与常规的总体最小二乘解法进行比较, 结果表明了算法的正确性和可靠性。

关键词：总体最小二乘,线性回归,平差模型,迭代解法,奇异值分解

参考文献

[1]GOLUB G H, VAN L C F.An Analysis of the Total Least Squares Problem[J].SIAM J Numer.Anal, 1980, 17:883~893.

[2]鲁铁定, 周世健.总体最小二乘的迭代解法[J].武汉大学学报 (信息科学版) , 2010, 35 (11) :1351~1354.

[3]孔建, 姚宜斌, 吴寒.整体最小二乘的迭代解法[J].武汉大学学报 (信息科学版) , 2010, 35 (6) :711~714.

[4]许超钤, 姚宜斌, 张豹等.基于整体最小二乘的参数估计新方法及精度评定[J].测绘通报, 2011, (10) :1~4.

[5]邱卫宁, 齐公玉, 田丰瑞.整体最小二乘求解线性模型的改进算法[J].武汉大学学报 (信息科学版) , 2010, 35 (6) 708~710.

[6]丁士俊, 姜卫平, 杨颜梅.整体最小二乘线性回归模型与算法[J].测绘通报, 2012, (12) :8~10.

[7]邱卫宁, 陶本藻, 姚宜斌等.测量数据处理理论与方法[M].武汉:武汉大学出版社, 2008.

[8]陈义, 陆珏, 郑波.总体最小二乘方法在空间后方交会中的应用[J].武汉大学学报 (信息科学版) , 2008, 33 (12) :1271~1274.

[9]陆珏, 陈义, 郑波.总体最小二乘方法在三维坐标转换中的应用[J].大地测量与地球动力学, 2008, 28 (5) :77~81.

非线性最小二乘篇4

二阶抛物问题的最小二乘扩展混合元格式

对二阶抛物问题提出了最小二乘扩展混合元数值模拟格式.数值算例验证了格式的.有效性.

作者：丁胜 DING Sheng 作者单位：山东师范大学数学科学学院,济南,250014刊名：科学技术与工程 ISTIC英文刊名：SCIENCE TECHNOLOGY AND ENGINEERING年，卷(期)：8(23)分类号：O241.82关键词：二阶抛物问题最小二乘扩展混合有限元方法

非线性最小二乘篇5

电子稳定程序(ESP)已成为现代车辆最重要的主动安全系统之一。ESP通过监测横摆角速度和质心侧偏角两个被控量来调整车辆姿态，从而避免或减轻行驶过程中的失稳现象[1]。其中，横摆角速度通常可以通过传感器直接测量，但质心侧偏角的获取需要通过估计等间接方法得到。因此，准确估计车辆的质心侧偏角对于ESP的性能有重要意义。

从估计方法来讲，应用车辆动力学和运动学模型来建立车辆状态估计器或观测器的方法被广泛使用[2,3,4]。其中，基于车辆动力学模型来估计质心侧偏角需要准确计算轮胎的侧偏力。然而，在不同附着系数的路面上，同一轮胎侧偏角对应的侧偏力却存在差异，尤其当侧偏角较大，侧偏力进入到侧偏角—侧偏力特性曲线的非线性区域时，这种差异更加明显。因此，寻找适应于不同附着系数的路面的轮胎模型是正确估计车辆质心侧偏角的关键。

从估计算法来看，自回归贝叶斯滤波算法(包括扩展卡尔曼滤波、无迹卡尔曼滤波和粒子滤波)，模糊逻辑[5]、滑膜变结构[6]、多模型交互[7]等算法均得到了应用，尤其是卡尔曼滤波算法应用得最为广泛，但卡尔曼滤波算法的收敛速度和去噪性能与卡尔曼增益和卡尔曼协方差矩阵初值的选择有很大关系[8]。而状态空间形式的自回归最小二乘算法(RLS)因原理简单、使用方便，特别适合于在线估计，从而为其在车辆动力学变量估计中的应用奠定了基础。

基于上述考虑，引入含变参数的轮胎侧偏力模型，在非线性单轨车辆模型的基础上，应用状态空间表达的广义自回归最小二乘算法建立了质心侧偏角估计器，通过高附着系数路面上的车辆蛇形试验和变附着系数路面上的双移线试验验证了估计方法的性能。

1 车辆动力学模型

1.1 车辆模型

不同于车辆性能仿真时需要建立复杂的多自由度模型，在面向状态变量估计问题时，需要协调计算效率即实时性和估计精度两方面的需求。参考文献[9]，将车辆模型简化为单轨模型，具有横摆运动和侧向平动两个自由度，如图1所示。其中，a、b分别为车辆质心到前后轴的距离，αf、αr和Fy1、Fy2分别为前后轮的侧偏角和侧偏力，δ为前轮转角，β为车辆质心侧偏角，vx、vy分别为车辆质心纵向、侧向速度。

车辆侧向平动的动力学方程为

式中，m为车辆质量。

车辆横摆运动的动力学方程为

式中，为车辆横摆角速度;Iz为垂直转动惯量。

车辆质心侧偏角定义为车辆侧向速度和纵向速度间的夹角，即

将式(1)和式(2)代入式(3)求导并忽略次要因素可得[9]

设状态矢量X=(vy，ωr，β)T，从而有状态方程：

或

1.2 非线性轮胎侧偏力模型

车辆侧向运动受轮胎侧偏力的支配。基于车辆动力学模型进行质心侧偏角估计，需要事先选择轮胎模型来计算侧偏力。然而，车辆行驶的道路附着系数却是时变的。应用一种适合于特定附着系数路面(如高附着系数路面)的轮胎侧偏力模型无法准确计算其他附着系数路面上(如低附着系数路面)的轮胎力，即轮胎侧偏力在不同的附着系数路面是有差异的，尤其当轮胎的侧偏角较大时，这种差异更加明显，如图2所示。

“魔术公式”(magic formula)是基于试验数据的轮胎模型，它表达简洁且具有较高的拟合精度。文献[10]给出的轮胎侧偏力公式为

其中，α为轮胎侧偏角，Fy为轮胎侧偏力。在不同附着系数的路面上，参数A1、A2、A3、A4分别对应不同的拟合值。考虑到不同附着系数路面上计算轮胎侧偏力的需要，对式(7)进行简化并引入动态参数，将其看作随路面附着系数作相应变化的时变参数，使原公式变为

式中，C(t)为随时间t变化的参数;Fz为轮胎垂直载荷;D、E分别为模型的常系数。

可以看到，式(8)中的侧偏力仍表示为垂直载荷的函数，并且减少了模型参数。

图3所示为Fz=3kN和Fz=5kN两种垂直载荷下高附着系数路面上魔术公式(式(7))和变换公式(式(8))侧偏力的对比，参考文献[10]，侧偏力曲线形状因子D取为1.3,E通过拟合取为13。随着参数C的不断变化(分别取0.92、0.94和0.96)，式(8)计算的侧偏力能够接近魔术公式所得曲线。在中、低附着系数路面，随着参数C的变化也可以得到相同的曲线族。

由式(8)可分别得到轮胎前后轮侧偏力(下标1、2分别表示前轮、后轮)：

2 状态空间自回归最小二乘算法

自回归最小二乘算法本质上是一种递推形式的算法，其基本原理是应用新测量的信息来不断修正先前估计的变量，因此，该算法适用于变量的在线估计。本文选择以状态空间形式表达的广义自回归最小二乘算法[11]，其基本的表达式为

式中，为第k步预测值的更新;为从第k-1步到第k步的预测;Kk为第k步的估计增益;yk为第k步的观测值;Hk为第k步的量测方程。

考虑到车辆侧向加速度(而非侧向速度)可直接测量，故将加速度作为状态变量更加方便，即状态矢量为X=(ay，ωr，β)T，变形的魔术公式中的参数C也作为估计参数，从而状态矢量扩展为X=(ay，ωr，β，C1,C2)T。

式(11)中，^Xk|k-1是通过状态方程得到的状态变量：

式中，Uk为第k步的输入量;Vk为第k步的状态噪声。

状态方程和观测方程中的输入变量为车辆的前轮转角δ和纵向车速vx。δ通过方向盘转角传感器测量并由转向系传动比折算成前轮转角，即

V为系统的过程噪声。对于状态方程的展开，编程中采用欧拉法来实现其线性化，从而式(12)表示为

其中，k、k-1分别代表时间序列上的第k、k-1步;ΔT为时间步长;Fy1、Fy2为中间变量，可通过式(9)和式(10)求得。

Kk为估计增益，其表达式为

其中，λk为遗忘因子。引入遗忘因子的作用是增大新测量数据即车辆横摆角速度和侧向加速度的权重，从而减小旧测量数据的权重。J为状态方程对扩展状态变量的雅可比矩阵：

式中，x1、x2、x3分别为车辆侧向加速度ay、横摆角速度ωr和质心侧偏角β;f1、f2、f3分别为三个变量的状态方程。

实际程序计算中，对于雅可比矩阵J的数值解，采用差分方法来实现，即

式中，i、j分别为状态方程的第i个分量和状态变量的第j个分量;ε为摄动常数，程序计算中取ε=10-5。

本文对车辆质心侧偏角的估计需要测量车辆质心的侧向加速度和横摆角速度，这两个变量也是ESP系统配置的传感器可直接得到的，即测量变量

其分量形式为

其中，,N=(n1,n2)T为系统的观测噪声。

综合上述各式，可以实现状态变量的递推更新和估计，即单轨车辆模型中扩展状态变量的更新是通过新得测量值y(车辆横摆角速度和侧向加速度)与测量估计值之差来修正状态变量的估计值，从而使状态变量蕴含新测量的信息。

3 试验验证和分析

为了验证状态空间RLS对质心侧偏角非线性估计的效果，分别通过高附着系数路面的实车场地蛇形试验和附着系数变化的虚拟仿真试验进行验证。虚拟试验在Carsim软件中实现，车辆模型参数见表1。用车辆速度、方向盘转角、横摆角速度及车辆侧向加速度作为真实传感器的测量值，用得到的车辆质心侧偏角值对估计结果进行验证。编程中的遗忘因子λ取0.98，状态矢量的初值X=(0,0,0,0.8,0.8)T。

3.1 试验1（高附着系数路面上的蛇形试验）

按照国标GB/T6323.1-94，在干燥路面上进行蛇形试验。车辆以45km/h的初始速度进入测试区进行蛇形试验。测得方向盘转角、车辆速度、横摆角速度如图4～图6所示。由图7所示的车辆侧向加速度变化范围-5.89～5.52m/s2可知，轮胎进入到了侧偏角一侧偏力特性曲线的非线性区。质心侧偏角估计结果如图8所示，由图8可知，在整体趋势上估计结果和试验曲线一致并较好地跟踪了测量值，达到较高的估计精度。轮胎模型中的参数C1和C2在估计过程中的变化如图9～图10所示，均发生了相应的变化。

3.2 试验2（变附着系数路面的移线试验）

变附着系数路面能全面检验估计的适应性。由于条件限制，变附着系数路面试验在车辆虚拟仿真平台Carsim中完成。路面附着系数依次为0.5、0.9和0.65，具体如图11所示，车辆在三段路面上均进行了不同程度的双移线试验。车辆的初始速度为80km/h，整个试验持续约24s。方向盘转角、车辆速度和横摆角速度分别如图12～图14所示。由图15可知，车辆在三种附着系数路面的最大侧向加速度分别为-5m/s2、-5.9m/s2和-5.6m/s2，轮胎均进入到了侧偏角—侧偏力特性曲线的非线性区。

质心侧偏角的估计结果如图16所示，可以看到，当车辆在低附着系数路面上运动时，其侧向运动较为剧烈，产生了较大的质心侧偏角，最大值达到0.148rad。整体来看，车辆质心侧偏角估计值较好地跟随了测量值。误差相对较大的地方出现在图形的“波峰”和“波谷”处。在这些地方，由于侧向加速度达到极值，轮胎侧偏相对严重，轮胎模型的参数C1和C2需要经历调整的过程，其相应变化过程分别如图17、图18所示，可以看到，它们的变化趋势随着路面附着系数的变化发生了相应变化，整体看来也经历了低-高-低的变化过程，与道路附着系数变化趋势相似。

作为对比，本文还应用扩展卡尔曼滤波(extended Kalman filter,EKF)算法对车辆质心侧偏角进行了估计。从图16中可以看到EKF的估计值同RLS估计值一样，较好地跟踪了测量值。定量的比较如表2所示，RLS和EKF估计值的最大误差分别为1.579°和1.889°，前者优于后者;估计均方根误差分别为0.1364°和0.0958°，后者估计精度稍微优于前者，但都在可以接受的范围内。从单位计算用时来看，在同样的计算机硬件条件下两者性能近乎相似，分别为0.893ms和0.917ms。轮胎模型参数C1和C2的变化见图17和图18。对比可知，基于EKF的参数变化和基于RLS的参数变化趋势一致，波动幅度小于RLS的估计值。但从使用条件来看，应用扩展卡尔曼滤波时要求事先假设噪声的统计服从高斯分布，而RLS却无此要求。

4 结论

(1)基于状态空间形式的自回归最小二乘算法建立了车辆质心侧偏角非线性估计器，对车辆动力学状态和轮胎参数同时进行了估计。

(2)通过实车高附着路面的蛇形试验和附着系数变化路面的双移线试验，对状态空间RLS算法的估计性能进行了验证。结果表明，即使车辆发生大侧偏现象使轮胎进入到侧偏角—侧偏力特性曲线的非线性区域，自回归最小二乘算法依然能够实现对单一附着路面和变附着系数路面上车辆质心侧偏角的估计，呈现出良好的适应性。

(3)状态空间RLS和EKF估计结果的对比表明，两者最大误差分别为1.579°和1.889°，均方根误差分别为0.1364°和0.0958°，单位计算耗时分别为1.193ms和1.217ms。从估计精度和计算效率来看，两者的估计性能相似，但EKF在使用过程中要求系统的噪声统计规律服从高斯分布，因此RLS具有更广的使用范围。

参考文献

[1]van Zanten A Z.Bosch ESP Systems:5 Years of Experience[C].SAE Paper,2000-01-1633.

[2]郭孔辉,付皓,丁海涛.基于扩展卡尔曼滤波的汽车质心侧偏角估计[J].汽车技术,2009,4(1):1-4.Guo Konghui,Fu Hao,Ding Haitao.Estimation of Vehicle Sideslip Angle Based on Extended Kalman Filter[J].Automobile Technology,2009,4(1):1-4.

[3]Fukada Y.Slip-angle Estimation for Vehicle Stability Control[J].Vehicle System Dynamics,1999,32(4):375-388.

[4]Piyabongkan D,Rajamani R,Grogg J,et al.Development and Experimental Evaluation of a Slip Angle Estimation for Vehicle Stability Control[J].IEEE Transactions on Control Systems Technology,2009,17(1):78-88.

[5]Cheli F,Sabbioni E,Pesce M,et al.A Methodology for Vehicle Sideslip Angle Identification:Comparison with Experimental Data[J].Vehicle System Dynamics,2007,45(6):549-563.

[6]Stephant J,Charara A,Meizel D.Evaluation of a Sliding Mode Observer for Vehicle Sideslip Angle[J].Control Engineering Practice,2007,15(7):803-812.

[7]Tsunashima H,Murakami M,Miyata J.Vehicle and Road State Estimation Using Interacting Multiple Model Approach[J].Vehicle System Dynamics,2006,44(S):750-758.

[8]Haykin S.Kalman Filtering and Neural Networks[M].New York:John Wiley&Sons,2001.

[9]Jazar R N.Vehicle Dynamics:Theory and Applications[M].New York:Springer,2008.

[10]Bakker E,Nyborg L,Pacejka H B.Tyre Modelling for Use in Vehicle Dynamics Studies[J].SAE Paper,870421,1987.

非线性最小二乘篇6

变电站接地网对变电站的安全运行具有重要的意义,是电力运行部门十分关注的问题,近年来也成为学术界研究的一个热点[1,2,3,4,5,6]。接地网导体埋于地下,由于化学腐蚀等原因使接地导体产生锈蚀极大地改变接地网的参数,从而影响接地网的接地性能。大型变电站接地网覆盖面积大,接地导体网格的数量也较大,因此接地导体的腐蚀程度是不均匀的,这给运行部门的维护检修带来相当大的困难。为了节约地网改造的成本,需要对接地导体网格的腐蚀情况有准确的评价,为接地网的改造提供准确的信息。

目前对接地网腐蚀程度的评价主要是进行测试,并通过一定的模型和算法诊断各个接地网格导体的实际电阻,与接地网的设计值对比即可知道接地导体的腐蚀程度,并在此基础上制定接地网改造方案。其中文献[3]利用接地网的可及点电压值检测接地网故障的判据,应用网络撕裂法概念建立接地网区域故障的诊断方程。文献[4]提出了一种基于禁忌搜索算法的接地网故障诊断方法,采用轮换激励位置和每处激励多处测量的方法,使可及节点得到更充分利用,观测信息显著增加。

本文提出了一种基于灵敏度分析的线性二乘优化模型。根据电路理论,电路元件参数的变化会引起电路响应参数的变化,当元件参数的变化量很小时,可以把电路响应参数的变化量近似表示成元件参数的变化量及响应参数对各个元件参数灵敏度的线性关系,从而建立基于最小线性二乘优化的诊断方程。通过求解线性最小二乘优化问题就可以诊断出腐蚀接地网各个接地导体的实际电阻值。

2 故障诊断模型

设接地网的支路电阻为Ri(i=1,2,…,b,其中b为支路数),测试参数(或响应参数)为各个可及测试节点的电压Vx(x=1,2,…,N,其中N为接地网可及节点数),若在接电网的两个可及测试节点i,i′间施加激励(可以为电流激励也可为电压激励),但考虑节点i,i′到电源间的引线电阻可能造成的影响,选择电流源作为接地网的激励源。接地网的测试框图如图1所示。

在如图1所示的测试条件下,故障条件接地网测试节点x的节点电压为V′x,非故障条件(所有接地支路电阻值都为标称值)节点x的节点电压为Vx,故障条件及非故障条件下Vx的变化量为ΔVx=V′x-Vx,Vx对各个支路电阻的灵敏度为Sundefined=∂Vx/∂Ri(i=1,2,…,b)。故障条件下每个支路电阻值为R′i,标称值为Ri,故障条件及非故障条件下各支路电阻的变化量为ΔRi=R′i-Ri(i=1,2,…,b),根据以上介绍的灵敏度理论,ΔVx可以表示如下:

undefined

这样即可把电路可及节点的节点电压变化量ΔVx表示成关于元件参数变化量ΔRi=R′i-Ri(i=1,2,…,b)及灵敏度Sundefined=∂Vx/∂Ri(i=1,2,…,b)的线性关系。

在同样的激励条件下(激励位置和大小相同),可以选择α不同的测试节点,得到一组多个测试节点电压的变化量与元件参数变化量及灵敏度的线性关系:

undefined

一般接地网的支路电阻数较多,为了通过解式(2)对各个支路电阻值进行识别,需要建立较多的形如式(1)的线性等式。若在相同的激励条件下选择过多的测试节点,方程(2)的条件会变坏。因此为了改善方程组(2)的条件,应该选择不同的节点i,i′施加激励,而且在同一激励条件下测试不同节点的节点电压。

3 灵敏度算法

接地网为纯电阻网络,根据电路理论中的节电电压法,电路的节电电压Vn向量与电路导纳矩阵Yn及电路施加的电流激励Is向量之间的关系可以表示为:

undefined

令:B=Yundefined。

undefined

节点电压对支路导纳的灵敏度为:

undefined

再由Yn·Yundefined=1得:

undefined

根据灵敏度关系,节点k的节点电压Vnk对电阻Rij的灵敏度为:

undefined

4 线性最小二乘优化模型

方程组(2)可以表示成如下的矩阵形式:

undefined

其中:ΔR=[ΔR1,ΔR2,…,ΔRb]T,d=[ΔVx1,ΔVx2,…,ΔVxα]T,C为式(2)中的灵敏度系数矩阵。

线性最小二乘优化模型为:

undefined

式(9)的线性最小二乘优化模型可以使用Matlab中的工具箱直接求解。

5 实例仿真

本文提出的的模型和方法在如图2所示的实际接地网上进行了数值仿真,数值仿真的结果如表1所示。

从表1可以看出实际电阻的诊断值与设定的故障值相当接近,误差不超过7.5%,取得了较好的诊断结果。

6 结语

本文在基本电路理论和灵敏度分析的基础上,提出了一种线性最小二乘优化模型。该方法充分利用可及节点所提供的信息,建立的模型可以反映节点电压变化量与接地网导体电阻变化量之间的关系,在此基础上建立了最小二乘优化模型。从仿真实例的结果可以看出该方法具有较高的故障诊断精度。

参考文献

[1]刘渝根,滕永禧,陈先禄,等.接地网腐蚀的诊断方法研究[J].高电压技术,2004,30(6):19-21.

[2]王新翠,彭敏放.基于对称复镜像法接地网接地电阻的算法[J].现代电子技术,2007,30(9):143-145.

[3]江修波.接地网故障诊断的一种新方法[J].福州大学学报:自然科学版,2005,33(6):749-752.

[4]程红丽,刘健,王森,等.基于禁忌搜索的接地网故障诊断[J].高电压技术,2007,33(5):139-142.

[5]刘渝根,吴立香,王硕.大中型接地网腐蚀优化诊断实用化分析[J].重庆大学学报:自然科学版,2008,31(4):417-420.

非线性最小二乘篇7

在多元线性回归模型中,如果解释变量之间存在着密切的线性相关关系,就称它们之间存在着多重共线性. 在出现多重共线性情形时,普通最小二乘估计不再适用; 回归参数的估计值方差会很大,从而影响自变量对因变量的解释;估计的精度会降低; 估计的效果也会变坏. 在实际经济问题的多元回归分析中,多重共线性的现象很多,这时我们就应该寻找另外的回归方法对参数进行估计.

二、方法介绍

如果在实际问题中出现了多重共线性的现象,我们可以选择用有偏回归方法———岭回归( RR) 和偏最小二乘回归( PLS) 来处理. 岭回归是利用岭估计( X'X + k I)- 1X' Y来替代普通最小二乘估计( X'X)- 1X' Y,从而消除了普通最小二乘估计中矩阵X'X无法求逆的问题. 偏最小二乘回归是先在自变量集和因变量集中分别提取第一潜在因子t1与u1,其中t1与u1分别是自变量与因变量的线性组合,要求t1与u1尽可能多地提取所在变量组的变异信息,且t1与u1的相关程度达最大,然后建立因变量与t1的回归方程,若回归方程不能达到满意的精度,则继续提取第二潜在因子,否则停止.

三、实例比较

根据理论及对现实情况的认识,拟建立以我国国民总收入( 单位: 亿元) 为因变量y,以就业人员数( 单位: 万人) 、财政收入( 单位: 亿元) 、能源生产总量 ( 单位: 万吨标准煤) 、国有单位工资总额( 单位: 亿元) 和城镇集体工资总额( 单位: 亿元) 分别为自变量x1,x2,x3,x4,x5的线性回归模型. 由《中国统计年鉴》查得相关数据如下:

在SAS软件上使用REG过程来建立最小二乘回归方程,所有自变量的方差膨胀因子都大于100,诊断出模型中存在非常严重的多重共线性问题. 用最小二乘法所得到的回归方程为y = - 431189 + 6. 13224x1- 0. 18088x2+0. 44051x3+ 5. 69125x4- 13. 63786x5.

可以看到方程中,自变量x2,x5的系数为负,这显然与事实不符,这正是由多重共线性所导致,因此最小二乘回归求出的回归方程不利于模型的解释,下面改用岭回归方法来建模.

用SAS软件中的REG过程,求解岭回归方程. 由岭迹图可以看出,当岭参数k≥0. 02后,岭迹曲线趋于稳定,因此,取k = 0. 02的岭回归估计来建立岭回归方程为

这时,回归系数的符号符合实际意义.

现在用偏最小二乘回归方法来进行处理,用SAS软件中的PLS过程建立偏最小二乘回归方程,用最常用的舍一交叉验证法来抽取偏最小二乘的成分,结果抽取了3个偏最小二乘成分,得到偏最小二乘回归方程为

这时,回归方程中的回归系数的符号也都符合实际意义.

根据前面得出的岭回归方程和偏最小二乘回归方程,计算出衡量模型拟合效果好坏的平均绝对百分误差和复测定系数,得到相应的数值如下:

四、总结

从上例可以看出,在多元线性回归模型中出现共线性问题时,最小二乘回归方法已经不再适用,而用岭回归和偏最小二乘回归这两种有偏回归方法都可以处理多重共线性问题,且从表2的结果可知,两种方法建立的回归方程拟合的效果都不错,而偏最小二乘回归方法相对岭回归方法要更优.

摘要：文章介绍了处理多元线性回归模型中多重共线性问题的有偏回归方法——岭回归和偏最小二乘回归,并通过实例比较了两种方法建立的回归方程的拟合效果,而偏最小二乘回归方法相对岭回归方法要更优.

非线性最小二乘篇8

在实验数据的处理、分析时, 数据拟合是经常采用的一种方法。数据拟合的目标是找到能反映变量之间关系的一种表达式, 使其在某种准则下最佳地接近已知数据。其原理有最小二乘法、契比雪夫法等, 且以最小二乘法最为常见和重要。

随着人类认识能力的不断进步以及计算技术的快速发展, 对于变量之间的未知关系, 应用曲线拟合的方法揭示其内在规律具有重要的理论和现实意义。针对最小二乘曲线拟合的有关问题以及相应的MATLAB实现进行探讨。

1 最小二乘曲线拟合

曲线拟合是指:已知n个数据点 (xi, yi) , i=1, 2, ……, n, 其中xi不全相同, 寻求函数f (x;a1, a2, …, am) 的待定参数a1, a2, …, am的一组取值, 使得在这组取值之下, 函数f (x;a1, a2, …, am) 与已知n个数据点整体上最为接近。

最小二乘曲线拟合方法根据已知数据, 首先构造出能够反映含有待定参数的函数f (x;a1, a2, …, am) 与n个数据点 (xi, yi) , i=1, 2……, n偏离程度的函数:

然后应用数学方法求函数J (a1, a2, …, am) 的最小值a1, am2, i…n, am=J (a1, a2, …, am) , 此时a1, a2, …, am的取值就是所求的待定值。这样一组取值使得函数f (x;a1, a2, …, am) 与n个数据点在二次平方和意义下最为接近。

2 最小二乘曲线拟合的MATLAB实现

2.1 线性最小二乘曲线拟合的MATLAB实现

线性最小二乘曲线拟合的含有待定参数的函数形式为:

其中rk (x) 为事先选定的一组已知函数, ak为待定系数, k=1, 2, …, n, m

该方程组可简化为RTRA=RTY, 其中:

当r1 (x) , …, rm (x) 线性无关时, RTR可逆, 方程组RTRA=RTY有唯一解。

根据线性最小二乘拟合理论可得关于待求参数构成的矩阵A的解, 利用MATLAB下命令A=RY, 可以直接求得待求参数, 下面举例说明。

例1下表给出了一组实测数据

已知数据 (xi, yi) 的函数原型为:

试用已知数据进行曲线拟合, 求出待定参数的值。

在MATLAB中输入以下程序:

运行程序, 输出待定参数结果为:

下面画出拟合得到的曲线及已知的数据散点图:

由图1可见, 曲线拟合效果非常好。事实上, 所给数据就是由已知曲线:

2.2 非线性最小二乘曲线拟合

在最小二乘曲线拟合中, 如果待定参数不能全部以线性形式出现, 如指数拟合函数f (x) =a0+a1e-a2·x, 这便是非线性最小二乘曲线拟合。

非线性最小二乘曲线拟合与线性最小二乘曲线拟合的原理没有什么区别, 但是最小二乘的解常常难以通过人的手算实现, 从而制约了该方法的应用。随着计算机技术的进步、专业软件的不断涌现, 这一问题的求解已不再困难。但是, 非线性曲线拟合中初值的选取是一个重要的问题, 目前为止还没有固定的理论或方法给出一般性的结论。

在MATLAB中实现非线性最小二乘曲线拟合有三个常用的命令。

(1) lsqcurvefit () 命令, 其使用格式为x=lsqcurvefit (fun, x0, xdata, ydata) , 其中fun是要拟合的非线性函数, x0是初始参数, xdata, ydata是拟合点的数据, 该函数最终返回系数矩阵。

(2) nlinfit () 命令, 其应用格式为beta=nlinfit (x, y, fun, beta0) , 其中x和y是拟合点数据, fun是被回归 (拟合) 的函数, beat0是初始参数。

(3) lsqnonlin () 命令, 其应用格式为x=lsqnonlin (fun, x0) , 其中fun的定义与lsqnonlin () 函数中fun的定义有差别, x0仍为初始参数向量, 将输出的系数结果放在变量x中。

例2假设已知

并已知该函数原型为y (x) =a1exp (a2x) +0.54exp (a3x) .sin (a4x) , 其中ai为待定系数。

在MATLAB中输入以下命令:

%建立函数原型, 则可以根据他来进行下面的求取系数的计算

运行结果为:

所求得的系数与原式中的系数相近。

如果想进一步提高精度, 则需修改最优化的选项, 函数的调用格式也将随之改变。

在MATLAB中输入以下命令:

%修改精度, 即是修改其限制条件

%两个空矩阵表示系数向量的上下限a', res运行结果为:

下面画出由拟合得到的曲线及已知的数据散点图, 如图2。

3 结论

首先介绍了最小二乘法的原理, 其次对线性最小二乘曲线拟合和非线性最小二乘曲线拟合应用MATLAB进行了具体实现, 使得相关曲线拟合理论更加生动易懂, 这使我们可以在今后的研究和工作中建立曲线模型对相应对象进行曲线拟合, 并应用MATLAB来实现, 从而找到更好更形象的反映变量之间关系的曲线, 也就找到最合适的模型了。

参考文献

[1]薛定宇, 陈阳泉.高等应用数学问题的MATLAB求解[M].北京:清华大学出版社, 2004.

[2]刘琼荪, 何中市, 傅鹂等.数学实验[M].北京:高等教育出版社, 2000.

[3]梁国业, 廖健平.数学建模[M].北京:冶金工业出版社, 2004.

非线性最小二乘篇9

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面。其中一种简单的情形是:根据相关影长数据确定参照物位置, 这一技术也称为太阳影子定位技术。在本文中, 根据2015 年全国大学生数学建模竞赛A题附件1 的数据[1], 确定出该地点的位置。

2 影长与时间的函数关系

2.1 影响影长的各个参数

太阳高度角是指太阳光的入射方向和地平面之间的夹角, 可以用太阳高度角来估计影子场电影。太阳高度角h与实测地的地理纬度θ, 太阳直射角纬度 φ, 时角 Ω 满足以下函数关系式[2]:

其中, 关于各个位置参数确定如下:

太阳直射点每时都在向西移动, 每小时移过15 度经度, 在地理题的计算中可粗略取每天移动0.25 度纬度[3]。据此求解出太阳直射角纬度与日期之间的关系, 进而通过MATLAB编程, 得到了每天太阳直射点的纬度, 部分表见表1, 在4 月18 日的太阳直射点纬度为 φ=7。

一个天体的时角表示该天体是否通过了当地的子午圈, 其数值则表示了该天体与当地子午圈的角距离, 以小时来计量 (1 时角=15 度) 。根据以上理论, 建立时角的关系表达式[4], 即:

其中Cr表示北京时;LC表示经度校正 (4min/度) , 经度大于120 的经度校正为正, 小于120 的经度校正为负;EQ为时差。

2.2 建立影长随时间变化函数

根据式 (1) , 我们得到了太阳高度角与时角、太阳直射点纬度以及实测地的地理纬度的函数关系式, 进而我们可以通过太阳高度角与影子长度之间的关系, 建立影长随时间变化的函数关系式, 其步骤如下:

根据太阳高度角的定义, 得到太阳高度角与影子长度之间满足以下关系:

其中L为影子长度;H为旗杆高度3m;h为太阳高度角;

根据式 (1) , 将各个参数值代入公式进行整理, 得到太阳高度角与各个参数之间的函数关系:

再根据式 (4) , 将 (5) 转化为:

其中:

3 利用最小二乘拟合确定直杆高度以及所处位置

3.1 模型的建立

已知对应时刻的影子长度 (ti, Li) , 影长与时间的函数关系, 见公式 (6) , 利用最小二乘原理, 建立如下模型:

其中c1, c2, c3的初始值设为1, 0.1, 10

3.2 参数及经度的确定

根据参数的初始值, 将21 组 (ti, Li) 代入拟合方程中进行拟合。在拟合之前, 本文根据c1, c2, c3的表达式做了以下猜想:发现c1, c2关于其自身的自变量的表达式十分简单, 而c3表达式较为复杂, 并且c3在拟合方程中的表达式也较为复杂, 所以关于c3的拟合效果相对于c1, c2来说, 拟合效果较差, 可能得不到精确的解。

而在实际拟合过程也证明了我们的猜想。在实际拟合过程中, 随着初始值的不断变化, c1, c2变化不大或者只在固定的值的范围内小幅度变化, 而c3的变化幅度较大, 且不能找到唯一确定的值, 因此通过求出的纬度以及直杆高度来确定经度。经求解:c1=2.1316m, c2=0.313或0.1287, 即H=2.1316m, φ=18.24N或7.3969N。

将H, φ 中的数据代入 (5) 中, 用MATLAB编程计算, 我们得到21个时间点所对应的经度位置, 我们发现, 这21 个经度变化都十分小, 故我们经度该位置的经度就选取求出的21 个经度的平均值, 得到4 个经度, 进而得到四个固定直杆可能所处的地理位置: (18.24N, 105.9528E) , (18.24N, 37.4972E) , (7.3969N, 107.1988E) , (7.3969N, 36.2512E) 。

4 地理位置的检验

根据测定的关于x, y的坐标, 我们发现随着时间的增加, 影子长度也在不断增加, 说明当地时间已经过了正午时间12:00, 而此时的北京时间为14:42-15:42, 结合北京时间以东八区120 度为标准, 每向西15 度, 当地时间往前1 小时的规律, 我们计算出在79.5 E时, 北京时间为14:52, 当地时间恰好为12:00, 因此我们所确定的经度必须大于79.5E。进而确定可能的地理位置为: (18.24N, 105.9528E) 和 (7.3969N, 107.1988 E) 。

根据地图确定地理位置为 (18.24N, 105.9528E) 的地点为越南, 地理位置为 (7.3969N, 107.1988E) 的地点为泰国湾 (为靠近马来西亚的一个海峡) , 故这个位置不符合我们的标准, 故我们确定的直杆所处位置大致在越南。

5 模型的检验

我们将所确定的地理位置代入关于太阳高度角的式子 (5) 中, 画出14:42-15:42的影子长度的变化曲线l1, 同时拟合出附件1中影子长度随时间变化曲线l2, 两者之间的误差, 可以这21个时刻的误差分别为:0.0004, 0.0005, 0.0005, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0005, 0.0005, 0.0005, 0.0005, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004。两条影子长度随时间变化的曲线 (见图1) 。

根据图1, 结合21 个点的误差, 可以发现两条曲线误差最大的为0.0006, 误差非常小, 认为该模型合理。

参考文献

[1]http://www.shumo.com/2015cumcm.html[OL].

[2]http://baike.baidu.com/view/86609.htm[OL].

[3]http://baike.baidu.com/link?url=_a Uu Zt IGSo3DBZgj Wi SGC2Nq H9NNBj Ec XX37JYu5MWs F6M22qio Ndju Y81e Nje P31Llpv JKTBhf_e PWF9g5CXa#1[OL].

[4]王国安, 米鸿涛, 邓天宏, 李亚男, 李兰霞.太阳高度角和日出日落时刻太阳方位角一年变化范围的计算[J].气象与环境科学, 2007, 30:161-163.

【非线性最小二乘】推荐阅读：