Optimization model definition problem

2 次查看（过去 30 天）
Paul AGAMENNONE
Paul AGAMENNONE 2023-4-20
回答： albara 2023-4-28
Hello,
I'm running an optimization model to minimize the system reliability of my problem. My variables are defined following Weibull distribution with unknown parameters, this what I try to optimize to determine my system reliability. I encounter two problems:
  • I got some errors in the definition of mvtcdf. I don't know if this function can replace the actual form of a multivariate cdf with Weibull.
  • I'm also stuck in defining my lower and upper bounds. I got my reliability model and I know how the form that sigma minimum should have, but I can't find the relation. (see note at the end)
I let my actual code and some pictures of my equations if anyone can help me to unblock this situation, it would help me a lot.
Thank you in advance
clearvars;
%Series system with 3 different components and 3 different loads
%(non-normal distribution)
%
%Information available to system designers
mu_L1 = 2000;
mu_L2 = 2000;
mu_L3 = 2000;
sig_L1 = 110;
sig_L2 = 110;
sig_L3 = 105;
R1 = 1.365*10^-4;
R2 = 3.983*10^-4;
R3 = 3.152*10^-4;
w1 = 1;
w2 = 0.8;
w3 = 0.65;
%
%Bounds of d
mu_Sr1_min = 0;
mu_Sr1_max = 10000;
mu_Sr2_min = 0;
mu_Sr2_max = 10000;
mu_Sr3_min = 0;
mu_Sr3_max = 10000;
sig_Sr1_min = sqrt(((((0.8-1)*w1*mu_L1)/wblinv(R1))^2)-(w1*sig_L1)^2);
sig_Sr1_max = sqrt(((((1.4-1)*w1*mu_L1)/wblinv(R1))^2)-(w1*sig_L1)^2);
sig_Sr2_min = sqrt(((((0.8-1)*w2*mu_L2)/wblinv(R2))^2)-(w2*sig_L2)^2);
sig_Sr2_max = sqrt(((((1.4-1)*w2*mu_L2)/wblinv(R2))^2)-(w2*sig_L2)^2);
sig_Sr3_min = sqrt(((((0.8-1)*w3*mu_L3)/wblinv(R3))^2)-(w3*sig_L3)^2);
sig_Sr3_max = sqrt(((((1.4-1)*w3*mu_L3)/wblinv(R3))^2)-(w3*sig_L3)^2);
%
%Optimization function
lb = [mu_Sr1_min,sig_Sr1_min,mu_Sr2_min,sig_Sr2_min,mu_Sr3_min,sig_Sr3_min];
ub = [mu_Sr1_max,sig_Sr1_max,mu_Sr2_max,sig_Sr2_max,mu_Sr3_max,sig_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);
disp(d);
%
function Rs = parameterfun(d,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,R1,R2,R3,w1,w2,w3)
%
mu_Sr1 = d(1)*gamma(1+(1/d(2)));
mu_Sr2 = d(3)*gamma(1+(1/d(4)));
mu_Sr3 = d(5)*gamma(1+(1/d(6)));
sig_Sr1 = sqrt(((d(1))^2)*gamma(1+(2/d(2)))-d(1)*gamma(1+(1/d(2))));
sig_Sr2 = sqrt(((d(3))^2)*gamma(1+(2/d(4)))-d(3)*gamma(1+(1/d(4))));
sig_Sr3 = sqrt(((d(5))^2)*gamma(1+(2/d(6)))-d(5)*gamma(1+(1/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);
%
%Mean vector & Covariance matrix
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));
%
%System reliability
Rs = 1-mvtcdf(y,Y_mean,Y_std);
%pf = 1-Rs;
%
end
%
function [c,ceq] = nonlcon(d,mu_Sr1,mu_Sr2,mu_Sr3,sig_Sr1,sig_Sr2,sig_Sr3,mu_L1,sig_L1,mu_L2,sig_L2,mu_L3,sig_L3,R1,R2,R3,w1,w2,w3)
%
mu_L1 = 2100;
mu_L2 = 2000;
mu_L3 = 2000;
sig_L1 = 120;
sig_L2 = 90;
sig_L3 = 105;
%
%Constraint functions
c(1) = 0.8 - ((d(1)*gamma(1+(1/d(2))))/mu_L1);
c(2) = 0.8 - ((d(3)*gamma(1+(1/d(4))))/mu_L1);
c(3) = 0.8 - ((d(5)*gamma(1+(1/d(6))))/mu_L1);
c(4) = ((d(1)*gamma(1+(1/d(2))))/mu_L1) - 1.4;
c(5) = ((d(3)*gamma(1+(1/d(4))))/mu_L1) - 1.4;
c(6) = ((d(5)*gamma(1+(1/d(6))))/mu_L1) - 1.4;
c(7) = 0.09 - (sqrt(((d(1))^2)*gamma(1+(2/d(2)))-d(1)*gamma(1+(1/d(2)))))/(d(1)*gamma(1+(1/d(2))));
c(8) = 0.09 - (sqrt(((d(3))^2)*gamma(1+(2/d(4)))-d(1)*gamma(1+(1/d(3)))))/(d(1)*gamma(1+(1/d(4))));
c(9) = 0.09 - (sqrt(((d(5))^2)*gamma(1+(2/d(6)))-d(1)*gamma(1+(1/d(5)))))/(d(1)*gamma(1+(1/d(6))));
c(10) = (sqrt(((d(1))^2)*gamma(1+(2/d(2)))-d(1)*gamma(1+(1/d(2)))))/(d(1)*gamma(1+(1/d(2)))) - 0.2;
c(11) = (sqrt(((d(3))^2)*gamma(1+(2/d(4)))-d(1)*gamma(1+(1/d(3)))))/(d(1)*gamma(1+(1/d(4)))) - 0.2;
c(12) = (sqrt(((d(5))^2)*gamma(1+(2/d(6)))-d(1)*gamma(1+(1/d(5)))))/(d(1)*gamma(1+(1/d(6)))) - 0.2;
%
ceq = [];
%
end
%
Considering the bounds, I've done it with normal distribution (see doc attached) but it was simplified as mu was replaced by other variables. In this case it doesn't work, so I'm stuck to determine it.

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

回答（1 个）

albara
albara 2023-4-28
It seems that you have two main issues: defining the mvtcdf function and determining the lower and upper bounds for the optimization problem.
Regarding the mvtcdf function, MATLAB does not have a built-in function to calculate the cumulative distribution function (CDF) of the multivariate Weibull distribution. You might need to implement a custom function or use a third-party package to compute the CDF of the multivariate Weibull distribution. Alternatively, you can use copulas to model the dependence structure between the Weibull-distributed random variables.
As for the bounds, I noticed that you're using the same mu values for mu_L1, mu_L2, and mu_L3 in the nonlcon function as you have defined in the main script. The mu values should be the same in both places to maintain consistency in the problem. If they are different, it might cause issues with the constraints and the optimization process.
Regarding your note on determining the bounds, it is unclear how you have derived the bounds for the sigma values. It might be helpful to revisit your problem formulation and the assumptions you have made in deriving the bounds. Understanding the physical constraints of your system and the properties of the Weibull distribution might help you come up with a reasonable range for the variables.
If you still have issues with the bounds, it might be helpful to provide more information on the problem and the assumptions you have made. This would allow for more specific guidance in determining the appropriate bounds for the optimization problem.
Important: There may be some mistakes in this answer Experts can tell if there are any mistakes

类别

在 Help Center 和 File Exchange 中查找有关 Exploration and Visualization 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by