Optimize a 2 variable function

I'm trying for a very long time to optimize this function:

V(a,b)=(26*(2000*cos((pi*a)/180) + (sin((pi*(a + 45))/180)*((2600000*cos((pi*(a + b))/180)*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180)))/(cos((pi*a)/180)^2*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180))) + (2600000*cos((pi*a)/180)*sin((pi*(a + 45))/180))/(cos((pi*(a + b))/180)*sin((pi*b)/180)*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180)))))/((1300*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180))^2)/cos((pi*a)/180)^3 + (1300*sin((pi*(a + 45))/180)^2)/(cos((pi*(a + b))/180)*sin((pi*b)/180)^2) + 8085710131136235/4398046511104)))/(3*cos((pi*(a + b))/180)*sin((pi*b)/180)) - (26*((2^(1/2)*((2600000*cos((pi*(a + b))/180)*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180)))/(cos((pi*a)/180)^2*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180))) + (2600000*cos((pi*a)/180)*sin((pi*(a + 45))/180))/(cos((pi*(a + b))/180)*sin((pi*b)/180)*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180)))))/(2*((1300*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180))^2)/cos((pi*a)/180)^3 + (1300*sin((pi*(a + 45))/180)^2)/(cos((pi*(a + b))/180)*sin((pi*b)/180)^2) + 8085710131136235/4398046511104)) + (cos((pi*(a + b))/180)*(2000*cos((pi*a)/180) + (sin((pi*(a + 45))/180)*((2600000*cos((pi*(a + b))/180)*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180)))/(cos((pi*a)/180)^2*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180))) + (2600000*cos((pi*a)/180)*sin((pi*(a + 45))/180))/(cos((pi*(a + b))/180)*sin((pi*b)/180)*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180)))))/((1300*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180))^2)/cos((pi*a)/180)^3 + (1300*sin((pi*(a + 45))/180)^2)/(cos((pi*(a + b))/180)*sin((pi*b)/180)^2) + 8085710131136235/4398046511104)))/sin((pi*b)/180)))/(3*cos((pi*a)/180)^2) - (539047342075749*((2600000*cos((pi*(a + b))/180)*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180)))/(cos((pi*a)/180)^2*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180))) + (2600000*cos((pi*a)/180)*sin((pi*(a + 45))/180))/(cos((pi*(a + b))/180)*sin((pi*b)/180)*(cos((pi*(a + b))/180)*sin((pi*a)/180) - sin((pi*(a + b))/180)*cos((pi*a)/180)))))/(43980465111040*((1300*(2^(1/2)/2 + (cos((pi*(a + b))/180)*sin((pi*(a + 45))/180))/sin((pi*b)/180))^2)/cos((pi*a)/180)^3 + (1300*sin((pi*(a + 45))/180)^2)/(cos((pi*(a + b))/180)*sin((pi*b)/180)^2) + 8085710131136235/4398046511104))

The variables i need to solve for are a and b which represent angles. I want the GLOBAL MAXIMUM of the function. The solutions i need must satisfy the following conditions: -the derivative of V with respect to a in the points determined as global maximum must be equal to 0; -likewise for the derivative of V with respect to b; -the sum of a and b must be less then 90 degrees(a+b<90);

I tried with Optimization toolbox but i have no ideea what I'm doing. Thank you very much

Answer 1

[Star Strider]

I needed to parameterise your V equation for it to work with the Optimization Toolbox routines:

% % % p(1) = a,    p(2) = b;
V = @(p) -(26*(2000*cos((pi*p(1))/180) + (sin((pi*(p(1) + 45))/180)*((2600000*cos((pi*(p(1) + p(2)))/180)*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180)))/(cos((pi*p(1))/180)^2*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180))) + (2600000*cos((pi*p(1))/180)*sin((pi*(p(1) + 45))/180))/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180)))))/((1300*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180))^2)/cos((pi*p(1))/180)^3 + (1300*sin((pi*(p(1) + 45))/180)^2)/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)^2) + 8085710131136235/4398046511104)))/(3*cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)) - (26*((2^(1/2)*((2600000*cos((pi*(p(1) + p(2)))/180)*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180)))/(cos((pi*p(1))/180)^2*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180))) + (2600000*cos((pi*p(1))/180)*sin((pi*(p(1) + 45))/180))/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180)))))/(2*((1300*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180))^2)/cos((pi*p(1))/180)^3 + (1300*sin((pi*(p(1) + 45))/180)^2)/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)^2) + 8085710131136235/4398046511104)) + (cos((pi*(p(1) + p(2)))/180)*(2000*cos((pi*p(1))/180) + (sin((pi*(p(1) + 45))/180)*((2600000*cos((pi*(p(1) + p(2)))/180)*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180)))/(cos((pi*p(1))/180)^2*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180))) + (2600000*cos((pi*p(1))/180)*sin((pi*(p(1) + 45))/180))/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180)))))/((1300*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180))^2)/cos((pi*p(1))/180)^3 + (1300*sin((pi*(p(1) + 45))/180)^2)/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)^2) + 8085710131136235/4398046511104)))/sin((pi*p(2))/180)))/(3*cos((pi*p(1))/180)^2) - (539047342075749*((2600000*cos((pi*(p(1) + p(2)))/180)*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180)))/(cos((pi*p(1))/180)^2*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180))) + (2600000*cos((pi*p(1))/180)*sin((pi*(p(1) + 45))/180))/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)*(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(1))/180) - sin((pi*(p(1) + p(2)))/180)*cos((pi*p(1))/180)))))/(43980465111040*((1300*(2^(1/2)/2 + (cos((pi*(p(1) + p(2)))/180)*sin((pi*(p(1) + 45))/180))/sin((pi*p(2))/180))^2)/cos((pi*p(1))/180)^3 + (1300*sin((pi*(p(1) + 45))/180)^2)/(cos((pi*(p(1) + p(2)))/180)*sin((pi*p(2))/180)^2) + 8085710131136235/4398046511104));

Since you want to optimise it and the Optimization Toolbox functions minimise it, negating it produced a maximum. (That explains the minus sign.)

The code for the optimisation is:

A = [1 1];                          % Constraints
b = 90;                             % Constraints

[angs, fval] = fmincon(V, [1 1]', A, b)

asum = angs(1) + angs(2)            % Check

The A and b arguments are the equivalent of: (a + b) <= 90. The asum statement simply checks to see that the constraint was met.

The optimisation results are:

angs =

     -93.3393e+000
     180.0000e+000

fval =

    -694.1977e+009

asum =

      86.6607e+000

I negated V to get the maximum, so in reality fval = 694.1977e+009.

Since the solver uses the Jacobian to find the minimum of the function and reported:

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.

the constraints on the derivatives were satisfied automatically.

[Pavel]

So, there is only one set of solutions? This is a surprise...

Anyway, thank you very much sir. I couldn't find them without your help

[Star Strider]

My pleasure!

There could be more than one solution, and the one fmincon found may not be the global maximum. However to explore for those you would need the Global Optimization Toolbox. If you don’t have that, the Genetic Algorithms Toolbox from the University of Sheffield looks like a good option. A genetic algorithm approach is almost guaranteed to find a global extremum of your description, but like any such search, will likely take a while to converge.

Depending on how the GA is written, you might need to remove the minus sign when you use the GA, since in some GA designs you can tell the GA to search for the global maximum. Some search for the global minimum by default (I believe the MATLAB GA does), so leaving it in would be necessary in those situations.

The more recent versions of MATLAB have trigonometric functions in degrees (sind, cosd, atan2d etc.) that you might want to substitute for the radian-argument functions. In my experience, for degree arguments, the degree functions are more accurate than multiplying by pi/180 to get radians. That may require your rewriting your function, but it could be worthwhile.

If possible (the function space permits), I plot functions over a range of argument values. The surface I plotted for your function plotted either as a plane, or a moonscape. There could be singularities.

You can write your constraints into your fitness function, however the derivative constraints will likely be automatically satisfied, because even though the GA doesn’t use a Jacobian in its calculations, the Jacobian will likely be at a minimum at the parameter estimates the genetic algorithm converges on.

[John D'Errico]

I looked at that function. It has MANY singularities that look like they go to +/- inf. So probably no global max (or min) that is less than inf.

If you can, ALWAYS plot your functions. Look at them. This will help you to understand what is happening.