Multiobjective optimization with constraints

Hi All,

I am trying to solve a game theoretic scenario where I have two players that each have three decision variables. They each have their own individual utility functions with their payoffs depending on what the other one does. As a result, I have to optimize these two equations simulataneously. I've looked into using the gamultiobj tool, and I've tried to follow an example from the MATLAB website, but obviously my scenario is a little bit more complicated. They just had one variable, whereas I have six variables in total, with additional constraints of all variables being nonnegative, and x(5) and x(6) being <=24. I'm not exactly sure where I went wrong, especially since my situation differs quite a bit from the example, so i would greatly appreciate any and all help. I've included my code below (the example, then mine). Thanks in advance!

% Example version
fitnessfcn = @(x)[sin(x),cos(x)];
nvars = 1;
lb = 0;
ub = 2*pi;
rng default % for reproducibility
x = gamultiobj(fitnessfcn,nvars,[],[],[],[],lb,ub)
plot(sin(x),cos(x),'r*')
xlabel('sin(x)')
ylabel('cos(x)')
title('Pareto Front')
legend('Pareto front')

%My version
funx=-0.5*x(1)^2+100*log(x(1)+x(2)-14)-x(5)/2+0.9*(-0.5*(1-x(5)/40)*x(3)^2+100*log(x(3)+x(4)-14))
funy=-0.5*x(2)^2+100*log(x(1)+x(2)-14)-x(6)/2+0.9*(-0.5*(1-x(5)/40)*x(4)^2+100*log(x(3)+x(4)-14))
fitnessfcn = @(x)[funx,funy];
nvars = 2;
lb = [0,0,0,0,0,0];
ub = [Inf,Inf,Inf,Inf,24,24];
rng default % for reproducibility
x = gamultiobj(fitnessfcn,nvars,[],[],[],[],lb,ub)
plot(funx,funy,'r*')
xlabel('funx')
ylabel('funy')
title('Pareto Front')
legend('Pareto front')

Answer 1

[Alan Weiss]

Your input data is inconsistent:

nvars = 2;
lb = [0,0,0,0,0,0];
ub = [Inf,Inf,Inf,Inf,24,24];

Your bounds are on your decision variables, not your objective functions. The nvars argument refers to exactly this, so you should have nvars = 6.

Your other errors relate to how to call functions and use ther resulting data. Your fitness functions need initial @x arguments:

funx = @(x)-0.5*x(1)^2+100*log(x(1)+x(2)-14)-x(5)/2+0.9*(-0.5*(1-x(5)/40)*x(3)^2+100*log(x(3)+x(4)-14));
funy = @(x)-0.5*x(2)^2+100*log(x(1)+x(2)-14)-x(6)/2+0.9*(-0.5*(1-x(5)/40)*x(4)^2+100*log(x(3)+x(4)-14));
fitnessfcn = @(x)[funx(x),funy(x)];

Your plot needs to pass the data to the objective functions:

plot(funx(x),funy(x),'r*')

When I ran this code I got other errors relating to the function being complex-valued. You also need to include linear constraints that keep the function from having complex values. But this should get you started.

Alan Weiss

MATLAB mathematical toolbox documentation

[Sabrina Chui]

Hi Alan,

Thanks, I input the corrections to the code you made and saw the error about the function being complex-valued. I didn't really think that my functions would need linear constraints, but when I thought about it, I guess having the log(x(1)+x(2)-14) might be a source of trouble since if x(1)+x(2)<=14, then it would be undefined. Same for x(3)+x(4). So I tried to put these two linear inequality constraints into the function, but I'm guessing my syntax wasn't correct?

funx = @(x)-0.5*x(1)^2+100*log(x(1)+x(2)-14)-x(5)/2+0.9*(-0.5*(1-x(5)/40)*x(3)^2+100*log(x(3)+x(4)-14));
funy = @(x)-0.5*x(2)^2+100*log(x(1)+x(2)-14)-x(6)/2+0.9*(-0.5*(1-x(5)/40)*x(4)^2+100*log(x(3)+x(4)-14));
fitnessfcn = @(x)[funx(x),funy(x)];
nvars = 6;
lb = [0,0,0,0,0,0];
ub = [Inf,Inf,Inf,Inf,24,24];
A=[A;-1 -1 0 0 0 0];
A=[A;0 0 -1 -1 0 0];
b=[b;-14.1 -14.1]

rng default % for reproducibility
x = gamultiobj(fitnessfcn,nvars,A,b,[],[],lb,ub)
plot(funx(x),funy(x),'r*')
xlabel('funx')
ylabel('funy')
title('Pareto Front')
legend('Pareto front')

I also thought I might eliminate the problem of complex values if just increased the lower bound for x(1),x(2),x(3) and x(4) to be 7, which did work in giving 70 solutions, but I'm not quite sure if they are valid, especially since these solutions change when I make the lower bounds 5, 4, or 8.

[Alan Weiss]

I would write the linear constraints this way:

A = [-1 -1 0 0 0 0
     0 0 -1 -1 0 0];
b = [-14.1;-14.1];

You might also want to set options to use the gaplotpareto plot function in order to have the plot drawn automatically while your optimization is underway.

Alan Weiss

MATLAB mathematical toolbox documentation

[Sabrina Chui]

Thanks for the suggestion Alan. In addition to your suggestion, I also took from the method shown here and set my problem up in the optimizer app, and ended up with a nice pareto front and some solutions. Since my example was quite complicated( it involved 6 variables, linear inequalities and upper and lower bounds), I wanted to test this method on a different simpler example that I already know what the answer should be. However, I ran into some issues and was wondering if you could tell me where I went wrong? This time I followed a similar tutorial, from the same youtuber. I have my problem set up like this:

function Output = objective_function_test(Input)
e1 = Input(1);
e2 = Input(2);

f1 = -(e1^(1/3)+(10-e1-e2)^(1/3));
f2 = -(e2^(1/3)+(10-e1-e2)^(1/3));

Output = [f1 f2]; 

and then put it into the optimization solver using the gamultiob option like this:

I end up with errors, which dont allow the pareto front to be found. The erro messages shown above are from me changing up the name of the fitness function, since that was the only thing different from the tutorial video. I'm not sure why this isn't working, because my first try with 6 variales worked. And it seems like in the second tutorial 2 variables aren't a problem. So is there something wrong with my set up? Thanks again. I just really want to be sure about this optimizer, and want to test it with an example that I already know the solutions for.

[Alan Weiss]

Please don't change the default tolerances unless you know exactly what you are doing. Leaving all options at their default values is the best way to ensure that your issues are not related to option choices, but are only related to your setup. In particular, your function tolerance and constraint tolerance choices are poor.

I don't know where you put your objective function file or what you called it. It seems that MATLAB cannot find it. I cannot help with that. Just put it somewhere MATLAB can find it, like the current folder, and name it properly.

Your objective function returns complex values for negative arguments or for arguments whose sum exceeds 10. To keep the sum bounded away from 10, put upper bounds of [10 10] and also set A = [1 1] and b = 10.

Alan Weiss

MATLAB mathematical toolbox documentation

[Sabrina Chui]

Hi Alan, I changed all the options back to default, changed the upper boundaries and that seemed to work. When I look through my set of solutions for e1 and e2, none of them even get close to the solutions that I worked out by hand. For the two objective functions I provided, f1 and f2, the nash equilibrium should be e1=3.3333 and e2=3.3333. When I ran this example, I think I was only expecting one solution on the pareto front, but it was confusing to see a long set of solutions. Is it because of the way the algorithm works, if I increase the iterations will I get closer to what I expected? Thank you!

[Alan Weiss]

For the problem as you have defined it, the solution has e1 going from 0 to 5, and e2 going from 5 to 0 as follows (in one run, these numbers are stochastic):

solution =

    0.0001    4.9999
    0.0001    4.9999
    0.0825    4.9803
    0.1373    4.5127
    0.0062    5.0064
    5.0000         0
    4.3021    0.0542
    4.9261    0.0507
    0.8771    3.0301
    3.8961    0.1244
    0.5075    3.3599
    0.9395    3.9216
    4.4089    0.3093
    0.5050    4.5467
    4.4039    0.6124
    2.2507    2.8033
    4.9980    0.0012
    0.3360    4.0666

There is a tradeoff between the objective functions; here are the associated values:

objectiveValue =

   -1.7589   -3.4199
   -1.7589   -3.4199
   -2.1381   -3.4105
   -2.2649   -3.4015
   -1.8919   -3.4192
   -3.4200   -1.7100
   -3.4068   -2.1589
   -3.4141   -2.0828
   -2.7837   -3.2735
   -3.3886   -2.3142
   -2.6281   -3.3282
   -2.7051   -3.3026
   -3.3813   -2.4178
   -2.5004   -3.3607
   -3.3472   -2.5573
   -3.0143   -3.1138
   -3.4198   -1.8165
   -2.4708   -3.3717

So I think that it is your interpretation of the solution that is at fault, not the solution.

Note that your expected solution is not even on the Pareto curve (there are better points). Your expected solution has value (-2.9876,-2.9876), but there are points (such as the third to last) with both values smaller than -3.

Alan Weiss

MATLAB mathematical toolbox documentation

[Sabrina Chui]

Thanks for clariying this Alan! I'll have to do a little bit more reading on the matter, but this is a great start!