Solving an optimization problem both graphically and numerically

Here are the problems, along with the solutions.

Problem 1

--

[x1, x2] = meshgrid(-3:.1:3);
z = -x1 .^ 2 - x2;
i = find(x1 .^ 2 + x2 .^ 2 > 9); z(i) = NaN;
i = find(x1 + x2 > 1); z(i) = NaN; 
surf(x1, x2, z); shading interp

Somehow from the graph we can see that solutions are x_1 = 0, x_2 = −3, and the maximum value of the function is 3.

Problem 2

--

[x1, x2] = meshgrid(0:0.02:1, 1:0.02:2);
z = x1 .^ 3 + x2 .^2 + 4 * x1 + 4;
ii = find(x1 - x2 + 2 < 0); z(ii) = NaN;
ii = find(-x1 .^ 2 + x2 - 1 < 0); z(ii) = NaN;
ii = find(x1 < 0); z(ii) = NaN;
ii = find(x2 < 0); z(ii) = NaN;
surf(x1, x2, z); shading interp

Somehow we can see that x_1 = 0, x_2 = 1.

Numerical solution:

function [c, ce] = exc6f1(x)
    ce = [];
    c = [x(1)^2 - x(2) + 1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = @(x) x(1) ^ 3 + x(2) ^ 2 + 4 * x(1) + 4;
A = [-1 1]; B = 2; Aeq = []; Beq = []; xm = [0;0];
x = fmincon(f, [0; 1], A, B, Aeq, Beq, xm, [], 'exc6f1')

My questions are:

1. How to determine arguments for the meshgrid function? For the first problem, my book says "from the given constraints, the initial square region [−3, 3] can be selected and grids can be made". Can you explain it step-by-step (for both problems)?

2. How do we read solutions from the graph?

3. How do we determine initial value for the numerical solution? I saw several different functions having this kind of argument (e.g. fsolve, ode45), but I don't know what to put here.

Answer 1

[Matt J]

1. How to determine arguments for the meshgrid function? For the first problem, my book says "from the given constraints, the initial square region [−3, 3] can be selected and grids can be made". Can you explain it step-by-step (for both problems)?

In problem 1, the first constraint says that (x1,x2) are bounded to a circle of radius 3. So by choosing a square region that covers this circle, you are certain to see at least the entire feasible set in the plot. In problem 2, you can analyze the feasible region more easily by re-arranging the constraints as

   1+x1^2 <= x2 <= x1+2

In other words, x2 is lower bounded by a quadratic function of x1 and upper bounded by a linear function of x1. If you plot these as a function of x1, or solve the equation 1+x1^2=x1+2 you will see that 0<=x1<=1 and 1<=x2<=2 are the only intervals where the linear bound lies above the quadratic bound, so choosing the meshgrid to cover these intervals again ensures that the plot will view the entire feasible set.

2. How do we read solutions from the graph?

The surface will reach a maximum or minimum at the solution, so just observe the highest peak or lowest valley on the surface plot. Keep in mind though that this is only an approximate solution. The plot only shows the values of the objective function at the discrete points that you've sampled in the meshgrid. In general, you can't expect the sample points to hit the true solution exactly, except in simple exercises like these that have been set up to have nice integer solutions.

3. How do we determine initial value for the numerical solution?

The idea is to initialize as close to the solution as you can. In general, there is no systematic way to choose this point. Insights specific to the problem must be used. However, in this case, your surface plots let you see approximately where the solution lies, so you can initialize somewhere close to that.

[Nikola Kostic]

1. Can MATLAB help me somehow here, or I have to solve this on paper?

2. As far I can see, highest peek is colored in bright yellow (and lowest valley in dark blue). My guess now is that if I rotate the graph to X-Y view I should be able to easily read the solution, I just check the coordinates of the brightest/darkest point? And I can easily see the function value by rotating to X-Z view? Something tells me that graphical method isn't very useful when there are more that two variables. Can it only be used for two?

[Matt J]

1. Can MATLAB help me somehow here, or I have to solve this on paper?

What do you mean? You have already used fmincon to obtain the solution.

I just check the coordinates of the brightest/darkest point? And I can easily see the function value by rotating to X-Z view?

Yes, or you can use the datatip tool to probe points on the surface.

[Nikola Kostic]

I'm talking about meshgird function arguments. You are suggesting that I can determine them by solving the system of inequalities?

[Matt J]

In general, yes, you have to do some sort of problem-specific pre-analysis to determine a rectangle containing the feasible set, although Matlab can help with that pre-analysis, e.g., with plotting.

An exception would be in the case when all your constraints are linear. In that case, and if the feasible region is bounded, then the feasible region will be a polygon whose vertices you can get numerically using lcon2vert (Download). A bounding rectangle can then be obtained by looking at the maximum and minimum vertex coordinates.

[Nikola Kostic]

For the second problem, I wouldn't really say that 0 <= x1 <= 1 and 1 <= x2 <= 2. I would more likely say that 0 <= x1 <= 1.618 and 1 <= x2 <= 3.618?

[Matt J]

Yes, I think you're right. But you can also see from the objective function that it is monotonic in x2 and so if you only care about seeing the minimum, the sub-region they gave is enough.

Answer 2

[Matt J]

2. How do we read solutions from the graph?

Note that you don't actually have to read the solution from the graph. You can also compute it from x1,x2, and z.

 optval=max(z(~isnan(z))); %or min()

 optlocs=(z==optval);

 solutions=[x1(optlocs), x2(optlocs)];