Weibull distribution reliability optimization model

Hello everyone,

I'm trying to determine the reliability of a system with unknown information, i.e. the distribution of random variables. To do so, I've used Weibull distribution to find the system reliability.

I have created a matlab function with unknown parameters called d that are optimized and at the end obtain the reliability. However, I'm not very good at Matlab and got some errors. I'll need some help to understand how to make my program works, and also know if what I'm doing is logic.

I got two mainly questions:

• How can I recalled my component reliability in the main function, everytime I got the error that one the variable is not recognized
• How can I optimize my Weibull parameters and get my system reliability? Matlab don't recognize the mean and standard deviation for Weibull distribution when I optimize them

I let my actual code, with a lot of mistake I know, to let you see.

Thank you very much for your help.

Paul

clearvars;
%Series system with 3 different components and 3 different loads
%(non-normal distribution)
%
%Information available to system designers
mu_L1 = 2100;
mu_L2 = 2000;
mu_L3 = 2000;
sig_L1 = 120;
sig_L2 = 90;
sig_L3 = 105;
w1 = 1;
w2 = 0.8;
w3 = 0.65;
%
%Optimization function
lb = [Sr1_min,Sr2_min,Sr3_min];

Unrecognized function or variable 'Sr1_min'.

ub = [Sr1_max,Sr2_max,Sr3_max];
A = [];
B = [];
Aeq = [];
Beq = [];
d0 = (lb+ub)/2;
fun = @(d) parameterfun(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,mu_Sr1,mu_Sr2,mu_Sr3,sig_Sr1,sig_Sr2,sig_Sr3,w1,w2,w3);
const = @(d) nonlcon(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,mu_Sr1,mu_Sr2,mu_Sr3,sig_Sr1,sig_Sr2,sig_Sr3,w1,w2,w3);
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options);
disp(d);
%
function Rs = parameterfun(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,mu_Sr1,mu_Sr2,mu_Sr3,sig_Sr1,sig_Sr2,sig_Sr3,R1,R2,R3,w1,w2,w3)
%
%Component reliability
R1 = wblpdf(d(1),d(4));
R2 = wblpdf(d(2),d(5));
R3 = wblpdf(d(3),d(6));
%
%Bounds of d
Sr1_min = sqrt(((((1.2-1)*w1*mu_L1)/norminv(R1))^2)-(w1*sig_L1)^2);
Sr1_max = sqrt(((((2-1)*w1*mu_L1)/norminv(R1))^2)-(w1*sig_L1)^2);
Sr2_min = sqrt(((((1.2-1)*w2*mu_L2)/norminv(R2))^2)-(w2*sig_L2)^2);
Sr2_max = sqrt(((((2-1)*w2*mu_L2)/norminv(R2))^2)-(w2*sig_L2)^2);
Sr3_min = sqrt(((((1.2-1)*w3*mu_L3)/norminv(R3))^2)-(w3*sig_L3)^2);
Sr3_max = sqrt(((((2-1)*w3*mu_L3)/norminv(R3))^2)-(w3*sig_L3)^2);
%
%Weibull parameters
[mu_Sr1, sig_Sr1] = wblstat(d(1),d(4));
[mu_Sr2, sig_Sr2] = wblstat(d(2),d(5));
[mu_Sr3, sig_Sr3] = wblstat(d(3),d(6));
%
Y1_mean = w1*mu_L1-mu_Sr1;
Y2_mean = w2*mu_L2-mu_Sr2;
Y3_mean = w3*mu_L3-mu_Sr3;
%
Y1_std = sqrt(sig_Sr1^2+(w1*sig_L1)^2);
Y2_std = sqrt(sig_Sr2^2+(w2*sig_L2)^2);
Y3_std = sqrt(sig_Sr3^2+(w3*sig_L3)^2);
%
Y_mean = [Y1_mean Y2_mean Y3_mean];
Y_std = [(Y1_std^2) (w1*w2*sig_L1*sig_L2) (w1*w3*sig_L1*sig_L3); (w2*w1*sig_L2*sig_L1) (Y2_std)^2 (w2*w3*sig_L2*sig_L3); (w3*w1*sig_L3*sig_L1) (w3*w2*sig_L3*sig_L2) (Y3_std)^2];
y = ones(size(d));
%
Rs = 1-mvncdf(y,Y_mean,Y_std);
%pf = 1-Rs;
%
end
%
function [c,ceq] = nonlcon(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,mu_Sr1,mu_Sr2,mu_Sr3,sig_Sr1,sig_Sr2,sig_Sr3,w1,w2,w3)
%
mu_L1 = 2100;
mu_L2 = 1600;
mu_L3 = 1750;
sig_L1 = 120;
sig_L2 = 230;
sig_L3 = 185;
%
c(1) = 1.2 - mu_Sr1/mu_L1;
c(2) = 1.2 - mu_Sr2/mu_L2;
c(3) = 1.2 - mu_Sr3/mu_L3;
c(4) = mu_Sr1/mu_L1 - 2;
c(5) = mu_Sr2/mu_L2 - 2;
c(6) = mu_Sr3/mu_L3 - 2;
c(7) = 0.06 - sig_Sr1/mu_Sr1;
c(8) = 0.09 - sig_Sr2/mu_Sr2;
c(9) = 0.09 - sig_Sr3/mu_Sr3;
c(10) = sig_Sr1/mu_Sr1 - 0.12;
c(11) = sig_Sr2/mu_Sr2 - 0.12;
c(12) = sig_Sr3/mu_Sr3 - 0.12;
%
ceq = [];
%
end
%

[Torsten]

I must admit that I don't have the slightest idea what you are trying to do in your code.

Could you add meaningful comments to the different sections of your program ?

[John D'Errico]

Its been far too many years since I worked with an expert in reliability to know how to fix your code. So I do see roughly what you are doing, but at the same time, if your code has bugs in it, then all we have is code that fails to work. I would strongly suggest that at the very least, you need to clearly explain what the code should do. That might help some one here, enough they could guide you to fix the code. But perhaps better is to consider posting this as a question on a site that has experts in reliability who visit it regularly.

You might consider the Cross Calidated site (under the stack exchange umbrella), but there are probably a few other places to search too.

[Paul AGAMENNONE]

Hello @Torsten and @John D'Errico,

Thank you for your answer, I must admit that myself don’t do what I’ve done with this code. I’ll try to explain more in details with a previous code that I made which is working.

In the code below, that is working, I optimized the reliability of a system composed of 3 components with random variables following normal distribution. In this case, I know the distribution of my variables.

%Series system with 3 different components and 3 different loads 
% 
%Information available to system designers
mu_L1 = 1400;
mu_L2 = 1500;
mu_L3 = 2000;
sig_L1 = 130;
sig_L2 = 150;
sig_L3 = 220;
R1 = 1-3.557*10^-4;
R2 = 1-2.4952*10^-4;
R3 = 1-7.1804*10^-4;
w1 = 0.5;
w2 = 1;
w3 = 0.35;
%
%Bounds of d
Sr1_min = sqrt(((((1.5-1)*w1*mu_L1)/norminv(R1))^2)-(w1*sig_L1)^2);
Sr1_max = sqrt(((((2.5-1)*w1*mu_L1)/norminv(R1))^2)-(w1*sig_L1)^2);
Sr2_min = sqrt(((((1.5-1)*w2*mu_L2)/norminv(R2))^2)-(w2*sig_L2)^2);
Sr2_max = sqrt(((((2.5-1)*w2*mu_L2)/norminv(R2))^2)-(w2*sig_L2)^2);
Sr3_min = sqrt(((((1.5-1)*w3*mu_L3)/norminv(R3))^2)-(w3*sig_L3)^2);
Sr3_max = sqrt(((((2.5-1)*w3*mu_L3)/norminv(R3))^2)-(w3*sig_L3)^2);
%
%Optimization function
lb = [Sr1_min,Sr2_min,Sr3_min];
ub = [Sr1_max,Sr2_max,Sr3_max];
A = [];
B = [];
Aeq = [];
Beq = [];
d0 = (lb+ub)/2;
fun = @(d) parameterfun(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,R1,R2,R3,w1,w2,w3);
const = @(d) nonlcon(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,R1,R2,R3,w1,w2,w3);
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options);

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           4    1.301109e-03     0.000e+00     1.000e+00     0.000e+00     5.929e-08  

Initial point is a local minimum that satisfies the constraints.

Optimization completed because at the initial point, the objective function is non-decreasing 
in feasible directions to within the value of the optimality tolerance, and 
constraints are satisfied to within the value of the constraint tolerance.

disp(d);

  191.8486  391.6546  199.2990

%
function Rs = parameterfun(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,R1,R2,R3,w1,w2,w3)
%
mu_Sr1 = w1*mu_L1+norminv(R1)*sqrt((w1*sig_L1)^2+(d(1))^2);
mu_Sr2 = w2*mu_L2+norminv(R2)*sqrt((w2*sig_L2)^2+(d(2))^2);
mu_Sr3 = w3*mu_L3+norminv(R3)*sqrt((w3*sig_L3)^2+(d(3))^2);
%
Y1_mean = w1*mu_L1-mu_Sr1;
Y2_mean = w2*mu_L2-mu_Sr2;
Y3_mean = w3*mu_L3-mu_Sr3;
%
Y1_std = sqrt((d(1))^2+(w1*sig_L1)^2);
Y2_std = sqrt((d(2))^2+(w2*sig_L2)^2);
Y3_std = sqrt((d(3))^2+(w3*sig_L3)^2);
%
Y_mean = [Y1_mean Y2_mean Y3_mean];
Y_std = [(Y1_std^2) (w1*w2*sig_L1*sig_L2) (w1*w3*sig_L1*sig_L3); (w2*w1*sig_L2*sig_L1) (Y2_std)^2 (w2*w3*sig_L2*sig_L3); (w3*w1*sig_L3*sig_L1) (w3*w2*sig_L3*sig_L2) (Y3_std)^2]; 
y = ones(size(d)); 
% 
Rs = 1-mvncdf(y,Y_mean,Y_std); 
%pf = 1-Rs; 
% 
end 
% 
function [c,ceq] = nonlcon(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,R1,R2,R3,w1,w2,w3)
%
mu_L1 = 1400;
mu_L2 = 1500;
mu_L3 = 2000;
sig_L1 = 130;
sig_L2 = 150;
sig_L3 = 220;
%
c(1) = 1.5 - ((w1*mu_L1+norminv(R1)*sqrt((d(1)^2)+((w1*sig_L1)^2)))/(w1*mu_L1));
c(2) = 1.5 - ((w2*mu_L2+norminv(R2)*sqrt((d(2)^2)+((w2*sig_L2)^2)))/(w2*mu_L2));
c(3) = 1.5 - ((w3*mu_L3+norminv(R3)*sqrt((d(3)^2)+((w3*sig_L3)^2)))/(w3*mu_L3));
c(4) = ((w1*mu_L1+norminv(R1)*sqrt((d(1)^2)+((w1*sig_L1)^2)))/(w1*mu_L1)) - 2.5;
c(5) = ((w2*mu_L2+norminv(R2)*sqrt((d(2)^2)+((w2*sig_L2)^2)))/(w2*mu_L2)) - 2.5;
c(6) = ((w3*mu_L3+norminv(R3)*sqrt((d(3)^2)+((w3*sig_L3)^2)))/(w3*mu_L3)) - 2.5;
c(7) = 0.09 - (d(1)/((w1*mu_L1+norminv(R1)*sqrt((d(1)^2)+((w1*sig_L1)^2)))));
c(8) = 0.09 - (d(2)/((w2*mu_L2+norminv(R2)*sqrt((d(2)^2)+((w2*sig_L2)^2)))));
c(9) = 0.09 - (d(3)/((w3*mu_L3+norminv(R3)*sqrt((d(3)^2)+((w3*sig_L3)^2)))));
c(10) = (d(1)/((w1*mu_L1+norminv(R1)*sqrt((d(1)^2)+((w1*sig_L1)^2))))) - 0.2;
c(11) = (d(2)/((w2*mu_L2+norminv(R2)*sqrt((d(2)^2)+((w2*sig_L2)^2))))) - 0.2;
c(12) = (d(3)/((w3*mu_L3+norminv(R3)*sqrt((d(3)^2)+((w3*sig_L3)^2))))) - 0.2; 
% 
ceq = []; 
% 
end 
%

For this different case, I’m working also with the same logic but with unknown information about variables distribution. It means that the component mean and standard deviation are unknown, as well as component reliability. Therefore, I used Weibull distribution to describe my unknown values. Then, I have to determine the unknown parameters of variables using my optimization model and be able to obtain the reliability. However, I’m not sure that mathematically I can follow the same principle used with normal distribution, this is why I tried the same code but I’m stuck with it.

I hope you understand a bit better what I hope to achieve, I tried to be as clear as possible.

Thank you for your help