vpasolve unable to find solution of 6 simultaneous equations

Hi there,

I am trying to use vpasolve to solve 6 simultaneous equations in which i need 7 variables solved for.

My matlab code is as follows:

Alt = 6562 * 0.3048; %Altitude (m)
M = 6300; %Aircraft mass (kg)
pos_cg = 0.29; %Centre of Gravity position 
FPA = 0 * (pi/180); %Fligth Path Angle (Radians)
g0 = 9.81; %Gravitational Constant (m/s^2)

%Valid up to 36000 ft
R = 287.5; %Gas constant (Nm/kgK)
Lr = -0.0065; %Lapse Rate (K/m)
T = 288.16 +(Lr * Alt); %Temperature (K)
Rho = 1.225 * ((T/288.16)^-((g0/(Lr * R))+1)); %Air Density (kg/m^3)
Sigma = Rho/1.225; %Density ratio

%Note that true airspeed is assumed unless otherwise stated
V_knots_vector = zeros(10,1);

for i=1:1:10
  
    V_knots = 100+(15 * i);%True airspeed range 
    V_knots_vector(i) = V_knots;

end

V_i = V_knots_vector * 0.515; %True airspeed
V_eas_i = V_knots_vector * sqrt(Sigma); %Equivalent airspeed

S = 25.08; %Wing area (m^2)
b = 15.85; %Wing span (m)
c_w = 1.716; %Wing mean chord (m)
Sweep = 0; %Sweep at 1/4Cw (deg)
Z_w = 0.45; %Z coordinate of 1/4Cw point above (-ve) or below (+ve) ox body axis (m)
alpha_wr = 1.0 * (pi/180); %Wing riggle angle (rads)

S_T = 7.79; %Tailplane area (m^2)
b_T = 6.6; %Tailplane span (m)
l_t = 6.184; %Tail arm, 1/4Cw to 1/4Ct (m)
Z_t = -1.435; %Z coordinate of 1/4Cw point above (-ve) or below (+ve) ox body axis (m)
n_T = 1.5 * (pi/180); %Tail setting angle (rad)
F_d = 1.981; %Fuselage diameter or width (m)

Z_tau = 0.312; %Thrust line z coordinate above (-ve) or below (+ve) ox body axis (m)
k = 0 * (pi/180); %Engine thrust line angle relative to ox body axis (+nose up) (rad)

%Insert values defined by the installed wing aerodynamic design 
a = 5.19; %Wing body CL-alpha (per rad)
CL_max = 1.37; %Maximum lift coefficient 
C_m0 = -0.0711; %Zero lift pitching moment
CD_0 = 0.03; %Zero lift drag coefficient
alpha_w0 = -2 * (pi/180); %Zero lift angle of attack (rad)
h_0 = -0.08; %Wing-body aero centre

%Insert values defined by the tailplane aerodynamic centre
a1 = 3.2; %Tailplane CL-alpha (per rad)
a2 = 2.414; %Elevator CL-n (per rad)
epsilon_0 = 2.0 * (pi/180); %Zero lift downwash angle (rad)

AR = b^2/S; %Aspect Ratio
s = b/2; %Wing semi-span (m)
l_T = l_t - (c_w * (pos_cg - 0.25)); %Tail arm, cg to 1/4ct (m)
V_T = (S_T * l_T)/(S * c_w); %Tail volume

x = l_t/b; %Tail position relative to wing (% of Span)
z = (Z_w - Z_t)/b; %Tail position relative to wing (% of Span)
fivector = zeros(81,1);

for fi=5:1:85

    indivfi = ((0.5 * cos((fi * pi)/180)^2)/(sqrt(x^2 + (0.5 * cos((fi * pi)/180))^2 + z^2))) * ((((x + sqrt(x^2 + (0.5 * cos((fi * pi)/180))^2 + z^2))/((0.5 * cos((fi * pi)/180))^2 + z^2)) + ((x)/(x^2 + z^2))) * pi/180);
    fivector(fi) = indivfi;

end

sumfi = sum(fivector);
d_epsilonalpha = (a/(((pi)^2) * AR))*sumfi;

%Drag Polar defined as CD = Cd_0 + KCL^2
s_d = 0.9998 + (0.0421*(F_d/b)) - (2.6286*(F_d/b)^2) + (2.000*(F_d/b)^3); %Fuselage drag factor 
K_D = (-3.333e-04 * Sweep^2) + (6.667e-05 * Sweep) + 0.38; %Emperical constant
e = 1/(pi * AR * K_D * CD_0 +(1/(0.99 * s_d))); %Oswald efficiency factor
K = 1/(pi * AR * e); %Induced drag factor

V_md = ((sqrt((2 * M * g0)/(Rho * S)))*(K/CD_0)^0.25) * (1/0.515); %Minimum drag speed (knots)
V_mdeas = V_md * sqrt(Sigma); % Equivalent minimum drag speed (knots)
V_stall = sqrt((2 * M * g0)/(Rho * S * CL_max))*(1/0.515); %Stall speed (knots)
V_stalleas = V_stall * sqrt(Sigma); %Equivalent stall speed (knots)
h_n = h_0 + (V_T * (a1/a) * (1 - d_epsilonalpha)); %Neutral point - controls fixed
K_n = h_n - pos_cg; %Static margin - controls fixed

%Trim computation finds the trim condition for each speed defined in the
%speed range table

syms V_i alpha_e C_tau C_D C_LT C_LW CL

Y = vpasolve((2*(M*g0)/Rho*(V_i^2)*S)*sin(alpha_e + FPA) == C_tau*cos(k) - C_D*cos(alpha_e) + CL*sin(alpha_e), (2*(M*g0)/Rho*(V_i^2)*S)*cos(alpha_e + FPA) == CL*cos(alpha_e) - C_D*sin(alpha_e) + C_tau*sin(k), 0 == (C_m0 + (pos_cg - h_0)*C_LW) - (V_T*C_LT) + (C_tau * (Z_tau/c_w)), CL == C_LW + (C_LT * (S_T)/S), C_D == CD_0 + K*(CL)^2, C_LW == a*(alpha_e + alpha_wr - alpha_w0), [CL,C_D,C_LT,C_LW,C_tau,V_i,alpha_e])

However i recieve an output of:

CL: [0x1 sym]

C_D: [0x1 sym]

C_LT: [0x1 sym]

C_LW: [0x1 sym]

C_tau: [0x1 sym]

V_i: [0x1 sym]

alpha_1: [0x1 sym]

It might be helpful to note that V_i is a 10 x 1 vector and i need to find the values of CL, C_D ect for each of the values of V.

Also apologies for the code length but i beleive that most of the terms are needed for the solving.

Any help/ alterations to my code would be highly appreciated. :)

[Torsten]

I count six equations for seven unknowns. So one unknown is a free variable ?

Answer 1

[Torsten]

A corrected version of the last part of your code would be

Alt = 6562 * 0.3048; %Altitude (m)
M = 6300; %Aircraft mass (kg)
pos_cg = 0.29; %Centre of Gravity position 
FPA = 0 * (pi/180); %Fligth Path Angle (Radians)
g0 = 9.81; %Gravitational Constant (m/s^2)
%Valid up to 36000 ft
R = 287.5; %Gas constant (Nm/kgK)
Lr = -0.0065; %Lapse Rate (K/m)
T = 288.16 +(Lr * Alt); %Temperature (K)
Rho = 1.225 * ((T/288.16)^-((g0/(Lr * R))+1)); %Air Density (kg/m^3)
Sigma = Rho/1.225; %Density ratio
%Note that true airspeed is assumed unless otherwise stated
V_knots_vector = zeros(10,1);
for i=1:1:10
  
    V_knots = 100+(15 * i);%True airspeed range 
    V_knots_vector(i) = V_knots;
end
V_i = V_knots_vector * 0.515; %True airspeed
V_eas_i = V_knots_vector * sqrt(Sigma); %Equivalent airspeed
S = 25.08; %Wing area (m^2)
b = 15.85; %Wing span (m)
c_w = 1.716; %Wing mean chord (m)
Sweep = 0; %Sweep at 1/4Cw (deg)
Z_w = 0.45; %Z coordinate of 1/4Cw point above (-ve) or below (+ve) ox body axis (m)
alpha_wr = 1.0 * (pi/180); %Wing riggle angle (rads)
S_T = 7.79; %Tailplane area (m^2)
b_T = 6.6; %Tailplane span (m)
l_t = 6.184; %Tail arm, 1/4Cw to 1/4Ct (m)
Z_t = -1.435; %Z coordinate of 1/4Cw point above (-ve) or below (+ve) ox body axis (m)
n_T = 1.5 * (pi/180); %Tail setting angle (rad)
F_d = 1.981; %Fuselage diameter or width (m)
Z_tau = 0.312; %Thrust line z coordinate above (-ve) or below (+ve) ox body axis (m)
k = 0 * (pi/180); %Engine thrust line angle relative to ox body axis (+nose up) (rad)
%Insert values defined by the installed wing aerodynamic design 
a = 5.19; %Wing body CL-alpha (per rad)
CL_max = 1.37; %Maximum lift coefficient 
C_m0 = -0.0711; %Zero lift pitching moment
CD_0 = 0.03; %Zero lift drag coefficient
alpha_w0 = -2 * (pi/180); %Zero lift angle of attack (rad)
h_0 = -0.08; %Wing-body aero centre
%Insert values defined by the tailplane aerodynamic centre
a1 = 3.2; %Tailplane CL-alpha (per rad)
a2 = 2.414; %Elevator CL-n (per rad)
epsilon_0 = 2.0 * (pi/180); %Zero lift downwash angle (rad)
AR = b^2/S; %Aspect Ratio
s = b/2; %Wing semi-span (m)
l_T = l_t - (c_w * (pos_cg - 0.25)); %Tail arm, cg to 1/4ct (m)
V_T = (S_T * l_T)/(S * c_w); %Tail volume
x = l_t/b; %Tail position relative to wing (% of Span)
z = (Z_w - Z_t)/b; %Tail position relative to wing (% of Span)
fivector = zeros(81,1);
for fi=5:1:85
    indivfi = ((0.5 * cos((fi * pi)/180)^2)/(sqrt(x^2 + (0.5 * cos((fi * pi)/180))^2 + z^2))) * ((((x + sqrt(x^2 + (0.5 * cos((fi * pi)/180))^2 + z^2))/((0.5 * cos((fi * pi)/180))^2 + z^2)) + ((x)/(x^2 + z^2))) * pi/180);
    fivector(fi) = indivfi;
end
sumfi = sum(fivector);
d_epsilonalpha = (a/(((pi)^2) * AR))*sumfi;
%Drag Polar defined as CD = Cd_0 + KCL^2
s_d = 0.9998 + (0.0421*(F_d/b)) - (2.6286*(F_d/b)^2) + (2.000*(F_d/b)^3); %Fuselage drag factor 
K_D = (-3.333e-04 * Sweep^2) + (6.667e-05 * Sweep) + 0.38; %Emperical constant
e = 1/(pi * AR * K_D * CD_0 +(1/(0.99 * s_d))); %Oswald efficiency factor
K = 1/(pi * AR * e); %Induced drag factor
V_md = ((sqrt((2 * M * g0)/(Rho * S)))*(K/CD_0)^0.25) * (1/0.515); %Minimum drag speed (knots)
V_mdeas = V_md * sqrt(Sigma); % Equivalent minimum drag speed (knots)
V_stall = sqrt((2 * M * g0)/(Rho * S * CL_max))*(1/0.515); %Stall speed (knots)
V_stalleas = V_stall * sqrt(Sigma); %Equivalent stall speed (knots)
h_n = h_0 + (V_T * (a1/a) * (1 - d_epsilonalpha)); %Neutral point - controls fixed
K_n = h_n - pos_cg; %Static margin - controls fixed
%Trim computation finds the trim condition for each speed defined in the
%speed range table

syms V_is alpha_e C_tau C_D C_LT C_LW CL
eqns = [2*M*g0/(Rho*V_is^2*S)*sin(alpha_e + FPA) == C_tau*cos(k) - C_D*cos(alpha_e) + CL*sin(alpha_e),...
        2*M*g0/(Rho*V_is^2*S)*cos(alpha_e + FPA) == CL*cos(alpha_e) + C_D*sin(alpha_e) + C_tau*sin(k),...
        0 == (C_m0 + (pos_cg - h_0)*C_LW) - (V_T*C_LT) + (C_tau * (Z_tau/c_w)),...
        CL == C_LW + (C_LT * S_T/S),...
        C_D == CD_0 + K*(CL)^2,...
        C_LW == a*(alpha_e + alpha_wr - alpha_w0)]

eqns = 

for i = 1:10
   sol{i} = vpasolve(subs(eqns,V_is,V_i(i)),[CL,C_D,C_LT,C_LW,C_tau,alpha_e],[0.7,0.02,0.1,0.5,0.4,0.1]);
end
sol{1}

ans = struct with fields:
         CL: 1.3736882716979893617039547243037
        C_D: 0.11409086583252592354570404895714
       C_LT: 0.37508019219059833968006975146431
       C_LW: 1.257186090790303513773023968167
      C_tau: 0.11617879574053315255268236817489
    alpha_e: 0.18987250987568138702235169855749

sol{2}

ans = struct with fields:
         CL: 1.0806515546858969868534927223957
        C_D: 0.08204078137203438620157669258433
       C_LT: 0.28138121010624787698783671680093
       C_LW: 0.99325284548622908566787677248027
      C_tau: 0.082839979599886595558852228776241
    alpha_e: 0.13901832002903891532579154524597

sol{3}

ans = struct with fields:
         CL: 0.87128570334868728148136370851425
        C_D: 0.063829347350396253038935338397756
       C_LT: 0.21520198222745664166451064099197
       C_LW: 0.80444266341440150641890207002432
      C_tau: 0.064167041145711935326839910167181
    alpha_e: 0.10263870883986211778263649660321

sol{4}

ans = struct with fields:
         CL: 0.71695898221471616548200459601452
        C_D: 0.052906622756685223827836853945781
       C_LT: 0.16681551359825705242118947185929
       C_LW: 0.6651450726879848082905745327855
      C_tau: 0.053058975039540426238551603000608
    alpha_e: 0.075799095983134430281995352626556

sol{5}

ans = struct with fields:
         CL: 0.6001210872992270641996464530754
        C_D: 0.046049082827406691891761933414645
       C_LT: 0.13040187348665011495621867784453
       C_LW: 0.55961747507988877091779065162369
      C_tau: 0.046120008794646441441601661130241
    alpha_e: 0.055466225538414577223655491131024

sol{6}

ans = struct with fields:
         CL: 0.50962278060260659047407069551502
        C_D: 0.041573643738985965218674062468029
       C_LT: 0.10232535392795833753009025351681
       C_LW: 0.47783990551892256139275478343783
      C_tau: 0.041606442819498990740642109523776
    alpha_e: 0.039709468397573303712087886470935

sol{7}

ans = struct with fields:
         CL: 0.43813998390860542567459309107445
        C_D: 0.038554570690162372329899563067032
       C_LT: 0.080226789612682063372447793971915
       C_LW: 0.41322105683193902720284794294681
      C_tau: 0.038568898986627967418968292595685
    alpha_e: 0.027258823178501331806518167686537

sol{8}

ans = struct with fields:
         CL: 0.38071407087365265742599012845424
        C_D: 0.036459074592921189740437247728118
       C_LT: 0.062523816276295533606510464746782
       C_LW: 0.36129379460601540835124066591925
      C_tau: 0.036464501938483477653475319872306
    alpha_e: 0.017253570341136472682894414887394

sol{9}

ans = struct with fields:
         CL: 0.33389672563068869193839754580555
        C_D: 0.034968173612232022878020538591145
       C_LT: 0.048124049689892834645375165653509
       C_LW: 0.31894910413610076602581889586772
      C_tau: 0.034969619823878641335193703984335
    alpha_e: 0.0090946704432723797492871374208174

sol{10}

ans = struct with fields:
         CL: 0.295232354057315935263822240204
        C_D: 0.033884188895691379596728113837506
       C_LT: 0.036254125283240502670705932366935
       C_LW: 0.28397160302237002155549691272639
      C_tau: 0.033884282878537438179818669745467
    alpha_e: 0.0023552675311855311037530366228743

However, vpasolve is not able to find a solution.

[Walter Roberson]

When you have a single output then the only real reason to pass a list of variables to vpasolve is if you continue on to pass a list of initial values. In cases such as this one if your list of variables does not include all of the variables in the equations then vpasolve will only be able to solve if the remaining variables drop out because equations are redundant.

[Torsten]

V_i is now substituted by a numerical value such that vpasolve should have a system of 6 equations in 6 unknowns. But it is not able to solve it.

I don't get what you mean with your comment.

[Torsten]

Doesn't seem to help ( see the complete code above ).

[Torsten]

I corrected the equations and it works now (see above).

Answer 2

[Walter Roberson]

You need to arrayfun or loop over elements of V_i as vpasolve and solve try to solve all of the passed equations simultaneously

[Joel Booth]

Thanks for the help, how exactly would i do this?

Im not familiar with the arrayfun function or as to how looping would work