how to tweak fmincon into yielding Global Minimum

I am using the fmincon function to derive the elements of a mxn matrix (X) which will maximize the following objective function subject to the non-negativity constraint and sum(xij)<100,

{f(x(1,1))+(f(x(1,2))+...+f(x(1,n))} + Beta*{f(x(2,1))+(f(x(2,2))+...+f(x(2,n))}+...+(Beta^(m-1))*{f(x(m,1))+(f(x(m,2))+...+f(x(m,n))}

where f(xij)=zeta*[((alpha*pf*(xij^delta))-(pw*xij))^gamma]

The code I used to do this was (fmincon was used, as the objective function is non-linear);

    function optim = optim(m,n)
    A= ones(1,m*n);
    b = 100;
    z = zeros(m,n);
    in = inf(m,n);
    X = ones(m,n);

    [x, bestval] = fmincon(@myfun2,X,A,b,[],[],z,in,[])

      function f = myfun2(x)
      Alpha = 5;
      %kappa = 5;
      zeta = 5;
      beta = 0.90909;
      delta = 0.4;
      gamma = 0.4;
      pf = 1;
      pw = 1;

      for i= 1:m
         f=0;
         sum(i)=0;
            for j=1:n
            sum(i) = sum(i) +((beta^(i-1))*(-1)*(zeta)*(((Alpha*pf.*((x(i,j)).^delta))- 
                     (pw.*x(i,j)))^gamma));

             end
         f = f+sum(i)
       end

     end
    end

When the code was run for a 5x5 matrix (optim(5,5)), the resulting solution was x =

    1.5439    1.5439    1.5439    1.5439    1.5439
    1.5439    1.5439    1.5439    1.5439    1.5439
    1.5439    1.5439    1.5439    1.5439    1.5439
    1.5439    1.5439    1.5439    1.5439    1.5439
    3.1748    3.1748    3.1748    3.1748    3.1748

    bestval =

    -31.8780

But this is not a global minimum; as the marginal conditions specify at the global minimum we would have; (for 1st row) f'(x11)=...=f'(x1n) and so on for each row.Also we would have for each column f'(x11)=beta*f'(x21)=...=(beta^m-1)*f'(xm1) (for the first column and so on). None of these conditions are satisfied by the resulting matrix. I have looked at the the related questions in stack overflow and also the documentation and I have no inkling of how to get better results. Is there a problem with the code? Can the code be tweaked to get better results?

Can I make the marginal conditions the stopping conditions and how could i go about doing that. or Can I use the Jacobian in some way? Any help would be appreciated. Thank You.

Answer 1

[Alan Weiss]

I think that the following two lines are out of order:

 for i= 1:m
    f=0;

The way you have it, f is reset for the first m - 1 calculations, and the only thing that affects the result is the line when i = m.

It is not at all clear to me why you have so many different values of x(i,j). It seems to me that the result depends only on the value of sum_j f(x(i,j)). So maybe you should make that the variable, and reduce on level of summation. And it is not clear to me why you have loops rather than a matrix multiply. But these are different issues, and it is possible that I misunderstand what you are really trying to do.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 2

[Matt J]

When I re-implement the code as below, with matrix multiplication as Alan suggested, and also essentially de-activating the TolX and TolFun stopping criteria, I get the following solution,

    x =

        0.6350    0.6350    0.6350    0.6350    0.6350
        0.6350    0.6350    0.6350    0.6350    0.6350
        0.6350    0.6350    0.6350    0.6350    0.6350
        0.6350    0.6350    0.6350    0.6350    0.6350
        0.6350    0.6350    0.6350    0.6350    0.6350

It is clearly the correct solution, since xij=0.635 is the location where every individual f(xij) term attains its unconstrained maximum. Since the constraints are inactive at this solution, it solves the constrained problem as well.

   function out = optim(m,n)
      A= ones(1,m*n);
      b = 100;
      z = zeros(m,n);
      in = z;
      X = ones(m,n);
      [x, bestval,exitflag,s,l,g] = ...
        fmincon(@myfun2,X,A,b,[],[],z,z+2.92,[],...
                optimset('TolFun',1e-20,'TolX',1e-20));

      x, bestval

                  function val = myfun2(x)

                          Alpha = 5;
                          %kappa = 5;
                          zeta = 5;
                          beta = 0.90909;
                          delta = 0.4;
                          gamma = 0.4;
                          pf = 1;
                          pw = 1;

                          betas=beta.^(0:m-1);

                           f= ( ((Alpha*pf)*x).^delta- pw*x ).^gamma; 
                           val=-zeta*sum(betas*f);

                  end
       end

   end

[Matt J]

A few additional remarks,

1. Do not use 'optim' for two different things like in the line "function optim = optim(m,n)"
2. It is in general a bad idea (or at least undesirable) to use 'sum' as a variable, since you then lose the MATLAB builtin sum() command.
3. Due to the symmetry of the problem, it is clear that x will always have a solution with identical columns. You may as well rewrite the problem in terms of just one of them.
4. You need additional constraints ensuring that ((alpha*pf*(xij^delta))-(pw*xij)) > 0. The objective function can be complex-valued and/or non-differentiable outside that region.

Answer 3

[Uday]

Thank You Alan and Matt for your responses and comments. To explain the problem better, think of n number of people using x over m time periods. The value they get from use is f(). What I am trying to do is maximize the sum of the total value they obtain from their use of x over all time i.e. sum of f(xij). But at this time the value gained from future use is discounted by a factor of B^i-1 at every time period i. The only constraints are that the use cannot be negative (therefore the lower bound) and x is in limited supply (provisionally at b=100). I cannot take f(xij) or sum(f(xij) as the variable being adjusted because I would also like to derive how much of x each individual j would use at each time i.

Thanks for pointing out my errors with optim, sum and f. The reason i cannot use an out and out matrix multiplication is that though currently f is the same for all js; I plan to tweak the parameters (alpha, delta, pf and gamma) as well as the functional form f() at a later part of time so that they are heterogeneous (but random) across the js. At which point a simple matrix multiplication would not help me much. Therefore I would like the method to be generalizable to a much larger matrix with heterogeneous parameters and functional form.

Also Matt regarding the results; given that the functional form is the same (for now) across js; the value of xij in each row should remain constant (identical columns as you pointed out). But as the contribution of each row to the total objective function gets discounted by a factor of B^i-1,the xij across all the rows cannot have the same value.Logically the xij values should decrease over time(i).

I will tweak the code a bit based on your suggestions and get back to you for your suggestions. Thank You.

[Matt J]

But as the contribution of each row to the total objective function gets discounted by a factor of B^i-1,the xij across all the rows cannot have the same value.Logically the xij values should decrease over time(i).

This is only true if the constraint sum(x(:))<=b is strong enough to have an impact on the solution. If applying or removing the constraint does not change the solution (as in your example, because b is so large) then the discounting has no impact at all.

If the constraints don't matter, then the solution can be found at the unconstrained maximum, which also happens to lie where each term beta^(i-1)*f(xij) term is maximized independently of the other terms. But the presence of the discount factor beta^(i-1) has no impact when you maximize each term independently. It is just a scale factor that doesn't change the location of the maximizing xij.

I plan to tweak the parameters (alpha, delta, pf and gamma) as well as the functional form f() at a later part of time so that they are heterogeneous (but random) across the js. At which point a simple matrix multiplication would not help me much.

Doesn't sound like it would make a difference to me. You still have a weighted sum of terms. A weighted sum of anything can always be written as a matrix multiplication.