第一次接触到1stOpt是因为N年前需要求解一组非常复杂的微分方程组，自己又懒得用matlab敲代码，于是就在网上搜索有没有更为轻松便捷的办法。一找就找到了1stOpt，虽然是个绿色和谐版，但基本上能满足自己需求，顺带把1stOpt的帮助文档都看了个遍。

看了1stOpt官网的介绍，据说1stOpt是一款不依赖初始值的数学优化综合工具，这一点连老牌的Matlab、Origin、Lingo、Gams等都是望洋兴叹的，这得益于其开发的通用全局优化算法。不依赖初始值对科学计算而言是极其有用的，比如求根、求解微分、优化等问题往往对初始值有高度依赖，初始值选择不当可能会得到高误差甚至是错误的结果，而1stOpt完美解决了此问题。现已被广泛用于工程、金融、决策、流通、生产等分析、计算、规划领域。揽括：模型自动优化率定；参数估算；任意模型公式线性，非线性拟合，回归；非线性连立方程组求解；常微方程（组）求解（初值、边值问题）；常微分方程（组）拟合求解；复数方程求解、复数非线性拟合；任意维函数，隐函数极值求解；隐函数根求解，作图，求极值；线性，非线性及整数规划；组合优化问题等。

下面就通过几个实例来看看1stOpt的实力：

特别说明：本次采用的1stOpt 5.5 英文试用版，实例全部来自于5.0官方帮助文档。

由于试用版有参数个数限制，计算参数不能超过4个，目前官方最新版为8.0，9.0版即将发布，下面就一起来看看吧。

1、二维针状全局最优的函数

该函数的函数图形如下所示，在(50,50)处取得全局最大值 1.1512，其第二极大值为1.12837，它是一个多峰值函数，采用传统优化方法几乎不能找到全局最优点。

1stOpt源代码：

Parameter x[0,100], y[0,100];
ConstStr r = sqrt((x-50)^2+(y-50)^2)+exp(1);
MaxFunction Sin(r)/r + 1;

2、隐函数优化

1stOpt源代码：

Parameter x[-1,7], y[-2,2];
Minimum = z;
Function z = sin((z*x-0.5)^2 + x*2*y^2-z/10)*exp(-((x-0.5-exp(-y+z))^2 + y^2-z/5+3));

3、非线性规划

1stOpt源代码：

Parameter x1[0,10], x2[0,10];
MaxFunction (sin(2*pi*x1))^3*sin(2*pi*x2)/(x1^3*(x1+x2));
x1^2-x2+1 <= 0;
1-x1+(x2-4)^2 <= 0;

4、整数规划

1stOpt源代码：

Parameter x(1:2)[0,100,0], x3[3,100,0], x4[0,100,0];
MinFunction -2*x1-x2-4*x3-3*x4;
2*x2+x3+4*x4<54;
3*x1+4*x2+5*x3-x4<62;
x1+2*x2+x3-2*x4<-3;

5、带等式约束的优化问题

1stOpt源代码：

Algorithm = DE1[100];
PenaltyFactor = 5;
ParameterDomain = [-50,50];
MinFunction exp(sin(50*x)) +sin(60*exp(y)) +sin(70*sin(x))
+sin(sin(80*y))-sin(10*(x+y))+(x^2+y^2)^sin(y)/4;
x-((cos(y))^x-x)^y<=0;

6、非线性拟合

美国国家标准与技术研究院（NIST）提供有一套 27 道非线性拟合测试题，世界上几乎所有著名的数据分析软件包都以能通过该套测试题集为验证标准。经对比测试，1stOpt 是目前唯一不依赖使用 NIST 提供的初始值， 而能以任意随机初始值就可求得全部最优解的软件包(如果使用 NIST 提供的初始值，则更可轻易求得最优解)。由于在实际应用当中，选择确定合理的初始值组是一件非常困难的事，尤其是在参数量比较多的情况下。从此意义而言1stOpt的实用能力达业界领先水平。 

示例代码：

Variables x, y;
Parameters a, b, c, d;
Function y = a-b*exp(-c*x^d);
Data;
0.05 0.13
0.15 0.13
0.25 0.19
0.35 0.34

7、求解带积分方程组

本题来源于小木虫，由于参数限制，改为4个参数：

http://muchong.com/t-4192342-1-authorid-8600892

Parameters a,b,c,d;
Function int(exp(a+b*x+c*x^2+d*x^3),x=0.5779,13.5883)=1;
int(x*exp(a+b*x+c*x^2+d*x^3),x=0.5779,13.5883)=6.86;
int(x^2*exp(a+b*x+c*x^2+d*x^3),x=0.5779,13.5883)=54.39;
int(x^3*exp(a+b*x+c*x^2+d*x^3),x=0.5779,13.5883)=478.14;

由于功能和参数限制，没法演示1stOpt更为强大的求解功能，暂且演示到此。

在少数铁杆粉丝的支持下，小编目前已拥有1stOpt 9.0最新版的正版版权。如有需要购买正版1stOpt的伙伴，加公众号任意Q联系群主，可能能为大家争取到专属优惠哦！

参考资料：www.7d-soft.com

图片来源：www.7d-soft.com

如需转载，请在公众号中回复“转载”获取授权，如未经授权擅自搬运抄袭的，本公众号将保留一切追责权利！

首先看抛物线函数：

现在我取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方法局部：

名称 灵敏度指数 灵敏度指数(%)

a3.916E1889.1113892365457

b4.125E179.38673341677096

c6.6E161.50187734668335

morris方法全局：

名称 灵敏度指数 灵敏度指数(%)

a2.00290323137447E1881.4890482414753

b2.10289523274038E178.55572692595605

c2.44687506069835E179.95522483256862

偏微分方法局部：

名称 灵敏度指数 灵敏度指数(%)

a3.99432000000567E1889.1113892365577

b4.20750000000213E179.38673341676365

c6.73199999998753E161.50187734667864

偏微分方法全局：

名称 灵敏度指数 灵敏度指数(%)

a3.98876054886661E1882.084202339751

b4.20542269921853E178.65428655186941

c4.50049450186404E179.26151110837963

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

1 %sobol 参数敏感性分析2 %算法参考：3 % csdn : https://blog.csdn.net/xiaosebi1111/article/details/46517409

4 % wiki： https://5 %运行环境 matlab2016b6 %作者 blzhu@buaa.edu.cn 2020年6月7日7 %%初始化8 clc;9 clear all;10 close all;11 %%设定：给定参数个数和各个参数的范围12

13 % % 1-自定义子函数114 % 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-自定义子函数221 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-自定义子函数328 % 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:numnPop40 [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 end46 %%绘图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','>','

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 on56 end57 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 on66 end67 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 train79 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 end89

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 end99

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 end108

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:nPop121 r=p(i,:);122 r=VarMin+r.*(VarMax-VarMin);123 R=[R; r];124 end125 % plot(R(:,1),'b*')126 %拆分为A B127 A=R(:,1:D);%每行代表一组参数，其中每列代表每组参数的一个参数；行数就代表共有几组参数128 B=R(:,D+1:end);129 %根据A B 产生矩阵AB130 AB=zeros(nPop,D,D);131 for i=1:D132 tempA=A;133 tempA(:,i)=B(:,i);134 AB(1:nPop,1:D,i)=tempA;135 end136 %%求各参数解137 YA=zeros(nPop,1);%解138 YB=zeros(nPop,1);139 YAB=zeros(nPop,D);%分别代表YAB1,YAB2,YAB3，YAB(:,D)就代表YABD140 for i=1:nPop141 YA(i)=myFunction(A(i,:));142 YB(i)=myFunction(B(i,:));143 for j=1:D144 YAB(i,j)=myFunction(AB(i,:,j));145 end146 end147 %%根据一阶影响指数公式：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:D155 for j=1:nPop156 VarX(i)=VarX(i)+YB(j)*(YAB(j,i)-YA(j));157 end158 VarX(i)=1/nPop*VarX(i);%蒙特卡罗估计量159 S(i)=VarX(i)/VarY;160 end161

162 %%总效应指数163 EX=zeros(D,1);164 ST=zeros(D,1);165 for i=1:D166 for j=1:nPop167 EX(i)=EX(i)+(YA(j)-YAB(j,i))^2;168 end169 EX(i)=1/(2*nPop)* EX(i);%蒙特卡罗估计量170 ST(i)=EX(i)/VarY;171 end172 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包分析

Method of Morris, including groups and optimal trajectories ([Morris 1991], [Campolongo et al. 2007])

Fourier Amplitude Sensitivity Test (FAST) ([Cukier et al. 1973], [Saltelli et al. 1999])

Delta Moment-Independent Measure ([Borgonovo 2007], [Plischke et al. 2013])

Derivative-based Global Sensitivity Measure (DGSM) ([Sobol and Kucherenko 2009])

Fractional Factorial Sensitivity Analysis ([Saltelli et al. 2008])

下面以sobol方法举例：

1 #https://salib.readthedocs.io/en/latest/basics.html#run-model

2 #-*- coding: utf-8 -*-

3 from SALib.sample importsaltelli4 from SALib.analyze importsobol5 from SALib.test_functions importIshigami6 importnumpy as np7 importmath8 from SALib.plotting.bar importplot as barplot9 importmatplotlib.pyplot as plot10

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, 1000)#不管用哪个方法计算y,这个要有

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 inenumerate(param_values):43 tarr=np.arange(-5,6,1);44 yerror=0.0;45 for t intarr: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 importsaltelli4 from SALib.analyze importsobol5 from SALib.analyze importmorris6 from SALib.test_functions importIshigami7 importnumpy as np8 importmath9 from SALib.plotting.bar importplot as barplot10 importmatplotlib.pyplot as plot11

12 #problem = {

13 #'num_vars': 3,

14 #'names': ['x1', 'x2', 'x3'],

15 #'bounds': [[-3.14159265359, 3.14159265359],

16 #[-3.14159265359, 3.14159265359],

17 #[-3.14159265359, 3.14159265359]]

18 #}

19

20 problem ={21 'num_vars': 3,22 'names': ['x1', 'x2', 'x3'],23 'bounds': [[1, 3],24 [2, 4],25 [3, 5]]26 }27

28

29 param_values = saltelli.sample(problem, 1000)#不管用哪个方法计算y,这个要有

30 np.savetxt("param_values.txt", param_values)#将参数变化保存，其实是各参数范围内的sobol抽样

31

32 ## 计算Y

33 ##1-自定义-1

34 #Y = np.zeros([param_values.shape[0]])

35 #A = 7

36 #B = 0.1

37 #for i, X in enumerate(param_values):

38 #Y[i] = math.sin(X[0]) + A * math.pow(math.sin(X[1]), 2) + \

39 #B * math.pow(X[2], 4) * math.sin(X[0])

40

41 #1-自定义-2

42 Y =np.zeros([param_values.shape[0]])43 for i, X inenumerate(param_values):44 tarr=np.arange(-5,6,1);45 yerror=0.0;46 for t intarr:47 ylab=2*t**2+3*t+4;48 ytheory=X[0]*t**2+X[1]*t+X[2];49 yerror=yerror+(ylab-ytheory)**2;50

51 Y[i] =math.sqrt(yerror);52

53 #2-load计算好的txt

54 #Y = np.loadtxt("outputs.txt", float)

55

56 #3-SALib自带测试函数

57 #Y = Ishigami.evaluate(param_values)

58 ## np.savetxt("outputs.txt", Y)#将因变量变化结果保存

59

60 #Si = sobol.analyze(problem, Y ,print_to_console=True)

61 Si = morris.analyze(problem, Y ,print_to_console=True)62 print()#自动输出S1(单个参数对因变量的影响)、ST(考虑各个变量相互影响)和S2(两两参数之间影响)，需有,print_to_console=True

63

64 print("all parameters first-order sensitivity indices S1:")65 print(Si['S1'])#一阶影响指数

66 print("all parameters second-order sensitivity indices S2:")67 print(Si['S2'])#二阶影响指数

68 print("all parameters total sensitivity indices ST:")69 print(Si['ST'])#总效应指数

70

71 #绘图 https://zhuanlan.zhihu.com/p/137953265

72 Si_df =Si.to_df()73 barplot(Si_df[0])74 plot.show()

python3.8.1+SALib1.3 use morris

Parameter S1 S1_conf ST ST_conf

x1 0.969397 0.072129 0.982232 0.062795

x2 0.007222 0.008033 0.009680 0.001303

x3 0.000848 0.009519 0.011699 0.001045

Parameter_1 Parameter_2 S2 S2_conf

x1 x2 0.000330 0.087924

x1 x3 0.010014 0.087743

x2 x3 -0.000129 0.013700

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]

小技巧：

SALib如果不便于将目标函数写为函数的形式，可以将代码：

np.savetxt("param_values.txt", param_values)# 将参数变化保存，其实是各参数范围内的sobol抽样

生成的抽样带入自己的系统，然后根据自己需要生成对应抽样的目标函数，将目标函数放入同目录下的 outputs.txt 文档中，一行一个结果，然后用这个语句代替上面求Y[i]:

Y = np.loadtxt("outputs.txt", float)

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

总结：

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