随机响应面法

随机响应面法

随机响应面法同神经网络法、Kriging法、响应面法、支持向量机等方法一样，都属于代理模型范畴。所谓代理模型，就是用于替代现实中难以掌握系统规律的一种模块的统称，比如，经济学中的供求关系与价格关系最为密切，但实际上仍有许多其他因素影响着供需关系，比如收入、城市级别等因素，这些因素共同影响着供需关系，假设X为需求，J为价格，S为收入，C为城市级别，则他们的关系为X=f(J,S,C)，此关系通常情况下是隐式关系，并不能用显式表达式表达。另外，对于边坡模型，内摩擦角、粘聚力、密度和渗透系数等都对安全系数有影响，但他们的关系是有限元计算软件这个“黑盒子”计算得出的，计算时间较长且不能直观看出各因素与安全系数结果之间的关系。故引入了随机响应面法等代理模型，用代理模型来代替这个“黑盒子”，毕竟多项式计算比有限元计算速度快的不是一点半点，但是代理模型也有很多不足之处，如果训练的不好，很可能会造成很大的误差。
本文将细致的讲解随机响应面法相关知识，并且利用随机matlab编写多项式的自动生成程序，以便随机变量和阶数较多情况下的使用。

响应面与随机响应面

不论是响应面法还是随机响应面法都属于数值拟合的范畴，不少朋友会把拟合与插值搞混，这里加以区分以下，插值是要求所有样本点都在所的函数曲线之上，换句话说，插值是将样本点用曲线“串”起来的过程。拟合则有些不同，它是指通过距离最小原则，即每个点到最终拟合好的函数之间的距离之和最小来表明拟合程度好坏的一种方式，可想而知，样本点不一定在拟合函数上。插值与拟合各有自己优势，感兴趣的朋友可以找本数值分析了解一下。

书归正传，响应面法和随机响应面法的区别有二，一、随机响应面所使用的多项式是Hermite随机多项式，而响应面法使用的是简单的多项式，且多用不含交叉项的二次多项式作为拟合多项式。二、随机响应面法更加适合拉丁超立方抽样作为其配点方法，而响应面法一般使用概率配点法。

响应面

说到这，大家一定觉得响应面法有些飘渺难以理解，没关系，只要从字面理解响应面三个字就可以了，“响应”是指系统对输入数据反馈，举个简单的例子，一个包含x,y的函数z=x+y，输入x=1,y=1，即可得到该系统的响应值z=2。下面理解“面”的概念，假设现在有一个二元函数，Z=X² +Y²，他在三维空间是什么样子的呢？我们通过Matlab来绘制一下把。

代码如下：

x=linspace(-2,2);
y=linspace(-2,2);
[x y]=meshgrid(x,y);
z=x.^2+y.^2;
surf(x,y,z)
shading interp

可见，绘制出的图像是三维空间的一个曲面，那么有人会问，假如这个函数是M=X²+Y²+Z²呢？它是不能在三维空间标识出来的，确实，所以此时我们称之为高维面，已经超出我们观察能力范围的面都视为高维面。既然都是响应系统，也都有面的性质，那么就能理解为什么叫“响应面”了。

响应面分析流程大致分为7步：

1. 选取含待定系数的响应面函数代替不能明确表达的实际功能函数，常用多项式功能函数。
2. 确定各个随机变量的概率分布形式和取样范围。
3. 根据功能函数和随机变量分布，选用一定数量的样本点。
4. 使用样本点数据，建立合适的模型（比如有限元计算）求解实际功能响应。使用N各样本点数据，可以得到系统在N个样本点处的响应值。
5. 将样本点及其功能响应数据代入响应面函数，建立线性方程组，求解待定系数，从而获得明确的响应面函数。
6. 基于响应面函数进行一些列操作
7. 如果想得到更加精确的响应面模型，则可以迭代计算，重新取样，求解新的响应面函数。

响应面法的存在有两大重要意义：
1.通过数据拟合出一个精度最高的响应面用来代替真实系统的响应关系。
2.寻找拟合后响应面山的极值点（寻优）

注意：
响应面过程中的自变量应限制在合理范围内，这样拟合后的响应面才符合真实情况

---

接下来以多项式响应面方法为例，介绍响应面拟合相关知识。
首先响应面是通过近似构造一个简单的显式函数来近似代替真实目标函数的过程，如果构造最高次为一次的多项式，就称为一次响应面法，如果构造的多项式最高次数为二次就成为二次响应面法，以此类推。
现假设有一复杂系统，它有许多因素影响，那么该系统的响应值与影响因素就可以表示为：

在这里插入图片描述
由上式可知，该系统有q个影响因素，现假设对每个影响因素进行N次取值，然后将N组影响因素分别带入系统求解，得到N个响应量，再通过一次响应面对上述系统进行拟合，则可以表示为下式：
式中x_ij为为第j个自变量第i次系统功能分析时的取值，z_i为第i次系统功能分析时的响应值λ_i为一次响应面函数的待定系数，ε_i是第i次分析时的误差项。上面提到的系统功能分析可以理解一次计算，比如我将边坡模型建立完成后，设定抗剪强度参数，进行一次计算，就是一次系统功能分析，当然系统的响应值就是这次分析后得到的安全系数。

为了提高计算精度，通常采用二次响应面法来进行拟合。二次响应面可以表示为：

当给定一组已知数据时，可通过上式得到含待定系数的线性方程组，通常，通常利用最小二乘法求解待定系数。待定系数确定后，就得到了响应面函数。
从上式还可以看出，，待定系数个数为2q+1+C²_q。当随机变量个数较多时，为确定各个待定系数需要取很多样本点进行分析，降低了该方法的效率，故很多情况下使用不含交叉项的二次多项式来拟合。此时待定系数个数为2q+1个。见下式：

随机响应面

为接近岩土工程实际，需要考虑岩土体应力-应变关系和渗流等作用的影响，采用有限元方法（或借助于有限元软件）分析边坡变形与稳定性，一般情况下所获得的边坡输出响应量（安全系数、位移、应力、孔隙水压力和渗透坡降等）是岩土体物理力学参数的非线性隐式函数。因此，在获得岩土体物理力学输入参数模型之后，为了便于可靠度分析，一般需要首先建立边坡输出响应量的代理模型，即建立输出响应量与输入参数之间的近似显示函数关系，再采用常规可靠度分析方法计算可靠度。本文采用Hermite随机多项式（PCE）建立输出输出响应量与输入参数间的近似显示函数关系。

Hermite随机多项式

根据Cameron-Martin定理，当输出响应量表示为正态分布随机变量时，Hermite随机多项式展开有关于阶数p的指数收敛速度。对于任意来自Hillbert空间的函数皆可由Hermite多项式函数来近似表示。
采用Hermite随机多项式展开建立边坡输出响应量F_p(ζ)与输入参数之间的近似显式表达式为：

式中a为待定系数;N为随机变量数目，也是参数随机场所离散的随机变量数目之和；ζ为独立标准正态随机变量。p为随机多项式展开阶数；Γ为自由度为p 的Hermite多项式：

Hermite多项式项数

Hermite随机多项式的项数是由其维度和阶次确定的，假设N_c为随机多项式展开的项数，N为维度，p为阶次，则随机多项式项数可由下面的公式求出：

Hermite随机多项式的推导

一维Hermite多项式

一维Hermite多项式既可以由微分方程求得，也可以使用递推公式求得：

1. 微分方程
   
2. 递推公式
   
   通过以上两种方法可以推导出一维Hermite多项式：
   

高维Hermite多项式

根据Hermite多项式的性质，N维Hermite随机多项式展开可表示为多个变量的一维Hermite多项式展开的乘积：

当然也可以通过微分方程得到多维高阶Hermite多项式了，我通过Matlab实现了四阶任意维度的Hermite多项式程序，可以实现四阶以下（包括四阶），任意维度的多项式生成，辅助matlab进行最小二乘回归。

随机多项式系数的计算

随机响应面法最终可以通过最小二乘法计算，计算出每一个多项式前的未知系数，从而得到显式功能函数。

当确定了输出响应量的随机多项式展开后，接着用随机配点的方法计算随机多项式展开的待定系数。其中会用到拉丁超立方抽样。将抽样得到的配点带入系统求出响应值，根据配点和响应值建立方程组

为准确求解 Hermite 随机多项式展开系数，般要求所需的配点数目 大于待定系数的数目。此外也可以基于下式进行最小二乘分析：

这种基于输入参数配点的随机多项式展开待定系数的求解，并建立显式表达功能函数的方法称为多项式展开法。

Pei_Number=140;
x1=[0.3790144	0.9450246	-0.183698	-0.5890345	-0.7252413	0.275119	-0.0286829	-0.4730099	0.4818152	1.5793378	-1.2062526	-2.6833556	0.2651452	-0.3656525	1.9877294	-2.2973643	-0.8909261	-0.6915144	0.07418	-0.2921757	0.2245247	0.1282181	-0.4385361	-1.5027866	-1.6420049	-0.1122994	-0.4632194	0.2369559	0.2002068	-0.0368794	0.0410539	0.6779487	0.9874637	-1.180399	0.5782359	-0.7572747	0.6114364	0.0683059	0.6980659	-0.9017496	1.7006964	-0.5224338	0.1965467	1.2848361	1.235287	1.8922828	-0.3197828	0.9543268	-0.699335	-2.0914089	0.172459	-1.0216256	-0.129145	0.106291	0.3583854	0.8909212	-0.423484	1.3407049	0.6311065	1.0543415	0.1124861	-1.6587895	-0.0903823	-0.2412594	-1.1038534	-0.1745446	0.806621	-0.6353967	-0.4979306	1.3981427	-1.5332693	0.7504411	-0.0833992	0.6713888	1.4386319	-0.3295978	-0.9381026	1.1111662	-1.2772358	-0.5349091	-0.2601284	-1.0842007	-0.769635	-1.7591731	2.0941481	2.3886093	0.8198773	-1.33246	1.489417	-0.2750108	0.5013085	-1.9596891	1.2773472	0.7806756	-0.5590123	0.9200276	0.859668	-0.9911279	3.1024445	-1.8382227	0.424901	1.724414	-0.8511768	-0.6191964	-1.038013	-0.2008088	-0.2260041	-0.3925501	0.0264446	-0.5701868	1.5407155	-0.6672858	-1.3166698	-1.1334171	1.1451734	0.0113391	0.5377827	0.4583812	-0.0151496	-0.1617397	0.1583899	1.0974508	0.3162461	0.5512921	-0.8330679	0.39712	-1.4586928	1.0200726	0.2985681