探索最优解的典型习题解析？如何解答这些优化相关的问题？

如何进行下列最优化问题的求解

请查阅数值分析

最小二乘法

在我们探讨两个变量(x, y)之间的相互作用时，通常可以获得一系列配对的数据(x1, y1、x2, y2... xm, ym)；将这些数据绘制在x-y直角坐标系中（如图1），若发现这些点靠近一条直线，可以设定这条直线方程为（式1-1）。

Y计= a0+ a1 X（式1-1）

其中：a0、a1是任意实数

为建立这条直线方程，需要确定a0和a1，应用《最小二乘法原理》，将实测值Yi与利用（式1-1）计算值(Y计=a0+a1X)的偏差(Yi-Y计)的平方和〔∑(Yi- Y计)2〕最小化作为“优化标准”。

令:φ=∑(Yi- Y计)2（式1-2）

将（式1-1）代入（式1-2）中得:

φ=∑(Yi- a0- a1 Xi)2（式1-3）

当∑(Yi-Y计)平方最小时，可用函数φ对a0、a1求偏导数，令这两个偏导数等于零。

（式1-4）

（式1-5）

亦即：

m a0+(∑Xi) a1=∑Yi（式1-6）

(∑Xi) a0+(∑Xi2) a1=∑(Xi, Yi)（式1-7）

得到的两个关于a0、 a1为未知数的两个方程组，解这两个方程组得出：

a0=(∑Yi)/ m- a1(∑Xi)/ m（式1-8）

a1= [∑Xi Yi-(∑Xi∑Yi)/ m]/ [∑Xi2-(∑Xi)2/ m)]（式1-9）

此时将a0、a1代入（式1-1）中，此时的（式1-1）就是我们回归的元线性方程即：数学模型。

在回归过程中，回归的关联式不可能完全通过每个回归数据点(x1, y1、 x2, y2...xm,ym)，为了判断关联式的好坏，可借助相关系数“R”，统计量“F”，剩余标准偏差“S”进行判断；“R”越接近于 1越好；“F”的绝对值越大越好；“S”越接近于 0越好。

R= [∑XiYi- m(∑Xi/ m)(∑Yi/ m)]/ SQR{[∑Xi2- m(∑Xi/ m)2][∑Yi2- m(∑Yi/ m)2]}（式1-10）＊

在（式1-1）中，m为样本容量，即实验次数；Xi、Yi分别任意一组实验X、Y的数值。微积分应用课题一最小二乘法

从前面的学习中，我们知道最小二乘法可以用来处理一组数据，可以从一组测定的数据中寻找变量之间的依赖关系，这种函数关系称为经验公式。本课题将介绍最小二乘法的精确定义及如何寻找与之间近似成线性关系时的经验公式。假定实验测得变量之间的n个数据，x1, y1, x2, y2, ..., xn, yn，则在平面上，可以得到n个点，这种图形称为“散点图”，从图中可以大致看出这些点大致散落在某直线近旁，我们认为x与y之间近似为一线性函数，下面介绍求解步骤。

考虑函数，其中a0和a1是待定常数。如果在一直线上，可以认为变量之间的关系为y=a0+a1x。但一般说来，这些点不可能在同一直线上。记，它反映了用直线来描述y时，计算值与实际值产生的偏差。当然要求偏差越小越好，但由于偏差可正可负，因此不能认为总偏差为0时，函数就很好地反映了变量之间的关系，因为此时每个偏差的绝对值可能很大。为了改进这一缺陷，就考虑用来代替。但是由于绝对值不易作解析运算，因此，进一步用来度量总偏差。因偏差的平方和最小可以保证每个偏差都不会很大。于是问题归结为确定中的常数a0和a1，使为最小。用这种方法确定系数的方法称为最小二乘法。

由极值原理得，即

解此联立方程得

()

（3）根据公式()，利用“Apply”函数及集合的有关运算编写一个小的程序，求经验公式；（程序编写思路为：任意给定两个集合A（此处表示温度）、B（此处表示得率），由公式()可定义两个二元函数（集合A和B为其变量）分别表示a0和a1。集合A元素求和：Apply[Plus,A]表示将加法施加到集合A上，即各元素相加，例如Apply[Plus,{1,2,3}]=6；Length[A]表示集合A元素的个数，即为n；A.B表示两集合元素相乘相加；AB表示集合A与B元素对应相乘得到的新的集合。）

（4）在同一张图中显示直线及散点图；

（5）估计温度为200时产品得率。

然而，不少实际问题的观测数据，x1, y1, ..., xn, yn的散点图明显地不能用线性关系来描述，但确实散落在某一曲线近旁，这时可以根据散点图的轮廓和实际经验，选一条曲线来近似表达x与y的相互关系。

问题 II下表是美国旧轿车价格的调查资料，今以表示轿车的使用年数，p表示相应的平均价格，求p与x之间的关系。

使用年数

1 2 3 4 5 6 7 8 9 10

平均价格

2651 1943 1494 1087 765 538 484 290 226 204

（1）利用“ListPlot”函数绘出数据的散点图，注意观察有何特征？

（2）令p，绘出数据的散点图，注意观察有何特征？

（3）利用“Line”函数，将散点连接起来，说明有何特征？

（4）利用最小二乘法，求p与x之间的关系；

（5）求p与x之间的关系；

（6）在同一张图中显示散点图及关于x的图形。

思考与练习

1.假设一组数据：x1, y1, ..., xn, yn，变量之间近似成线性关系，试利用集合的有关运算，编写一简单程序：对于任意给定的数据集合，通过求解极值原理所包含的方程组，不需要给出a0、a1的计算表达式，立即得到a0、a1的值，并就本课题 I/(3)进行实验。

注：利用Transpose函数可以得到数据A的第一个分量的集合，命令格式为：

先求A的转置，然后取第一行元素，即为数据A的第一个分量集合，例如

(A即为矩阵)

(A即为数据A的第一个分量集合)

(A即为数据A的第二个分量集合)

B-C表示集合B与C对应元素相减所得的集合，如B-C。

2.最小二乘法在数学上称为曲线拟合，请使用拟合函数“Fit”重新计算a0、a1的值，并与先前的结果作一比较。

1. 最小二乘法在数学上称作曲线拟合，请采用拟合函数“适配”重新计算与的值，并与先前结果进行对比。

注：适配函数使用格式：

设变量为x，对数据A进行线性适配，例如对题1中的A适配函数为：

最小二乘优化反演方法

1. 对称四极装置电阻率测深一维线性反演的基本思路

· 第1步：构建目标函数

电阻率测深数据的反演方法可以归结为使以下目标函数趋于最小：

地球物理数据处理基础

这里，NS是供电极距数；ρai是第i个极距的实测视电阻率；M（Mj，j＝1，2…，NM）是模型参数（地层的电阻率或厚度），NM是预测模型参数个数；ρci是由预测模型正演计算所得的第i个极距的理论视电阻率；α是范数，当α＝2时即为最小二乘法。

由于地球物理反问题存在多解性，因此在求解时，为减少反演的多解性，常在目标函数式（9－48）中对模型加入先验信息，如加入模型先验信息的目标函数变为

地球物理数据处理基础

式中：Mb为背景模型；C为平滑度矩阵；λ为拉格朗日系数。即要求预测模型正演结果与实测数据最接近，而且要求预测模型与给定的背景模型最接近并且平滑。显然拉格朗日

系数λ可以调节对预测模型的要求偏重于哪边。

· 第2步：线性化

对式（9－48）利用泰勒展开，将模型的解在初始模型M处展开，并忽略二次及以上高次项，得

地球物理数据处理基础

要上式趋于最小，则对于各供电极距i（i＝1，2，…，NS）要满足下面的线性方程：

地球物理数据处理基础

· 第3步：计算视电阻率对模型参数的偏导数得到偏导数矩阵。

· 第4步：解线性方程组（9－50），得到预测模型M的修正量ΔM，建立新的预测模型。

· 第5步：计算新模型的正演理论视电阻率，与实测视电阻率进行对比，如果精度满足要求，则新的预测模型即为反演结果；否则回到第3步重新计算偏导数矩阵，重新计算模型修正量，直到精度满足要求。

上述反演过程中有三个问题需要解决：①偏导数的求解；②模型参数的处理，模型参数中有不同量纲的层电阻率和厚度，尤其电阻率参数变化范围很大，如果直接由上述方法求解，不但会导致方程（9－50）严重病态，而且电阻率和厚度参数的修正也不会正确，从而会导致反演方法不收敛；③线性方程组（9－50）的求解，线性方程组（9－50）往往是不相容方程组，求解方法不同，收敛速度和收敛性也不同。下面将对这三个问题的处理方法逐一进行介绍。

2. 偏导数的计算

对称四极装置电阻率测深的视电阻率对模型参数的偏导可以从数值滤波正演过程中直接给出。但一般反演中，偏导数通常用差分方法来计算。这里介绍用差分法来求对称四极装置电阻率测深的视电阻率对模型参数的偏导

取ΔMj＝0.1Mj，则

地球物理数据处理基础

3. 模型参数的无量纲化处理

为了解决视电阻率数据和电阻率模型参数变化范围大的问题，不少曾用它们的对数值来进行计算，但是这并没有解决层厚度参数与电阻率参数量纲不同的问题，因此，在反演过程中，出现了视电阻率测深曲线拟合很好，但预测模型与真实模型相差很大的情况，所以必须对方程（9－50）无量纲化。

将式（9－51）的偏导数代入方程（9－50）中可得

地球物理数据处理基础

其中：i＝1，2，…，NS。

上式两端同时除以ρci（M1，M2，…，Mj，…，MNM）可得

地球物理数据处理基础

这样方程组（9－52）可写为矩阵形式

地球物理数据处理基础

上式方程组中，左端系数矩阵A中各系数Aij、未知向量x中各变量xj以及右端向量B中各系数Bi都无量纲，从而解决了参数变量量纲不同的问题。

解线性方程组（9－54），则新模型参数Μ｛Mj，j＝1，2，…，NM｝为

地球物理数据处理基础

另外为了防止模型参数修改过量，实际过程中又作如下规定：xj>1.5时取xj＝1.5，xj<－0.5时取xj＝－0.5，即每次修改量不超过原有模型参数值的一半，以保证收敛稳定。

4. 线性方程组求解

线性方程组（9－54）的求解问题，由于A为NS×NM阶矩阵，且对地球物理反演问题而言A多半是接近奇异的，所以不能用常规方法求解。解这个非方阵线性方程组最常用的方法有：

（1）奇异值分解法：它的基本思想是建立在奇异值分解定理上，即任意NS×NM阶矩阵A均可分解为A＝UWVT，这里U为NS×NS阶正交阵和V为NM×NM阶正交阵，

地球物理数据处理基础

其中：δ1，…，δr为矩阵A的奇异值，r是矩阵的秩。当A非奇异时，奇异值较大，方程组（9－54）有广义逆解X＝VW－1UT，这里

地球物理数据处理基础

当A接近奇异时，有的奇异值就较小，此时由于W－1中系数过大，上述解的误差就较大。为了解决这个问题，维根斯（Wiggins）提出用最接近的矩阵R来代替A，R＝UWeVT，而

地球物理数据处理基础

W中小的奇异值在这里便被零代替了。因此有较精确的广义逆解X＝VW－1eUT。

（2）最小二乘法：在方程组（9－54）两端左乘AT，则方程变为ATAx＝ATB，当A非奇异时，有唯一解X＝（ATA）－1ATB；当A接近奇异时，用上式求解误差就会很大，导致反演不收敛。为了解决这个问题，马奎特提出给系数矩阵ATA的对角元素上加上一个常数，即ATAx＋λI＝ATB，以改善方程的条件，这就是著名的马奎特法又称阻尼最小二乘法。

分析奇异值（维根斯）分解方法和阻尼最小二乘（马奎特）方法，可以认为奇异值分解（维根斯）方法是将导致小的奇异值的方程取消来改善奇异值，而马奎特方法则是通过在ATA的对角元素上加常数来增大A的奇异值，两者的思想刚好相反。