Maximizing a function in Matlab (use of fsolve and derivatives)

Question

I would like to find the values for vector p (i.e p(1),p(2),p(3),...) which maximizes the function A(p).

I am using MATLAB to do that and I found fsolve which I thought it could help me. So I made function A:

function A = myfun(p)

    R  = 0.1; 
    u1 = 500;
    u2 = 400;
    u3 = 300;

    A = ( (p(1)+p(2)+p(3)) * (1/u1+1/u2+1/u3)) * ...
        (1 + R*(p(1)^2+p(2)^2+p(3)^2) * (1/u1+1/u2+1/u3) );

And then I need to solve a system of equations which will be:

diff(A,p(1))==0

diff(A,p(2))==0

diff(A,p(3))==0

where the resulting p vector will be the solution to my problem.

How could fsolve solve this system of equations (p0=[1 1 1]) ?

Answer 1

Here's an example of how to do it using the symbolic toolbox:

%// (there's probably a way to generalize this as well)
dAdp = cellstr(char(...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p1')) ';'],...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p2')) ';'],...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p3')) ';']))

%// convert to proper vector equation
dAdp = regexprep(dAdp, 'p([0-9])', 'p\($1\)');

%// convert to function handle
F = str2func( strcat('@(p) [', dAdp{:}, ']') );

but I would not recommend doing it this way (it's completely unreadable and bug-prone).

You could write a proper function and evaluate dAdp above using the symbolic toolbox, but I also would not recommend doing it that way (it's just very slow).

I'm a numerical guy, because I trust computers only with what they are built for: calculations. Derivations I do on paper, except for the excessively simple, yet long and tedious ones (you could argue that this is such a case, I still like doing it on paper, because I also like to exercise my brain :).

I propose you do this, and/or checking yourself continuously using the symbolic toolbox. IMHO, you should use it as an aid, not as the main engine.

So, here we go. Here's your function again, this time in a different form, with slightly more brain applied:

function A = myfun(p)    
    R = 0.1; 
    u = [500; 400; 300];     
    C = sum(1./u);
    B = C * sum(p) * (1 + R*C*sum(p.^2));

So, you want to solve ∇A(p) = 0. Best way is to apply more brain. You may verify that the vectorial derivative is equal to:

function F = Aprime(p)
    R = 0.1;
    u = [500; 400; 300];
    C = sum(1./u);
    F = C*( C*R*sum(p.^2) + 2*C*R*p*sum(p) + 1 );

which you can either solve with more brain:

C·( CR·Σp² + 2CR·p·Σp + 1 ) = 0
                Σp² + 2p·Σp = -1/(CR)

vector == scalar: that means all elements in p are equal.
Substituting q = p1 = p2 = p3, then

 3q² + 2q·3q = -1/(CR)
         9q² = -1/(CR)

⇒ q = ⅓·√(-1/(CR))

(which indicates trouble), or in fsolve like so:

fsolve(@Aprime, [1 1 1])

which will immediately result in

fsolve completed because the vector of function values at the initial point is near zero as measured by the default value of the function tolerance, and the problem appears regular as measured by the gradient.

which indeed indicates trouble.

Since you have indicated that you're interested in the maximum, and the function does not appear to have stationary points, the only conclusion is that the function has no maximum and no minimum, unless you impose bounds on p. Indeed, if you reduce the dimensionality to 2, and make a plot:

R = 0.1;
u = [500; 400; 300];
C = sum(1./u);
B = @(p) C * sum(p) .* (1 + R*C*sum(p.^2));

[p1,p2] = meshgrid(-10:0.1:10);
surf(p1,p2,reshape(B([p1(:) p2(:)].'), size(p1)), 'edgecolor', 'none')

...there is no extremum.