在水电站优化运行中，水轮机的动力特性曲线，需要根据离散的特征点拟合得到。拟合方法有多项式拟合、径向基函数神经网络拟合等，其中多项式最小二乘拟合，是最经典、最简单的一种方法。本文在水电站优化运行的应用背景下，研究了多项式最小二乘拟合的算法，给出了 Delphi/Pascal 实现。

1 多项式最小二乘拟合

2 Cholesky分解法解方程组

3 数值试验

参考文献

源程序：（基于 Delphi 实现）

多项式最小二乘拟合

Cholesky分解法解线性方程组

1 多项式最小二乘拟合

算法 1 多项式最小二乘拟合

2 Cholesky分解法解方程组

算法2 Cholesky分解法解方程组Ax=b, 参见[1, P60]

输入：对称正定矩阵A，向量b

3 数值试验

采用文献[1,P304]中例3进行测试，PolyFit计算得结果为(80.34762, -0.12429)，而[1]中结果为(80.3463，-0.12427)，有较大出入。为确定正确结果，采用通用科学计算自由软件 GNU Octave 计算，得到结果如下图所示，与PolyFit结果相同，从而算法实现的正确性得到了初步验证。要进一步测试本程序的稳定性可以找更多有代表性的例子，与通用计算软件进行比较。

参考文献

[1] 郑慧娆，陈绍林，莫忠息，黄象鼎. 数值计算方法[M]. 武汉大学出版社，2002.

源程序：（基于 Delphi 实现）

多项式最小二乘拟合

unit PolyFit;interfaceusesVectorMatrix;typeTPolyCurveFit = classprivateFN: Integer;    // 拟合多项式阶数；FM: Integer;    // 离散数据点数；FX: TVector;    // 离散点横坐标；FY: TVector;    // 离散点纵坐标；FPolyCoef: TVector;  // 拟合多项式系数；publicconstructor Create(aOrder, aDataNum: Integer; aDataX, aDataY: TVector);procedure Compute;publishedproperty Order: Integer read FN;property PolyCoef: TVector read FPolyCoef;end;implementationusesMath,Cholesky;constructor TPolyCurveFit.Create(aOrder, aDataNum: Integer; aDataX, aDataY: TVector);
beginFN := aOrder;FM := aDataNum;SetLength(FX, FM);SetLength(FY, FM);FX := Copy(aDataX);FY := Copy(aDataY);
end;procedure TPolyCurveFit.Compute;
varA, AA: TMatrix;   // A, A^T * A;Ab: TVector;      // A^T * b, b即向量FY;AAElem: TVector;  // AA 的所有元素；I, J, K, L: Integer;
begin// 计算矩阵 A;SetLength(A, FM, FN + 1);for I := 0 to FM - 1 dobeginA[I, 0] := 1;for J := 1 to FN dobeginA[I, J] := A[I, J-1] * FX[I];end;end;// 计算矩阵 AA；SetLength(AAElem, 2*FN + 1);for L := 0 to 2*FN dobeginI := Floor(L / 2.0);J := L - I;AAElem[L] := 0;for K := 0 to FM-1 doAAElem[L] := AAElem[L] + A[K, I] * A[K, J]end;SetLength(AA, FN + 1, FN + 1);for I := 0 to FN dofor J := 0 to FN  doAA[I, J] := AAElem[I + J];// 计算向量 Ab；SetLength(Ab, FN + 1);for I := 0 to FN dobeginAb[I] := 0;for  K := 0 to FM - 1 doAb[I] := Ab[I] + A[K, I] * FY[K];end;// Cholesky 分解法解线性方程组 AA * x = Ab；SetLength(FPolyCoef, FN + 1);FPolyCoef := Copy(CholeskySolve(FN + 1, AA, Ab));
end;end.

Cholesky分解法解线性方程组


unit Cholesky;interfaceusesVectorMatrix;function CholeskySolve(aN: Integer; A: TMatrix; b: TVector): TVector;implementation{ Cholesky分解解线性方程组 Ax=b ；输入：aN：阶数；A：系数矩阵，为 aN*aN 阶正定矩阵；b：aN阶向量；返回：aN阶解向量，TVector类型，即数组；
}
function CholeskySolve(aN: Integer; A: TMatrix; b: TVector): TVector;
varX: TVector;I, J, K: Integer;TmpNum: Double;
beginif aN < 2 thenExit;// Cholesky分解 A = G * G^T, G 保存在A的下三角部分;// 计算 G 首列（第0列）；if A[0,0] <= 0  thenExit;A[0, 0] := Sqrt(A[0, 0]);for I:= 1 to aN - 1 doA[I, 0] := A[I, 0] / A[0, 0];// 计算 G 其他列；for K := 1 to aN - 1 dobeginTmpNum := 0;for J := 0 to K-1 doTmpNum := TmpNum + A[K, J]* A[K, J];TmpNum := A[K, K] - TmpNum;if TmpNum <= 0 thenExit;A[K, K] := Sqrt(TmpNum);if K < aN-1 thenfor I := K+1 to aN-1 dobeginTmpNum := 0;for J := 0 to K-1 doTmpNum := TmpNum + A[I, J] * A[K, J];{ 2012-04-15，hll严重bug修正, PolyFit测试程序中二阶方程测试成功为巧合；}//A[K, K] := (A[I, K] - TmpNum) / A[K, K];A[I, K] := (A[I, K] - TmpNum) / A[K, K];end;end;// 求解下三角方程组 G * y = b；SetLength(X, aN);X[0] := b[0] / A[0, 0];for I := 1 to aN-1 dobeginTmpNum := 0;for J := 0 to i-1 doTmpNum := TmpNum + A[I, J] * X[J];X[I] := (b[I] - TmpNum) / A[I, I];end;// 求解上三角方程组  G^T * x  = y;X[aN-1] := X[aN-1] / A[aN-1, aN-1];for I := aN-2 downto 0 dobeginTmpNum := 0;for J := I + 1 to aN - 1 doTmpNum := TmpNum + A[J, I] * X[J];X[I] := (X[I] - TmpNum) / A[I, I];end;Result := X;
end;end.

源码和文档（因为比较小，可执行程序也包含在内，还可以运行，但发现结果不正确。时过境迁，原因不明），放在百度网盘，供参考。

下载地址：https://pan.baidu.com/s/1DsuzGnCx1zqyIqrCrgugNg

提取码：3az1

多项式最小二乘拟合算法实现相关推荐

【多项式最小二乘拟合实验】
一.实验目的编制以函数{xk}k=0,-,n:为基的多项式最小二乘拟合程序.了解在matlab环境下曲线拟合问题的思想,掌握最小二乘拟合多项式拟合法和曲线拟合法,区分插值于拟合的不同之处. 二.实验 ...
Zernike多项式波前重建算法
首先,干涉图需要进行图像处理.通常通过干涉仪所获得的干涉图是光强分布,并非干涉图像的相位分布. 采用傅立叶变换空间载波法从强度分布得到相位分布,即从空间上呈正弦分布的光强信息,恢复出波面的相位信息,流 ...
最小二乘拟合n阶多项式【Matlab】
最小二乘拟合问题的求解附录4 最小二乘拟合n阶多项式程序Ⅰ 附录5 最小二乘拟合n阶多项式程序Ⅱ 附录5所示程序运行得到的拟合曲线等结果,与运行附录4所示程序结果相同,也相互验证了附录4.5所示 ...
最小二乘拟合的三种算法推导及Python代码
目录 1 最小二乘拟合(方法一) 1.1 数学推导 1.2 算例 1.3 Python 代码 2.最小二乘拟合(方法二) 2.1 数学推导 2.2 算例 2.3 Python 代码 3 最小二乘法拟合 ...
基于粒子群算法与最小二乘拟合函数参数
前言今天更新较晚主要还是学业繁忙,学习素材也不是很好找,可能很多同学们都在做数学建模以及应用统计时都会涉及到函数参数拟合的问题,一般最常用的方法是最小二乘法,但是当函数参数很多时,往往去普通最小二乘 ...
直线拟合算法（续：加权最小二乘）
直线拟合算法(续:加权最小二乘) 在此之前,我写过两篇文章介绍直线拟合算法: https://blog.csdn.net/liyuanbhu/article/details/50866802 http ...
多项式、正交多项式最小二乘拟合
最小二乘法求解矛盾方程组矛盾方程组:方程个数多于未知数个数,不能得到精确解析解. 使用最小二乘拟合得近似解. 误差函数: L=∑i=1n[∑j=1maijxj−bi]2L = \sum_{i = 1 ...
数学建模拟合(最小二乘拟合，多项式拟合，自定义函数拟合)
文章目录 matlab拟合工具箱最小二乘拟合理论推导用最小二乘法求解线性回归的k,b 怎么评价拟合的精度一个例子另一个例子,薄膜渗透率题目,最小二乘拟合溶液浓度变化多项式拟合自定义函数拟 ...
数学建模4 拟合算法
1.插值与拟合插值算法中,得到的多项式f(x)要经过所有样本点.但是如果样本点太多,那么这个多项式次数过高,会造成龙格现象.尽管我们可以选择分段的方法避免这种现象,但是更多时候我们倾向于得到一个确定 ...

多项式最小二乘拟合算法实现

1 多项式最小二乘拟合

2 Cholesky分解法解方程组

3 数值试验

参考文献

源程序：（基于 Delphi 实现）

多项式最小二乘拟合

Cholesky分解法解线性方程组

多项式最小二乘拟合算法实现相关推荐

最新文章

热门文章