Optimization mvncdf and integral problem

1 次查看（过去 30 天）
Paul AGAMENNONE
Paul AGAMENNONE 2022-12-13
评论： Torsten 2022-12-14
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

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-12-13
编辑：Torsten 2022-12-13
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 个评论
显示 4更早的评论
Paul AGAMENNONE
Paul AGAMENNONE 2022-12-14
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
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Genomics and Next Generation Sequencing 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by