How to solve this optimization problem (with matlab optimization toolbox)?

Hi

How to solve this problem?

[John D'Errico]

Since basic mathematical reasoning gives the answer, you don't solve it with the optimization toolbox. Think of it another way. Is the answer unique?

[Brede roberson]

Thanks

How do you think it can be solved (minimized) in a different way?(n=4,l=d=h=2)

[Torsten]

https://www.wolframalpha.com/input/?i=FindMinimum%5B%7BSqrt%5B1%2B(y%2Fx)%5E2%2B(z%2Fx)%5E2%5D%2BSqrt%5B1%2B((20-y)%2F(16-x))%5E2%2B((14-z)%2F(16-x))%5E2%5D+%26%26+x%3E%3D0+%26%26+x%3C%3D16+%26%26+y%3E%3D0+%26%26+y%3C%3D20+%26%26+z%3E%3D14+%26%26+z%3C%3D14%7D,%7Bx,y,z%7D%5D

Best wishes

Torsten.

Answer 1

[John D'Errico]

No other answer yet? This is probably homework, or some sort of contest puzzle.

The optimum does not exist under the existing constraints. Essentially, there is no optimum, because that optimum lies at a point where a singularity exists, and outside of the boundaries of a partially open set. You an go arbitrarily close to that point, but you cannot achieve the true optimum.

As such, you cannot pose this as a problem for the optimization toolbox, because it does not allow strict inequality constraints. In fact, this is not really even a MATLAB question, because the optimum lies at a point of singularity. Instead, some thought experiments allow you to identify the solution.

Inspection is often a good way to think these things out on a problem like this. Try a few examples, so build some intuition about what is happening.

The problem is essentially to find an array XYZ of size nx3. Thus x is the first column, etc.

We are given constraints, but only a few. We are told that all of the elements of xyz are STRICTLY positive. zero is not allowed. Also,

sum(x) = L
sum(y) = d
sum(z) = h

Or, we can write it in pseudo-MATLAB form, given a vector LDH of size 1x3, as

sum(xyz,1) == LDH

or as

ones(1,n)*xyz == LDH

Now, we want to minimize a form that could be written in MATLAB as:

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))

At least, that applies as long as you have a current MATLAB release. So lets try it out, and see what happens for some numbers. Again, it helps to build intuition here. Large mathematical leaps tend to be confusing to someone who does not see where it is going.

n = 5;
LDH = [2 3 4]
LDH =
     2     3     4

xyz = rand(n,3);
xyz = xyz./sum(xyz,1).*LDH
xyz =
      0.64278      0.29857      0.75635
      0.19223      0.73254      0.71042
      0.73354      0.56279       1.1853
      0.27627      0.41816      0.36939
      0.15519      0.98794      0.97854

sum(xyz,1)
ans =
            2            3            4

Ok, so I picked some sets of numbers for xyz, scaled them to sum correctly. What is our objective function for that set?

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       20.333

It is almost always true that problems like this have a solution that lies at a boundary of some sort. If a minimum exists, since none of the rows of xyz are any different from the rest in context of the constraint or for the objective, lets consider the case where all the rows of xyz are the same. In fact, this candidate solution does lie at a boundary.

xyz = ones(n,1)*(LDH./n)
xyz =
          0.4          0.6          0.8
          0.4          0.6          0.8
          0.4          0.6          0.8
          0.4          0.6          0.8
          0.4          0.6          0.8

sum(xyz,1)
ans =
     2     3     4

Clearly, it satisfies the requirements. And it yields a significantly lower objective function than the random choice I tried.

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       13.463

Does this mean it is the solution? Well, no. Look at the objective carefully. Move the x_i into the sqrt. Remember that x is both positive and greater than zero. We can write it as a sum of squared ratios.

sum( sqrt( 1 + (y_i/x_i)^2 + (z_i/x_i)^2 ) )

How can we make this as small as possible? Hmm. If we start at the point where x,y,z were all constant, then head towards the array that I will describe as

xyz = [L-(n-1)*eps), D-(n-1)*realmin, H-(n-1)*realmin;
       eps, realmin, realmin;
       ...
       eps, realmin, realmin]

So for all rows of xyz, make y and z zero, or as small as possible. Now, eps is avery small number in MATLAB, but realmin is essentially immeasurably small compared to eps.

eps
ans =
   2.2204e-16

realmin
ans =
  2.2251e-308

I'll create xyz as

xyz = [LDH - (n-1)*[eps realmin realmin];ones(n-1,1)*[eps realmin realmin]]
xyz =
            2            3            4
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308

You can see that it clearly satisfies the requirements.

sum(xyz,1)
ans =
     2     3     4

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       6.6926

As well, it has a significantly smaller objective function than the previous attempt.

Note that I could have permuted the rows of xyz in any of factorial(n) ways and gotten the same value for the objective.

Could I have moved even closer to the boundary whee y(2:n) would be zero? In theory, yes. Go as close as you want to zero there, as long as you don't also divide by zero. You want to make x(1) as large as possible, so that (y(1)/x(1))^2 is as small as possible. You cannot reach that point though, because you will run into a divide by zero.

Could you have made two of the rows of xyz non-zero, whicle leaving the remainder as TINY numbers in the same way?

xyz = [ones(2,1)*[LDH/2 - (n-2)*[eps realmin realmin]];ones(n-2,1)*[eps realmin realmin]]
xyz =
            1          1.5            2
            1          1.5            2
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       8.3852

As you can see in the example, this is considerably worse.

Is this the solution? Well, since your question is probably a homework problem, no I won't prove that this is the global solution. And I refuse to do homework. At most, all I've done here is point you in a direction that should give you the optimum, in case this is not homework.

Again, you need to recognize that MATLAB cannot in fact solve this problem using strictly numerical methods, because the optimization toolbox cannot accept strict inequality constraints. The feasible set of candidate solutions forms a partially open set. You can keep moving towards the boundary, but you will never reach it. Infinity is a long way out, but you will never get there. The same thing happens as you approach zero. In fact, I would argue this is why Wolfram Alpha also fails to provide a solution, saying that no solution exists. In fact, the solution lies just over the horizon, but you can never reach that particular horizon.

Finally, are the actual values of L,D,H or n material to the problem? My intuition tells me that they are not, that you could have chosen other values, yet the solution will still lie as unattainably in that corner of the solution space that you cannot reach. Proof of all this is left to the student of course, since I refuse to do homework.

[John D'Errico]

I just noticed the comment with a specific set for L,D,H,n. The same set of conclusions apply.

LDH = [2 2 2];
n = 4;

xyz = rand(n,3);
xyz= xyz./sum(xyz,1).*ldh
xyz =
      0.84791     0.046955      0.22018
      0.42798      0.46198      0.96409
      0.63877      0.67815      0.79556
     0.085333      0.81292     0.020173

sum(xyz,1)
ans =
     2     2     2

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       15.225

xyz = ones(n,1)*(LDH./n)
xyz =
          0.5          0.5          0.5
          0.5          0.5          0.5
          0.5          0.5          0.5
          0.5          0.5          0.5

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       6.9282

But a better solution lies near the corner that we cannot get all the way into:

xyz = [LDH - (n-1)*[eps realmin realmin];ones(n-1,1)*[eps realmin realmin]]
xyz =
            2            2            2
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308
   2.2204e-16  2.2251e-308  2.2251e-308

sum(sqrt(sum(xyz.^2,2))./xyz(:,1))
ans =
       4.7321

Again, I won't go any further.