MATLAB: Optimize a 2 variable function

optimization

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

Best Answer

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.

Related Solutions

MATLAB: How to solve a non-linear equation system as function handle with fsolve

qini = zeros(1,6);             % Initial guess
ci = fsolve(@funec,qini)       % fsolve with the function handle, unknown variables and initial guess
                               % done many variations here, none of them worked
function ec = funec(x);
q1 = x(1);
q2 = x(2);
q3 = x(3);
q4 = x(4);
q5 = x(5);
q6 = x(6);
ec =[...
    sin(q6)*(cos(q1)*cos(q2)*sin(q3)*sin(q4) - cos(q4)*sin(q1) + cos(q1)*cos(q3)*sin(q2)*sin(q4)) - cos(q6)*(cos(q2 + q3)*cos(q1)*sin(q5) + cos(q5)*sin(q1)*sin(q4) + cos(q1)*cos(q2)*cos(q4)*cos(q5)*sin(q3) + cos(q1)*cos(q3)*cos(q4)*cos(q5)*sin(q2)) - 3097002709122497/18014398509481984
    cos(q6)*(cos(q1)*cos(q4) + cos(q2)*sin(q1)*sin(q3)*sin(q4) + cos(q3)*sin(q1)*sin(q2)*sin(q4)) + sin(q6)*(cos(q2 + q3)*sin(q1)*sin(q5) - cos(q1)*cos(q5)*sin(q4) + cos(q2)*cos(q4)*cos(q5)*sin(q1)*sin(q3) + cos(q3)*cos(q4)*cos(q5)*sin(q1)*sin(q2)) - 6124274932042799/36028797018963968
    - sin(q2 + q3)*cos(q5) - cos(q2 + q3)*cos(q4)*sin(q5) - 139701203459709/140737488355328
    70*cos(q1) + 360*cos(q1)*sin(q2) - 65*sin(q1)*sin(q4)*sin(q5) + 380*cos(q1)*cos(q2)*cos(q3) - 380*cos(q1)*sin(q2)*sin(q3) + 65*cos(q1)*cos(q2)*cos(q3)*cos(q5) - 65*cos(q1)*cos(q5)*sin(q2)*sin(q3) - 65*cos(q1)*cos(q2)*cos(q4)*sin(q3)*sin(q5) - 65*cos(q1)*cos(q3)*cos(q4)*sin(q2)*sin(q5) + 19869/50
    sin(q6)*(cos(q1)*cos(q4) + cos(q2)*sin(q1)*sin(q3)*sin(q4) + cos(q3)*sin(q1)*sin(q2)*sin(q4)) - cos(q6)*(cos(q2 + q3)*sin(q1)*sin(q5) - cos(q1)*cos(q5)*sin(q4) + cos(q2)*cos(q4)*cos(q5)*sin(q1)*sin(q3) + cos(q3)*cos(q4)*cos(q5)*sin(q1)*sin(q2)) + 2202239090516653/2251799813685248
    cos(q2 + q3)*cos(q6)*sin(q4) - sin(q6)*(sin(q2 + q3)*sin(q5) - cos(q2 + q3)*cos(q4)*cos(q5)) + 7591943296130745/288230376151711744
    cos(q2 + q3)*cos(q1)*cos(q5) - sin(q5)*(sin(q1)*sin(q4) + cos(q1)*cos(q2)*cos(q4)*sin(q3) + cos(q1)*cos(q3)*cos(q4)*sin(q2)) - 102014544375725/18014398509481984
    70*sin(q1) + 360*sin(q1)*sin(q2) - 380*sin(q1)*sin(q2)*sin(q3) + 380*cos(q2)*cos(q3)*sin(q1) + 65*cos(q1)*sin(q4)*sin(q5) + 65*cos(q2)*cos(q3)*cos(q5)*sin(q1) - 65*cos(q5)*sin(q1)*sin(q2)*sin(q3) - 65*cos(q2)*cos(q4)*sin(q1)*sin(q3)*sin(q5) - 65*cos(q3)*cos(q4)*sin(q1)*sin(q2)*sin(q5) - 1563/20
    cos(q6)*(sin(q2 + q3)*sin(q5) - cos(q2 + q3)*cos(q4)*cos(q5)) + cos(q2 + q3)*sin(q4)*sin(q6) - 2129854249224409/18014398509481984
    cos(q6)*(cos(q1)*cos(q2)*sin(q3)*sin(q4) - cos(q4)*sin(q1) + cos(q1)*cos(q3)*sin(q2)*sin(q4)) + sin(q6)*(cos(q2 + q3)*cos(q1)*sin(q5) + cos(q5)*sin(q1)*sin(q4) + cos(q1)*cos(q2)*cos(q4)*cos(q5)*sin(q3) + cos(q1)*cos(q3)*cos(q4)*cos(q5)*sin(q2)) - 8872946535830847/9007199254740992
    cos(q2 + q3)*cos(q5)*sin(q1) - sin(q5)*(cos(q2)*cos(q4)*sin(q1)*sin(q3) - cos(q1)*sin(q4) + cos(q3)*cos(q4)*sin(q1)*sin(q2)) - 8718732974408059/72057594037927936
    360*cos(q2) - 380*cos(q2)*sin(q3) - 380*cos(q3)*sin(q2) - 65*cos(q2)*cos(q5)*sin(q3) - 65*cos(q3)*cos(q5)*sin(q2) - 65*cos(q2)*cos(q3)*cos(q4)*sin(q5) + 65*cos(q4)*sin(q2)*sin(q3)*sin(q5) - 39571/100];
end

MATLAB: Solving non linear equation using fsolve

You are not passing theta_n and other parameters to root2d(). One fix is to make root2d a Nested Function. Also, you must choose a different initial point. phi0 = [10,10] worked for me.

function solve_it()
    beta_a = 10;
    alpha_n = 10 ;
    eta= 18;
    i = 18;
    theta_i = 10;
    theta_n = 10 ;
    
    fun = @root2d;
    phi0 = [10,10];
    phi = fsolve(fun,phi0);
    phi_n = phi(1),
    phi_i = phi(2),
    function F = root2d(phi)
    F = zeros(2,1);
    F(1) = (pi*sin((pi*(phi(1) + theta_n))/180)*cos((pi*phi(2))/180)*(cos((pi*theta_n)/180) + tan((pi*i)/180)*tan((pi*theta_i)/180)))/(180*sin((pi*phi(1))/180)*(cos((pi*(phi(1) + theta_n))/180)*cos((phi(2)*pi)/180) + sin((phi(2)*pi)/180)*tan((theta_i*pi)/180))^2) - (pi*cos((pi*phi(1))/180)*(cos((pi*theta_n)/180) + tan((pi*i)/180)*tan((pi*theta_i)/180)))/(180*sin((phi(1)*pi)/180)^2*(cos((pi*(phi(1) + theta_n))/180)*cos((pi*phi(2))/180) + sin((pi*phi(2))/180)*tan((pi*theta_i)/180)));
    F(2) = ((cos((pi*theta_n)/180) + tan((pi*i)/180)*tan((pi*theta_i)/180))*((pi*cos((pi*(phi(1) + theta_n))/180)*sin((pi*phi(2))/180))/180 - (pi*cos((pi*phi(2))/180)*tan((pi*theta_i)/180))/180))/(sin((pi*phi(1))/180)*(cos((pi*(phi(1) + theta_n))/180)*cos((phi(2)*pi)/180) + sin((phi(2)*pi)/180)*tan((theta_i*pi)/180))^2);
    end
end

Best Answer

Related Solutions

MATLAB: How to solve a non-linear equation system as function handle with fsolve

MATLAB: Solving non linear equation using fsolve

Related Question