猪冰龙

导航

1stopt、matlab和python用morris、sobol方法实现参数敏感性分析（sensitivity analysis）

首先看抛物线函数：

 

 

 现在我取a=2,b=3,c=4 ,得到如下函数：

 

 

x或t都是指自变量，就不改了，一个意思。 

 

问题是，我想知道对于此数据和模型，参数a,b,c的敏感性，也就是y的改变量与a、b、c的改变量的比值关系。

 

  • 首先用1stopt分析：

 

 1 NewCodeBlock"SA";
 2 Parameter a=[1,3],b=[2,4],c=[3,5];//要优化的参数及其范围
 3 Variable t,y;//变量
 4 Function y=a*t^2+b*t+c;
 5 Data;
 6 //t    y 变量顺序要和Variable后变量对应
 7 -5    39
 8 -4    24
 9 -3    13
10 -2    6
11 -1    3
12 0    4
13 1    9
14 2    18
15 3    31
16 4    48
17 5    69

 

 

 

 

 以下分别为各种参数敏感性方法（包括morris局部和全局，偏微分方法局部和全局）：

 

 

 

 

 

 

 

 

 

 

 

 以上就是4中方法的结果。用的目标函数是SSR,点的范围我用的是上下浮动50%，正好布满整个给定空间：

可以看morris局部敏感度分析具体数据：

a：

 

 

 

 

 

 

b：

 

 

 

 

 

c：

 

 

 

然后将灵敏度指数平均得到morris局部方法（本质是OAT方法，一次改变一个）分析结果：

 

 

 

 

全局方法（本质是通过拉丁抽样实现同时考虑多个参数的影响）就不是这个思路了，看一下超拉丁抽样：

 

 

 

 

 

 

 

总结一下1stopt中4种方法参数灵敏度结果：

 

morris方法局部：

 

 

名称 灵敏度指数 灵敏度指数(%)
a 3.916E18 89.1113892365457
b 4.125E17 9.38673341677096
c 6.6E16 1.50187734668335

morris方法全局：

名称 灵敏度指数 灵敏度指数(%)
a 2.00290323137447E18 81.4890482414753
b 2.10289523274038E17 8.55572692595605
c 2.44687506069835E17 9.95522483256862

 

 

偏微分方法局部：

名称 灵敏度指数 灵敏度指数(%)
a 3.99432000000567E18 89.1113892365577
b 4.20750000000213E17 9.38673341676365
c 6.73199999998753E16 1.50187734667864

 

 

偏微分方法全局：

 

 

名称 灵敏度指数 灵敏度指数(%)
a 3.98876054886661E18 82.084202339751
b 4.20542269921853E17 8.65428655186941
c 4.50049450186404E17 9.26151110837963

 

  • 下面用matlab写sobol方法进行分析：

 

  1 % sobol 参数敏感性分析
  2 %算法参考：
  3 % csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409
  4 % wiki： https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis   具体参考文献这里也有
  5 %运行环境 matlab2016b
  6 %作者 blzhu@buaa.edu.cn 2020年6月7日
  7 %% 初始化
  8 clc;
  9 clear all;
 10 close all;
 11 %% 设定：给定参数个数和各个参数的范围
 12 
 13 % % 1-自定义子函数1
 14 % D=3;% 维度3,几个参数
 15 % nPop=4:500:5000;% 采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢
 16 % VarMin=[-pi -pi -pi ];%各个参数下限  SALib ：S1: [ 0.30797549  0.44776661 -0.00425452] ；ST: [0.56013728 0.4387225  0.24284474]
 17 % VarMax=[pi pi pi];%各个参数上限
 18 % myFunction=@(x) Ishigami(x);%目标函数，也可以是个黑盒子
 19 
 20 % % 1-自定义子函数2
 21 D=3;% 维度3,几个参数
 22 nPop=4:50:1000;% 采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢
 23 VarMin=[1 2 3 ];%各个参数下限
 24 VarMax=[3 4 5];%各个参数上限
 25 myFunction=@(x) myx2(x);
 26 
 27 % % 1-自定义子函数3
 28 % D=4;% 维度3,几个参数
 29 % nPop=4:50:1000;% 采样点个数,也就是参数水平数 ,取大了好，比如4000，但慢
 30 % VarMin=[1 2 3 4];%各个参数下限
 31 % VarMax=[3 4 5 5];%各个参数上限
 32 % myFunction=@(x) myx3(x);
 33 
 34 
 35 %% 开始计算
 36 numnPop=numel(nPop);
 37 SAll=zeros(numnPop,D+1);%分别是各参数的敏感度，最后一列是各参数敏感度之和
 38 STAll=zeros(numnPop,D+1);
 39 for i=1:numnPop
 40     [S,ST]=sobol(D,nPop(i),VarMin,VarMax,myFunction);
 41     SAll(i,1:D)=S';
 42     SAll(i,D+1)=sum(SAll(i,1:D));
 43     STAll(i,1:D)=ST';
 44     STAll(i,D+1)=sum(STAll(i,1:D));
 45 end
 46 %% 绘图
 47 color=[1 0 0;0 1 0;0 0 1;0.5 0.1 0;0 0.3 0.4;0.6 0.7 0.2;0.5 0.8 0.9;0 0.2 0.1;0.1 0.5 0;0.1 0 0.5;0.5 0 0.1];%12种颜色 一般颜色不一样
 48 marker=['o','+','*','.','x','s','d','^','v','>','<','p','h'];%13种标记 一般标记也不一样; 字符数组，每个字符占1个位置
 49 linestyle=[string('-'),string('--'),string(':'),string('-.'),string('none')];
 50 useL=1;
 51 
 52 figure(1)
 53 for i=1:D+1
 54 plot(nPop,SAll(:,i),'Marker',marker(i),'LineStyle',char(linestyle(useL)),'Color',color(i,:),'LineWidth',1);
 55 hold on
 56 end
 57 title('Sobol-S');
 58 whichPara=sprintfc('%g',repmat(1:D+1,1,2));%把数字数组转化成字符串数组
 59 legend(whichPara,'Location','bestoutside');%加图例
 60 
 61 
 62 figure(2)
 63 for i=1:D+1
 64 plot(nPop,STAll(:,i),'Marker',marker(i),'LineStyle',char(linestyle(useL)),'Color',color(i,:),'LineWidth',1);
 65 hold on
 66 end
 67 title('Sobol-ST');
 68 whichPara=sprintfc('%g',repmat(1:D+1,1,2));%把数字数组转化成字符串数组
 69 legend(whichPara,'Location','bestoutside');%加图例
 70 
 71 disp('一阶影响指数（左方向收敛于1）Sobol-S：');
 72 disp(S);
 73 disp('总效应指数（大于等于1，且仅当myfun是纯相加时取等号）Sobol-ST：');
 74 disp(ST);
 75 disp(datetime);
 76 disp('parameter sensitive analyse success use sobol method!');
 77 %% 火车声音提示已经算完了
 78 load train
 79 sound(y,Fs)
 80 
 81 
 82 
 83 %% -------------------------子函数 matlab2016之前不支持子函数写在同一个m文档中----------------------------
 84 % 1-自定义子函数1（3个参数）Ishigami  https://www.sfu.ca/~ssurjano/ishigami.html
 85 function y=Ishigami(x)
 86 y=sin(x(1))+7*(sin(x(2)))^2+0.1*x(3)^4*sin(x(1));% SALib用的这个 
 87 % y=sin(x(1))+7*(sin(x(2)))^2+0.05*x(3)^4*sin(x(1));
 88 end
 89 
 90 % 1-自定义子函数2 （3个参数）
 91 function y=myx2(x)
 92 t=-5:1:5;%与此处有t范围和步距有关系
 93 % t=-5:0.1:5;%与此处有t范围和步距有关系
 94 ylab=2*t.^2+3*t+4;
 95 ytheory=x(1)*t.^2+x(2)*t+x(3);
 96 y=sum((ylab-ytheory).^2);%残差平方和SSR作为目标函数
 97 % y=sum((ylab-ytheory).^2)/numel(t);%各参数灵敏度与上式相同
 98 end
 99 
100 
101 % 1-自定义子函数3（4个参数）
102 function y=myx3(x)
103 t=-5:1:5;
104 ylab=2*t.^3+3*t.^2+4*t+5;
105 ytheory=x(1)*t.^3+x(2)*t.^2+x(3)*t+x(4);
106 y=sum((ylab-ytheory).^2);
107 end
108 
109 
110 
111 %% 2-求sobol敏感度
112 function [S,ST]=sobol(D,nPop,VarMin,VarMax,myFunction)
113 M=D*2;% 
114 %% 产生所需的各水平参数
115 VarMin=[VarMin,VarMin];
116 VarMax=[VarMax,VarMax];
117 p= sobolset(M);% https://www.cnblogs.com/zhubinglong/p/12260292.html
118 % R=p(1:nPop,:);% 我只用前nPop个
119 R=[];
120 for i=1:nPop
121     r=p(i,:);
122     r=VarMin+r.*(VarMax-VarMin);
123     R=[R; r];
124 end
125 % plot(R(:,1),'b*')
126 % 拆分为A B
127 A=R(:,1:D);% 每行代表一组参数，其中每列代表每组参数的一个参数；行数就代表共有几组参数
128 B=R(:,D+1:end);
129 % 根据A B 产生矩阵AB
130 AB=zeros(nPop,D,D);
131 for i=1:D
132     tempA=A;
133     tempA(:,i)=B(:,i);
134     AB(1:nPop,1:D,i)=tempA;
135 end
136 %% 求各参数解
137 YA=zeros(nPop,1);% 解
138 YB=zeros(nPop,1);
139 YAB=zeros(nPop,D);%分别代表YAB1,YAB2,YAB3，YAB(:,D)就代表YABD
140 for i=1:nPop
141     YA(i)=myFunction(A(i,:));
142     YB(i)=myFunction(B(i,:));
143     for j=1:D
144         YAB(i,j)=myFunction(AB(i,:,j));
145     end
146 end
147 %%  根据一阶影响指数公式：
148 VarX=zeros(D,1);% S的分子
149 S=zeros(D,1);
150 
151 % 0: 估算基于给定样本的方差(EXCEL var.p) ;   1:计算基于给定的样本总体的方差(EXCEL var.p())
152 % var([2.091363878    1.110366059    3.507651769    1.310950363    2.091363878    3.507651769    1.110366059    1.7066512],1);
153 VarY=var([YA;YB],1,'omitnan');% S的分母。 计算基于给定的样本总体的方差(EXCEL var.p())
154 for i=1:D
155     for j=1:nPop
156          VarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j));
157     end
158      VarX(i)=1/nPop*VarX(i);% 蒙特卡罗估计量
159      S(i)=VarX(i)/VarY;
160 end
161 
162 %% 总效应指数
163 EX=zeros(D,1);
164 ST=zeros(D,1);
165 for i=1:D
166     for j=1:nPop
167          EX(i)=EX(i)+(YA(j)-YAB(j,i))^2;
168     end
169       EX(i)=1/(2*nPop)* EX(i);% 蒙特卡罗估计量
170      ST(i)=EX(i)/VarY;
171 end
172 end

 

 

我分别取了不同个数的样点 4:50:1000 ，结果如下，可见1000个样点基本稳定了。

各参数的灵敏度：

一阶影响指数（左方向收敛于1）Sobol-S：
0.9728
0.0030
0.0001

总效应指数（大于等于1，且仅当myfun是纯相加时取等号）Sobol-ST：
0.9860
0.0031
0.0155

当然，也可以在matlab的fileexchange上下载各种工具箱，但这个根据csdn和wiki上写的算法相对简单些，便于魔改。

  • 用python的SALib包分析

python中也有参数灵敏度分析的包SALib（https://salib.readthedocs.io/en/latest/index.html），这个包支持以下算法：

下面以sobol方法举例：

 

 1 #https://salib.readthedocs.io/en/latest/basics.html#run-model
 2 # -*- coding: utf-8 -*-
 3 from SALib.sample import saltelli
 4 from SALib.analyze import sobol
 5 from SALib.test_functions import Ishigami
 6 import numpy as np
 7 import math
 8 from SALib.plotting.bar import plot as barplot
 9 import matplotlib.pyplot as plot
10 
11 # problem = {
12 #     'num_vars': 3,
13 #     'names': ['x1', 'x2', 'x3'],
14 #     'bounds': [[-3.14159265359, 3.14159265359],
15 #                [-3.14159265359, 3.14159265359],
16 #                [-3.14159265359, 3.14159265359]]
17 # }
18 
19 problem = {
20     'num_vars': 3,
21     'names': ['x1', 'x2', 'x3'],
22     'bounds': [[1, 3],
23                [2, 4],
24                [3, 5]]
25 }
26 
27 
28 param_values = saltelli.sample(problem, 1024)# 不管用哪个方法计算y,这个要有.个数建议为2^n
29 np.savetxt("param_values.txt", param_values)# 将参数变化保存，其实是各参数范围内的sobol抽样
30 
31 ##  计算Y
32 # #1-自定义-1
33 # Y = np.zeros([param_values.shape[0]])
34 # A = 7
35 # B = 0.1
36 # for i, X in enumerate(param_values):
37 #     Y[i] = math.sin(X[0]) + A * math.pow(math.sin(X[1]), 2) + \
38 #             B * math.pow(X[2], 4) * math.sin(X[0])
39 
40 #1-自定义-2
41 Y = np.zeros([param_values.shape[0]])
42 for i, X in enumerate(param_values):
43     tarr=np.arange(-5,6,1);
44     yerror=0.0;
45     for t in tarr:
46         ylab=2*t**2+3*t+4;
47         ytheory=X[0]*t**2+X[1]*t+X[2];
48         yerror=yerror+(ylab-ytheory)**2;
49 
50     Y[i] = math.sqrt(yerror);
51 
52 # 2-load计算好的txt
53 # Y = np.loadtxt("outputs.txt", float)
54 
55 # 3-SALib自带测试函数
56 # Y = Ishigami.evaluate(param_values)
57 ## np.savetxt("outputs.txt", Y)#将因变量变化结果保存
58 
59 Si = sobol.analyze(problem, Y ,print_to_console=True)
60 print()#自动输出S1（单个参数对因变量的影响）、ST（考虑各个变量相互影响）和S2(两两参数之间影响)，需有,print_to_console=True
61 
62 print("all parameters first-order sensitivity indices   S1:")
63 print(Si['S1'])# 一阶影响指数
64 print("all parameters second-order sensitivity indices   S2:")
65 print(Si['S2'])#二阶影响指数
66 print("all parameters total  sensitivity indices   ST:")
67 print(Si['ST'])# 总效应指数
68 
69 #绘图 https://zhuanlan.zhihu.com/p/137953265
70 Si_df = Si.to_df()
71 barplot(Si_df[0])
72 plot.show()

 

 

 

输出结果：

Parameter S1 S1_conf ST ST_conf
x1 0.969397 0.069624 0.982232 0.058139
x2 0.007222 0.009055 0.009680 0.001325
x3 0.000848 0.008468 0.011699 0.001033

Parameter_1 Parameter_2 S2 S2_conf
x1 x2 0.000330 0.070070
x1 x3 0.010014 0.069389
x2 x3 -0.000129 0.013788

all parameters first-order sensitivity indices   S1:
[9.69397123e-01 7.22243327e-03 8.47690887e-04]
all parameters second-order sensitivity indices   S2:
[[        nan  0.00032978  0.01001386]
 [        nan         nan -0.00012935]
 [        nan         nan         nan]]
all parameters total  sensitivity indices   ST:
[0.9822323  0.00968001 0.01169928]
 
也可以用morris方法：
只需要导入 
from SALib.analyze import morris
然后用 
Si = morris.analyze(problem, Y ,print_to_console=True)  代替 
Si = sobol.analyze(problem, Y ,print_to_console=True)

 

 1 # https://salib.readthedocs.io/en/latest/basics.html#run-model
 2 # -*- coding: utf-8 -*-
 3 from SALib.sample import saltelli
 4 from SALib.analyze import sobol
 5 from SALib.analyze import morris
 6 from SALib.test_functions import Ishigami
 7 import numpy as np
 8 import math
 9 import SALib
10 from SALib.plotting.bar import plot as barplot
11 import matplotlib.pyplot as plt
12 
13 # problem = {
14 #     'num_vars': 3,
15 #     'names': ['x1', 'x2', 'x3'],
16 #     'bounds': [[-3.14159265359, 3.14159265359],
17 #                [-3.14159265359, 3.14159265359],
18 #                [-3.14159265359, 3.14159265359]]
19 # }
20 
21 problem = {
22     'num_vars': 3,
23     'names': ['x1', 'x2', 'x3'],
24     'bounds': [[1, 3],
25                [2, 4],
26                [3, 5]]
27 }
28 
29 
30 param_values = saltelli.sample(problem, 1024)# 产生参数样本； 不管用哪个方法计算y,这个要有
31 np.savetxt("param_values.txt", param_values)# 将参数变化保存，其实是各参数范围内的sobol抽样
32 
33 ##  计算Y
34 # #1-自定义-1
35 # Y = np.zeros([param_values.shape[0]])
36 # A = 7
37 # B = 0.1
38 # for i, X in enumerate(param_values):
39 #     Y[i] = math.sin(X[0]) + A * math.pow(math.sin(X[1]), 2) + \
40 #             B * math.pow(X[2], 4) * math.sin(X[0])
41 
42 #1-自定义-2
43 Y = np.zeros([param_values.shape[0]])
44 for i, X in enumerate(param_values):
45     tarr=np.arange(-5,6,1);
46     yerror=0.0;
47     for t in tarr:
48         ylab=2*t**2+3*t+4;
49         ytheory=X[0]*t**2+X[1]*t+X[2];
50         yerror=yerror+(ylab-ytheory)**2;
51 
52     Y[i] = math.sqrt(yerror);
53 
54 # 2-load计算好的txt
55 # Y = np.loadtxt("outputs.txt", float)
56 
57 # 3-SALib自带测试函数
58 # Y = Ishigami.evaluate(param_values)
59 ## np.savetxt("outputs.txt", Y)#将因变量变化结果保存
60 
61 # Si = sobol.analyze(problem, Y ,print_to_console=True)
62 
63 X = np.loadtxt("param_values.txt", float)# X = morris.sample(problem, 1024)
64 Si = morris.analyze(problem,X,Y,print_to_console=True)
65 print()#自动输出S1（单个参数对因变量的影响）、ST（考虑各个变量相互影响）和S2(两两参数之间影响)，需有,print_to_console=True
66 # morris.analyze源代码：https://salib.readthedocs.io/en/latest/_modules/SALib/analyze/morris.html?highlight=morris.analyze#

 

结果：

mu mu_star sigma mu_star_conf
x1 -0.027094 31.955401 41.773095 1.105048
x2 0.007385 35.174581 39.713491 0.719772
x3 0.031841 18.773557 27.977986 0.828840


 
小技巧：
SALib如果不便于将目标函数写为函数的形式，可以将代码：
np.savetxt("param_values.txt", param_values)# 将参数变化保存，其实是各参数范围内的sobol抽样
生成的抽样带入自己的系统，然后根据自己需要生成对应抽样的目标函数，将目标函数放入同目录下的 outputs.txt 文档中，一行一个结果，然后用这个语句代替上面求Y[i]:
Y = np.loadtxt("outputs.txt", float)

 

就是说提供了参数变化以及目标函数变化，用SALib就可以求参数灵敏度了。

 

 

 总结：

matlab我用的sobol生成的抽样和别人的不一样，不知道为什么，这个是造成与python计算不一样的一个原因吧。但差别不大。

灵敏度分析的概念在此处没有详细讲，可以参考：https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis

 

 

 

 

 

 

 

 

 

posted on 2020-06-09 22:30  猪冰龙  阅读(16981)  评论(2)  编辑  收藏  举报