Question

哪位大师帮忙看看下面这个MATLAB程序咋回事？纠结了好久了。

Mu=[10, 2.5]; ％输入平均值。 %主函数部分

Sigama=[2, 0.375]; ％输入标准差。
X0=[10, 2.5]; ％输入初始迭代点。
A=[];b=[];Aeq=[];beq=[];
lb=[];ub=[];
[X,fval,exitflag,output]=fmincon(@bata2,X0,A,b,Aeq,beq,lb,ub,@st);
bata=sqrt(fval); ％计算可靠指标β。
Pf=cdf('norm',-bata,0.1); ％计算失效概率。

function CC=bata2(X)％目标函数子函数。
global Mu Sigama ％定义全局变量。
CC=((X(1)-Mu(1))/Sigama(1))^2+((X(2)-Mu(2))/Sigama(2))^2 ％X均为正态分布

function [c,ceq]=st(X) ％约束条件子函数。
c=[]; ％非线性不等式约束。
ceq=xym(X) ％非线性等式约束

function Z=xym(X) ％响应面函数值子函数。
A1=[1 10 2.5 100 6.25
1.0 16.0 2.5 256 6.25
1.0 4.0 2.5 16 6.25
1.0 10.0 3.625 100.0 13.1406
1.0 10.0 1.375 100.0 1.8906];
B1=[13.6728;10.8005;16.5451;16.8897;-10.31363];
C=A1\B1; %求解二次响应面函数各项系数。
Z=C(1)+C(2)*X(1)+C(3)*X(2)+C(4)*X(1)+C(5)*X(2); ％计算功能函数值

Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you have
specified. FMINCON will use the active-set algorithm instead. For information on applicable algorithms, see Choosing
the Algorithm in the documentation.
> In fmincon at 504
In rsm1 at 10
Warning: Your current settings will run a different algorithm (interior-point) in a future release.
> In fmincon at 509
In rsm1 at 10
Attempted to access Mu(1); index out of bounds because numel(Mu)=0.

Error in bata2 (line 3)
CC=((X(1)-Mu(1))/Sigma(1))^2+((X(2)-Mu(2))/Sigma(2))^2

Error in fmincon (line 635)
initVals.f = feval(funfcn{3},X,varargin{:});

Error in rsm1 (line 10)
[X,fval,exitflag,output]=fmincon(@bata2,X0,A,b,Aeq,beq,lb,ub,@st);

Caused by:
Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.

Answer 1

你好

由于你的目标函数里面有全局变量，所以你的主函数部分必须也要定义全局变量。

fmincon函数提示用默认的算法找不到最小，提示你使用active-set算法。

所以进行这两个修改后就能运行。主函数修改如下：

global Mu Sigama%加一句，必须要
Mu=[10, 2.5]; 
Sigama=[2, 0.375]; 
X0=[10, 2.5]; 
A=[];b=[];Aeq=[];beq=[];
lb=[];ub=[];
options=optimset('Algorithm','active-set');%加一句，改算法，不会有警告
[X,fval,exitflag,output]=fmincon(@bata2,X0,A,b,Aeq,beq,lb,ub,@st,options);
bata=sqrt(fval); 
Pf=cdf('norm',-bata,0.1);

希望对你有帮助

追问

非常感谢，已经解决了，这50分给你了。那在此基础上，如何修改一下程序，以循环调用响应面函数语句，最终使bata（k）-bata（k-1）的绝对值小于给定的精度0.001？如果可以，再追加50分，

追答

修改主函数如下，其他子函数不用修改。

clear,clc
global Mu Sigama%加一句，必须要
Mu=[10, 2]; 
Sigama=[2, 0.375]; 

A=[];b=[];Aeq=[];beq=[];
lb=[];ub=[];
options=optimset('Algorithm','active-set');%加一句，改算法，不会有警告
X0=[10, 2.5]; 
Max=10;%最大求解次数，可以设置一个较大的值
bata=zeros(1,Max);
    [X,fval,exitflag,output]=fmincon(@bata2,X0,A,b,Aeq,beq,lb,ub,@st,options);
    bata(1)=sqrt(fval);

for k=2:Max%循环调用响应面函数语句
    X0=X;%以上次解出的X为下次求解的初值。
    [X,fval,exitflag,output]=fmincon(@bata2,X0,A,b,Aeq,beq,lb,ub,@st,options);
    bata(k)=sqrt(fval);
    es=abs(bata(k)-bata(k-1));%判断bata（k）-bata（k-1）的绝对值小于给定的精度0.001
    if es<0.001
        break;
    end
end
Pf=cdf('norm',-bata,0.1);

Answer 2

结果是不是这样的：
CC =

0

ceq =

44.4431

CC =

5.5511e-015

ceq =

44.4431

CC =

9.8686e-015

ceq =

44.4431

CC =

6.9618

ceq =

4.2533e-007

CC =

6.9618

ceq =

3.5392e-007

CC =

6.9618

ceq =

1.4363e-006

CC =

6.9462

ceq =

-1.5684e-009

CC =

6.9462

ceq =

-7.4008e-008

CC =

6.9462

ceq =

1.0104e-006

CC =

6.9410

ceq =

2.2053e-009

CC =

6.9410

ceq =

-7.1262e-008

CC =

6.9410

ceq =

1.0152e-006