# Optimization mvncdf and integral problem

3 Ansichten (letzte 30 Tage)
Paul AGAMENNONE am 13 Dez. 2022
Kommentiert: Torsten am 14 Dez. 2022
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
##### 0 Kommentare-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Torsten am 13 Dez. 2022
Bearbeitet: Torsten am 13 Dez. 2022
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
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 Kommentare5 ältere Kommentare anzeigen5 ältere Kommentare ausblenden
Torsten am 14 Dez. 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.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Find more on Point Cloud Processing in Help Center and File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by