Optimization mvncdf and integral problem
Show older comments
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
Accepted Answer
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
6 Comments
Paul AGAMENNONE
on 13 Dec 2022
Torsten
on 13 Dec 2022
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
on 13 Dec 2022
Paul AGAMENNONE
on 14 Dec 2022
Torsten
on 14 Dec 2022
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.
More Answers (0)
Categories
Find more on Fit Postprocessing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
