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
采纳的回答
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 个评论
Paul AGAMENNONE
2022-12-13
Torsten
2022-12-13
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
2022-12-13
Paul AGAMENNONE
2022-12-14
Torsten
2022-12-14
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Fit Postprocessing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
