Optimization mvncdf and integral problem

Hello,

I'm running an optimzation model to minimize the system reliability of my problem. The model works correctly as I'm looking for the right values of d with mvncdf and using the simple method of the integral with the pdf. However, I can't find the same reliability at the end of my problem. With mvncdf, I get the good value, but with the integrals method I find 1 (or 0 depending if I take the probability of failure or reliability).

I don't know where the problem comes from, so I ask for help.

Here is my code: change the value of the flag to get both methods

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)

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)
disp(d);
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/sqrt((2*pi)^(3/2).*det_Y_std).*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

Answer 1

[Torsten]

There were some errors in the function supplied to compute the MND.

Use

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).*sqrt(abs(det_Y_std))).*exp(-(1/2).*([x(i,j,k),y(i,j,k),z(i,j,k)]-Y_mean)*inv_Y_std*([x(i,j,k),y(i,j,k),z(i,j,k)]-Y_mean).');
        end
    end
  end
end

instead.

You can also try the call

Rs = 1-integral3(@(x,y,z)arrayfun(@(X,Y,Z)fy(X,Y,Z,det_Y_std,inv_Y_std,Y_mean),x,y,z),-Inf,0,-Inf,0,-Inf,0)

together with the function

function value = fy(x,y,z,det_Y_std,inv_Y_std,Y_mean)
  value = 1/((2*pi)^(3/2).*sqrt(abs(det_Y_std))).*exp(-(1/2).*([x,y,z]-Y_mean)*inv_Y_std*([x,y,z]-Y_mean).');
end

[Paul AGAMENNONE]

Oh I see, but now my Matlab is running without giving me any results. Did it take times to you to compute the whole problem?

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)

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)arrayfun(@(X,Y,Z)fy(X,Y,Z,det_Y_std,inv_Y_std,Y_mean),x,y,z),-Inf,0,-Inf,0,-Inf,0)
else 
Rs = 1 - mvncdf(zeros(1,3),Y_mean,Y_std)
disp(d);
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 value = fy(x,y,z,det_Y_std,inv_Y_std,Y_mean)
  value = 1/((2*pi)^(3/2).*sqrt(abs(det_Y_std))).*exp(-(1/2).*([x,y,z]-Y_mean)*inv_Y_std*([x,y,z]-Y_mean).');
end

[Torsten]

I tested integral3 and mvncdf for the standard 3d normal distribution and both gave the same results.

I didn't test integral3 for your input matrices and determinants.

integral3 is a general-purpose integrator over 3d domains that is not especially optimized for the gaussian probability density function, but mvncdf is. So why don't you want to stick to a program that is best suited for your purpose ?

[Paul AGAMENNONE]

I wanted to try with the classic method but I have seen that it is too complicated. However I still have a question. With mvncdf I minimizing the function and thus finding the minimum reliability. I've seen that for maximizing f(x) I need to minimize -f(x) but it has no effect with mvncdf, how should I do?

[Torsten]

Maybe your problem has only one feasible point ?

I don't know, but that's an observation that would prevent me from trusting in the results.

[Paul AGAMENNONE]

I found the solution, it comes from the definition of y in my example with mvncdf

y = zeros(size(d));

Thanks for the help

[Torsten]

But there is no such setting in the code above and nevertheless you get the same solution when you try to minimize the objective as you get when you try to maximize it.