matlab solver - grid search

edited June 2007 in Science & Tech
Hi,

I have a problem with matlab.
I need a solver that can change three variables as follows:
x1 = x1+0.001
x2 = x2+0.001
x3 = x3+0.001
until another parameter y < 0.1

I now made it with 3 for loops, but this is way to slow for the interval i have to variate the variables.

Does anyone know a quicker way, like a solver, because i dont understand the matlab help concerning solvers.

Thanks in advance.

Comments

  • shwaipshwaip bluffin' with my muffin Icrontian
    edited June 2007
    can you post your code please, it'll make it easier for me to understand what you want to do.

    please copy and paste it in between [php]
    CODE HERE!
    
    [/php] tags!
  • edited June 2007
    [SIZE=2][FONT=Courier New][COLOR=#228b22]%%Constants%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/FONT][/SIZE]
    [FONT=Courier New][SIZE=2]rho0 = 1.225;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]labda = -0.0065;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]T0 = 288.15;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]g0 = 9.80665;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]R = 287.05;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Aircraft 1 parameters (index 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]S_1 = 154;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]W_1 = [0 0 -17463.49];[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Cd0_1 = 0.01232;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]AR_1 = 23;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]e_1 = 0.8765;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Aircraft 2 parameters (index 2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]S_2 = 154;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]W_2 = [0 0 -17463.49];[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Cd0_2 = 0.01232;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]AR_2 = 23;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]e_2 = 0.8765;[/FONT][/SIZE]
     
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Cable parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]sigma_ultimate = 1.1E9;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]j = 2;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]rho_c = 700;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Inputs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [SIZE=2][FONT=Courier New][COLOR=#228b22]%Jet stream inputs[/COLOR][/FONT][/SIZE]
    [FONT=Courier New][SIZE=2]H_delta = 1978.50;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]W_delta = 16.48;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]windratio = 0.7;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Aircraft 1 inputs[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]AOA_1 = 9;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Vy_1 = [0 29 0];[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Vz_1 = [0 0 0];[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]H_cruise_1 = 11000;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Kite inputs[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Vz_2 = [0 0 0];[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Cable inputs[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]d_c = 5;[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Atmosphere calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]rho_1 = rho0*(1+labda*H_cruise_1/T0)^(-(g0/(R*labda)+1));[/SIZE][/FONT]
     
    [SIZE=2][FONT=Courier New]H_cruise_2 = H_cruise_1 - H_delta;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]rho_2 = rho0*(1+labda*H_cruise_2/T0)^(-(g0/(R*labda)+1));[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Aircraft calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [SIZE=2][FONT=Courier New][COLOR=#228b22]%Lift and drag coefficient[/COLOR][/FONT][/SIZE]
    [FONT=Courier New][SIZE=2]Cl_1 = AOA_1*pi/180*2*pi*AR_1/(2+sqrt(4+AR_1^2));[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Cd_1 = Cd0_1 + Cl_1^2/(pi*AR_1*e_1);[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Speed[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Vx_1 = [W_delta*windratio 0 0];[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]V_1 = Vx_1 + Vy_1 + Vz_1;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]V_1_size = sqrt(sum(V_1.^2));[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%LOOP: Change phi_1 until Fy_1 = 0%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]phi_1 = 0;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]T(2) = 1;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]while[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2](abs(T(2)) > 0.01) [/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]phi_1 = phi_1+0.0001;[/FONT][/SIZE]
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
     
    [SIZE=2][FONT=Courier New][COLOR=#228b22]%Lift and drag[/COLOR][/FONT][/SIZE]
    [FONT=Courier New][SIZE=2]L_1_size = 0.5*rho_1*V_1_size^2*Cl_1*S_1;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]L_1 = [-L_1_size*sin(phi_1*pi/180)*Vy_1(2)/V_1_size L_1_size*sin(phi_1*pi/180)*Vx_1(1)/V_1_size L_1_size*cos(phi_1*pi/180)];[/FONT][/SIZE]
     
    [SIZE=2][FONT=Courier New]D_1_size = 0.5*rho_1*V_1_size^2*Cd_1*S_1;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]D_1 = [-D_1_size*Vx_1(1)/V_1_size -D_1_size*Vy_1(2)/V_1_size 0];[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Cable calculations 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [SIZE=2][FONT=Courier New][COLOR=#228b22]%Cable height 1[/COLOR][/FONT][/SIZE]
    [FONT=Courier New][SIZE=2]H_cable_1 = H_delta*windratio;[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Maximum cable tension[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]T_max = sigma_ultimate/j*(d_c/1000)^2*pi/4;[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Cable drag 1[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Dx_c1 = -1/3*1/2*rho_1*Vx_1.^2*d_c/1000*H_cable_1;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Dy_c1 = -1/3*1/2*rho_1*Vy_1.^2*d_c/1000*H_cable_1;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Dz_c1 = -1/3*1/2*rho_1*Vz_1.^2*d_c/1000*H_cable_1;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]D_c1 = Dx_c1 + Dy_c1 + Dz_c1;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Tension force equilibrium%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Fx_1 = [L_1(1)+D_1(1)+D_c1(1) 0 0];[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Fy_1 = [0 L_1(2)+D_1(2)+D_c1(2) 0];[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Fz_1 = [0 0 L_1(3)+ D_1(3)+ D_c1(3)+W_1(3)];[/FONT][/SIZE]
     
    [SIZE=2][FONT=Courier New]T = Fx_1 + Fy_1 + Fz_1;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]T_size = sqrt(sum(T.^2));[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%cable angle[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]beta_1 = atan(Fz_1(3)/Fx_1(1))*180/pi;[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]end[/COLOR][/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Kite calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
     
    [SIZE=2][FONT=Courier New][COLOR=#228b22]%%LOOP: vary AOA_2, Vy_2, phi_2 so that T1 = T2%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/FONT][/SIZE]
    [FONT=Courier New][SIZE=2]Vy_2 = [0 0 0];[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2] AOA_2 = -2:0.1:14[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2] i = -40:0.1:40[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Vy_2(2) = i;[/FONT][/SIZE]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2] phi_2 = -90:0.1:90[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Lift and drag coefficient[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Cl_2 = AOA_2*pi/180*2*pi*AR_2/(2+sqrt(4+AR_2^2));[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Cd_2 = Cd0_2 + Cl_2^2/(pi*AR_2*e_2);[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Speed[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Vx_2 = [W_delta-Vx_1(1) 0 0];[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]V_2 = Vx_2 + Vy_2 + Vz_2;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]V_2_size = sqrt(sum(V_2.^2));[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Lift and drag[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]L_2_size = 0.5*rho_2*V_2_size^2*Cl_2*S_2;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]L_2 = [-L_2_size*sin(phi_2*pi/180)*Vy_2(2)/V_2_size L_2_size*sin(phi_2*pi/180)*Vx_2(1)/V_2_size L_2_size*cos(phi_2*pi/180)];[/FONT][/SIZE]
     
    [SIZE=2][FONT=Courier New]D_2_size = 0.5*rho_2*V_2_size^2*Cd_2*S_2*sign(V_2_size);[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]D_2 = [-D_2_size*Vx_2(1)/V_2_size -D_2_size*Vy_2(2)/V_2_size 0];[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Cable calculations 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Cable height 1[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]H_cable_2 = H_delta - H_cable_1;[/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Cable drag 1[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Dx_c2 = -1/3*1/2*rho_2*Vx_2.^2*d_c/1000*H_cable_2;[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Dy_c2 = -1/3*1/2*rho_2*Vy_2.^2*d_c/1000*H_cable_2;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Dz_c2 = -1/3*1/2*rho_2*Vz_2.^2*d_c/1000*H_cable_2;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]D_c2 = Dx_c2 + Dy_c2 + Dz_c2;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%%Tension force equilibrium%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]Fx_2 = [L_2(1)+D_2(1)+D_c2(1) 0 0];[/SIZE][/FONT]
    [SIZE=2][FONT=Courier New]Fy_2 = [0 L_2(2)+D_2(2)+D_c2(2) 0];[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]Fz_2 = [0 0 L_2(3)+ D_2(3)+ D_c2(3)+W_2(3)];[/FONT][/SIZE]
     
    [SIZE=2][FONT=Courier New]T_2 = Fx_2 + Fy_2 + Fz_2;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]T_size_2 = sqrt(sum(T.^2));[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#228b22]%Cable angle[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2]beta_2 = atan(Fz_2(3)/Fx_2(1))*180/pi;[/SIZE][/FONT]
     
    [SIZE=2][FONT=Courier New]dFx = Fx_1(1) - Fx_2(1);[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]dFy = Fy_1(2)^2 + Fy_2(2)^2;[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]dFz = Fz_1(3) + Fz_2(3);[/FONT][/SIZE]
    [SIZE=2][FONT=Courier New]dF = dFx^2 + dFy^2 +dFz^2;[/FONT][/SIZE]
     
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2] dF < 1[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]break[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2];[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]end[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]end[/COLOR][/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2] dF < 1[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]break[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2];[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]end[/COLOR][/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]end[/COLOR][/SIZE][/FONT]
     
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2] dF < 1[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]break[/COLOR][/SIZE][/FONT][FONT=Courier New][SIZE=2];[/SIZE][/FONT]
    [FONT=Courier New][SIZE=2][COLOR=#0000ff]end[/COLOR][/SIZE][/FONT]
    [SIZE=2][FONT=Courier New][COLOR=#0000ff]end[/COLOR][/FONT][/SIZE]
    
  • shwaipshwaip bluffin' with my muffin Icrontian
    edited June 2007
    Soo...do you know if the optimization problem is convex? Because there are much better ways to do convex optimization than a linear search. I don't think there's a really quick way to do the linear search in matlab...you might be able to do it a bit faster with matricies instead of for loops, but i dunno how much that'll help.

    one easy way to do this is to change the step size on your loops (make them like .5 instead of .1). If you find a value that works, then great. If you don't, then take the space that corresponds to the smallest value, and do the search with the smaller step size.
  • edited June 2007
    I don't really know if it's convex, because I never heard about convex optimization. Is such a optimization easy to do?

    Increasing the stepsize is not really the option, because it will not find a solution then.
  • shwaipshwaip bluffin' with my muffin Icrontian
    edited June 2007
    Does the code that you gave me have a solution? I want to see if there's an error in my code.

    On thing to think about is if you really want to search for the answer over the entire space - phi is probably some angle (correct) - would you expect to get an answer that had |phi| > 80?

    This code does what is called a grid search. It steps over the space at a large step size (.5), finds the best value at that step size. Then, it looks at a range near that optimum at a smaller step size(.1), and finds the best value at that step size. This is repeated at (.01).

    There is an error in the printf statement (at least on my version of matlab), but it is working.

    [php]
    %%Constants%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    rho0 = 1.225;
    labda = -0.0065;
    T0 = 288.15;
    g0 = 9.80665;
    R = 287.05;

    %%Aircraft 1 parameters (index 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    S_1 = 154;
    W_1 = [0 0 -17463.49];
    Cd0_1 = 0.01232;
    AR_1 = 23;
    e_1 = 0.8765;

    %%Aircraft 2 parameters (index 2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    S_2 = 154;
    W_2 = [0 0 -17463.49];
    Cd0_2 = 0.01232;
    AR_2 = 23;
    e_2 = 0.8765;


    %%Cable parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    sigma_ultimate = 1.1E9;
    j = 2;
    rho_c = 700;

    %%Inputs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Jet stream inputs
    H_delta = 1978.50;
    W_delta = 16.48;
    windratio = 0.7;

    %Aircraft 1 inputs
    AOA_1 = 9;
    Vy_1 = [0 29 0];
    Vz_1 = [0 0 0];
    H_cruise_1 = 11000;

    %Kite inputs
    Vz_2 = [0 0 0];

    %Cable inputs
    d_c = 5;

    %%Atmosphere calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    rho_1 = rho0*(1+labda*H_cruise_1/T0)^(-(g0/(R*labda)+1));

    H_cruise_2 = H_cruise_1 - H_delta;
    rho_2 = rho0*(1+labda*H_cruise_2/T0)^(-(g0/(R*labda)+1));

    %%Aircraft calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Lift and drag coefficient
    Cl_1 = AOA_1*pi/180*2*pi*AR_1/(2+sqrt(4+AR_1^2));
    Cd_1 = Cd0_1 + Cl_1^2/(pi*AR_1*e_1);

    %Speed
    Vx_1 = [W_delta*windratio 0 0];
    V_1 = Vx_1 + Vy_1 + Vz_1;
    V_1_size = sqrt(sum(V_1.^2));

    %%LOOP: Change phi_1 until Fy_1 = 0%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi_1 = 0;
    T(2) = 1;
    while(abs(T(2)) > 0.01)
    phi_1 = phi_1+0.0001;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Lift and drag
    L_1_size = 0.5*rho_1*V_1_size^2*Cl_1*S_1;
    L_1 = [-L_1_size*sin(phi_1*pi/180)*Vy_1(2)/V_1_size L_1_size*sin(phi_1*pi/180)*Vx_1(1)/V_1_size L_1_size*cos(phi_1*pi/180)];

    D_1_size = 0.5*rho_1*V_1_size^2*Cd_1*S_1;
    D_1 = [-D_1_size*Vx_1(1)/V_1_size -D_1_size*Vy_1(2)/V_1_size 0];

    %%Cable calculations 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Cable height 1
    H_cable_1 = H_delta*windratio;

    %Maximum cable tension
    T_max = sigma_ultimate/j*(d_c/1000)^2*pi/4;

    %Cable drag 1
    Dx_c1 = -1/3*1/2*rho_1*Vx_1.^2*d_c/1000*H_cable_1;
    Dy_c1 = -1/3*1/2*rho_1*Vy_1.^2*d_c/1000*H_cable_1;
    Dz_c1 = -1/3*1/2*rho_1*Vz_1.^2*d_c/1000*H_cable_1;
    D_c1 = Dx_c1 + Dy_c1 + Dz_c1;

    %%Tension force equilibrium%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Fx_1 = [L_1(1)+D_1(1)+D_c1(1) 0 0];
    Fy_1 = [0 L_1(2)+D_1(2)+D_c1(2) 0];
    Fz_1 = [0 0 L_1(3)+ D_1(3)+ D_c1(3)+W_1(3)];

    T = Fx_1 + Fy_1 + Fz_1;
    T_size = sqrt(sum(T.^2));

    %cable angle
    beta_1 = atan(Fz_1(3)/Fx_1(1))*180/pi;
    end
    %kite calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%LOOP: vary AOA_2, Vy_2, phi_2 so that T1 = T2%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Vy_2 = [0 0 0];

    minA=0;
    minI=0;
    minPhi=0;
    minVal=inf;

    ARange =[-2 14];
    iRange =[-40 40];
    phiRange = [-90 90]

    for step = [.5 .1 .01]
    for AOA_2 = ARange(1):step:ARange(2)
    for i = iRange(1):step:iRange(2)
    Vy_2(2) = i;
    for phi_2 = phiRange(1):step:phiRange(2)
    %Lift and drag coefficient
    Cl_2 = AOA_2*pi/180*2*pi*AR_2/(2+sqrt(4+AR_2^2));
    Cd_2 = Cd0_2 + Cl_2^2/(pi*AR_2*e_2);

    %Speed
    Vx_2 = [W_delta-Vx_1(1) 0 0];
    V_2 = Vx_2 + Vy_2 + Vz_2;
    V_2_size = sqrt(sum(V_2.^2));

    %Lift and drag
    L_2_size = 0.5*rho_2*V_2_size^2*Cl_2*S_2;
    L_2 = [-L_2_size*sin(phi_2*pi/180)*Vy_2(2)/V_2_size L_2_size*sin(phi_2*pi/180)*Vx_2(1)/V_2_size L_2_size*cos(phi_2*pi/180)];

    D_2_size = 0.5*rho_2*V_2_size^2*Cd_2*S_2*sign(V_2_size);
    D_2 = [-D_2_size*Vx_2(1)/V_2_size -D_2_size*Vy_2(2)/V_2_size 0];

    %%Cable calculations 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Cable height 1
    H_cable_2 = H_delta - H_cable_1;

    %Cable drag 1
    Dx_c2 = -1/3*1/2*rho_2*Vx_2.^2*d_c/1000*H_cable_2;
    Dy_c2 = -1/3*1/2*rho_2*Vy_2.^2*d_c/1000*H_cable_2;
    Dz_c2 = -1/3*1/2*rho_2*Vz_2.^2*d_c/1000*H_cable_2;
    D_c2 = Dx_c2 + Dy_c2 + Dz_c2;

    %%Tension force equilibrium%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Fx_2 = [L_2(1)+D_2(1)+D_c2(1) 0 0];
    Fy_2 = [0 L_2(2)+D_2(2)+D_c2(2) 0];
    Fz_2 = [0 0 L_2(3)+ D_2(3)+ D_c2(3)+W_2(3)];

    T_2 = Fx_2 + Fy_2 + Fz_2;
    T_size_2 = sqrt(sum(T.^2));

    %Cable angle
    beta_2 = atan(Fz_2(3)/Fx_2(1))*180/pi;

    dFx = Fx_1(1) - Fx_2(1);
    dFy = Fy_1(2)^2 + Fy_2(2)^2;
    dFz = Fz_1(3) + Fz_2(3);
    dF = dFx^2 + dFy^2 +dFz^2;

    if dF < 1
    printf('done');
    break;
    elseif dF < minVal
    minA=AOA_2;
    minI=i;
    minPhi=phi_2;
    minVal = dF

    end
    end

    % if dF < 1
    % printf('done');
    % break;
    % end
    end

    % if dF < 1
    % printf('done');
    % break;
    % end
    end
    ARange = [minA-4*step minA+4*step]
    iRange = [minI-4*step minI+4*step]
    phiRange = [minPhi-4*step minPhi+4*step]
    end
    [/php]
  • edited June 2007
    hey thanks very much for your time. It is kind of working. The solution is there if the components of T are equal to the components of T_2. This is not yet the case within a certain accuracy. Decreasing the stepsize helps a bit, but not enough, and it increases calculating time a lot.
  • edited June 2007
    Ii just checked, i need at least an accuracy in stepsize of 0.00001 to get to the right solution. Do you know a way in which to get to a solution faster?
Sign In or Register to comment.