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

  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 %作者 [email protected] 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);%残差平方和s-s-r作为目标函数
 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

 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()