Error in Fmincon, objective function must return scalar

Hello everyone,

I'm trying to run an optimization model with Weibull distribution and unknown distribution parameters about components. To do so, I'm trying to optimize my functions to determine the unknown parameters and the system reliability. I'm using fmincon but I got an error. I know it musts return a scalar value but I don't know where I am wrong in the definition of my functions.

I let my code for you to see, if someone can enlighten me it would be a great help.

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;
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 = 0;
sig_Sr1_max = 1000;
sig_Sr2_min = 0;
sig_Sr2_max = 1000;
sig_Sr3_min = 0;
sig_Sr3_max = 1000;

%
%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);
const = @(d) nonlcon(d,mu_L1,mu_L2,mu_L3,mu_Sr1,sig_Sr1,mu_Sr2,sig_Sr2,mu_Sr3,sig_Sr3);
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_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))-(gamma(1+1/d(2)))^2));
sig_Sr2 = sqrt(d(3)^2*(gamma(1+2/d(4))-(gamma(1+1/d(4)))^2));
sig_Sr3 = sqrt(d(5)^2*(gamma(1+2/d(6))-(gamma(1+1/d(6)))^2));
%
F_S1 = wblpdf(d,mu_Sr1,sig_Sr1);
F_S2 = wblpdf(d,mu_Sr2,sig_Sr2);
F_S3 = wblpdf(d,mu_Sr3,sig_Sr3);
%
%Define joint probability function
F_jpdf = @(d) F_S1.*F_S2.*F_S3;
%
%System reliability
Rs = 1-integral(F_jpdf,0,1,'ArrayValued',true);
%pf = 1-Rs;
%
end
%
function [c,ceq] = nonlcon(d,mu_L1,mu_L2,mu_L3,mu_Sr1,sig_Sr1,mu_Sr2,sig_Sr2,mu_Sr3,sig_Sr3)
%
%Constraint functions
c(1) = 0.8 - ((d(1)*gamma(1+(1/d(2))))/(w1*mu_L1));
c(2) = 0.8 - ((d(3)*gamma(1+(1/d(4))))/(w2*mu_L2));
c(3) = 0.8 - ((d(5)*gamma(1+(1/d(6))))/(w3*mu_L3));
c(4) = ((d(1)*gamma(1+(1/d(2))))/(w1*mu_L1)) - 1.4;
c(5) = ((d(3)*gamma(1+(1/d(4))))/(w2*mu_L2)) - 1.4;
c(6) = ((d(5)*gamma(1+(1/d(6))))/(w3*mu_L3)) - 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(3)*gamma(1+(1/d(3)))))/(d(3)*gamma(1+(1/d(4)))) - 0.2;
c(12) = (sqrt(((d(5))^2)*gamma(1+(2/d(6)))-d(5)*gamma(1+(1/d(5)))))/(d(5)*gamma(1+(1/d(6)))) - 0.2;
%
ceq(1) = wblcdf(d,mu_Sr1,sig_Sr1);
ceq(2) = wblcdf(d,mu_Sr2,sig_Sr2);
ceq(3) = wblcdf(d,mu_Sr3,sig_Sr3);
%
end
%

Thank you in advance

Answer 1

[Torsten]

F_S1 = wblpdf(d,mu_Sr1,sig_Sr1)
F_S2 = wblpdf(d,mu_Sr2,sig_Sr2)
F_S3 = wblpdf(d,mu_Sr3,sig_Sr3)

d is a vector with six values. Thus F_S1,F_S2 and F_S3 are vectors with six values each. Thus F_S1.*F_S2.*F_S3 is a vector with six (constant) values.

Defining

F_jpdf = @(d) F_S1.*F_S2.*F_S3;

thus makes no sense because F_S1.*F_S2.*F_S3 is just a constant vector (calculated with the vector for d transfered to "parameterfun" by "fmincon"), not a function that changes value when called with different d's.

It follows that

Rs = 1-integral(F_jpdf,0,1,'ArrayValued',true);

is simply

Rs = 1 - F_S1.*F_S2.*F_S3

which is also a vector with six values (and not a scalar as required).

I have the impression that you are far off with your code with respect to what you are really trying to do.

[Paul AGAMENNONE]

Hello @Torsten,

In reality, I'm trying to determine the probability of failure of the system, with components following Weibull distribution. Moreover, I don't know the values of mean and standard deviation of the components following Weibull distribution. This is what I call "d" in my code. I got the optimization model, and tried to code my problem, but I face some errors. You can see with the photos what I'm talking about.

Mean and standard deviation for Weibull distribution

Thank you for the help