Using Mathematica and webMathematica in Engineering Research:
Symbolic Evaluation for Statistical Distribution Modelling

Arvind Rajan
PhD Candidate, Monash University, Malaysia

Reliability and Robustness in Engineering Design

In the increasingly complex and competitive modern technological environment, engineers from a wide variety of fields constantly strive toward finding the most economical designs by using the least amount of resources.

For high quality designs, two additional factors are considered during the process:

Robustness; and

Reliability.

Where are they used?

Distributions

Knowing the distribution of the quantity of interest is the key to the assessment of robustness and reliability of a system-under-design

Mean

Provides the expected value of a quantity

Plot[Table[PDF[NormalDistribution[μ,1.5],x],{μ,{-3,3,10}}]//Evaluate,{x,-10,15},Filling->Axis]

Standard Deviation

Directly related to the robustness - Smaller standard deviation means more robust

Plot[Table[PDF[NormalDistribution[0,σ],x],{σ,{.8,1,2}}]//Evaluate,{x,-6,6},Filling->Axis]

Other Parameters

Skewness - the measure of the asymmetry of a distribution

Kurtosis - the measure of the “peak” and “tail” of a distribution

Research at Monash University, Malaysia

One of the intensive research area at Monash University, Malaysia is reliability analysis of complex engineering systems:

Finding the parameters of the distribution using mathematical tools

The parameters are then used to reconstruct the distribution to study the robustness and reliability of the system-under-design

System Representation using Mathematical Notations

A system Y is conventionally expressed using a mathematical function g as

Y=g((Subscript[X, 1]),(Subscript[X, 2]),... ,(Subscript[X, N]))

Example:

Y=(Subscript[X, 1])+(Subscript[X, 2])

Distributions (and their parameters) of the inputs X are usually known; and distribution of Y is sought after

Research Aim

The aim of our research is to automate the process of finding the distribution parameters of Y in terms of the distribution parameters of the inputs X

Mathematica: Automated Framework

Algorithm

(*first stage, which is to be the input of function file later*)
NoVar = 2; (*number of variables*)

(*initialisation & declaration of the symbolic terms*)
ex={x1,x2}; (*initialise variable*)

(*building the struct*)
Symmetry = {"symmetrical","symmetrical"};
Distribution = {"normal","normal"};
Pmeter = {{m1,s1},{m2,s2}};

(*writing down the equation*)
P =(x1^2)+(x2^2);(*from the user input box*)

M={1,2,3,4}; (*number of Moment*)
(*----------------------------------------------------------------------------------------*)
(*derivation starts*)
(*----------------------------------------------------------------------------------------*)
V=Table[v[i],{i,NoVar}]; (*Untranslated standard variable*)
EP=Table[i/i, {i,Length[M]}];(*Mth moment of input polynomial*)
PolyError = 0;

(*Main loop begins*)
For[
ptM=1,ptM≤Length[M],ptM++,
(*for loop begins*)
PIn =Expand[P^(M[[ptM]])];

f[HoldPattern[+z__]]:={z};
f[else_]:={else};
segment = f[PIn];

B = Table[i/i,{i,Length[segment]}];
T = Table[i/i,{i,Length[segment]}];

For[(*For loop to find the monomials and coefficients*)
pos = 1,pos≤Length[segment],pos++,
ExpVar = Rationalize[Exponent[segment[[pos]],ex]];

For[ (*For loop to extract the monomials out*)
part=1,part≤Length[ExpVar],part++,
T[[pos]] = T[[pos]]*ex[[part]]^ExpVar[[part]];

];

B[[pos]]=segment[[pos]]/T[[pos]] ;(*Finding corresponding coefficient*)
];

ET=Table[i/i, {i,Length[T]}]; (*Expectation of each monomial in P*)

(*For each monomial*)
For[
index=1,index≤Length[T],index++,
(*for loop begins*)
SymList = Variables[T[[index]]]; (*Decouple product terms*)
(*For finding moment of each independent variable*)

For[
jndex = 1,jndex≤Length[SymList],jndex++,
(*for loop begins*)
y = SymList[[jndex]]; (*Isolate independent variable*)
m = Exponent[T[[index]],y];(*Finding the degree of the variable*)

ptx = Position[ex,y][[1]]; (*Identify variable being isolated*)
SymType=Symmetry[[ptx[[1]]]]; (*Finding the symmetry type*)
Param=Pmeter[[ptx[[1]]]]; (*Finding the parameters of the variable*)

If[ (*For uniform distribution not 1 sided*)!(((StringMatchQ[(Distribution[[ptx[[1]]]]) ,"uniform"]∨StringMatchQ[(Distribution[[ptx[[1]]]]) ,"triangular"]∨StringMatchQ[(Distribution[[ptx[[1]]]]),"trapezoidal"]) && (StringMatchQ[SymType,"left sided"]∨StringMatchQ[SymType,"right sided"]))∨(StringMatchQ[(Distribution[[ptx[[1]]]]),"beta"]&&m≤(-1))),

(*For un-symmetrical 2-sided*)
If[SymType==3,EvalCount=2,EvalCount=1];

For[(*kndex*)
kndex = 1,kndex≤EvalCount,kndex++,
Q = Expand[(V[[ptx[[1]]]]+Param[[1]])^m];
tegment=MonomialList[Q];
Test = Quiet[Check[coef=CoefficientList[Q,V[[ptx[[1]]]]],1]];
If[Test == 1, PolyError = PolyError+Test];
vm = tegment/coef;

If[kndex == 1,
Evm=Table[EVM[i],{i,Length[vm]}];,
Evm;];


For[(*lndex*)
lndex = 1,lndex≤Length[vm],lndex++,
Switch[vm[[lndex]],
1,
ni=0;,
_,
ni = Exponent[vm[[lndex]],Variables[vm[[lndex]]][[1]]];
];(*End of Switch case*)

(*Select formula according to distribution type*)
Switch[ToLowerCase[(Distribution[[ptx[[1]]]])],
"uniform",
EvmTemp = Param[[2]]^ni/(((ni+1)*2));,
"trapezoidal",
Shape = Param[[3]];
EvmTemp = Param[[2]]^ni/(ni^2+3*ni+2)*(1-Shape^(ni+2))/(1-Shape^2);,
"triangular",
EvmTemp=Param[[2]]^ni/(ni^2+3*ni+2);,
"beta",
Shape1 = Param[[3]];
Shape2 = Param[[4]];
gain = 1/2;
For[ptm= 1,ptm≤m,ptm++,gain = gain*(Shape1+ptm-1)/(Shape1+Shape2+ptm-1);]; (*End of ptm For loop*)
EvmTemp = (Param[[2]]^ni)*gain;,
"normal",
If[Mod[ni,2]==0,EvmTemp=(Param[[2]]^ni)*Product[i,{i,Range[1,ni-1,2]}]/2;,EvmTemp=(Param[[2]]^ni)*Product[i,{i,Range[1,(ni-1)/2,1]}]*2^(ni/2-1)/Sqrt[Pi];](*End of If statement*),
"lognormal",
EvmTemp=Exp[((ni*Param[[2]])^2)/2]/2;(*lognormal here*),
"student",
Shape=Param[[3]];
EvmTemp=(Param[[2]]^ni)*(Shape^(ni/2)*Gamma[(ni+1)/2]*Gamma[((Shape-ni)/2)]/(2*Sqrt[Pi]*Gamma[Shape/2]));,
"laplace",
EvmTemp=(Param[[2]]^ni)*Factorial[ni]/2;,
"gamma",
Shape=Param[[3]];
EvmTemp=(Param[[2]]^ni)*Product[i,{i,Range[Shape,Shape+ni-1,1]}]/2;,
"weibull",
Shape=Param[[3]];
EvmTemp=(Param[[2]]^ni)*Gamma[ni/Shape+1]/2;,
"maxwell",
If[Mod[ni,2]==0,EvmTemp=(Param[[2]]^ni)*Product[i,{i,Range[1,ni+1,2]}];,EvmTemp=(Param[[2]]^ni)*(2^(ni/2+1)/Sqrt[Pi])*Factorial[(ni+1)/2];](*End of If statement*),
_,
(*Default*)
Print["Unknown distribution type!"];
Abort[];
](*End of Switch case*)


(*Adjust the moment according to the stated symmetry*)
Switch[SymType,
"symmetrical"(*Symmetrical Distribution*),
If[Mod[ni,2]==0,
Evm[[lndex]]=2*EvmTemp;,Evm[[lndex]]=0;](*End of If statement*),
"right sided"(*One Sided Distribution (right)*),
Evm[[lndex]]=2*EvmTemp;,
"left sided"(*One Sided Distribution (left)*),
Evm[[lndex]]=((-1)^ni)*2*EvmTemp;
](*End of Switch case*)
](*End of If statement (lndex)*)
];(*End of If statement (kndex)*)

ET[[index]]=Total[ET[[index]]*coef*Evm];,
(*#### ELSE STATEMENT HERE!!! ####*)
(*SPECIAL ROUTINE TO CALCULATE THE MOMENT OF FINITE DISTRIBUTION*)

Switch[ToLowerCase[Distribution[[ptx[[1]]]]],
"uniform",
mean = Param[[1]];
scale = Param[[2]];
If[m==(-1),Exm=Log[(mean+scale)/(mean-scale)]/(2*scale),Exm=Expand[(mean+scale)^(m+1)-(mean-scale)^(m+1)]/(2*scale*(m+1))];,
"beta",
mean = Param[[1]];
scale = Param[[2]];
scale1 = Param[[3]];
scale2 = Param[[4]];
Exm=(mean^m)*(Hypergeometric2F1[scale1,-m,scale1+scale2,(-scale/mean)]);,
"triangular",
mean = Param[[1]];
scale = Param[[2]];
Switch[m,
(-1),
Exm=((Log[(mean+scale)/mean]-Log[(mean-scale)/mean])/scale)+((mean*(Log[(mean+scale)/mean]-Log[(mean-scale)/mean]))/scale^2);,
(-2),
Exm=(-Log[(mean+scale)/mean]-Log[(mean-scale)/mean])/scale^2;,
_,
Exm=(((mean+scale)^m-(mean-scale)^m)/(m^2+3*m+2))-((scale*(2*mean^(m+1)-2*mean*(mean+scale)^m+m*mean^(m+1))-((mean^2)*(mean+scale)^m)+mean^(m+2))/((scale^2)*(m^2+3*m+2)))+((scale*(2*mean^(m+1)-2*mean*(mean-scale)^m+m*mean^(m+1))+((mean^2)*(mean-scale)^m)-mean^(m+2))/((scale^2)*(m^2+3*m+2)));
];(*end of switch case*),
"trapezoidal",
mean = Param[[1]];
scale = Param[[2]];
beta1 = Param[[3]];
Switch[m,
(-1),
Exm=(Log[(mean+beta1)/(mean-beta1)]/(scale*(1+beta1)))+((beta1-scale*(Log[(mean+beta1)/(mean+scale)]+1)-(mean*Log[beta1+mean])+(mean*Log[scale+mean]))/((scale^2)*(1+(beta1^2))))-((beta1-scale*(Log[(mean-beta1)/(mean-scale)]+1)+(mean*Log[mean-beta1])-(mean*Log[mean-scale]))/((scale^2)*(1+(beta1^2))));,
(-2),
Exm=((scale-beta1*(Log[(mean+scale)/(mean+beta1)]+1)+(mean*Log[beta1+mean])-(mean*Log[scale+mean]))/((scale^2)*(1+(beta1^2))*(beta1-mean)))+((scale-beta1*(Log[(mean-scale)/(mean-beta1)]+1)-(mean*Log[mean-beta1])+(mean*Log[mean-scale]))/((scale^2)*(1+(beta1^2))*(beta1-mean)))-(((1/(beta1+mean))+(1/(beta1-mean)))/(scale*(beta1+1)));,
_,
Exm=(((((mean-beta1)^(m+1)-(mean-scale)^(m+1))/(m+1)))/(scale*(beta^2+1)))+((((mean-beta1)^(m+2)-(mean-scale)^(m+2))/(m+2))/((scale^2)*(beta1^2+1)))-(((((beta1+mean)^(m+1)-(mean+scale)^(m+1))/(m+1)))/(scale*(beta^2+1)))+((((mean+beta1)^(m+2)-(mean+scale)^(m+2))/(m+2))/((scale^2)*(beta1^2+1)))-((mean*((((mean-beta1)^(m+1)-(mean-scale)^(m+1))/(m+1))))/((scale^2)*(beta1^2+1)))-((mean*((((beta1+mean)^(m+1)-(mean+scale)^(m+1))/(m+1))))/((scale^2)*(beta1^2+1)))-(((mean-beta1)^(m+1)-(beta1+mean)^(m+1))/(scale*(beta1+1)*(m+1)));
];(*end of switch case*),
_,
(*Default*)
Print["Unknown distribution type!"];
Abort[];
];(*End of Switch case*)

ET[[index]]=ET[[index]]*Exm; (*changed here*)
](*End of If statement*)
](*End of loop of finding moment of each independant variable*)
];(*End of loop for each monomial*)
EP[[ptM]]=Total[B*ET];
];(*End of main loop*)

(*CALCULATION FOR UNCERTAINTY STARTS HERE*)
VarP=DeleteDuplicates@Cases[P,_Symbol,Infinity];

If[SubsetQ[ex,VarP]==False,UPSQR="check equation typed in";,UPSQR = EP[[2]]-(EP[[1]])^2;];
If[PolyError>0,Print["Bad"]];

InputForm[Simplify[EP[[2]]-EP[[1]]^2]]

webMathematica

A webpage was created to promote the automated framework

In the webpage, a graphical user interface (GUI) was built

The GUI interacts with the user and Mathematica kernel via webMathematica

The tool also recognized and being promoted by the IEEE Instrumentation and Measurement Society, Technical Committee 32 - Fault Tolerance Measurement Systems

http://aue.erc.monash.edu.au/webMathematica/Asue/Evaluate.jsp

A Case Study: Car Side Impact Crashworthiness Study

A complex real-world design optimization problem with reliability analysis

10 constraints and 11 design variables

Aim is to build a safe car by minimising the vehicle weight subject to 90% reliability

Distributions of the constraints are found using the tool published using webMathematica

Results

Achievements

Five high-quality journal publications in the Institute of Electrical and Electronics Engineers (IEEE) and the Institution of Engineering and Technology (IET)

Five conference publications in international IEEE conferences

Future Project

Reliability analysis on structural designs subject to loading from earthquake

In collaboration with Massey University, New Zealand

Acknowledgement

OEMS IntiPakar Corporation Sdn Bhd

Thank you
&
Have a nice day