最小二乘问题

2024-10-22

最小二乘问题（精选10篇）

最小二乘问题篇1

考虑下列二阶抛物问题

${\begin{cases} u_{t} (x, t) - \nabla \cdot (a (x) \nabla u (x, t)) = f (x, t), \\ (x, t) \in Ω \times J, \\ u (x, t) = 0, (x, t) \in Γ \times J, \\ u (x, 0) = u_{0} (x), x \in Ω ‚ \end{cases} (1)$

(1)式中J=[0,T],Ω⊂Rn为有界凸域,边界Γ是Lipschitz连续的。假定a(x)有界,即存在常数α*和α*,使得

0≤α*≤a(x)≤α*,∀x∈Ω (2)

对上述问题人们提出了许多数值模拟格式,如有限元格式[1],混合有限元格式[2],扩展混合有限元格式[3],最小二乘混合有限元格式[4]等。现将最小二乘思想和扩展混合有限元方法相结合,提出并分析了最小二乘扩展混合有限元格式。该格式继承了最小二乘和扩展混合元方法的优点,即格式能同时高精度逼近未知函数,未知函数的梯度以及流体的通量,所形成的有限元方程组是对称正定的,有限元空间无须满足LBB相容性条件。数值算例说明了该格式的有效性。

1 最小二乘扩展混合元格式

定义空间Λ=(L2(Ω))n,V=H $_{0}^{1}$ (Ω),W=H(div;Ω)={q∈(L2(Ω))n,∇·q∈L2(Ω)},相应范数为 $∥ μ ∥_{0, Ω} = (μ, μ)^{\frac{1}{2}}, ∥ v ∥_{1, Ω} = (∥ v ∥_{0, Ω}^{2} + ∥ \nabla v ∥_{0, Ω}^{2})^{\frac{1}{2}}, ∥ q ∥_{Η (d i v; Ω)} = (∥ q ∥_{0, Ω}^{2} + ∥ \nabla \cdot q ∥_{0, Ω}^{2})^{\frac{1}{2}}$ ,其中(·,·)表示空间(L2(Ω))n中的内积。作时间剖分τ=T/N,N为正的常数,在时间tn=nτ,n=1,2,…,N处逼近方程的解。 $Ω = \underset{Κ \in τ_{h}}{\cup} Κ$ 是Ω的一个剖分,h=max{diam(K):K∈Th}。设Vh,Λh和Wh分别是V,Λ和W的有限维子空间,具有通常逼近性质。Vh,Λh和Wh的一个标准的选法是k,r和s次的分片多项式空间。

记λ=-∇u,σ=aλ,并用一阶差分近似ut,则式(1)可变成离散形式I

${\begin{cases} \bar{\partial}_{t} u^{n} = \nabla \cdot σ^{n} = f^{n} + R_{1}^{n}, x \in Ω, \\ λ^{n} + \nabla u^{n} = 0, x \in Ω, \\ a λ^{n} - σ^{n} = 0, x \in Ω, \\ u^{n} (x) = 0, x \in Ω, \\ u^{0} (x) = u_{0} (x) x \in Ω \end{cases} (3)$

(3)式中 $\bar{\partial}_{t} u^{n} = (u^{n} - u^{n - 1}) / τ, R_{1}^{n} = \bar{\partial}_{t} u^{n} - u_{t}^{n}$ 。

对于(u,λ,σ)∈V×Λ×W和(v,u,q)∈V×Λ×W,定义双线性形式

$\begin{array}{l} a_{Ν} (u, λ, σ; υ, μ, q) = (u + τ \nabla \cdot σ, υ + τ \nabla \cdot q) + \\ τ (λ + \nabla u, μ + \nabla υ) + τ (a^{- 1} (a λ - σ), a μ - q) (4) \end{array}$

根据离散形式式(3),省略时间截断误差定义最小二乘扩展混合元格式

格式Ⅰ 给定初始值u $_{h}^{0}$ ∈Vh,对于n=1,2,…,N,求(u,λ,σ)∈Vh×Λh×Wh,使得

$\begin{array}{l} a_{Ν} (u_{n}^{n}, λ_{n}^{n}, σ_{n}^{n}; υ_{h}, μ_{h}, q_{h}) = (u_{h}^{n - 1} + τ f^{n}, υ_{h} + \\ τ \nabla \cdot q_{h}), \forall (υ_{h}, μ_{h}, q_{h}) \in V_{h} \times Λ_{h} \times W_{h} 。 (5) \end{array}$

2 解的存在唯一性

为证明格式Ⅰ解的唯一性,首先说明双线性形式aN(·;·)是强制的。

引理1 对任意的(v,μq)∈V×A×W存在c>0,使得aN(v,μ,q;v,μq)≥c(‖v‖ $_{1, n}^{2}$ +‖μ‖ $_{0, n}^{2}$ +‖q‖H(div n))。 (6)

证明引入适当选取的常数β>0,式(4)可整理为

$\begin{array}{l} a_{Ν} (υ, μ, q; υ, μ, q) = (τ \nabla \cdot q + (1 - β) υ, \\ τ \nabla \cdot q + (1 - β) υ) + τ (q - a (μ + β υ), \\ a^{- 1} (q - a (μ + β υ))) + \\ 2 τ (μ + (1 - a β) \nabla υ, μ + (1 - a β) \nabla υ) + \\ (2 β - β^{2}) (υ, υ) 。 (7) \end{array}$

由于υ(x)=0,x∈Γ,所以由Poincaré-Friedrichs不等式知,存在正常数CF,使得

‖v‖ $_{0, n}^{2}$ ≤CF‖∇v‖ $_{0, n}^{2}$ (8)

结合式(2)有

$\begin{array}{l} a_{Ν} (υ, μ, q; υ, μ, q) \geq β (2 (τ α_{*} + C_{F}) - \\ β (τ α^{* 2} + τ α^{*} + C_{F})) ∥ \nabla υ ∥_{0, Ω}^{2} 。 (9) \end{array}$

取β=τα*/(τα*2+τα*+CF),则

$a_{Ν} (υ, μ, q; υ, μ, q) \geq C ∥ \nabla υ ∥_{0, Ω}^{2} \geq C ∥ υ ∥_{1, Ω}^{2} 。 (10)$

从而由式(10)和式(4),知

$\begin{array}{l} ∥ μ ∥_{0, Ω}^{2} \leq 2 ∥ μ + \nabla υ ∥_{0, Ω}^{2} + 2 ∥ \nabla υ ∥_{0, Ω}^{2} \leq C a_{Ν} \\ (υ, μ, q; υ, μ, q), (11) \\ ∥ q ∥_{0, Ω}^{2} \leq 2 ∥ a μ - q ∥_{0, Ω}^{2} + 2 ∥ a μ ∥_{0, Ω}^{2} \leq C a_{Ν} \\ (υ, μ, q; υ, μ, q), (12) \\ ∥ \nabla \cdot q ∥_{0, Ω}^{2} \leq \frac{2}{τ} ∥ υ + τ \nabla \cdot q ∥_{0, Ω}^{2} + \frac{2}{τ} ∥ υ ∥_{0, Ω}^{2} \leq C a_{Ν} \\ (υ, μ, q; υ, μ, q), (13) \end{array}$

结合式(10)—式(13)便证得式(6)成立。

定理1 设f∈L2(Ω),则格式Ⅰ存在唯一解(μ $_{h}^{n}$ ,λ $_{h}^{n}$ ,σ $_{h}^{n}$ )∈Vh×Λh×Wh。

证明考虑空间V×Λ×W,其范数定义为‖υ‖1,Ω+‖μ‖0,Ω+‖Q‖H(div;Ω)。由引理1知双线民生形式αN(·;·)在空间Vh×Λh×Wh中是强制的;而且由Cauchy不等式易证αN(·;·)在Vh×Λh×Wh中是有界的。因此由Lax-Milgram引理知格式Ⅰ解存在唯一。

3 数值算列

为验证所提格式的有效性,在区域[0,1]×[0,1]上考虑方程

${\begin{cases} u_{t} - (a u_{x})_{x} = f, (x, t) \in (0, 1) \times (0, 1], \\ u (0, t) = u (1, t) = 0, t \in [0, 1], \\ u (x, 0) = \sin π x, x \in [0, 1] \end{cases} (14)$

(14)式中 $a = \frac{1}{100}, f = (1 + \frac{π^{2}}{100}) e^{t} \sin π x$ 。

容易验证方程的真解为 $u = e^{t} \sin π x, λ = - π e^{t} \cos π x, σ = - \frac{1}{100} π e^{t} \cos π x$ 。取Λh为分片常数多项式空间,Vh和Wh为分片线性多项式空间,用MATLAB编程计算近似解uh,λh和σh,绘出真解和最小二乘扩展混合元解的比较图,见图1—图3。

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

关键词：二阶抛物问题,最小二乘,扩展混合有限元方法

参考文献

[1]Thome V.Galerkin finite element methods for parabolic problems.Lecture Notes in Mathematics,Springer,1984

[2]Johnson C V.Thome,Error estimates for some mixed finite element methods for parabolic type problem.RAIRO Numer Anal,1981;15:41—78

[3] Guo Ling,Chen Huanzhen.An expanded characteristic mixed finite element method for a convection-dominated transport problem.Journal of Coputational Mthematics,2005;23:479—490

[4] Yang Danping.Analysis of least-squares mixed finite element methods for nonlinear nonstationary convection-diffusion problems.Math Comp 1999;69:929—963

最小二乘问题篇2

基于模糊判断矩阵的对数最小二乘排序方法

基于决策者提供方案偏好信息的排序方法研究是决策分析的一个重要问题.不同类型的偏好信息有不同的排序方法和理论依据.根据模糊判断矩阵完全一致性的.概念,从最优化角度提出了模糊判断矩阵的对数最小二乘排序方法,并研究了其性质,给出了相应的一致性检验指标.最后进行了算例分析.

作者：和媛媛周德群王强 HE Yuan-yuan ZHOU De-qun WANG Qiang 作者单位：南京航空航天大学经济与管理学院,江苏,南京,210016刊名：中国管理科学 ISTIC PKU CSSCI英文刊名：CHINESE JOURNAL OF MANAGEMENT SCIENCE年，卷(期)：15(z1)分类号：C934关键词：模糊判断矩阵排序方法完全一致性对数最小二乘法

最小二乘问题篇3

L 1+ae-bx

中参数L，a，b的估计往往是基于一定的经济意义或生物统计的背景，但从纯粹的数学方程角度来看，当L，a，b的取值范围不加任何限制时，是否也能估计出L，a，b呢？文章证明了当（L，a，b）

∈R3

时，y=

L 1+ae-bx

中参数L，a，b不能用最小二乘法估计，并给出了其参数在一定的限制条件下，可以用最小二乘法估计.

[关键词] 逻辑曲线参数最小二乘估计

[中图分类号] G633.6 [文献标识码] A [文章编号] 1674 6058（2016）17 0049

一、引言

基于残差平方和最小的思想，最小二乘法已在社会生产和实践中得到了广泛的应用.在解决非线性回归模型的参数估计问题时，我们往往将其转化为线性回归模型，进而用最小二乘法估计出其参数，但有的模型并不能通过直观的观察和简单的计算就能转化为线性回归模型，也就谈不上用最小二乘法估计其参数了.这就引发了我们对非线性回归模型的参数是否能用最小二乘法估计的沉思，如常见的逻辑曲线y=

L 1+ae-bx

（（L，a，b∈R3）），其参数L，a，b能否用最小二乘法的思路估计出或至少能被确定属于一定的范围呢？这就是本文所要研究的重点.

二、猜想

任意取n个不同的样本（xi，yi），i=1，2，…，n.只要样本中不含用（0，y（0））；那么逻辑曲线y= L 1+ae-bx （（L，a，b）∈R3）的参数L，a，b就不能用最小二乘法估计出且不能估出（L，a，b）属于R3中某一真子空间，即无法确定（L，a，b）∈H，HR3.

三、证明

假设用最小二乘法可以估计出（L，a，b），尽管（L，a，b）并不一定唯一，即（L，a，b）也可属于R3中一真子空间.

第一步，观察y= L 1+ae-bx ，

∵（xi，yi）可取n个不同的样本，

∴可推断出L的估计，L≠0.

第二步，（Ⅰ）y= L 1+ae-bx （*）两边同时对x求导，得

dy dx =by（1-

y L

）

，其中，y（0）=

L 1+a

（**），

∴参数（L，a，b）的估计若满足（*）式，则也满足（**）.

（Ⅱ）考查微分方程y= L 1+ae-bx ，y（0）= L 1+a （**）.

令f（x，y）=by（

1- y L ）

，取R3平面上一区域D：|x| ≤h1|y|≤h2.

不难发现：f（x，y）在D上连续，且

|f（x1，y1）-f（x2，y2）|= 1 L |bL

（y1-y2）

-b（y21-y22）|≤

b|y1-y2|+2

bh2 L

|y1-y2|=（b+2 bh2 L ）|y1-y2|.

故f（x，y）在D上满足Lipschitz条件.

∴由常系数微分方程的解存在唯一性定理知，（**）式满足y（0）= L 1+a 的解在|x|≤h上存在且唯一.

这里的h=min{h1，h2}，M=max|f（x，y）|，（x，y∈D）.

（Ⅲ）现在解（**）式：

1）观察f（x，y），得：y=0或y=L为（**）的解.

2）在f（x，y）中，令z= 1 y ，则

dz dx =

dz dy × dy dx =

- 1 y2 · b L y（L-y）

=-bz+

b L .

（1）

再令z=c（x）e-bx，则

dz dx =

dc（x） dx e-bx-

be-bxc（x）

= dc（x） dx e-bx

-bz. （2）

比较（1）（2）得： dc（x） dx = b L ebx.

利用变量分离法，求得：c（x）= 1 L ebx

+C，C为任一常数.

∴y= 1 z =

1 c（x）

ebx=

L ebx+CL ebx=

L 1+CLe-bx .

（3）

将y（0）=

L 1+a

代入（3），得：C= a L .

∴y= L 1+ae-bx .

综上，（**）式的解为y=0或L或 L 1+ae-bx .

但由（Ⅰ）知，样本（xi，yi）对（**）式中参数L，a，b的估计问题同样适用.故实际操作中，由于n个样本（xi，yi）是不同的，所以y=0，L要排除.

综合（Ⅱ）（Ⅲ）知，在整个实数域 R 上，（**）式的解都是存在参数（L，a，b）作最小二乘估计.

第三步，多元线性回归方程的一般形式是y=β0+β1x1+…+βpxp，不失一般性，总可以设β0=0.因为可以形式上引入自变量x0=1，所以在理论研究时，可以假设多元线性回归模型的形式为y=β1x1+…+βpxp.

由假设知，y= L 1+ae-bx 中参数a，b可作最小二乘估计，故y= L 1+ae-bx 必可划为线性回归形式：g（x，y）=p（L）u（x，y）+q（a）u（x，y）+h（b）w（x，y）.（4）

由于我们是先估计出p（L），q（a），h（b）后，利用L=p（L），

q（a）=q（a），h（b）=h（b），进而估计出，L，a，b的，所以（4）式与y= L 1+ae-bx 实质上对L，a，b的估计并没有任何区别.

故不妨设y= L 1+ae-bx 可划为y=Lx1+ax2+bx3（5）.

（5）式中，残差θi=

（yi-Lxi1-axi2-bxi3）=（L，a，b）xx′（Lab）-2y′x（Lab）+y′y.

其中，x=x11x22…xnn，y=y1y2…yn，

∴2xx′（Lab）-2x′y=0.

由假设知，L，a，b可用最小二乘估计出：若（L，a，b）唯一，则r（xx′）=3；若（L，a，b）不唯一，则（L，a，b）属于R3中一真子空间，此时r（xx′）=1或2.

综上，由假设得到r（xx′）>0.

第四步，由第二步知，对y= L 1+ae-bx 中参数L，a，b的估计等价于对 dy dx

=by（1- y L ），y（0）= L 1+a 中参数L，a，b的估计.

分析后者，我们发现通过样本（xi，yi）至多只能直接估出b，L，而a与样本（xi，yi）并没有直接关系，它只与y（0）有直接联系，a= L y（0） -1（∵由第一步知，L≠0，∴y（0）= L 1+a ≠0）.

但实际的样本数据中，由于（0，y（0））并未给出，故y（0）的取值范围为R＼{0}.

现任意取定L，b，一方面残差θ退化为只含a的二次函数，即r，s，t∈R1→θ=θ（a）=ra2+sa+t，

则a的估计值必须满足θ′（a）=0.即2ra+s=0. （6）

另一方面，a=

L y（0）

-1，由L≠0，y（0）的取值范围为R＼{0}，得到a的取值范围为R＼{-1}，故R＼{-1}中任一点都可作为a，显然也要满足（6）式.

∴r=s=0，

∴θ（a）=t为常数，故残差θ与a无关.

同理，任意取定a，b，也可得到残差θ与L无关.

故残差θ至多与b有关.

如果θ与b有关，则θ=θ（b），即y= L 1+ae-bx 可转化为线性回归形式，g（x，y）=h（b）w（x，y）；但前式含有3个参数L，a，b，后式只含有1个参数b，显然，这种转化不成立.

故残差θ=t为常数.

反馈到第三步，可知此时的xx′=03×3，即r（xx′）=0，与假设得到的r（xx′）>0矛盾.

∴原假设不成立，故猜想成立.

证毕.

四、总结与归纳

尽管L，a，b在非限制条件下，y= L 1+ae-bx 中的参数L，a，b不能用最小二乘法估计，但在实际经济模型或生物统计模型中，L，a，b往往具有特殊的经济含义和现实意义.我们可以用其他方法先估计出L或b等，然后再利用最小二乘法估计剩余的两个参数.比如，著名的Logistic人口发展模型y= L 1+ae-bx 中的参数L代表自然资源和环境条件所能容纳的最大人口数量，我们可以先根据生态学的知识预测出L，然后将y= L 1+ae-bx 转化为线性回归形式，即

-bx+lna=ln（ L y -1）

，再利用不同年份的人口统计数据和最小二乘法估计出a，b.

[ 参考文献 ]

[1] 王高雄.周之铭等.常微分方程[M].北京：高等教育出版社，2007.

[2] 北大数学系几何与代数教研室前代数小组.高等代数[M].北京：高等教育出版社，2003.

[3] 陈家鼎，孙山泽等.数量统计学讲义[M].北京：高等教育出版社，2011.

[4] 吴梅村.数理统计学基本原理和方法[M].成都：西南财经大学出版社，2006.

最小二乘问题篇4

根据电力系统负荷预测的层级化特点,文献[1]提出了多级负荷预测的概念和研究框架。基于这个基本框架,后续研究的多级负荷预测的基本协调模型[2]和关联协调模型[3]能够较好地解决各种具有直接加和特性的一维两级和多维多级预测协调问题。这2种协调模型都是借鉴电力系统状态估计思想[4]建立的,量测量可以通过量测方程用状态量表示出来,量测方程就是反映电网物理特性的潮流等网络方程[4,5]。多级负荷预测的协调问题也存在类似特点。对于多级负荷预测,只需要部分负荷值就可以反映全部负荷状况,本文称这组负荷值为状态负荷。多级负荷预测的各预测值都是量测量,它们也可以通过一组量测方程用状态负荷表示出来,本文称这组方程为负荷量测方程。引入状态负荷和负荷量测方程的概念后,多级负荷预测协调问题即可表示为典型的状态估计的形式[4,5]。

正如文献[1]所指出的,在负荷预测中总需求和子需求的对照几乎无处不在,预测问题呈现明显的“多维多级”特征。这里的“多维多级”包括多个时间周期、多种空间划分、多个行政级别、多种属性分类,等等。姑且从预测主体角度看,国家电网、区域电网、省网以及地区电网等都需要进行相应的负荷预测,而由于电网级别、负荷水平、负荷构成、负荷预测方法以及负荷预测人员等的不同,负荷预测精度不同,负荷预测结果的可信程度一般也不相同。由此可知,对于多级负荷预测协调问题,无论是基本协调模型,还是关联协调模型,都涉及一个重要的物理量——预测可信度。如何确定合理的预测可信度成为电力系统多级负荷预测及其协调问题的重要组成部分。既然多级预测协调模型是借鉴状态估计思想建立的,那么预测可信度是否也可以采用类似状态估计权重的形式?或者有其他更好的可信度方式?必须对此问题展开探讨。

本文首先分析多级负荷预测协调模型的特征以及加权最小二乘(WLS)准则下的可信度,然后针对WLS准则下可信度的不足提出约束最小二乘可信度优化模型,并分析了基于约束最小二乘可信度的多级协调模型提高预测精度的机理及其精度评价。

1 加权最小二乘准则下的可信度

考虑如下所示的两维两级协调问题:

虽然该问题总共有(m+1)×(n+1)个负荷变量,但只需要m×n个状态负荷值就可以完全反映此两维两级协调的负荷情况,其他负荷值可以通过负荷量测方程表示。对于状态负荷,其负荷量测方程就是其自身;对于非状态负荷,其负荷量测方程为其相关负荷的某种组合。

如果取各子需求负荷为状态负荷,则非状态负荷(各总需求负荷)的负荷量测方程为其相关状态负荷的加和,即关联协调模型[3]的约束条件。用数学方程表示为:

${\begin{cases} x_{i, j} = x_{i, j} i = 1, 2, \dots, m; j = 1, 2, \dots, n \\ x_{i, 0} = \sum_{j = 1}^{n} x_{i, j} i = 1, 2, \dots, m \\ x_{0, j} = \sum_{i = 1}^{m} x_{i, j} j = 0, 1, \dots, n \end{cases} (2)$

式(2)所示的所有方程统称为负荷量测方程,其向量记为h(X),其中X为m×n维状态负荷向量。式(2)中第1个式子看似多余,但其必须包含在h(X)中,才能在下文导出式(3),即典型的状态估计形式。引入状态负荷和负荷量测方程概念后,多级负荷预测协调问题可表示为典型的状态估计形式[4,5]。

$\min J (X) = (Ζ - h (X))^{Τ} R^{- 1} (Ζ - h (X)) (3)$

式中:Z为(m+1)×(n+1)维负荷预测值列向量;h(X)为对应的负荷量测方程向量;R-1=Q=diag(W)/(diag(Z))2为权重。

式(3)所示的数学模型相当于矩阵形式的关联协调模型[3]中把约束条件代入目标函数。因此,其解与关联协调模型相同。另外,由于基本协调模型是关联协调模型的简化[3],因此,基本协调模型也可以转换为式(3)所示的形式。

式(3)所示的模型是状态估计器中最典型的WLS估计器。它的特点是:对量测误差为理想正态分布时,估计具有最优一致无偏等优良统计特性[6,7,8]。模型(3)中的权重一般取为方差的倒数[4,5,9],这样选取可以使WLS具有最小方差性质[8,10]。因此,对多级负荷预测可以采用类似的权重形式,但预测误差的方差是未知的,需要根据历史预测误差估计。而负荷预测,特别是中长期负荷预测相邻2次预测的时间差较大,相邻2次预测值的差别较大,所以多级协调中各预测值的误差方差估计不能简单统计历史预测绝对误差的方差,而应根据更有可比性的历史相对误差的方差来估计,考虑到相对误差与绝对误差的关系以及方差的统计公式,得出如下的多级协调可信度权重形式(以基本协调模型为例):

$q_{i} = \frac{w_{i}}{z_{i}^{2}} = \frac{1}{z_{i}^{2} σ_{i}^{2}} i = 0, 1, \dots, n (4)$

式中:σ $_{i}^{2}$ 为第i项电力需求的预测相对误差的方差。

由式(4)可得预测可信度为:

$w_{i} = \frac{1}{σ_{i}^{2}} i = 0, 1, \dots, n (5)$

为更便于比较各可信度的大小,可将式(5)所示的可信度归一化:

$w_{i} = \frac{\frac{1}{σ_{i}^{2}}}{\max_{0 \leq i \leq n} \frac{1}{σ_{i}^{2}}} i = 0, 1, \dots, n (6)$

WLS准则下的可信度能够反映预测精度,但基于此可信度的多级负荷预测协调模型却有可能出现调整过大的现象,即多级负荷预测协调模型的相对调整系数[2]有可能小于-1。其详细数学推导参见附录A。

2 约束最小二乘可信度

为了保持加权最小二乘可信度反映预测精度的优点,同时还要避免可能出现的调整过大的现象,本文考虑对其施加一定的约束,构造约束最小二乘可信度。本文主要分析基本协调模型下的约束最小二乘可信度,关联协调模型可以以此类推。

基本协调模型下,最简单的约束最小二乘可信度思路是让子需求相对调整系数都不小于-1,但这样的做法不可行。假设一种极端的情况:子需求预测值相差不多,总需求预测值大于各子需求预测值之和,总需求历史预测相对误差方差非常小,子需求预测相对误差方差相差不多。这种情况下,基本协调模型将类似于自上而下式协调[2]。如果不对各预测项的可信度作调整,由于总需求预测值大于各子需求预测值之和,子需求相对调整系数应该小于-1,所以在这种假设情况下,把子需求相对调整系数限制在不小于-1是不合理的,即这样的一组可信度值不存在。

因此,为了保证基本协调模型不调整过大,基本协调方法对各子需求预测的调整量应限制在不大于自上而下式协调对各子需求预测的调整量[2]。即相对调整系数:

$δ_{i}^{r} = - \frac{\frac{1}{w_{i}} \frac{z_{i}}{z_{0}}}{\sum_{i = 0}^{n} \frac{1}{w_{i}} (\frac{z_{i}}{z_{0}})^{2}} \geq - \frac{1}{\sum_{i = 1}^{n} \frac{z_{i}}{z_{0}}} i = 1, 2, \dots, n (7)$

如对第k个可信度,式(7)可变形为:

$w_{k} \geq \frac{z_{k} (\sum_{i = 1}^{n} z_{i} - z_{k})}{\sum_{i = 0, i \neq k}^{n} \frac{z_{i}^{2}}{w_{i}}} k = 1, 2, \dots, n (8)$

由式(8)可见,某个可信度的临界值是其他可信度的函数。如果增大基本协调模型中的某个可信度,其他可信度的临界值也将相应增大,这就有可能违背了其他可信度的约束条件。所以,若有可信度不满足其约束条件,必须做到各可信度协调调整。为此,构造如下的可信度优化调整模型:

${\begin{cases} \min f = \sum_{k = 1}^{n} (w_{k} - \frac{\frac{1}{σ_{k}^{2}}}{\max_{0 \leq k \leq n} \frac{1}{σ_{k}^{2}}})^{2} \\ s . t . w_{k} \geq \frac{z_{k} (\sum_{i = 1}^{n} z_{i} - z_{k})}{\sum_{i = 0, i \neq k}^{n} \frac{z_{i}^{2}}{w_{i}}} k = 1, 2, \dots, n \\ w_{k} \geq 0 k = 1, 2, \dots, n \end{cases} (9)$

数学模型(9)必然有可行解吗?回答是肯定的。简要分析如下:当总需求预测完全不可信,即所有子需求预测的可信度远大于总需求预测时,数学模型(9)的约束条件必然成立,而这种可信度形式正是传统的自下而上式协调;当子需求预测的可信度与其基数规模成正比,总需求预测的可信度远大于各子需求预测时,数学模型(9)的约束条件刚好处于边界上,而这种可信度形式正是传统的自上而下式协调。因此,传统的自下而上式协调和自上而下式协调的可信度形式,分别是约束最小二乘可信度的上下界,边界及界内的所有可信度组合都是式(9)的可行解。其最优解可以采用罚函数法求解[11],此处不再赘述。

3 模型机理与精度评价

3.1 模型机理

基于约束最小二乘可信度的多级负荷预测协调模型能够通过合理调整原始预测值提高预测精度,主要有以下3方面的原因:

1)多级负荷预测协调模型建立在多级电力需求关联性的基础上,多维多级电力需求的冗余量和相互参照作用是预测精度提高的内在机制。这一点与电力系统状态估计的思想完全一致。

2)约束最小二乘可信度反映了历史预测精度,精度越高则调整越小。基于预测精度的可信度形式能够在概率意义上发现导致多级电力需求不平衡的原因并对其调整,从而使得对各预测值的调整有章可循,是预测精度提高的重要保障。

3)约束最小二乘可信度能够避免对某些预测项调整过大,使得对预测值的调整更具理性,在一定程度上防止矫枉过正,进一步保证预测精度的提高。

基于约束最小二乘可信度的多级负荷预测协调模型从以上3个层面的内在机理上确保了其对多级负荷预测精度的提高。

3.2 多级负荷预测的精度评价

本文采用如下2个指标整体评价多级负荷预测问题的预测精度。

1)总体相对误差:

是指所有预测项误差绝对值之和与总需求实际值的比,用于评价多级负荷预测问题的总体误差状况。其计算公式为:

$E_{Τ} = \frac{\sum_{i = 0}^{n} | z_{i} - x_{i} |}{x_{0}} \times 100 ％ (10)$

式(10)还可变形为:

$E_{Τ} = \sum_{i = 0}^{n} (\frac{x_{i}}{x_{0}} \frac{| z_{i} - x_{i} |}{x_{i}}) \times 100 ％ (11)$

由式(11)可见,总体相对误差是所有预测项相对误差的加权和,其权重为各预测项实际值与总需求实际值的比。对于2个相对误差相同的预测项而言,实际值越大的预测项对总体相对误差的贡献越大。从这个意义上讲,总体相对误差也反映了考虑预测项基数规模的多级负荷预测问题的相对误差状况。

2)平均相对误差:

是指所有预测项相对误差的平均值,用于评价多级负荷预测问题的平均误差状况。其中没有考虑各预测项的基数规模。其计算公式为:

$E_{A} = \frac{1}{n + 1} \sum_{i = 0}^{n} \frac{| z_{i} - x_{i} |}{x_{i}} \times 100 ％ (12)$

4 算例分析

以某省2000年—2004年的数据作为基础数据,对2005年进行负荷预测。全年和各月全社会用电量构成基本协调问题,2005年各物理量的实际值、预测值和预测误差如表1所示。

对2002年—2004年进行虚拟预测得到的历史预测误差也列于表1中。

以历史预测误差方差的倒数为可信度(即WLS准则下的可信度)的基本协调结果如表2所示。由表2可见,最小二乘可信度下的基本协调结果中,2月预测值的调整百分比达到了-5.8%,使得其预测误差由-4.5%变为了-10.1%,调整过大,远大于自上而下式协调的调整百分比-2.9%。同样,3月和10月的调整量也偏大。

为了解决上述问题,采用本文提出的约束最小二乘可信度代替最小二乘可信度,其协调结果也列于表2。约束最小二乘可信度优化方法在尽量少调整原有可信度的原则下适当增大了2月、3月和10月的可信度,从而减小了其预测值调整百分比。另外,还可以发现,在调整这3个月的可信度的同时,其他月份的可信度也进行了相应的调整。从图1中能更直观地看出约束最小二乘可信度的优化效果。

对本算例各种协调方式的精度评价如表3所示。由表3可见,最小二乘可信度和约束最小二乘可信度下的基本协调都能明显减小基本协调问题的总体相对误差和平均相对误差;而且约束最小二乘可信度在不出现调整过大现象的同时,其协调精度也优于最小二乘可信度;而传统自上而下式协调虽然能达到总体平衡,但总体相对误差和平均相对误差都比原始预测大。

5 结语

本文分析了多级负荷预测协调问题与状态估计的相似性,由此推导出了WLS准则下的多级协调预测可信度形式。由于WLS准则下的可信度形式使得多级负荷预测协调可能出现调整过大现象,为此,本文研究了约束最小二乘可信度优化模型。在此基础上,分析了基于约束最小二乘可信度的多级协调模型提高预测精度的机理。最后,通过实例验证了基本协调模型约束最小二乘可信度的有效性及其提高多级预测精度的实效。约束最小二乘可信度,既能反映预测精度,又可以有效避免对预测项的调整过大现象,为多级负荷预测协调模型提供了一种有效的可信度确定方式。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

参考文献

[1]康重庆,牟涛,夏清.电力系统多级负荷预测及其协调问题:(一)研究框架.电力系统自动化,2008,32(7):34-38.KANG Chongqing,MU Tao,XI A Qing.Power system multilevel load forecasting and coordinating:Part one research framework.Automation of Electric Power Systems,2008,32(7):34-38.

[2]牟涛,康重庆,夏清,等.电力系统多级负荷预测及其协调问题:(二)基本协调模型.电力系统自动化,2008,32(8):14-18.MU Tao,KANG Chongqing,XI A Qing,et al.Power system multilevel load forecasting and coordinating:Part two basic coordinating model.Automation of Electric Power Systems,2008,32(8):14-18.

[3]牟涛,康重庆,夏清,等.电力系统多级负荷预测及其协调问题:(三)关联协调模型.电力系统自动化,2008,32(9):20-24.MU Tao,KANG Chongqing,XI A Qing,et al.Power system multilevel load forecasting and coordinating:Part three correlative coordinating model.Automation of Electric Power Systems,2008,32(9):20-24.

[4]于尔铿.电力系统状态估计.北京:水利电力出版社,1985.

[5]于尔铿,刘广一,周京阳,等.能量管理系统.北京:科学出版社,1998.

[6]李碧君,薛禹胜,顾锦汶,等.电力系统状态估计问题的研究现状和展望.电力系统自动化,1998,22(11):53-60.LI Bijun,XUE Yusheng,GU Jinwen,et al.Status quo and prospect of power system state estimation.Automation of Electric Power Systems,1998,22(11):53-60.

[7]SCHWEPPE FRED C,WILDES J.Power system static-state estimation:PartⅠexact model.IEEE Trans on Power Apparatus and Systems,1970,89(1):120-125.

[8]夏天昌.系统辨识——最小二乘法.北京:国防工业出版社,1984.

[9]兰华,李积捷.电力系统状态估计算法的研究现状和展望.继电器,2007,35(10):78-82.LAN Hua,LI Jijie.Status quo and prospect of algorithm of power systemstate estimation.Relay,2007,35(10):78-82.

[10]国家测绘总局测绘研究所.近代最小二乘法:译文集.北京:测绘出版社,1980.

最小二乘问题篇5

基于偏最小二乘回归的飞机维修保障费用预测

分析了影响飞机维修保障费用的.参数,提出用偏最小二乘回归方法来预测飞机维修保障费用.该方法对变量进行主成分分析、典型相关分析和多元线性回归,在处理存在多重线性相关的小样本多元数据方面效果很好.实例证明,与传统普通多元线性回归方法相比,偏最小二乘回归在飞机维修保障费用预测中精度更高.

作者：郭风张恒喜李寿安张琦 GUO Feng ZHANG Heng-xi LI Shou-an ZHANG Qi 作者单位：空军工程大学,工程学院,陕西,西安,710038刊名：空军工程大学学报(自然科学版) ISTIC PKU英文刊名：JOURNAL OF AIR FORCE ENGINEERING UNIVERSITY (NATURAL SCIENCE EDITION)年，卷(期)：20056(3)分类号：V37关键词：维修保障费用预测偏最小二乘回归

最小二乘问题篇6

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

二、方法介绍

如果在实际问题中出现了多重共线性的现象,我们可以选择用有偏回归方法———岭回归( 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的结果可知,两种方法建立的回归方程拟合的效果都不错,而偏最小二乘回归方法相对岭回归方法要更优.

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

最小二乘问题篇7

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

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

最小二乘问题篇8

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面。其中一种简单的情形是:根据相关影长数据确定参照物位置, 这一技术也称为太阳影子定位技术。在本文中, 根据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.

广义逆的计算与最小二乘估计篇9

在系统理论、信号处理和矩阵方程求解等许多应用领域中,往往需要给出形如:

的矩阵的某种类型的广义逆,其中表示样本数目的k通常都很大,而且随时间增长而增大,因而计算量也迅速增加,这在实时应用中是应该避免的。

本文给出了某些类型的广义逆的算法,并用所得结果,给出了求解由多个矩阵方程所构成的系统的最小二乘估计,从而推广了[1]的工作,并简化了[2]中有关求多个方程公共解的算法,文中使用广义逆的记号均引自[3]。

1广义逆的计算

设k≥1为一整数,为具有相同行数的矩阵且满足

再设的某种广义逆,我们将证明,存在着满足

的矩阵Mk和Gk,使得

为的同种广义逆。

1°,如果(·)-表{1}-逆;

2°,如果(·)-表{1,2}-逆;

1°的证明:注意到Qk-1为对称幂等矩阵且Ak=(AkTAk)(1)AkT=(AkTAk)(1)HkTQk-1[3],我们有

2°的证明:注意到对任意A(1,3)∈A{1,3},总存在着矩阵Z,使得A(1,3)=A-+(I-A-A)Z,因此AA-=AA+,于是由(1),(2)和(3)可知

为对称阵,另外由1°可知:

3°的证明:由2°,我们只需再进一步证明

引理3[5]如果Gk如下确定:

从引理1到引理3知,为获得的广义逆,必须求出

容易看出,以上两式的计算量将随着k的增大而迅速增大,为避免这一困难,注意到

或等价地,

显然,(5)～(8)的计算量与k无关.

综上所述,我们证明了以下两定理。

定理1设Mk和Zk由(1),(2)和(3)所定义且对k=0,1,…,令

Q-1=I,Ak=Qk-1Hk,Qk=Qk-1-AkGk

则

定理2设Mk和Zk由(1),(2)和(3)所定义且对k=0,1,…,令

下面,我们将定理1与定理2应用于

其中和Hk为具有相同列数的矩阵,且注意到如下事实[6]

我们有

定理3设由(9)所定义,且令

则

定理4设和Zk由(9)和(10)所定义且对k=0,1,…,令

则

2最小二乘估计

考虑观测系统

其中Yk∈Rm×p,Hk∈Rm×n为已知,Vk∈Rm×p为零均值随机矩阵,X为待估未知矩阵,对任何k,称X赞k为X的最小二乘估计(LSE),如果X赞k为最优化问题

的最小范数解,而

及加权矩阵wi∈Rm×n皆正定,i=0,1,…,k.

首先考虑无加权情形,即的情形。

为使估计误差的方差最小[6,7],我们有

其中的Moore-Penrose逆.

为得到(14)的递推式,注意到定理4以及(9)、(10)和(13),我们有

定理5系统(11)的X的无加权的极小范数最小二乘估计可如下递推地给出

其中Gk由定理4递推而得。

证明:

如果m=p=1,则定理5便退化为[1]的结果。定理5还有如下简单的推论:

若(11)中的Vk=0,坌k,则定理5给出多个方程所构成的系统

的公共解,亦即我们有

则(16)的公共解可按

递推地求得。

该推论中要求计算Ak+,考虑广义逆的性质和定理3,如果Gk=Qk-1Ak(1),则(17)给出(16)的公共解,从而只需计算A(1)∈A{1}而使计算量大为减少。类似地,如果Gk=Ak(1,4),则(17)给出(16)的具有极小范数的公共解。

现在再考虑一般加权最小二乘估计,因

故可通过无加权问题(18)的右端来处理加权问题。其中

或等价地,

应用定理5于无加权问题(18),便有

定理6令P-1=0,Q-1=I,且对k=0,1,…令

H′kTCk,则系统(11)的加权极小范数二乘估计可如下递推地求得:

为使(11)的加权最小二乘估计是最优估计[8],其中加权矩阵应如何选择?答案如下:

设而其方差设为

则由[4]中的推论6.4.4知,当时,为(11)在给定Yo,…,Yk时的最优最小二乘估计。

关于定理3和定理4的更多应用,可进一步参考[6,9,10]等。

3结论

应用定理1、2、3的结果,将[1]中关于单变量情形下的最小二乘估计推广到了多变量的情形,给出了后者的最优最小二乘估计的算法,同时还简化了[2]中求多个相容矩阵方程公共解的算法。

摘要：给出了某些广义逆的计算方法,将这些结果应用于观测系统的参数估计,给出了系统无加权和加权最小二乘的最新估计结果,推广了[1]的工作,同时简化了[2]中求多个相容矩阵方程公共解的算法。

关键词：算法,最小二乘估计,广义逆,方差,加权

参考文献

[1]Albert A,Sittle R.W.A method for computing leastsquares estimators that keep up with data SIAM[J].Control,1965,3:384-417.

[2]Morris G.L,Odell P.L.Common solutions for a matrixequations with applications[J].Assoc comput.math,1968(15):272-274.

[3]A.B.Israel,T.N.E.Greville.Generalized inverses:Theory and Applications[M].New York:John Wiley,1974.

[4]Campbell S.L,Meyer,C.D.Jr.Generalized Inverses of lin-ear transformations[M].London:Pitman publishing Ltd,1979:104-115.

[5]Cline R.E.Representations for the generalized inverses ofa partitioned matrix[J].Soc.Indust.Appl.Math.1964(12):588-600.

[6]X.Liu(刘轩黄).Recursive computation of generalizedinverses with applications to optimal state estimation[J].Control theory and advanced technology.1995(10):1485-1497.

[7]Rao G.R.Linear statistical inference 2nd its applications(2nd edtion)[M].New York:John Wiley&Sons,1973.

[8]刘轩黄.无初始状态统计知识时多变量系统的最优平滑与预报[J].海南大学学报,1998,16(3):210-214.

[9]X.Liu(刘轩黄).Optimal and minimum energy optimaltracking problems[J].苏州科技学院学报,2005,22(1):1-9;22(2):9-18.

最小二乘问题篇10

空间直线拟合在工程中具有重要应用[1,2,3,4,5,6,7,8],如在计算机视觉[4]、三坐标测量问题[5]、机器人焊缝中心三维定位[6]和空间管线网络干线的确定[7]中。

由于最小二乘准则下问题求解比较容易,一些文献在该准则下给出确定空间直线的算法。文献[1]将直线理解为两个平面的交集,通过求取两个平面方程,从而确定直线。文献 [2,3]根据最佳平方逼近原理,证明直线必过数据中心,并用逐步迭代给定方向向量,从而确定点向式直线方程。文献 [5]则通过两次降维,先按全最小二乘准则拟合二维平面,将空间中的点投影到该平面上,再将该平面上数据利用全最小二乘准则拟合直线,从而确定空间直线。文献[6]先确定直线经过的点,再求导得到一个非线性方程组,解方程组确定方向向量,从而最终确定空间直线。这些算法比较复杂,本文给出一种基于主成分分析的简便算法,通过计算机仿真验证了算法的有效性。

由于最小一乘准则下确定空间直线是一个不可微优化问题,难度较大。文献 [7]提出了空间中按全加权最小一乘准则下确定直线问题,但没有给出具体的算法,文献[8,9]给出了平面上的求解方法,没有对空间中全加权最小一乘准则下确定直线问题进行讨论。本文给出一种基于无约束优化的搜索算法,并进行了计算机仿真。

1 空间直线拟合问题

空间直线拟合问题可以描述为:在三维空间给定n个点P1(x1,y1,z1),P2(x2,y2,z2),…,Pn(xn,yn,zn),求空间直线L: $\frac{x - x_{0}}{m} = \frac{y - y_{0}}{n} = \frac{z - z_{0}}{p}$ ,其中M(x,y,z)为直线上任意点,M0(x0,y0,z0)为直线上确定的一点,非零向量 $\vec{s} = (m, n, p)$ 为直线的方向向量,使得目标函数:

$\sum_{i = 1}^{n} d^{2} (Ρ_{i}, L) (1)$

为最小,其中 $d (Ρ_{i}, L) = \frac{| \vec{Μ_{0} Ρ_{i}} \times \vec{s} |}{| \vec{s} |}$ 表示点Pi到直线L的距离[10] ,其中|·|表示向量的模。式(1)实际上是按全最小二乘拟合空间直线,等价于文献[1,2,3,4,5,6]中拟合空间直线的准则。

如在上述假设下,改为求目标函数:

$\sum_{i = 1}^{n} d (Ρ_{i}, L) (2)$

最小,则是按全最小一乘拟合空间直线,这样得到的直线不易受奇异点(outliers)的影响。

如在上述假设下,改为求目标函数:

$\sum_{i = 1}^{n} w_{i} d (Ρ_{i}, L) (3)$

最小,其中wi(i=1,2,…,n)为权重,则是按加权全最小一乘拟合空间直线,这样得到的直线可以用于空间中管道的主干线确定问题。当wi=1(i=1,2,…,n)时,式(3)变为式(2)。

2 拟合空间直线的算法

2.1 全最小二乘准则下的拟合算法

将空间中三维分布的点集拟合为一维的直线,其本质就是降维,因此可以按照主成分的思想,对三维分布的空间数据进行主成分分析,提取第一主成分向量作为直线方向向量。文献[3]证明在最小二乘准则下,空间直线必通过数据中心点,因此可以由点向式确定空间直线方程。按此思想给出求问题式(1)中直线的算法如下:

第1步输入三维空间中点坐标的n个数据点矩阵Xn×3,求出均值 $\bar{X} = (\bar{x}, \bar{y}, \bar{z})$ ,拟合直线必过数据中心 $(\bar{x}, \bar{y}, \bar{z})$ ;

第2步求协方差矩阵 $S = \frac{1}{n} (X - \vec{1}_{n} \bar{X})^{Τ} \cdot (X - \vec{1}_{n} \bar{X})$ ,其中1n为n行1列元素全为1的n维列向量;

第3步求出S的最大特征值λ1和对应特征向量 $\vec{v}_{1}$ ,以 $\vec{s} = \vec{v}_{1} = (a, b, c)$ 作为直线方向向量;

第4步由点向式确定直线 $\frac{x - \bar{x}}{a} = \frac{y - \bar{y}}{b} = \frac{z - \bar{z}}{c}$ ,其中(x,y,z)为直线上任意一点。

由主成分意义,点到直线距离平方和为n(λ2+λ3),其中λ2与λ3为S的其余两个特征值,与由式(1)得到的结果相同。

2.2 加权全最小一乘准则下的拟合算法

由于问题式(2)是式(3)的特殊情形,因此讨论式(3)。在问题式(3)中目标函数中含有绝对值,是一个不可微函数,因而求解比较困难。

要在空间中确定一条直线,无论采用确定两个相交平面还是点向式方程都至少需要确定6个未知参数。

将问题式(3)中的目标函数看作是一个含有6个自变量的不可微无约束优化问题,可以采用文献[11]中不依赖于导数的Nelder-Mead的单纯形算法与Hooke-Jeeves的模式搜索法求解。MATLAB中有基于单纯形算法的fminsearch()命令可以直接调用。

本文采用点向式方程确定直线,约定6个自变量按下面方式排列 $\vec{v} = (x_{0}, y_{0}, z_{0}, a, b, c)$ ,给出求问题式(3)中直线的算法如下:

第1步输入三维空间中点坐标的n个数据点矩阵Xn×3,给出初始值,给出 $\vec{v} = (x_{0}, y_{0}, z_{0}, a, b, c)$ 的初始值。初始值可以选为全最小二乘准则下给出的值。

第2步利用初始值和MATLAB中fminsearch()对目标函数 $\sum_{i = 1}^{n} w_{i} d (Ρ_{i}, L)$ 进行迭代求解,给出收敛值,从而确定直线经过的点(x0,y0,z0)和方向向量(a,b,c)。收敛的最终结果(x0,y0,z0)不唯一,方向向量(a,b,c)经标准化并指定方向后唯一。

第3步由点向式确定直线 $\frac{x - x_{0}}{a} = \frac{y - y_{0}}{b} = \frac{z - z_{0}}{c}$ ,其中(x,y,z)为直线上任意一点。

2.3 空间点到拟合直线的投影

为了在计算机仿真中绘制空间中点到拟合直线的垂线,这里给出点到直线的垂足坐标确定方法。

由于空间拟合直线已经求出,故可以确定空间中每一个点Pi(xi,yi,zi)在直线上的垂直投影点P′i(x′i,y′i,z′i)。由空间解析几何知识[10]可知投影点P′i(x′i,y′i,z′i)为直线在经过点Pi(xi,yi,zi)且与直线垂直平面上的垂足,联立直线的参数式方程与直线垂直的平面方程。

${\begin{cases} x = \bar{x} + a t y = \bar{y} + b t z = \bar{z} + c t \\ a (x - x_{i}) + b (y - y_{i}) + c (z - z_{i}) = 0 \end{cases} (4)$

可以解出 $t = \frac{a (x_{i} - \bar{x}) + b (y_{i} - \bar{y}) + c (z_{i} - \bar{z})}{a^{2} + b^{2} + c^{2}}$ ,从而代入参数式方程解出(x,y,z),即为垂足P′i(x′i,y′i,z′i)。

3 计算机仿真

为了检验本文算法的正确性和有效性,给出下面例子。

仿真环境为 Intel酷睿2双核 P8400、2G内存、Windows Vista OS、MATLAB7.5环境。

例1 设直线方程为x=1+t,y=1+t,z=1+t,t为参数。

参数t取值区间为[-6,6],先将该区间等分为30个小区间,t取区间端点,从而可在直线上确定31个点,可以确定31×3的数据矩阵。对直线上每个点进行随机扰动,让x,y,z分别加上随机数种子为rand(‘twister’,10),rand(‘twister’,20),rand(‘twister’,30) 的服从[-1,1]上均匀分布的n维随机列向量。绘出这些点三维散点图,如图1所示。

表1给出不同准则和算法得到的数值结果。其中方向向量已经标准化并约定向上(c>0),加权最小一乘准则中取权重全为1,即最小一乘,优化使用MATLAB自带的fminsearch()。

按最小二乘准则,绘制出拟合直线和点到拟合直线的垂线,如图2所示。

按最小一乘准则绘制出拟合直线和点到拟合直线的垂线,如图3所示。

4 结论

本文讨论了空间中的直线拟合问题,针对全最小二乘和全加权最小一乘准则下,分别给出了基于主成分分析和无约束优化的空间直线拟合方法,这些方法在计算机上实现起来比较简便,计算机仿真验证了算法的有效性。这对使用空间直线拟合方法的科研人员具有一定的参考价值。

空间中曲线拟合比较复杂,未见到统一的方法,然而特殊情况比较容易。如空间中拟合圆周,可以先拟合平面,将三维数据投影到拟合平面上(降维到2),再根据平面法向量的方向角对平面上投影数据旋转使得一个分量全部为0,将旋转后的数据拟合圆周,得到圆周对应的圆心和半径。最后进行旋转逆变换,将球面方程与平面方程联立即得空间圆周方程。

参考文献

[1]袭杨.空间直线拟合的一种方法[J].齐齐哈尔大学学报,2009,25(2):64-68.

[2]许海涛,高彩霞.直线拟合算法[J].电脑知识与技术,2009,5(4):864-867.

[3]杜明芳.空间直线拟合[J].北京印刷学院学报,1996,4(2):27-31.

[4]李畅,张剑清,邹松柏,等.自动识别街道立面凹凸边界[J].计算机工程与应用,2009(15):6-8.

[5]胡学军,王俊亮,王刚.全最小二乘法在三坐标测量中的应用[J].武汉工程大学学报,2008,30(4):112-113.

[6]郭学军,王喜平,刘贺平.机器人焊缝中心三维定位的解析方法[J].数学的实践与认识,2009,39(9):123-127.

[7]梁怡,吴可法.一类非光滑最优化问题的有限步解法[C].中国运筹学会第六届学术交流会论文集(章祥荪等主编),Global-Link出版公司,香港,2000,801-811.

[8]林诒勋,尚松蒲.平面上的点—线选址问题[J].运筹学学报,2002,6(3):61-68.

[9]冯守平,杨桂元.关于加权全最小一乘法[J].应用概率统计,2009,25(2):135-142.

[10]同济大学应用数学系.高等数学(第五版上册)[M].北京:高等教育出版社,2002.

【最小二乘问题】推荐阅读：