Matrix size and scalar problem using fmincon

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

Answer 1

[William Rose]

@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.

[Paul AGAMENNONE]

Hi @William Rose

First thank you for your answer I understand more about the problem.

But what should be the form of the calculation I need to add to convert Rs to a scalar? I'm a beginner in Matlab so I don't know exatcly what should be the form.

Thanks

[Torsten]

Shouldn't this be matrix multiplication instead of elementwise multiplication ?

transpose(y-Y_mean).*inv_Y_std.*(y-Y_mean)

[William Rose]

@Paul AGAMENNONE,

I think @Torsten is right. I'll get to that.

To answer your question about how to convert Rs from a 3x3 array to a scalar, it would help me if you would describe Rs in words, or with an equation (other than the equation in your code). This would include a description in words of y and fy. Why do you integrate fy from y=-Inf to y=0?

fy includes the factor

transpose(y-Y_mean).*inv_Y_std.*(y-Y_mean)

where y-yMean is 3x1 and inv_Y_std is 3x3. The result of these element-wise multiplications is a 3x3 matrix. What is the meaning of this matrix?

I think @Torsten is correct: you probably want to do standard matrix multiplication:

transpose(y-Y_mean)*inv_Y_std*(y-Y_mean)

which will result in a scalar. I think the other factors in fy are scalars, so you can remove "'arrayValued',true" from the code. Then Rs will be a scalar.

[Paul AGAMENNONE]

@William Rose

Rs is the system reliability of a system with 3 components. In term of equation, Rs is equal to the cumulative derivative function (cdf) of the joint probability density function (which is fy in my case) of the components

I intregate from -inf to 0 because I want the probability that Y <0.

Also if I transform fy like this,

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

Matlab gives me an error "Arrays have incompatible sizes for this operation."

[Torsten]

@Paul AGAMENNONE

I think you want

Rs = 1 - mvncdf(zeros(1,3),Y_mean,Y_std)

[Paul AGAMENNONE]

@Torsten

I used mvncdf already (and it works) but I wanted to see if it would work with the general expression.

I was not sure if mvncdf was the same as my expression

[Torsten]

mvncdf works with the usual multivariate normal distribution. Isn't your fy the density function of the multivariate normal distribution with mean = Y_mean and Sigma = Y_std ?

Ah, I see, you multiply by sqrt(det_Y_std) and don't divide by this term. But I think it's just a coding error on your side.

[Paul AGAMENNONE]

yes it is

Ok I see I don't need to worry too much though my code is working with mvncdf

[Torsten]

I think you should stick to MATLAB's predefined solution with mvncdf.

But if you search for

Rs = Pr(Y<0)

,

Rs = mvncdf(zeros(1,3),Y_mean,Y_std)

, not

Rs = 1 - mvncdf(zeros(1,3),Y_mean,Y_std)

[Paul AGAMENNONE]

Ok I see

I put the 1-mvncdf because I was looking for the probability of failure, Rs is just the name I gave to not change every calculation.

Thanks for the help

[Torsten]

Isn't it a failure if Y<0 ?

[Torsten]

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

[William Rose]

Impressive, @Torsten.

[Paul AGAMENNONE]

@Torsten

I've used you code which is really clear, but why can't I find the same values when I use mvncdf? There are suppose to be the same but with different methods right?

[Torsten]

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

[Paul AGAMENNONE]

Hi @Torsten

Everything is working, thanks for your help.

Just one last question, what is exactly the use of flag in your coding? And why do I find the good result when the flag=2? I'm not familiar with this

[Torsten]

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]

Ok I see thank you