Matrix size and scalar problem using fmincon

Question

Paul AGAMENNONE 2022-11-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1860793-matrix-size-and-scalar-problem-using-fmincon

评论： Paul AGAMENNONE 2022-11-29

Hello,

I'm trying to run an optimization reliability problem using fmincon but I got a size problem when I integrate my function to search for the global reliability and thus fmincon cannot return a scalar value.

I don't know exactly where my problem stands so I ask for help.

Here is my code

muL = 2000;
sigL = 200;
R1 = 1-9.92*10^-5;
R2 = 1-1.2696*10^-4;
R3 = 1-3.87*10^-6;
Sr1_min = sqrt(((((1.5-1)*muL)/norminv(R1))^2)-(sigL)^2);
Sr1_max = sqrt(((((2.5-1)*muL)/norminv(R1))^2)-(sigL)^2);
Sr2_min = sqrt(((((1.5-1)*muL)/norminv(R2))^2)-(sigL)^2);
Sr2_max = sqrt(((((2.5-1)*muL)/norminv(R2))^2)-(sigL)^2);
Sr3_min = sqrt(((((1.5-1)*muL)/norminv(R3))^2)-(sigL)^2);
Sr3_max = sqrt(((((2.5-1)*muL)/norminv(R3))^2)-(sigL)^2);
lb = [Sr1_min,Sr2_min,Sr3_min];
ub = [Sr1_max,Sr2_max,Sr3_max];
A = [];
B = [];
Aeq = [];
Beq = [];
d0 = (lb+ub)/5;
fun = @(d) parameterfun(d,muL,sigL,R1,R2,R3);
const = @(d) nonlcon(d,muL,sigL,R1,R2,R3);
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options);
Warning: Inf or NaN value encountered.
Error using fmincon
Supplied objective function must return a scalar value.
function Rs = parameterfun(d,muL,sigL,R1,R2,R3)
%
mu_Sr1 = muL+norminv(R1)*sqrt((sigL)^2+(d(1))^2);
mu_Sr2 = muL+norminv(R2)*sqrt((sigL)^2+(d(2))^2);
mu_Sr3 = muL+norminv(R3)*sqrt((sigL)^2+(d(3))^2);
%
Y1_mean = muL-mu_Sr1;
Y2_mean = muL-mu_Sr2;
Y3_mean = muL-mu_Sr3;
%
Y1_std = sqrt((d(1))^2+(sigL)^2);
Y2_std = sqrt((d(2))^2+(sigL)^2);
Y3_std = sqrt((d(3))^2+(sigL)^2);
%
Y_mean = [Y1_mean Y2_mean Y3_mean];
Y_std = [(Y1_std^2) (sigL)^2 (sigL)^2; (sigL)^2 (Y2_std)^2 (sigL)^2; (sigL)^2 (sigL)^2 (Y3_std)^2];
inv_Y_std = inv(Y_std);
det_Y_std = det(Y_std);
fy = @(y) (1/(2*pi)^(3/2).*(det_Y_std)^0.5).*exp(-(1/2).*transpose(y-Y_mean).*inv_Y_std.*(y-Y_mean));
%
Rs = 1-integral(fy,-inf,0,'ArrayValued',true);
%pf = 1-Rs;
%
end
function [c,ceq] = nonlcon(d,muL,sigL,R1,R2,R3)
muL = 2000;
sigL = 200;
c(1) = 1.5 - ((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))/muL);
c(2) = 1.5 - ((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))/muL);
c(3) = 1.5 - ((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))/muL);
c(4) = ((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))/muL) - 2.5;
c(5) = ((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))/muL) - 2.5;
c(6) = ((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))/muL) - 2.5;
c(7) = 0.08 - (d(1)/((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))));
c(8) = 0.08 - (d(2)/((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))));
c(9) = 0.08 - (d(3)/((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))));
c(10) = (d(1)/((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2))))) - 0.2;
c(11) = (d(2)/((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2))))) - 0.2;
c(12) = (d(3)/((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2))))) - 0.2;
ceq = [];
end

Thank you in advance.

Paul

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2022-11-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1860793-matrix-size-and-scalar-problem-using-fmincon#answer_1109713

@Paul AGAMENNONE,

It appears that the value Rs returned by parameterfun() is a vector (or array). Rs is a vector because function fy(), inside the integral on the right hand side of Rs, is a vector (or array). Add another calculation inside parameterfun(), after Rs is computed, to convert the vector Rs to a scalar, which is a global measure of reliability. parameterfun() should return this scalar, instead of the vector Rs.

fmincon() minimizes a function. If you want to maximize global reliability, insert a minus sign, or a reciprocal, somewhere.

18 个评论
显示 16更早的评论隐藏 16更早的评论

Torsten 2022-11-23

编辑：Torsten 2022-11-23

在 MATLAB Online 中打开

Here is the code how you could implement fy on your own.

muL = 2000;
sigL = 200;
R1 = 1-9.92*10^-5;
R2 = 1-1.2696*10^-4;
R3 = 1-3.87*10^-6;
Sr1_min = sqrt(((((1.5-1)*muL)/norminv(R1))^2)-(sigL)^2);
Sr1_max = sqrt(((((2.5-1)*muL)/norminv(R1))^2)-(sigL)^2);
Sr2_min = sqrt(((((1.5-1)*muL)/norminv(R2))^2)-(sigL)^2);
Sr2_max = sqrt(((((2.5-1)*muL)/norminv(R2))^2)-(sigL)^2);
Sr3_min = sqrt(((((1.5-1)*muL)/norminv(R3))^2)-(sigL)^2);
Sr3_max = sqrt(((((2.5-1)*muL)/norminv(R3))^2)-(sigL)^2);
lb = [Sr1_min,Sr2_min,Sr3_min];
ub = [Sr1_max,Sr2_max,Sr3_max];
A = [];
B = [];
Aeq = [];
Beq = [];
d0 = (lb+ub)/5;
fun = @(d) parameterfun(d,muL,sigL,R1,R2,R3);
const = @(d) nonlcon(d,muL,sigL,R1,R2,R3);
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options)
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           4    1.000000e+00     3.244e-02     1.000e+00     0.000e+00     0.000e+00  
    1           8    1.000000e+00     4.236e-03     1.000e+00     1.414e+02     8.643e+01  
    2          12    1.000000e+00     1.074e-04     1.000e+00     2.334e+01     5.799e+00  
    3          16    1.000000e+00     7.345e-08     1.000e+00     6.230e-01     3.254e-02  
    4          20    1.000000e+00     3.401e-14     1.000e+00     4.268e-04     4.464e-06  
    5          24    1.000000e+00     0.000e+00     1.000e+00     1.977e-10     4.372e-13  

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
d = 1×3
  261.2446  259.2901  284.3916
fval = 1
function Rs = parameterfun(d,muL,sigL,R1,R2,R3)
%
mu_Sr1 = muL+norminv(R1)*sqrt((sigL)^2+(d(1))^2);
mu_Sr2 = muL+norminv(R2)*sqrt((sigL)^2+(d(2))^2);
mu_Sr3 = muL+norminv(R3)*sqrt((sigL)^2+(d(3))^2);
%
Y1_mean = muL-mu_Sr1;
Y2_mean = muL-mu_Sr2;
Y3_mean = muL-mu_Sr3;
%
Y1_std = sqrt((d(1))^2+(sigL)^2);
Y2_std = sqrt((d(2))^2+(sigL)^2);
Y3_std = sqrt((d(3))^2+(sigL)^2);
%
Y_mean = [Y1_mean Y2_mean Y3_mean];
Y_std = [(Y1_std^2) (sigL)^2 (sigL)^2; (sigL)^2 (Y2_std)^2 (sigL)^2; (sigL)^2 (sigL)^2 (Y3_std)^2];
inv_Y_std = inv(Y_std);
det_Y_std = det(Y_std);
%fy = @(X,Y,Z) arrayfun(@(x,y,z) 1/((2*pi)^(3/2).*(det_Y_std)^0.5).*exp(-(1/2).*([x,y,z]-Y_mean)*inv_Y_std*([x,y,z]-Y_mean).'),X,Y,Z);
%Rs = 1 - integral3(fy,-Inf,0,-Inf,0,-Inf,0);
Rs = 1 - integral3(@(x,y,z)fy(x,y,z,det_Y_std,inv_Y_std,Y_mean),-inf,0,-inf,0,-inf,0);
%Rs = 1 - mvncdf(zeros(1,3),Y_mean,Y_std);
%pf = 1-Rs;
%
end
function [c,ceq] = nonlcon(d,muL,sigL,R1,R2,R3)
muL = 2000;
sigL = 200;
c(1) = 1.5 - ((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))/muL);
c(2) = 1.5 - ((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))/muL);
c(3) = 1.5 - ((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))/muL);
c(4) = ((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))/muL) - 2.5;
c(5) = ((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))/muL) - 2.5;
c(6) = ((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))/muL) - 2.5;
c(7) = 0.08 - (d(1)/((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))));
c(8) = 0.08 - (d(2)/((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))));
c(9) = 0.08 - (d(3)/((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))));
c(10) = (d(1)/((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2))))) - 0.2;
c(11) = (d(2)/((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2))))) - 0.2;
c(12) = (d(3)/((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2))))) - 0.2;
ceq = [];
end
function values = fy(x,y,z,det_Y_std,inv_Y_std,Y_mean)
  values = zeros(size(x));
  for i=1:size(x,1)
    for j=1:size(x,2)
        for k=1:size(x,3)
            values(i,j,k) = 1/((2*pi)^(3/2).*(det_Y_std)^0.5).*exp(-(1/2).*([x(i),y(j),z(k)]-Y_mean)*inv_Y_std*([x(i),y(j),z(k)]-Y_mean).');
        end
    end
  end
end

Torsten 2022-11-25

编辑：Torsten 2022-11-25

在 MATLAB Online 中打开

I get the same constants for both methods with the code above.

But the integrals are just switched - the integral with mvncdf seems to be 1 - the integral with integral3. I don't know why. You should test both variants with inputs for mu and Sigma where you know the result.

muL = 2000;
sigL = 200;
R1 = 1-9.92*10^-5;
R2 = 1-1.2696*10^-4;
R3 = 1-3.87*10^-6;
Sr1_min = sqrt(((((1.5-1)*muL)/norminv(R1))^2)-(sigL)^2);
Sr1_max = sqrt(((((2.5-1)*muL)/norminv(R1))^2)-(sigL)^2);
Sr2_min = sqrt(((((1.5-1)*muL)/norminv(R2))^2)-(sigL)^2);
Sr2_max = sqrt(((((2.5-1)*muL)/norminv(R2))^2)-(sigL)^2);
Sr3_min = sqrt(((((1.5-1)*muL)/norminv(R3))^2)-(sigL)^2);
Sr3_max = sqrt(((((2.5-1)*muL)/norminv(R3))^2)-(sigL)^2);
lb = [Sr1_min,Sr2_min,Sr3_min];
ub = [Sr1_max,Sr2_max,Sr3_max];
A = [];
B = [];
Aeq = [];
Beq = [];
d0 = (lb+ub)/5;
fun = @(d) parameterfun(d,muL,sigL,R1,R2,R3);
const = @(d) nonlcon(d,muL,sigL,R1,R2,R3);
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
format long
flag = 1;
fun = @(d) parameterfun(d,muL,sigL,R1,R2,R3,flag);
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options)
Rs = 
     1
Rs = 
     1
Rs = 
     1
Rs = 
     1
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           4    1.000000e+00     3.244e-02     1.000e+00     0.000e+00     0.000e+00  
Rs = 
     1
Rs = 
     1
Rs = 
     1
Rs = 
     1
    1           8    1.000000e+00     4.236e-03     1.000e+00     1.414e+02     8.643e+01  
Rs = 
     1
Rs = 
     1
Rs = 
     1
Rs = 
     1
    2          12    1.000000e+00     1.074e-04     1.000e+00     2.334e+01     5.799e+00  
Rs = 
     1
Rs = 
     1
Rs = 
     1
Rs = 
     1
    3          16    1.000000e+00     7.345e-08     1.000e+00     6.230e-01     3.254e-02  
Rs = 
     1
Rs = 
     1
Rs = 
     1
Rs = 
     1
    4          20    1.000000e+00     3.401e-14     1.000e+00     4.268e-04     4.464e-06  
Rs = 
     1
Rs = 
     1
Rs = 
     1
Rs = 
     1
    5          24    1.000000e+00     0.000e+00     1.000e+00     1.977e-10     4.372e-13  

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
d = 1×3
   2.612446122521141   2.592901015034655   2.843915659497657
fval = 
     1
flag = 2;
fun = @(d) parameterfun(d,muL,sigL,R1,R2,R3,flag);
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options)
Rs = 
     2.261964601025790e-04
Rs = 
     2.261964601985023e-04
Rs = 
     2.261964602014999e-04
Rs = 
     2.261964601206756e-04
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           4    2.261965e-04     3.244e-02     1.000e+00     0.000e+00     3.381e-08  
Rs = 
     2.288078562057150e-04
Rs = 
     2.288078562411311e-04
Rs = 
     2.288078562413531e-04
Rs = 
     2.288078562107110e-04
    1           8    2.288079e-04     4.236e-03     1.000e+00     1.414e+02     8.643e+01  
Rs = 
     2.289743038106362e-04
Rs = 
     2.289743038415004e-04
Rs = 
     2.289743038417225e-04
Rs = 
     2.289743038146330e-04
    2          12    2.289743e-04     1.074e-04     1.000e+00     2.334e+01     5.799e+00  
Rs = 
     2.289783317231953e-04
Rs = 
     2.289783317539484e-04
Rs = 
     2.289783317540595e-04
Rs = 
     2.289783317270810e-04
    3          16    2.289783e-04     7.345e-08     1.000e+00     6.230e-01     3.254e-02  
Rs = 
     2.289783344753271e-04
Rs = 
     2.289783345060803e-04
Rs = 
     2.289783345063023e-04
Rs = 
     2.289783344793239e-04
    4          20    2.289783e-04     3.401e-14     1.000e+00     4.268e-04     4.471e-06  
Rs = 
     2.289783344751051e-04
Rs = 
     2.289783345058582e-04
Rs = 
     2.289783345060803e-04
Rs = 
     2.289783344791019e-04
    5          24    2.289783e-04     0.000e+00     1.000e+00     1.616e-08     8.017e-09  

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
d = 1×3
   2.612446122081119   2.592901014603244   2.843915659497657
fval = 
     2.289783344751051e-04
function Rs = parameterfun(d,muL,sigL,R1,R2,R3,flag)
%
mu_Sr1 = muL+norminv(R1)*sqrt((sigL)^2+(d(1))^2);
mu_Sr2 = muL+norminv(R2)*sqrt((sigL)^2+(d(2))^2);
mu_Sr3 = muL+norminv(R3)*sqrt((sigL)^2+(d(3))^2);
%
Y1_mean = muL-mu_Sr1;
Y2_mean = muL-mu_Sr2;
Y3_mean = muL-mu_Sr3;
%
Y1_std = sqrt((d(1))^2+(sigL)^2);
Y2_std = sqrt((d(2))^2+(sigL)^2);
Y3_std = sqrt((d(3))^2+(sigL)^2);
%
Y_mean = [Y1_mean Y2_mean Y3_mean];
Y_std = [(Y1_std^2) (sigL)^2 (sigL)^2; (sigL)^2 (Y2_std)^2 (sigL)^2; (sigL)^2 (sigL)^2 (Y3_std)^2];
inv_Y_std = inv(Y_std);
det_Y_std = det(Y_std);
%fy = @(X,Y,Z) arrayfun(@(x,y,z) 1/((2*pi)^(3/2).*(det_Y_std)^0.5).*exp(-(1/2).*([x,y,z]-Y_mean)*inv_Y_std*([x,y,z]-Y_mean).'),X,Y,Z);
%Rs = 1 - integral3(fy,-Inf,0,-Inf,0,-Inf,0);
if flag == 1
Rs = 1 - integral3(@(x,y,z)fy(x,y,z,det_Y_std,inv_Y_std,Y_mean),-Inf,0,-Inf,0,-Inf,0)
else 
Rs = 1 - mvncdf(zeros(1,3),Y_mean,Y_std)
end
%pf = 1-Rs;
%
end
function [c,ceq] = nonlcon(d,muL,sigL,R1,R2,R3)
muL = 2000;
sigL = 200;
c(1) = 1.5 - ((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))/muL);
c(2) = 1.5 - ((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))/muL);
c(3) = 1.5 - ((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))/muL);
c(4) = ((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))/muL) - 2.5;
c(5) = ((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))/muL) - 2.5;
c(6) = ((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))/muL) - 2.5;
c(7) = 0.08 - (d(1)/((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2)))));
c(8) = 0.08 - (d(2)/((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2)))));
c(9) = 0.08 - (d(3)/((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2)))));
c(10) = (d(1)/((muL+norminv(R1)*sqrt((d(1)^2)+(sigL^2))))) - 0.2;
c(11) = (d(2)/((muL+norminv(R2)*sqrt((d(2)^2)+(sigL^2))))) - 0.2;
c(12) = (d(3)/((muL+norminv(R3)*sqrt((d(3)^2)+(sigL^2))))) - 0.2;
ceq = [];
end
function values = fy(x,y,z,det_Y_std,inv_Y_std,Y_mean)
  values = zeros(size(x));
  for i=1:size(x,1)
    for j=1:size(x,2)
        for k=1:size(x,3)
            values(i,j,k) = 1/((2*pi)^(3/2).*(det_Y_std)^0.5).*exp(-(1/2).*([x(i),y(j),z(k)]-Y_mean)*inv_Y_std*([x(i),y(j),z(k)]-Y_mean).');
        end
    end
  end
end

Torsten 2022-11-29

I wanted to compare the results from fmincon when using integral3 and mvncdf in one run.

So I called fmincon first with flag = 1 to compute with integral3 and called fmincon thereafter with flag = 2 to compute with mvncdf.

Paul AGAMENNONE 2022-11-29

Ok I see thank you

请先登录，再进行评论。

Matrix size and scalar problem using fmincon

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

18 个评论
显示 16更早的评论隐藏 16更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Matrix size and scalar problem using fmincon

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

18 个评论 显示 16更早的评论隐藏 16更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

18 个评论
显示 16更早的评论隐藏 16更早的评论