Reliability based robust optimization using fmincon

I am trying to solve the following problem from a text:

minimize f=(x1.*x2.*cos(x1))+(x1.^2)-((1/4).*(x2.^2))-exp(x2)

subject to: g1=((x1-1).^2)+(x2.^2)-x1-6 g2=((3/7).*(x1.^2))-((1/10).*(x2))+((x2-1).^2)-5 -2.0<=x1<=2.0; -2.0<=x2<=2.0;

The design variables x1 and x2 have the normal distribution with standard deviation of 0.1 The objective function used for robust optimization is (0.5*mean(f))+(0.5*std(f))

Code using fmincon for robust optimization is:

main.m:

options = optimoptions('fmincon','Display','iter');
x0=[0.4,0.4];
[xopt,fval]=fmincon('objective1',x0,[],[],[],[],[-2,-2],[2,2],'constraint1',options)

objective1.m:

function fn=objective1(x)
x1=x(1)+normrnd(0,0.1,1e6,1);
x2=x(2)+normrnd(0,0.1,1e6,1);
f=(x1.*x2.*cos(x1))+(x1.^2)-((1/4).*(x2.^2))-exp(x2);% objective function
fn=(0.5*mean(f))+(0.5*std(f));
end

constraint1.m:

function [c,ceq]=constraint1(x)
c(1)=((x(1)-1).^2)+(x(2).^2)-x(1)-6; %inequality constraint 1
c(2)=((3/7).*(x(1).^2))-((1/10).*(x(2)))+((x(2)-1).^2)-5; %inequality constraint 2
ceq = []; 
end

I am not getting the robust value of x = [0.265, -0.593] given in the text.

Code using fmincon for reliability based robust optimization (target probability of failure being 0.1) is:

main.m:

options = optimoptions('fmincon','Display','iter');
x0=[0.4,0.4];
[xopt,fval]=fmincon('objective1',x0,[],[],[],[],[-2,-2],[2,2],'constraint1',options)

objective1.m:

function fn=objective1(x)
x1=x(1)+normrnd(0,0.1,1e6,1);
x2=x(2)+normrnd(0,0.1,1e6,1);
f=(x1.*x2.*cos(x1))+(x1.^2)-((1/4).*(x2.^2))-exp(x2);% objective function
fn=(0.5*mean(f))+(0.5*std(f));
end

constraint1.m:

function [c,ceq]=constraint1(x)

pt=0.1; %target probability of failure
x1=x(1)+normrnd(0,0.05,1e5,1);
x2=x(2)+normrnd(0,0.05,1e5,1);

g1=((x1-1).^2)+(x2.^2)-x1-6; %inequality constraint 1
pf1=mean(g1<=0); % Prob of failure 1
c(1)=pf1-pt; % checking for probability of failure of 1st constraint;

g2=((3/7).*(x1.^2))-((1/10).*(x2))+((x2-1).^2)-5; %inequality constraint 2
pf2=mean(g2<=0); % Prob of failure 2
c(2)=pf2-pt; % checking for probability of failure of 2nd constraint;

ceq = []; 
end

Is the code correct? Please help me to get the solution. Thanks.

Answer 1

[Matt J]

I'm not completely sure what you're trying to do, but it makes no sense to put any random number generation inside your objective function, as you are doing here,

x1=x(1)+normrnd(0,0.1,1e6,1);
x2=x(2)+normrnd(0,0.1,1e6,1);

The objective function has to be a deterministic function of the design variables x(i) if the iterative algorithm used by fmincon to minimize it is to have any hope of converging. When randomness is present, then the location of the minimum is also changing randomly, and the iterations of the algorithm end up trying to converge to a randomly moving target.

[Raghavi B]

Matt, thanks for your reply. I can understand it. But how do I consider x1 and x2 to be random following normal distribution then?

I tried making the random variable as global and kept it same for every iteration as shown in the code below. But this also did not work.

main.m:

options = optimoptions('fmincon','Display','iter');
global r;
r=normrnd(0,0.1,100,1); %random variable following noramal distribution
x0=[0.4,0.4];
[xopt,fval]=fmincon('objective1',x0,[],[],[],[],[-2,-2],[2,2],'constraint1',options)

objective1.m

function fn=objective1(x)
global r;
x1=x(1)+r;
x2=x(2)+r;
f=(x1.*x2.*cos(x1))+(x1.^2)-((1/4).*(x2.^2))-exp(x2);% objective function
fn=(0.5*mean(f))+(0.5*std(f));
end

constraint1.m:

function [c,ceq]=constraints(x)
c(1)=((x(1)-1).^2)+(x(2).^2)-x(1)-6; %inequality constraint 1
c(2)=((3/7).*(x(1).^2))-((1/10).*(x(2)))+((x(2)-1).^2)-5; %inequality constraint 2
ceq = []; 
end

[Matt J]

But how do I consider x1 and x2 to be random following normal distribution then?

The only thing that makes sense to me is that you are supposed to be first simulating noisy measurements of f

   xsim = [0.265, -0.593];

   x1=xsim(1)+normrnd(0,0.1,1e6,1);
   x2=xsim(2)+normrnd(0,0.1,1e6,1);

   fmeas=(x1.*x2.*cos(x1))+(x1.^2)-((1/4).*(x2.^2))-exp(x2);

And then, your objective function is really supposed to be

function fn=objective(x, fmeas)

    f=(x1.*x2.*cos(x1))+(x1.^2)-((1/4).*(x2.^2))-exp(x2);

    Error=(fmeas-f);

    fn=(0.5*mean(Error).^2)+(0.5*std(Error).^2);

end